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

Question

我正在运行一个程序，该程序需要重复使用二维点集合中的成对距离计算，然后计算向量。这最终导致我的运行时出现瓶颈，因此我尝试将代码从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

我能够得到一个速度vs matlab使用朱莉娅0.5多线程。

在我的机器上（i5有4个核心）我得到以下时间：
matlab R2012a-8.5秒
julia 0.5单线程-fo3（）（见下文）-18.5秒
julia 0.5多线程-fo4（）（见下文）-4.5秒

也就是说，我能够使julia单线程函数的运行速度比matlab慢一倍，但多线程函数的运行速度比matlab快一倍。

很抱歉这是一个非常冗长的回答——我认为最好是全面的。我在下面发布了使用的每个内部函数以及主要函数-fo3（）和foo（4）。

1.单Thread：

下面的展开函数的目的是避免不必要的内存分配，并利用数组的对称性。从tim的回答来看，这其中的大部分似乎可以用0.6中带有点符号的单行来处理。

function pdist2!(pd, a)
    m = size(a, 1)
    for col in 1:m
        for row in (col + 1):m
            s = 0.0
            for i in 1:2
                @inbounds s += abs2(a[col, i] - a[row, i])
            end
            @inbounds pd[row, col] = pd[col, row] = sqrt(s)
        end
    end
end

function dotminustranspose!(xgv, ygv, a)
    m = size(a, 1)
    for col in 1:m
        @inbounds for row in (col + 1):m
            xgv[row, col] = a[col, 1] - a[row, 1]
            ygv[row, col] = a[col, 2] - a[row, 2]
            xgv[col, row] = - xgv[row, col]
            ygv[col, row] = - ygv[row, col]
        end
    end
end

function atan2!(th, ygv, xgv)
    for i in eachindex(ygv)
        @inbounds th[i] = atan2(ygv[i], xgv[i])
    end
end

function invpdp1!(fv, pd)
    for i in eachindex(pd)
        @inbounds fv[i] = 1 / (pd[i] + 1)
    end
end

function fv_times_cos_th!(xfv, fv, th)
    for i in eachindex(th)
        @inbounds xfv[i] = fv[i] * cos(th[i])
    end
end

function fv_times_sin_th!(yfv, fv, th)
    for i in eachindex(th)
        @inbounds yfv[i] = fv[i] * sin(th[i])
    end
end

function adsum2!(a, xfv, yfv)
    n = size(a, 1)
    for j in 1:n
        for i in 1:n
            @inbounds a[i, 1] += xfv[i, j]
            @inbounds a[i, 2] += yfv[i, j]
        end
    end
end

function foo3()
    a = reshape(sin(1:300), 150, 2)
    histarray = zeros(5000, 150, 2)
    pd = zeros(size(a, 1), size(a, 1))
    xgv = zeros(pd)
    ygv = zeros(pd)
    th = zeros(pd)
    fv = zeros(pd)
    xfv = zeros(pd)
    yfv = zeros(pd)
    for i in 1:5000
        pdist2!(pd, a)
        dotminustranspose!(xgv, ygv, a)
        atan2!(th, ygv, xgv)
        invpdp1!(fv, pd)
        fv_times_cos_th!(xfv, fv, th)
        fv_times_sin_th!(yfv, fv, th)
        adsum2!(a, xfv, yfv)

        histarray[i, :, :] = view(a, :)
    end
    return histarray
end

时间：

@time histarray = foo3()
17.966093 seconds (24.51 k allocations: 13.404 MB)

1.多线程：

elementwise trig函数可以使用@threads宏进行多线程处理。这给了我大约4倍的加速。这仍然是实验性的，但我测试了输出，它们是相同的。

using Base.Threads

function atan2_mt!(th, ygv, xgv)
    @threads for i in eachindex(ygv)
        @inbounds th[i] = atan2(ygv[i], xgv[i])
    end
end

function fv_times_cos_th_mt!(xfv, fv, th)
    @threads for i in eachindex(th)
        @inbounds xfv[i] = fv[i] * cos(th[i])
    end
end

function fv_times_sin_th_mt!(yfv, fv, th)
    @threads for i in eachindex(th)
        @inbounds yfv[i] = fv[i] * sin(th[i])
    end
end

function foo4()
    a = reshape(sin(1:300), 150, 2)
    histarray = zeros(5000, 150, 2)
    pd = zeros(size(a, 1), size(a, 1))
    xgv = zeros(pd)
    ygv = zeros(pd)
    th = zeros(pd)
    fv = zeros(pd)
    xfv = zeros(pd)
    yfv = zeros(pd)
    for i in 1:5000
        pdist2!(pd, a)
        dotminustranspose!(xgv, ygv, a)
        atan2_mt!(th, ygv, xgv)
        invpdp1!(fv, pd)
        fv_times_cos_th_mt!(xfv, fv, th)
        fv_times_sin_th_mt!(yfv, fv, th)
        adsum2!(a, xfv, yfv)

        histarray[i, :, :] = view(a, :)
    end
    return histarray
end

时间：

@time histarray = foo4()
4.569416 seconds (54.51 k allocations: 14.320 MB, 0.20% gc time)

Answer 2

摆脱trig重构背后的思想是sin（atan（x，y））==y/sqrt（x^2y^2）。函数hypot方便地计算平方根分母inv用于摆脱慢速分割。守则：

# a constant input matrix to allow foo2/foo3 comparison
a = randn(150,2)

# calculation using trig functions
function foo2(b,n)
  a = copy(b)
  historyarray = zeros(n,size(a,1),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:n
      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

# helper function to handle annoying Infs from self interaction calc
nantoone(x) = ifelse(isnan(x),1.0,x)
nantozero(x) = ifelse(isnan(x),0.0,x)

# calculations using Pythagoras
function foo3(b,n)
  a = copy(b)
  historyarray = zeros(5000,size(a,1),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:n
      pairwise!(pd, Euclidean(),a.')
      xgv .= a[:,1].' .- a[:,1]
      ygv .= a[:,2].' .- a[:,2]
      th .= inv.(hypot.(ygv,xgv))
      fv .= inv.(pd.+1)
      xfv .= nantoone.(fv.*xgv.*th)
      yfv .= nantozero.(fv.*ygv.*th)
      a[:,1:1] .+= sum!(tmp, xfv)
      a[:,2:2] .+= sum!(tmp, yfv)
      historyarray[i,:,:] = a
  end
end

还有一个长凳：

julia> @time foo2(a,5000)
  9.698825 seconds (809.51 k allocations: 67.880 MiB, 0.33% gc time)

julia> @time foo3(a,5000)
  2.207108 seconds (809.51 k allocations: 67.880 MiB, 1.15% gc time)

A.

另一个收获是可以添加到Base的NaN到某物函数的便利性（类似于来自SQL世界的coalesce和nvl）。

Answer 3

您是正确的，调用sin，cos和atan2是Julia代码中的瓶颈。然而，大量拨款意味着仍有优化的潜力。

在即将发布的Julia版本中，您的代码可以通过使用改进的点广播语法a轻松重写，以避免不必要的分配f、 （b，c）。这相当于广播！（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%

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