arXiv:1901.07653v2 [quant-ph] 15 Mar 2019
Quantum Imaginary Time Evolution, Quantum Lanczos, and Quantum Thermal Averaging
Mario Motta,
1,
∗
Chong Sun,
1
Adrian T. K. Tan,
2
Matthew J. O’Rourke,
1
Erika Ye,
2
Austin J. Minnich,
2
Fernando G. S. L. Brand˜ao,
3, 4
and Garnet Kin- Lic Chan
1,
†
1
Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, CA 91125, USA
2
Division of Engineering and Applied Sciences, California Institute of Technology, Pasadena, CA 91125, USA
3
Institute for Quantum Information and Matter, California Institute of Technology, Pasadena, CA 91125, USA
4
Google LLC, Venice, CA 90291, USA
An efficient way to compute Hamiltonian ground-states on a quantum computer stands to impact many prob-
lems in the physical and computer sciences, fr om quantum simulation to machine learning. Existing techniques,
such as phase estimation and variational algorithms, display potential disadvantages, including requirements for
deep circuits with ancillae and high-dimensional optimization. Here we describe the quantum imaginary time
evolution and quantum Lanczos algorithms, analogs of classical algorithms for ground (and excited) states, but
with exponentially reduced space and time requirements per iteration, and avoiding deep circuits with ancillae
and high-dimensional optimization. We discuss quantum imaginary time evolution as a natural subroutine to
generate Gibbs averages through an analog of minimally entangled typical thermal states. We implement these
algorithms with exact classical emulation and prototype circuits on the Rigetti quantum virtual machine and
Aspen-1 quantum processing unit, demonstrating the power of quantum elevations of classical algorithms.
An important application for a quantum computer is to
compute the ground-state Ψ of a Hamiltonian
ˆ
H [
1, 2]. This
arises in simulations, for example, of the electronic structure
of molecules and materials, [
3–6] as well as in more general
optimization problems. While efficient ground-state determi-
nation cannot be guaranteed for all Hamiltonians, as this is a
QMA-hard problem [
7], several heuristic quantum algorithms
have been proposed, inc luding adiabatic state pre paration with
quantum phase estimation [
8, 9] (QPE) a nd quantum-classical
variational algorithms, such as the quantum approximate op-
timization algorith m [
10–12] and variational quantum e igen-
solver [13–15]. Despite many advances, these algorith ms also
have potential disadvantages, especially in the context of near-
term quantum computing architectures with limited quantum
resources. For example, phase estimation produces a ne arly
exact eigenstate, but appears impractical without error cor-
rection, while variational algorithms, although somewhat ro-
bust to cohere nt errors, are limited in accuracy for a fixed
Ansatz, and involve a high-dimensional noisy classical opti-
mization [
16].
In classical simulations, different strategies are employed to
numerically determine nearly exact ground-states. One pop-
ular approac h is imaginary-time evolution, which expresses
the ground-state as the long-time limit of the imaginary-
time Schr¨odinger equation −∂
β
|Φ(β)i =
ˆ
H|Φ(β)i, |Ψi =
lim
β→∞
|Φ(β)i
kΦ(β)k
(for hΦ(0)|Ψi 6= 0). Unlike variational al-
gorithms with a fixed An satz, imaginary-time evolution al-
ways converges to the ground-state, as distinguished fro m
imaginary-time Ansatz optimization [
17]. Another comm on
algorithm is the iterative La nczos algorithm [
18] and its vari-
ants. The Lanczos iteration constructs the Hamiltonian ma-
trix H in a Krylov subspace {|Φi,
ˆ
H|Φi,
ˆ
H
2
|Φi. . .}; diag-
onalizing H yields a variational estimate of the ground-state
which tends to |Ψi for a large number of iterations. For an N-
qubit Ha miltonian, the classical complexity of imaginary time
evolution and La nczos algorithm scales as ∼ exp (O(N)) in
space and time. Exponential space comes from sto ring Φ(β)
or the Lanczos vector, wh ile exponential time comes from the
cost of Hamiltonian multiplicatio n
ˆ
H|Φi, as well as, in princi-
ple, though not in practice, the N-dependence of the number
of propagation steps or Lanczos iterations. Thus it is natu-
ral to consider quantum versions of these algor ithms that can
overcome the expo nential bottlenecks.
Here we descr ibe the quantum imaginary time evolution
(QITE) and the quan tum Lanc z os (QL anczos) algorithms to
determine g round-states (and excited states in the case of
QLanczos) on a quantum compu ter. As we show, under well
defined assumptions, these use exponentially reduced spa ce
and time per propagation step or iteration compared to their
direct classical counterparts. They also offer ad vantages over
existing ground-state quantum algorithms as they do not use
deep circuits and ar e guaranteed to converge to the ground-
state without non-linear optimization. We further d escribe
inexact QITE and QLanczos algorithms that present a hierar-
chy of approxim ations to apply within a limited computational
budget. A crucial common component is the efficient imple-
mentation of the non-H e rmitian operation of an imaginary-
time step e
−∆τ
ˆ
H
(for sma ll ∆τ ) assuming a finite correlation
length in th e state. Non-Hermitian operations are not natural
on a quantum computer and are usually achieved using an-
cillae and p ostselection, but w e describe how to implement
imaginary time evolution on a given state without these re-
sources. The lack of ancillae and complex circuits make QITE
and QLanczos potentially suitable for near-term qua ntum ar-
chitectures. Using the QITE algorithm, we show how to sam-
ple from thermal (Gibbs) states, also w ithout deep c ircuits or
ancillae as is usually the case, via a quantum analo g of the
minimally entangled typical thermal states (QMETTS) algo-
rithm [
19, 20]. We demonstrate the algorithms on spin and
fermion ic Hamiltonians (short- and lon g-range spin and Hub-
bard models, MAXCUT optimization, and dihy drogen m in-
imal molecular m odel) using exact classical emulation, and
demonstra te pro of-of-co ncept implementations on the Rigetti
quantum virtual machine (QVM) and Aspen-1 quantu m pro-
2
cessing units (QPUs).
Quantum Imaginary-Time Evolution. Define a geometric
k-local Hamiltonian
ˆ
H =
P
m
ˆ
h
m
(where each term
ˆ
h
m
acts
on at most k neighbouring qubits on an underlying graph) and
a Trotter decomposition of the corresponding imaginary-time
evolution,
e
−β
ˆ
H
= (e
−∆τ
ˆ
h
1
e
−∆τ
ˆ
h
2
. . .)
n
+ O(∆τ) ; n =
β
∆τ
(1)
applied to a state |Ψi. After a single Trotter step, we have
|Ψ
′
i = e
−∆τ
ˆ
h
m
|Ψi. (2)
The b a sic idea is th a t the norm alized state |
¯
Ψ
′
i = |Ψ
′
i/kΨ
′
k
is generate d from |Ψi by a unitary operator e
−i∆τ
ˆ
A[m]
acting
on a neighbourhood of the qu bits ac te d on by
ˆ
h
m
, where
ˆ
A[m]
can be determined from tomography of |Ψi in this neighbour-
hood up to controllable errors. This is illustrated by the sim-
ple example wher e |Ψi is a produ ct state. The squared norm
c = kΨ
′
k
2
can b e calculated from the exp ectation value of
ˆ
h
m
, requiring me a surements over k qubits,
c = hΨ|e
−2∆τ
ˆ
h[m]
|Ψi = 1 − 2∆τhΨ|
ˆ
h
m
|Ψi + O(∆τ
2
)
(3)
Because |Ψi is a product state, |Ψ
′
i is obtain ed applying the
unitary operator e
−i∆τ
ˆ
A[m]
also o n k qubits.
ˆ
A[m] can be
expanded in terms of a n operator basis, e.g. the Pauli basis
{σ
i
} on k qubits,
ˆ
A[m] =
X
i
1
i
2
...i
k
a[m]
i
1
i
2
...i
k
σ
i
1
σ
i
2
. . . σ
i
k
. (4)
Up to O(∆τ), the co efficients a[m]
i
1
i
2
...i
k
are defined by the
linear system Sa[m] = b where th e elements of S and b are
expectation values over k qubits,
S
i
1
i
2
...i
k
,i
′
1
i
′
2
...i
′
k
= hΨ|σ
†
i
1
σ
†
i
2
. . . σ
†
i
k
σ
i
′
1
σ
i
′
2
. . . σ
i
′
k
|Ψi
b
i
1
i
2
...i
k
= −i c
−
1
2
hΨ|σ
†
i
1
σ
†
i
2
. . . σ
†
i
k
ˆ
h[m]|Ψi (5)
In general, S ha s a null space; to ensure a[m] is real, we min-
imize kc
−1/2
Ψ
′
− (1 − i∆τ
ˆ
A[m])Ψk
2
w.r.t. real variations
in a[m] (see SI). Becau se the solution is d etermined from a
linear problem, there are no local minima.
In this simple case, the normalized result of the imaginary
time evolution step could be represented by a unitary update
over k qubits, because |Ψi had correlation length z ero. After
the initial step, this is no longer the case. However, for a more
general |Ψi with finite correlations over at most C qubits (i.e.
correlations between observables separated by distance L are
bounded by exp(−L/C)), |Ψ
′
i can be generate d by a unitary
acting on a do main of width at m ost O(C) qubits surrounding
the qubits acted on by
ˆ
h
m
(this follows f rom Uhlmann’s the-
orem [
21]; see SI). The unitary e
−i∆τA[m]
can then be deter-
mined by measurements and solving the least squares problem
in this domain (Fig. 1). For example, for a ne a rest-neighbo r
local Hamiltonian on a d-dimension cubic lattice, the domain
size D is bounded by O(C
d
). In many p hysical systems, we
expect the maximum correlation len gth throughout the Trot-
ter steps to increase with β and saturate for C
max
≪ N [
22].
Fig.
1 shows the mutual information between qubits i and j
as a function of imaginary time in the 1D an d 2D ferromag-
netic transverse field Ising models computed by tensor net-
work simulation (see SI), dem onstrating a monotonic increase
and clear saturation.
The above replacement of ima ginary time evolution steps
by unitary up dates ca n be extended to more general Hamilto-
nians, such as ones with long-range interactions and fermionic
Hamiltonians. For example, for a Hamiltonian with long-
range pairwise terms, the action of e
−∆τ
ˆ
h[m]
(if
ˆ
h[m] acts o n
qubits i and j) can be emulated by a unitary constructed in the
neighborhoods of i and j, over a domain of (2C log(1/δ))
k
sites (see SI). Th e assumption of finite correlation length,
however, is less natural fo r such Hamiltonians. For fermions,
the lo c ality of the corresponding qubit Ham iltonian depends
on the spin mapping. In principle, a geom etric k-local
fermion ic Ham iltonian can be mapped to a geometric local
qubit Hamiltonian [
23], allowing the above techniques to be
directly applied. Alternatively, we conjecture that by using a
fermion ic unitary, where the Pauli basis in Eq. (
4) is replaced
by the fermio nic operator basis {1, ˆa, ˆa
†
, ˆa
†
ˆa}, the unitary up-
date can be constructed over a domain size D ∼ O(C
d
) where
C is the fermio nic correlation len gth.
Cost of QITE. The number of measureme nts and classical
storage at a given time step ( star ting propagation from a prod-
uct state) is bounded by exp(O(C
d
)) (with C the co rrelation
length at that time step), since each unitary at that step acts
on at most O(C
d
) sites; classical solution of the least squares
problem has a similar scaling exp(O(C
d
)), as does the syn-
thesis and application as a quantum circuit (comp osed of two-
qubit gate s) of the unitary e
−i∆τA[m]
. Thus, space and time
requirements are bounded by expon entials in C
d
, but are poly-
nomial in N when one is interested in a local approx imation
of the state (or quasi-polynomial for a global ap proximation);
the polynomial in N comes from the number of terms in H;
see SI for details).
The exponential dependence on C
d
can be greatly reduced
in many cases. Suppose the Hamiltonian A[m] of the unitary
update has a loc a lity structure, i.e. it is (approximately) a p-
local Hamilton ian (i.e. in E q. (
4), all a[m]
i
1
...i
k
coefficients
are zero except for those where at most p of the σ
i
operators
are different from the identity). Then the cost of tomo graphy
becomes only C
O(d p)
, while the cost of finding and imple-
menting the unitary is O(pC
d
T
e
), with T
e
the cost to com-
pute one entry of A[m] [
24]. If we assume further that A[m]
is geometric local, the cost of tomogr aphy is reduced further
to O(pC
d
). However, even if C is too large to construc t the
unitaries exactly, we can still run the algorithm a s a h e uris-
tic, truncating the unitary updates to domain sizes that fit the
computational budget. This gives the inexact QITE algorithm,
described further below.
Comparison to classical implementations. Compared to a di-
3
FIG. 1: (color online) (a) Schematic of the QITE algorithm. Top: i maginary-time evolution under a geometric k-local operator
ˆ
h[m] can be
reproduced by a unitary operation acting on D > k qubits. Bottom: exact imaginary-time evolution starting from a product state requires
unitaries acting on a domain D that grows with correlations. (b,c) Left: mutual information I(i, j) between qubits i, j as a function of distance
d(i, j) and imaginary time β, for a 1D (b) and a 2D (c) FM transverse-field Ising model, wi th h = 1.25 (1D) and h = 3.5 (2D). I(i, j) saturates
at longer times. Right: relative error in the energy ∆E and fidelity F = |hΦ(β)|Ψi|
2
between the finite-time state Φ(β) and infinite-time state
Ψ as a function of β. The noise in the 2D fidelity error at large β arises from the approximate nature of the algorithm used. See SI for details.
rect classical implementation of ima ginary time evolution, the
cost of a QITE time-step (for boun ded c orrelation length C) is
linear in N in spa ce and polynomial in N in time, thus giving
an exponential reduction in space and time. We can also com-
pare to other classical algorithms. As QITE defines a quantum
circuit for the imaginary time evolution, we could attempt to
use it for a faster classical simulation. If we are only interested
in local observables, we can apply the circuit in the Heisen-
berg picture in a classical emulation. However, this gives an
extra exponential dependence on th e number of previous time-
steps: After the unitaries associated to (e
−∆τ
ˆ
h
1
e
−∆τ
ˆ
h
2
. . .)
l
have been applied, the cost of applying the next unitar y scales
as exp(O(lD)), with D the domain size of the unitar ie s, in-
stead of exp(O(D)) in QITE . Alternatively, if |Ψi is rep-
resented by a ten sor network in a classical simulation, then
e
−∆τ
ˆ
h[m]
|Ψi can be represented as a classical tensor network
with increased bond dimension [
25, 26]. However, the bond
dimension will scale as exp(O(lD)). Apart from the extra
exponential dependence on l, a further potential drawback in
this approach is that we cannot guarantee contracting the re-
sulting classical te nsor network for an observable is efficient;
it is a #P-hard problem in the worst case in 2D (and even in
the average case for Gaussian distributed tensors) [
27, 28].
Finally, we can compare QITE with bounded C with the clas-
sical heuristic of truncating the problem size a t the correlation
length C
0
of the ground-state a nd solving by exact diagonal-
ization, which can be do ne in time exp(O(C
d
0
)) in d spatial di-
mensions. While this is a competitive strategy in many cases,
it may not converge to the correct ground-state when there is
frustration in the Hamilton ia n, for example in glassy models.
Inexact QITE. Given limited resou rces, for example on near-
term devices, we can choose to measure and construc t the uni-
tary over a domain D smaller than induced b y correlations, to
fit the computational budget. For example, if D = 1, this
gives a mean-field approximation of the im aginary tim e evo-
lution. While the unitary is n o longer an exact representation
of th e imagina ry time evolution, there is no issue of a loca l
minimum in its construction, although the energy is no longer
guaran teed to decrease at every step. In this case, one can
apply inexact imaginary time evolution until the energy stops
decreasing; the energy will still be a variational upper bound.
One ca n also use the quantum Lanczos algorithm, de scribed
later.
QITE experiments. To illustra te the QITE a lgorithm, we
have carried o ut exact classical emulations (assuming per-
fect expec tation values and perfect gates) for several Hamil-
tonians: short-range 1D Heisenberg; 1D AFM transverse-
field Ising; long-range 1D Heisenberg with spin-spin coupling
J
ij
= (|i − j| + 1)
−1
; i 6= j; 1D Hubbard at half-filling
(mapped b y Jordan-Wigner transform ation to a spin model) ;
a 6-qubit MAXCUT [
10–12] instance, and a minimal basis 2-
qubit dihydrogen molecular Hamiltonian [
29]. To assess the
feasibility of implem e ntation on near-term quantum devices,
we have also carried out n oisy classical emulation (sampling
expectation values and with an error model) using the Rigetti
quantum virtual machine (QVM) and a physical simulation
using the Rigetti Aspen-1 QPUs, for a single qubit field model
(2
−1/2
(X + Z))[
30] and a 1D AFM transverse-field Ising
model. We carry out QITE using different fixed domain sizes
D for the unitary or ferm ionic u nitary (see SI for descriptions
of simulations and mod els).
Figs.
2 and 3 show the energy obtained by QITE as a func-
tion of β and D for the various models. As w e increase D, the
asymptotic ( β → ∞) energies rapidly converge to the exact
4
-18
-16
-14
-12
-10
E(β)
(a)
0.00 2.00 4.00 6.00 8.00
β
0.00
0.25
0.50
0.75
1.00
F (Φ
β
, Ψ)
(b)
D = 2
D = 4
D = 6
D = 8
exact
-36
-32
-28
-24
-20
E(β)
(c)
D = 2
D = 2, QLanczos
exact
0.00 1.00 2.00 3.00 4.00
β
-36
-32
-28
-24
-20
E(β)
(d)
D = 4
D = 4, QLanczos
exact
-1.2
-0.7
-0.2
0.3
0.8
E(β)
(e)
QPUs, QITE
QPUs, QLanczos
exact
0.00 0.75 1.50 2.25 3.00
β
-2.0
-1.0
0.0
1.0
2.0
E(β)
(f)
QVM, QITE
QPUs, QITE
QVM, QLanczos
QPUs, QLanczos
exact
FIG. 2: Left: QITE energy E(β) (a) and fidelity F (b) between finite-time state Φ(β) and exact ground state Ψ as function of imaginary time
β, for a 1D 10-site Heisenberg model, showing the convergence wi th increasing unitary domains of D = 2 − 8 qubits. Middle: QITE ( dashed
red, dot-dashed green lines) and QLanczos (solid red, solid green lines) energies as function of imaginary time β, for a 1D Heisenberg model
with N = 20 qubits, using domains of D = 2 (c) and 4 qubits ( d), showing improved convergence of QLanczos over QITE. Black line is the
exact ground-state energy/fidelity. Right: QITE and QLanczos energy E(β) as a function of imaginary time β for (e) 1-qubit field model using
the QVM and QPU (qubit 14 on Aspen-1), (f) 2-qubit A FM transverse field Ising model using the QVM and QPU (qubit 14, 15 on Aspen-1).
Black line is the exact ground-state energy (see SI for details).
ground-state. For small D, the inexact QITE tracks the exact
QITE for a time until the correlation length exceeds D. Af-
terwards, it may go down or up. The non-monotonic behavior
is strongest f or small domains; in the MAXCUT example, the
smallest domain D = 2 gives an oscillating energy; the first
point at which the energy stops decreasing is a reasonable es-
timate of the ground-state energy. In all models, increasing
D past a maximum value (less than N) no longer affects the
asymptotic e nergy, showing that the c orrelations have satu-
rated (this is true even in the MAXCUT instance).
Figs.
2e and 2f show the results of running the QITE al-
gorithm on Rigetti’s QVM and Aspen-1 QPUs for 1- and 2-
qubits, respec tively. The e rror bars are due to gate, readout,
incohere nt and cross-talk errors. Sufficient samples were used
to ensure that sampling error is negligib le . Encouragingly for
near-term simulations, despite these errors it is possible to
converge to a ground-state energy close to th e exact energy for
the 1-qubit case. This result reflects a robustness that is some-
times informally observed in imaginary time evolution algo-
rithms in which the ground state energy is appro ached even
if the imaginary time step is not perfectly implemented. In
the 2-qubit case, although the QITE energy converges, there
is a systematic shift which is reproduced on the QVM using
ava ilable noise parame te rs for readout, d ecoherence and de-
polarizing noise [
31]. Remaining discrepancies between the
emulator and hardware are likely attributable to cross-talk be-
tween parallel gates not included in the noise model (see SI).
However, reducing decoherence and depolarizing errors in the
QVM o r using different sets of qubits with improved noise
characteristics (see SI) all lead to improved convergence to
the exact ground-state e nergy.
Quantum Lanczos algorithm. G iven the QITE subroutine,
0 1 2 3 4
β
-5
-4
-3
-2
-1
E(β)
(a)
D = 2
D = 4
D = 6
exact
0 2 4 6 8
β
(b)
D = 2
D = 4
exact
0 4 8 12 16
β
0.00
0.25
0.50
0.75
1.00
P (C = C
max
)
(c)
D = 2
D = 4
D = 6
exact
0.0 1.0 2.0 3.0
R [
˚
A]
(d)
β = 0
β = 1
β = 2
β = 3
exact
-5
-4
-3
-2
-1
E(β)
-1.2
-0.9
-0.6
-0.3
0.0
0.3
E(β) [E
Ha
]
FIG. 3: (a) QITE energy E(β) as a f unction of imaginary time β
for a 6-site 1D long-range Heisenberg model, for unitary domains
D = 2 − 6; (b) a 4-site 1D Hubbard model with U/t = 1, for
unitary domains D = 2, 4. (c) Probability of MAXCUT detection,
P (C = C
max
) as a function of imaginary time β, for the 6-site
graph in the panel. (d) QIT E energy for the H
2
molecule in the STO-
6G basis as a function of bond-length R and β. Black line is the
exact ground-state energy/probability of detection.
we now consider how to formulate a quantum Lanczos al-
gorithm, which is an especially economic al realizatio n of a
quantum subspace method [32, 33]. An im portant practical
motivation is that th e Lanczos alg orithm typically converges
much more quickly than imaginary time evolution, and often
in physical simulations only tens of iterations are needed to
converge to good precision. In addition, Lanczos provides a
5
0.0 1.0 2.0 3.0 4.0
β
-1.80
-1.35
-0.90
-0.45
0.00
h
ˆ
Hi(β)
(a)
D = 2
D = 3
D = 4
D = 5
exact
1.0 2.0 3.0 4.0
β
-1.10
-1.00
-0.90
-0.80
-0.70
h
ˆ
Hi(β)
(b)
QVM, noiseless
QVM, with noise
QPUs
exact
1.0 2.0 3.0 4.0
β
-1.80
-1.60
-1.40
-1.20
-1.00
h
ˆ
Hi(β)
(c)
Emulator, noise
Emulator, noise
exact
FIG. 4: Left: Thermal (Gibbs) average h
ˆ
Hi at temperature β from QMETTS for a 1D 6-site Heisenberg model (exact emulation). Black line
is the exact thermal average without sampling error. Middle, Right: Thermal average h
ˆ
Hi at temperature β from QMETTS for (b) a 1 qubit
field model using QVMs and QPUs, and (c) 2 qubit AFM transverse field Ising model using QVM.
natural way to compute excited states. Con sider the sequence
of imaginary time vectors |Φ
l
i = e
−l∆τ
ˆ
H
|Φi, l = 0 , 1, . . . n,
where c
l
= kΦ
l
k. In QLanczos, we consider the vectors af-
ter even numbers of time steps |Φ
0
i, |Φ
2
i. . . to form a ba-
sis fo r the ground- state. (SI describes the equiva le nt treat-
ment in terms of normalized imaginary time vectors). These
vectors define an overlap matrix whose elements can be com-
puted entirely f rom norms, S
ll
′
= hΦ
l
|Φ
l
′
i = c
2
(l+l
′
)/2
, where
c
(l+l
′
)/2
is the norm of another integer tim e step vector, and
the overlap matrix elements for n/2 vecto rs can be accumu-
lated for free after n steps of time evolution. The Hamilto-
nian matrix elements satisfy the identity H
ll
′
= hΦ
l
|
ˆ
H|Φ
l
′
i =
hΦ
(l+l
′
)/2
|
ˆ
H|Φ
(l+l
′
)/2
i. Although the Hamiltonian h a s ∼ n
2
matrix elem ents, there are only ∼ n unique elements, and im-
portantly, each is a simple expectation value of the energy
during the imaginary time evolution. This ec onomy of ma-
trix elements is a property shared with the classical Lanczos
algorithm . Whereas the c la ssical Lanczos iteration builds a
Krylov space in powers of
ˆ
H, QLanczos builds a Krylov sp a ce
in powers of e
−2∆τ
ˆ
H
; in the limit of small ∆τ these Krylov
spaces ar e identical. Diagonaliza tion of th e QLancz os Hamil-
tonian matrix is guaranteed to give a groun d-state energy that
is lower tha n that of the last imaginary time vector Φ
n
(while
higher roots approximate excited states). Thus, as long as one
is willing to take measur e ments of the energy during th e imag-
inary time evolution process, one can use QLanczos to gener-
ate an im proved ground state (or excited states).
With a limited computational budget, we can use inexact
QITE to generate Φ
l
, Φ
′
l
. However, in th is case the above
expressions for S
ll
′
and H
ll
′
in terms of expectation values are
no longer exactly satisfied which can create numerical issues
(e.g. the overlap may no longer be positive). To handle this as
well as errors due to noise and sampling in real experiments,
the QLanczos algorithm needs to be stabilized by ensuring
that successive vectors are not nearly linearly depen dent (see
SI).
We demonstrate the QLanczo s algorithm using classical
emulation on the 1D Heisenberg Hamilto nian, as used for
the QITE algorithm in Fig.
2 (see SI). Using exact QITE
(large domains) to generate matrix elements, quan tum Lanc-
zos converges much more rapidly than imaginary time evo-
lution. Using inexact QITE (sm a ll domains), co nvergence is
usually faster and also reaches a lower en ergy. We also as-
sess the feasibility of QL a nczos in presence of noise, using
emulated noise on the Rigetti QVM as well as on the Rigetti
Aspen-1 QPUs. In Fig.
2, we see that QLanczos also provides
more rapid convergence than QITE with both noisy classical
emulation as well as on the physical d evice for 1 an d 2 q ubits.
Quantum thermal averages. The QITE subr outine can be
used in a range of other alg orithms. For example , we discuss
how to compute thermal averages Tr
ˆ
Oe
−β
ˆ
H
/Tr
e
−β
ˆ
H
us-
ing imaginary time evolution. Several procedure s have been
proposed for qu antum thermal averaging, ranging from gen-
erating the finite-temperature state explicitly by eq uilibration
with a bath [
34], to a quantum analog o f Metropolis sam-
pling [35] that re lies on phase estimation, as we ll as meth-
ods based on ancilla based Hamilto nian simulation with post-
selection [
36] and approaches based on r ecovery maps [37].
However, given a method for imaginary time evolution, one
can generate thermal averages of observables without any
ancillae or d e ep circu its. This can be done by adapting to
the quantum setting the classical minimally entangled typical
thermal state (METTS) algorithm [
19, 20], which gen e rates
a Markov cha in from which the thermal average can be sam-
pled. The Q METTS algorithm can be car ried out as follows (i)
start from a product state, carry out imaginary-time evolution
(using QITE) up to time β (ii) measure the expectation value
of th e observable that one wants to produce a thermal aver-
age for (iii) measure a product operator such as
ˆ
Z
1
ˆ
Z
2
. . .
ˆ
Z
N
,
to collapse back onto a random product state (iv) rep e at (i).
Note th at in step (iii) on e can measure in any product basis,
and randomizing the product basis can be used to reduce the
autocorrelation time and avoid ergodicity problems in sam-
pling.
In Fig.
4 we show the results of quantum METTS (using
exact classical emulation) for the thermal average h
ˆ
Hi as a
function of temperature β, f or the 6-site Heisenberg model
for several temperatu res and domain sizes; sufficiently large
D converges to the exact thermal average at each β; err or bars
reflect only finite QM ETTS samples. We also show an im-
plementation of quantum METTS o n the Aspen-1 QPU and
QVM with a 1-qubit field model (Fig.
4b), and using the QVM
for a 2-qubit AFM transverse field Ising m odel (Fig. 4c).
Conclusions. We have introduced quantum an alogs of
![FIG. 1: (color online) (a) Schematic of the QITE algorithm. Top: imaginary-time evolution under a geometric k-local operator ĥ[m] can be reproduced by a unitary operation acting on D > k qubits. Bottom: exact imaginary-time evolution starting from a product state requires unitaries acting on a domainD that grows with correlations. (b,c) Left: mutual information I(i, j) between qubits i, j as a function of distance d(i, j) and imaginary time β, for a 1D (b) and a 2D (c) FM transverse-field Ising model, with h = 1.25 (1D) and h = 3.5 (2D). I(i, j) saturates at longer times. Right: relative error in the energy∆E and fidelity F = |〈Φ(β)|Ψ〉|2 between the finite-time state Φ(β) and infinite-time state Ψ as a function of β. The noise in the 2D fidelity error at large β arises from the approximate nature of the algorithm used. See SI for details.](/figures/fig-1-color-online-a-schematic-of-the-qite-algorithm-top-mao60ks0.webp)
