向量化这个循环

Posted 2023-02-16

【中文标题】向量化这个循环

【英文标题】：Vectorize this loop

【发布时间】：2015-10-30 18:39:49

【问题描述】：

我在 MATLAB 中有以下循环：

n = 20000
rho=0.9;
sigma=[0.1 0.2 0.3];
epsilon = normrnd(0,1, 3, n);
z =  NaN(1,n);
z(1,1) = 0;
for i=2:n
   z(1,i) = rho * z(1,i-1) + sigma* epsilon(:,i);
end

我尝试通过以下方式对其进行矢量化：

z(1,2:end) = rho * z(1,1:end-1) + sigma * epsilon

没有用。我知道问题在于这一点：z(1,2:end) = rho * z(1,1:end-1) 不是递归的。

我该如何解决这个问题？

【问题讨论】：

问题是每个元素都依赖于前一个元素，因此需要循环。也许bsxfun 可以解决它，但是对于递归函数，我总是使用循环。

这很难向量化，而且 for 循环非常快。在我的系统上，该示例耗时不到 0.02 秒。除非您的实际问题更大，否则我认为不值得花时间。

为什么要矢量化它？

写一个 Mex 文件。如果您不知道如何操作，那么它将是您技能的一个很好的补充。如果你知道，那么这应该不会超过半小时。

哇！很多答案的方式，所有这些都很有用！不确定接受哪个。 @Daniel，我的程序要大得多。

【参考方案1】：

在充满 crazy fast parallel processors and almost-free memory latency machines 的后世界末日世界中，bsxfun 和 matrix multiplication 等矢量化工具可以轻松跨 20,000 个内核生成，一个迷失的灵魂尝试将此类问题矢量化可能会冒险尝试这样的事情 -

parte1 = tril(rho.^(bsxfun(@minus,(0:n-2)',0:n-2)));
parte2 = sigma*epsilon(:,2:end);
z = [0 parte2*parte1.']

说真的，它不应该是内存效率高的，因此不太可能很快.. 现在，而是为了演示一种跟踪递归并应用它的方法用于矢量化解决方案。

【讨论】：

您可以使用 parte1 = tril(rho.^[0:n-2]'*rho.^(0:-1:-n+2)); 来避免 n^2 次幂计算，这是代码中最慢的部分。

或parte1 = toeplitz(rho.^[0:n-2],[1,zeros(1,n-2)]);

@Daniel 谢谢！可以使用这些，还有一个类似的问题是在this one 中创建像parte1 这样的矩阵。这篇文章更多地是作为跟踪此类递归的指南。为了提高性能，您的解决方案中的部分矢量化方法将是首选解决方案。

【参考方案2】：

对代码进行部分矢量化可将您的示例表单在我的系统上的执行时间从 0.015 秒降低到 0.0005 秒。我只是使用单个矩阵乘法预先计算了sigma* epsilon(:,i)：

n = 20000
rho=0.9;
sigma=[0.1 0.2 0.3];
epsilon = normrnd(0,1, 3, n);
se=sigma*epsilon;
z =  NaN(1,n);
z(1,1) = 0;
for i=2:n
   z(1,i) = rho * z(1,i-1) + se(i);
end

【讨论】：

我想你的意思是se(:,i)？

哦，好的。 sigma 是1x3，我错过了。

在我的机器上很有趣，原始代码需要 0.04 秒，而您的替代代码需要 0.1 秒。

【参考方案3】：

除非您以不使用递归的方式重新编写此特定程序，否则您无法对其进行矢量化。

很难进行涉及递归的向量化计算，因为大多数时候您无法为您尝试求解的函数提供封闭形式的表达式（实际上在任何给定时间点，您都没有有限数量的操作；并且即使你确实有一个有限的数字会呈指数增长）。

我认为您可能需要考虑使用filter 重写您的函数。它专为此类关系而设计。

编辑：

如前所述，您描述的是一个基本过滤器。假设您的设置与@Daniel 相同：

...
epsilon = normrnd(0,1, 3, n);
se=sigma*epsilon;

a = [1 -rho];
b = 1; 
z = filter(b, a, [0  se(2:n)] );

我相信您会始终发现这是所提出的三种解决方案中速度更快的解决方案。对我来说，这似乎也是最自然的，因为您的问题陈述描述了一种过滤算法。