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

Question

我现在想把一个很大的稀疏矩阵（大约100万行和20万列）和它的转置矩阵相乘。结果矩阵的值会是浮点数。

• 我尝试用scipy的稀疏矩阵来加载这个矩阵，然后把第一个矩阵的每一行和第二个矩阵相乘。这个乘法花了大约2个小时才完成。

有没有更高效的方法来完成这个乘法呢？因为我发现了计算中的一些规律。

• 这个矩阵很大，而且是稀疏的。
• 矩阵和它的转置相乘，所以结果矩阵会是对称的。

我想知道有哪些库可以更快地完成这个计算。可以是Python、R、C、C++或者其他语言的库。

Answer 1

你遇到的问题主要是因为你原始矩阵的大小。如果 A 有大约 1e6 行，那么 A*A.T 的结果会有大约 1e12 个元素！在你的桌面电脑上处理这么大的矩阵是完全不现实的，也不可能快。此外，我认为问题不在于使用 Python，因为 Scipy 中的关键循环是用 C 或 Fortran 实现的（你可以在 这里查看）。

话虽如此，我自己也遇到过一个较小的类似问题（大约 20 万行，3000 列），我发现如果你小心使用稀疏矩阵，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 矩阵相乘会更快。

请注意，结果还会依赖于你稀疏矩阵的密度。

Answer 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

当然可以！请把你想要翻译的内容发给我，我会帮你把它变得更简单易懂。

Answer 3

我想你主要是想节省内存。首先，当你把一个矩阵和它的转置相乘时，其实不需要为转置分配内存：转置矩阵的所有元素都可以直接通过原矩阵访问（tA[i,j] = A[j,i]）。这样可以节省将近三分之一的内存。

我也明白计算时间也是很重要的。因为结果矩阵是对称的，你只需要计算一半的内容，然后直接存储另一半。这样可以节省将近一半的计算时间。

如果你确定你的初始矩阵是稀疏的，也就是说大部分元素都是零，那么你可以直接把结果存储在一个scipy的稀疏矩阵中，使用COO格式：只需要三个列表来存储非零值。

不过……我不知道有什么库可以做到这一点，你可能需要用你喜欢的编程语言自己编写代码（可能是python，因为你提到了scipy）。

下面是一个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))