在 MATLAB 中循环添加多个稀疏矩阵的最快方法答案

【问题标题】：Fastest way to add multiple sparse matrices in a loop in MATLAB在 MATLAB 中循环添加多个稀疏矩阵的最快方法
【发布时间】：2015-04-15 19:21:19
【问题描述】：

我有一个代码在循环中重复计算稀疏矩阵（准确地说，它执行此计算 13472 次）。这些稀疏矩阵中的每一个都是唯一的。

每次执行后，它将新计算的稀疏矩阵添加到原来的稀疏零矩阵中。

当所有 13742 个矩阵都添加完毕后，代码退出循环，程序终止。

代码瓶颈出现在添加稀疏矩阵中。我制作了一个虚拟版本的代码，它与我的真实代码表现出相同的行为。它由一个 MATLAB 函数和下面给出的脚本组成。

（1）生成稀疏矩阵的函数：

function out = test_evaluate_stiffness(n)
ind = randi([1 n*n],300,1);
val = rand(300,1);
[I,J] = ind2sub([n,n],ind);
out = sparse(I,J,val,n,n);
end

(2) 主脚本（程序）

% Calculate the stiffness matrix
n=1000;
K=sparse([],[],[],n,n,n^2);
tic
for i=1:13472
    temp=rand(1)*test_evaluate_stiffness(n);
    K=K+temp;
end
fprintf('Stiffness Calculation Complete\nTime taken = %f s\n',toc)

我对稀疏矩阵运算不是很熟悉，所以我可能在这里遗漏了一个关键点，这可能会大大加快我的代码速度。

我是否在代码中以合理的方式处理刚度矩阵的更新？是否有另一种方法可以让我使用 sparse 来获得更快的解决方案？

下面还提供了探查器报告：

【问题讨论】：

您需要 13742 个临时矩阵还是它们的总和？
中只需要总和。

标签： matlab sparse-matrix

【解决方案1】：

如果您只需要这些矩阵的总和，而不是单独构建所有矩阵然后求和它们，只需连接向量 I、J 和 vals 并只调用一次 sparse。如果[I,J]中有重复行[i,j]，则对应的值S(i,j)会自动求和，所以代码是完全等价的。由于调用sparse 涉及对排序算法的内部调用，因此您节省了 13742-1 个中间排序，并且只需要一个即可。

这涉及将test_evaluate_stiffness 的签名更改为输出[I,J,val]：

function [I,J,val] = test_evaluate_stiffness(n)

并删除out = sparse(I,J,val,n,n);这一行。

然后您将其他功能更改为：

n = 1000;
[I,J,V] = deal([]);
tic;
for i = 1:13472
    [I_i, J_i, V_i] = test_evaluate_stiffness(n);
    nE = numel(I_i);
    I(end+(1:nE)) = I_i;
    J(end+(1:nE)) = J_i;
    V(end+(1:nE)) = rand(1)*V_i;
end
K = sparse(I,J,V,n,n);
fprintf('Stiffness Calculation Complete\nTime taken = %f s\n',toc);

如果您提前知道test_evaluate_stiffness 的输出长度，您可以通过使用适当大小的zeros 矩阵预先分配数组I、J 和V 并设置他们使用类似的东西：

I((i-1)*nE + (1:nE)) = ...
J((i-1)*nE + (1:nE)) = ...
V((i-1)*nE + (1:nE)) = ...

【讨论】：

这听起来是个好主意。我会试试这个，然后报告可能的效率。每个评估的 I,J,val 的大小不一定是固定的，尽管它是有界的。因此，用零“填充”较小的 I、J、val 向量可能是值得的，这样我就可以不断使用预先分配的“全局”I、J、V 向量。我会尝试这两种方法！
@Kobye：根据我的经验，使用I(end+(1:nE)) 附加值应该执行得很好。也许你可以通过一些预分配技巧来加速它。上面的代码为您提供了大约27 = 85s/3.11s 的加速。
我只是想确认一下，在实际代码中，我能够立即将我的 K 矩阵的计算从 146 秒提高到 23 秒，这是一个巨大的进步！剩下的最大计算时间为 11 秒，是对最终 I、J、V 向量的稀疏运算，所以我认为我们已经将其简化为最基本的计算。非常感谢您的好建议。
@Kobye：很高兴听到这个消息！
@Kobye：顺便说一句：除了通常的 sparse 函数之外，还有性能更好的替代方法：blogs.mathworks.com/loren/2007/03/01/…

【解决方案2】：

剩余最大的计算，耗时 11 秒，是稀疏运算在最后的I,J,V 向量上，所以我想我们已经把它简化了骨头。

几乎...但最后一个技巧：如果您可以创建向量以便J 按升序排序，那么您将大大提高sparse 调用的速度，根据我的经验，大约是 4 倍。

（如果对I 进行排序更容易，则创建转置矩阵sparse(J,I,V) 并在之后取消转置。）

【讨论】：