5

使用 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}]

所以我搜索了一个 fold 函数，我发现的只是saccumulate上的一个方法numpy.ufunc。老实说，我知道我可能注定要失败，因为尝试

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

如numpy 邮件列表中所述产生错误

Reduction not defined on ufunc with signature

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

提前致谢。

4

1 回答 1

4

值得深思的是，这里有 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索引表达式看起来像一个操作序列，但它实际上是作为一个 3 轴求和的 5d 乘积来计算的。在 C 代码中这是用 annditer和 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不能采用任意数量的ops. 此外，迭代空间随着每个操作和索引的增加而增加。但我可以想象在 Cython 级别或从 Python 反复调用它。我认为它有可能比repeat(dot...)许多小型阵列更快。

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]
于 2015-01-17T05:27:18.533 回答