如何折叠/累积一个numpy矩阵乘积(点)? [英] How to fold/accumulate a numpy matrix product (dot)?

问题描述

使用python库numpy，可以使用函数 cumprod 评估累积产品，例如

With using python library numpy, it is possible to use the function cumprod to evaluate cumulative products, e.g.

a = np.array([1,2,3,4,2])
np.cumprod(a)

给予

array([ 1,  2,  6, 24, 48])

确实可以仅在一个轴上应用此功能.

It is indeed possible to apply this function only along one axis.

我想对矩阵(表示为numpy数组)做同样的事情，例如如果我有

I would like to do the same with matrices (represented as numpy arrays), e.g. if I have

S0 = np.array([[1, 0], [0, 1]])
Sx = np.array([[0, 1], [1, 0]])
Sy = np.array([[0, -1j], [1j, 0]])
Sz = np.array([[1, 0], [0, -1]])

和

b = np.array([S0, Sx, Sy, Sz])

然后我想要一个类似cumprod的函数，该函数可以提供

then I would like to have a cumprod-like function which gives

np.array([S0, S0.dot(Sx), S0.dot(Sx).dot(Sy), S0.dot(Sx).dot(Sy).dot(Sz)])

(这是一个简单的示例，实际上，我可能会在n维网格上评估大型矩阵，因此我寻求评估该问题的最简单和方法.)

(This is a simple example, in reality I have potentially large matrices evaluated over n-dimensional meshgrids, so I seek for the most simple and efficient way to evaluate this thing.)

例如我会使用

FoldList[Dot, IdentityMatrix[2], {S0, Sx, Sy, Sz}]

所以我搜索了折叠函数，发现的只是numpy.ufunc上的accumulate方法.老实说，我知道我可能因为尝试

so I searched for a fold function, and all I found is an accumulate method on numpy.ufuncs. To be honest, I know that I am probably doomed because an attempt at

np.core.umath_tests.matrix_multiply.accumulate(np.array([pauli_0, pauli_x, pauli_y, pauli_z]))

如一个numpy邮件列表中所述的

错误

as mentioned in a numpy mailing list yields the error

Reduction not defined on ufunc with signature

您知道如何(有效)进行这种计算吗?

Do you have an idea how to (efficiently) do this kind of calculation ?

谢谢.

推荐答案

值得深思的是，这里有3种方法来评估3种连续点积:

As food for thought, here are 3 ways of evaluating the 3 sequential dot products:

使用普通的Python reduce(也可以写成循环)

With the normal Python reduce (which could also be written as a loop)

In [118]: reduce(np.dot,[S0,Sx,Sy,Sz])
array([[ 0.+1.j,  0.+0.j],
       [ 0.+0.j,  0.+1.j]])

等效于einsum

In [119]: np.einsum('ij,jk,kl,lm',S0,Sx,Sy,Sz)

einsum索引表达式看起来像一个操作序列，但实际上被评估为3d轴求和的5d乘积.在C代码中，这是通过nditer完成并大步前进的，但效果如下:

The einsum index expression looks like a sequence of operations, but it is actually evaluated as a 5d product with summation on 3 axes. In the C code this is done with an nditer and strides, but the effect is as follows:

In [120]: np.sum(S0[:,:,None,None,None] * Sx[None,:,:,None,None] *
    Sy[None,None,:,:,None] * Sz[None,None,None,:,:],(1,2,3))

In [127]: np.prod([S0[:,:,None,None,None], Sx[None,:,:,None,None],
    Sy[None,None,:,:,None], Sz[None,None,None,:,:]]).sum((1,2,3))

---

前一段时间，我从np.einsum创建补丁时，将C代码转换为Python，还编写了Cython产品和功能.这段代码在github上

---

A while back while creating a patch from np.einsum I translated that C code to Python, and also wrote a Cython sum-of-products function(s). This code is on github at

https://github.com/hpaulj/numpy-einsum

einsum_py.py是Python einsum，具有一些有用的调试输出

einsum_py.py is the Python einsum, with some useful debugging output

sop.pyx是Cython代码，已编译为sop.so.

sop.pyx is the Cython code, which is compiled to sop.so.

这是将其用于部分问题的方法.我跳过了Sy数组，因为我的sop没有为复数编码(但是可以更改).

Here's how it could be used for part of your problem. I'm skipping the Sy array since my sop is not coded for complex numbers (but that could be changed).

import numpy as np
import sop
import einsum_py    

S0 = np.array([[1., 0], [0, 1]])
Sx = np.array([[0., 1], [1, 0]])
Sz = np.array([[1., 0], [0, -1]])

print np.einsum('ij,jk,kl', S0, Sx, Sz)
# [[ 0. -1.] [ 1.  0.]]
# same thing, but with parsing information
einsum_py.myeinsum('ij,jk,kl', S0, Sx, Sz, debug=True)
"""
{'max_label': 108, 'min_label': 105, 'nop': 3, 
 'shapes': [(2, 2), (2, 2), (2, 2)], 
 'strides': [(16, 8), (16, 8), (16, 8)], 
 'ndim_broadcast': 0, 'ndims': [2, 2, 2], 'num_labels': 4,
 ....
 op_axes [[0, -1, 1, -1], [-1, -1, 0, 1], [-1, 1, -1, 0], [0, 1, -1, -1]]
"""    

# take op_axes (for np.nditer) from this debug output
op_axes = [[0, -1, 1, -1], [-1, -1, 0, 1], [-1, 1, -1, 0], [0, 1, -1, -1]]
w = sop.sum_product_cy3([S0,Sx,Sz], op_axes)
print w

如sum_product_cy3所述，不能取任意数量的ops.再加上迭代空间随每个操作和索引而增加.但是我可以想象在Cython级别或从Python重复调用它.我认为对于许多小型阵列，它有可能比repeat(dot...)更快.

As written sum_product_cy3 cannot take an arbitrary number of ops. Plus the iteration space increases with each op and index. But I can imagine calling it repeatedly, either at the Cython level, or from Python. I think it has potential for being faster than repeat(dot...) for lots of small arrays.

Cython代码的精简版本是:

A condensed version of the Cython code is:

def sum_product_cy3(ops, op_axes, order='K'):
    #(arr, axis=None, out=None):
    cdef np.ndarray[double] x, y, z, w
    cdef int size, nop
    nop = len(ops)
    ops.append(None)
    flags = ['reduce_ok','buffered', 'external_loop'...]
    op_flags = [['readonly']]*nop + [['allocate','readwrite']]

    it = np.nditer(ops, flags, op_flags, op_axes=op_axes, order=order)
    it.operands[nop][...] = 0
    it.reset()
    for x, y, z, w in it:
        for i in range(x.shape[0]):
           w[i] = w[i] + x[i] * y[i] * z[i]
    return it.operands[nop]

这篇关于如何折叠/累积一个numpy矩阵乘积(点)?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！