每列的 MATLAB 矩阵乘法与 for 循环答案

【问题标题】：MATLAB matrix multiplication vs for loop for each column每列的 MATLAB 矩阵乘法与 for 循环
【发布时间】：2011-02-06 20:52:41
【问题描述】：

当两个矩阵相乘时，我尝试了以下两个选项：

res = X*A;

for i = 1:size(A,2)
    res(:,i) = X*A(:,i);
end

我在两者中都为 res 预先分配了内存。令人惊讶的是，我发现选项 2 更快。

谁能解释一下这是怎么回事？

编辑：我试过了

K=10000;
clear t1 t2
t1=zeros(K,1);
t2=zeros(K,1);

for k=1:K
    clear res
    x = rand(100,100);
    a = rand(100,100);
    tic
    res = x*a;
    t1(k) = toc;
end

for k=1:K
    clear res2
    res2 = zeros(100,100);
    x = rand(100,100);
    a = rand(100,100);
    tic
    for i = 1:100
        res2(:,i) = x*a(:,i);
    end
    t2(k) = toc;
end

【问题讨论】：

在 OSX 2010b 上，我得到的版本一的中位时间为 0.0001，版本二的中位时间为 0.0008。换句话说，矢量化版本快了大约 7 倍。
为了消除一些可变性，最好在开始第二次运行之前重置随机生成器。无论如何，我还发现第一个选项要快几倍。

标签： matlab matrix-multiplication

【解决方案1】：

我在循环中运行这两个代码 1000 次。平均而言（但不总是），第一个矢量化代码的速度要快 3-4 倍。我在启动计时器之前清除了结果变量并进行了预分配。

x = rand(100,100);
a = rand(100,100);

K=1000;
clear t1 t2
t1=zeros(K,1);
t2=zeros(K,1);

for k=1:K
    clear res
    tic
    res = x*a;
    t1(k) = toc;
end

for k=1:K
    clear res2
    res2 = zeros(100,100);
    tic
    for i = 1:100
        res2(:,i) = x*a(:,i);
    end
    t2(k) = toc;
end

所以，永远不要根据单次运行得出时间结论。

【讨论】：

我修改了您的代码并编辑了问题。在我的机器上，如果我使用不同的 x 和 a，在每次迭代中，第二个版本仍然更快。

【解决方案2】：

我相信我可以了解这两种方法在时间上的差异，以及人们获得不同相对速度的原因。

在 Matlab 版本 2008a（或该版本附近的版本）之前，for 循环在任何 Matlab 代码中都会受到重大影响，因为解释器（非常易读的脚本和代码的较低级别实现之间的一层）必须重新- 每次通过for循环解释代码。

自该版本发布以来，解释器逐渐变得更好，因此，在运行现代版本的 Matlab 时，解释器可以查看您的代码并说“啊哈！我知道他在做什么，让我稍微优化一下" 并通过重新解释代码来避免可能受到的影响。

我希望执行矩阵乘法的两种方法能够在相同的时间内进行评估，为什么 for 循环实现运行得更快是因为解释器优化中的一些细节，我们普通人是不知道的。

我们应该从中吸取的一个广泛教训是，并非所有版本都是平等的。我确实使用两个 Matlab 附加组件 SimBiology 和并行计算工具箱处理了几个前沿案例，这两者（特别是如果您希望它们一起工作）在执行速度方面取决于版本，并且不时其他稳定性问题。因此，我保留了三个最新版本的 Matlab，将测试我从每个版本中得到相同的答案，如果我发现某些功能存在问题，我偶尔会回滚到早期版本。对于大多数人来说，这可能有点过头了，但可以让您了解版本差异。

希望这会有所帮助。

编辑：

为了澄清，代码向量化仍然很重要。但是给定一个像这样的脚本：

x_slow = zeros(1,1e5);
x_fast = zeros(1,1e5);


tic;
for i=1:1e5
    x_slow(i) = log(i);
end
time_slow = toc; % evaluates for me in .0132 seconds

tic;
x_fast = log(1:1e5);
time_fast = toc; % evaluates for me in .0055 seconds

在过去的几个版本中，基于解释器的改进，time_slow 和 time_fast 之间的差异已经减少。我看到的例子我相信是在 2000a 与 2008b 上的，但这取决于我的回忆。

Oli 和 Yuk 可能还解决了其他一些问题。 time_1 和 time_2 之间经常存在差异：

tic; x = log(1:1e5); time_1 = toc
tic; x = log(1:1e5); time_2 = toc

因此，100 万次评估与一次评估的测试是有价值的，这取决于内存 x 的位置（在缓存中或没有）。

希望这对您有帮助。

【讨论】：

真的吗？每次迭代都会重新解释源代码？我觉得很难相信！你有引用吗？
参考，恐怕不行。这是我在两年前的一次 HPC 讲座中回忆的。我已经编辑了帖子以澄清我的意思。

【解决方案3】：

这很可能是缓存的影响。 a 在您执行第二个版本时已经在缓存中，因此它具有优势。尝试创建一组独立的输入以使其公平。此外，测量时间可能会更好。对此进行 100 万次迭代，以消除由于外部影响导致的典型变化。

【讨论】：

【解决方案4】：

在我看来，您没有正确地乘以矩阵，您需要将 X 矩阵的第 i 行和 A 矩阵的第 j 列的所有乘积相加，这可能是一个原因。查看here 看看它是如何完成的。

【讨论】：