将大型稀疏矩阵与其转置相乘的最佳方法是什么？答案

【问题标题】：Which is the best way to multiply a large and sparse matrix with its transpose?将大型稀疏矩阵与其转置相乘的最佳方法是什么？
【发布时间】：2014-08-25 08:10:36
【问题描述】：

我目前想将一个大型稀疏矩阵（~1M x 200k）与其转置相乘。结果矩阵的值将是浮点数。

我尝试将矩阵加载到 scipy 的稀疏矩阵中，并将第一个矩阵的每一行与第二个矩阵相乘。乘法需要大约 2 小时才能完成。

实现这种乘法的有效方法是什么？因为我在计算中看到了一个模式。

矩阵大而稀疏。
矩阵与其转置的乘法。因此，生成的矩阵将是对称的。

我想知道哪些库可以更快地实现计算。它可以是 Python、R、C、C++ 或任何其他语言。

【问题讨论】：

你试过matrix.dot(matrix.transpose())吗？
@user2357112 使用 Scipy？是的.. 由于内存空间，它给出了分段错误。由于生成的矩阵很大，其中包含浮点值...
stackoverflow.com/questions/7477733/… 重复？
你能发布代码吗？你把数据放在什么格式？ scipy.sparse.csr_matrix?
@Jblasco：不是重复的。将矩阵乘以其转置是一种特殊但常见的情况。例如。将向量乘以其转置得到其长度的平方。

标签： python c++ linear-algebra sparse-matrix matrix-multiplication

【解决方案1】：

我想您的主要需求是节省内存。首先，当您将矩阵与其转置相乘时，您不需要任何用于转置的内存：它的所有单元格都可以通过第一个矩阵 (tA[i,j] = A[j,i]) 直接访问。节省了近 1/3 的内存。

我可以看到计算时间也不容忽视。由于生成的矩阵将是对称的，因此您只能计算一半并直接存储另一半。节省了近一半的计算时间。

如果你确定你的初始矩阵真的很稀疏，那么希望得到的结果也是如此，你可以直接将结果存储在 scipy 稀疏矩阵中，COO 格式：只有三个列表来存储非空值。

但是......我不知道有任何库可以做到这一点，你必须自己用你喜欢的语言编写代码（可能是你所说的 scipy 的 python）。

Python 代码示例（矩阵 = A[M][N]）

I = []
J = []
V = []
for i in range(M):
    for j in range(i:M) :
        X = 0.0
        for k in range(N):
            X += A[i ][k] * A[k][j]
        if X != 0.0 # or abs (X) > epsilon if floating point accuracy is a concern ... 
            I.append (i )
            J.append(j)
            V.append(X)
            I.append (j )
            J.append(i)
            V.append(X)

而 I、J、V 是 scipy COO 稀疏矩阵所需要的：

RESULT = sparse.coo_matrix((V,(I,J)),shape=(N, N))

【讨论】：

是的，这确实有效...但是如果生成的矩阵不是稀疏的并且大小为 100000 x 100000 怎么办？
@sravan_kumar 答案就是答案，无论你如何表示

【解决方案2】：

def URMdist(URM):
    NLin=URM.shape[0]
    NCol=URM.shape[1]
    URMt=URM.T
    Result = lil_matrix((NLin,NLin))
    for Lin in range(0,NLin):
        X = 0.0
        for Col in range(Lin,NLin):
            X = URM[Col,:].dot(URMt[:,Lin])
            if X != 0.0: 
                Result[Lin,Col] = Result[Col,Lin] = X
    return Result

【讨论】：

你能仔细检查你的代码吗？我相信你需要在这个函数的某个地方使用 NCol，否则它只有在 NLin == NCol 时才有效，这仅适用于方阵
好的，所以这只有在你这样做时才有效：Result[Lin,Col] = Result[Col,Lin] = X.toarray()[0][0]
我已经有一段时间没有发布这个了，但我很确定我已经测试过了。矩阵是对称的，所以因为我正在计算巨大的矩阵，所以我只计算了一半。是的，这个解决方案只适用于方阵。

【解决方案3】：

您的问题是原始矩阵的剪切尺寸。如果A 具有1e6 行的顺序，则A*A.T 具有1e12 条目的顺序！您绝对不可能在台式计算机中使用这样的矩阵并期望它很快。另外，我认为问题不在于使用 Python，因为 Scipy 中的关键循环是用 C 或 Fortran 实现的（请参阅here）。

话虽如此，我自己也遇到了较小版本的问题（约 200k 行，约 3k cols），我发现如果您小心并且在任何地方都使用稀疏矩阵，Python 可以很好地工作。（我打破了稀疏性，因为我将稀疏矩阵除以向量。不要这样做。）

对于稀疏矩阵，您可以使用正常的重载运算符 * 进行乘法运算。

import scipy.sparse as sp
import numpy as np

# ... in my case my sparse matrix is obtained from real data
# ... but we get something similar below

A = sp.rand(int(2e5), int(3e3), density=1e-3, format='csr')

# If used in the jupyter notebook
%time A * A.T      # Wall time: 2.95 s

# Another option
At = A.T.tocsr()
%time A * At       # Wall time: 2.89 s

我发现很有趣的是，第二个选项并没有明显更快，尽管文档说将两个 csr 矩阵相乘会更快。

请注意，结果还取决于稀疏矩阵的密度。

【讨论】：