элементарное внешнее произведение с разреженными матрицами - PullRequest
Новая
       94

элементарное внешнее произведение с разреженными матрицами

3 голосов
/ 07 июня 2019

Я хочу сделать поэлементное внешнее произведение трех (или четырех) больших двумерных массивов в python (значения float32 округлены до 2 десятичных знаков).Все они имеют одинаковое количество строк «n», но разное количество столбцов «i», «j», «k».
Полученный массив должен иметь форму (n, i * j * k).Затем я хочу суммировать каждый столбец результата, чтобы получить массив 1D формы (i * j * k). 

np.shape(a) = (75466, 10)
np.shape(b) = (75466, 28)
np.shape(c) = (75466, 66)

np.shape(intermediate_result) = (75466, 18480)
np.shape(result) = (18480)

Благодаря ruankesi и divakar , яполучил фрагмент кода, который работает: 

# Multiply first two matrices
first_multi = a[...,None] * b[:,None]
# could use np.einsum('ij,ik->ijk',a,b), which is slightly faster
ab_fills = first_multi.reshape(a.shape[0], a.shape[1]*b.shape[1])

# Multiply the result with the third matrix
second_multi = ab_fills[..., None] * c[:,None]
abc_fills = second_multi.reshape(ab_fills.shape[0], ab_fills.shape[1] * c.shape[1])

# Get the result: sum columns and get a 1D array of length 10*28*66 = 18 480
result = np.sum(abc_fills, axis = 0)

Задача 1 : Производительность

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

Задача 2 : Мои матрицы редки

Действительно, например, «a» содержит 70% от 0.Я пытался играть с scipy csc_matrix, но действительно не мог получить рабочую версию.(чтобы получить поэлементный внешний продукт здесь, я иду через преобразование в трехмерную матрицу, которая не поддерживается в scipy sparse_matrix)

Проблема 3 : использование памяти

Если я пытаюсь также работать с 4-й матрицей, у меня возникают проблемы с памятью.


Я предполагаю, что преобразование этого кода в sparse_matrix сэкономит много памяти и ускорит вычисления, игнорируя многочисленные значения 0.Это правда?Если да, может ли кто-нибудь мне помочь?
Конечно, если у вас есть предложения по улучшению реализации, я также очень заинтересован.Мне не нужны никакие промежуточные результаты, только конечный результат 1D.
Уже несколько недель я застрял в этой части кода, я схожу с ума!

Спасибо!



Редактировать после ответа Дивакара

Подход № 1:
Очень хороший лайнер, но на удивление медленнее, чем оригинальный подход (?).
В моем тестовом наборе данных подход № 1 занимает 4,98 с ± 3,06 мс на цикл (без ускорения с optimize = True)
Первоначальный разложенный подход занимал 3,01 с ± 16,5 мс на цикл


Подход № 2:
Абсолютно здорово, спасибо!Какое впечатляющее ускорение!
62,6 мс ± 233 мкс на цикл


Что касается Numberxpr, я стараюсь максимально избегать требований к внешним модулям и не планирую использовать многоядерные потоки / потоки.Это «смущающая» распараллеливаемая задача с сотнями тысяч объектов для анализа, я просто распределю список по доступным процессорам во время производства.Я попробую оптимизировать память.
В качестве краткой попытки по Numberxpr с ограничением на 1 поток, выполняя 1 умножение, я получаю время выполнения 40 мс без Numberxpr и 52 мс с Numexpr.
Еще раз спасибо !!

1 Ответ

0 голосов
/ 07 июня 2019

Подход № 1

Мы можем использовать np.einsum для уменьшения суммы за один раз - 

result = np.einsum('ij,ik,il->jkl',a,b,c).ravel()

Кроме того, поэкспериментируйте с флагом optimize в np.einsum, установив его как True для использования BLAS.

Подход № 2

Мы можем использовать broadcasting, чтобы сделать первый шаг, также упомянутый в размещенном коде, а затем использовать тензор-матричное умножение с np.tensordot - 

def broadcast_dot(a,b,c):
    first_multi = a[...,None] * b[:,None]
    return np.tensordot(first_multi,c, axes=(0,0)).ravel()

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

import numexpr as ne

def numexpr_broadcast_dot(a,b,c):
    first_multi = ne.evaluate('A*B',{'A':a[...,None],'B':b[:,None]})
    return np.tensordot(first_multi,c, axes=(0,0)).ravel()

Время для случайных данных с плавающей точкой с заданными размерами набора данных - 

In [36]: %timeit np.einsum('ij,ik,il->jkl',a,b,c).ravel()
4.57 s ± 75.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [3]: %timeit broadcast_dot(a,b,c)
270 ms ± 103 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [4]: %timeit numexpr_broadcast_dot(a,b,c)
172 ms ± 63.8 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

Просто чтобы дать ощущение улучшения с numexpr - 

In [7]: %timeit a[...,None] * b[:,None]
80.4 ms ± 2.64 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [8]: %timeit ne.evaluate('A*B',{'A':a[...,None],'B':b[:,None]})
25.9 ms ± 191 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Это должно быть существенно при расширении этого решения на большее количество входов.

...