Quantum inverse iteration algorithm for programmable quantum simulators

TL;DR: In this article, a quantum inverse power iteration algorithm is proposed to estimate the ground state properties of a programmable quantum device, where the sequential application of the Hamiltonian inverse to an initial state prepares an approximate groundstate.

Abstract: We propose a quantum inverse iteration algorithm which can be used to estimate the ground state properties of a programmable quantum device. The method relies on the inverse power iteration technique, where the sequential application of the Hamiltonian inverse to an initial state prepares an approximate groundstate. To apply the inverse Hamiltonian operation, we write it as a sum of unitary evolution operators using the Fourier approximation approach. This allows to reformulate the protocol as separate measurements for the overlap of initial and propagated wavefunction. The algorithm thus crucially depends on the ability to run Hamiltonian dynamics with an available quantum device. We benchmark the performance using paradigmatic examples of quantum chemistry, corresponding to molecular hydrogen and beryllium hydride. Finally, we show its use for studying the ground state properties of relevant material science models which can be simulated with existing devices, considering an example of the Bose-Hubbard atomic simulator.

Summary

INTRODUCTION

• Quantum computing offers drastic speed up for certain computational problems, and has evolved as a unique direction in the theoretical information science.
• The field of experimental quantum computing is yet at its infancy.
• The algorithms of ever-increasing complexity were implemented on different platforms, with circuit depth exceeding a thousand gates, 2, 3 ultimately allowing for quantum supremacy demonstration.
• 19, 20, 24 VQE can use a chemically inspired ansatz, 24 Hamiltonian variational ansatz, 25 or rely on the variational imaginary time evolution.
• Finally, when applied to the Bose-Hubbard quantum simulator, it allows to study its ground state properties, showing promise as a protocol for analog quantum simulators and near-term quantum devices.

RESULTS

• The authors start by considering a generic interacting system, which can be described by the second quantized Hamiltonian.
• In the fermionic case, Hamiltonian (1) can describe the full configuration interaction problems in quantum chemistry, with operator âj corresponding to molecular orbital j.
• Using the existing mappings between fermionic and spin-1/2 systems, one can rewrite Eq. (1) in the form of a local Hamiltonian Ĥ for interacting qubits, which involves strings of Pauli operators.

Inverse iteration

• The authors propose the procedure which can be seen as a quantum version of the inverse power iteration algorithm for finding the dominant eigenvalue of the matrix, represented by the inverse of Hamiltonian matrix ĤÀ1 , which is treated as a dimensionless matrix in this section.
• Given that Ĥ is invertible, the order of eigenvalues is reversed and, with the appropriate shift of the diagonal to make eigenvalues positive, the power iteration allows to find GSE.

Sequential energy estimation

• The authors final goal is to estimate system observables, provided that the approximate ground state is prepared.
• This allows to examine the protocol for a system of higher complexity (N = 4), comparable to lithium hydrate four-qubit simulation considered in ref.
• As the propagation phase grows, the approximation λ ðaÞ k starts to resemble the idealized iteration procedure.
• The performance of the quantum inverse iteration procedure is further analyzed in Fig. 2c, d where energy distance to ground state and trace distance are shown as a function ϕ max for several fixed iteration steps (k = 2, 4, 7).
• For this the authors have performed the analysis including relevant dephasing processes, which influence the estimate for overlaps (see details in the Supplemental Material).

DISCUSSION

• The authors have presented the algorithm for the ground state energy (GSE) estimation of a quantum Hamiltonian.
• Targeting near-term quantum simulators, the authors described the protocol as a separate estimation of GSE contributions from the wavefunction overlap measurements.
• Both digital and noisy operation was considered, and found to be sufficient for a GSE calculation with chemical accuracy.
• 46 For instance, the analog-type fermionic quantum chemistry simulator 47 would be much valued for the task.
• This poses the question of connection between the measurementbased cooling scheme and the dynamic protocol described in the current study.

Scaling and fault-tolerant implementation

• In this section the authors consider the scaling for the quantum inverse iteration algorithm.
• Moreover, since the authors also require simulation of Hamiltonian dynamics for expðÀiϕ k;'.
• The latest represents conceptually the closest algorithm to the one described in the paper, and thus will serve as benchmark.
• 56 This enlarges the actual circuit depth (while being polynomial), and precludes the implementation of U ¼ expðÀiϕ ĤÞ in analog fashion.

