Vectorizing ODE in Octave / Matlab

Question

If I have an ode and wrote it in two ways, like here:

function re=rabdab()
x=linspace(0,2000,2000)';
tic;
[T,Y] = ode45(@fun,[x],[0 1 1]);
[T,Y] = ode45(@fun,[x],[0 1 1]);
[T,Y] = ode45(@fun,[x],[0 1 1]);
toc;

tic;
[A,B] = ode45(@fun2,[x],[0 1 1]);
[A,B] = ode45(@fun2,[x],[0 1 1]);
[A,B] = ode45(@fun2,[x],[0 1 1]);
toc;



function dy = fun(t,y)
dy = zeros(3,1);    % a column vector
dy = [y(2) * y(3);...
-y(1) * y(3);...
-0.51 * y(1) * y(2);];

function dy = fun2(t,y)
dy = zeros(3,1);    % a column vector
dy(1) = y(2) * y(3);
dy(2) = -y(1) * y(3);
dy(3) = -0.51 * y(1) * y(2);

There is almost no difference in time. One takes just as long as the other. But I thought that fun is the vectorized version of fun2. Or am I mistaken here? The purpose is to speed up my code a little. The example is taken from the matlab webpage. I think I haven't really understood what "vectorized" means. If this is already vectorized, what would a non-vectorized code look like?

Answer 1

Vectorization is a concept which is closely related to functional programming. In MATLAB it means performing operations on arrays (vectors or matrices) without implicitly writing a loop.

For example, if you were to compute the function f(x) = 2x for every integer x between 1 and 100, you could write:

for x = 1:100
    f(x) = 2 * x;
end

which is not vectorized code. The vectorized version is:

x = 1:100; %// Declare a vector of integer values from 1 to 100
f = 2 * x; %// Vectorized operation "*"

or even shorter:

f = 2 * (1:100);

MATLAB is an interpreted language, so obviously the interpreter translates this into some kind of loop "under the hood", but it's optimized and is usually much faster than actually interpreting a loop (read this question for reference). Well, sort of -- it's been like that until the recent releases of MATLAB, where JIT acceleration has been integrated (read here).

Now getting back to your code: what you have here is two vectorized versions of your code: one that concatenates vertically three values and one that directly assigns these values into a column vector. That's basically the same thing. Should you have done it with an explicit for loop, this would not have been "vectorized". Regarding the actual performance gain from "vectorizing" a loop (that is, converting a for loop into a vectorized operation), this depends on the how fast the for loop actually was due to JIT acceleration in the first place.

It doesn't seem that there's much to be done to speed up your code. Your functions are pretty basic, so it boils down to the internal implementation of ode45, which you cannot modify.

If you're interested in further reading about vectorization and writing faster MATLAB code in general, here's an interesting article: Loren on the Art of MATLAB: "Speeding Up MATLAB Applications".

Happy coding!

Answer 2

While the previous answer is correct in general terms, vectorized in the context of ODEs means something more specific. In short, a function f(t,y) is vectorized iff f(t,[y1 y2 ...]) returns [f(t,y1) f(t,y2) ...], for y1,y2 column vectors. As per the documentation [1], "this allows the solver to reduce the number of function evaluations required to compute all the columns of the Jacobian matrix."

Functions fun3 and fun4 below are properly vectorized in the ODE sense:

function re=rabdab()
x=linspace(0,20000,20000)';
opts=odeset('Vectorized','on');

tic;
[T,Y] = ode45(@fun,[x],[0 1 1]);
[T,Y] = ode45(@fun,[x],[0 1 1]);
[T,Y] = ode45(@fun,[x],[0 1 1]);
toc;

tic;
[A,B] = ode45(@fun2,[x],[0 1 1]);
[A,B] = ode45(@fun2,[x],[0 1 1]);
[A,B] = ode45(@fun2,[x],[0 1 1]);
toc;

tic;
[A,B] = ode45(@fun3,[x],[0 1 1],opts);
[A,B] = ode45(@fun3,[x],[0 1 1],opts);
[A,B] = ode45(@fun3,[x],[0 1 1],opts);
toc;

tic;
[A,B] = ode45(@fun4,[x],[0 1 1],opts);
[A,B] = ode45(@fun4,[x],[0 1 1],opts);
[A,B] = ode45(@fun4,[x],[0 1 1],opts);
toc;


function dy = fun(t,y)
dy = zeros(3,1);    % a column vector
dy = [y(2) * y(3);...
-y(1) * y(3);...
-0.51 * y(1) * y(2);];

function dy = fun2(t,y)
dy = zeros(3,1);    % a column vector
dy(1) = y(2) * y(3);
dy(2) = -y(1) * y(3);
dy(3) = -0.51 * y(1) * y(2);

function dy = fun3(t,y)
dy = zeros(size(y)); % a matrix with arbitrarily many columns, rather than a column vector
dy = [y(2,:) .* y(3,:);...
-y(1,:) .* y(3,:);...
-0.51 .* y(1,:) .* y(2,:);];

function dy = fun4(t,y)
dy = [y(2,:) .* y(3,:);... % same as fun3()
-y(1,:) .* y(3,:);...
-0.51 .* y(1,:) .* y(2,:);];

(As a side remark: omitting the unnecessary memory allocation with zeros lets fun4 run slightly faster than fun3.)

• Q: What about vectorizing with respect to the first argument?
• A: For the ODE solvers, the ODE function is vectorized only with respect to the second argument. The boundary value problem solver bvp4c, however, does require vectorization with respect to the first and second arguments. [1]

The official documentation [1] provides further details on ODE-specific vectorization (see section "Description of Jacobian Properties").

[1] https://www.mathworks.com/help/releases/R2015b/matlab/ref/odeset.html