带
我在使用scipy.linal.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乘积；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)) )
~anaconda3libsite-packagesscipylinalgmatfuncs.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
~anaconda3libsite-packagesscipysparselinalgmatfuncs.py in expm(A)
590 [ 0. , 0. , 20.08553692]])
591 """
--> 592 return _expm(A, use_exact_onenorm='auto')
593
594
~anaconda3libsite-packagesscipysparselinalgmatfuncs.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.
~anaconda3libsite-packagesscipysparselinalgmatfuncs.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
~anaconda3libsite-packagesscipysparsebsr.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是一种更常见的格式，它也没有这个(。
克朗可能制造bsr
查看CCD_ 5代码，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>
随着CCD_ 22的增长，除大小外，CCD_。
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是一个复杂的计算，它能做的最好的(数值(就是对它进行近似。对于小矩阵，这种近似更容易。这段代码背后有很多数学理论。