矢量化這個循環

我在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

矢量化這個循環

回答

相關問題