Exact diagonalization : tutorial 1
Lanczos algorithm
Lanczos algorithm
Power method and Krylov basis¶
1/ Give an estimate of the number of non-zero matrix elements for a spin chain Hamiltonian of length $L$ and nearest neighbour interactions. Comment.
We do not want to store the Hamiltonian in a full matrix to take advantage of its sparse structure. There are two solutions: either use a sparse storage scheme (only non-zero elements are kept) or we can even recompute matrix elements on-the-fly. Yet, we need an algorithm that allows us to find the ground-state based on the mere operation $$ H|{\psi}\rangle $$ and basics linear algebra for vectors (and not matrices).
A variational argument¶
Consider an Hamiltonian $H$ on a finite-size Hilbert space, with $\{|{\psi_n}\rangle,E_n\}$ its eigenstates/energies. The mean-energy (Rayleigh quotient) of a given vector $|{x}\rangle$ reads: $$E(x) = \frac{\langle{x}|{H}|{x}\rangle}{\langle{x}|{x}\rangle}$$
2/ Show that $E_0 \leq E(x) \leq E_{\text{max}}$, with $E_0$ the ground-state energy and $E_{\text{max}}$ the maximal energy. By considering the gradient of $E(x)$, prove that equality is attained if and only if $|{x}\rangle$ corresponds to the extremal eigenstates.
Therefore, given in inital wave-vector $|{x}\rangle$, we may look for the ground-state in the space spanned by vectors $\{|{x}\rangle, H|{x}\rangle, H^2|{x}\rangle, \ldots, H^p|{x}\rangle \}$ with $p$ a number of iterations. This set of vectors is called the Krylov space.
The power method¶
A first possible algorithm is called the power method and amounts to filter the ground-state $|{\psi_0}\rangle$ by starting from an initial random and normalized state $|{x}\rangle$, and by iteratively applying $H$.
from quantumising import *
L, hx = 6, 0.5
HamH = GeneralOp([ [ -hx, Sx(i)] for i in range(L) ])
HamIsing = GeneralOp([ [ 1.0, SzSz(i,(i+1)%L)] for i in range(L) ])
Ham = HamIsing + HamH
hilbert = IsingHilbert(L=L)
x = randomWf(hilbert)
for i in range(101):
x = Ham*x
x.normalize()
if not i%10: print(i,x.measure(Ham))
# to cross-check results
Ham.createMatrix(hilbert)
Ham.diagonalize(with_eigvec=False)
E0 = np.min(Ham.En)
print("expected E0=",E0)
3/ What is the condition on $|{x}\rangle$ for the algorithm to work. Why the algorithm above is not in general working ?
4/ Let is introduce a constant $C$ to circumvent this issue. How must $C$ be chosen in the algorithm below to ensure convergence ? Is there a better choice for $C$ ? Is that useful in practice?
E0, Emax = np.min(Ham.En), np.max(Ham.En)
C = -Emax
x = randomWf(hilbert)
#x = basisvectorWf(hilbert,3)
for i in range(201):
x = Ham*x - C*x
x.normalize()
if not i%10: print(i,x.measure(Ham),E0,Emax)
The Lanczos algorithm¶
The Lanczos algorithm builds a tridiagonal representation of the Hamiltonian in an orthonormalized basis of the Krylov space. It possesses remarkable properties, such as exponential convergence to the extremal eigenvalues and numerical statibility.
Recursion relations¶
The Lanczos vectors are denoted by $\{|{i}\rangle\}_{i=0,p-1}$ with $p$ the number of iterations. They statisfy the following recursion relations : \begin{align} |{i+1'}\rangle &= H|{i}\rangle - a_i|{i}\rangle - b_{i-1} |{i-1}\rangle \\ |{i+1}\rangle &= |{i+1'}\rangle/b_{i} \\ \end{align} where $b_i = \||{i+1'}\rangle\|$ and $b_{-1}=0$ by convention.
5/ Give the expression of the $a_i$ coefficient as a function of known data at iteration $i$. Show that the matrix representation of the Hamiltonian in the Lanczos subbasis is tridiagonal and gives the expression of its elements.
6/ Implement the Lanczos algorithm to compute the spectrum of the approximate representation of the Hamiltonian.
def lanczosSpectrum(a,b):
tridiag = np.diag(a)
tridiag+= np.diag(b,k=-1)
return np.linalg.eigvalsh(tridiag)
numiteration = 20
numprint = 3
a, b = [], []
u2 = randomWf(hilbert)
for i in range(numiteration):
a.append(u2.measure(Ham))
print(i+1,lanczosSpectrum(a,b)[:numprint],E0)
u3 = Ham*u2 - a[i]*u2
if i>0: u3 -= b[i-1]*u1
b.append(u3.norm())
u3 /= b[i]
u1 = u2
u2 = u3
Computing the ground-state¶
7/ We want to compute the ground-state without keeping in memory all the intermediate Lanczos vectors. How to do that ? Implement the corresponding algorithm.
def lanczosSpectrum(a,b):
tridiag = np.diag(a)
tridiag+= np.diag(b,k=-1)
return np.linalg.eigvalsh(tridiag)
def lanczosGS(a,b):
tridiag = np.diag(a)
tridiag+= np.diag(b,k=-1)
ev, vec = np.linalg.eigh(tridiag)
return vec[:,0]
numiteration = 20
numprint = 3
a, b = [], []
u2 = randomWf(hilbert)
gs = copy(u2) # save the initial random state
for i in range(numiteration):
a.append(u2.measure(Ham))
print(i+1,lanczosSpectrum(a,b)[:numprint],E0)
u3 = Ham*u2 - a[i]*u2
if i>0: u3 -= b[i-1]*u1
b.append(u3.norm())
u3 /= b[i]
u1 = u2
u2 = u3
print("Now computing the ground-state")
# start again from the SAME initial random WF
u2 = copy(gs)
gs[:] = 0.0
gs_krylov = lanczosGS(a,b[:-1])
a,b = [],[]
for i in range(numiteration):
gs += gs_krylov[i]*u2 # recomputing u2 on the fly
a.append(u2.measure(Ham))
u3 = Ham*u2 - a[i]*u2
if i>0: u3 -= b[i-1]*u1
b.append(u3.norm())
u3 /= b[i]
u1 = u2
u2 = u3
Egs = gs.measure(Ham)
VarEgs = scalar(gs,Ham*(Ham*gs))-Egs**2
print("Checking energy and variance of GS : ", Egs, VarEgs)
Practical use¶
Play with the Lanczos algorithm to observe the exponential convergence and the presence of ghosts states due to lack of numerical precision.
We have written a new module quantumising_sparse.py that implements the lanczos calculations above, the sparse storage of the Hamiltonian matrix and also the Arnoldi algorithm which preserves orthogonality (using scipy).
Check that the Arnoldi algorithm cures the ghosts states issue.
from quantumising_sparse import *
L, h = 10, 0.3
def IsingHam(L=4,h=0.0,J=1.0):
HamH = GeneralOp([ [ -h, Sx(i)] for i in range(L) ])
HamIsing = GeneralOp([ [ -J, SzSz(i,(i+1)%L)] for i in range(L) ])
return HamIsing + HamH
hilbert = IsingHilbert(L=L)
# without storing the matrix
Ham = IsingHam(L,h)
En = Ham.Energies(hilbert,num=2,algorithm="Lanczos")
print(En)
# with storing the matrix
Ham.createMatrix(hilbert,sparse=True)
En = Ham.Energies(hilbert,num=2,algorithm="Lanczos")
print(En)
Fulldiag or Arnoli¶
Ham = IsingHam(L,h)
En = Ham.Energies(hilbert,num=2,algorithm="Fulldiag")
print(En)
Ham = IsingHam(L,h)
En = Ham.Energies(hilbert,num=2,algorithm="Arnoldi") # recommended, default algorithm
print(En)
# without storing the matrix
Ham = IsingHam(L,h)
psi0 = Ham.Groundstate(hilbert,algorithm="Lanczos")
print(psi0.measure(Ham))
# with storing the matrix
Ham.createMatrix(hilbert,sparse=True)
psi0 = Ham.Groundstate(hilbert,algorithm="Lanczos")
print(psi0.measure(Ham))
Fuldiag or Arnoldi¶
Ham = IsingHam(L,h)
psi0 = Ham.Groundstate(hilbert,algorithm="Fulldiag")
print(psi0.measure(Ham))
Ham = IsingHam(L,h)
psi0 = Ham.Groundstate(hilbert,algorithm="Arnoldi") # recommended
print(psi0.measure(Ham))
Downloading the modules¶
modules written in Python 3¶
- http://lptms.u-psud.fr/membres/groux/Test/ED/quantumising.py
- http://lptms.u-psud.fr/membres/groux/Test/ED/quantumising_sparse.py