julia - 欧几里得距离矩阵之间的最小距离

Question

我有一些代码可以计算一个矩阵中的每个笛卡尔坐标与另一个矩阵中的每个其他坐标之间的距离。对于每个坐标，最小距离将与产生最小值的坐标的索引位置一起返回。

function MED3D(m1, m2)
    n1::Int = size(m1,1)
    Dist = SharedArray{Float64}((n1,3))
    @sync @distributed for k in 1:n1
        Dist[k,:] = MD3D(m1[k,:], m2, k)
    end
    return Dist
end

@everywhere function MD3D(v1, m2, k)
    dsum::Float64 = Inf
    dtemp::Float64 = Inf
    i = 0
    for j in 1:size(m2,1)
        @inbounds dtemp = sqrt((v1[1] - m2[j,1]) * (v1[1] - m2[j,1]) + (v1[2] - m2[j,2]) * (v1[2] - m2[j,2]) + (v1[3] - m2[j,3]) * (v1[3] - m2[j,3]))
        if dtemp < dsum
            dsum = dtemp
            i = j
        end
    end
    return [dsum, k, i]
end

m1 = rand(10,3)
m2 = rand(15,3)
results = MED3D(m1,m2)

虽然这对于较小的 3D 点云来说相当不错，但我希望通过使用基于 GPU 的分析来提高大型点云的性能。然而，在 Julia 中使用更典型的方法进行矩阵运算似乎是不可能的，因为我必须返回索引位置和最小距离。我已经尝试了几种不同的方法来为这个任务采用 CUarrays，但到目前为止，它们都失败了，而没有使用实际的 for 循环。此外，由于将距离矩阵存储在内存中，许多实现它的方法似乎效率极低，对于我的特定数据集，它很快超过了 128gb 的 ram。

有人可以帮助我如何在 Julia 中正确实现它以在 GPU 上运行吗？ CUarrays 甚至是正确的方法，还是考虑到除了距离之外我还返回索引，它是否过于抽象？我尝试使用乘积和点来计算 L2 范数，但它并没有完全满足我的要求。

更新:

这是我使用广播对内循环进行 GPU 化的失败尝试。

using CuArrays
function difff(m1,m2)
    n1 = size(m1,1)
    Dist = Array{Float64}(undef, n1,3)
    m2 = CuArray(m2)
    m1 = CuArray(m1)
    for z in 1:size(m1)
        v1 = transpose(m1[z,:])
        i = 0
        dsum::Float64 = Inf
        mi = v1 .- m2
        mi = mi .* mi
        mi = sum(mi, dims=2)
        mi = mi .^ 0.5
        mi = findmin(mi)
        i = mi[2][1]
        dsum = mi[1]
        @inbounds Dist[z,:] = [dsum,z,i]
    end
end

更新:

失败的尝试 #2。我试图计算最小距离而忘记了指数。这对我的应用程序来说并不理想，但我可以忍受。但是，这仅在第一个数组只有一行时才能正常工作。我试图通过使用 map 切片来解决这个问题，但它不起作用。

using CuArray
a = rand(1,3)
b = rand(3,3)

a = CuArray(a)
b = CuArray(b)

function GK(m1, m2)
    reduce(min, sum((m1 .- m2) .^ 2,dims=2) .^ 0.5)
end

mapslices(GK(b), a, 2)

更新:

通过使用外循环取得进展，但肯定有更好的方法来做到这一点？

using CuArray
using BenchmarkTools
aa = rand(2,3)
bb = rand(5000000,3)

a = CuArray(aa)
b = CuArray(bb)

function GK(m1, m2)
    reduce(min, sum((m1 .- m2) .^ 2,dims=2) .^ 0.5)
end

function D(a,b)
    Dist = Array{Float64}(undef,size(a,1),1)
    for i in 1:size(a,1)
        Dist[i] = GK(a[i,:]',b)
    end
    return Dist
end

@benchmark test = D(a,b)
@benchmark test = D(aa,bb)

更新:

我之前的分布式版本、修改后的分布式版本、GPU 版本和串行版本之间的一些基准测试。 ??

using Distributed
using SharedArrays
using CuArrays
using BenchmarkTools

aa = rand(4,3)
bb = rand(500000,3)
a = CuArray(aa)
b = CuArray(bb)

function MED3D(m1, m2)
    n1::Int = size(m1,1)
    Dist = SharedArray{Float64}((n1,1))
    @sync @distributed for k in 1:n1
        Dist[k] = MD3D(m1[k,:]', m2)
    end
    return Dist
end

@everywhere function MD3D(v1, m2)
    dsum::Float64 = Inf
    dtemp::Float64 = Inf
    for j in 1:size(m2,1)
        @inbounds dtemp = sqrt((v1[1] - m2[j,1]) * (v1[1] - m2[j,1]) + (v1[2] - m2[j,2]) * (v1[2] - m2[j,2]) + (v1[3] - m2[j,3]) * (v1[3] - m2[j,3]))
        if dtemp < dsum
            dsum = dtemp
        end
    end
    return dsum
end

function MED3DGK(m1, m2)
    n1::Int = size(m1,1)
    Dist = SharedArray{Float64}((n1,1))
    @sync @distributed for k in 1:n1

        @inbounds Dist[k] = GK(m1[k,:]',m2)
    end
    return Dist
end

@everywhere function GK(m1, m2)
    reduce(min, sum((m1 .- m2) .^ 2,dims=2) .^ 0.5)
end

function D(a,b)
    Dist = Array{Float64}(undef,size(a,1),1)
    for i in 1:size(a,1)
        @inbounds Dist[i] = GK(a[i,:]',b)
    end
    return Dist
end

@benchmark test = D(a,b)
@benchmark test = D(aa,bb)
@benchmark test = MED3D(aa,bb)
@benchmark test = MED3DGK(aa,bb)

更新:

使用 NearestNeighbors.jl 实现分布式处理。关于如何使这更快的任何想法？:

function MED3D(m1, m2)
    m2 = Matrix(m2')
    kdtree = KDTree(m2)
    n1::Int = size(m1,1)
    Dist = SharedArray{Float64}((n1,1))
    Ind = SharedArray{Float64}((n1,1))
    @sync @distributed for k in 1:n1
        Ind[k,:], Dist[k,:] = knn(kdtree, m1[k,:], 1)
    end
    return [Ind,Dist]
end

Answer 1

我不确定它对你的情况是否有帮助，但是当切片时 m1[k,:]默认情况下，julia 会复制该内存(尽管这可能取决于 knn 函数对该切片的作用。
如果你把它改成 knn(kdtree, @view m1[k,:], 1) 有什么改善吗？