Overlap measurement

• This nicely fits the task of GSE estimation for the fermionic Hamiltonian, as its Hilbert space includes a vacuum state with no fermions present (unless space reduction procedure was performed).
• The HF state can be prepared from the reference jψ R i (vacuum or other product state) using the product of local operators, and the authors note that these states are orthogonal.
• The authors note that in principle these are two related ways to estimate the real part of the overlap, corresponding to Eqs. ( 11) and ( 13).
• For the larger system the measurement is generalized to GHZ-type, and consequently requires increasing number of two qubit operators, which depends on how much an initial HF state is different from the reference state.
• Finally, the authors note that direct measurement is possible in atomic setups, where many-body interferometry is applied to two copies.

Molecular hydrogen Hamiltonian

• The fermionic Hamiltonian ĤH2 is first written in the form of Eq. (1) (main text), where coefficients v ij and V ijkl are calculated by conventional quantum chemistry methods.
• Here, the authors exploited the OpenFermion package for Python, 58 which allows to extract the interfermionic interactions for four Gaussian orbitals fit via STO-3G method and perform the fermions-to-qubits transformation.
• For the small N = 4 system the authors have chosen to use the Jordan-Wigner transformation, although other options may be used as the system size increases.
• The energy scale J for the actual H 2 Hamiltonian corresponds to Hartree units, while for the quantum simulator J corresponds to the effective qubit coupling.
• This shall be done within the chemical precision ϵ, which is equal to ϵ = 0.0016 Hartree, and thus defines the relevant cutoff for the iteration procedure.

Figures

Fig. 4 Bose-Hubbard ground state energy estimation for a chain of five atoms. a Sketch of the system, which depicts coherent tunneling processes with the rate J and the nonlinear interaction of strength U. Correlations are measured for a set of distances r= [0, 1, 2]. b Energy difference for the true GSE and quantum inverse iteration estimate are shown for different iteration steps k (log scale). Several tunneling values are chosen, J∕U= 0.01, 0.05, 0.1, 0.2, capturing the transition between insulating and superfluid behaviour. c–e Correlation functions of the form hâycþr âci are calculated for the approximate ground state at increasing iteration step k, and for several values of J∕U (log scale). Correlations for true ground state are shown in blue (large dots and top curves).

Fig. 3 Beryllium hydride (H–Be–H) ground state energy estimation (N= 8 qubits). a Energy difference for true GSE and quantum inverse iteration shown for different iteration steps k and maximal propagation phases ϕmax (log scale). The red line corresponds to an ideal iteration procedure. b Same energy difference shown as a function of maximal phase at k= 2, 4, 7. In both panels the blue shaded area corresponds to the chemically precise estimate.

Full text

ORE Open Research Exeter
TITLE
Quantum inverse iteration algorithm for programmable quantum simulators
AUTHORS
Kyriienko, O
JOURNAL
npj Quantum Information
DEPOSITED IN ORE
23 January 2020
This version available at
http://hdl.handle.net/10871/40549
COPYRIGHT AND REUSE
Open Research Exeter makes this work available in accordance with publisher policies.
A NOTE ON VERSIONS
The version presented here may differ from the published version. If citing, you are advised to consult the published version for pagination, volume/issue and date of
publication

