加快稀疏矩陣計算

是否有可能通過例如加速大型稀疏矩陣計算來加速稀疏矩陣計算。最佳地放置偏食？加快稀疏矩陣計算

我問的是：我可以通過強制Matlab以指定的順序執行操作（例如「從右到左」或類似的東西）來加速以下代碼嗎？

我有一個稀疏正方形對稱矩陣H，先前已被因式分解，並用長度等於H的尺寸的稀疏矢量M我想要做的是以下內容：

編輯：一些額外的信息：H通常是4000x4000。 z和c的計算大約需要4000次，而dVa和dVaComp的計算每4000個循環完成10次，因此總計爲40000次。（迭代地求解dVa和dVaComp，其中P_mis被更新）。

這裏M*c*M'將變成具有4個非零元素的稀疏矩陣。在Matlab中：

[L U P] = lu(H);     % H is sparse (thus also L, U and P) 
% for i = 1:4000    % Just to illustrate 
M = sparse([bf bt],1,[1 -1],n,1); % Sparse vector with two non-zero elements in bt and bf 
z = -M'*(U \ (L \ (P * M)));  % M^t*H^-1*M = a scalar 
c = (1/dyp + z)^-1;    % dyp is a scalar 
    % while (iterations < 10 && ~=converged) 
    dVa = - (U \ (L \ (P * P_mis))); 
    dVaComp = (U \ (L \ (P * M * c * M' * dVa))); 
    % Update P_mis etc. 
    % end while 
% end for

並記錄：儘管我多次使用H的倒數，但預先計算它並不會更快。

[L,U,P,Q] = lu(H);

額外的置換矩陣Q重新排序的列增加：

來源

2013-05-12 Stewie Griffin

有幾個問題：你說「'P'是一個預定義的完整向量」，但是你將它生成爲一個與'H'大小相同的稀疏矩陣。另外，在「數學」「z」和實現的「z」之間，有一個負號不同。這使得你不清楚你想要完成什麼;你能否糾正/澄清這些問題？ – 2013-05-13 08:17:39

感謝您指出Rody的錯誤。「預定義」P現在變爲P_mis，並且是在找到dVa和dVaComp時針對每次迭代更新的失配向量。不要與置換矩陣P混淆。它應該是z前面的負號。 – 2013-05-13 08:31:33

有一點讓你感到困惑，你使用t一次換位，一次使用變量。順便說一句。我記得黑暗時，如果你計算M * M'，Matlab有一些優化 - 因爲c是標量，你應該可以在dVaComp中切換它。 – bdecaf 2013-05-13 13:05:16

有幾件事情並不完全清楚對我說：

命令M = sparse([t f],[1 -1],1,n,1);不可能是正確的;你是說在行t,f和列1,-1應該有一個1;列-1顯然不能正確。
結果dVaComp是一個完整的矩陣，由於乘以P_mis，而你說它應該是稀疏的。

離開這些問題放在一邊，有一些小的優化我看到：

您使用inv(H)*M兩次，所以你可以預先計算這一點。
否定的dVa可以移出循環。
如果您不需要明確地指定dVa，那麼也不要指定給變量。
標量的反轉意味着將該標量除以1（計算c）。

實施的變化，並試圖相當比較（我只用40次迭代，以保持總時間小）：

%% initialize 
clc 
N = 4000; 

% H is sparse, square, symmetric 
H = tril(rand(N)); 
H(H<0.5) = 0; % roughly half is empty 
H = sparse(H+H.'); 

% M is sparse vector with two non-zero elements. 
M = sparse([1 N],[1 1],1, N,1); 

% dyp is some scalar 
dyp = rand; 

% P_mis = full vector 
P_mis = rand(N,1); 


%% original method 

[L, U, P] = lu(H); 

tic    

for ii = 1:40 

    z = -M'*(U \ (L \ (P*M))); 
    c = (1/dyp + z)^-1; 

    for jj = 1:10   
     dVa = -(U \ (L \ (P*P_mis))); 
     dVaComp = (U \ (L \ (P*M * c * M' * dVa)));  
    end 

end 

toc 


%% new method I 

[L,U,P,Q] = lu(H);  

tic    

for ii = 1:40 

    invH_M = U\(L\(P*M)); 

    z = -M.'*invH_M; 
    c = -1/(1/dyp + z); 

    for jj = 1:10   
     dVaComp = c * (invH_M*M.') * (U\(L\(P*P_mis))); 
    end 

end 

toc

這得出以下結果：

Elapsed time is 60.384734 seconds. % your original method 
Elapsed time is 33.074448 seconds. % new method

來源

2013-05-13 12:17:58

你可能想融通（稀疏）矩陣H時使用擴展語法lu嘗試

感謝=）因子的稀疏性L,U（而排列矩陣P僅重排行部分轉軸）。

的具體結果取決於H稀疏圖案，但在許多情況下，使用良好的列置換顯著減少了在因式分解非零的數目，減少了存儲器的使用和增加的速度。

您可以閱讀更多關於lu語法here的信息。

來源

2013-05-13 09:56:17

加快稀疏矩陣計算

回答

相關問題