我正在嘗試執行具有以下方案的矩陣乘法：
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，這會導致一個 N x M 矩陣，其中 M>>N 并與對角矩陣相乘，它是一個 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 個元素的瘋狂大矩陣來進行乘法運算。
你有什么想法和/或解釋為什么會出現這個問題嗎？
uj5u.com熱心網友回復：
如果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)
雙matmul：
In [10]: A@[email protected]
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*[email protected],A@[email protected])
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@[email protected]
7.26 μs ± 20.4 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
In [23]: timeit A*[email protected]
8.17 μs ± 12.6 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
因為einsum“濃縮”版本更快。使用matmul完整的對角線稍微快一些。
但是對于創建完整的大型陣列B可能會出現問題，使用b可能會更快。在其他 SO 中也觀察到，由于更好的記憶體處理，較小陣列上的迭代可以更快。
np.array([A[i,:]*[email protected] for i in range(3)])
uj5u.com熱心網友回復：
你這樣做...過于復雜。這是一個簡單的路徑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那么就是你想要的陣列。你們的中間產品都不是M x M。如果您仍然有記憶體問題，您需要獲得更多記憶體，或者您需要將您的問題分成m軸上的較??小部分并分段計算。
轉載請註明出處，本文鏈接：https://www.uj5u.com/net/454652.html