create a vector whose elements depend on its previous elements without a loop

Hi,

As loops take considerable time, I'm trying without loops.

I want to create a vector whose elements depend on its previous elements, for example:

a(i) = 3*a(i-1) + 2

a(1) = 5

where a(i) is the i-th value of vector A.

Can it be done?

[Stephen23]

"As loops take considerable time..."

Loops do not take "considerable time". Badly written code inside a loop might be slow, or an algorithm might require slow code or many iterations, but loops in themselves are fast.

Have you read the documentation on how to write efficient MATLAB code?:

https://www.mathworks.com/help/matlab/matlab_prog/techniques-for-improving-performance.html

In particular, did you preallocate any output arrays before the loop?

[Galgool]

I know I can write the loops much more efficiantly, and that in my code I won't get rid of all of them and I need the outer loop to do a few million iterations...

But probably I can get rid of the inner loops with Element-wise calculations. Specificaly, I need to create a vector whose elements depend on its previous elements.

Can This be done?

[Guillaume]

Yes, it can be done. Haven't you seen my answer?

Have you proven yet that the loop is a bottleneck? If not, why are you trying to optimise it? You're probably focusing on the wrong part of your code.

Note that while it can be done without a loop, testing on my computer shows that it is signficantly faster with a loop:

>> filter_vs_loop(1e6)
Comparing array creation of 1000000 elements, with a loop or filter

With loop:   0.00579613 s
With filter: 0.0146097 s

With loop:   0.00560415 s
With filter: 0.0146441 s

With loop:   0.00564641 s
With filter: 0.0145454 s

With loop:   0.00556884 s
With filter: 0.0147495 s

With loop:   0.00566302 s
With filter: 0.0146408 s

[Galgool]

Thanks a lot. Took me a while to understand filter, but it works beautifuly.

Is there a way to use filter (or any other way) for:

a(i) = f(i) * a(i-1) + 2

(ofcourse a(0)==0)

where the coefficient f(i) is also a vector.

[Guillaume]

filter can't be used for that and I doubt there's anything faster than a loop.

[Jan]

This might be 10% faster than fillwithfilter:

function a = fillwithfilter2(nelem)
    q          = zeros(1, nelem);
    q(1)       = 5;
    q(2:nelem) = 10;
    a = filter(1, [1, -3], q);
end

But this is of academic interest only, if the loop is much faster.

Answer 1

[Guillaume]

Doing this will a loop shouldn't be slow as long as you preallocate the array beforehand:

a = zeros(1, 100);  %preallocation
a(1) = 5;
for i = 2:numel(a)
    a(i) = 3 * a(i-1) + 2;
end

It can indeed be done without a loop, with the filter function, the most difficult is to figure out the inputs. The more about section is the most helpful.

a = filter(1, [1 -3], [5, 2*ones(1, 99)]);  % 1*y(n) = 1*x(n) - (-3)*y(n-1), with x = [5, 2, 2, ...] and y(0) = 0

Frankly, the loop is clearer so you should prefer it.

[Jan]

+1: The filter method is faster, most likely. So you can save some milliseconds, but you loose a minute or 10 for writing and testing.

By the way, mulitplications take more cpu cycles than addition. So why is 2*ones(1, n) not measurably slower than 1+ones(1, n) ?

[Guillaume]

If we're trying to eke out a few cpu cycles from 2*ones vs 1+ones, we're doing something wrong!. With that said,

a = filter(2, [1, -3], [2.5, ones(1, 99)])

would possibly gain you an extra cycle. (but you would have probably spend more time figuring it out than you'd ever gained over the lifetime of the program).

[Guillaume]

Turns out that filter is actually about 3 times slower than a loop (R2019a)

[Jan]

@Guillaume: Even if we do something wrong (a kind of pre-mature optimization), CPUs are deterministic machines and an arithmetic operation, which needs less cycles, should consume less time.

x = ones(1, 1e6);
timeit(@() x * 2)   % 1.2651e-04
timeit(@() x + 1)   % 1.2923e-04

By the way, the timings are equivalent for x * 2.22 and x + 2.22, so this is not a magic effect of the JIT performing some smart repelacement for the multiplication by 2.

Nevertheless, this does not concern the question of this thread. I'm going to discuss this somewhere else.