ARTICLE
OPEN
Quantum inverse iteration algorithm for programmable
quantum simulators
Oleksandr Kyriienko
1,2,3
*
We propose a quantum inverse iteration algorithm, which can be used to estimate ground state properties of a programmable
quantum device. The method relies on the inverse power iteration technique, where the sequential application of the Hamiltonian
inverse to an initial state prepares the approximate ground state. To apply the inverse Hamiltonian operation, we write it as a sum
of unitary evolution operators using the Fourier approximation approach. This allows to reformulate the protocol as separate
measurements for the overlap of initial and propagated wavefunction. The algorithm thus crucially depends on the ability to run
Hamiltonian dynamics with an available quantum device, and can be used for analog quantum simulators. We benchmark the
performance using paradigmatic examples of quantum chemistry, corresponding to molecular hydrogen and beryllium hydride.
Finally, we show its use for studying the ground state properties of relevant material science models, which can be simulated with
existing devices, considering an example of the Bose-Hubbard atomic simulator.
npj Quantum Information (2020) 6:7 ; https://doi.org/10.1038/s41534-019-0239-7
INTRODUCTION
Quantum computing offers drastic speed up for certain
computational problem s, and has evolved as a unique direction
in the theoretical information sc ience.
1
However, the field of
experimental quantum computing is yet at its infanc y. The typical
size of quantum chips for the reliable gate based quantum
computation ranges from one to several tens of physical qubits,
with the main limits posed by decoherence. Despite the
imperfections, the algorithms of ever-increasing complexity were
implemented on different platforms, with circuit depth exceedin g
a thousand gates,
2,3
ultimately allowing for quantum supremacy
demonstration.
4
At the same time, the vox populi of quantum engineers says that
while experimental setups are developed and mastered rapidly,
the theorists in the field lag behind. Whereas by now textbook
examples of quantum algorithms with exponential and quadratic
speed up for factoring and search serve as a great motivation,
1
the
estimates of gate counts are daunting, making them distant goals
for the future fault-tolerant quantum computers.
5
Recent devel-
opments in this fast evolving field call for new short depth
algorithms which can solve useful problems in the era of noisy
intermediate scale quantum (NISQ) devices,
6
and in future lead to
quantum advantage.
One of the most promising directions for quantum computation
is the field of quantum chemistry and materials.
5,7
Targeting the
access to ground state properties of molecules and strongly
correlated matter, it can offer huge gain for various technological
applications, for instance helping to find a catalyst for the nitrogen
fixation.
8
To date, different quantum theoretical protocols were
developed, and several proof-of-principle experiments on various
platforms were performed in the simplest cases. Examples include
simulation of molecular hydrogen with the linear optical setup,
9
superconducting circuits,
10–13
and trapped ions.
14
Finally, the
variational simulation for larger molecules (LiH and BeH
2
) were
reported recently.
11
From the material science perspective, the use
of cold atom quantum simulators has shown great promise, where
simulations of Fermi-Hubbard lattice dynamics,
15
large scale
quantum Rydberg chain
16
and Ising model,
17
and two-
dimensional many-body localization
18
have been performed.
However, in the latter cases the analog approach to simulation
is taken, given an access to unitary dynamics, while precluding the
study of ground state properties.
To access the ground state properties of quantum chemical
Hamiltonian, several routines can be used (see refs
19,20
for the
review). First option corresponds to the quantum phase estima-
tion algorithm (PEA),
21
which exploits unitary dynamics of the
system controlled by register qubits. Although this algorithm is
efficient, giving logarithmic error and polynomial gate scaling, its
implementation requires substantial circuit depth for currently
available circuits.
10
Moreover, the controlled type of operations
require the digitization of the circuit, thus complicating the use of
analog quantum simulators for PEA. Another approach is adiabatic
quantum computing, which was already applied to quantum
chemistry problems.
22
However, the required adiabaticity of
dynamics typically results in the effectively long circuit depth.
Finally, an alternative route to quantum chemistry and materials is
offered by hybrid-classical variational approaches, which were
proposed recently.
23
They rely on term-by-term energy measure-
ment for the prepared trial quantum state (ansatz) with
consequent classical optimization, and are referred to as Varia-
tional Quantum Eigensolvers (VQE).
19,20,24
VQE can use a
chemically inspired ansatz,
24
Hamiltonian variational ansatz,
25
or
rely on the variational imaginary time evolution.
26
The search for a
simple and efficient ansatz represents an important ongoing
research direction.
27,28
In this case the depth of the quantum
circuit is greatly reduced, though at the expense of increased
number of measurements, being favorable strategy for NISQ
devices. For VQE the number of variational parameters scales as
O[(3N)
k
], where N is a number of qubits and k represents an
approximation order.
24
While for quantum chemistry applications
k = 2 suffices to give useful results, these approaches are yet to be
tested for larger system sizes, where the multi-variable
1
Department of Physics and Astronomy, University of Exeter, Stocker Road, Exeter EX4 4QL, UK.
2
ITMO University, Kronverkskiy prospekt 49, Saint Petersburg 197101, Russia.
3
NORDITA, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, SE-106 91 Stockholm, Sweden. *email: kyriienko@ukr.net
www.nature.com/npjqi
Published in partnership with The University of New South Wales
1234567890():,;

