两个巨大的密集矩阵的乘法 Hadamard-乘以稀疏矩阵答案

【问题标题】：Multiplication of two huge dense matrices Hadamard-multiplied by a sparse matrix两个巨大的密集矩阵的乘法 Hadamard-乘以稀疏矩阵
【发布时间】：2021-04-24 09:58:00
【问题描述】：

我有两个稠密矩阵A 和B，每个矩阵都有3e5x100 的大小。另一个稀疏二进制矩阵C，大小为3e5x3e5。我想找到以下数量：C ∘ (AB')，其中∘ 是 Hadamard 产品（即元素方面），B' 是 B 的转置。显式计算AB' 将要求大量内存（~500GB）。由于最终结果不需要整个AB'，因此只计算乘法A_iB_j' 就足够了，其中C_ij != 0，其中A_i 是矩阵i 的列A 和C_ij 是矩阵C 的位置(i,j) 处的元素。建议的方法类似于以下算法：

result = numpy.initalize_sparse_matrix(shape = C.shape)
while True:
 (i,j) = C_ij.pop_nonzero_index() #prototype function returns the nonzero index and then points to the next nonzero index
 if (i,j) is empty:
   break
 result(i,j) = A_iB_j'

然而，这个算法需要太多时间。无论如何使用LAPACK/BLAS 算法来改进它？我正在用 Python 编码，所以我认为 numpy 可以成为 LAPACK/BLAS 的更人性化的包装器。

【问题讨论】：

你的 C 有多稀疏？获取 C 中 1 的所有坐标需要多长时间？我不记得有任何算法可以减少这种计算

标签： python numpy matrix sparse-matrix lapack

【解决方案1】：

假设C 存储为scipy.sparse 矩阵，您可以使用以下方法进行此计算：

C = C.tocoo()
result_data = C.data * (A[C.row] * B[C.col]).sum(1)
result = sparse.coo_matrix((result_data, (row, col)), shape=C.shape)

这里我们展示了结果与一些较小输入的朴素算法相匹配：

import numpy as np
from scipy import sparse

N = 300
M = 10

def make_C(N, nnz=1000):
  data = np.random.rand(nnz)
  row = np.random.randint(0, N, nnz)
  col = np.random.randint(0, N, nnz)
  return sparse.coo_matrix((data, (row, col)), shape=(N, N))


A = np.random.rand(N, M)
B = np.random.rand(N, M)
C = make_C(N)

def f_naive(C, A, B):
  return C.multiply(np.dot(A, B.T))

def f_efficient(C, A, B):
  C = C.tocoo()
  result_data = C.data * (A[C.row] * B[C.col]).sum(1)
  return sparse.coo_matrix((result_data, (C.row, C.col)), shape=C.shape)

np.allclose(
    f_naive(C, A, B).toarray(),
    f_efficient(C, A, B).toarray()
)
# True

在这里我们看到它适用于完整的输入大小：

N = 300000
M = 100

A = np.random.rand(N, M)
B = np.random.rand(N, M)
C = make_C(N)

out = f_efficient(C, A, B)

print(out.shape)
# (300000, 300000)

print(out.nnz)
# 1000

【讨论】：