4

括弧を最適に配置するなどして、大規模な疎行列の計算を高速化することは可能ですか?

私が求めているのは、Matlab に指定された順序 (たとえば、「右から左」など) で操作を実行させることで、次のコードを高速化できますか?

以前に因数分解されたスパース正方対称行列 H と、H の次元に等しい長さのスパース ベクトル M があります。やりたいことは次のとおりです。

編集:追加情報: H は通常 4000x4000 です。z と c の計算は約 4000 回実行されますが、dVa と dVaComp の計算は 4000 ループごとに 10 回実行されるため、合計で 40000 回実行されます。(dVa と dVaComp は反復的に解決され、P_mis が更新されます)。

ここM*c*M'で は、4 つの非ゼロ要素を持つスパース行列になります。マトラブでは：

[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 の逆関数を何度も使用していますが、事前に計算する方が高速ではありません。

ありがとう =)

4

2 に答える 2

1

私には完全に明確ではないことがいくつかあります。

  • コマンドM = sparse([t f],[1 -1],1,n,1);は正しくありません。t,f行と列には;1,-1が必要だと言っています。1列-1は明らかに正しくありません。
  • 結果dVaCompは による乗算による完全な行列P_misですが、スパースである必要があります。

今のところこれらの問題は脇に置いておきますが、いくつかの小さな最適化が見られます。

  • 2 回使用する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-13T12:17:58.443 に答える