我刚开始使用 Julia，并将我的 MATLAB 代码翻译成 Julia（基本上是逐行）。我注意到 Julia 代码要慢得多（比如 50x）。最初的问题是一个动态编程问题，我在其中插值函数——插值是代码花费大部分时间的地方。因此，我尝试制作一个显示性能差异的最小示例代码。需要注意的重要事项是，它是插值的样条近似，并且网格最好是不规则的，即不等间距。MATLAB 代码：
tic
spacing=1.5;
Nxx = 300;
Naa = 350;
Nalal = 200;
sigma = 10;
NoIter = 500;
xx=NaN(Nxx,1);
xmin = 0.01;
xmax = 400;
xx(1) = xmin;
for i=2:Nxx
xx(i) = xx(i-1) + (xmax-xx(i-1))/((Nxx-i+1)^spacing);
end
f_U = @(c) c.^(1-sigma)/(1-sigma);
W=NaN(Nxx,1);
W(:,1) = f_U(xx);
xprime = ones(Nalal,Naa);
for i=1:NoIter
W_temp = interp1(xx,W(:,1),xprime,'spline');
end
toc
Elapsed time is 0.242288 seconds.
朱莉娅代码：
using Dierckx
function performance()
const spacing=1.5
const Nxx = 300
const Naa = 350
const Nalal = 200
const sigma = 10
const NoIter = 500
xx=Array(Float64,Nxx)
xmin = 0.01
xmax = 400
xx[1] = xmin
for i=2:Nxx
xx[i] = xx[i-1] + (xmax-xx[i-1])/((Nxx-i+1)^spacing)
end
f_U(c) = c.^(1-sigma)/(1-sigma)
W=Array(Float64,Nxx,1)
W[:,1] = f_U(xx)
W_temp = Array(Float64,Nalal,Naa)
interp_spline = Spline1D(xx,W[:,1])
xprime = ones(Nalal,Naa)
for i=1:NoIter
for j=1:Naa
for k=1:Nalal
W_temp[k,j] = evaluate(interp_spline, xprime[k,j])
end
end
end
end
@time(performance())
4.200732 seconds (70.02 M allocations: 5.217 GB, 4.35% gc time)
我制作了 W 一个二维数组，因为在原始问题中它是一个矩阵。我对不同的插值包进行了一些研究，但对于不规则网格和样条曲线没有太多选择。MATLAB 的 interp1 显然不可用。
我的问题是我在想，如果我编写 Julia 代码并给出与 MATLAB 相同的结果，那么 Julia 通常应该更快。但显然情况并非如此，因此您需要注意您的编码。我不是程序员，当然我会尽力而为，但我想知道我是否在这里犯了一些可以轻松修复的明显错误，或者它是否会（太）经常发生，我必须注意我的 Julia 编码——因为那样我可能不值得学习它。同样，如果我可以在这里让 Julia 更快（我很确定我可以，例如分配看起来有点大），我也可以让 MATLAB 更快。我对 Julia 的希望是——对于类似的代码——它会比 MATLAB 运行得更快。
在对时间进行了一些评论之后，我还运行了以下代码：
using Dierckx
tic()
const spacing=1.5
const Nxx = 300
const Naa = 350
const Nalal = 200
const sigma = 10
const NoIter = 500
xx=Array(Float64,Nxx)
xmin = 0.01
xmax = 400
xx[1] = xmin
for i=2:Nxx
xx[i] = xx[i-1] + (xmax-xx[i-1])/((Nxx-i+1)^spacing)
end
f_U(c) = c.^(1-sigma)/(1-sigma)
W=Array(Float64,Nxx,1)
W[:,1] = f_U(xx)
W_temp = Array(Float64,Nalal,Naa)
interp_spline = Spline1D(xx,W[:,1])
xprime = ones(Nalal,Naa)
for i=1:NoIter
for j=1:Naa
for k=1:Nalal
W_temp[k,j] = evaluate(interp_spline, xprime[k,j])
end
end
end
toc()
elapsed time:
7.336371495 seconds
甚至慢得多，嗯...
另一个编辑：在这种情况下，消除一个循环实际上使它更快，但仍然无法与 MATLAB 相比。编码：
function performance2()
const spacing=1.5
const Nxx = 300
const Naa = 350
const Nalal = 200
const sigma = 10
const NoIter = 500
xx=Array(Float64,Nxx)
xmin = 0.01
xmax = 400
xx[1] = xmin
for i=2:Nxx
xx[i] = xx[i-1] + (xmax-xx[i-1])/((Nxx-i+1)^spacing)
end
f_U(c) = c.^(1-sigma)/(1-sigma)
W=Array(Float64,Nxx,1)
W[:,1] = f_U(xx)
W_temp = Array(Float64,Nalal,Naa)
interp_spline = Spline1D(xx,W[:,1])
xprime = ones(Nalal,Naa)
for i=1:NoIter
for j=1:Naa
W_temp[:,j] = evaluate(interp_spline, xprime[:,j])
end
end
end
@time(performance2())
1.573347 seconds (700.04 k allocations: 620.643 MB, 1.08% gc time)
另一个编辑，现在通过循环迭代相同次数：
function performance3()
const spacing=1.5
const Nxx = 300
const Naa = 350
const Nalal = 200
const sigma = 10
const NoIter = 500
xx=Array(Float64,Nxx)
xmin = 0.01
xmax = 400
xx[1] = xmin
for i=2:Nxx
xx[i] = xx[i-1] + (xmax-xx[i-1])/((Nxx-i+1)^spacing)
end
f_U(c) = c.^(1-sigma)/(1-sigma)
W=Array(Float64,Nxx,1)
W[:,1] = f_U(xx)
W_temp = Array(Float64,Nalal,Naa)
W_tempr = Array(Float64, Nalal*Naa)
interp_spline = Spline1D(xx,W[:,1])
xprime = ones(Nalal,Naa)
xprimer = reshape(xprime, Nalal*Naa)
for i=1:NoIter
W_tempr = evaluate(interp_spline, xprimer)
end
W_temp = reshape(W_tempr, Nalal, Naa)
end
tic()
performance3()
toc()
elapsed time:
1.480349334 seconds
不过，与 MATLAB 完全没有可比性。顺便说一句，在我的实际问题中，我轻松地运行了 50k 次循环，并且我正在访问更大的 xprime 矩阵，尽管不确定该部分是否有所作为。