Блог

У меня есть матрица 16000 * 16000, и я хочу найти минимальную запись. Эта матрица является матрицей расстояний, поэтому она симметрична относительно диагонали. Чтобы получить ровно один минимум за каждый раз, я установил нижний треугольник и диагональ np.inf равными . Ниже приведен пример матрицы 5 * 5:

 inf a0  a1  a2  a3
inf inf a4  a5  a6
inf inf inf a7  a8
inf inf inf inf a9
inf inf inf inf inf
 

Я хочу найти индекс минимальной записи только в верхнем треугольнике. Однако, когда я использую np.argmin() , он все равно будет проходить через всю матрицу. Есть ли какой-нибудь способ «игнорировать» нижний треугольник и увеличить скорость?

Я перепробовал много методов, таких как:

1. Используйте массив с маской
2. Используется triu_indices() для извлечения верхнего треугольника, а затем для нахождения минимального
3. Установите для записей в нижнем треугольнике и диагонали значение None вместо np.inf , затем используйте np.nanargmin() для поиска минимального

Однако все методы, которые я пробовал, медленнее np.argmin() , чем прямое использование.

Спасибо за ваше время, я был бы признателен, если бы вы могли мне помочь.

ОБНОВЛЕНИЕ 1: Некоторая предыстория моей проблемы

Фактически, я внедряю модифицированную версию агломеративной кластеризации с нуля. Исходный набор данных равен 16000 * 64 (у меня 16000 точек, каждая из которых 64-мерная). Сначала я создаю 16000 кластеров, и каждый содержит ровно одну точку. На каждой итерации я нахожу ближайшие 2 кластера и объединяю их, пока не будет выполнено условие завершения.

Чтобы избежать повторного вычисления расстояний, я сохраняю расстояния в матрице расстояний 16000 * 16000. Я установил диагональный и нижний треугольник np.inf равными . На каждой итерации я найду наименьшую запись в матрице расстояний, а индекс этой записи соответствует 2 ближайшим кластерам, скажем c_i , и c_j . После этого в матрице расстояний я заполняю 2 строки и 2 столбца, соответствующие c_i и c_j np.inf, что означает, что эти 2 кластера объединены и больше не существуют. Затем я вычислю массив расстояний между новым кластером и всеми другими кластерами, затем помещу массив в 1 строку и 1 столбец, соответствующие c_i .

Позвольте мне пояснить: во всем процессе размер матрицы расстояний никогда не меняется. На каждой итерации, поскольку 2 строки и 2 столбца соответствуют 2 ближайшим кластерам, которые я нашел, я заполняю 1 строку и 1 столбец np.inf и помещаю массив расстояний нового кластера в другие 1 строку и 1 столбец.

Теперь узким местом производительности является нахождение наименьшей записи в матрице расстояний, которая занимает 0,008 с. Время выполнения всего алгоритма составляет около 40 минут.

ОБНОВЛЕНИЕ 2: как я вычисляю матрицу расстояний

Ниже приведен код, который я использовал для генерации матрицы расстояний:

 from sklearn.metrics import pairwise_distances

dis_matrix = pairwise_distances(dataset)

for i in range(num_dim):
    for j in range(num_dim):
        if i >= j or (cluster_list[i].contain_reference_point and cluster_list[j].contain_reference_point):
            dis_matrix[i][j] = np.inf
 

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

Если мы создадим резервную копию шага, предполагая, что матрица расстояний симметрична и основана на (i, n) фигурном массиве с i точками в n измерениях, а метрика расстояния декартова, это можно сделать очень эффективно со KDTree структурой данных:

 i = 16000
n = 3
points = np.random.rand(i, n) * 100

from scipy.spatial import cKDTree
tree = cKDTree(points)
close = tree.sparse_distance_matrix(tree, 
                                    max_distance = 1, #can tune for your application
                                    output_type  = "coo_matrix") 
close.eliminate_zeros()
ix = close.data.argmin()
i, j = (close.row[ix], close.col[ix])
 

Это довольно быстро, но это зависит от вашего приложения и метрики расстояния, если это полезно для вас.

Если вам вообще не нужна матрица расстояний (и нужны только индексы), вы можете сделать:

 d, ix = tree.query(points, 2)
j, i = ix[d[:, 1].argmin()]
 

РЕДАКТИРОВАТЬ: это плохо работает для данных с высокой размерностью. Поскольку вы столкнулись с проклятием размерности, вам, вероятно, потребуется грубая сила. Я рекомендую scipy.spatial.distance.pdist для этого:

 from scipy.spatial.distance import pdist
D = pdist(points, metric = 'seuclidean')  # this only returns the upper diagonal
ix = np.argmin(D)

def ix_to_ij(ix, n):
    sorter = np.arange(n-1)[::-1].cumsum()
    j = np.searchsorted(sorter, ix)
    i = ix - sorter[j]
    return i, j

ix_to_ij(ix, 16000)
 

Не полностью протестировано, но я думаю, что это должно сработать.

Одна вещь, о которой я могу подумать, которая может дать вам толчок, — это использование numba.njit :

 @njit
def upper_min(m):
    x = np.inf
    for r in range(0, m.shape[0] - 1):
        for c in range(r   1, m.shape[1]   1):
            if x < m[r, c]:
                x = m[r, c]
 

Не указывайте время при первом запуске. Компиляция идет медленно.

Другим способом может быть использование разреженных матриц каким-либо образом.

Вы можете выбрать верхний треугольник массива путем маскирования, простой пример:

 import numpy as np
arr = np.array([[0, 1], [2, 3]])
# Mask of upper triangle
mask = np.array([[True, True],[False, True]])
# Masking returns only upper triangle as 1D array
min_val = np.min(arr[mask]) # Equal to np.min([0, 1, 3])
 

Таким образом, вместо того, чтобы создавать нижний треугольник as inf , вам нужно сгенерировать маску, в которой находится нижний треугольник False , а верхний треугольник True , и применить маскировку arr[mask] , которая возвращает 1D массив верхнего треугольника, затем вы применяете min