Exact diagonalization : tutorial 1
Lanczos algorithm
Lanczos algorithm
Power method and Krylov basis¶
1/ Give an estimate of the number of non-zero matrix-element 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 in which only non-zero elements are stored or we can even recompute matrix elements on-the-fly on demand. 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$ in a finite-size Hilbert space and the mean-energy in a given vector $|{x}\rangle$: $$E(x) = \frac{\langle{x}|{H}|{x}\rangle}{\langle{x}|{x}\rangle}$$
2/ Show that minimizing $E(x)$ suggests to look for the ground-state in a certain direction.
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 starting from a random initial state $|{x}\rangle$ and iteratively applying $H$.
from quantumising import *
L, hx = 4, 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(201):
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.
1/ 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.
2/ 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):
...
...
...
...
Computing the ground-state¶
3/ We want to compute the ground-state without keeping in memory all the intermediate Lanczos vectors. How to do that ? Implementation 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 = deepcopy(u2) # save the initial random state
...
...
...
...
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¶
- http://lptms.u-psud.fr/membres/groux/Test/ED/quantumising.py
- http://lptms.u-psud.fr/membres/groux/Test/ED/quantumising_sparse.py