CODEWITHSUNDEEP
performancevectorizationoctavesparse-matrix
Mattia
Mattia

Reputation: 13

Sparse matrix to speed up octave

I have a loop where "i" depends on "i-1" value, so I cannot vectorize it. I've read that I can use a sparse matrix in order to vectorize it and so to speed up my code, but I don't understand how this work. Any help? Thanks

Upvotes: 1

Views: 224

Answers (1)

Tasos Papastylianou
Tasos Papastylianou

Reputation: 22225

You are referring to this technique, as referenced from this (rather old) how to speed up octave article.

I'll rephrase the gist here in case the link dies in the future.

Suppose you have the following loop:

p1(1) = 0;

for i = 2 : N
    t     = t + dt;
    p1(i) = p1(i - 1) + dt * 2 * t;
endfor

You note here that, purely from a mathematical point of view, the last step in the loop could be rephrased as:

-1 * p1(i - 1) + 1 * p1(i) = dt * 2 * t

This makes it possible to recast the problem as a sparse matrix solve, by thinking of p1 as the vector of unknowns, and each iteration of the loop as a row in a (sparse) system of equations. E.g.:

Given that t is a known vector, this makes the above a straightforward problem that can be solved via a simple matrix division operation, which is guaranteed to be fast.

Having said that, presumably this 'trick' is only useful if you are able to recast the problem in this manner in the first place. Presumably this will only be the case for linear problems of your unknown. I don't think this can necessarily be used for more complicated loops.

Also, as Cris has mentioned in the comments, if this method does not work for you, there's a chance you can optimize your loop in other ways (or even that the loop solution may not necessarily be slow in the first place).

By the way, in theory, Octave provides jit-speedup like matlab does, though unlike matlab you need to enable it explicitly (in the sense that you need to compile your octave with jit options, which tends not to be the default), and my personal experience is that this is mostly experimental and may not do much except in the simplest of loops (see this post).

Upvotes: 3

Related Questions