optimization may raise problems for the genuine ground state
estimation.
29
In the following, we propose the quantum inverse iteration
algorithm for the estimation of the ground state energy (GSE) of a
quantum system. It is inspired by the classical inverse power
iteration algorithm for finding the dominant eigenstate of the
matrix, where the computationally demanding part of matrix
inversion and multiplication is performed by a quantum circuit.
Previously, a direct iteration approach was considered as a general
purpose quantum algorithm,
30
aiming for large scale fault-
tolera nt implementation. Here, we present the protocol of the
hybrid quantum-classical nature. It relies on performing quantum
evolution for different propagation times and classical post-
processing of the measured observables. The approach is applied
to quantum chemistry examples (H
2
and BeH
2
molecules),
showing favourable scaling with system parameters. Finally, when
applied to the Bose-Hubbard quantum simulator, it allows to study
its ground state properties, showing promise as a protocol for
analog quantum simulators and near-term quantum devices.
RESULTS
We start by considering a generic interacting system, which can
be described by the second quantized Hamiltonian. It can be
written as sum of two-body and four-body parts
^
H¼
X
ij
v
ij
^
a
y
i
^
a
j
þ
X
ijkl
V
ijkl
^
a
y
i
^
a
y
j
^
a
k
^
a
l
;
(1)
where
^
a
y
i
(
^
a
i
) can correspond to fermionic or bosonic creation
(annihilation) operators, and cover broad range of models. In the
fermionic case, Hamiltonian (1) can describe the full configuration
interaction problems in quantum chemistry, with operator
^
a
j
corresponding to molecular orbital j. Using the existing mappings
between fermionic and spin-1/2 systems, one can rewrite Eq. (1)in
the form of a local Hamiltonian
^
H for interacting qubits, which
involves strings of Pauli operators. The task is then to find the
lowest eigenvalue of large matrix
^
H, corresponding to GSE.
Inverse iteration
We propose the procedure which can be seen as a quantum
version of the inverse power iteration algorithm for finding the
dominant eigenvalue of the matrix, represented by the inverse of
Hamiltonian matrix
^
H
1
, which is treated as a dimensionless
matrix in this section. Given that
^
H is invertible, the order of
eigenvalues is reversed and, with the appropriate shift of the
diagonal to make eigenvalues positive, the power iteration allows
to find GSE. Namely, starting with an initial state jψ
0
i, which has
nonzero overlap with the sought ground state jψ
gs
i, by repetitive
application of the inverted matrix one can prepare (unnormalized)
state,
e
ψ
k



¼ð
^
H
1
Þ
k
ψ
0
ji
, such that ψ
gs
jψ
k





2
< ϵ for sufficiently
large number of iterations k ≥ K,
31
where ψ
k
ji¼
e
ψ
k



= k
e
ψ
k



k
(Fig. 1). We note that generally this method has favorable
logarithmic complexity in the iteration depth, being K = log[ϵ
sin
−2
(θ
0
)∕(λ
1
− λ
n
)]∕[2 log(λ
2
∕λ
1
)], where λ
1
, λ
2
, and λ
n
correspond
to dominant, sub-dominant, and smallest eigenvalue of
^
H
1
. Here
sin
2
θ
0
parametrizes the overlap between jψ
0
i and jψ
gs
i, marking
that convergence of the procedure depends on the initial guess,
and generally can be made nonzero taking jψ
0
i as a random state.
While classical power iteration methods generally have good
convergence in the number of iterations, the main caveat comes
from the complexity scaling with the system size N. The
requirement for K matrix multiplications leads to O[K2
2N
]
operations (for dense matrices), yielding exponential scaling. Even
worse situation is for the inverse matrix algorithm, where an
overhead comes from the
^
H
1
calculation, requiring extra O[2
N
]
operations.
Fourier approximation
In the following we show that we can exploit the iterative
procedure with logarithmic iteration depth in ϵ, while providing
exponential speed up for the inverse Hamiltonian multiplication
process. The latter comes from the approximation theory,
32
observing that the inverse can be represented as an integral
x
1
¼
R
þ1
0
expðxyÞdy, which by applying the trapezoidal rule
can be written as a sparse sum of exponents. For quantum
systems a similar idea was proposed in ref.,
33
where Fourier
approximation of the Hamiltonian inverse was presented as a
double integral of the unitary propagator. This was further used to
design an efficient solver for the quantum linear equation system
problem.
33,34
Here we extend the Fourier approximation to the k-
th power of the inverse (k ≥ 1), which formally reads
^
H
k
¼
i N
k
ffiffiffiffiffiffi
2π
p
Z
þ1
0
dy
Z
þ1
1
dzðzy
k 1
Þexpðz
2
=2Þexpðiyz
^
HÞ;
(2)
and N
k
is a normalization factor. The integral can be then
discretized as
^
H
k

