Vectorization recursive vector operation

Hi! I need some help trying to vectorize a code. I have a loop that creates a vector, using the last element of the same matrix. An example code using a for loop is very simple. For example:

p(1) = 1;
for i = 1:6
  p2(i+1) = 2*p2(i);
end

the problem comes when I try to vectorize this kind of code. I've used things like this:

p(1) = 1;
p(2:6) = 2*p(1:5);

but it does not update the vector each time, so it is not the correct solution... I suppose there is a very simple way to do it, but I don't know how.

Thank you!

[Kareem Elgindy]

What if the relation was like this: p(i+2) = 2*p(i+1) - p(i) with p(1) = 1 and p(2) = 5? Can we still use filter?

Answer 1

[Jan]

Do not overestimate a vectorization.

n = 100000;
tic
p = 1;
for i = 1:n-1
  p(i+1) = 1.001 * p(i);
end
toc

tic
p = zeros(1, n);  % Pre-allocation
p(1) = 1;
for i = 1:n-1
  p(i+1) = 1.001 * p(i);
end
toc

tic
x = [1, zeros(1, n-1)];
p = filter(1, [1 -1.001], x);
toc

tic;
p = repmat(1.001, 1, n);
p(1) = 1;
p = cumprod(p);
toc

tic;
p = 1.001 .^ (0:n-1);
toc

 Elapsed time is 8.463157 seconds.  % Loop without pre-allocation
 Elapsed time is 0.001182 seconds.  % Loop with pre-allocation
 Elapsed time is 0.001233 seconds.  % filter
 Elapsed time is 0.000507 seconds.  % cumprod
 Elapsed time is 0.011692 seconds.  % power

The filter command is efficient here, cumprod also, but the main problem was the missing pre-allocation.

Note: These are only rough timings. Prefer timeit for more accurate values.

[mathango]

This is very nice and useful example that your have demonstrated for vectorization of the recursive problem. I wonder if it can be done with y=y+x, y=sin(y) and etc.?

Note that the second example where you employed p = zeros(1, n); % Pre-allocation did not work for my example. I had to use y=0.0;

Answer 2

[Guillaume]

As long as the recursive relationship is linear you can use filter:

x = [1, zeros(1, 6)]; %starting value, + how many elements desired
p2 = filter(1, [1 -2], x) %creates a(1)*p2(n) = b*x(n) - a(2)*p2(n-1)=> p2(1) = x(1); p2(n) = 2*p2(n-1)

If the relationship is non linear then you don't have a choice but to use a loop.