Є фрагмент коду, де обчислюється попарна евклідова відстань в
V_pd[k][i], де k — компонета координати x,y,z,.. i — номер точки.
dist = []
for i in range(0, 10):
for j in range(i+1, 10):
su=0
for k in range(0, D):
mih=V_pd[k][i]-V_pd[k][j]
su +=mih*mih
dist.append(su)
Як можна пришвидшити цей обрахунок?
как вариант, кстати, если все данные не получается за раз засунуть в один вызов pdist(), можно соорудить minibatch из pdist() и cdist().
А какая общая задача? Может можно обойтись и без попарных дистанций.
Насправді треба знайти, всі попарні відстані коротші певної заданої довжини.
простая оптимизация: отсортируй точки по одной из координат, сразу же можно будет скипать пары, у которых дистанция больше
Для цього методу написано, These approaches are not efficient for dimension d>2, в мене d>20
Полная фраза: Computing either the Delaunay triangulation or the Voronoi diagram takes O(n log n) time. These approaches are not efficient for dimension d>2, while the divide-and-conquer algorithm can be generalized to take O(n log n) time for any constant value of d..
По каждой размерности (или по
en.wikipedia.org/wiki/K-d_tree — тоже прикольная штука, но для 10 точек в
pdist — это очень часто необходимость и в данном случае разумнее заюзать какую либу для математики. Тот же MKL, но я питона настолько не знаю, возможно numpy, scipy это уже и делают.
Или просто реализовать в С, С++ и прилинковать.
Ну мало ли, вдруг можно на кластеры побить и работать как с единым объектом, либо отбрасывать их пачками. Вон даже N-Body задачу научились за O(n log n) решать вместо квадрата.
всі попарні відстані коротші певної заданої довжини.— хороший намек. Для десяти точек можно и в лоб, если больше, то Quad, Oct и так далее tree.
Вепроятнее всего — это маленькая часть какой-то большой задачи. Просто это «бутылочное горлышко» в его случае. И самый простой способ — это заюзать готовую либу, что считает pdist максимально быстро на его железе (кластере).
Например в матлаб такой код тоже долго будет считься. Но можно тупо зарепматить (repmat) и сделать все вычисления в виде (A-B).^2, это если памяти хватает. Или вообще сделать блочный алгоритм и в параллель на нескольких ядрах запустить (в кластере). Путей море.
Можно матлабовский pdist заюзать.
Но выбор решения зависит от конкретной задачи. У него маленький массив, так что я бы просто тупо repmat заюзал, если бы родной pdist не удовлетворил.
scipy.spatial.distance.pdist обчислює дуже швидко, проблема в пам’яті. Наприклад, як заставити його зберігати в пам’яті відстані в int а не в float? Або зберігати лише відстані менші граничної? Я не знайшов в параметрах.
Або зберігати лише відстані менші граничної?По определению он этого делать не будет. Тут уже самому ручками на сях писать надо.
зберігати в пам’яті відстані в int а не в float?int нынче 4 или 8 байт, float-4, double 8 или 10. Но переход в int может просто повысить скорость вычислений (в теории), но ты наткнешьс яна переполнение, которое надо будет учитывать. Но для float и double есть вся мощь SSE (всех).
Но переход в int может просто повысить скорость вычислений (в теории), но ты наткнешьс яна переполнениеНе зрозумів, якщо я редукую з double 8 на int 4, я отримую виграв по пам’яті в двічі, те що треба. Швидкість в принципі значення не має.
Это я не знаю, можно-ли питон считать в даблах или вообще с ними работать. Матлаб нельзя, точнее можно, но считать он будет всё одно в даблах. Это уже возможности С и С++. Да напиши ты либу на сях и заюзай ее в Питон, в ней ты уже сможешь всё что хочешь делать..
Не зрозумів, якщо я редукую з double 8 на int 4, я отримую виграв по пам’яті в двічіа ти точно data scientist?
Відстані можна масштабувати, діапазон для задачі від 0 до 10000 більш ніж достатній, так що за межі 2**31-1 не вийде.
Интересно, а чем обусловлен отступ во второй строке ?
for i in range(0, 10):
нічим, його не повинно бути, опечатка
Не підходить, бо ціною високої продуктивності та використання JIT-компіляції є більше споживання пам’яті. А пам’яті обмаль.
нада перевести координаты точек в вид
a = np.array((2,2,2,2,2))
b = np.array((4,4,4,4,4))
тогда можна делать так tmp = a — b
ну и решение переводится в вид
dist = []
for i in range(0, points):
for j in range(i+1, points):
su = np.linalg.norm(V_pd1[i]-V_pd1[j])
dist.append(su)
su это уже корень из суммы квадратов
+1 за numpy и по скорости, и по памяти. можно еще использовать numpy.float16 или какой-нибудь int для экономии памяти.
а, кстати, точно нужно вычислять все попарные расстояния, или достаточно, скажем, найти 10 ближайших точек? тогда можно сделать что-то вроде
Sergiy Matusevych
Alex Koshterek
Alex Koshterek
Viktor Chyzhdzenka
Andy Wilk
ScorpZ Lapshov
anonymous
Ванкель Ванкелевич
mobile mobile
Sergiy Matusevych
import numpy as np
import heapq
arr_size = (10000, 1000)
points = np.zeros(arr_size, dtype=np.float16)
points[:] = np.random.standard_normal(arr_size)
def closest_n(n, points):
return heapq.nsmallest(n, (
(np.linalg.norm(points[i] - points[j]), i, j)
for i in xrange(len(points))
for j in xrange(i + 1, len(points))))
подойдет такое?
о, нашел. все, что вам нужно — scipy.spatial.distance.pdist():
import numpy as np
import scipy.spatial.distance as dist
points = np.random.randn(10000, 20)
distances = dist.pdist(points, metric='euclidean')
считает все расстояния за 1.3 сек на моем лаптопе. batteries included, что называется :)
В точку, те що шукав, величезне спасибі.
Для великих матриць dist.pdist жере дуже багато пам’яті. Можливо знаєте як відстані зберігати в int а не в дефолтному double?
Сорри я сейчас без компа, но попробуй сделать dtype=float16 для входных данных, и/или оставлять только n наименьших расстояний — и то, и другое есть в моём первом примере
Я робив dtype=float16, але об’єм використаної памяті той самий.
В принципі зменшення необхідної пам’яті в двічі, вирішило би проблему. Мені достатньо int діапазону від 0 до 10000 вище голови.
ну, по крайней мере, для входных данных можно объем памяти уменьшить:
arr_size = (10000, 20)
points = np.zeros(arr_size, dtype=np.float16)
points[:] = np.random.standard_normal(arr_size)
print sys.getsizeof(points)
400112

