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

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

给

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

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

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

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的功能

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

(这是一个简单的例子,实际上我有可能通过n维网格物体评估大型矩阵,所以我寻求最简单有效的方法来评估这个东西.)

在例如Mathematica我会用

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

所以我搜索了一个折叠函数,我发现的是numpy.ufuncs上的累积方法.说实话,我知道我可能注定要失败,因为我是在尝试

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

如a numpy mailing list中所述,产生错误

Reduction not defined on ufunc with signature

你知道如何(有效地)进行这种计算吗？

提前致谢.

解决方法:

作为思考的食物,这里有3种评估3个连续点产品的方法：

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

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索引表达式看起来像一系列操作,但它实际上被评估为5d乘积,在3轴上求和.在C代码中,这是通过nditer和strides完成的,但效果如下：

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上

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

einsum_py.py是Python einsum,带有一些有用的调试输出

sop.pyx是Cython代码,编译为sop.so.

以下是它如何用于部分问题.我正在跳过Sy数组,因为我的sop没有编码为复数(但可以改变).

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不能采取任意数量的操作.此外,迭代空间随着每个操作和索引而增加.但我可以想象,无论是在Cython级别还是从Python中重复调用它.我认为对于许多小型阵列来说,它有可能比重复(点…)更快.

Cython代码的精简版本是：

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]