Vectorization recursive vector operation

4 次查看（过去 30 天）
Alejandro Castilla
评论： mathango 2021-3-29
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!
  1 个评论
Kareem Elgindy
Kareem Elgindy 2020-7-16
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?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Jan
Jan 2017-5-2
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.
  1 个评论
mathango
mathango 2021-3-29
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;

请先登录，再进行评论。

Guillaume
Guillaume 2017-5-2
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.

类别

在 Help Center 和 File Exchange 中查找有关 Performance and Memory 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by