help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Vectorizing for loops .. can't get the hang of it


From: Jordi Gutiérrez Hermoso
Subject: Re: Vectorizing for loops .. can't get the hang of it
Date: Tue, 6 Nov 2012 13:08:01 -0500

On 6 November 2012 11:31, Joza <address@hidden> wrote:
> I understand I should REALLY try use vectors instead of for loops ... anyone
> care to explain how to do it? I come from a CS programming background so
> traditional for loops always seem the easiest to me.

There is nothing traditional about loops in CS. If you have a wide
enough programming background, you should be familiar with procedural,
object-oriented and functional programming, to mention a few
programming paradigms. An example of very loopless code is functional
programming, and Octave supports some basic functional constructs with
arrayfun, structfun, and cellfun.

You can also avoid loops in object-oriented programming, since you can
*implement* loops with objects, quoting Steve Yegge:

    // print even numbers from 2 to 100
    LoopObject loop = new LoopObject();
    loop.setStartingValue(2);
    loop.setEndingValue(10);
    loop.setCounterIncrement(2);
    loop.setOperation(new LoopBody() {
        public void loopBody(int currentIndex) {
            System.out.println(currentIndex);
        }
    });
    loop.start();

(from
https://sites.google.com/site/steveyegge2/language-trickery-and-ejb)

My point here isn't to write that Java monstrosity to not write loops
in Java, but rather, that loops aren't the only hammer in your toolbox
just because you consider yourself a programmer.

> A piece of code I've used recently is below. I couldn't figure out
> how to vectorize it. Where do i start?

The first step is to attempt to not think in terms of indices. In
fact, you initially do not have to think of indices, you think of
"first column" or "diagonal" or "previous column". For writing a loop,
you have to jump from this interpretation to giving an subscript pair
(i,j) to that column, row, or diagonal. A step to thinking vectorially
is to avoid this step of assigning indices to begin with.

Another thing you should do is recognise matrix multiplication, since
y = A*x is

    y = zeros(rows(A), columns(x));
    for i = 1:rows(A)
        for j = 1:columns(A) ## where columns(A) == rows(x)
            y(i) += A(i,j)*x(j);
        endfor
    endfor

So, to unloop your code...

> for i=1:n
>                 sum = 0.0;
>
>                 if i==1
>                         for j=1:2
>                                 sum = sum + A(i,j)*x(j);
>                         end
>                 elseif i==n
>                         for j=n-1:n
>                                 sum = sum + A(i,j)*x(j);
>                         end
>                 else
>                         for j=i-1:i+1
>                                 sum = sum + A(i,j)*x(j);
>                         end
>                 end
> end

the same expression as matrix multiplication occurs here, A(i,j)*x(j),
but I see that you've special-cased your loops in order to especially
handle zero entries, since you're using sparse matrices. If just do

    A = sparse (A);
    x = sparse (x);
    y = A*x;

assuming A and x weren't sparse to begin with, you can just let A's
zero entries handle the special cases you branched into ifs.

Your code, however, is attempting to sum all of the entries into a
single "sum" variable (which, btw, overwrites a core function called
"sum" too), so your loop and branching, assuming A is tridiagonal as
you've frequently been stating, should be instead

    s = sum (y);

If A is not tridiagonal but you only wanted to extract the tridiagonal
sum as in your original loopy implementation,

    ## extract just the tridiagonal part of A
    A = sparse (triu (tril (A,1),-1));

    ## make x sparse in case it  wasn't already.
    x = sparse(x);

    ## Compute the sum of the entries of the matrix-vector product
    s = sum (A*x);

HTH,
- Jordi G. H.


reply via email to

[Prev in Thread] Current Thread [Next in Thread]