CODEWITHSUNDEEP
pythonnumpymatrixscipylinear-algebra
LittleByBlue
LittleByBlue

Reputation: 454

scipy.linalg.expm of hermitian is not special unitary

If I have a (real) hermitian matrix, for instance 

H = matrix([[-2. ,  0.5,  0.5,  0. ],
        [ 0.5,  2. ,  0. ,  0.5],
        [ 0.5,  0. ,  0. ,  0.5],
        [ 0. ,  0.5,  0.5,  0. ]])

(This matrix is hermitian; it is the Hamiltonian of a 2-spin Ising chain with coupling to an external field.)

Then there exists a special orthogonal transformation O (preserves length of column and row vectors of a matrix) s.t.

H = O.transpose() @ D @ O

where D is diagonal. For the matrix exponential this leads to

T = expm(1j * H) = O.transpose() @ expm(1j * D) @ O

so all column/row vectors of T must have length 1.

If I use scipy.linalg.expm this property is violated:

In [1]: import numpy as np                                                      

In [2]: from numpy import matrix                                                

In [3]: from scipy.linalg import expm                                           

In [4]: H = matrix([[-2. ,  0.5,  0.5,  0. ], 
   ...:         [ 0.5,  2. ,  0. ,  0.5], 
   ...:         [ 0.5,  0. ,  0. ,  0.5], 
   ...:         [ 0. ,  0.5,  0.5,  0. ]])                                      

In [5]: T = expm(1j * H)                                                        

In [6]: np.sum(np.abs(T[0]))                                                    
Out[6]: 1.6099093263121051

In [7]: np.sum(np.abs(T[1]))                                                    
Out[7]: 1.609909326312105

In [8]: np.sum(np.abs(T[2]))                                                    
Out[8]: 1.7770244703003222

In [9]: np.sum(np.abs(T[3]))                                                    
Out[9]: 1.7770244703003222

Is this a bug in expm or am I making a mistake here?

Upvotes: 0

Views: 294

Answers (1)

Guilhem L.
Guilhem L.

Reputation: 421

you're using the wrong norm. Use

np.sqrt( np.sum( np.abs(T[0])**2 ) )

Or even in a shorter way

np.linalg.norm( T[0] )

Upvotes: 2

Related Questions