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

這很有趣。我可以說的一件事是，該問題是專門針對np.matrix子類的。例如，下面的工作正常：

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

---

挖得更深入回溯，異常被升高在_fragment_2_1在scipy.sparse.linalg.matfuncs.py，具體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被假定爲任一維或一列向量，但對角線的np.matrix對象是一個行向量。這突出了一個更普遍的觀點，即np.matrix通常不如np.ndarray支持得好，所以在大多數情況下最好使用後者。

一個可能的解決辦法是使用np.ravel()拼合diag_T成一維np.ndarray：

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

這似乎解決您遇到的問題，雖然有可能是與np.matrix我的避風港等問題尚未發現。

---

我已經打開了拉請求here。