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?

The script doesn't look too complex, so I'll hazard a guess you will have reasonable performance with an optimized Matlab implementation. Matlab is also pretty forgiving for a newbie - unless best possible performance is an absolute necessity, I feel you should stick with it for now. Especially since you've implemented your simulation there now :) — mikkola, Dec 05 '15 at 17:13

flawr · Accepted Answer · 2015-12-05T17:27:42.507

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!

Nice work. I believe for the simulation it is necessary to draw multiple samples with `normrnd` rather than one. This is achieved with the syntax `normrnd(mu, sigma, M, N, ...)` where `M`, `N` etc. are the desired dimensions of the output array. And the [usual caveats of using `i` or `j` as loop variables](http://stackoverflow.com/questions/14790740/using-i-and-j-as-variables-in-matlab) apply. — mikkola, Dec 05 '15 at 17:26
@mikkola Thank you for letting me know, I just noticed this too and added a paragraph. — flawr, Dec 05 '15 at 17:28
It works glorious, thank you! And yes, correction (your PS) was needed, since I wanted matrix of simulated paths. — Photon Light, Dec 05 '15 at 17:39
@PhotonLight And another thing, this code now works only because `tau` is exactly 1. — flawr, Dec 05 '15 at 17:40
Yes, I see, that's another issue, but for now I can assume it's true - in reality it represents time step, for simple cases it can be 1. I will talk with my prof. if it can be solved, maybe another implementation is needed, but now I know more about optimization, so should be easier! — Photon Light, Dec 05 '15 at 17:46

Speeding up simulation of the Levy motion algorithm

1 Answers1