Эффективное распределение многопоточных и плотных матриц Numpy/Scipy

Я работаю над реализацией следующего уравнения:

X =(Y.T * Y + Y.T * C * Y) ^ -1

Y - матрица (n x f), а C - (n x n) - диагональная; n составляет около 300 тыс., а f будет варьироваться от 100 до 200. В рамках процесса оптимизации это уравнение будет использоваться почти 100 миллионов раз, поэтому его нужно обрабатывать очень быстро.

Y инициализируется случайным образом, а C - очень разреженная матрица, и только несколько чисел из 300k по диагонали будут отличаться от 0. Поскольку диагональные функции Numpy создают плотные матрицы, я создал C как разреженную матрицу csr. Но при попытке решить первую часть уравнения:

r = dot(C, Y)

Компьютер выходит из строя из-за ограничений памяти. Затем я решил преобразовать Y в csr_matrix и выполнить ту же операцию:

r = dot(C, Ysparse)

и этот подход занял 1,38 мс. Но это решение несколько "сложно", так как я использую разреженную матрицу для хранения плотной, интересно, насколько это эффективно.

Итак, мой вопрос в том, есть ли способ умножения разреженного C и плотного Y без необходимости превращать Y в разреженный и улучшать производительность? Если бы какой-то С мог быть представлен как диагональный плотный, не потребляя тонны памяти, возможно, это привело бы к очень эффективной работе, но я не знаю, возможно ли это.

Я ценю вашу помощь!

Ответ 1

Причина, по которой dot-продукт сталкивается с проблемами памяти при вычислении r = dot (C, Y), заключается в том, что функция numpy dot не имеет встроенной поддержки обработки разреженных матриц. То, что происходит, это numpy, думает о разреженной матрице C как о объекте python, а не о массиве numpy. Если вы осмотрите мелкий масштаб, вы можете увидеть проблему из первых рук:

>>> from numpy import dot, array
>>> from scipy import sparse
>>> Y = array([[1,2],[3,4]])
>>> C = sparse.csr_matrix(array([[1,0], [0,2]]))
>>> dot(C,Y)
array([[  (0, 0)    1
  (1, 1)    2,   (0, 0) 2
  (1, 1)    4],
  [  (0, 0) 3
  (1, 1)    6,   (0, 0) 4
  (1, 1)    8]], dtype=object)

Ясно, что приведенное выше не является тем результатом, который вас интересует. Вместо этого вы должны вычислить с помощью функции scipy sparse.csr_matrix.dot:

r = sparse.csr_matrix.dot(C, Y)

или более компактно

r = C.dot(Y)

Ответ 2

Try:

import numpy as np
from scipy import sparse

f = 100
n = 300000

Y = np.random.rand(n, f)
Cdiag = np.random.rand(n) # diagonal of C
Cdiag[np.random.rand(n) < 0.99] = 0

# Compute Y.T * C * Y, skipping zero elements
mask = np.flatnonzero(Cdiag)
Cskip = Cdiag[mask]

def ytcy_fast(Y):
    Yskip = Y[mask,:]
    CY = Cskip[:,None] * Yskip  # broadcasting
    return Yskip.T.dot(CY)

%timeit ytcy_fast(Y)

# For comparison: all-sparse matrices
C_sparse = sparse.spdiags([Cdiag], [0], n, n)
Y_sparse = sparse.csr_matrix(Y)
%timeit Y_sparse.T.dot(C_sparse * Y_sparse)

Мои тайминги:

In [59]: %timeit ytcy_fast(Y)
100 loops, best of 3: 16.1 ms per loop

In [18]: %timeit Y_sparse.T.dot(C_sparse * Y_sparse)
1 loops, best of 3: 282 ms per loop

Ответ 3

Во-первых, вы действительно уверены, что вам нужно выполнить полную матричную инверсию в своей проблеме? Большую часть времени действительно нужно вычислить x = A ^ -1 y, что намного проще решить.

Если это действительно так, я бы рассмотрел вычисление аппроксимации обратной матрицы вместо полной инверсии матрицы. Поскольку инверсия матрицы действительно дорогостоящая. См., Например, алгоритм Lanczos для эффективного приближения обратной матрицы. Аппроксимация может храниться редко в качестве бонуса. Кроме того, для этого требуются только операции с матричным вектором, поэтому вам даже не нужно хранить полную матрицу для инверсии.

В качестве альтернативы, используя pyoperators, вы также можете использовать метод .todense для вычисления матрицы для инверсии с использованием эффективных векторных векторных операций. Существует специальный разреженный контейнер для диагональных матриц.

Для реализации алгоритма Ланцоша вы можете взглянуть на pyoperators (отказ от ответственности: я являюсь одним из соавторов этого часть программного обеспечения).