4

是否可以通过优化放置括号来加速大型稀疏矩阵计算？

我要问的是：我可以通过强制 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 的倒数，预计算它也并不快。

谢谢 =)

4

2 回答 2

1

有几件事我并不完全清楚：

  • 命令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-13T12:17:58.443 回答
1

lu在分解（稀疏）矩阵时，您可能想尝试使用扩展语法H：

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

额外的置换矩阵Q对列进行重新排序以增加因子的稀疏性L,U（而置换矩阵P仅对部分旋转的行进行重新排序）。

具体结果取决于 的稀疏模式H，但在许多情况下，使用良好的列排列会显着减少因式分解中非零的数量，从而减少内存使用并提高速度。

您可以在此处阅读有关lu语法的更多信息。

于 2013-05-13T09:56:17.800 回答