为什么矩阵指数不能超过一定的大小？

Question

我在使用 scipy.linalg.expm 进行矩阵指数计算时遇到问题。

The code is like this:
import numpy as np    
import scipy.sparse as sp
import scipy.linalg as linA 
from scipy.sparse.linalg import expm 

linA.expm(sp.kron(np.arange(N), np.identity(2)))

其中，N 可以是任意整数，number 是具有对角元素 diag{0,1,2, ... N-1} 的 NxN 对角矩阵，sp.kron 是 scipy.sparse 中的 kronecker 积。这仅适用于 N

linA.expm(sp.kron(np.arange(6), np.identity(2)))

这应该是一个相当简单的代码，但我不知道为什么它会给出以下错误：

NotImplementedError                       Traceback (most recent call last)
<ipython-input-168-b0f10db2dd69> in <module>
----> 1 linA.expm(sp.kron(number(6), identity(2)) )

~\anaconda3\lib\site-packages\scipy\linalg\matfuncs.py in expm(A)
    253     # Input checking and conversion is provided by sparse.linalg.expm().
    254     import scipy.sparse.linalg
--> 255     return scipy.sparse.linalg.expm(A)
    256 
    257 

~\anaconda3\lib\site-packages\scipy\sparse\linalg\matfuncs.py in expm(A)
    590             [  0.        ,   0.        ,  20.08553692]])
    591     """
--> 592     return _expm(A, use_exact_onenorm='auto')
    593 
    594 

~\anaconda3\lib\site-packages\scipy\sparse\linalg\matfuncs.py in _expm(A, use_exact_onenorm)
    675     if structure == UPPER_TRIANGULAR:
    676         # Invoke Code Fragment 2.1.
--> 677         X = _fragment_2_1(X, h.A, s)
    678     else:
    679         # X = r_13(A)^(2^s) by repeated squaring.

~\anaconda3\lib\site-packages\scipy\sparse\linalg\matfuncs.py in _fragment_2_1(X, T, s)
    811             lam_1 = scale * diag_T[k]
    812             lam_2 = scale * diag_T[k+1]
--> 813             t_12 = scale * T[k, k+1]
    814             value = _eq_10_42(lam_1, lam_2, t_12)
    815             X[k, k+1] = value

~\anaconda3\lib\site-packages\scipy\sparse\bsr.py in __getitem__(self, key)
    313 
    314     def __getitem__(self,key):
--> 315         raise NotImplementedError
    316 
    317     def __setitem__(self,key,val):

NotImplementedError: 

Answer 1

根据回溯T[k, k+1] 不起作用，因为T 是bsr 格式的稀疏矩阵，它不实现索引。 （coo 是更常见的格式，也没有这个）。

kron 可以制造 bsr

查看sp.kron代码，https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.kron.html

if (format is None or format == "bsr") and 2*B.nnz >= B.shape[0] * B.shape[1]:
    # B is fairly dense, use BSR
    A = csr_matrix(A,copy=True)
    ...
    return bsr_matrix((data,A.indices,A.indptr), shape=output_shape)

所以在sp.kron(number(N), np.identity(2))

In [251]: B = sparse.coo_matrix(np.identity(2))
In [252]: B
Out[252]: 
<2x2 sparse matrix of type '<class 'numpy.float64'>'
    with 2 stored elements in COOrdinate format>
In [253]: 2*B.nnz >= B.shape[0]*B.shape[1]
Out[253]: True

In [254]: sparse.kron(np.arange(4).reshape(2,2), np.identity(2))
Out[254]: 
<4x4 sparse matrix of type '<class 'numpy.float64'>'
    with 12 stored elements (blocksize = 2x2) in Block Sparse Row format>

测试

In [258]: lg.expm(sparse.kron(np.identity(6), np.identity(2)))
/usr/local/lib/python3.8/dist-packages/scipy/sparse/linalg/dsolve/linsolve.py:144: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format
  warn('spsolve requires A be CSC or CSR matrix format',
/usr/local/lib/python3.8/dist-packages/scipy/sparse/linalg/dsolve/linsolve.py:215: SparseEfficiencyWarning: spsolve is more efficient when sparse b is in the CSC matrix format
  warn('spsolve is more efficient when sparse b '
Out[258]: 
<12x12 sparse matrix of type '<class 'numpy.float64'>'
    with 12 stored elements in Compressed Sparse Column format>

更改为 csc 以避免此警告：

In [265]: lg.expm(sparse.kron(np.identity(6), np.identity(2)).tocsc())
Out[265]: 
<12x12 sparse matrix of type '<class 'numpy.float64'>'
    with 12 stored elements in Compressed Sparse Column format>

因此，简单地向expm 提供bsr 不会导致您的错误。看起来我们必须检查expm 还会发生什么。几年前我看过这个函数（和 MATLAB 的）。它使用包含inv（即spsolve(I,A)）的pade 近似值。这是一个复杂的函数，它尝试不同的事情，包括不同的Pade 订单。

所以你必须告诉我们更多关于这个number 和kron() 结果的性质。我的猜测都没有重现您的错误。

https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.expm.html#scipy.sparse.linalg.expm

上三角形

更正，回溯告诉我们它检测到您的矩阵是upper triangular：

if structure == UPPER_TRIANGULAR:
    # Invoke Code Fragment 2.1.
    X = _fragment_2_1(X, h.A, s)

所以有更多的代码可以追踪。

无论如何，在将矩阵传递给expm 之前进行tocsc 转换可能会解决问题：

lg.expm(sp.kron(...).tocsc())

测试小上三角阵列

In [268]: A = np.array([[1,2,3],[0,4,5],[0,0,6]])
In [269]: M = sparse.bsr_matrix(A)
In [270]: M
Out[270]: 
<3x3 sparse matrix of type '<class 'numpy.int64'>'
    with 6 stored elements (blocksize = 1x1) in Block Sparse Row format>

你的错误：

In [271]: lg.expm(M)
/usr/local/lib/python3.8/dist-packages/scipy/sparse/linalg/dsolve/linsolve.py:144: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format
  warn('spsolve requires A be CSC or CSR matrix format',
/usr/local/lib/python3.8/dist-packages/scipy/sparse/linalg/dsolve/linsolve.py:215: SparseEfficiencyWarning: spsolve is more efficient when sparse b is in the CSC matrix format
  warn('spsolve is more efficient when sparse b '
Traceback (most recent call last):
  File "<ipython-input-271-d1b1437dc466>", line 1, in <module>
    lg.expm(M)
  File "/usr/local/lib/python3.8/dist-packages/scipy/sparse/linalg/matfuncs.py", line 592, in expm
    return _expm(A, use_exact_onenorm='auto')
  File "/usr/local/lib/python3.8/dist-packages/scipy/sparse/linalg/matfuncs.py", line 677, in _expm
    X = _fragment_2_1(X, h.A, s)
  File "/usr/local/lib/python3.8/dist-packages/scipy/sparse/linalg/matfuncs.py", line 813, in _fragment_2_1
    t_12 = scale * T[k, k+1]
  File "/usr/local/lib/python3.8/dist-packages/scipy/sparse/bsr.py", line 315, in __getitem__
    raise NotImplementedError
NotImplementedError

使用 csc 校正：

In [272]: lg.expm(M.tocsc())
Out[272]: 
<3x3 sparse matrix of type '<class 'numpy.float64'>'
    with 6 stored elements in Compressed Sparse Column format>

与np.diag(np.arange(N))

In [303]: sparse.kron(np.diag(np.arange(3)), np.identity(2)).A
Out[303]: 
array([[0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0.],
       [0., 0., 1., 0., 0., 0.],
       [0., 0., 0., 1., 0., 0.],
       [0., 0., 0., 0., 2., 0.],
       [0., 0., 0., 0., 0., 2.]])
In [304]: sparse.kron(np.diag(np.arange(5)), np.identity(2))
Out[304]: 
<10x10 sparse matrix of type '<class 'numpy.float64'>'
    with 16 stored elements (blocksize = 2x2) in Block Sparse Row format>
In [305]: sparse.kron(np.diag(np.arange(6)), np.identity(2))
Out[305]: 
<12x12 sparse matrix of type '<class 'numpy.float64'>'
    with 20 stored elements (blocksize = 2x2) in Block Sparse Row format>

kron 结果没有显着差异，只是随着N 的大小增加。

In [308]: lg.expm(sparse.kron(np.diag(np.arange(6)), np.identity(2)))
 ...
    t_12 = scale * T[k, k+1]
  File "/usr/local/lib/python3.8/dist-packages/scipy/sparse/bsr.py", line 315, in __getitem__
    raise NotImplementedError
NotImplementedError

在kron 中指定csc 格式可以避免该错误（我们可以忽略此效率警告）：

In [309]: lg.expm(sparse.kron(np.diag(np.arange(6)), np.identity(2),'csc'))
/usr/local/lib/python3.8/dist-packages/scipy/sparse/_index.py:82: SparseEfficiencyWarning: Changing the sparsity structure of a csc_matrix is expensive. lil_matrix is more efficient.
  self._set_intXint(row, col, x.flat[0])
Out[309]: 
<12x12 sparse matrix of type '<class 'numpy.float64'>'
    with 23 stored elements in Compressed Sparse Column format>

为什么N=6 给出这个警告，而不是小的N 可能与它必须尝试的Pade 命令有关。请记住，expm 是一个复杂的计算，它所能做的（数字上）最好的就是近似它。对于小矩阵，这种近似更容易。这段代码背后有很多数学理论。