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

Question

Я пытаюсь реализовать следующее уравнение, используя редкий пакет scipy:

W = x[:,1] * y[:,1].T + x[:,2] * y[:,2].T + ...

, где x & y - это nxm csc_matrix.В основном я пытаюсь умножить каждый столбец x на каждый столбец y и сложить полученные nxn матрицы вместе.Затем я хочу сделать все ненулевые элементы равными 1.

Это моя текущая реализация:

    c = sparse.csc_matrix((n, n))
    for i in xrange(0,m):
        tmp = bam.id2sym_thal[:,i] * bam.id2sym_cort[:,i].T
        minimum(tmp.data,ones_like(tmp.data),tmp.data)
        maximum(tmp.data,ones_like(tmp.data),tmp.data)

        c = c + tmp

Эта реализация имеет следующие проблемы:

1. Использование памяти, кажется, взрывается.Насколько я понимаю, память должна увеличиваться только по мере того, как c становится менее разреженным, но я вижу, что цикл начинает поглощать> 20 ГБ памяти с = 10 000, m = 100 000 (в каждой строке x & y есть только около 60ноль элементов).

2. Я использую цикл Python, который не очень эффективен.

Мой вопрос: есть ли лучший способсделай это?Моя первая задача - контролировать использование памяти, но было бы здорово сделать это быстрее!

Спасибо!

Answer 1

Обратите внимание, что сумма внешних произведений описанным вами способом просто равна умножению двух матриц вместе.Другими словами,

sum_i X[:,i]*Y[:,i].T == X*Y.T

Так что просто умножьте матрицы вместе.

Z = X*Y.T

Для n = 10000 и m = 100000, где каждый столбец имеет один ненулевой элемент в X и Y, он почти мгновенно вычисляется на моем ноутбуке.

Answer 2

С точки зрения памяти и производительности, это может быть основным кандидатом для использования Cython .

Есть раздел следующей статьи, описывающий его использование с разреженными матрицами:

http://folk.uio.no/dagss/cython_cise.pdf