Как складывать/накапливать матричный продукт (точка)?

Question

С использованием функции numpy библиотеки python можно использовать функцию 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}] 

так что я искал кратное функции, и все, что я нашел это accumulate метод на numpy.ufunc с. Честно говоря, я знаю, что я, вероятно, обречены, потому что попытка

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 

Есть ли у вас представление о том, как (эффективно) делают этот вид расчет?

Заранее спасибо.

Answer 1

В пищу для размышлений, вот 3 способа оценки 3 последовательных точечных продуктов:

с нормальным Python уменьшить (которые также могут быть записаны в виде петли)

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 остроумия продукта h суммирования по 3 осям. В коде C это делается с nditer и махов, но эффект заключается в следующем:

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, а также написал a Cython функция (-ы) суммы продуктов.Этот код на GitHub в

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

einsum_py.py является einsum Python, с некоторой полезной отладочный вывод

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]