与Matlab2016b相比,如何在Julia 0.5中获得更快的成对距离和矢量计算？

Question

我正在运行一个程序,需要在2D中的一组点上重复使用成对距离计算,然后计算向量.这最终导致我的运行时间出现瓶颈,因此我尝试将我的代码从Matlab重新编写为Julia,以利用其更快的速度.然而,我遇到的问题是我在Julia中编写的函数实际运行速度比我的Matlab实现慢五倍.鉴于朱莉娅的声誉是一种快得多的语言,我假设我做错了什么.

我写了一个简单的例子来说明我所看到的.

朱莉娅代码:

using Distances
function foo()
  historyarray = zeros(5000,150,2)
  a = randn(150,2)
  for i in 1:5000
    pd = pairwise(Euclidean(),a.')
    xgv = broadcast(-,a[:,1].',a[:,1])
    ygv = broadcast(-,a[:,2].',a[:,2])
    th = atan2(ygv,xgv)
    fv = 1./(pd+1)
    xfv = fv*cos(th)
    yfv = fv*sin(th)
    a[:,1]+= sum(xfv,2)
    a[:,2]+= sum(yfv,2)

    historyarray[i,:,:] = copy(a)
  end
end

Matlab代码:

function foo
histarray = zeros(5000,150,2);
a = randn(150,2);
for i=1:5000

    pd = pdist2(a,a);
    xgv = -bsxfun(@minus,a(:,1),a(:,1)');
    ygv = -bsxfun(@minus,a(:,2),a(:,2)');
    th = atan2(ygv,xgv);
    fv = 1./(pd+1);
    xfv = fv.*cos(th);
    yfv = fv.*sin(th);
    a(:,1) = a(:,1)+sum(xfv,2);
    a(:,2) = a(:,2)+sum(yfv,2);

    histarray(i,:,:)=a;
end

结束

当我测试Julia代码的速度(多次考虑编译时间)时,我得到:

@time foo()
16.110077 seconds (2.65 M allocations: 8.566 GB, 6.30% gc time)

另一方面,Matlab函数的性能是:

tic
foo
toc
Elapsed time is 3.339807 seconds. 

当我在Julia代码上运行配置文件查看器时,花费最多时间的组件是第9,11和12行.三角函数是否可能发生奇怪的事情？

Answer 1

你是正确的调用sin,cos以及atan2在你的朱莉娅代码的瓶颈.然而,大量的分配意味着仍有可能进行优化.

在即将推出的Julia版本中,您可以轻松地重写代码,以避免使用改进的点广播语法进行不必要的分配a .= f.(b,c).这相当于broadcast!(f, a, b, c)并更新a到位.此外,rhs上的几个点广播呼叫会自动融合为一个.最后,@views宏将所有切片操作(如a[:,1]复制)转换为视图.新代码如下所示:

function foo2()
    a = rand(150,2)
    historyarray = zeros(5000,150,2)
    pd = zeros(size(a,1), size(a,1))
    xgv = similar(pd)
    ygv = similar(pd)
    th = similar(pd)
    fv = similar(pd)
    xfv = similar(pd)
    yfv = similar(pd)
    tmp = zeros(size(a,1))
    @views for i in 1:5000
        pairwise!(pd, Euclidean(),a.')
        xgv .= a[:,1].' .- a[:,1]
        ygv .= a[:,2].' .- a[:,2]
        th .= atan2.(ygv,xgv)
        fv .= 1./(pd.+1)
        xfv .= fv.*cos.(th)
        yfv .= fv.*sin.(th)
        a[:,1:1] .+= sum!(tmp, xfv)
        a[:,2:2] .+= sum!(tmp, yfv)
        historyarray[i,:,:] = a
    end
end

(我xfv .= fv.*cos.(th)在Matlab代码中使用逐元素乘法而不是矩阵乘法.)

对新代码进行基准测试显示分配的内存大幅减少:

julia> @benchmark foo2()
BenchmarkTools.Trial: 
  memory estimate:  67.88 MiB
  allocs estimate:  809507
  --------------
  minimum time:     7.655 s (0.06% GC)
  median time:      7.655 s (0.06% GC)
  mean time:        7.655 s (0.06% GC)
  maximum time:     7.655 s (0.06% GC)
  --------------
  samples:          1
  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%

(大部分可以在0.5上实现,但需要更多输入)

但是,这仍然是您的Matlab版本的两倍.分析表明,大部分时间都花在三角函数上.

我只是为了好玩而尝试:

const atan2 = +
const cos = x->2x
const sin = x->2x

得到了:

julia> @benchmark foo2()
BenchmarkTools.Trial: 
  memory estimate:  67.88 MiB
  allocs estimate:  809507
  --------------
  minimum time:     1.020 s (0.69% GC)
  median time:      1.028 s (0.68% GC)
  mean time:        1.043 s (2.10% GC)
  maximum time:     1.100 s (7.75% GC)
  --------------
  samples:          5
  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%

我想,三角函数缓慢的一个原因可能是我使用了预构建的二进制文件,并且没有libmJulia使用的库的自编译版本.因此,libm代码未针对我的处理器进行优化.但我怀疑在这种情况下,这将使Julia比Matlab快得多.对于这种算法,Matlab代码似乎已经接近最优.