Esponenziazione della matrice in Python

Question

Sto tentando di esponenziare una matrice complessa in Python e sto incontrando qualche problema. Sto usando la funzione scipy.linalg.expm, e sto avendo un messaggio di errore piuttosto strano quando provo il seguente codice:

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] 

L'errore durante l'esecuzione del secondo esperimento è:

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 

Questo sembra davvero strano dal momento che tutto ciò che ho cambiato era la gamma di t che stavo usando. È perché l'Hamiltoniana è diagonale? In generale, gli Hamiltoniani non lo saranno, ma voglio anche che funzioni per quelli diagonali. Non conosco davvero i meccanismi di expm, quindi qualsiasi aiuto sarebbe molto apprezzato.

Answer 1

Questo è interessante. Una cosa che posso dire è che il problema è specifico della sottoclasse np.matrix. Ad esempio, il seguente funziona bene:

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

---

Scavando un po 'più in profondità la traceback, l'eccezione viene sollevata in _fragment_2_1 in scipy.sparse.linalg.matfuncs.py, in particolare 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] 

Il messaggio di errore

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

mi suggerisce che exp_diag[k] dovrebbe essere uno scalare, ma restituisce invece un vettore (e non è possibile assegnare un vettore a X[k, k], che è uno scalare).

Impostazione di un punto di interruzione ed esaminando le forme di queste variabili conferma questo:

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

Il problema fondamentale è che exp_diag si presume essere o 1D o un vettore colonna, ma la diagonale di un oggetto np.matrix è un vettore di riga. Ciò evidenzia un punto più generale del fatto che lo standard np.matrix è generalmente meno supportato di np.ndarray, quindi nella maggior parte dei casi è preferibile utilizzare quest'ultimo.

Una possibile soluzione potrebbe essere quella di utilizzare np.ravel() per appiattire diag_T in una 1D np.ndarray:

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

Questo sembra risolvere il problema incontrato, anche se ci possono essere altri problemi relativi alla np.matrix che mi porto ancora avvistato.

---

Ho aperto una richiesta pull here.