iN
k
ffiffiffiffiffiffi
2π
p
X
M
y
1
j
y
¼0
Δ
y
ðj
y
Δ
y
Þ
k 1
X
M
z
j
z
¼M
z
Δ
z
ðj
z
Δ
z
Þexp½j
2
z
Δ
2
z
=2exp½iðj
y
Δ
y
Þðj
z
Δ
z
Þ
^
H;
(3)
where Δ
y,z
correspond to the discretization steps for integration
variables, and M
y,z
represent cutoffs for integration. Notably, once
applied to the physical Hamiltonian inverse, the discretization
variable Δ
z
remains dimensionless, while Δ
y
has the units of
inverse energy, serving akin to discrete time variable. The success
of approximation (3) depends on the condition number of the
Hermitian matrix
^
H, given by the ratio of its largest to smallest
eigenvalue, κ ¼ λ
^
H
n
=λ
^
H
1
. Finally, Eq. (3) can be conveniently
redefined as
^
H
k
¼
X
L
k
‘¼1
c
k ;‘
expði ϕ
k ;‘
^
HÞ 
^
H
k
a
;
(4)
where we have rewritten the double summation in Eq. (3) using
the superindex ℓ(j
y
, j
z
), ϕ
k,ℓ
= ( j
y
Δ
y
)( j
z
Δ
z
) is a phase of evolution for
parameters chosen to discretize k-th inverse, and L
k
= M
y
(2M
z
+ 1).
Here c
k,ℓ
represent purely imaginary coefficients for the series, and
∣c
k,ℓ
∣ define corresponding weights. The required size and number
of discretization steps Δ
y,z
and M
y,z
depends on ϵ and κ (see
Fig. 1 Flowchart of the quantum inverse iteration algorithm. First, the initial product state is prepared and the inverse Hamiltonian
operation is represented as a sum unitary evolution operators. Next, the iterated wavefunction can be formally obtained by applying the
inverse. Finally, physical quantities (e.g., energy) are estimated as expectation values of corresponding operator, and recast as a sum of
wavefunction overlaps.
O. Kyriienko
2
npj Quantum Information (2020) 7 Published in partnership with The University of New South Wales
1234567890():,;

Methods, section A, for the details). Importantly, they set the
maximal evolution phase ϕ
max
= (M
y
Δ
y
)(M
z
Δ
z
), which serves as an
equivalent of the total gate count for analog quantum simulation.
As we need to prepare the approximate ground state by
applying generally nonunitary operator
^
H
K
to the initial state
jψ
0
i, we shall either introduce an ancillary register to perform it, or
properly account for the normalization of the resulting wavefunc-
tion. The former option is an excellent strategy for the future fault-
tolerant devices, and has beneficial scaling (Methods, Section A). It
has deep connection to linear composition of unitaries (LCU)
methods and duality quantum computing.
35
The latter is more
suitable for programmable quantum simulators. In the following
we present the strategy, which can be applied to estimate the
ground state properties by sequential evaluation of terms in the
series. Similarly to VQE approaches, this relies on performing large
number of measurements, and thus adds an extra complexity as
compared to the generic implementation of the inverse operator.
At the same time, term-by-term readout offers better resilience to
errors where even imperfect procedure can yield reasonable GSE
estimate for quantum simulators.
Sequential energy estimation
Our final goal is to estimate system observables, provided that the
approximate ground state is prepared. For any operator
^
A it can
be retrieved from the measurement A ¼hψ
gs
j
^
Ajψ
gs
i=hψ
gs
jψ
gs
i,
where the normalization is accounted for explicitly. In particular,
we are interested in calculating the ground state energy λ
gs
≈ λ
k
,
choosing the operator
^
A as
^
H. This amounts to measurement of
Hamiltonian expectation value for ψ
k
ji¼
^
H
k
ψ
0
jiin the form
λ
k
¼
hψ
k
j
^
Hjψ
k
i
hψ
k
jψ
k
i
:
(5)
We proceed by considering each propagated wavefunction
separately, such that λ
k
can be related to wavefunction overlaps
(see flowchart in Fig. 1). This is motivated by the Hamiltonian
averaging procedure
36
used in VQE to reduce the circuit depth at
the expense of larger number of sequential measurements. Using
Fourier expansion of the inverse Hamiltonian (4), the estimated
energy reads
λ
ðaÞ
k
¼
P
‘;‘
0
hψ
k ;‘
0
j
^
Hjψ
k ;‘
i
P
‘;‘
0
hψ
k ;‘
0
jψ
k ;‘
i
¼
P
‘;‘
0
c

