Matrix Multiply> for loop> bsxfun - 奇数速度结果

Question

我有一个n乘n矩阵A,一组n个系数k(n乘1),以及一个叫做row(1乘n)的矩阵.我的目标是从第i行的第i行中减去第k行中第i个系数加权的行.我有三个问题:为什么我的for-loop表现优于内置矩阵乘法,一般来说,是什么解释了每种方法的优越性接下来,还有比我提出的三个方法更好的方法吗？

% 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循环来增加时间.

Answer 1

简短的回答:代码的更快版本是这样的:

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 = alpha x y'+ A是GER函数) .问题很可能是MATLAB不能识别它,而是运行代码,因为它理解它,即显式执行包括k*row在内的所有操作.这正是这个解决方案的问题.外部产品必须分配额外的内存大小等于矩阵的大小,这本身需要时间.考虑一下:

• 内存分配 - 因为我们不知道MATLAB如何管理内存分配,这可能意味着很多开销.
• 读矢量k,行
• 写结果矩阵(KR)
• 读取KR从A中减去它
• 阅读和写作A.

两个1000*1000矩阵是16 MB - 我怀疑你有这么多缓存.这就是为什么这个版本不是最好的原因,并且实际上可能比'内存低效'循环慢,这取决于可用的内存带宽和CPU缓存大小.

无需分配KR矩阵并将值明确存储在内存中 - 您可以在循环中计算所需的产品.因此,您有效地减少了一半的内存带宽需求并消除了内存分配开销!

• 读矢量k,行
• 读写A.

假设一个向量适合缓存(它做了 - 1000*8字节不多)你只从内存中读取k和行一次.现在,具有循环列的算法非常有意义(这可能是BLAS实现此计算的方式)

• 读取k和行并将它们保存在缓存中
• 流A具有到CPU的全内存带宽,减去k*行产品,流回内存

现在最后的效率考虑.采取我的循环结构.在每次迭代中,我们读取和写入A,并读取向量.这是每次迭代移动的16MB数据.1000次迭代总共移动了16 GB的数据.计算结果所需的两秒钟可提供8GB/s的有效内存带宽.当使用2个CPU核心时,我的系统具有16GB/s的流带宽,使用时,大约为11-12 GB/s.因此,这个顺序循环在一个核心上以60-70%的效率运行.还不错,考虑到这是一个MATLAB循环 :)

另请注意,至少在我的计算机上,基于列的循环版本比(更多一点)快了两倍,比Ak 行实现快两倍(2s vs 4.37s).这强烈表明k行确实是由MATLAB显式执行和构造的,并且总内存流量是环路版本的两倍.因此表现差了两倍.

编辑您仍然可以尝试通过直接调用相应的BLAS函数在第一个算法中执行此操作.看看这个FEX的贡献.它允许您直接从MATLAB调用BLAS和LAPACK函数.可能有用..