使用矢量化求解多个线性系统答案

【问题标题】：Solving multiple linear systems using vectorization使用矢量化求解多个线性系统
【发布时间】：2011-06-14 14:05:19
【问题描述】：

对不起，如果这很明显，但我搜索了一段时间并没有找到任何东西（或错过了）。

我正在尝试解决 Ax=B 形式的线性系统，其中 A 是 4x4 矩阵，B 是 4x1 向量。

我知道对于单个系统我可以使用mldivide 来获取x：x=A\B。

但是，我正在尝试解决大量系统（可能 > 10000）并且我不愿意使用 for 循环，因为我被告知在许多 MATLAB 问题中它比矩阵公式要慢得多。

然后我的问题是：有没有办法使用 A 4x4x N 和 BAx=B /em> 一个矩阵 4x N ?

PS：我不知道它是否重要，但 B 向量对于所有系统都是相同的。

【问题讨论】：

使用矩阵在 MATLAB 中提高性能的整个概念称为 vectorization，它并不像“如果它是矩阵形式 -> 它总是更快”那么简单从答案中可以看出。
@Jacob：是的，我的意思是矢量化，感谢您的编辑。我知道有时您无法避免循环，但这只是从我学习 Matlab 的方式来看，并考虑到所有面向向量化的本机函数，我想知道该站点的有经验的人是否已经遇到过这个问题.感谢@Tom 和@Amro 的回答，感谢@Rasman 和@eat 的替代方案。
我感觉这是一个常见的场景（Matlab 开发人员可以改进）。我正在运行类似的 for 循环数百万次。

标签： matlab linear-algebra

【解决方案1】：

您应该使用 for 循环。如果A 保持不变而B 发生变化，则预先计算分解并重用它可能会有好处。但是对于 A 发生变化而 B 保持不变的问题，除了求解 N 个线性系统之外别无选择。

您也不应该太担心循环的性能成本：MATLAB JIT compiler 意味着循环在最新版本的 MATLAB 上通常可以同样快。

【讨论】：

【解决方案2】：

我认为您无法进一步优化它。正如@Tom 所解释的，由于A 是一个变化的，所以事先考虑各种A 没有任何好处......

考虑到您提到的尺寸，除了循环解决方案非常快：

A = rand(4,4,10000);
B = rand(4,1);          %# same for all linear systems

tic
X = zeros(4,size(A,3));
for i=1:size(A,3)
    X(:,i) = A(:,:,i)\B;
end
toc

经过的时间是 0.168101 秒。

【讨论】：

0.091238 秒...比我预期的要快得多:)。 A 的小维度很好，因为它显着降低了问题的复杂性

【解决方案3】：

问题来了：您正在尝试对 3d 矩阵执行 2D 操作 (mldivide)。无论您如何看待它，您都需要按索引引用矩阵，这是时间惩罚开始的地方......问题不是 for 循环，而是人们如何使用它们。如果您可以以不同的方式构建问题，那么也许您可以找到更好的选择，但现在您有几个选择：

1 - 墨西哥

2 - 并行处理（写一个parfor循环）

3 - CUDA

【讨论】：

@Jacob：不，parfor 在这里实际上会更慢。根据 Amro 的时间（以及评论中的 Rasman），每个任务花费的平均时间约为 10^-6 秒。调度程序与工作人员通信所花费的时间比这要慢，或者可能是相同的数量级。所以考虑到这一点，它会更慢。在我的机器上，不带parfor 总共需要 ~0.1 秒，带 parfor 和 8 个工作人员总共需要 ~0.5 秒。
@yoda 如果除法是唯一的操作，那么你是对的，否则，它可能有用也可能没用
@yoda 和 @Jacob：差异很小，但是可以节省打开池的时间，并且我的 parfor 运行系统的速度比 for 循环快一点（当您增加迭代次数，例如 100000)
@yoda：parfor 开销在很大程度上取决于处理器。较新的 Intel 多核显着降低了通信开销。

【解决方案4】：

这是一个相当深奥的解决方案，它利用了 MATLAB 的特殊优化。构建一个巨大的 4k x 4k sparse 矩阵，其中 4x4 块在对角线上。然后同时解决所有问题。

在我的机器上，这获得了与 @Amro/Tom 的 for-loop 解决方案相同的单精度解决方案，但速度更快。

n = size(A,1);
k = size(A,3);
AS = reshape(permute(A,[1 3 2]),n*k,n);
S = sparse( ...
  repmat(1:n*k,n,1)', ...
  bsxfun(@plus,reshape(repmat(1:n:n*k,n,1),[],1),0:n-1), ...
  AS, ...
  n*k,n*k);
X = reshape(S\repmat(B,k,1),n,k);

随便举个例子：

For k = 10000
For loop: 0.122570 seconds.
Giant sparse system: 0.032287 seconds.

如果您知道您的 4x4 矩阵是正定矩阵，那么您可以在 S 上使用 chol 来提高准确性。

这愚蠢。但是，即使使用 JIT，matlab 的 for 循环 still 在 2015 年也是如此缓慢。当 k 不太大时，此解决方案似乎找到了一个最佳点，因此所有内容仍适合内存。

【讨论】：

我很好奇——如果B 的大小从B = rand(4,1) 更改为B = rand(4,4)，同样的策略会起作用吗？对于Nx=32，我的测试数据参考A = rand(Nx+1,Nx+1,Nx); % size (33,33,32) 和B = rand(Nx+1,Nx) % size (33,32)。 X = reshape(S\repmat(B,k,1),n,k); 行不再有效，因为 S\repmat(B,k,1) 的大小为 1056 x 32 而不是 1056 x 1。
如果将 B 的列放在求解的 rhs 列中，它应该可以工作。
您能否澄清一下您所说的解决方案的 RHS 是什么意思？
我认为我不了解您的矩阵的维度。在原始帖子中，A 是 4x4xN，因此 A(:,:,1) 是 4x4 矩阵。 B 是 4xN，因此 B(:,1) 是 4x1。这样 A(:,:,1) \ B(:,1) 才有意义。如果您说 B 现在应该是 4x4xN ，那么 B(:,:,1) 是 4x4 并且 A(:,:,1) \ B(:,:,1) 仍然有意义。然后类似的解决方案是沿着 S 的对角线撒上 As 并将 Bs 重塑为垂直堆栈，这样您就可以做到 Xs = As \ Bs 并且 Xs 将是 4*N x 4 的解决方案堆栈。

【解决方案5】：

我知道这篇文章已经有好几年了，但无论如何我都会贡献我的两分钱。您可以将所有 A 矩阵放入一个更大的块对角矩阵中，其中一个大矩阵的对角线上将有 4x4 块。右手边将是你所有的 b 向量一遍又一遍地堆叠在一起。一旦你设置了它，它就被表示为一个稀疏系统，并且可以使用 mldivide 选择的算法有效地解决。这些块在数字上是解耦的，所以即使那里有奇异块，当您使用 mldivide 时，非奇异块的答案应该是正确的。在 MATLAB Central 上有一段代码采用了这种方法：

http://www.mathworks.com/matlabcentral/fileexchange/24260-multiple-same-size-linear-solver

我建议尝试看看该方法是否比循环更快。我怀疑它可能更有效，尤其是对于大量小型系统。特别是，如果在 N 个矩阵中 A 的系数有很好的公式，您可以使用 MATLAB 向量运算（无循环）构建完整的左侧，这可以为您节省额外的成本。正如其他人所指出的，矢量化操作并不总是更快，但根据我的经验，它们通常是这样。

【讨论】：