矢量化这个循环

我在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)不是递归的。

我该如何解决这个问题？

来源

2015-10-30 phdstudent

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

这是非常难以矢量化，除了for循环非常快。在我的系统上，这个例子少于0.02s。除非你的实际问题更大，否则我认为这是不值得的。 – Daniel

你为什么要引导它？ – IKavanagh

在充满crazy fast parallel processors and almost-free memory latency machines，其中矢量化工具，如bsxfun和matrix multiplication可以轻松地在20000个核产卵世界末日后的世界，一个丢了魂拼命向量化这样的问题可能会冒险进入这样的事情 -

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

在所有的严重性，它不应该是内存的效率，因此，不太可能要快.. 现在，但为了演示的方式来跟踪递归和应用它的AV ectorized解决方案。

来源

2015-10-30 19:57:40 Divakar

你可以使用'parte1 = tril（rho。^ [0：n-2]'* rho。是代码中最慢的部分。 – Daniel

or'parte1 = toeplitz（rho。^ [0：n-2]，[1，zeros（1，n-2）]）;' – Daniel

@Daniel谢谢！这些可以使用，也有类似的问题来创建一个矩阵，如''这个one''中的'parte1'（http://stackoverflow.com/questions/29209361/replicate-vector-and-shift-each-copy-通过-1-行下不换循环）。这篇文章更多的是作为跟踪这些递归的指南。为了提高性能，您的解决方案中的部分矢量化方法将成为解决方案。 – Divakar

部分向量化代码将我的系统上的示例表单执行时间从0.015s降低到0.0005s。使用单个矩阵乘法我只是预先计算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

来源

2015-10-30 19:19:01 Daniel

我想你是指'se（：，i）'？ –

'se'大小为[1，n]'，两者都是正确的。 – Daniel

哦，好的。西格玛是'1x3'，我错过了。 –

除非你在不使用递归的方式重新写这个特殊的程序，你不能向量化它。

很难有涉及递归的向量化计算，因为大多数情况下，您无法为您试图解决的函数设置一个封闭的表达式（有效地在任何给定的时间点，您不需要有限的操作次数;即使你确实有一个数量有限的指数增长）。

我想你可能想考虑用filter重写你的函数。它是专门为这种关系而设计的。

编辑：

如前所述，你的描述是一个基本的过滤器。假设为@Daniel您只需使用相同的设置：

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

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

我相信你会一直觉得这是提出了三种解决方案的速度更快。对我来说，它似乎也是最自然的，因为你的问题陈述描述了一个过滤算法。

来源

2015-10-30 19:22:09

矢量化这个循环

回答

相关问题