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

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

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)

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

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

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 и серийной версией. РЕДАКТИРОВАНИЕ: после масштабирования до 100 миллиардов сравнений версия 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
Источник

1 ответ

Я не уверен, что это поможет в вашем случае, но когда вы берете кусочек m1[k,:]по умолчанию julia копирует эту память (хотя, возможно, это зависит от того, что knnфункция делает с этим фрагментом.

Улучшится ли что-нибудь, если вы измените его на knn(kdtree, @view m1[k,:], 1)

Источник
Другие вопросы по тегам