欧几里得距离矩阵之间的最小距离答案

【问题标题】：Minimum distances among a Euclidean distance matrix欧几里得距离矩阵之间的最小距离
【发布时间】：2019-10-21 23:26:27
【问题描述】：

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

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 中使用更典型的方法进行矩阵运算似乎是不可能的，因为我必须返回索引位置和最小距离。我已经尝试了几种不同的方法来为这项任务采用 CUarray，但到目前为止，在没有使用实际 for 循环的情况下，它们都失败了。此外，由于将距离矩阵存储在内存中，许多实现它的方法似乎效率极低，对于我的特定数据集，这很快超过了 128gb 的内存。

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

更新：

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

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 失败。我试图计算最小距离并忘记了索引。这对我的应用程序来说并不理想，但我可以接受它。但是，这只有在第一个数组只有一行时才能正常工作。我试图通过使用 maplices 来解决这个问题，但它不起作用。

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 版本和串行版本之间进行了一些基准测试。编辑：在扩展到 1000 亿次比较之后，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

【问题讨论】：

@kevbonham 我已经查看了 Distances.jl。我没有使用它的原因是缺乏并行支持。鉴于我必须分析多少数据，使用串行方法是不现实的。此外，据我所见，他们倾向于返回所有数据而不是丢弃不需要的距离，这导致我的特定情况下内存分配显着增加。返回整个距离矩阵对我来说简直太疯狂了。
您是否考虑过构建 NNTrees 并使用 NearestNeighbors.jl 中的 knn？我认为性能很难被击败。
NearestNeighbors 包会产生索引和精确距离 - 您不是将其与不同的算法混为一谈吗？
好吧，我绝对错了！ NearestNeighbor 软件包即使在扩展到数十亿次比较并且速度提高约 10,000 倍时也产生了相同的结果。我一直在做这一切都错了。
你能获得如此巨大的加速真是太好了！我想知道是否有办法以有用且可回答的方式重新表述这个问题。拥有所有更新是很好的文档，但确实取决于问题和答案。这有点像一段旅程（也许在 discourse.julialang.org 上重现会很好？

标签： julia gpu julia-gpu

【解决方案1】：

我不确定它是否会对您的情况有所帮助，但是在获取切片时 m1[k,:] 默认情况下，julia 会复制该内存（尽管这可能取决于 knn 函数对该切片的作用。

如果你把它改成knn(kdtree, @view m1[k,:], 1)，有什么改善

【讨论】：