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#
2import sys, os
3
4os.environ["KMP_DUPLICATE_LIB_OK"] = (
5 "True" # uncomment this line if omp error occurs on OSX for python 3
6)
7os.environ["OMP_NUM_THREADS"] = "1" # set number of OpenMP threads to run in parallel
8os.environ["MKL_NUM_THREADS"] = "1" # set number of MKL threads to run in parallel
9#
10
11#######################################################################
12# example 20 #
13# This example shows how to use the `Lanczos` submodule of the #
14# `tools` module to compute the time evolution of a quantum state #
15# and how to find ground states of hermitian Hamiltonians. #
16#######################################################################
17from quspin.basis import spin_basis_1d
18from quspin.operators import hamiltonian
19from scipy.sparse.linalg import expm_multiply
20from quspin.tools.lanczos import lanczos_full, lanczos_iter, lin_comb_Q_T, expm_lanczos
21import numpy as np
22import scipy.linalg as sla
23
24#
25np.random.seed(17) # set random seed, cf initial state below
26#
27##### Heisenberg model
28L = 20 # system size
29dt = 0.1 # unitary evolution time step
30# basis object
31basis = spin_basis_1d(L, m=0, kblock=0, pblock=1, zblock=1, pauli=False)
32print("\nHilbert space dimension: {}.\n".format(basis.Ns))
33# Heisenberg Hamiltonian
34J_list = [[1.0, i, (i + 1) % L] for i in range(L)]
35static = [[op, J_list] for op in ["xx", "yy", "zz"]]
36H = hamiltonian(static, [], basis=basis, dtype=np.float64)
37#
38##### Lanczos time evolution calculation
39#
40m_evo = 20 # Krylov subspace dimension
41#
42# initial states
43v0 = np.random.normal(0, 1, size=basis.Ns)
44v0 /= np.linalg.norm(v0)
45# make copies to test the lanczos routines independently
46v_expm_multiply = v0.copy()
47v_lanczos_full = v0.copy()
48v_lanczos_iter = v0.copy()
49#
50print("\nchecking lanczos matrix exponential calculation:\n")
51for i in range(100):
52 # compute Lanczos decomposition
53 E_full, V_full, Q_full = lanczos_full(
54 H, v_lanczos_full, m_evo
55 ) # all Lanczos vectors at once
56 E_iter, V_iter, Q_iter = lanczos_iter(
57 H, v_lanczos_iter, m_evo
58 ) # Lanczos vectors as an iterator
59 # evolve quantum state using different routines
60 v_expm_multiply = expm_multiply(
61 -1j * dt * H.static, v_expm_multiply
62 ) # cf tools.expm_multiply_parallel with OMP speedup
63 v_lanczos_full = expm_lanczos(E_full, V_full, Q_full, a=-1j * dt)
64 v_lanczos_iter = expm_lanczos(E_iter, V_iter, Q_iter, a=-1j * dt)
65 # test results against each other
66 np.testing.assert_allclose(v_lanczos_full, v_expm_multiply, atol=1e-10, rtol=0)
67 np.testing.assert_allclose(v_lanczos_iter, v_expm_multiply, atol=1e-10, rtol=0)
68 #
69 print("finished unitary evolution step: {0:d}.".format(i))
70#
71print("\ntime evolution complete.\n")
72#
73###### Lanczos ground state calculation
74#
75# compute exact GS data
76E_GS, psi_GS = H.eigsh(k=1, which="SA")
77psi_GS = psi_GS.ravel()
78#
79###### apply Lanczos
80# initial state for Lanczos algorithm
81v0 = np.random.normal(0, 1, size=basis.Ns)
82v0 /= np.linalg.norm(v0)
83#
84m_GS = 50 # Krylov subspace dimension
85#
86# Lanczos finds the largest-magnitude eigenvalues:
87E, V, Q_T = lanczos_full(H, v0, m_GS, full_ortho=False)
88#
89# check GS energy convergence
90try:
91 # compute difference to exact GS energy value
92 dE = np.abs(E[0] - E_GS[0])
93 assert dE < 1e-10
94except AssertionError:
95 raise AssertionError(
96 "Energy failed to converge |E_lanczos-E_exact| = {} > 1e-10".format(dE)
97 )
98#
99# compute ground state vector
100psi_GS_lanczos = lin_comb_Q_T(V[:, 0], Q_T)
101# check ground state convergence
102try:
103 # compute fidelity of being in exact GS
104 F = np.abs(np.log(np.abs(np.vdot(psi_GS_lanczos, psi_GS))))
105 assert F < 1e-10
106except AssertionError:
107 raise AssertionError(
108 "wavefunction failed to converge to fidelity = {} > 1e-10".format(F)
109 )
110#
111print("\nground state calculation complete.\n")