Speeding up simulation of the Levy motion algorithm

Question

Here is my little script for simulating Levy motion:

clear all;
clc; close all;
t = 0; T = 1000; I = T-t;
dT = T/I; t = 0:dT:T; tau = T/I;
alpha = 1.5;
sigma = dT^(1/alpha);
mu = 0; beta = 0;
N = 1000;
X = zeros(N, length(I));
for k=1:N
    L = zeros(1,I);
    for i = 1:I-1
       L( (i + 1) * tau ) = L(i*tau) + stable2( alpha, beta, sigma, mu, 1);
    end
    X(k,1:length(L)) = L;
end

q = 0.1:0.1:0.9;
quant = qlines2(X, q, t(1:length(X)), tau);
hold all
for i = 1:length(quant)
    plot( t, quant(i) * t.^(1/alpha), ':k' );
end

Where stable2 returns a stable random variable with given parameters (you may replace it with normrnd(mu, sigma) for this case, it's not crucial); qlines2 returns quantiles needed for plotting.

But I don't want to talk about math here. My problem is that this implementation is pretty slow, and I would like to speed it up. Unfortunately, computer science is not my main field - I heard something about methods like memoization, vectorization and that there is a lot of other techniques, but I don't know how to use them.
For example, I'm pretty sure I should replace this filthy double for-loop somehow, but I'm not sure what to do instead.
EDIT: Maybe I should use (and learn...) another language (Python, C, any functional one)? I always though that Matlab/OCTAVE is designed for numerical computation, but if change, then for which one?

Answer 1

The crucial bit is, as you said, the for loops, Matlab does not like those, so vectorization is indeed the keyword. (Together with preallocating the space.

I just altered you for loop section somewhat so that you do not have to reset L over and over again, instead we save all Ls in a bigger matrix (also I elimiated the length(L) command).

L = zeros(N,I);
for k=1:N
    for i = 1:I-1
       L(k,(i + 1) * tau ) = L(k,i*tau) + normrnd(mu, sigma);
    end
    X(k,1:I) = L(k,1:I);
end

Now you can already see that X(k,1:I) = L(k,1:I); in the loop is obsolete and that also means that we can switch the order of the loops. This is crucial, because the i-steps are recursive (depend on the previous step) that means we cannot vectorize this loop, we can only vectorize the k-loop.

Now your original code needed 9.3 seconds on my machine, the new code still needs about the same time)

L = zeros(N,I);
for i = 1:I-1
    for k=1:N
       L(k,(i + 1) * tau ) = L(k,i*tau) + normrnd(mu, sigma);
    end
end
X = L;

But now we can apply the vectorization, instead of looping throu all rows (the loop over k) we can instead eliminate this loop, and doing all rows at "once".

L = zeros(N,I);
for i = 1:I-1
   L(:,(i + 1) * tau ) = L(:,i*tau) + normrnd(mu, sigma); %<- this is not yet what you want, see comment below
end
X = L;

This code need only 0.045 seconds on my machine. I hope you still get the same output, because I have no idea what you are calculating, but I also hope you could see how you go about vectorizing code.

PS: I just noticed that we now use the same random number in the last example for the whole column, this is obviously not what you want. Instad you should generate a whole vector of random numbers, e.g:

L = zeros(N,I);
for i = 1:I-1
   L(:,(i + 1) * tau ) = L(:,i*tau) + normrnd(mu, sigma,N,1);
end
X = L;

PPS: Great question!