4

当我有两个非稀疏矩阵A和 我有存储在此处指定的 CSC 格式的所需索引。BC=A.T.dot(B)CC

4

3 回答 3

2

如果您事先知道您想要 C 的哪些部分并且其中一些部分是连续的矩形区域*，那么您可以使用与分区矩阵乘法 ( 1 ) 或块矩阵乘法 ( 2 ) 相关的矩阵代数规则来加速这些计算中的一些。因此，例如，您可以使用与@GaryBishop 相同的基本逻辑，但是您可以使用由 i_start、i_end 和 j_start、j_end 组成的四元组的列表（或数组）来代替“i”和“j”元素的列表定义 C 的子矩阵，那么您可以使用这些索引（尽管在这些链接中建立规则）来找出您需要求解所需的 C 块的 A 和 B 的子矩阵。

举个简单的例子，假设你只想要 C 的中间块，所以我们将 C 按行划分为 C1、C2 和 C3，我们只关心 C2。如果 A^{T} 同样被划分为三组行 A1、A2、A3 则 C2 = A2 * B。这个想法推广到任何形状的矩形，它只需要 A 和 B 的不同划分来计算。这个想法是一样的。

  • - 这是微不足道的，但如果区域大于单个元素，您只会节省时间。
于 2012-12-21T12:54:36.037 回答
2

您可以让 numpy 进行循环，而不是使用 Python（GaryBishop 的回答）对坐标进行迭代，这构成了显着的加速（时间如下）：

def sparse_mult(a, b, coords) :
    rows, cols = zip(*coords)
    rows, r_idx = np.unique(rows, return_inverse=True)
    cols, c_idx = np.unique(cols, return_inverse=True)
    C = np.dot(a[rows, :], b[:, cols])
    return C[r_idx, c_idx]

>>> A = np.arange(12).reshape(3, 4)
>>> B = np.arange(15).reshape(3, 5)
>>> np.dot(A.T, B)
array([[100, 112, 124, 136, 148],
       [115, 130, 145, 160, 175],
       [130, 148, 166, 184, 202],
       [145, 166, 187, 208, 229]])
>>> sparse_mult(A.T, B, [(0, 0), (1, 2), (2, 4), (3, 3)])
array([100, 145, 202, 208])

sparse_mult返回您作为元组列表提供的坐标处的值的扁平数组。我对稀疏矩阵格式不是很熟悉，所以我不知道如何从上述数据中定义 CSC，但以下工作：

>>> coords = [(0, 0), (1, 2), (2, 4), (3, 3)]
>>> sparse.coo_matrix((sparse_mult(A.T, B, coords), zip(*coords))).tocsc()
<4x5 sparse matrix of type '<type 'numpy.int32'>'
    with 4 stored elements in Compressed Sparse Column format>

这是各种替代方案的时间安排：

>>> import timeit
>>> a = np.random.rand(2000, 3000)
>>> b = np.random.rand(3000, 5000)
>>> timeit.timeit('np.dot(a,b)[[0, 0, 1999, 1999], [0, 4999, 0, 4999]]', 'from __main__ import np, a, b', number=1)
5.848562187263569
>>> timeit.timeit('sparse_mult(a, b, [(0, 0), (0, 4999), (1999, 0), (1999, 4999)])', 'from __main__ import np, a, b, sparse_mult', number=1)
0.0018596387374678613
>>> np.dot(a,b)[[0, 0, 1999, 1999], [0, 4999, 0, 4999]]
array([ 758.76351111,  750.32613815,  751.4614542 ,  758.8989648 ])
>>> sparse_mult(a, b, [(0, 0), (0, 4999), (1999, 0), (1999, 4999)])
array([ 758.76351111,  750.32613815,  751.4614542 ,  758.8989648 ])
于 2013-01-07T23:23:15.600 回答
1

忽略 CSC 业务，也许回答一个比你问的更简单的问题。在给定 C 索引值的元组列表的情况下，我将如何计算 C 元素的子集。

由于您正在评估 C=ATdot(B) ，因此您将 A 的列乘以 B 的列。所以，

for i, j in indexList:
    C[i, j] = np.dot(A[:,i], B[:,j])

我猜这不是你要找的，但我发现简单的答案有时有助于澄清问题。

于 2012-12-07T16:50:32.437 回答