Matrix Expression in Python

Question

I am trying to express a complex matrix in Python and I am having problems. I use the scipy.linalg.expm function and I get a rather strange error message when I try to execute the following code:

 import numpy as np from scipy import linalg hamiltonian = np.mat('[1,0,0,0;0,-1,0,0;0,0,-1,0;0,0,0,1]') # This works t_list = np.linspace(0,1,10) unitary = [linalg.expm(-(1j)*t*hamiltonian) for t in t_list] # This doesn't t_list = np.linspace(0,10,100) unitary = [linalg.expm(-(1j)*t*hamiltonian) for t in t_list] 

Error during the second experiment:

 This works! Traceback (most recent call last): File "matrix_exp.py", line 11, in <module> unitary_t = [linalg.expm(-1*t*(1j)*hamiltonian) for t in t_list] File "/usr/lib/python2.7/dist-packages/scipy/linalg/matfuncs.py", line 105, in expm return scipy.sparse.linalg.expm(A) File "/usr/lib/python2.7/dist- packages/scipy/sparse/linalg/matfuncs.py", line 344, in expm X = _fragment_2_1(X, A, s) File "/usr/lib/python2.7/dist- packages/scipy/sparse/linalg/matfuncs.py", line 462, in _fragment_2_1 X[k, k] = exp_diag[k] TypeError: only length-1 arrays can be converted to Python scalars 

This seems really strange, since all I changed is the range of t that I used. Is it because the Hamiltonian is diagonal? In general, there will be no Hamiltonian, but I also want it to work on diagonal ones. I really don't know the mechanics of expm , so any help would be greatly appreciated.

Answer 1

It is interesting. I can only say that the problem is with the np.matrix subclass. For example, the following works just fine:

 h = np.array(hamiltonian) unitary = [linalg.expm(-(1j)*t*h) for t in t_list] 

---

Digging a little deeper into the trace, the exception is expressed in _fragment_2_1 in scipy.sparse.linalg.matfuncs.py , specifically these lines :

 n = X.shape[0] diag_T = T.diagonal().copy() # Replace diag(X) by exp(2^-s diag(T)). scale = 2 ** -s exp_diag = np.exp(scale * diag_T) for k in range(n): X[k, k] = exp_diag[k] 

Error message

  X[k, k] = exp_diag[k] TypeError: only length-1 arrays can be converted to Python scalars 

offers me that exp_diag[k] should be a scalar, but instead returns a vector (and you cannot assign a vector X[k, k] , which is a scalar).

Setting a breakpoint and examining the forms of these variables confirms this:

 ipdb> l 751 # Replace diag(X) by exp(2^-s diag(T)). 752 scale = 2 ** -s 753 exp_diag = np.exp(scale * diag_T) 754 for k in range(n): 755 import ipdb; ipdb.set_trace() # breakpoint e86ebbd4 // --> 756 X[k, k] = exp_diag[k] 757 758 for i in range(s-1, -1, -1): 759 X = X.dot(X) 760 761 # Replace diag(X) by exp(2^-i diag(T)). ipdb> exp_diag.shape (1, 4) ipdb> exp_diag[k].shape (1, 4) ipdb> X[k, k].shape () 

The main problem is that exp_diag is considered either 1D or a column, but the diagonal of the np.matrix object is a row vector. This emphasizes the more general point that np.matrix is generally less well supported than np.ndarray , so in most cases it is better to use the latter.

One possible solution would be to use np.ravel() to smooth diag_T into 1D np.ndarray :

 diag_T = np.ravel(T.diagonal().copy()) 

This seems to np.matrix issue you are facing, although there may be other issues related to np.matrix that I haven't noticed yet.

---

I opened the transfer request here .