Эффективно вычислять только уникальные элементы матрицы симметричного внешнего произведения - PullRequest
Новая
       5

Эффективно вычислять только уникальные элементы матрицы симметричного внешнего произведения

0 голосов
/ 04 февраля 2019

У меня есть большой набор из N d-мерных векторов (находящихся в матрице), которые я поднимаю, беря собственное внешнее произведение (т.е. каждый вектор i умножается на себя).Для каждого вектора получается симметричная матрица, в которой (d + 1) выбираются 2 уникальные записи.Для всех данных это тензор N xdxd.Я хотел бы вычислить только уникальные (d + 1), выбрать 2 записи из нижней диагонали каждого тензорного среза и сохранить их в векторе.Я хочу сделать это с минимальным использованием памяти и как можно быстрее в Python - включая использование привязок Си.

Если вы делаете это, используя стандартные методы numpy, он выделяет всю матрицу целиком.Это примерно вдвое увеличивает сложность памяти по сравнению с тем, что действительно требуется.

Для определения масштаба рассмотрим случай, когда N = 20k и d = 20k.Тогда N * d ^ 2 * ~ 8 байт на элемент = (2 * 10 ^ 4) ^ 3 * 8 байт = 64 терабайта.

Если мы вычисляем только векторы, которые кодируют уникальныйзаписей, у нас есть (20001 выбор 2) * 20k * 8 = 200010000 * 20000 * 8 байт = 32 терабайта.

Есть ли быстрый способ сделать это, не прибегая к медленным методам (таким как кодирование)мой собственный внешний продукт в python)?

Редактировать: я отмечу, что аналогичный вопрос был задан в Создать массив внешних продуктов в numpy

Я уже знаю, какчтобы вычислить это с помощью einsum (как в приведенном выше вопросе).Тем не менее, не было получено ответа об этом без дополнительных (d выбрать 2) вычислений и распределений

Редактировать 2: Этот поток Как использовать симметрию во внешнем продукте в Numpy (или других решениях Python)? задает связанный вопрос, но не затрагивает сложность памяти.Верхний ответ по-прежнему будет выделять массив adxd для каждого внешнего продукта.

Этот поток Numpy Performance - внешний продукт вектора с его транспонированием также учитывает вычислительные соображения собственного внешнего продукта, но делаетне достичь решения, эффективного для использования памяти.

Редактировать 3: Если кто-то хочет выделить весь массив, а затем извлечь элементы, np.tril_indices или scipy.spatial.distance.squareform добьются цели.

1 Ответ

0 голосов
/ 05 февраля 2019

Не совсем точно, как вы хотите, чтобы ваш вывод, но всегда есть возможность использовать Numba: 

import numpy as np
import numba as nb

# Computes unique pairwise products
@nb.njit(parallel=True)
def self_outer_unique(a):
    n, d = a.shape
    out = np.empty((n, (d * d + d) // 2), dtype=a.dtype)
    for i in nb.prange(n):
        for j1 in range(d):
            for j2 in range(j1, d):
                idx = j1 * (2 * d - j1 + 1) // 2 + j2 - j1
                out[i, idx] = a[i, j1] * a[i, j2]
    return out

Это даст вам массив со всеми уникальными продуктами в каждой строке (то есть сплющенный верхний треугольникполного вывода). 

import numpy as np

a = np.arange(12).reshape(4, 3)
print(a)
# [[ 0  1  2]
#  [ 3  4  5]
#  [ 6  7  8]
#  [ 9 10 11]]
print(self_outer_unique(a))
# [[  0   0   0   1   2   4]
#  [  9  12  15  16  20  25]
#  [ 36  42  48  49  56  64]
#  [ 81  90  99 100 110 121]]

С точки зрения производительности это быстрее, чем вычисление полного внешнего продукта с помощью NumPy, хотя воссоздание полного массива из этого занимает больше времени. 

import numpy as np

def np_self_outer(a):
    return a[:, :, np.newaxis] * a[:, np.newaxis, :]

def my_self_outer(a):
    b = self_outer_unique(a)
    n, d = a.shape
    b_full = np.zeros((n, d, d), dtype=a.dtype)
    idx0 = np.arange(n)[:, np.newaxis]
    idx1, idx2 = np.triu_indices(d)
    b_full[idx0, idx1, idx2] = b
    b_full += np.triu(b_full, 1).transpose(0, 2, 1)
    return b_full

n, d = 1000, 100
a = np.arange(n * d).reshape(n, d)
print(np.all(np_self_outer(a) == my_self_outer(a)))
# True

%timeit np_self_outer(a)
# 24.6 ms ± 248 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit self_outer_unique(a)
# 6.32 ms ± 69.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit my_self_outer(a)
# 124 ms ± 770 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
...