Wykładanie macierzy w Pythonie

Question

Próbuję wykładać skomplikowaną macierz w Pythonie i mam problemy. Używam funkcji scipy.linalg.expm, i mam dość dziwny komunikat o błędzie, gdy próbuję następujący kod:

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] 

Błąd gdy drugi eksperyment prowadzony jest:

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 

Wydaje się naprawdę dziwne, ponieważ zmieniłem tylko zakres t, którego używałem. Czy to dlatego, że hamiltonian ma przekątną? Ogólnie rzecz biorąc, Hamiltonianie nie będą, ale ja też chcę, żeby działał na przekątnych. Naprawdę nie znam mechaniki expm, więc każda pomoc byłaby bardzo doceniana.

Answer 1

To jest interesujące. Jedno mogę powiedzieć, że problem jest specyficzny dla podklasy np.matrix. Na przykład, następujący działa dobrze:

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

---

Kopanie trochę głębiej w traceback wyjątek jest podniesione w _fragment_2_1 w scipy.sparse.linalg.matfuncs.py konkretnie 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] 

Komunikat o błędzie

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

sugeruje, że exp_diag[k] powinien być skalarem, ale zamiast tego zwraca wektor (i nie można przypisać wektora do X[k, k], który jest skalarem).

Ustawianie przerwania i sprawdzenie kształty tych zmiennych potwierdza:

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

Podstawowym problemem jest to, że exp_diag zakłada się albo 1D lub wektor kolumny, ale przekątna obiektu np.matrix jest wektor wiersza. Podkreśla to ogólniejszą zasadę, że np.matrix jest generalnie mniej obsługiwany niż np.ndarray, więc w większości przypadków lepiej jest używać tego drugiego.

Możliwym rozwiązaniem byłoby wykorzystanie np.ravel() spłaszczyć diag_T do 1D np.ndarray:

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

To wydaje się rozwiązać problem jesteś napotykają, chociaż mogą istnieć inne kwestie związane z np.matrix że haven Jeszcze nie zauważyłem.

---

Mam otwarte żądanie pobrania here.