Lanczos module: time-evolution and ground state search¶
This example demonstrates how to use the Lanczos submodule of the tools module to do time evolvution and ground state search in the Heisenberg model:
where \(S_j\) is the spin-1/2 operator on lattice site \(j\); we use periodic boundary conditions.
The Lanczos decomposition for the \(n\times n\) Hamiltonian matrix is defined as
for a real-valued, symmetric tridiagonal \(m\times m\) matrix \(T=Q^\dagger HQ\), and (in general) a complex-valued \(n\times m\) matrix \(Q\) containing the orthonormal Lanczos vectors in the rows. Here \(m\) is the number of states kept in the Krylov subspace which controls the quality of the “Lanczos compression” of \(H\). We further apply the eigenvalue decomposition \(T=V \mathrm{diag}(E) V^T\) and compute the eigenvectors \(V\) of \(T\) (note that this is computationally cheap for \(m\ll n\)).
Time evolution¶
With this information, we can compute an approximation to the matrix exponential, applied to a state \(|\psi\rangle\) as follows:
If we use \(|\psi\rangle\) as the (nondegenerate) initial state for the Lanczos algorithm, then \(\sum_{j,k}V^T_{ij}Q^\dagger_{jk}\psi_k = \sum_{j}V_{ji}\delta_{0,j} = V_{i,0}\) [by construction, \(\psi_k\) is the zero-th row of \(Q\) and all the rows are orthonormal], and the expression simplifies further. Notice that these lines of thought apply to any matrix function, not just the matrix exponential.
The convergence of this method depends heavily on the function that is being approximated as well as the structure of the matrix. For the matrix exponential there is some literature in math/computer science that discuss error a-priori and a-posteriori bounds for krylow methods. The a-priori bounds are typically far from saturation when performing numerical experiments and the a-posteriori bounds are often impracticle to implement. One can easily check convergence by calculating the lanczos basis of size m but performing the calculation with m and m-1 basis vectors and comparing the results, or by comparing the results of the lanczos calculation to some other method that is implemented, e.g. expm_multiply or expm_multiply_parallel.
In the case for expm_multiply_parallel the convergence is always guaranteed to be machine precision. The tolerance can be slightly controlled by switching between single and double precision floating point types which can often speed up the calculation by a factor of about 1.5. That being said, it is possible to get faster code; however, this requires more memory to store the lanczos vectors in memory during the calculation and often one has to experiment quite a bit to find the optimal time-step and number of lanczos vectors required to beat expm_multiply_parallel.
Ground State Search¶
One of the major uses of the Lanczos method is to find the ground state of a given matrix. It is important to remember that the Lanczos iteration projects out the eigenstates with the largest magnitude eigenvalues of the operator. As such, depending on which eigenvalues one is targeting one might have to transform the operator to make sure that the Lanczos algorithm targets that particular eigenvalue. In the case of the ground state, one either shift the operator by a constant as to make the magitude of the ground state the largest, however, in a lot of cases the ground state already has one of the largest magnitude eigenvalues.
After creating the lanczos basis, QuSpin will return the eigenvalues and vectors of the Krylov sub-space matrix \(T\). If the operator has been transformed to create the Lanczos basis, one should perform the inverse transform of the eigenvalues to get the eigenvalues of the original operator. In the example below the ground state energy is the largest magnitude eigenvalue, hence we do not need to transform the Hamiltonian and likewise, the eigenvalues. The eigenvectors of the Hamiltonian can be constructed by taking linear combinations of the Lanczos basis with coefficients given by the eigenvectors of \(T\).
Script¶
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 | from __future__ import print_function, division
#
import sys,os
os.environ['KMP_DUPLICATE_LIB_OK']='True' # uncomment this line if omp error occurs on OSX for python 3
os.environ['OMP_NUM_THREADS']='1' # set number of OpenMP threads to run in parallel
os.environ['MKL_NUM_THREADS']='1' # set number of MKL threads to run in parallel
#
quspin_path = os.path.join(os.getcwd(),"../../")
sys.path.insert(0,quspin_path)
#######################################################################
# example 20 #
# This example shows how to use the `Lanczos` submodule of the #
# `tools` module to compute the time evolution of a quantum state #
# and how to find ground states of hermitian Hamiltonians. #
#######################################################################
from quspin.basis import spin_basis_1d
from quspin.operators import hamiltonian
from scipy.sparse.linalg import expm_multiply
from quspin.tools.lanczos import lanczos_full,lanczos_iter,lin_comb_Q_T,expm_lanczos
import numpy as np
import scipy.linalg as sla
#
np.random.seed(17) # set random seed, cf initial state below
#
##### Heisenberg model
L = 20 # system size
dt= 0.1 # unitary evolution time step
# basis object
basis = spin_basis_1d(L,m=0,kblock=0,pblock=1,zblock=1,pauli=False)
print("\nHilbert space dimension: {}.\n".format(basis.Ns))
# Heisenberg Hamiltonian
J_list = [[1.0,i,(i+1)%L] for i in range(L)]
static = [[op,J_list] for op in ["xx","yy","zz"]]
H = hamiltonian(static,[],basis=basis,dtype=np.float64)
#
##### Lanczos time evolution calculation
#
m_evo=20 # Krylov subspace dimension
#
# initial states
v0 = np.random.normal(0,1,size=basis.Ns)
v0 /= np.linalg.norm(v0)
# make copies to test the lanczos routines independently
v_expm_multiply = v0.copy()
v_lanczos_full = v0.copy()
v_lanczos_iter = v0.copy()
#
print("\nchecking lanczos matrix exponential calculation:\n")
for i in range(100):
# compute Lanczos decomposition
E_full,V_full,Q_full = lanczos_full(H,v_lanczos_full,m_evo) # all Lanczos vectors at once
E_iter,V_iter,Q_iter = lanczos_iter(H,v_lanczos_iter,m_evo) # Lanczos vectors as an iterator
# evolve quantum state using different routines
v_expm_multiply = expm_multiply(-1j*dt*H.static,v_expm_multiply) # cf tools.expm_multiply_parallel with OMP speedup
v_lanczos_full = expm_lanczos(E_full,V_full,Q_full,a=-1j*dt)
v_lanczos_iter = expm_lanczos(E_iter,V_iter,Q_iter,a=-1j*dt)
# test results against each other
np.testing.assert_allclose(v_lanczos_full,v_expm_multiply,atol=1e-10,rtol=0)
np.testing.assert_allclose(v_lanczos_iter,v_expm_multiply,atol=1e-10,rtol=0)
#
print("finished unitary evolution step: {0:d}.".format(i))
#
print("\ntime evolution complete.\n")
#
###### Lanczos ground state calculation
#
# compute exact GS data
E_GS,psi_GS = H.eigsh(k=1,which="SA")
psi_GS = psi_GS.ravel()
#
###### apply Lanczos
# initial state for Lanczos algorithm
v0 = np.random.normal(0,1,size=basis.Ns)
v0 /= np.linalg.norm(v0)
#
m_GS=50 # Krylov subspace dimension
#
# Lanczos finds the largest-magnitude eigenvalues:
E,V,Q_T = lanczos_full(H,v0,m_GS,full_ortho=False)
#
# check GS energy convergence
try:
# compute difference to exact GS energy value
dE = np.abs(E[0]-E_GS[0])
assert(dE < 1e-10)
except AssertionError:
raise AssertionError("Energy failed to converge |E_lanczos-E_exact| = {} > 1e-10".format(dE))
#
# compute ground state vector
psi_GS_lanczos = lin_comb_Q_T(V[:,0],Q_T)
# check ground state convergence
try:
# compute fidelity of being in exact GS
F = np.abs(np.log(np.abs(np.vdot(psi_GS_lanczos,psi_GS))))
assert(F < 1e-10)
except AssertionError:
raise AssertionError("wavefunction failed to converge to fidelity = {} > 1e-10".format(F))
#
print("\nground state calculation complete.\n")
|