k ;‘
0
c
k ;‘
hψ
0
je
iðϕ
k;‘
ϕ
k;‘
0
Þ
^
H
^
Hjψ
0
i
P
‘;‘
0
c

k ;‘
0
c
k ;‘
hψ
0
je
iðϕ
k;‘
ϕ
k;‘
0
Þ
^
H
jψ
0
i
:
(6)
Note that expression (6) now includes overlaps between initial and
evolved wavefunction for the fixed phase, which shall be
calculated separately for the numerator (“energy”) and denomi-
nator (“norm”). Finally, we note that the wavefunction overlap can
be inferred using different approaches. One option corresponds to
using the SWAP test,
37,38
which represents a common quantum
measurement strategy and requires system doubling. It can be
conveniently realized in some near-term setups, being successfully
demonstrated for cold atom lattices by many-body interferometry
of two copies of a quantum state.
39
The overlap measurement
schemes continue to improve.
40
Additionally, in the Methods,
section B, we describe an alternative approach, which does not
require extra qubits and relies on the measurements of
observables, once the reference state for the system is chosen.
In the previous section we described the general algorithm and
discussed its key properties, namely the scaling and sequential
operation. To show its use for the ground state estimation and
characterize the required resources for realistic problems, we
apply it to quantum chemistry.
Applications: molecular hydrogen
We start with by now the standard example of molecular
hydrogen, H
2
. As a test task we consider the spinful case. This
allows to examine the protocol for a system of higher complexity
(N = 4), comparable to lithium hydrate four-qubit simulation
considered in ref.
11
The details of mapping of quantum chemical
structure into qubits are presented in Methods, section C. In the
following we work with four-qubit molecular hydrogen Hamilto-
nian
^
H
H
2
, with all eigenenergies shifted to positive values. Starting
from the Hartree-Fock (HF) energy λ
0
, the task is to estimate GSE
λ
gs
, using the protocol described in the preceding section. This
shall be done within the chemical precision ϵ, which is equal to
ϵ = 0.0016 Hartree, and thus defines the relevant cutoff for the
iteration procedure.
We start by benchmarking the inverse power procedure in its
general form, and define how many iteration steps one needs to
come close to the ground state. For this, we first perform the
inverse Hamiltonian iteration in the ideal setting, assuming that an
exact inverse is known. Then, we compare it to the quantum
inverse iteration, which uses the Fourier approximation (4). GSE is
estimated using the measurement of propagated and initial
wavefunction overlaps. To quantify the performance two char-
acteristics are employed. The first, and the most natural one,
corresponds to the difference between estimated energy value λ
k
and true GSE λ
gs
, being Δλ∕J ≡ (λ
k
− λ
gs
)∕J. It allows to observe the
convergence and provides an indication of how well the
procedure works for a given system. The second quantity
corresponds to the trace distance between an idealized inverse
iteration matrix
^
H
k
and its approximation
^
H
k
a
,defined as a half
of trace norm for the difference of two matrices. It reveals the
actual success of mimicking the ideal inverse in full generality. At
the same time, this is the quantity which cannot be straightfor-
wardly observed in the experiment, and only serves for the
analysis.
The results of the inverse power iteration for molecular
hydrogen Hamiltonian
^
H
H
2
are shown in Fig. 2a as a function of
iteration step k. The ideal version of inverse iteration is plotted in
red and reveals exponential convergence to GSE. The chemical
precision is achieved already at the second iteration step, as
depicted by the blue shaded area starting at Δλ = 1.6 × 10
−3
J. The
idealized case is then compared to the quantum inverse iteration
procedure with combined measurement of wavefunction overlaps
as stated in Eq. (6). Here, we assumed that the genuine unitary
evolution with Hamiltonian
^
H
H
2
is run in the analog simulation
fashion. The case of digital evolution with associated Trotterization
technique and its benchmarking is considered in the Supple-
mental Material, where we also present the circuit scheme for
digital evolution. The approximation was performed using equal
number of steps M
z
= M
y
= 30, and the discretization values Δ
z
=
Δ
y
J were adjusted to match the maximal propagation phases of
ϕ
max
∕2π = J(M
y
Δ
y
)(M
z
Δ
z
)∕2π = {0.3, 0.35, 0.6, 0.95, 1.35} (here, the
propagation phase is taken to be dimensionless by absorbing
energy unit prefactor J from the Hamiltonian). The corresponding
curves show the improvement of the quantum power iteration
estimation for increasing number of iteration steps. The conver-
gence rate also depends on the maximal phase of the
propagation. For small phases (top curves in Fig. 2a), the initial
estimator does not give successful convergence, but comes closer
to GSE for large k. As the propagation phase grows, the
approximation λ
ðaÞ
k
starts to resemble the idealized iteration
procedure. However, this only happens up to a certain value of k
past which the approximate energy grows, thus deviating from
the ideal solution. From the point of view of process fidelity, the
trace distance Tr½
^
H
k
;
^
H
k
a
 between ideal and approximate
