如何做總和的平方和的總和？

Question

我有一筆我想要加快的款項。在一種情況下，它是：

S_ {X，Y，K，L} Fu_ {ク} Fv_ {LV} Fx_ {KX} Fy_ {LY}

在其他情況下，它是：

S_ {X，Y}（S_ {K，L} Fu_ {ク} Fv_ {LV} Fx_ {KX} Fy_ {LY}）^ 2

注：S_ {指數}：是對這些總和指數

第一個案例我已經想出瞭如何使用numpy的einsum，它導致了驚人的加速〜x160。

此外，我曾試圖擴大平方，但不會是一個殺手，因爲我需要總結x，y，k，l，k，l而不是x，y，k，l ？

這是一個實現，它演示了我與einsum的區別和解決方案。

Nx = 3 
Ny = 4 
Nk = 5 
Nl = 6 
Nu = 7 
Nv = 8 
Fx = np.random.rand(Nx, Nk) 
Fy = np.random.rand(Ny, Nl) 
Fu = np.random.rand(Nu, Nk) 
Fv = np.random.rand(Nv, Nl) 
P = np.random.rand(Nx, Ny) 
B = np.random.rand(Nk, Nl) 
I1 = np.zeros([Nu, Nv]) 
I2 = np.zeros([Nu, Nv]) 
t = time.time() 
for iu in range(Nu): 
    for iv in range(Nv): 
     for ix in range(Nx): 
      for iy in range(Ny): 
       S = 0. 
       for ik in range(Nk): 
        for il in range(Nl): 
         S += Fu[iu,ik]*Fv[iv,il]*Fx[ix,ik]*Fy[iy,il]*P[ix,iy]*B[ik,il] 
       I1[iu, iv] += S 
       I2[iu, iv] += S**2. 
print time.time() - t; t = time.time() 
# 0.0787379741669 
I1_ = np.einsum('uk, vl, xk, yl, xy, kl->uv', Fu, Fv, Fx, Fy, P, B) 
print time.time() - t 
# 0.00049090385437 
print np.allclose(I1_, I1) 
# True 
# Solution by expanding the square (not ideal) 
t = time.time() 
I2_ = np.einsum('uk,vl,xk,yl,um,vn,xm,yn,kl,mn,xy->uv', Fu,Fv,Fx,Fy,Fu,Fv,Fx,Fy,B,B,P**2) 
print time.time() - t 
# 0.0226809978485 <- faster than for loop but still much slower than I1_ einsum 
print np.allclose(I2_, I2) 
# True 

如圖所示我已經成功地做到I1_與我已經找到了如何做到以上einsum爲I1。

編輯：

我加了如何通過擴大正方形，但速度可達做I2_是有點令人失望，也值得期待...〜x3.47增速相比，X160〜

編輯2：

加速似乎並不一致，我得到了x40和x1.2之前，但現在得到不同的數字。無論哪種方式，差異和問題依然存在。

EDIT3： 我試圖簡化實際上我不過之後搞砸及以上的總和允許通過@ user5402提供的出色答卷的總和。

我已經編輯上面的代碼以證明其是下面的總和：

I1 = S_ {X，Y，K，L} Fu_ {ク} Fv_ {LV} Fx_ {KX} Fy_ {LY } P_ {xy} B_ {kl}

I2 = S_ {x，y} {S_ {k，l} Fu_ {ku} Fv_ {lv} Fx_ {kx} Fy_ {ly} P_ {xy} B_ { kl}）^ 2

Answer 1

由於問題已更改，我將開始一個新答案。

試試這個：

E = np.einsum('uk, vl, xk, yl, xy, kl->uvxy', Fu, Fv, Fx, Fy, P, B) 
E1 = np.einsum('uvxy->uv', E) 
E2 = np.einsum('uvxy->uv', np.square(E)) 

我發現它運行得一樣快的時間I1_。

這是我的測試代碼：http://pastebin.com/ufwy7cLy

Answer 2

（更新：跳轉到最後查看錶示爲幾個矩陣乘法的結果。）

我想你可以通過使用身份大大簡化了計算：

例如，

S_{k,l} Fu_{ku} Fv_{lv} Fx_{kx} Fy_{ly} 
    = S_{k,l} Fu_{ku} Fx_{kx} Fv_{lv} Fy_{ly}   -- rearrange the factors 
      \___ A ____/ \___ B ____/ 
    = (S_k Fu_{ku} Fx_{kx}) * (S_l Fv_{lv} Fy_{ly}) -- from the identity 
    = A_{ux}    * B_{vy} 

其中A_{ux}只取決於u和x和B_{vy}只取決於v和y。

對於平方和，我們有：

S_{xy} A_{ux} * B_{vy} 
    = S_x A_{ux} * S_y B_{vy}      -- from the identity 
    = C_u  * D_v 

然後終於總結了u和v：

繼續總和當過 x和 y發生

S_k [ S_l Fu_{ku} Fv_{lv} Fx_{kx} Fy_{ly} ]^2 
    = S_k Fu_{ku} Fx_{kx} * [ S_l Fv_{lv} Fy_{ly} ]^2 
    = S_k Fu_{ku} Fx_{kx} * B_{vy}^2     -- B is from the above calc. 
    = B_{vy}^2 * S_k Fu_{ku} Fx_{kx}     -- B_vy is free of k 
    = B_{vy}^2 * A_{ux}        -- A is from the above calc. 

類似的削減

S_{uv} C_u D_v = (S_u C_u) * (S_v D_v)    -- from the identity 

希望這有助於。

更新：我剛剛意識到也許平方和你想計算 [ S_k S_l ... ]^2在這種情況下，你可以這樣進行：

[ S_k S_l Fu_{ku} Fv_{lv} Fx_{kx} Fy_{ly} ]^2 
    = [ A_{ux}    * B_{vy} ]^2 
    = A_{ux}^2 * B_{vy}^2 

所以，當我們總結過了兩個變量，我們得到：

S_{uvxy} A_{ux}^2 B_{vy}^2 
    = S_{uv} (S_{xy} A_{ux}^2 B_{vy}^2) 
    = S_{uv} (S_x A_{ux}^2) * (S_y B_{vy}^2)  -- from the identity 
    = S_{uv}  C_u   *  D_v 
    = (S_u C_u) * (S_v D_v)       -- from the identity 

更新2：這可以歸結爲幾個矩陣乘法。

A和的定義B：

A_{uv} = S_k Fu_{ku} Fx_{kx} 
B_{vy} = S_l Fv_{lv} Fy_{ly} 

也可以寫成矩陣形式：

A = (transpose Fu) . Fx    -- . = matrix multiplication 
B = (transpose Fv) . Fy 

和C和d的定義：

C_u = S_x A_{ux} 
D_v = S_y B_{vy} 

我們看到向量C只是A的行和，向量D只是B的行和。由於整個求和的答案（未平方）是：

total = (S_u C_u) * (S_v D_v) 

我們看到總畢竟是一個時代的矩陣元素的簡單相加所有B.

這裏的矩陣元素的總和是numpy的代碼：

from numpy import * 
# ... set up Fx, Fv, Fu, Fy as above... 

A = Fx.dot(Fu.transpose()) 
B = Fv.dot(Fy.transpose()) 
sum1 = sum(A) * sum(B) 

A2 = square(A) 
B2 = square(B) 
sum2 = sum(A2) * sum(B2) 

print "sum of terms:", sum1 
print "sum of squares of terms:", sum2