我有一个n乘n的矩阵A，一组n个系数k(n乘1)，以及一个称为行(1乘n)的矩阵。我的目标是从A的第I行减去由k中的第I个系数加权的行。我有三个问题：为什么我的for循环比内置矩阵乘法表现得更好，一般来说，是什么解释了每种方法比下一种方法的优越性，有没有比我提出的三种方法更好的方法？
% Define row & coefficients used in each method
k = (1:1000).';
row = 1:1000;
% Method 1 (matrix multiply) ~15 seconds
A = magic(1e3);
tic;
for z = 1:1e3
A = A - k*row;
end
toc;
% Method 2 (for loop) ~11 seconds
B = magic(1e3);
tic;
for z = 1:1e3
for cr = 1:1000
B(cr,:) = B(cr,:) - k(cr)*row;
end
end
toc;
% method 3 (bsxfun) ~ 4 seconds
C = magic(1e3);
tic;
for z = 1:1e3
C = C - bsxfun(@times, k, row);
end
toc
isequal(A,B)
isequal(A,C)
注意，我在一个算法中做这些行减法。我简化了一些代码，创建了这个玩具测试用例，但计算的关键仍然存在。此外，为了避免混淆，使用带有z的for循环来增加时间。
简短的回答：代码的更快版本是：
tic;
for z = 1:1e3
for cr = 1:1000
B(:,cr) = B(:,cr) - k*row(cr);
end
end
toc;
你可能想看看我之前对这个问题的回答。简而言之，循环是基于行的，而MATLAB是基于列的。这意味着列在内存中是连续的。原始循环在行上迭代，这是低效的。
我电脑上的执行时间：
% A - k*row
Elapsed time is 4.370238 seconds.
% B(cr,:) = B(cr,:) - k(cr)*row;
Elapsed time is 9.537338 seconds.
% C = C - bsxfun(@times, k, row);
Elapsed time is 3.039836 seconds.
B(:,cr) = B(:,cr) - k*row(cr);
Elapsed time is 2.028186 seconds.
解释。您的第一个版本不是矩阵乘法，而是两个向量的外积，从而产生大小为1000 x 1000的矩阵。空穴计算是BLAS2秩1更新(a=alphaxy’+a是GER函数)。最有可能的问题是MATLAB没有识别它，而是按照它所理解的方式运行代码，即显式执行包括k*行在内的所有操作。这正是这个解决方案的问题所在。外积必须分配与矩阵大小相等的额外内存，这本身需要时间。考虑一下：
- 内存分配-由于我们不知道MATLAB是如何管理内存分配的，这可能意味着大量开销
- 读取矢量k，行
- 结果矩阵(KR)的写入
- 读取KR以将其从A中减去
- A的读写
两个1000*1000矩阵是16 MB-我怀疑您是否有这么多缓存。这就是为什么这个版本不是最好的，而且实际上可能比"内存低效"循环慢，这取决于可用的内存带宽和CPU缓存大小。
不需要分配KR矩阵并将值显式存储在内存中——您可以在循环中计算所需的乘积。因此，您可以有效地将内存带宽需求减半，并消除内存分配开销！
- 读取矢量k，行
- 读写A
假设一个向量适合缓存(它确实这样做了-1000*8字节并不多)，则只从内存中读取k和行一次。现在，在列上循环的算法是完全合理的(这可能是BLAS实现这种计算的方式)
- 读取k和row并将它们保存在缓存中
- 流A具有到CPU的全部内存带宽，减去k*行乘积，流回到内存
现在是最后的效率考虑因素。以我的循环结构为例。在每次迭代中，我们读写A，并读取向量。这是每次迭代移动的16MB数据。1000次迭代总共移动了16GB的数据。计算结果所需的两秒钟可提供8GB/s的有效内存带宽。当使用两个CPU核时，我的系统具有16GB/s的流带宽，当使用一个时，大约为11-12GB/s。因此，这种顺序循环在一个核心上以60-70%的效率运行。不错，考虑到这是一个MATLAB循环：)
还需要注意的是，至少在我的计算机上，基于列的循环版本比a-k行实现(2s比4.37s)快(多一点)。这有力地表明，k行确实是由MATLAB显式执行和构建的，并且总内存流量是循环版本的两倍。因此性能差了两倍。
编辑您仍然可以通过直接调用相应的BLAS函数来尝试像第一个算法中那样进行编辑。看看这个FEX贡献。它允许您直接从MATLAB中调用BLAS和LAPACK函数。可能有用。。