Re: help with vectorization of update rule

[Søren Hauberg]

tor, 21 10 2010 kl. 08:48 -0700, skrev grg:
> 
> Søren Hauberg wrote:
> > 
> > Hi,
> > 
> > You can do something like this
> > 
> >         ## Problem parameters
> >         A = 2;
> >         T = 10;
> >         x0 = rand (20, 1);
> >         
> >         ## Loop implementation
> >         x = x0;
> >         for it=1:T
> >             x(:,it+1)=A*x(:,it);
> >         endfor
> >         
> >         ## Vectorised implementation
> >         y = repmat (x0, 1, T+1);
> >         y = repmat (x0, 1, T+1) .* repmat (A.^(0:T), rows (x0), 1);
> >         
> >         ## Check (should return 1)
> >         isequal (x, y)
> >         
> > Søren
> > 
> > 
> > 
> > 
> 
> hi Søren,
> 
> Thanks a lot for your prompt response.
> 
> I tried your code and checked the execution times:
> 
>         ## Problem parameters
>         A = 2;
>       T = 10;
>         x0 = rand (20, 1);
> 
>       tic
>       ## Loop implementation
>         x = x0;
>         for it=1:T
>             x(:,it+1)=A*x(:,it);
>         endfor
>       toc
> 
>       tic
>       ## Vectorised implementation
>         y = repmat (x0, 1, T+1);
>         y = repmat (x0, 1, T+1) .* repmat (A.^(0:T), rows (x0), 1);
>       toc
> 
>       ## Check (should return 1)
>         isequal (x, y)
> 
> 
> The output was:
> 
>     octave:11> test_vec
>     Elapsed time is 0.00022 seconds.
>     Elapsed time is 0.0016 seconds.
>     ans =  1
>     octave:12> 
> 
> So, in this case the vectorized version of the code is significantly slower
> than the for loop.  Is this normal?

Try the following (more simple solution):

    z = x0 * (A.^(0:T));

That seems to be faster than the loopy solution even for small problems.

Søren