I have the following loop in 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
I tried vectorizing it by doing:
z(1,2:end) = rho * z(1,1:end-1) + sigma * epsilon
It didn't work. I understand that the problem is that this bit: z(1,2:end) = rho * z(1,1:end-1) is not recursive.
How can I solve this?
In a post apocalyptic world filled with crazy fast parallel processors and almost-free memory latency machines, where vectorizing tools like bsxfun and matrix multiplication could easily spawn across 20,000 cores, one lost soul desperately trying to vectorize such a problem could venture into something like this -
parte1 = tril(rho.^(bsxfun(@minus,(0:n-2)',0:n-2)));
parte2 = sigma*epsilon(:,2:end);
z = [0 parte2*parte1.']
In all seriousness, it's not supposed to be memory-efficient and thus, isn't likely to be fast .. now, but is meant to demonstrate a way to trace the recursion and apply it for a vectorized solution.
Partially vectorizing the code lowered the execution time for your example form 0.015s to 0.0005s on my system. I simply pre-calculated sigma* epsilon(:,i) using a single matrix multiplication:
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
Unless you re-write this particular program in a way that does not use recursion you cannot vectorize it.
It is very hard to have vectorized calculations involving recursion because most of the time you cannot have a closed form expression for the function you try to solve (effectively at any given time-point you do not a finite number of operations; and even if you do have a finite number that will grow exponentially).
I think you might want to consider re-writing your function using filter. It is specifically designed for such relations.
EDIT:
As mentioned, what you describe is an elementary filter. Assuming the same setting as @Daniel you simply have:
...
epsilon = normrnd(0,1, 3, n);
se=sigma*epsilon;
a = [1 -rho];
b = 1;
z = filter(b, a, [0 se(2:n)] );
I believe that you will consistently find this to be the faster of the three solutions proposed. To me it also appears to be the most natural as your problem statement describes a filtering algorithm.
