如何在 Matlab 中对这个 4-D 矩阵运算进行矢量化?

计算科学 matlab 矢量化
2021-12-24 14:47:39

我需要在 Matlab 中进行大量矩阵运算,但其中一个矩阵是 4-D,我对正确矢量化的正常直觉令我失望。现在我正在使用循环。也许有人可以修复我的愚蠢大脑?

for k=1:kNum
   for l=1:lNum
      L(k,l,:) = squeeze(M(k,l,:,:))\squeeze(P(k,l,:));
   end
end
1个回答

我认为不太可能对这个表达式进行矢量化。通常,向量化与对向量的赋值以及为每个项目输出单独结果的函数一起使用,例如

for i = 1:N
   A(i) = sin(i)
end

等于

i = 1:N
A(i) = sin(i)

因为 sin 对元素进行元素正弦。使用 mldivide/“反斜杠运算符”求解线性方程组时,解决方案取决于向量中的所有元素 - 除非矩阵具有特殊属性,否则不能仅按元素进行 mldivide。

为了优化这个循环,你可以考虑使用 parfor-loops(循环变量是正常的数字加一,所以它应该像做“parfor”而不是“for”一样简单)。这需要并行计算工具箱。你也可以如果索引的编号有些可预测,则展开循环以获得更好的性能。

或者,您可以尝试编写一个 .mex 文件来执行相同的操作。在 C/Fortran 中,循环与 MATLAB 向量化一样快/快。

要考虑的另一件事是重新排序内存访问:由于矩阵只是使用 [i + j*dim1 + k*dim1*dim2] 模式访问的长数组,因此确保内部循环对应于顺序存储的内存很重要. 对于二维矩阵,至少,做

for k=1:kNum
   for l=1:lNum
       L(k,l,:) = squeeze(M(k,l,:,:))\squeeze(P(k,l,:));

for l= 1:lNum
   for k=1:kNum
       L(k,l,:) = squeeze(M(k,l,:,:))\squeeze(P(k,l,:));

因为与 k-index 对应的元素是按顺序存储的。Mathworks 有一些关于它的信息,但它并不是 MATLAB 本身所特有的,任何关于 C 或 FORTRAN 数组的高级文本都应该详细介绍。这里还有一些缓存的东西。

编辑:虽然在这里没有用,因为每个解决方案的右侧和左侧都不同,如果你正在做类似的事情

for i = 1:N
    solution(:,i) = A\b(:,i)
end

你也可以

solution = A\b