gpt4 book ai didi

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

转载 作者:行者123 更新时间:2023-12-02 16:55:06 27 4
gpt4 key购买 nike

我有两个矩阵,a1a2a1 是 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 矩阵划分。

最佳答案

通过将 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 Manswer ,使用一个大方程组(使用 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

关于matlab - 提高许多子矩阵左除运算的性能(mldivide,\),我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/56825113/

27 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com