Джулия: БЫСТРЫЙ способ расчета наименьших расстояний между двумя наборами точек
У меня есть 5000 3D-точек в матрице A и еще 5000 3D-точек в матрице B.
Для каждой точки в A я хочу найти наименьшее расстояние до точки в B. Эти расстояния должны храниться в массиве из 5000 записей.
Пока у меня есть это решение, работающее примерно 0.145342 seconds (23 allocations: 191.079 MiB). Как я могу улучшить это дальше?
using Distances
A = rand(5000, 3)
B = rand(5000, 3)
mis = @time minimum(Distances.pairwise(SqEuclidean(), A, B, dims=1), dims=2)
Ответы 2
Генерация матрицы больших расстояний с использованием Distances.pairwise(SqEuclidean(), A, B, dims=1) неэффективна, потому что основная память в настоящее время довольно медленная по сравнению с кешем ЦП и вычислительной мощностью современных ЦП, и в ближайшее время ситуация не улучшится (см. «стену памяти»). Быстрее вычислить минимум на лету, используя два базовых вложенных цикла for. Кроме того, можно использовать несколько ядер для более быстрого вычисления с использованием нескольких потоков.
function computeMinDist(A, B)
n, m = size(A, 1), size(B, 1)
result = zeros(n)
Threads.@threads for i = 1:n
minSqDist = Inf
@inbounds for j = 1:m
dx = A[i,1] - B[j,1]
dy = A[i,2] - B[j,2]
dz = A[i,3] - B[j,3]
sqDist = dx*dx + dy*dy + dz*dz
if sqDist < minSqDist
minSqDist = sqDist
end
end
result[i] = minSqDist
end
return result
end
mis = @time computeMinDist(A, B)
Обратите внимание, что интерпретатор Julia по умолчанию использует 1 поток, но это можно настроить с помощью переменной среды JULIA_NUM_THREADS=auto или просто запустив его с помощью флага --threads=auto. См. документацию по многопоточности для получения дополнительной информации.
Результаты производительности
Вот результаты производительности на моей машине i5-9600KF с 6 ядрами (с двумя матрицами 5000x3):
Initial implementation: 93.4 ms
This implementation: 4.4 ms
Таким образом, эта реализация в 21 раз быстрее. Результаты одинаковы для нескольких ULP.
Обратите внимание, что код, безусловно, может быть дополнительно оптимизирован с помощью разбиения циклов и, возможно, путем перестановки A и B, чтобы JIT мог генерировать более эффективную реализацию с использованием SIMD-инструкций.
@Shayan Хорошая идея. Я не знал этого обозначения. Я узнаю Джулию, отвечая на подобные вопросы :) . Спасибо.
Спасибо! Ваш ответ можно изменить, добавив @inbounds for j = 1:m ...
@Benvorth Действительно, это почти в два раза быстрее :)! И тебе спасибо.
Это стандартный способ сделать это, поскольку он будет иметь лучшую временную сложность (особенно для больших данных):
using NearestNeighbors
nn(KDTree(B'; leafsize = 10), A')[2] .^ 2
Два комментария:
- по умолчанию вычисляется евклидово расстояние (поэтому я возводил его в квадрат)
- по умолчанию NearestNeigbors.jl предполагает, что наблюдения хранятся в столбцах (поэтому мне нужны
B'иA'в решении; если бы ваши исходные данные были перенесены, это не понадобилось бы; причина, по которой он разработан таким образом, заключается в том, что Джулия использует хранение основной матрицы столбцов )
О, ничего себе, это решение примерно на 50% быстрее, чем (уже БЫСТРО), предоставленное @Jérôme Richard
Хороший. Это действительно интересное решение для больших матриц. Это немного быстрее, чем два основных вложенных цикла на входе OP. С оптимизацией SIMD/кеша 2 цикла должны работать быстрее. При этом я думаю, что это решение использует только 1 поток, и оно также имеет преимущество в простоте.
Вы можете легко сделать его многопоточным, разбив матрицу A на куски. Однако реальная разница будет заключаться в том, что каждая из этих матриц будет иметь примерно 1 миллион строк. Тогда подход kdtree занимает 1,3 секунды (без многопоточности), а другое решение займет не менее нескольких минут (быстрая оценка на моем компьютере — 30 минут).
Я думаю, было бы неплохо избежать скобок в size(A)[1], size(A)[1], просто заменив его на size(A, 1), size(A, 1).