matlab - 提高许多子矩阵左除运算的性能(mldivide，\)

Question

我有两个矩阵，a1 和 a2。 a1 是 3x12000，a2 是 3x4000。我想创建另一个 3x4000 数组，它是 a1 的 3x3 子矩阵的左矩阵除法 (mldivide, \) 和a2 的 3x1 子矩阵。您可以使用 for 循环轻松完成此操作:

for ii = 1:3:12000
  a = a1(:,ii:ii+2)\a2(:, ceil(ii/3));
end

但是，我想知道是否有更快的方法来做到这一点。

编辑:我知道预分配会提高速度，我只是出于视觉目的展示它。

Edit2:删除了数组的迭代增加。看来我的问题被误解了一点。我主要想知道是否有一些矩阵运算可以实现我的目标，因为这可能比 for 循环更快，即将 a1 reshape 为 3x3x4000 矩阵和 a2到 3x1x4000 矩阵和左矩阵一次性划分每个级别，但是，您不能将左矩阵与 3D 矩阵划分。

Answer 1

通过将 a1 的子矩阵放在 12000x12000 矩阵的对角线上，您可以创建一个包含许多独立“子系统”的方程组，如下所示:

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 寻找求解线性方程组的最佳方法，因此如果您求解每个子系统独立地，将为每个子系统选择最佳方法，但如果您将主题作为一个系统来解决，Matlab 将为所有子系统一起找到最佳方法 - 而不是为每个子系统找到最佳方法。

注意:如图Stefano M的 answer ，使用一个大方程组(使用 kron 和稀疏矩阵)比使用 for 循环(带预分配)更快，仅适用于非常小的方程子系统(在我的 PC 上，用于数字of equation <= 7) 对于更大规模的方程子系统，使用 for 循环更快。

比较不同的方法

我编写并运行了一个代码来比较 4 种不同的方法来解决这个问题:

1. for 循环，无预分配
2. for 循环，带预分配
3. 克朗
4. 细胞乐趣

测试:

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