何时不应该对Matlab进行向量化处理？

Question

我正在处理一些Matlab代码，用于处理大型（但不是巨大的）数据集：10,000个784元素向量（非稀疏），并计算存储在10,000 x 10稀疏矩阵中的信息。为了让代码正常工作，我使用循环处理了其中一些较棘手的部分，对这10k个项目进行了循环处理，并对稀疏矩阵中的10个项目进行了清理。

我的处理过程最初需要73次迭代（因此大约循环了730k次），并且运行时间约为120秒。这还不错，但由于这是Matlab，所以我试图将其向量化以加快速度。

最终，我有了一个完全向量化的解决方案，它得到了与初始解决方案相同的答案（因此它是正确的，或者至少与我的初始解决方案一样正确），但运行时间为274秒，几乎慢了一半！

这是我第一次遇到向量化后比迭代更慢的Matlab代码。是否有任何经验法则或最佳实践可用于确定何时可能/可能出现这种情况？

我很想分享代码以获得反馈，但它是当前开放的学校任务，所以我现在真的不能分享。如果它最终成为那些“哇，那很奇怪，你可能做错了什么”的事情之一，我可能会在一两周内重新审视它，看看我的向量化是否有误。

- Donnie

---

对你的代码进行性能分析，并将慢的部分发布在这里。只要你不会因速度而被评分，那就不会违反任何学术政策。 - user57368

Answer 1

在Matlab中，矢量化通常意味着分配更多的内存（通过托尼技巧等方式制作更大的数组以避免循环）。最近版本中循环的JIT编译得到了改进，可能意味着您矢量化解决方案所需的内存分配没有优势，但是不看代码很难说。Matlab具有出色的逐行分析器，可以帮助您查看哪些特定部分的矢量化版本需要时间。

- robince

---

是的，大的重复矩阵是向量化版本中花费所有时间的地方，我怀疑@las3rjock有正确的想法，即我的向量化解决方案可能在某些版本上更快，但在这个大小上更慢，我想我可以做他建议的绘图来检查。 - Donnie

有时候你可以通过使用bsxfun等方式来重写代码，避免使用缓慢的repmat函数。 - Amro

@Amro：你觉得bsxfun和内置循环一起工作吗？我一直想知道repmat和bsxfun之间的内存使用比较。我尝试了以下代码:A = rand(3,5e7);A = bsxfun(@minus,A,[1 2 3]');通常情况下，原地操作应该会防止分配更多的内存（就像在repmat中），但是我遇到了"内存不足"的错误！看起来我唯一（而且非常慢）可以做的方式是:for i=1:size(A,2),A(:,i)=A(:,i)-[1 2 3]';end - Jacob

@Jacob：就内存使用而言，正如@thrope所说，bsxfun必须为输出分配内存（就像MATLAB中返回结果的任何函数一样）。因此，除非您正在分配接近最大容量的数据（请参见memory函数），否则仍然最好使用bsxfun来提高速度。顺便说一句，如果您将A缩小一半，我想上面的方法也可以解决问题（在我的机器上，我甚至无法先分配A！）。但是就CPU时间而言，bsxfun要快得多...看看这个比较：http://blogs.mathworks.com/loren/2008/08/04/comparing-repmat-and-bsxfun-performance/ - Amro

我认为答案是你向量化得不好。_repmat_非常慢，我总是尽量避免使用它，特别是用于扩展标量维度。正如答案中所述，分析器是你的朋友。 - Nzbuu

Answer 2

您尝试过将执行时间作为问题规模的函数进行绘图吗（可以是每个向量的元素数量[目前为784]或向量数量[目前为10,000]）？当我对Gram-Schmidt正交化算法进行向量化时，也遇到了类似的异常情况。结果发现，在问题增长到一定大小时，向量化版本会更快，而在这一点之后，迭代版本实际上更快，如下图所示: 以下是两个实现和基准测试脚本:

clgs.m

function [Q,R] = clgs(A)
% QR factorization by unvectorized classical Gram-Schmidt orthogonalization

[m,n] = size(A);

R = zeros(n,n);     % pre-allocate upper-triangular matrix

% iterate over columns
for j = 1:n
    v = A(:,j);

    % iterate over remaining columns
    for i = 1:j-1
        R(i,j) = A(:,i)' * A(:,j);
        v = v - R(i,j) * A(:,i);
    end

    R(j,j) = norm(v);
    A(:,j) = v / norm(v);   % normalize
end
Q = A;

clgs2.m

function [Q,R] = clgs2(A)
% QR factorization by classical Gram-Schmidt orthogonalization with a
% vectorized inner loop

[m,n] = size(A);
R = zeros(n,n);     % pre-allocate upper-triangular matrix

for k=1:n
    R(1:k-1,k) = A(:,1:k-1)' * A(:,k);
    A(:,k) = A(:,k) - A(:,1:k-1) * R(1:k-1,k);
    R(k,k) = norm(A(:,k));
    A(:,k) = A(:,k) / R(k,k);
end

Q = A;

benchgs.m

n = [300,350,400,450,500];

clgs_time=zeros(length(n),1);
clgs2_time=clgs_time;

for i = 1:length(n)
    A = rand(n(i));
    tic;
    [Q,R] = clgs(A);
    clgs_time(i) = toc;

    tic;
    [Q,R] = clgs2(A);
    clgs2_time(i) = toc;
end

semilogy(n,clgs_time,'b',n,clgs2_time,'r')
xlabel 'n', ylabel 'Time [seconds]'
legend('unvectorized CGS','vectorized CGS')

- las3rjock

Answer 3

要回答“什么时候不需要将MATLAB代码向量化”的问题，更一般的说法是：

如果向量化不直接，并且使代码非常难以阅读，请不要对代码进行向量化。这是基于以下假设：

1.除您之外的其他人可能需要阅读和理解它。

2.未向量化的代码对您所需的速度已足够快。

- Hannes Ovrén

Answer 4

这不会是一个非常具体的答案，但我处理非常大的数据集（4D心脏数据集）。

有时我需要执行涉及多个4D集的操作。我可以创建一个循环或向量化的操作，它基本上在连接的5D对象上工作。 （例如，作为一个微不足道的例子，假设您想要获得平均4D对象，则可以创建一个收集行走平均值的循环，或在第5维中连接，并在其上使用平均函数）。

根据我的经验，在抛开首先创建5D对象所需的时间之后，由于执行计算时涉及的庞大大小和内存访问跳跃，通常更快地采用仍然很大但更易管理的4D对象的循环。

我要指出的另一个“微优化”技巧是matlab是“列优先顺序”。 意思是，对于我的微不足道的例子，我认为沿着第1个维度进行平均值计算比沿着第5个维度进行计算更快，因为前者涉及内存中连续的位置，而后者涉及巨大的跳跃，可以这么说。 因此，如果有意义的话，将您的超级数组存储在以您将要操作的数据为第一维的尺寸顺序中可能是值得的。

这是一个简单的例子，展示了在处理行和列时的差异：

>> A = randn(10000,10000);
>> tic; for n = 1 : 100; sum(A,1); end; toc
Elapsed time is 12.354861 seconds.
>> tic; for n = 1 : 100; sum(A,2); end; toc
Elapsed time is 22.298909 seconds.

- Tasos Papastylianou