这是我模拟 Levy 运动的小脚本:
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
stable2
返回 stable random variable使用给定的参数(对于这种情况,您可以将其替换为 normrnd(mu, sigma)
,这并不重要); qlines2
返回绘图所需的分位数。
但我不想在这里谈论数学。我的问题是这个实现很慢,我想加快它的速度。不幸的是,计算机科学不是我的主要领域——我听说过记忆化、矢量化等方法,还有很多其他技术,但我不知道如何使用它们。
例如,我很确定我应该以某种方式替换这个肮脏的双 for 循环,但我不确定该怎么做。
编辑:也许我应该使用(并学习......)另一种语言(Python,C,任何功能性语言)?我一直认为 Matlab/OCTAVE 是为数值计算而设计的,但如果改变,那是为哪个?
最佳答案
正如您所说,关键是 for 循环,Matlab 不喜欢这些,所以向量化确实是关键字。 (连同预分配空间。
我只是稍微改变了你的循环部分,这样你就不必一遍又一遍地重置 L
,而是将所有 L
保存在一个更大的矩阵中(我还消除了 length(L)
命令)。
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
现在您已经可以看到循环中的 X(k,1:I) = L(k,1:I);
已过时,这也意味着我们可以切换顺序循环。这是至关重要的,因为 i
-steps 是递归的(取决于前面的步骤),这意味着我们不能向量化这个循环,我们只能向量化 k
-loop。
现在你的原始代码在我的机器上需要 9.3 秒,新代码仍然需要大约相同的时间)
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;
但是现在我们可以应用向量化,而不是遍历所有行(k
上的循环)我们可以改为消除这个循环,并在“一次”。
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;
这段代码在我的机器上只需要 0.045 秒。我希望你仍然得到相同的输出,因为我不知道你在计算什么,但我也希望你能看到你是如何对代码进行矢量化的。
PS:我刚刚注意到我们现在在上一个示例中对整个列使用相同的随机数,这显然不是您想要的。相反,您应该生成一个完整的随机数向量,例如:
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:好问题!
关于performance - 加速 Levy 运动算法的仿真,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/34108255/