MATLAB的bsxfun是最好的嗎？ Python的numpy.einsum？

我有一個非常大的乘法和求和操作，我需要儘可能高效地實現。到目前爲止，我已經找到了最好的方法是在MATLAB，在那裏我制訂了問題，因爲bsxfun：MATLAB的bsxfun是最好的嗎？ Python的numpy.einsum？

L = 10000; 
x = rand(4,1,L+1); 
A_k = rand(4,4,L); 
tic 
for k = 2:L 
    i = 2:k; 
    x(:,1,k+1) = x(:,1,k+1)+sum(sum(bsxfun(@times,A_k(:,:,2:k),x(:,1,k+1-i)),2),3); 
end 
toc

注意L將在實踐中更大。有更快的方法嗎？奇怪的是，我需要首先將單身人士維度添加到x，然後sum，但我不能讓它工作，否則。

它仍然比我試過任何其他方法快得多，但還不足以對我們的應用程序。我聽說過有關Python函數numpy.einsum可能更高效的傳聞，但在考慮移植我的代碼之前，我想先問一下。

我使用MATLAB R2017b。

來源

2017-10-19 Alex

for循環的許多次比bsxfun快。但你似乎已經嘗試過......？ – Adiel

矩陣乘法比'bsxfun'快，但這似乎並不能很容易地借給自己是 –

@LuisMendo我認爲改造問題到稀疏矩陣乘法，但它是一個有點嚇人......你以爲這將是更快？ – Alex

相信無論你的求和可以被刪除，但我只取出暫時就越容易之一。在第二尺寸的總和是微不足道的，因爲它僅影響A_k數組：

B_k = sum(A_k,2); 
for k = 2:L 
    i = 2:k; 
    x(:,1,k+1) = x(:,1,k+1) + sum(bsxfun(@times,B_k(:,1,2:k),x(:,1,k+1-i)),3); 
end

利用該單個更改運行時從〜8秒〜我的筆記本電腦2.5秒降低。

第二求和也可以通過轉化倍+總結成矩陣矢量乘積被去除，。它需要一些單身擺弄來獲得尺寸正確的，但如果你定義一個輔助陣列是B_k與第二維逆轉，您可以生成餘額爲〜x*C_k這種輔助陣列C_k，或採取了幾個電話給reshape 。

所以仔細一看後，我意識到，我原來的評估過於樂觀：你必須在你剩餘期限兩個維度乘法，所以它不是一個簡單的矩陣產品。無論如何，我們可以將該術語改寫爲矩陣產品的對角線。這意味着，我們計算了一堆不必要的矩陣元素，但這似乎仍比bsxfun方法稍快，我們可以擺脫討厭的單維度的過於：

L = 10000; 
x = rand(4,L+1); 
A_k = rand(4,4,L); 
B_k = squeeze(sum(A_k,2)).'; 

tic 
for k = 2:L 
    ii = 1:k-1; 
    x(:,k+1) = x(:,k+1) + diag(x(:,ii)*B_k(k+1-ii,:)); 
end 
toc

這個運行在〜在我的筆記本電腦上2.2秒，比之前獲得的約2.5秒稍快。

來源

2017-10-19 18:56:07

哇！這也是我筆記本電腦的4倍。謝謝！ – Alex

@Alex只是確保它是正確的，我只是匆匆一瞥，並沒有嚴格檢查，因爲所有'NaN'和'Inf's。我想我可以擺脫你的第二筆錢，這將進一步提高性能，我現在沒有時間:) –

@Alex所以我刪除了第二筆款項，它只是比以前的版本略快，但它可能對更大的問題更好地擴展。 –

由於您使用Matlab的一個新版本，你可以嘗試廣播/ implicit expansion代替bsxfun：

x(:,1,k+1) = x(:,1,k+1)+sum(sum(A_k(:,:,2:k).*x(:,1,k-1:-1:1),3),2);

我也改變了求和的順序，並取消了進一步改善i變量。在我的機器上，使用Matlab R2017b，這個速度比L = 10000快了25％。

來源

2017-10-19 18:50:49 horchler

@LuisMendo：絕對。 25％的改善中有一半是由於改變求和順序。 – horchler

看起來，隱式版本比'bsxfun'稍快。 – Alex

@horchler哦，我沒有注意到那部分 –

MATLAB的bsxfun是最好的嗎？ Python的numpy.einsum？

回答

相關問題