points = np.random.standard_normal(arr_size) # dtype=float64
print sys.getsizeof(points)
1600112

впрочем, 1.6M против 400K — не думаю, что это может вызвать какие-либо проблемы.
А pdist(), похоже, всегда возвращает float64:
distances = dist.pdist(points)
print sys.getsizeof(points)
399960096
результат можно сжать во float16, используя тот же прием:
n = len(points)
distances = np.zeros(int(n * (n - 1) / 2), dtype=np.float16)
distances[:] = dist.pdist(points)
print sys.getsizeof(points)
99990096
объем памяти, кстати, не такие уж большой — ~400M максимум. может, у вас там какие-то лишние временные объекты память отъедают?
Мені
print sys.getsizeof(points)завжди повертає число 80. Я використовую
ну, points я мерял просто за компанию. понятно, что distances ~ O(points^2). кстати, странно, что у вас sys.getsizeof() дает 80 — попробуйте тогда points.nbytes, это точно должно работать.
в общем, если памяти все равно не хватает, то minibatch вам в помощь. разбиваем points на несколько сегментов, и для каждого считаем pdist() между точками внутри сегмента, и cdist() между точками этого сегмента и остальных. будет, конечно, не так быстро, как один pdist(), но наверняка намного быстрее обработки каждой пары отдельно.
а, сейчас напишу :) там все просто должно быть, главное, с индексами не ошибиться — ну да вы и потестируете, если что :)
def minibatch_distances(points, chunk_size):
for i in xrange(0, len(points), chunk_size):
chunk1 = points[i:i + chunk_size, :]
d = dist.pdist(chunk1)
k = 0
for m in xrange(len(chunk1) - 1):
for n in xrange(m + 1, len(chunk1)):
yield (d[k], m + i, n + i)
k += 1
for j in xrange(i + chunk_size, len(points), chunk_size):
chunk2 = points[j:j + chunk_size, :]
for (m, d) in enumerate(dist.cdist(chunk1, chunk2)):
for n in xrange(len(chunk2)):
yield (d[n], m + i, n + j)
якось так..
Дякую, так я розбиваю матрицю на 3×3 частини, це допомагає економити пам’ять. В мене питання, чому просунутий Python не дозволяє редукувати розмір змінних, наприклад мені би int16 вистачило з головою?
Потому что это не С и не С++. С такими требованиями пишешь функцимю на С или С++, делаешь либу и юзаешь из Питона.
Ви пропонуєте pdist написати на С/С++?
Давно уже, кроме того можно просто заюзать ATLAS, MKL и т.п.
pdist для pd.DataFrame на С/С++ не осилю.
Ну, тут я тебе ничем помочь не могу. В другое время написал бы для тебя за деньги, но сейчас сильно занят. выложи задачу на какой фрилансерской бирже, тебе сделают.
pdist() и так написан на С, и там всегда возвращается
Ніби heapq входить в Python Standard Library, але помилка:
NameError: global name ’heapq’ is not defined
D = 1000;
points = 100000;
у меня дает мемори эррор, сколько у вас точек ?
В мене d=20 і points = 10000, ерору пам’яті нема.
можна записать строку еще так , это дает гдето 15 процентов на моем компе
su += math.pow(V_pd[k][i]-V_pd[k][j], 2)
но там вверху я описал другое решение, оно дает гдето то 20%
на 20% малувато, а чи можна якось використати pandas ?
по сравнению с оригинальм кодом будет гдето 30 процентов и вы сразу получаете корень а не сумму квардратов. а кто это пандас ? на С++ можна распараллелить на количество ядер процессора, но на питоне я такое не разу не делал.
pandas-docs.github.io/pandas-docs-travis
Маю на увазі pandas Data Frames.
дык где numpy, там и pandas.. а у вас в векторе есть еще какие-то фичи? если нет, то используйте numpy и не парьтесь
В таком случае, попробуйте заменить range на xrange.
range(0, 9)
У Вас в самом деле 10 точек или массив все же «слегка» подлиннее будет?
Если сильно подлиннее я бы переставил местами индексы — так как сейчас координаты одной точки лежат не подряд в памяти — при больших массивах падает эффективность работы кэша.
Ну и тогда можно будет использовать sum, map и т.д. вместо циклов.
Насправді набагато більше ніж 10. Переставити місцями це
mih=V_pd[k][i]-V_pd[k][j] -----> mih=V_pd[i][k]-V_pd[j][k]
да, только и заполнить массив данными именно таким образом.
А то у Вас сейчас получается [ [x0,x1,x2...], [y0,y1,y2...], [z0,z1,z2...],...], а так как во внутреннем цикле Вы перебираете координаты точки, то логичней сделать так:
[[x0,y0,z0...], [x1,y,1,z1,...], [x2,y2,z2,...]...] — чтобы данные, который обрабатываются подряд — и в памяти находились подряд. Не знаю, как в питоне — но в том же С++ при этом компилятор задействует потоковые команды (SSEx и др), что значительно ускоряет вычисления. Ну и кэш помогает.
Спробую це, бо після
dequeі
xrangeрахує на 8% довше.
mobile mobile
Max M
Max M
Volodymyr Sydorenko
Volodymyr Sydorenko
58 комментариев
Добавить комментарий Подписаться на комментарииОтписаться от комментариев