Pythonでの行列の指数化

Question

私はPythonで複雑な行列を指数化しようとしていますが、何か問題があります。私はscipy.linalg.expm機能を使用しています、と私は、次のコードをしようとすると、かなり奇妙なエラーメッセージを持っています：

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] 

第二の実験が実行されたエラーは次のとおりです。

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 

これは本当にそうです私が変更したすべてがtの範囲だったので奇妙です。ハミルトニアンは対角であるからでしょうか？一般的に、ハミルトニアンはそうではありませんが、私はそれが対角線のものでも働きたいと思っています。私は実際にexpmの仕組みを知らないので、どんな助けでも大歓迎です。

Answer 1

これは興味深いことです。私が言うことの1つは、問題がnp.matrixサブクラスに固有であるということです。たとえば、次のように正常に動作します：トレースバックに少し深く掘り

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

---

は、例外がscipy.sparse.linalg.matfuncs.pyに_fragment_2_1に提起され、特に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] 

エラーメッセージ

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

は私にそれを示唆していますexp_diag[k]はスカラーでなければなりませんが、代わりにベクトルを返します（ベクトルをX[k, k]に割り当てることはできません。これはスカラーです）。

ブレークポイントを設定し、これらの変数の形状を調べ、これを確認：

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 
() 

根本的な問題はexp_diagが1D又は列ベクトルのいずれかであると仮定されていることであるが、np.matrixオブジェクトの対角線は、Aであります行ベクトル。これは、np.matrixが一般的にnp.ndarrayよりよくサポートされていないというより一般的な点を強調しています。したがって、ほとんどの場合、後者を使用する方が良いです。

1D np.ndarrayにdiag_Tを平らにするnp.ravel()を使用する一つの可能​​な解決策は、次のようになります。

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

np.matrix I避難所に関連する他の問題があるかもしれませんが、これは、あなたが遭遇している問題を解決しているようですまだ見つかりませんでした。

---

私はプル要求hereを開きました。