O. Kyriienko
3
Published in partnership with The University of New South Wales npj Quantum Information (2020) 7

inverse operators increases monotonically with k (Fig. 2b). The
increase of ϕ
max
allows to reduce Tr½
^
H
k
;
^
H
k
a
 at each k.
The performance of the quantum inverse iteration procedure is
further analyzed in Fig. 2c, d where energy distance to ground
state and trace distance are shown as a function ϕ
max
for several
fixed iteration steps (k = 2, 4, 7). Calculations were performed
accounting for two different ways of arranging the phase. First, the
approximation grid was fixed setting M
y
= M
z
= 30 while changing
Δ
z
= Δ
y
J (solid curves in Fig. 2c, d). In the second case the fixed
step size Δ
z
= Δ
y
J = 0.05 was combined with the increment of M
z,y
(dashed curves in Fig. 2c, d). For both energy distance (Fig. 2c) and
trace distance (Fig. 2d) we observe no difference between two
approximation procedures, but clear indication of the importance
of maximal propagation phase (time). For Δλ one sees a non-
monotonic dependence on ϕ
max
, which starts with a decrease of
the energy difference for increasing maximal phase (ϕ
max
∕2π <
0.4). At larger phases the dependence experiences pronounced
dips (note the log scale), which are more visible for many
iterations. Overall the difference remains well-within chemical
precision and experiences saturation. When the trace distance is
considered, one sees that success of the approximation mono-
tonically improves with ϕ
max
. At the same time, for fixed
approximation parameters {M
z,y
, Δ
z,y
} it is more difficult to
represent inverse iteration operator faithfully, in-line with scaling
analysis discussed in Methods. Finally, the comparison of results in
Fig. 2c, d allows to suggest that nonmonotonicity in the
spectroscopic signatures can come from the particular structure
of the Hamiltonian and the initial state, where certain phases
might be preferable (i.e., not all elements of the Hamiltonian
matrix contribute equally to the inverse iteration procedure).
To decide on the optimal way to approximate the inverse, we
consider different discretization steps for y and z auxiliary
variables, characterized by the skewness parameter Δ
y
J∕Δ
z
. The
calculation is done for M
z
= M
y
= 30 with the maximal phase fixed
to ϕ
max
∕2π = 0.92. The results are shown in Fig. 2e, f as a function
of skew. The energy difference parameter shows that for
approximating the inverse for small iteration numbers (k = 2
curve in Fig. 2e) larger skew factors are preferable, with z variable
requiring finer approximation. However, for increased iteration
number the optimum flows to Δ
y
J∕Δ
z
~ 1 values, suggesting close-
to-equal spacing can work well for varied k. Examining the trace
distance, we see that in unbiased setting the skew ratio of Δ
y
J∕Δ
z
~
2 is preferable.
Finally, the very important issue to address is an influence of
noise on the operation of quantum inverse iteration protocol. For
this we have performed the analysis including relevant dephasing
processes, which influence the estimate for overlaps (see details in
the Supplemental Material). Although noise makes the estimation
of energies at large iteration step k less reliable, it is possible to
estimate the energy within chemical precision using simple noise
mitigation techniques.
Applications: beryllium hydride
To test the scalability of the approach, we consider a molecule of
bigger size, which requires larger Hilbert space simulation. For this,
we choose to simulate beryllium hydride (BeH
2
) in the full spinful
version using N = 8 qubits (see Methods, section D, for the details).
We proceed in the same manner as for H
2
molecule, and
quantify the operation of the quantum inverse iteration procedure
for BeH
2
. The approximation parameters were chosen as Δ
y
J = Δ
z
=
0.05, with the number of discretization points M
y,z
adjusted
accordingly to maintain maximal propagation phase. The results
of the simulation are shown in Fig. 3. The first plot (Fig. 3a) shows
that ideal iteration works well for the beryllium hydride, with
chemically precise GSE obtained already at k = 1 iteration step.
The Fourier approximation for the inverse at small phases does
not reach required accuracy, while for the increased iteration step
number and ϕ
max
∕2π > 1 chemically accurate ground state
Fig. 2 Molecular hydrogen (H–H) example for benchmarking the quantum inverse power iteration (N = 4 qubits). a Energy difference
between the exact ground state and quantum inverse iteration estimate, shown as a function of the iteration step k for different maximal
phases of Fourier approximation (log scale). Solid red line shows the result for the ideal inverse iteration. Blue shaded area corresponds to the
chemically precise estimate (same in c, e). b Trace distance Tr½
^
H
k
;
^
H
k
a
 between the ideal and approximate inverse operators shown as a
