如何在 Python 中处理大矩阵和矩阵乘法

Question

我正在尝试执行具有以下方案的矩阵乘法：

C = np.dot(np.dot(sparse.csr_matrix(np.double(A).transpose()),sparse.spdiags(B,0,Ngrid,Ngrid)), sparse.csr_matrix(np.double(A)))

因此，我想对矩阵 A 进行转置，得到 M>>N 的 N x M 矩阵，然后乘以对角矩阵，即 M x M 矩阵。 B是“主对角线”。生成的矩阵 (N x M) 应与矩阵 A (M x N) 相乘并得出 N x N 矩阵 C.

出现的错误如下：

<2000x921600 sparse matrix of type '<class 'numpy.float64'>'
    with 1843066024 stored elements in Compressed Sparse Row format>

由于最终的矩阵是N x N，我想把这个矩阵作为一个numpy数组。你怎么看，我试着把中间的矩阵做成一个效果很好的稀疏对角矩阵。但是，我不明白为什么 Python 需要这个具有 1843066024 个元素的疯狂大矩阵来进行乘法运算。

你有什么想法吗and/or解释一下为什么会出现这个问题？

Answer 1

你这样做...太复杂了。这是 M >> N 的直接路径（您对此不一致）。

import numpy as np

B = sparse.spdiags(B,0,Ngrid,Ngrid)) # [M x M] sparse
A = np.ndarray(..., dtype=float)     # [M x N] dense

C = A.T @ B   # [N x M] dense
C = C @ A     # [N x N] dense

C 就是你想要的数组。 None 的中间产品是 M x M。如果您仍然有内存问题，您要么需要获得更多内存，要么需要将问题分成 m 轴上的更小部分并分段计算。

Answer 2

如果B是对角线，就不需要用稀疏来节省内存。您可以使用一维数组进行计算，只是对角线值。

我将用小尺寸进行演示，制作完整的 B 不会造成问题。其他人可以用大尺寸测试这个。

In [5]: A = np.arange(24).reshape(3,8)
In [7]: B = np.diag(np.arange(8))
In [8]: B.shape
Out[8]: (8, 8)

双矩阵相乘：

In [10]: A@B@A.T
Out[10]: 
array([[  784,  1904,  3024],
       [ 1904,  4816,  7728],
       [ 3024,  7728, 12432]])

相当于einsum：

In [12]: np.einsum('ij,jk,lk',A,B,A)
Out[12]: 
array([[  784,  1904,  3024],
       [ 1904,  4816,  7728],
       [ 3024,  7728, 12432]])

B的一维对角线：

In [15]: b = np.diag(B)

广播乘法与 matmul 做同样的事情：

In [17]: np.allclose(A*b,A@B)
Out[17]: True
In [18]: np.allclose(A*b@A.T,A@B@A.T)
Out[18]: True

用einsum表示：

In [19]: np.einsum('ij,j,lj',A,b,A)
Out[19]: 
array([[  784,  1904,  3024],
       [ 1904,  4816,  7728],
       [ 3024,  7728, 12432]])

一些对比时间：

In [20]: timeit np.einsum('ij,j,lj',A,b,A)
10.6 µs ± 22.1 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
In [21]: timeit np.einsum('ij,jk,lk',A,B,A)
15.2 µs ± 13.3 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
In [22]: timeit A@B@A.T
7.26 µs ± 20.4 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
In [23]: timeit A*b@A.T
8.17 µs ± 12.6 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

对于 einsum，'condensed' 版本更快。使用 matmul 时，整个对角线稍微快一些。

但是对于创建完整 B 的大型数组可能会出现问题，使用 b 可能会更快。在其他 SO 中也观察到，由于更好的内存处理，较小数组上的迭代可以更快。

np.array([A[i,:]*b@A.T for i in range(3)])

如何在 Python 中处理大矩阵和矩阵乘法

How to handle large matrices and matrix multiplication in Python

python

numpy

matrix

sparse-matrix

matrix-multiplication