Speeding up Levy's motion algorithm simulation

Question

Here is my little script to simulate Levy's movement:

clear all; clc; close all; t = 0; T = 1000; I = Tt; 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 the given parameters (you can replace it with normrnd(mu, sigma) for this case, it doesn’t matter); qlines2 returns the quantiles needed to plot the graph.

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

Answer 1

The key bit is, as you said, for loops, Matlab doesn't like it, so vectorization is really a keyword. (Together with the preliminary distribution of space.

I just changed a little for the loop section so that you do not have to reset L again and again, instead we save all L in a larger matrix (I also deleted 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 a cycle is outdated, which also means that we can switch the order of the cycles. This is important because i -steps are recursive (depending on the previous step), which means that we cannot vectorize this loop, we can only vectorize k -loop.

Now your source code took 9.3 seconds on my machine, the new code is still needed at 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 vector, instead of looping all the lines (loop through k ), we can exclude this loop and do all the rows "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; 

Only 0.045 seconds are required for this code on my machine. Hope you still get the same result, because I have no idea what you are calculating, but I also hope that you can see how you go about vectorizing the code.

PS: I just noticed that now we use the same random number in the last example for the entire column, this is clearly not what you want. Instad, you must create a whole vector of random numbers, for example:

 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!