加快稀疏矩阵计算答案

【问题标题】：Speed up sparse matrix calculations加快稀疏矩阵计算
【发布时间】：2013-05-12 16:05:14
【问题描述】：

是否可以通过例如加速大型稀疏矩阵计算以最佳方式放置花括号？

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

谢谢 =)

【问题讨论】：

有几个问题：你说“P 是一个预定义的全向量”，但是你将它生成为一个与H 大小相同的稀疏矩阵。此外，在“数学”z 和实现的z 之间，有一个减号不同。这让人有点不清楚你想要完成什么；您能否更正/澄清这些问题？
感谢罗迪指出错误。 “预定义” P 现在更改为 P_mis，并且是一个失配向量，当找到 dVa 和 dVaComp 时，每次迭代都会更新。不要与置换矩阵P混淆。它应该是z前面的减号。
将 t 用于转置一次，将其用作变量有点令人困惑。顺便提一句。我记得很清楚，如果你计算 M*M'，Matlab 会有一些优化 - 因为 c 是标量，你应该能够在 dVaComp 中切换它。

标签： performance matlab sparse-matrix

【解决方案1】：

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

M = sparse([t f],[1 -1],1,n,1); 命令不正确；你是说在行t,f 和列1,-1 上应该有一个1；专栏-1 显然不可能是对的。
由于乘以P_mis，结果dVaComp 是一个完整的矩阵，而您说它应该是稀疏的。

暂时把这些问题放在一边，我看到了一些小的优化：

您使用了两次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

【讨论】：

【解决方案2】：

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

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

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

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

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

【讨论】：