Минимальные расстояния среди евклидовой матрицы расстояний - PullRequest
Новая
       24

Минимальные расстояния среди евклидовой матрицы расстояний

10 голосов
/ 22 октября 2019

У меня есть некоторый код, который будет вычислять расстояния между каждой декартовой координатой в одной матрице и любой другой координатой в другой. Для каждой координаты будет возвращено минимальное расстояние вместе с позициями индекса для координат, которые дали минимум. 

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)

Хотя это работает достаточно хорошо с небольшими трехмерными облаками точек, я стремлюсь повысить производительностьдля больших облаков точек с помощью анализа на основе графического процессора. Тем не менее, использование более типичных способов выполнения матричных операций в Джулии кажется невозможным, так как мне приходится возвращать позиции индекса и минимальное расстояние. Я пробовал несколько разных способов использования CUarrays для этой задачи, но до сих пор все они потерпели неудачу без использования фактических циклов for. Кроме того, многие способы его реализации кажутся исключительно неэффективными из-за хранения матрицы расстояний в памяти, которая быстро превышает 128 ГБ ОЗУ для моего конкретного набора данных.

Может кто-нибудь помочь мне с тем, как правильно реализовать это в Джулии для работы на GPU? Является ли CUarrays даже правильным подходом, или это слишком абстрактный уровень, учитывая, что я возвращаю индексы в дополнение к расстоянию? Я пытался вычислить норму 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)

ОБНОВЛЕНИЕ:

Некоторый сравнительный анализ между моей предыдущей распределенной версией, модифицированной распределенной версией, версией графического процессора и последовательной версией. EDIT: после масштабирования до 100 миллиардов сравнений, версия графического процессорабольше не превосходит мою предыдущую распределенную версию ... Любые мысли о том, почему это ???? 

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
...