关于Julia和Matlab之间时序比较的奇怪观察

Question

这个实验的目的是将Matlab和Julia的速度与下面的一小段代码进行比较.

首先是Matlab代码:

>> t = 5000; n = 10000;  x = 1:t*n;
>> x = reshape(x, t, n);
>> tic(); y1 = sum(x(:) .* x(:)); toc()
Elapsed time is 0.229563 seconds.
>> y1

    y1 =
         4.1667e+22
>> tic(); y2 = trace(x * x'); toc()
Elapsed time is 15.332694 seconds.

>> y2
y2 =
       4.1667e+22

与朱莉娅对决

julia> t = 5000; n = 10000; x = 1: t*n;    
julia> x = reshape(x, t, n);
julia> tic(); y1 = sum(x[:].* x[:]); toc();
elapsed time: 1.235170533 seconds
julia> y1
-4526945843202100544
julia> tic();y2 = trace(x*x'); toc();

第二个没有在1分多钟内完成这项工作.那么朱莉娅在这里有什么问题？这段代码碰巧在Julia中运行得更慢而且溢出？我的风格有问题吗？我认为从Matlab转换为Julia的一个原因是速度,我曾经认为Julia默认处理大数字算术运算.现在,看起来这些都不正确.谁能解释一下？

Answer 1

这里有几件事情.

首先,与Matlab不同,它x是一个机器整数数组,而不是浮点值.这似乎是速度的主要差异,因为它无法使用BLAS例程进行线性代数.

你需要做任何一件事

x = 1.0:t*n

或明确转换它

x = float(x)

这也是答案不同的原因:Julia使用机器算术,对于整数将在溢出时回绕(因此为负数).浮点值不会出现此问题.

(Julia确实有任意精度的整数,但默认情况下不提升:相反,你需要自己推广它们big(n))

这会给你一些加速,但你可以做得更好.Julia确实需要与Matlab略有不同的编程风格.请参阅以上评论中Isaiah提供的链接.

关于你的具体例子

sum(x[:].* x[:])

这很慢,因为它创建了3个中间向量(2个副本x[:],但希望将来会改变,以及它们的产品).

同样的,

trace(x*x')

计算中间矩阵(其中大部分未在跟踪中使用).

我的建议是使用

dot(vec(x),vec(x))

vec(x)只是重塑x成一个向量(所以没有复制),并且dot是通常的总和产品.你也可以"滚动你自己":

function test(x)
    s = zero(eltype(x)) # this prevents type-instability
    for xi in x
        s += xi*xi
    end
    s
end

test(x)

(这需要在JIT编译器的功能中发挥作用).它应该相当快,尽管可能仍然没有那么快dot,它使用下面的BLAS调用.