Сумма строк точечного произведения для огромной матрицы в python

У меня есть 2 матрицы 100кх200 и 200кх100к если бы они были маленькой матрицей, я бы просто использовал numpy dot product 

sum(a.dot(b), axis = 0)

Однако матрица слишком велика, а также я не могу использовать циклы есть ли разумный способ сделать это?

3 14
python numpy matrix-multiplication

3 ответа:

Возможной оптимизацией является 

>>> numpy.sum(a @ b, axis=0)
array([  1.83633615,  18.71643672,  15.26981078, -46.33670382,  13.30276476])
>>> numpy.sum(a, axis=0) @ b
array([  1.83633615,  18.71643672,  15.26981078, -46.33670382,  13.30276476])

Вычисление a @ b требует 10k×200×10k операций, в то время как суммирование строк сначала уменьшит умножение до 1×200×10k операций, давая 10k× улучшение.

Это в основном связано с распознаванием 

   numpy.sum(x, axis=0) == [1, 1, ..., 1] @ x
=> numpy.sum(a @ b, axis=0) == [1, 1, ..., 1] @ (a @ b)
                            == ([1, 1, ..., 1] @ a) @ b
                            == numpy.sum(a, axis=0) @ b

Аналогично для другой оси.

>>> numpy.sum(a @ b, axis=1)
array([  2.8794171 ,   9.12128399,  14.52009991,  -8.70177811, -15.0303783 ])
>>> a @ numpy.sum(b, axis=1)
array([  2.8794171 ,   9.12128399,  14.52009991,  -8.70177811, -15.0303783 ])

(Примечание: x @ y эквивалентно x.dot(y) для 2D матриц и 1D векторов на Python 3.5+ с numpy 1.10.0+)

$ INITIALIZATION='import numpy;numpy.random.seed(0);a=numpy.random.randn(1000,200);b=numpy.random.rand(200,1000)'

$ python3 -m timeit -s "$INITIALIZATION" 'numpy.einsum("ij,jk->k", a, b)'
10 loops, best of 3: 87.2 msec per loop

$ python3 -m timeit -s "$INITIALIZATION" 'numpy.sum(a@b, axis=0)'
100 loops, best of 3: 12.8 msec per loop

$ python3 -m timeit -s "$INITIALIZATION" 'numpy.sum(a, axis=0)@b'
1000 loops, best of 3: 300 usec per loop

Иллюстрация: 

In [235]: a = np.random.rand(3,3)
array([[ 0.465,  0.758,  0.641],
       [ 0.897,  0.673,  0.742],
       [ 0.763,  0.274,  0.485]])

In [237]: b = np.random.rand(3,2)
array([[ 0.303,  0.378],
       [ 0.039,  0.095],
       [ 0.192,  0.668]])
Теперь, если мы просто сделаем a @ b, нам понадобится 18 операций умножения и 6 операций сложения. С другой стороны, если мы сделаем np.sum(a, axis=0) @ b, нам понадобится только 6 умножений и 2 операции сложения. Улучшение в 3 раза, потому что у нас было 3 строки в a. Что касается случая опа, это должно дать 10k-кратное улучшение по сравнению с простым вычислением a @ b, так как у него есть 10k строк в a.
18

Происходит два sum-reductions - Один из марикса-мультиликации с np.dot, а затем с явным sum.

Мы могли бы использовать np.einsum чтобы сделать и то и другое за один раз, как так - 

np.einsum('ij,jk->k',a,b)

Пробный прогон - 

In [27]: a = np.random.rand(3,4)

In [28]: b = np.random.rand(4,3)

In [29]: np.sum(a.dot(b), axis = 0)
Out[29]: array([ 2.70084316,  3.07448582,  3.28690401])

In [30]: np.einsum('ij,jk->k',a,b)
Out[30]: array([ 2.70084316,  3.07448582,  3.28690401])

Тест времени выполнения - 

In [45]: a = np.random.rand(1000,200)

In [46]: b = np.random.rand(200,1000)

In [47]: %timeit np.sum(a.dot(b), axis = 0)
100 loops, best of 3: 5.5 ms per loop

In [48]: %timeit np.einsum('ij,jk->k',a,b)
10 loops, best of 3: 71.8 ms per loop

К сожалению, не похоже, что мы делаем что-то лучше с np.einsum.

Для изменения на np.sum(a.dot(b), axis = 1), просто замените нотацию выходной строки там - np.einsum('ij,jk->i',a,b), например так - 

In [42]: np.sum(a.dot(b), axis = 1)
Out[42]: array([ 3.97805141,  3.2249661 ,  1.85921549])

In [43]: np.einsum('ij,jk->i',a,b)
Out[43]: array([ 3.97805141,  3.2249661 ,  1.85921549])
5

Несколько быстрых тестов времени, используя идею, которую я добавил к ответу Дивакара: 

In [162]: a = np.random.rand(1000,200)
In [163]: b = np.random.rand(200,1000)

In [174]: timeit c1=np.sum(a.dot(b), axis=0)
10 loops, best of 3: 27.7 ms per loop

In [175]: timeit c2=np.sum(a,axis=0).dot(b)
1000 loops, best of 3: 432 µs per loop

In [176]: timeit c3=np.einsum('ij,jk->k',a,b)
10 loops, best of 3: 170 ms per loop

In [177]: timeit c4=np.einsum('j,jk->k', np.einsum('ij->j', a), b)
1000 loops, best of 3: 353 µs per loop

In [178]: timeit np.einsum('ij->j', a) @b
1000 loops, best of 3: 304 µs per loop

einsum На самом деле быстрее, чем np.sum!

In [180]: timeit np.einsum('ij->j', a)
1000 loops, best of 3: 173 µs per loop
In [181]: timeit np.sum(a,0)
1000 loops, best of 3: 312 µs per loop

Для больших массивов преимущество einsum уменьшается 

In [183]: a = np.random.rand(100000,200)
In [184]: b = np.random.rand(200,100000)
In [185]: timeit np.einsum('ij->j', a) @b
10 loops, best of 3: 51.5 ms per loop
In [186]: timeit c2=np.sum(a,axis=0).dot(b)
10 loops, best of 3: 59.5 ms per loop
2