function of iteration step number for different phases. c Energy difference Δλ vs maximal propagation phase at different k (log scale). d Trace
distance Tr½
^
H
k
;
^
H
k
a
 vs phase for k = 2, 4, 7. e, f Energy difference (e) and trace distance (f) plotted for the fixed maximal phase of ϕ
max
∕2π =
0.92, but different arrangement of the approximation grid defined by the skew parameter Δ
y
J∕Δ
z
. Several iteration steps k = 2, 4, 7 are
depicted.
O. Kyriienko
4
npj Quantum Information (2020) 7 Published in partnership with The University of New South Wales

Frequently asked questions

Q1. What are the contributions mentioned in the paper "Quantum inverse iteration algorithm for programmable quantum simulators" ?

In this paper, the authors proposed a hybrid-classical variational approach for the estimation of ground state properties of molecules and strongly correlated matter. 

Q2. What is the main caveat of the inverse power iteration algorithm?

While classical power iteration methods generally have good convergence in the number of iterations, the main caveat comes from the complexity scaling with the system size N. The requirement for K matrix multiplications leads to O[K22N] operations (for dense matrices), yielding exponential scaling. 

Q3. What is the energy scale for the actual H2 Hamiltonian?

The energy scale J for the actual H2 Hamiltonian corresponds to Hartree units, while for the quantum simulator J corresponds to the effective qubit coupling. 

Q4. How many discretization points were chosen to maintain the maximum propagation phase?

The approximation parameters were chosen as ΔyJ= Δz= 0.05, with the number of discretization points My,z adjusted accordingly to maintain maximal propagation phase. 

Q5. What was the original goal of the inverse iteration approach?

a direct iteration approach was considered as a general purpose quantum algorithm,30 aiming for large scale faulttolerant implementation. 

Q6. What is the performance of the quantum inverse iteration procedure?

The performance of the quantum inverse iteration procedure is further analyzed in Fig. 2c, d where energy distance to ground state and trace distance are shown as a function ϕmax for several fixed iteration steps (k = 2, 4, 7). 

Q7. How many bits of precision is required to approach an error of = 2m?

The complexity of IPEA was discussed in ref., 55 showing the requirement of O [log(ϵ)log(log(ϵ)∕ϵ)] phase iterations to approach an error of ϵ= 2−m (energy is rescaled such that k Ĥ k < 2π, and m is the number of relevant bits of precision, typically limited to <20 for quantum chemistry applications). 

Q8. How is the model solvable in the so called Mott insulating regime?

In particular, the model was shown to be easily solvable in the socalled Mott insulating regime where U ≫ J and ground state corresponds to the product state of one atom per site, jψMotti = ∏i j1ii. 

Q9. What is the way to implement the sum of unitary operators?

One possible option here is the amplitude amplification approach,51 which addresses the task of implementing the sum of unitary operators, of the same type as the one in Eq. (4).