Check Levi9 best QA positions to Backbase team!
×Закрыть

Оптимізація коду Python

Є фрагмент коду, де обчислюється попарна евклідова відстань в D-мірному просторі між 10-ма точками. Координати
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)

Як можна пришвидшити цей обрахунок?

👍НравитсяПонравилось0
В избранноеВ избранном0
LinkedIn

Похожие топики

Допустимые теги: blockquote, a, pre, code, ul, ol, li, b, i, del.
Ctrl + Enter
Допустимые теги: blockquote, a, pre, code, ul, ol, li, b, i, del.
Ctrl + Enter

как вариант, кстати, если все данные не получается за раз засунуть в один вызов pdist(), можно соорудить minibatch из pdist() и cdist().

А какая общая задача? Может можно обойтись и без попарных дистанций.

Насправді треба знайти, всі попарні відстані коротші певної заданої довжини.

простая оптимизация: отсортируй точки по одной из координат, сразу же можно будет скипать пары, у которых дистанция больше

en.wikipedia.org/...st_pair_of_points_problem

Для цього методу написано, 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..

По каждой размерности (или по 1-2) можно массово отфильтровывать явно неподходящие пары и не считать их дистанции. Если граничная дистанция 2, а разница по координате 4, то эта пара явно не подходит.

en.wikipedia.org/wiki/K-d_tree — тоже прикольная штука, но для 10 точек в 20-мерном пространстве выигрыша не даст.

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?
int має межі від −2**31 до 2**31-1, або від 0 до 2**32-1 (але в python нема unsigned int)
double має набагато більші межі, точно не скажу, але він принципово заточений на величезні числа, правда, за рахунок втрати точності «знизу» — my_very_long_double + 1 == my_very_long_double
.
тому, якщо ти переведеш double -> int, то для кожної відстані > 2**31-1 python переведе число в long, який принципово безрозмірний, але займає негарантовану кількість пам’яті (явно більше 4 байт) — ось тут описано stackoverflow.com/...um-value-for-long-integer
.
або якщо це якась ліба, яка працює іменно з int 4, то ти отримаєш старе добре переповнення — всі числа більші за 2**31-1, стануть від’ємними або обріжуться до 32 біт

Відстані можна масштабувати, діапазон для задачі від 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 ближайших точек? тогда можно сделать что-то вроде

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. Я використовую
print memory_usage_psutil()
Для dtype=np.float16 дає 182.91 МВ, для float64 448.69 МВ. Основний об’єм пам’яті займає distances, чому ви міряєте для points ?

ну, 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() и так написан на С, и там всегда возвращается 64-bit double. так что питон тут ни при чем — это проблема конкретной библиотеки, которая легко обходится. в моем примере выше можно обрабатывать исходные данные по кускам и результаты сразу преобразовывать в нужный формат

Ніби 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 процентов и вы сразу получаете корень а не сумму квардратов. а кто это пандас ? на С++ можна распараллелить на количество ядер процессора, но на питоне я такое не разу не делал.

дык где numpy, там и pandas.. а у вас в векторе есть еще какие-то фичи? если нет, то используйте numpy и не парьтесь

Какой версии Python?

В таком случае, попробуйте заменить range на xrange.

Після deque і xrange рахує на 8% довше.

range(0, 9)

у Вас вроде как 9 точек перебирает (0..8) — если память не изменяет.
А с точки зрения скорости вначале вместо list используйте deque — подозреваю, что сейчас самая большая задержка при добавлении в конец массива. А дальше по идее можно попробовать лишние квадраты не считать, но не думаю, что тут это критично.

У Вас в самом деле 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% довше.

Подписаться на комментарии