稀疏矩阵乘以对角矩阵

Question

我想实现以下操作： B = diag(u) @ K @ diag(v) - 一些矩阵 K 从左到右乘以对角矩阵。显式创建密集矩阵和 2 个矩阵矩阵乘积非常昂贵。如果K 是np.ndarray，我想出了以下解决方案：

u[:, None] * K * v[None, :]

效果很好，因为广播发生了，但是如果 K 是 scipy.sparse.csr_matrix 我会遇到尺寸不匹配的问题：

ValueError: dimension mismatch

我试过了

u[:, None] * K.multiply(v[None, :])

计算了左操作数，但乘法也失败了。应该如何对稀疏矩阵实现相同的操作？

编辑

我对较小的数据进行了测试，在这种情况下没有这样的错误，但np.ndarray 和csr_matrix 的结果不同：

import numpy as np
from scipy.sparse import csr_matrix

A = np.random.randint(low = 0, high = 10, size = (5, 5))
A_sp = csr_matrix(A)
x = np.random.randint(low = 1, high = 6, size = 5)
y = np.random.randint(low = 1, high = 6, size = 5)

x[None, :] * A * y[:, None]
# output
array([[ 72,  48, 128,  60,  48],
       [ 45,  20, 140, 175,  60],
       [ 18,  24,  32,  15,  28],
       [ 42,  24,   0,  70,  16],
       [ 30,  40,  32,  60,   0]])

x[None, :] * A_sp * y[:, None]
# output
array([[436, 288, 400, 432, 236],
       [545, 360, 500, 540, 295],
       [109,  72, 100, 108,  59],
       [218, 144, 200, 216, 118],
       [218, 144, 200, 216, 118]])

Answer 1

像你这样的 3 个数组：

In [33]: A,x,y
Out[33]: 
(array([[7, 4, 2, 9, 0],
        [6, 2, 5, 7, 4],
        [6, 6, 4, 3, 1],
        [5, 6, 5, 1, 2],
        [4, 8, 6, 5, 6]]),
 array([5, 3, 3, 3, 4]),
 array([5, 4, 4, 2, 2]))

并制作对角二维数组：

In [34]: X = np.diag(x); Y = np.diag(y)

稀疏等价物：

In [35]: A_sp
Out[35]: 
<5x5 sparse matrix of type '<class 'numpy.int64'>'
    with 24 stored elements in Compressed Sparse Row format>

从密集的对角线生成稀疏的对角线：

In [36]: x_sp = sparse.csr_matrix(X); y_sp = sparse.csr_matrix(Y)
In [37]: x_sp
Out[37]: 
<5x5 sparse matrix of type '<class 'numpy.int64'>'
    with 5 stored elements in Compressed Sparse Row format>

或直接：

In [38]: x_sp = sparse.diags(x); y_sp = sparse.diags(y)
In [39]: x_sp
Out[39]: 
<5x5 sparse matrix of type '<class 'numpy.float64'>'
    with 5 stored elements (1 diagonals) in DIAgonal format>

您的原始矩阵乘积：

In [40]: X@A@Y
Out[40]: 
array([[175,  80,  40,  90,   0],
       [ 90,  24,  60,  42,  24],
       [ 90,  72,  48,  18,   6],
       [ 75,  72,  60,   6,  12],
       [ 80, 128,  96,  40,  48]])

这也适用于稀疏的A_sp：

In [41]: X@A_sp@Y
Out[41]: 
array([[175,  80,  40,  90,   0],
       [ 90,  24,  60,  42,  24],
       [ 90,  72,  48,  18,   6],
       [ 75,  72,  60,   6,  12],
       [ 80, 128,  96,  40,  48]])

请注意，第一个 @ 会生成一个二维密集数组：

In [42]: X@A_sp
Out[42]: 
array([[35, 20, 10, 45,  0],
       [18,  6, 15, 21, 12],
       [18, 18, 12,  9,  3],
       [15, 18, 15,  3,  6],
       [16, 32, 24, 20, 24]])

稀疏对角线也一样：

In [43]: x_sp@A_sp
Out[43]: 
<5x5 sparse matrix of type '<class 'numpy.float64'>'
    with 24 stored elements in Compressed Sparse Row format>
In [44]: x_sp@A_sp@y_sp
Out[44]: 
<5x5 sparse matrix of type '<class 'numpy.float64'>'
    with 24 stored elements in Compressed Sparse Row format>
In [45]: _.A
Out[45]: 
array([[175.,  80.,  40.,  90.,   0.],
       [ 90.,  24.,  60.,  42.,  24.],
       [ 90.,  72.,  48.,  18.,   6.],
       [ 75.,  72.,  60.,   6.,  12.],
       [ 80., 128.,  96.,  40.,  48.]])

广播的等价物：

In [47]: x[:,None]*A*y[None,:]
Out[47]: 
array([[175,  80,  40,  90,   0],
       [ 90,  24,  60,  42,  24],
       [ 90,  72,  48,  18,   6],
       [ 75,  72,  60,   6,  12],
       [ 80, 128,  96,  40,  48]])

但广播不适用于稀疏乘法

In [49]: y[None,:]*A_sp         # same as @
Out[49]: array([[101,  80,  68,  97,  36]])

我们可以使用带有稀疏元素乘法的向量。广播在这里确实有效，但一定要使用sparse 方法（无运算符）：

In [50]: A_sp.multiply(y[None,:])
Out[50]: 
<5x5 sparse matrix of type '<class 'numpy.int64'>'
    with 24 stored elements in COOrdinate format>
In [51]: _.A
Out[51]: 
array([[35, 16,  8, 18,  0],
       [30,  8, 20, 14,  8],
       [30, 24, 16,  6,  2],
       [25, 24, 20,  2,  4],
       [20, 32, 24, 10, 12]])
In [52]: A_sp.multiply(y[None,:]).multiply(x[:,None])
Out[52]: 
<5x5 sparse matrix of type '<class 'numpy.int64'>'
    with 24 stored elements in COOrdinate format>
In [53]: _.A
Out[53]: 
array([[175,  80,  40,  90,   0],
       [ 90,  24,  60,  42,  24],
       [ 90,  72,  48,  18,   6],
       [ 75,  72,  60,   6,  12],
       [ 80, 128,  96,  40,  48]])

我会让你看看相对速度。这里的数组很小，A 并不稀疏。