提高许多子矩阵左除运算的性能（mldivide，\）
Improving the performance of many sub-matrix left division operations (mldivide, \)
我有两个矩阵，a1 和 a2。 a1 是 3x12000，a2 是 3x4000。我想创建另一个 3x4000 数组，它是 a1 的 3x3 子矩阵和 [=12 的 3x1 子矩阵的左矩阵除法（mldivide，\） =].您可以使用 for 循环轻松完成此操作：
for ii = 1:3:12000
a = a1(:,ii:ii+2)\a2(:, ceil(ii/3));
end
但是，我想知道是否有更快的方法来做到这一点。
编辑：我知道预分配会提高速度，我只是出于视觉目的展示它。
Edit2：删除了数组的迭代增加。看来我的问题被误解了一点。我主要想知道是否有一些矩阵运算可以实现我的目标，因为这可能比 for 循环更快，即重塑 a1 到 3x3x4000 矩阵和 a2 到 3x1x4000 矩阵并离开matrix divide 每一级一次，然而，你不能用 3D 矩阵离开 matrix divide.
Marginal最大的改进是预分配输出矩阵，而不是增长它：
A1 = reshape(A1, 3, 3, []);
a = zeros(size(A2));
for i = 1:size(A1,3)
a(:,i) = A1(:,:,i) \ A2(:,i);
end
有了预分配数组，如果Parallel Toolbox可用，你可以试试parfor
编辑
这个答案不再相关，因为 OP 改写了问题以避免增加结果数组，这是最初的主要瓶颈。
这里的问题是要解决4000个独立的3x3线性系统。矩阵是如此之小，以至于 ad hoc 解决方案可能很有趣，特别是如果有人有一些关于矩阵属性的信息（对称或不对称，条件数等）。然而，坚持使用 \ matlab 运算符，加速计算的最佳方法是明确利用并行性，例如通过 parfor 命令。
Eliahu Aaron的另一个的稀疏矩阵解法确实很聪明，但是它的速度优势不是一般的，而是要看具体问题的大小。
使用此功能，您可以探索不同规模的问题：
function [t2, t3] = sotest(n, n2)
a1 = rand(n,n*n2);
a2 = rand(n,n2);
tic
a1_reshape = reshape(a1, n, n, []);
a_method2 = zeros(size(a2));
for i = 1:size(a1_reshape,3)
a_method2(:,i) = a1_reshape(:,:,i) \ a2(:,i);
end
t2 = toc;
tic
a1_kron = kron(speye(n2),ones(n));
a1_kron(logical(a1_kron)) = a1(:);
a_method3 = a1_kron\a2(:);
a_method3 = reshape(a_method3, [n n2]);
t3 = toc;
assert ( norm(a_method2 - a_method3, 1) / norm(a_method2, 1) < 1e-8)
确实对于 n=3 稀疏矩阵方法明显优越，但是对于增加 n 它变得不那么有竞争力了
上图是用
得到的
>> for i=1:20; [t(i,1), t(i,2)] = sotest(i, 50000); end
>> loglog(1:20, t, '*-')
我最后的评论是，使用密集 \ 运算符的显式循环确实很快；稀疏矩阵公式的准确性稍差，在边缘情况下可能会出现问题；并且可以肯定的是，稀疏矩阵解决方案的可读性不是很好。如果要解决的系统数量 n2 非常大 (>1e6)，那么也许应该探索 ad hoc 解决方案。
您正在求解一系列 N 个系统，每个系统有 m 个线性方程，N 个系统的形式为
Ax = b
您可以将它们转换为单一的 Nm 线性方程组：
|A1 0 0 ... 0 | |x1| |b1|
|0 A2 0 ... 0 | |x2| |b2|
|0 0 A3 ... 0 | |x3| = |b3|
|. . . ... . | |. | |. |
|0 0 0 ... AN| |xN| |bN|
然而，求解这一方程组比求解所有小方程组要昂贵得多。通常，成本是 O(n^3)，所以你从 O(N m^3) 到 O((Nm)^3)。巨大的悲观情绪。 （，显然可以利用矩阵的稀疏性。）
降低计算成本是可以的，但是你需要对数据提供保证。例如，如果矩阵 A 是正定的，则可以更便宜地求解系统。尽管如此，考虑到你正在处理 3x3 矩阵，那里的收益将很小，因为这些是非常简单的系统来解决。
如果您问这个问题是因为您认为 MATLAB 中的循环天生就很慢，那么您应该知道情况已不再如此，自从 MATLAB 大约 15 年前获得 JIT 以来就不再是这种情况了。今天，许多矢量化尝试导致同样快速的代码，并且通常导致更慢的代码，尤其是对于大数据。 （如果有必要，我可以找出我在 SO 上发布的一些时间来证明这一点。）
我认为一次解决所有系统可以减少每次调用运算符 \ 时 MATLAB 执行的检查次数。也就是说，hard-coding 问题的大小和类型可能会自始至终得到改善。但这样做的唯一方法是写一个 MEX-file.
您可以通过将 a1 的 sub-matrices 放在 12000x12000 矩阵的对角线上来创建一个包含许多独立 'sub-systems' 方程的方程组，如下所示：
a1(1,1) a1(1,2) a1(1,3) 0 0 0 0 0 0
a1(2,1) a1(2,2) a1(2,3) 0 0 0 0 0 0
a1(3,1) a1(3,2) a1(3,3) 0 0 0 0 0 0
0 0 0 a1(1,4) a1(1,5) a1(1,6) 0 0 0
0 0 0 a1(2,4) a1(2,5) a1(2,6) 0 0 0
0 0 0 a1(3,4) a1(3,5) a1(3,6) 0 0 0
0 0 0 0 0 0 a1(1,7) a1(1,8) a1(1,9)
0 0 0 0 0 0 a1(2,7) a1(2,8) a1(2,9)
0 0 0 0 0 0 a1(3,7) a1(3,8) a1(3,9)
然后左除以 a2(:)。
这可以使用 kron 和像这样的稀疏矩阵 (source) 来完成：
a1_kron = kron(speye(12000/3),ones(3));
a1_kron(logical(a1_kron)) = a1(:);
a = a1_kron\a2(:);
a = reshape(a, [3 12000/3]);
优点 - 速度：这比我 PC 上带有预分配的 for 循环快 3-4 倍。
缺点：使用这种方法必须考虑一个缺点：使用左除法时，Matlab 寻找求解线性方程组的最佳方法，因此如果您独立解决每个 sub-system，将为每个 sub-system 选择最佳方法，但如果您将主题作为一个系统解决，Matlab 将一起为所有 sub-system 找到最佳方法 - 不是每个 sub-system.
的最佳选择
注意：如Stefano M's 所示，使用一个大方程组（使用kron和稀疏矩阵）比使用for循环更快（使用预分配）仅适用于非常小的 sub-systems 方程（在我的 PC 上，方程数 <= 7）对于更大尺寸的 sub-systems 方程，使用 for 循环更快.
比较不同的方法
我写了 运行 一个代码来比较 4 种不同的方法来解决这个问题：
- for 循环，无预分配
- for 循环，带预分配
kroncellfun
测试：
n = 1200000;
a1 = rand(3,n);
a2 = rand(3,n/3);
disp('Method 1: for loop, no preallocation')
tic
a_method1 = [];
for ii = 1:3:n
a_method1 = [a_method1 a1(:,ii:ii+2)\a2(:, ceil(ii/3))];
end
toc
disp(' ')
disp('Method 2: for loop, with preallocation')
tic
a1_reshape = reshape(a1, 3, 3, []);
a_method2 = zeros(size(a2));
for i = 1:size(a1_reshape,3)
a_method2(:,i) = a1_reshape(:,:,i) \ a2(:,i);
end
toc
disp(' ')
disp('Method 3: kron')
tic
a1_kron = kron(speye(n/3),ones(3));
a1_kron(logical(a1_kron)) = a1(:);
a_method3 = a1_kron\a2(:);
a_method3 = reshape(a_method3, [3 n/3]);
toc
disp(' ')
disp('Method 4: cellfun')
tic
a1_cells = mat2cell(a1, size(a1, 1), repmat(3 ,1,size(a1, 2)/3));
a2_cells = mat2cell(a2, size(a2, 1), ones(1,size(a2, 2)));
a_cells = cellfun(@(x, y) x\y, a1_cells, a2_cells, 'UniformOutput', 0);
a_method4 = cell2mat(a_cells);
toc
disp(' ')
结果：
Method 1: for loop, no preallocation
Elapsed time is 747.635280 seconds.
Method 2: for loop, with preallocation
Elapsed time is 1.426560 seconds.
Method 3: kron
Elapsed time is 0.357458 seconds.
Method 4: cellfun
Elapsed time is 3.390576 seconds.
比较四种方法的结果，你可以看到使用方法3 - kron给出的结果略有不同：
disp(['sumabs(a_method1(:) - a_method2(:)): ' num2str(sumabs(a_method1(:)-a_method2(:)))])
disp(['sumabs(a_method1(:) - a_method3(:)): ' num2str(sumabs(a_method1(:)-a_method3(:)))])
disp(['sumabs(a_method1(:) - a_method4(:)): ' num2str(sumabs(a_method1(:)-a_method4(:)))])
结果：
sumabs(a_method1(:) - a_method2(:)): 0
sumabs(a_method1(:) - a_method3(:)): 8.9793e-05
sumabs(a_method1(:) - a_method4(:)): 0
我有两个矩阵，a1 和 a2。 a1 是 3x12000，a2 是 3x4000。我想创建另一个 3x4000 数组，它是 a1 的 3x3 子矩阵和 [=12 的 3x1 子矩阵的左矩阵除法（mldivide，\） =].您可以使用 for 循环轻松完成此操作：
for ii = 1:3:12000
a = a1(:,ii:ii+2)\a2(:, ceil(ii/3));
end
但是，我想知道是否有更快的方法来做到这一点。
编辑：我知道预分配会提高速度，我只是出于视觉目的展示它。
Edit2：删除了数组的迭代增加。看来我的问题被误解了一点。我主要想知道是否有一些矩阵运算可以实现我的目标，因为这可能比 for 循环更快，即重塑 a1 到 3x3x4000 矩阵和 a2 到 3x1x4000 矩阵并离开matrix divide 每一级一次，然而，你不能用 3D 矩阵离开 matrix divide.
Marginal最大的改进是预分配输出矩阵，而不是增长它：
A1 = reshape(A1, 3, 3, []);
a = zeros(size(A2));
for i = 1:size(A1,3)
a(:,i) = A1(:,:,i) \ A2(:,i);
end
有了预分配数组，如果Parallel Toolbox可用，你可以试试parfor
编辑
这个答案不再相关，因为 OP 改写了问题以避免增加结果数组，这是最初的主要瓶颈。
这里的问题是要解决4000个独立的3x3线性系统。矩阵是如此之小，以至于 ad hoc 解决方案可能很有趣，特别是如果有人有一些关于矩阵属性的信息（对称或不对称，条件数等）。然而，坚持使用 \ matlab 运算符，加速计算的最佳方法是明确利用并行性，例如通过 parfor 命令。
Eliahu Aaron的另一个
使用此功能，您可以探索不同规模的问题：
function [t2, t3] = sotest(n, n2)
a1 = rand(n,n*n2);
a2 = rand(n,n2);
tic
a1_reshape = reshape(a1, n, n, []);
a_method2 = zeros(size(a2));
for i = 1:size(a1_reshape,3)
a_method2(:,i) = a1_reshape(:,:,i) \ a2(:,i);
end
t2 = toc;
tic
a1_kron = kron(speye(n2),ones(n));
a1_kron(logical(a1_kron)) = a1(:);
a_method3 = a1_kron\a2(:);
a_method3 = reshape(a_method3, [n n2]);
t3 = toc;
assert ( norm(a_method2 - a_method3, 1) / norm(a_method2, 1) < 1e-8)
确实对于 n=3 稀疏矩阵方法明显优越，但是对于增加 n 它变得不那么有竞争力了
上图是用
得到的>> for i=1:20; [t(i,1), t(i,2)] = sotest(i, 50000); end
>> loglog(1:20, t, '*-')
我最后的评论是，使用密集 \ 运算符的显式循环确实很快；稀疏矩阵公式的准确性稍差，在边缘情况下可能会出现问题；并且可以肯定的是，稀疏矩阵解决方案的可读性不是很好。如果要解决的系统数量 n2 非常大 (>1e6)，那么也许应该探索 ad hoc 解决方案。
您正在求解一系列 N 个系统，每个系统有 m 个线性方程，N 个系统的形式为
Ax = b
您可以将它们转换为单一的 Nm 线性方程组：
|A1 0 0 ... 0 | |x1| |b1|
|0 A2 0 ... 0 | |x2| |b2|
|0 0 A3 ... 0 | |x3| = |b3|
|. . . ... . | |. | |. |
|0 0 0 ... AN| |xN| |bN|
然而，求解这一方程组比求解所有小方程组要昂贵得多。通常，成本是 O(n^3)，所以你从 O(N m^3) 到 O((Nm)^3)。巨大的悲观情绪。 （
降低计算成本是可以的，但是你需要对数据提供保证。例如，如果矩阵 A 是正定的，则可以更便宜地求解系统。尽管如此，考虑到你正在处理 3x3 矩阵，那里的收益将很小，因为这些是非常简单的系统来解决。
如果您问这个问题是因为您认为 MATLAB 中的循环天生就很慢，那么您应该知道情况已不再如此，自从 MATLAB 大约 15 年前获得 JIT 以来就不再是这种情况了。今天，许多矢量化尝试导致同样快速的代码，并且通常导致更慢的代码，尤其是对于大数据。 （如果有必要，我可以找出我在 SO 上发布的一些时间来证明这一点。）
我认为一次解决所有系统可以减少每次调用运算符 \ 时 MATLAB 执行的检查次数。也就是说，hard-coding 问题的大小和类型可能会自始至终得到改善。但这样做的唯一方法是写一个 MEX-file.
您可以通过将 a1 的 sub-matrices 放在 12000x12000 矩阵的对角线上来创建一个包含许多独立 'sub-systems' 方程的方程组，如下所示：
a1(1,1) a1(1,2) a1(1,3) 0 0 0 0 0 0
a1(2,1) a1(2,2) a1(2,3) 0 0 0 0 0 0
a1(3,1) a1(3,2) a1(3,3) 0 0 0 0 0 0
0 0 0 a1(1,4) a1(1,5) a1(1,6) 0 0 0
0 0 0 a1(2,4) a1(2,5) a1(2,6) 0 0 0
0 0 0 a1(3,4) a1(3,5) a1(3,6) 0 0 0
0 0 0 0 0 0 a1(1,7) a1(1,8) a1(1,9)
0 0 0 0 0 0 a1(2,7) a1(2,8) a1(2,9)
0 0 0 0 0 0 a1(3,7) a1(3,8) a1(3,9)
然后左除以 a2(:)。
这可以使用 kron 和像这样的稀疏矩阵 (source) 来完成：
a1_kron = kron(speye(12000/3),ones(3));
a1_kron(logical(a1_kron)) = a1(:);
a = a1_kron\a2(:);
a = reshape(a, [3 12000/3]);
优点 - 速度：这比我 PC 上带有预分配的 for 循环快 3-4 倍。
缺点：使用这种方法必须考虑一个缺点：使用左除法时，Matlab 寻找求解线性方程组的最佳方法，因此如果您独立解决每个 sub-system，将为每个 sub-system 选择最佳方法，但如果您将主题作为一个系统解决，Matlab 将一起为所有 sub-system 找到最佳方法 - 不是每个 sub-system.
的最佳选择注意：如Stefano M's kron和稀疏矩阵）比使用for循环更快（使用预分配）仅适用于非常小的 sub-systems 方程（在我的 PC 上，方程数 <= 7）对于更大尺寸的 sub-systems 方程，使用 for 循环更快.
比较不同的方法
我写了 运行 一个代码来比较 4 种不同的方法来解决这个问题：
- for 循环，无预分配
- for 循环，带预分配
kroncellfun
测试：
n = 1200000;
a1 = rand(3,n);
a2 = rand(3,n/3);
disp('Method 1: for loop, no preallocation')
tic
a_method1 = [];
for ii = 1:3:n
a_method1 = [a_method1 a1(:,ii:ii+2)\a2(:, ceil(ii/3))];
end
toc
disp(' ')
disp('Method 2: for loop, with preallocation')
tic
a1_reshape = reshape(a1, 3, 3, []);
a_method2 = zeros(size(a2));
for i = 1:size(a1_reshape,3)
a_method2(:,i) = a1_reshape(:,:,i) \ a2(:,i);
end
toc
disp(' ')
disp('Method 3: kron')
tic
a1_kron = kron(speye(n/3),ones(3));
a1_kron(logical(a1_kron)) = a1(:);
a_method3 = a1_kron\a2(:);
a_method3 = reshape(a_method3, [3 n/3]);
toc
disp(' ')
disp('Method 4: cellfun')
tic
a1_cells = mat2cell(a1, size(a1, 1), repmat(3 ,1,size(a1, 2)/3));
a2_cells = mat2cell(a2, size(a2, 1), ones(1,size(a2, 2)));
a_cells = cellfun(@(x, y) x\y, a1_cells, a2_cells, 'UniformOutput', 0);
a_method4 = cell2mat(a_cells);
toc
disp(' ')
结果：
Method 1: for loop, no preallocation
Elapsed time is 747.635280 seconds.
Method 2: for loop, with preallocation
Elapsed time is 1.426560 seconds.
Method 3: kron
Elapsed time is 0.357458 seconds.
Method 4: cellfun
Elapsed time is 3.390576 seconds.
比较四种方法的结果，你可以看到使用方法3 - kron给出的结果略有不同：
disp(['sumabs(a_method1(:) - a_method2(:)): ' num2str(sumabs(a_method1(:)-a_method2(:)))])
disp(['sumabs(a_method1(:) - a_method3(:)): ' num2str(sumabs(a_method1(:)-a_method3(:)))])
disp(['sumabs(a_method1(:) - a_method4(:)): ' num2str(sumabs(a_method1(:)-a_method4(:)))])
结果：
sumabs(a_method1(:) - a_method2(:)): 0
sumabs(a_method1(:) - a_method3(:)): 8.9793e-05
sumabs(a_method1(:) - a_method4(:)): 0