0

I am following this notes to tridiagonalize a symmetric matrix using Householder transformations. The above note has an example worked out in page two. Let me paste my code written in python below.

Although this code did the job for a particular Hamiltonian I was trying, it doesn't work in general. I hope the following code is easy to follow with the added comments. If the matrix is not symmetric, one will get an upper triangular matrix, as I understand.

import numpy as np
from numpy.linalg import matrix_power,norm,eigvals

"Lets first generate a random matrix of dimension d"
d=4
a=np.random.rand(d,d)
"Making it symmetric, my interest is in Hermitian matrices in physics"
h=(a+a.T)/2
ei,ev=np.linalg.eig(h)
hp=h.copy()

for i in range(d-2):  
    "storing the norm of the elements we want to put to zero in the matrix element next to   diagonal"
    if norm(h[i+2:,i])>0: 
        h[i][i+1]=norm(h[i+2:,i])
       
    for j in range(i+2,d):
        h[i][j]=0
   
   
    w= (hp[i]-h[i])
    
    if norm(w)==0:
        w=w
    else:
        w=w/norm(w)
    p=(np.eye(d) - 2*np.outer(w, w))
   
    h=(p@hp@p.T).round(4)
    
    hp=h.copy()
    

print(h)
sreeram pg
  • 1

1 Answers1

0

I could tridiagonalise by following this article in Wikipedia. Let me paste it here since it will be useful for someone else.

import numpy as np
from numpy.linalg import norm,eigvals
def sgn(x):
    if x>=0:
        return 1
    else:
        return -1
"Lets first generate a random matrix of dimension d"
d=10
a=np.random.rand(d,d)
"Making it symmetric, since my interest is in Hermitian matrices in physics"
H=(a+a.T)/2
for i in range(d-1):  
    alpha= -sgn(H[i+1][i])*norm(H[i+1:,i])
    r= np.sqrt(0.5*(alpha**2- H[i+1][i]*alpha))
    v= np.zeros((d,1))
    v[i+1]= (H[i+1][i]-alpha)/(2*r)
    for j in range(i+2,d):
        v[j]= H[j][i]/(2*r)
    w=np.eye(d)-2*v@v.T
    H= w@H@w
print(H)
sreeram pg
  • 1