Vectorize этот цикл

Question

Я следующий цикл в 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) не является рекурсивным.

Как я могу это решить?

Answer 1

В пост апокалиптическом мире, заполненном 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.'] 

Со всей серьезностью, это не должно быть эффективно использует память и, таким образом, скорее всего, не будет быстро .. Теперь, но призван продемонстрировать способ проследить рекурсию и применять его для пр екторизованное решение.

Answer 2

Частичная прорисовка кода снизила время выполнения для вашей формы от 0,015 до 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 

Answer 3

Если вы повторно не написать эту конкретную программу таким образом, чтобы не использовать рекурсию, вы не можете векторизации его.

Очень сложно иметь векторизованные вычисления с рекурсией, потому что большую часть времени у вас не может быть выражения закрытой формы для функции, которую вы пытаетесь решить (эффективно в любой момент времени вы не имеете конечного числа операций; и даже если у вас есть конечное число, которое будет экспоненциально расти).

Возможно, вы захотите рассмотреть возможность повторной записи своей функции с помощью filter. Он специально разработан для таких отношений.

EDIT:

Как уже упоминалось, что вы описываете элементарный фильтр. Предполагая, что один и тот же параметр, как @Daniel вы просто:

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

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

Я считаю, что вы будете последовательно найти, что это будет быстрее из трех предложенных решений. Для меня это также кажется наиболее естественным, так как ваш оператор проблемы описывает алгоритм фильтрации.