提高稀疏矩阵的幂

Question

我有一个 10001 行 + 10001 列（有很多 0）的稀疏矩阵，

我正在尝试提高这个稀疏矩阵的能力

即

A = [[1,1],[1,0]]
AS = sparse.csr_matrix(A)
AS

def matrixMul(AS, n):
    if(n <= 1):
        return AS
    else:
        return np.matmul(matrixMul(AS, n-1), AS)

matrixMul(AS, 10)

如果我将 AS 提高到 2 的幂，预期的结果应该是 [[2, 1] [1, 1]]

我想找 AS ^ 10

我应该调用什么函数？我已经尝试了上面的代码，但收到了这个错误。

谢谢。

ValueError Traceback（最近调用 最后）在（） 9 返回 np.matmul(matrixMul(AS, n-1), AS) 10 ---> 11 矩阵Mul(AS, 10)

matrixMul(AS, n) 中的 8 帧 7 返回 AS 8 其他： ----> 9 返回 np.matmul(matrixMul(AS, n-1), AS) 10 11 矩阵Mul(AS, 10)

ValueError: matmul: 输入操作数 0 没有足够的维度 （有 0，带有签名 (n?,k),(k,m?)->(n?,m?) 的 gufunc 核心需要 1)

Answer 1

In [3]: from scipy import sparse
In [4]: A = [[1,1],[1,0]]
   ...: AS = sparse.csr_matrix(A)
In [5]: AS
Out[5]: 
<2x2 sparse matrix of type '<class 'numpy.int64'>'
    with 3 stored elements in Compressed Sparse Row format>
In [6]: AS.A
Out[6]: 
array([[1, 1],
       [1, 0]])
In [7]: (AS@AS@AS).A
Out[7]: 
array([[3, 2],
       [2, 1]])

@ 映射到 AS.__matmul__，但 np.matmul 不以这种方式委派它。

In [8]: np.matmul(AS,AS)
Traceback (most recent call last):
  File "<ipython-input-8-37972b025121>", line 1, in <module>
    np.matmul(AS,AS)
ValueError: matmul: Input operand 0 does not have enough dimensions (has 0, gufunc core with signature (n?,k),(k,m?)->(n?,m?) requires 1)

更正：

In [9]: def matrixMul(AS, n):
   ...:     if(n <= 1):
   ...:         return AS
   ...:     else:
   ...:         return matrixMul(AS, n-1)@ AS
   ...: 
In [10]: matrixMul(AS,3)
Out[10]: 
<2x2 sparse matrix of type '<class 'numpy.int64'>'
    with 4 stored elements in Compressed Sparse Row format>
In [11]: _.A
Out[11]: 
array([[3, 2],
       [2, 1]])
In [12]: matrixMul(AS,10).A
Out[12]: 
array([[89, 55],
       [55, 34]])

如 cmets 所示，** 电源工作：

In [15]: (AS**10).A
Out[15]: 
array([[89, 55],
       [55, 34]])

稀疏矩阵以np.matrix 为模型。 * 是矩阵乘法，** 是矩阵幂。

此幂次与此半智能显式乘法大致相同：

def foo(AS):
    AS2=AS@AS
    return AS2@AS2@AS2@AS2@AS2
In [33]: timeit foo(AS).A
865 µs ± 378 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [34]: timeit AS**10
767 µs ± 189 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Answer 2

您可以使用幂运算符：AS ** 10 与普通数学符号中的AS ^ 10 相同。

def matrixMul(AS, n):
    if(n <= 1):
        return AS
    else:
        return AS ** n