python求幂_Python中的矩阵求幂

I'm trying to exponentiate a complex matrix in Python and am running into some trouble. I'm using the scipy.linalg.expm function, and am having a rather strange error message when I try 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]

The error when the second experiment is run is:

This works!

Traceback (most recent call last):

File "matrix_exp.py", line 11, in

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 was the range of t I was using. Is it because the Hamiltonian is diagonal? In general, the Hamiltonians won't be, but I also want it to work for diagonal ones. I don't really know the mechanics of expm, so any help would be greatly appreciated.

解决方案

That is interesting. One thing I can say is that the problem is specific to the np.matrix subclass. For example, the following works fine:

h = np.array(hamiltonian)

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

Digging a little deeper into the traceback, the exception is being raised 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]

The error message

X[k, k] = exp_diag[k]

TypeError: only length-1 arrays can be converted to Python scalars

suggests to me that exp_diag[k] ought to be a scalar, but is instead returning a vector (and you can't assign a vector to X[k, k], which is a scalar).

Setting a breakpoint and examining the shapes 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 underlying problem is that exp_diag is assumed to be either 1D or a column vector, but the diagonal of an np.matrix object is a row vector. This highlights a more general point that np.matrix is generally less well-supported than np.ndarray, so in most cases it's better to use the latter.

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

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

This seems to fix the problem you're encountering, although there may be other issues relating to np.matrix that I haven't spotted yet.

I've opened a pull request here.