Comments 50

Так изначальный код на питоне считал в один поток, я правильно понял листинг?

Так как основное тут матричное умножение, а 20 Гига сравнений могут выполниться на CPU быстро, можно просто сравнить Гигафлопсы на рынке и получить разницу 10-100 раз между CPU/GPU.
Немного устаревший график
image

Можно, но 20 секунд и 2 дня — это не 100 раз.
Даже не 1000, а ближе к 10k, что намекает, что не в CPU дело.

Мне в аналогичой ситуации очень помог HNSW (было на Хабре: https://habr.com/en/company/mailru/blog/338360/), позиций было порядка 10^6, около 5000 признаков. Индекс строился около суток, но поиск потом был очень быстрым (и, более того, настраиваемо быстрым в обмен на качество). Формально корректные результаты, правда, были так себе, потому что исходные данные были так себе.

Есть еще библиотека annoy, которая строит индекс и позволяет быстро находить ближайших соседей с разными метриками. Работает с размерностью вектора до 1000, желательно чтобы индекс помещался в оперативку. Естественно, поиск соседей приближенный, это не прямой перебор, используется лес из kd-деревьев.
Мда, функция top_3_similar_items(X,mindex) сходу ускоряется раза в 3. См. строчку, где вычисляется tmp. Ну нафига каждый раз длины векторов пересчитывать? Кстати, та же фигня и в последнем варианте кода.

Там много такого. Он бы мог сразу отнормировать свои вектора вместо того, чтобы в цикле каждый раз нормировать. А по нормированным векторам можно сразу вычислять матрицу косинусных дистанций хорошо оптимизированным матричным умножением.
Но сам посыл статьи эти ляпы не отменяют.


Прочитал ваш комментарий еще раз — вы как раз и имели ввиду под пересчетом длины векторов нормировку.

Матрицу — фиг, слишком большой размер. Даже если половину хранить. Но, скорее всего, на симметрии все равно можно двойку по времени выиграть. Только придется промежуточные результаты для всех товаров хранить и корректировать по ходу дела.

Да, матрица с результатом целиком не влезет в память.
Но у меня была похожая задача, я кусками по 5000 считал, отбирал каждый раз 1000 наиболее похожих. Наверное, тоже не самый оптимальный подход, но по сравнению с алгоритмом из этой статьи раз в 20 быстрее.

Неплохая статья как вводная.
А карточка уже была своя или ее специально для задачи фирма купила или в облаке использовали? Ценник не так что бы дружелюбен для «купить что бы попробовать».
А карточка уже была своя или ее специально для задачи фирма купила или в облаке использовали?
Так в самом конце статьи написано?

Этими алгоритмами можно пользоваться только если функция расстояния является метрикой, по крайней мере должно выполняться неравенство треугольника. Скалярное произведение (делённое на длины векторов) этому не удовлетворяет, например на двумерных векторах (0, 1), (0.5, 1), (0.5, 0). Если бы автор взял функцию, являющуюся метрикой, то действительно можно было бы обойтись без GPU.

Если векторы нормировать, их скалярное произведение монотонно зависит от квадрата евклидового расстояния между ними, а значит, задача нахождения ближайших соседей по скалярному произведению сводится к задаче нахождения ближайших соседей по евклидовой метрике.
Код был хорошо оптимизирован с применением операций и массивов Numpy. Отмечу, что все эти операции последовательно выполняются на CPU.

Скалярное произведение на 64 признака — numpy, двойной цикл по 105 индексов — на питоне, код хорошо оптимизирован, да.


Сдаётся, что с JIT-компиляцией без GPU вполне бы решалось всё за несколько десятков секунд, а автор в своём восхвалении либо лукавит, либо реально не понимает, в чём там дело.

У меня тоже после слов «мой код на питоне исполнялся долго, но я...» в голове зазвучала песенка «Write in C».
Ничего так ситуация, дали 10^5 строк данных, а еще совершенно случайно Volta V100 откуда-то взялась, можно сказать ждала своего часа.
Человек, собравший приложение на CUDA, с использованием CUDA Toolkit, передал это приложение мне. Могу ли я запустить его без установки CUDA Toolkit (по-сути SDK), имея лишь CUDA-совместимое железо, и, возможно, «драйвер» — либу CUDA?
Смотрите зависимости приложения. Для «чистой» CUDA должно быть достаточно драйвера видеокарты (nvidia-xxx) и cuda-runtime соответствующей версии.
Задача нахождения косинусного расстояния сводится к перемножению матриц. Имеется масса стандартных библиотек, которые решают эту операцию как на CPU так и на GPU.
И отличие между CPU и GPU при этом будет на более десяти раз.

К сожалению, матрица 10^5 на 10^5 чисел с плавающей запятой будет занимать около 80 ГБ, так что просто перемножить матрицы нельзя.

Это не так. Существуют методы, которые эффективно работают и с такими большими матрицами (пример реализации умножения матриц, в нем добавляется еще один уровень иерархии — и можно умножать матрицы превышающие размер оперативной памяти).

И сколько же тогда займёт перемножение, если данные будут постоянно дёргаться с диска?

Вся соль блочного подхода в том, что при правильном порядке умножения матриц (вычислительная сложность: O^3, зависимость по данным: O^2) скорость всегда упирается в вычислительные ресурсы. При этом последовательно задействуются: регистры, кэш процессора 1 уровня, 2-го уровня, 3-его, оперативная память, SSD, HDD и т.п.

Можете подсказать библиотеки (желательно для Пайтона), которые реализовывают подобное умножение для стогигабайтных матриц через диск? Потому что писать это всё руками, сами понимаете, такой себе вариант.

На одной моей прошлой работе были доступны сервера с 384-512 гигами памяти, и это было году в 2016-м.

"Программист на Фортране на любом языке может писать на Фортране".


У меня максимально "влобным" и неэффективным подходом без JITa ушло 20 сек на 10^4 элементов (то есть где-то в 60 раз быстрее, чем у автора). Цикл for в distance_func тривиально параллелизируется, то есть если есть CPU на ~30 потоков (или какой там Тредриппер можно купить за цену одной Вольты), то никаких "двух суток" и близко не нужно.


Заголовок спойлера
items = np.random.randint(1, 1000, size=(10**4, 64))

# normalize
items_norm = items / np.sqrt(np.sum(items**2, axis=1, keepdims=True))

def three_closest(X, item) -> (int, int, int):
    distances = np.sum(X[item, :] * X, axis=1)

    first = np.argmax(distances)
    distances[first] = 0
    second = np.argmax(distances)
    distances[second] = 0
    third = np.argmax(distances)
    distances[third] = 0

    return first, second, third

def distance_func(X):
    all_items = {}

    for i in range(len(X)):        
        all_items[i] = three_closest(X, i)

    return all_items

distance_func(items_norm)

Прошу прощения, прошляпил. В three_closest переменная first будет всегда иметь значение item, исправляется добавлением distances[item] = 0 перед первым argmax.

Тоже удивлен двумя сутками. Мы получали ускорение где-то в 100+ раз на Titan Black по сравнению с однопоточным режимом на CPU (i7, чистый Си без ++). Но это был расчет интегралов, а интегралы очень хорошо параллелятся.

Показательная статья как делать не надо.


Код мало того что не "хорошо оптимизирован", он вообще написан ужасно с точки зрения производительности. Одно и тоже вычисляется в цикле несколько раз и т.д. и т.п.


Не нужно использовать циклы для работы с матрицами. Все современные процессоры давно поддерживают Advanced Vector Extensions, которые позволяют существенно ускорить работу с матрицами.


Если уж в статье упоминается про NumPy, то именно её и надо было использовать для вычислений. Под капотом она использует высокооптимизированный код, написанный на C, и, насколько я знаю, используется тот самый AVX. Все можно было свести к паре строчек кода и работало бы оно максимум раз в 10-15 медленее чем на видеокарте.


Ну и как написали выше, есть более быстрые алгоритмы.

Есть библиотека rapids от Nvidia: rapids.ai которая реализует на GPU практически полный интерфейс pandas и, как бонус, поиск ближайших соседей в памяти GPU.
Считается именно 20 секунд, а не 0.5 секунд (на самом деле 0.5 мс будет выведено, перезапустил код локально). Автор посчитал затраченное время на асинхронный вызов.
Также хочется отметить, число колонок 59 — очень неудачное. Если добавить ещё 5 (например, нулевых) до 64, то выполняться на gpu будет быстрее примерно в два раза.

И, конечно, стоило бы воспользоваться готовой библиотекой для этой задачи. Например, при помощи faiss те же 10^5 векторов можно обработать за 25 секунд на TR 1920x (против 50 секунд для кода из статьи на 1080Ti). Результаты совпадают.
Код
import numpy as np
import faiss


n = 10**5
k = 59
batch_size = 18

items = np.random.rand(n, k).astype(np.float32)
items = items / np.sqrt(np.sum(items**2, axis=1, keepdims=True))
bounds = np.linspace(0, n, n // batch_size + 1, dtype=np.int64)

index = faiss.IndexFlatIP(59)
index.add(items)

distance = np.zeros((n, 4), dtype=np.float32)
index_val = np.zeros((n, 4), dtype=np.int64)
for lower, upper in zip(bounds[:-1], bounds[1:]):
    distance[lower:upper], index_val[lower:upper] = index.search(items[lower:upper], 4)
assert np.sum(index_val[:, 0] != np.arange(n)) == 0
distance = distance[:, 1:]
index_val = index_val[:, 1:]

Забавно, но перенесённой на Torch и немного оптимизированной версии решения от Physmatik уже требуется всего 0.5 секунды на GPU (1080Ti) и 10.5 секунд на CPU (TR 1920X) на обработку тех же 10^5 векторов. Учитывая возможность эффективного использования FP16 на карте V100, то на ней будет ещё быстрее.
Код на Torch
import torch

n = 10**5
k = 59
batch_size = 10**3
device = 'cuda:0' # 'cpu'
nn_number = 3

items = np.random.rand(n, k).astype(np.float32)
items = items / np.sqrt(np.sum(items**2, axis=1, keepdims=True))

def three_closest(X, indices):
    sample_n = indices.shape[0]
    index = torch.zeros((sample_n, nn_number), dtype=torch.int64, device=device)
    distances = X[indices, :] @ X.T
    ind_range = torch.arange(sample_n, device=device)
    distances[ind_range, indices] = 0
    for j in range(nn_number):
        max_index = torch.argmax(distances, dim=1)
        index[:, j] = max_index
        distances[ind_range, max_index] = 0
    return index

def distance_func(X):
    X = torch.from_numpy(X).to(device=device)
    all_items = torch.zeros((n, nn_number), dtype=torch.int64, device=device)
    for indices in torch.split(torch.arange(n, device=device), batch_size):        
        all_items[indices] = three_closest(X, indices)
    return all_items.cpu().numpy()

nearest_torch = distance_func(items)

Надеюсь в скором будущем поддержка вычислений на GPU будет добалена в ОС.
А зачем CUDA операционке? Что ей нужно такое большое обсчитывать?
Казалось бы. Но ускорение 2D происходит, насколько я понимаю, как одной из поверхностей 3D, т.е. чистых 2D ускорителей не делают. А работать с отключённым ускорением (например, при отсутствии драйвера видеокарты, т.е. «стандартный адаптер vga») тяжело.
Возможно так и происходит, но тем не менее это означает лишь что ей нужен 3D-ускоритель, а не 3D вообще :)
Хотя вроде в видеокартах с развитием 3D не отказывались и от 2D-операций…
Точно не знаю, но как-то нерационально иметь железо и не использовать его. Будет возможность-наверняка и напишут софт для этого. Хотя бы джипеги жать:)
Возможность-то уже давно есть, CUDA не вчера появилась :) А вот необходимости в ней у ОС нет :)
10 секунд переносим данные к gpu, полсекунды считаем, 10 секунд переносим данные обратно.
нвидия, это же какой-то каменный век.
Это не нвидия каменный век, а нужно понимать что именно делается инструментом. Там выше уже написали про асинхронность вызовов, по сути вычисления лениво выполняются только тогда, когда хост их забирает, т.е. 20 секунд и есть время выполнения.

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

У вольты вполне себе нвлинк может быть. Уже менее тормозной. А что делать? Это фундаментальная проблема — чем дальше данные перемещать, тем медленнее это будет (и дороже с точки зрения энергии).
Only those users with full accounts are able to leave comments. Log in, please.