为什么矩阵指数不能超过一定大小？
Why does matrix exponential not work beyond certain size?
我在使用 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是一个NxN的对角矩阵，对角元素为diag{0,1,2, ... N-1}，sp.kron是[中的克罗内克积=21=]。这仅适用于 N<6。当我尝试 运行 N=6
的代码时
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:
根据回溯，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() 结果的性质。 None 我的猜测重现了您的错误。
上三角
更正，回溯告诉我们它检测到您的矩阵是 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>
除了大小随着 N 的增长外，kron 结果没有显着差异。
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 是一个复杂的计算，它可以做的最好的（数字）是近似值。对于小矩阵，这种近似更容易。这段代码背后有很多数学理论。
我在使用 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是一个NxN的对角矩阵，对角元素为diag{0,1,2, ... N-1}，sp.kron是[中的克罗内克积=21=]。这仅适用于 N<6。当我尝试 运行 N=6
的代码时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:
根据回溯，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() 结果的性质。 None 我的猜测重现了您的错误。
上三角
更正，回溯告诉我们它检测到您的矩阵是 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>
除了大小随着 N 的增长外，kron 结果没有显着差异。
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 是一个复杂的计算，它可以做的最好的（数字）是近似值。对于小矩阵，这种近似更容易。这段代码背后有很多数学理论。