CODEWITHSUNDEEP
pythonscipyexponentiation
TheDorkSide
TheDorkSide

Reputation: 37

Why does matrix exponential not work beyond certain size?

I'm having trouble with the matrix exponential calculation using 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)))

where, N can be any integer, number is an NxN diagonal matrix with diagonal elements diag{0,1,2, ... N-1}, and sp.kron is a kronecker product in scipy.sparse. This works fine only for N<6. When I try to run the code with N=6

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

This should be a fairly simple code, but I don't know why it is giving the following error:

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:

Upvotes: 0

Views: 538

Answers (1)

hpaulj
hpaulj

Reputation: 231385

According to the traceback T[k, k+1] doesn't work because T is a bsr format sparse matrix, which does not implement indexing. (coo is a more common format that doesn't have this either).

kron may make bsr

Looking at the sp.kron code, 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)

So in 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>

testing

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>

changing to csc to get avoid this warning:

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>

So simply providing a bsr to expm doesn't cause your error. Looks like we'd have to examine expm what else happening. I looked at this function (and MATLAB's) years ago. It uses a pade approximation that includes a inv (which is spsolve(I,A)). It's a complex function that tries different things, including different Pade orders.

So you'll have to tell us more about this number and the nature of the kron() result. None of my guesses reproduce your error.

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

upper triangular

Correction, the traceback tells us that it detect that your matrix is upper triangular:

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

so there's more code to trace through.

In any case doing the tocsc conversion before passing the matrix to expm might solve the problem:

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

test small upper triangular array

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>

Your error:

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

with csc correction:

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>

with 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>

No significant difference in the kron result except size as N grows.

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

Specifying a csc format in kron avoids that error (we can ignore this efficiency warning):

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>

Why N=6 gives this warning and not small N probably has to do with which Pade order it had to try. Keep in mind that expm is a complex calculation, and the best it can do (numerically) is approximate it. That approximation is easier with small matrices. There's a lot of math theory behind this code.

Upvotes: 3

Related Questions