# State preparation and evolution in quantum computing: a perspective from
Hamiltonian moments

Joseph C. Aulicino  Pritzker Molecular Engineering, University of Chicago,
5640 S Ellis Ave, Chicago, IL 60637, United States of America  Trevor Keen
Department of Physics and Astronomy, University of Tennessee, Knoxville,
Tennessee 377996, United States of America  Bo Peng<EMAIL_ADDRESS>Physical
Sciences and Computational Division, Pacific Northwest National Laboratory,
Richland, WA 99354, United States of America

###### Abstract

Quantum algorithms on noisy intermediate-scale quantum (NISQ) devices are
expected to soon have the ability to simulate quantum systems that are
classically intractable, demonstrating a quantum advantage. However, the non-
negligible gate error present on NISQ devices impedes the implementation of
conventional quantum algorithms. Practical strategies usually exploit hybrid
quantum-classical algorithms to demonstrate potentially useful applications of
quantum computing in the NISQ era. Among the numerous hybrid quantum-classical
algorithms, recent efforts highlight the development of quantum algorithms
based upon quantum computed Hamiltonian moments,
$\langle\phi|\hat{\mathcal{H}}^{n}|\phi\rangle$ ($n=1,2,\cdots$), with respect
to quantum state $|\phi\rangle$. In this tutorial review, we will give a brief
review of these quantum algorithms with focuses on the typical ways of
computing Hamiltonian moments using quantum hardware and improving the
accuracy of the estimated state energies based on the quantum computed
moments. Furthermore, we will present an example to show how we can measure
and compute the Hamiltonian moments of a four-site Heisenberg model, and
compute the energy and magnetization of the model utilizing the imaginary time
evolution on the real IBM-Q NISQ hardware. Along this line, we will further
discuss some practical issues associated with these algorithms. We will
conclude this tutorial review by discussing some possible developments and
applications in this direction in the near future.

quantum computing, state preparation and evolution, Hamiltonian moments,
hybrid quantum classical algorithm, quantum simulation

###### Contents

  1. I Introduction
  2. II Quantum computation of Hamiltonian moments
    1. II.1 Term-by-term measurement
    2. II.2 Expansion by Chebyshev polynomials via quantum walk
    3. II.3 Linear combination of unitary time propagator
  3. III Classical computation of target states from quantum computed moments
    1. III.1 Lanczos approach
    2. III.2 Real time evolution
    3. III.3 Imaginary time evolution
    4. III.4 Variational simulation of real and imaginary time dynamics
  4. IV Quantum simulation of a model Hamiltonian
    1. IV.1 Computation preparation
    2. IV.2 Quantum measurement
    3. IV.3 Estimated ground state energy and magnetism
    4. IV.4 Discussions
      1. IV.4.1 In comparison with VQE
      2. IV.4.2 Grouping and readout error mitigation
      3. IV.4.3 Approximation in ITE
  5. V Conclusion and outlook
  6. VI Acknowledgement
  7. VII Funding Information
  8. VIII Data Availability Statement

##  I Introduction

The original idea of solving many-body quantum mechanical problems on a
quantum computer dates back to R. Feynman’s insight almost four decades ago
Feynman (1982), when high-performance classical computing was at its early
state and quantum computing was only a thought experiment. In the last forty
years, developments in classical quantum approaches and high performance
computing technology yielded adequate tools to perform highly accurate small-
to medium-size quantum simulations and qualitative large-scale quantum
simulations. Despite these exciting developments and achievements, the
classical quantum approach is gradually approaching the computing limit
attributed to the unfavorable scaling along with the expanded system size. For
instance, consider the well-known examples of full Configuration Interaction
(full CI) approach David Sherrill and Schaefer (1999) and Coupled Cluster (CC)
approach Bartlett and Musiał (2007). Classically, in order to achieve the
chemical accuracy ($\sim$1 kcal/mol) for the energy evalulation, the runtime
in the full CI and the “golden standard” CCSD(T) scale as
$\mathcal{O}\Big{(}\begin{subarray}{c}N\\\ n\end{subarray}\Big{)}$ and
$\mathcal{O}(N^{7})$, respectively, with respect to the system size $N$ (and
number of electrons $n$), which impedes their routine applications for large
systems. On the other hand, quantum computing (simulation performed on a
quantum computer) is a next-generation computing technology that bears the
hope of outperforming high-performance classical computing techniques to solve
a wide class of problems encountered in scientific research and industrial
applications. Quantum computation enables the encoding of an exponential
amount of information about many-body quantum systems into a polynomial number
of qubits, which provides a more efficient way to expand and explore the
computational state-space in polynomial time Sugisaki  _et al._ (2019).

Since the first application of quantum computing algorithm for computational
chemistry problems in 2005 Aspuru-Guzik  _et al._ (2005), numerous quantum
algorithms have been proposed aiming at demonstrating the potential quantum
speedup, in comparison to the classical approaches, for solving problems in
domain science such as chemistry, nuclear physics, quantum field theory, and
high energy physics O’Malley  _et al._ (2016); Linke  _et al._ (2017); Kandala
_et al._ (2017); Dumitrescu  _et al._ (2018); Klco  _et al._ (2018); Colless
_et al._ (2018); McCaskey  _et al._ (2019). However, many proposed quantum
algorithms featuring favorable scaling typically have deep circuit demands,
and are not suitable for the near-term noisy intermediate-scale quantum (NISQ)
devices Preskill (2018). In the search for quantum algorithms suitable running
on the NISQ devices, one should mention hybrid quantum-classical Variational
Quantum algorithm (VQA) Peruzzo  _et al._ (2014); McClean  _et al._ (2016);
Romero  _et al._ (2018); Shen  _et al._ (2017); Kandala  _et al._ (2017,
2019); Colless  _et al._ (2018); Huggins  _et al._ (2020), quantum approximate
optimization algorithm (QAOA) Farhi, Goldstone, and Gutmann (2014), quantum
annealing Bharti  _et al._ (2021); Albash and Lidar (2018), gaussian boson
sampling Aaronson and Arkhipov (2011), analog quantum simulation Trabesinger
(2012); Georgescu, Ashhab, and Nori (2014), iterative quantum assisted
eigensolver McArdle  _et al._ (2019); Motta  _et al._ (2020); Parrish and
McMahon (2019); Kyriienko (2020), and many others, which bear the hope to
capitalize the near-term NISQ quantum devices. Nevertheless, these NISQ
quantum algorithms have their own limitations. For example, VQAs face
challenges such as Barren plateau McClean  _et al._ (2018); Cerezo  _et al._
(2021); Cerezo and Coles (2021), ansatz searching Herasymenko and O’Brien
(2019), and measurement overhead Bonet-Monroig, Babbush, and O’Brien (2020)
which stimulate the emerging of some many recent improvements. Numerous
comprehensive reviews on the NISQ quantum algorithms and techniques have
appeared in recent years to pave the foundation of and provide guidance to
future studies (see e.g. Refs. 35; 21 for more recent ones). However, due to
the rapid development of this emerging field, it becomes challenging for the
reviews to comprehensively cover many newly developed quantum algorithms and
techniques in a timely manner. Furthermore, it is usually more efficient, in
particular for non-specialists, to choose a thread to understand and digest
the basic ideas and the interconnection of these algorithms by going through
some tutorials. With these concerns in mind, instead of giving another
comprehensive review covering some newly developed quantum algorithms, here we
present a tutorial review on some recently developed quantum algorithms based
on the quantum computation of Hamiltonian moments. The Hamiltonian moment here
is referring to the expectation value of Hamiltonian powers with respect to a
given state $|\phi\rangle$ for the systems of interest
$\langle\phi|\hat{\mathcal{H}}^{n}|\phi\rangle$ ($n=1,2,\cdots$), and can be
considered as one of the building blocks for performing variational or
perturbative calculations solving the energy and property of a many-body
quantum system. We will introduce the state preparation and evolution, and
hybrid algorithms associated with quantum Hamiltonian moments. With the
examples and discussions, we hope to provide a relatively easy-to-read review
with a clear demonstration of the connections and differences between some
relevant quantum algorithms.

This tutorial review is organized as follows. Section II reviews the ways of
quantum computing Hamiltonian moments. Section III reviews some recently
proposed hybrid approaches based on quantum Hamiltonian moments (see Fig. 1
for a typical workflow of the hybrid algorithm of this kind). Section IV gives
a tutorial of how to target the ground state and magnetization of a
generalized four-site Heisenberg model employing quantum Hamilonian moments in
the imaginary time evolution approach on IBM-Q quantum hardware. Section IV
also includes the discussions of some practical issues associated with
performing quantum computing on the NISQ device with a focus on the impact of
these issues on the results of quantum simulations. We will conclude this
tutorial review by pointing out possible developments and applications in this
direction in the near future.

Figure 1: The workflow of a typical hybrid quantum classical algorithm based
on quantum computed Hamiltonian moments $\langle\hat{\mathcal{H}}^{r}\rangle$
($r=1,2,\cdots$).

##  II Quantum computation of Hamiltonian moments

Given a system Hamiltonian $\hat{\mathcal{H}}$ of interest and a trial state
$|\phi\rangle$, the $n$-th order ($n\geq 1$) Hamiltonian moment with respect
to $|\phi\rangle$ is defined as

$\displaystyle\langle\hat{\mathcal{H}}^{n}\rangle\equiv\langle\phi|\hat{\mathcal{H}}^{n}|\phi\rangle.$
(1)

In classical computations, $\langle\hat{\mathcal{H}}^{n}\rangle$’s are
obtained by power iteration, i.e. a repeated multiplication of
$\hat{\mathcal{H}}$ and $|\phi\rangle$. However, as the dimension of the
Hilbert space grows exponentially with the system size, and the number of
terms that constitutes $\langle\hat{\mathcal{H}}^{n}\rangle$ would quickly
blow up as the power $n$ becomes large, the classical computations of
$\langle\hat{\mathcal{H}}^{n}\rangle$’s would quickly become intractable. On
the other hand, quantum computing allows for a more straightforward encoding
of the quantum states defined in a Hilbert space of potentially large and
classically intractable dimensions. Therefore, we consider (a) how to utilize
quantum resources for the direct computation of
$\langle\hat{\mathcal{H}}^{n}\rangle$, and (b) how to exploit the noisy
$\langle\hat{\mathcal{H}}^{n}\rangle$’s computed on NISQ devices to accurately
evaluate ground/excited states and dynamics of quantum many-body systems. In
this section, we try to answer the first question by briefly reviewing some of
the methods that have been proposed for directly computing
$\langle\hat{\mathcal{H}}^{n}\rangle$ on quantum devices. We leave the
discussion of the second question to Section III.

To directly compute $\langle\hat{\mathcal{H}}^{n}\rangle$ on quantum devices
for a quantum many-body system, a first step is to encode the many-body
Hamiltonian in terms of qubits. There have been many ways proposed to encode
the system Hamiltonian. For example, in the Jordan-Wigner Jordan and Wigner
(1928) and Bravyi-Kitaev methods Bravyi and Kitaev (2002); Seeley, Richard,
and Love (2012); Tranter  _et al._ (2015), the $N_{q}$-qubit system
Hamiltonian is encoded as a linear combination of tensor products of Pauli
qubit operators,

$\displaystyle\hat{\mathcal{H}}=\sum_{m=1}^{M}h_{m}P_{m},$ (2)

where $h_{m}$ is a scalar and
$P_{m}\equiv{\prod_{q=1}^{N_{q}}}^{\otimes}\sigma_{i}(q)$ is a Pauli string
with $\sigma_{i}(q)$ being either $2\times 2$ identity matrix or one of the
three Pauli matrices

$\displaystyle\sigma_{x}=\left(\begin{array}[]{cc}0&1\\\
1&0\end{array}\right),~{}~{}\sigma_{y}=\left(\begin{array}[]{cc}0&-i\\\
i&0\end{array}\right),~{}~{}\sigma_{z}=\left(\begin{array}[]{cc}1&0\\\
0&-1\end{array}\right)$ (9)

for the $q$-th qubit. From (2), a naive way to directly compute
$\langle\hat{\mathcal{H}}^{n}\rangle$ is to plug (2) into (1), and evaluate
the expectation values of the products of Pauli strings that constitute
$\langle\hat{\mathcal{H}}^{n}\rangle$. Apparently, without any simplification
or reduction, the naive way would require evaluating $\mathcal{O}(M^{n})$
terms and would quickly become intractable as $M$ and/or $n$ scale up.
Therefore, major efforts in this direction focus on mitigating this evaluation
overhead as much as possible by introducing some simplifications and
reductions based on the properties of the Pauli strings.

###  II.1 Term-by-term measurement

One straightforward way is to apply Pauli reduction and commutativity to
reduce the number of terms, and then do term-by-term measurements. The Pauli
reduction simply utilizes the commutation and anti-commutation relations
between Pauli matrices

$\displaystyle\left\\{\begin{array}[]{cl}[\sigma_{i},\sigma_{j}]&=2i\epsilon_{ijk}\sigma_{k}\\\
\\{\sigma_{i},\sigma_{j}\\}&=2\delta_{ij}I\end{array}\right.\Rightarrow\sigma_{i}\sigma_{j}=i\epsilon_{ijk}\sigma_{k}+\delta_{ij}I,$
(12)

where $\sigma_{i},\sigma_{j},\sigma_{k}$ are Pauli matrices, $\epsilon_{ijk}$
is the structure constant following the Levi-Civita symbol Tyldesley (1975),
and $\delta_{ij}$ is the Kronecker delta. From (12) it is straightforward to
see the product of Pauli strings is another Pauli string with an appropriate
phase factor. Therefore, in comparison to the evaluation of
$\langle\hat{\mathcal{H}}\rangle$, the evaluation of
$\langle\hat{\mathcal{H}}^{n}\rangle$ does not increase the circuit depth.
Furthermore, we notice that the measurement of a single Pauli string could be
used to determine classes of contributions to Hamiltonian moments of arbitrary
order. For example, since every Pauli string is unitary, the evaluation of
$\langle P_{i}\rangle$ can simultaneously contribute all the $\langle
H^{2k+1}\rangle$’s ($k\geq 0,k\in\mathbb{N}$) Kowalski and Peng (2020).
Suchsland et al. systematically analyzed the actual number of Pauli strings
corresponding to $\hat{\mathcal{H}}^{n}$ ($n=1,2,3$) for a collection of 11
molecules, including Hn ($n=2,4,6,8,10$), LiH, NaH, H2O, NH3, and N2, that
require up to $N_{q}=18$ qubits Suchsland  _et al._ (2021). It was found that
the actual numbers of Pauli strings in these systems after applying (12)
exhibit at most $\mathcal{O}(N_{q}^{3})$, $\mathcal{O}(N_{q}^{5})$, and
$\mathcal{O}(N_{q}^{7})$ scalings in comparison to the predicted
$\mathcal{O}(N_{q}^{4})$, $\mathcal{O}(N_{q}^{8})$, and
$\mathcal{O}(N_{q}^{12})$ scalings for $\hat{\mathcal{H}}^{n}$ ($n=1,2,3$)
respectively.

Figure 2: A pictorial QWC grouping procedure.

On the top of the above simplification, a more significant term reduction can
be achieved by utilizing the commutativity of the Pauli strings. To see this
method works mathematically, consider operators $\hat{\mathcal{A}}$ and
$\hat{\mathcal{B}}$ that commute, i.e.
$[\hat{\mathcal{A}},\hat{\mathcal{B}}]=0$, then there exists a complete
orthonormal eigenbasis ${|u_{i}\rangle}$ ($i=0,1,\cdots$) that simultaneously
diagonalizes both operators, and we can write

$\displaystyle\hat{\mathcal{A}}=\sum_{i}\lambda_{i}|u_{i}\rangle\langle
u_{i}|,~{}~{}\hat{\mathcal{B}}=\sum_{i}\kappa_{i}|u_{i}\rangle\langle u_{i}|,$
(13)

with $\\{\lambda_{i},i=0,1,\cdots\\}$ and $\\{\kappa_{i},i=0,1,\cdots\\}$
being the eigenvalues of $\hat{\mathcal{A}}$ and $\hat{\mathcal{B}}$,
respectively. Now the expectation values of $\hat{\mathcal{A}}$ and
$\hat{\mathcal{B}}$ with respect to a trial state $|\phi\rangle$ can be
expressed as

$\displaystyle\langle\hat{\mathcal{A}}\rangle$
$\displaystyle=\langle\phi|\hat{\mathcal{A}}|\phi\rangle=\sum_{i}\lambda_{i}|\langle
u_{i}|\phi\rangle|^{2},$ $\displaystyle\langle\hat{\mathcal{B}}\rangle$
$\displaystyle=\langle\phi|\hat{\mathcal{B}}|\phi\rangle=\sum_{i}\kappa_{i}|\langle
u_{i}|\phi\rangle|^{2}.$ (14)

Therefore, if we measure each $|\langle u_{i}|\phi\rangle|^{2}$’s, we can
deduce $\langle\hat{\mathcal{A}}\rangle$ and $\langle\hat{\mathcal{B}}\rangle$
based on (14).

The simplest commutativity among Pauli strings is qubitwise commutativity
(QWC) Kandala  _et al._ (2017); McClean  _et al._ (2016), where two Pauli
strings qubitwise commute if two single-qubit Pauli matrices at each index
commute, i.e.

$\displaystyle[I,I]=[I,\sigma_{x}]=[I,\sigma_{y}]=[I,\sigma_{z}]=0,$
$\displaystyle[\sigma_{x},\sigma_{x}]=[\sigma_{y},\sigma_{y}]=[\sigma_{z},\sigma_{z}]=0.$
(15)

For example, $\\{\sigma_{x}\otimes\sigma_{y},\sigma_{x}\otimes
I,I\otimes\sigma_{y},I\otimes I\\}$ is a QWC class, since any pair of Pauli
strings in this class commute (a more pictorial QWC grouping procedure is
shown in Fig. 2). For the observables grouped in one QWC class, one can easily
find a shared eigenbasis that simultaneously diagonalizes all the observables.
Take the previous example, if we want to measure the QWC class in the
computational basis, utilizing the Clifford representation of single Pauli
matrices $\sigma_{x}$ and $\sigma_{y}$

$\displaystyle\sigma_{x}=H\sigma_{z}H,$
$\displaystyle\sigma_{y}=S\sigma_{x}S^{\dagger}=SH\sigma_{z}HS^{\dagger},$
(16)

we will only need perform the following measurement

which gives the probabilities of the variational circuit rotated into the
shared eigenbasis, and the expectation values of each Pauli string can then be
recovered by classically combining these probabilities with the eigenvalues of
these Pauli strings.

Now, given a large amount of Pauli strings that constitute
$\hat{\mathcal{H}}^{n}$, how should we find an optimal grouping of these Pauli
strings such that the total number of measurements is minimized?
Unfortunately, it turns out that finding the optimal grouping is NP-hard Karp
(2010); Gokhale  _et al._ (2020), and can be converted to minimum clique cover
or minimum graph coloring problem, which in practice can be approximately
solved through heuristic approaches that scale quadratically with the number
of the Pauli strings Yen, Verteletskyi, and Izmaylov (2020). In terms of term
reduction, for small molecules the QWC grouping can bring roughly 75% saving
for evaluating $\langle\hat{\mathcal{H}}\rangle$, and even bigger savings for
for $\langle\hat{\mathcal{H}}^{n}\rangle$ ($n>1$). Early studies of applying
QWC bases to evaluating $\langle\hat{\mathcal{H}}^{4}\rangle$ for the
Heisenberg models represented by 10 to 40 qubits exhibits a sub-linear scaling
of the number of measurements in the number of the qubits Vallury  _et al._
(2020). Furthermore, for some small molecular systems, it has been shown that
the number of QWC bases grouping the Pauli strings that constitute
$\hat{\mathcal{H}}^{n}$ eventually reach a plateau, regardless of the power
$n$ Claudino  _et al._ (2021).

Essentially, the QWC grouping has its roots in mutually unbiased bases (MUB)
Schwinger (1960); Klappenecker and Rötteler (2004) from quantum information
theory associated with maximizing the information learned from a single
measurement. In general, MUBs group the potentially $4^{N_{q}}-1$
$N_{q}$-qubit Pauli strings (excluding identity string) into commuting groups
of maximal size, of which using QWCs would allow reducing this requirement to
$3^{N_{q}}$ Altepeter, James, and Kwiat (2004). Further, it is known that
there in principle exists an MUB grouping for $N_{q}$-qubit that can further
reduce this requirement to even $2^{N_{q}}+1$ with maximum $2^{N_{q}}-1$
distinct Pauli strings in each commuting group, however, the entanglement in
the MUBs Lawrence, Brukner, and Zeilinger (2002) makes MUB quantum state
tomography challenging. In practice, beside QWC, other grouping techniques
have alternatively been proposed, including general commutativity Gokhale  _et
al._ (2020); Yen, Verteletskyi, and Izmaylov (2020), unitary partitioning
Izmaylov  _et al._ (2020), and Fermionic basis rotation grouping Huggins  _et
al._ (2021). Generally speaking, the applications of these grouping rules for
evaluating $\langle\hat{\mathcal{H}}\rangle$ have shown that, at a cost of
introducing additional one-/multi-qubit unitary transformation before the
measurement, the total number of terms can be significantly reduced from
$\mathcal{O}(N_{q}^{4})$ to $\mathcal{O}(N_{q}^{2\sim 3})$. In particular, for
simpler cases where the entire Hamiltonian or its power could even be
transformed by a single unitary, the evaluation of the corresponding
expectation values can be done in a single set of measurements.

It is worth mentioning that the above discussion is limited to the number of
terms that can be efficiently reduced through the groupings that feature the
commutativity of Pauli strings, while the total number of measurements
required from the number of groups also critically depends on the covariance
between Pauli strings, $\text{cov}\big{(}P_{i},P_{j}\big{)}$, and the desired
precision $\epsilon$. Given that we can always write the Hamiltonian powers as
linear combination of Pauli strings similar to Eq. (2), the total number of
measurements for evaluating the moments can then be expressed as Gonthier  _et
al._ (2020); Rubin, Babbush, and J. (2018); Wecker, Hastings, and Troyer
(2015)

$\displaystyle\text{\\# of Measurements}=$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\Bigg{(}\displaystyle\frac{\sum_{G}\sqrt{\sum_{i,j,\in
G}h_{i}h_{j}\text{cov}\big{(}P_{i},P_{j}\big{)}}}{\epsilon}\Bigg{)}^{2}$ (17)

with $G$ indexing the groups. Therefore, it is likely that the number of
groups decreases at the cost of introducing larger covariances that could
essentially increase the total number of measurements required to achieve a
desired precision.

There are also many other advanced measurement schemes proposed recently. One
example is to simultaneously obtain expectation values of multiple observables
by randomly measuring and projecting the quantum state into classical shadows
Huang, Kueng, and Preskill (2020, 2021); Aaronson (2018); Struchalin  _et al._
(2021); Chen  _et al._ (2021); Zhao, Rubin, and Miyake (2021); Acharya, Saha,
and Sengupta (2021); Hadfield (2021); Hillmich  _et al._ (2021); Zhang  _et
al._ (2021). The algorithm in principle enables the measurements of $M$ low-
weight observables using only $\mathcal{O}(\log_{2}M)$ samples. The practical
performance of the algorithm for model and molecular Hamiltonians on NISQ
device, in terms of accuracy and efficiency, is still under intense study.

###  II.2 Expansion by Chebyshev polynomials via quantum walk

Another way to exactly implement $\hat{\mathcal{H}}^{n}$ is by Chebyshev
polynomial expansion. We can exactly express a monomial $x^{n}$ on
$\left[-1,1\right]$ as a finite sum of Chebyshev polynomials of the first
kind, i.e.

$\displaystyle x^{n}=\sum\limits_{k=0}^{n}C_{nk}T_{k}(x)$ (18)

where

$\displaystyle
C_{nk}=\begin{cases}\dfrac{1}{2^{n-1}}~{}\binom{n}{(n-k)/2}&\text{if
}(n-k)\text{ is even, }j>0\\\\[10.00002pt]
\dfrac{1}{2^{n}}~{}\binom{n}{n/2}&\text{if }n\text{ is even,
}j=0\\\\[10.00002pt] 0&\text{else}\end{cases}$ (19)

and $T_{k}(x)$ is the $k^{th}$ degree Chebyshev polynomial of the first kind
Subramanian, Brierley, and Jozsa (2019). This expansion is useful because a
“quantum walk” can exactly produce the effect of Chebyshev polynomials in
$\hat{\mathcal{H}}/d$, where $d$ is the sparsity of the matrix. A quantum walk
is a process in which the Hilbert space is enlarged, and then a walk operator
that works in the enlarged Hilbert space is repeatedly implemented. For
details on quantum walks and their implementations, see e.g. Refs. 68; 69.
With amplitude amplification to obtain a constant success probability of
preparing a quantum state and a flag qubit to indicate if the preparation was
successful, the query complexity of this method is
$\mathcal{O}\left(nd^{n-1}\lambda_{0}^{-n}\right)$ Subramanian, Brierley, and
Jozsa (2019). Using this method without amplitude amplification, the query
complexity of this implementation is $\mathcal{O}\left(n\right)$. Here,
$\lambda_{0}$ denotes the smallest (magnitude) eigenvalue of $H$, with the
assumption that $\hat{\mathcal{H}}$ does not have 0 as an eigenvalue. However,
this method has its own drawbacks. The implementation of the quantum walk can
be expensive, requiring a number of ancillary qubits that grows linearly with
the system size Subramanian, Brierley, and Jozsa (2019). The Chebyshev
polynomial expansion’s potential application for computing the imaginary time
evolution of an unnormalized quantum state, and in general
$e^{-\beta\hat{\mathcal{H}}}$ for evaluating partition function has been
discussed in Ref. 70, where they use the recursion relation for Chebyshev
polynomials.

###  II.3 Linear combination of unitary time propagator

Recently, Seki and Yunoki proposed a quantum power method to evaluate
$\hat{\mathcal{H}}^{n}|\phi\rangle$ Seki and Yunoki (2021), where the
Hamiltonian power $\hat{\mathcal{H}}^{n}$ is first expressed as the $n$-th
time derivative of $\hat{\mathcal{U}}(t)=e^{i\hat{\mathcal{H}}t}$ at $t=0$,
then the time derivative is approximated through the central finite difference
(CFD) which can be alternatively expressed as a linear combination of
$\hat{\mathcal{U}}(\Delta_{t})$’s at different time variables $\Delta_{t}$’s
in close proximity to $t=0$, i.e.

$\displaystyle\hat{\mathcal{H}}^{n}$
$\displaystyle=i^{n}\frac{\partial^{n}\hat{\mathcal{U}}(t)}{\partial
t^{n}}\bigg{|}_{t=0}\approx
i^{n}\left[\hat{\mathcal{U}}(\frac{\Delta_{t}}{2})-\hat{\mathcal{U}}(-\frac{\Delta_{t}}{2})\right]^{n}/\Delta_{t}^{n}$
$\displaystyle=\sum_{k=0}^{n}\left[\frac{i^{n}}{\Delta_{t}^{n}}(-1)^{k}\Big{(}\begin{subarray}{c}n\\\
k\end{subarray}\Big{)}\right]\left[\hat{\mathcal{U}}(\frac{\Delta_{t}}{2})\right]^{n-2k}.$
(20)

To implement the approach to quantum circuit, each
$\hat{\mathcal{U}}(\frac{\Delta_{t}}{2})$ is further decomposed using
symmetric Suzuki-Trotter decomposition (SSTD) Trotter (1959); Suzuki (1990).
There are two sources of error in this approach, the CFD error and SSTD error.
The systematic CFD error scales quadratically with $\Delta_{t}$, i.e.
$\mathcal{O}(\Delta_{t}^{2})$, while the $2m$-th order SSTD brings the
systematic error of $\mathcal{O}(\Delta_{t}^{2m})$. However, there are several
advantages in the proposed approach. First, both aforementioned errors can be
systematically suppressed through Richardson extrapolation Richardson and
Gaunt (1927); Temme, Bravyi, and Gambetta (2017) at the cost of involving more
terms in the linear combination in Eq. (20), which on the other hand offers
the space for using low order SSTD. Second, the number of gates required for
approximate $\hat{\mathcal{H}}^{n}$ for an $N_{q}$-qubit Hamiltonian
$\hat{\mathcal{H}}$ scales as $\mathcal{O}(nN_{q})$. In principle, the quantum
power method can be combined with many classical approaches, which will be
elaborated in Section III for targeting the ground/excited states and
properties of a quantum many-body system. As shown in Ref. 71, numerical
noiseless simulations employing the quantum power method in Krylov-subspace
diagonalization for targeting ground state of model systems demonstrate
systematically improved accuracy over the conventional VQE.

Inspired by the classical inverse power iteration approach for finding the
dominant eigenstate of a given hermitian matrix with a more favorable
logarithmic complexity in the iteration depth Sachdeva and Vishnoi (2014),
Kyriienko recently proposed a quantum inverse iteration algorithm Kyriienko
(2020) for approximately computing the inverse of $n$-th power of the
Hamiltonian, $\hat{\mathcal{H}}^{-n}$. The approach extends the discretized
Fourier approximation of the Hamiltonian inverse Childs, Kothari, and Somma
(2017) to the $n$-th power,

$\displaystyle\hat{\mathcal{H}}^{-n}\approx\frac{iN_{n}}{\sqrt{2\pi}}\sum_{j=0}^{J}y_{j}^{n-1}\Delta_{y}\sum_{k=-K}^{K}z_{k}\Delta_{z}e^{-z_{k}^{2}/2}\hat{\mathcal{U}}_{j,k},$
(21)

where $N_{n}$ is a normalization factor, $y_{j}\equiv j\Delta_{y}$ and
$z_{k}\equiv k\Delta_{z}$ with the integers $j\in[0,J]$ and $k\in[-K,K]$ and
intervals $\Delta_{y,z}\in\mathbb{R}$ defining a discretization grid space,
and unitary $\hat{\mathcal{U}}_{j,k}$ is given by
$e^{-iy_{j}z_{k}\hat{\mathcal{H}}}$. Then
$\langle\phi|\hat{\mathcal{H}}^{-n}|\phi\rangle$ is evaluated through the SWAP
test Ekert  _et al._ (2002); Higgott, Wang, and Brierley (2019), overlap
measurement Mitarai and Fujii (2019), or Bell-type-like measurement Kyriienko
(2020). However, this method requires a number of calls to a time-evolution
oracle that is dependent on the condition number of the system, as well as a
requirement that the Hamiltonian has been shifted such that all of the
eigenvalues are positive. This latter requirement requires the user to have
some knowledge of the ground state energy of the system a priori. A suitable
guess can be found, for example, using the minimum label finding algorithm of
Ref. 81.

Alternatively, if considering a fault-tolerant implementation with favorable
resource scaling $\hat{\mathcal{H}}^{n}|\phi\rangle$ or
$\hat{\mathcal{H}}^{-n}|\phi\rangle$ can be done through amplitude
amplification approach Brassard  _et al._ (2002), Hamiltonian simulation Berry
_et al._ (2015), qubitization Low and Chuang (2017), or the direct block-
encoding Gilyén  _et al._ (2019) methods at the cost of introducing deeper
circuit and implementing controlled-$\hat{\mathcal{U}}$ operations which,
however, are still challenging to be implemented in the NISQ devices. Some of
these methods have deep connection to linear combination of unitaries (LCU)
technique Childs and Wiebe (2012). The basic idea of LCU can be briefly
demonstrated from the following toy example. Suppose we want to implement
$\hat{\mathcal{U}}_{0}+\hat{\mathcal{U}}_{1}$ on $|\phi\rangle$. We can start
by introducing an ancilla qubit and preparing $|+\rangle|\phi\rangle$ with
$|+\rangle\equiv\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$. Then by performing
the controlled unitary $|0\rangle\langle
0|\otimes\hat{\mathcal{U}}_{0}+|1\rangle\langle 1|\hat{\mathcal{U}}_{1}$ on
$|+\rangle|\phi\rangle$ we are able to obtain the state
$\frac{1}{\sqrt{2}}(|0\rangle\hat{\mathcal{U}}_{0}|\phi\rangle+|1\rangle\hat{\mathcal{U}}_{1}|\phi\rangle)$.
Now, if we measure the ancilla qubit in the $\\{|+\rangle,|-\rangle\\}$ basis
and obtain the outcome associated with $|+\rangle$, then we obtain a state
proportional to $(\hat{\mathcal{U}}_{0}+\hat{\mathcal{U}}_{1})|\phi\rangle$.

##  III Classical computation of target states from quantum computed moments

Now supposing we already have necessary Hamiltonian moments computed from
quantum simulation, how can we proceed to accurately evaluate the energy
and/or properties of the Hamiltonian? In this section, we will try to give our
answers from some classical approaches.

###  III.1 Lanczos approach

One of the most common approaches that can be exploited is the Lanczos
approach Lanczos (1950). In the Lanczos approach, an orthonormal set of
quantum states $\\{|\phi_{i}\rangle,i=1,2,\cdots\\}$ is generated through a
three-term recursion

$\displaystyle\hat{\mathcal{H}}|\phi_{i+1}\rangle=\beta_{i}|\phi_{i}\rangle+\alpha_{i+1}|\phi_{i+1}\rangle+\beta_{i+1}|\phi_{i+2}\rangle,$
(22)

where $\alpha_{i}=\langle\phi_{i}|\hat{\mathcal{H}}|\phi_{i}\rangle$,
$\beta_{i}=\langle\phi_{i}|\hat{\mathcal{H}}|\phi_{i+1}\rangle$. The three-
term recursion essentially transforms the Hamiltonian $\hat{\mathcal{H}}$ into
a tridiagonal form

$\displaystyle\hat{\mathcal{T}}_{i}=\left[\begin{array}[]{cccc}\alpha_{1}&\beta_{1}&&\\\
\beta_{1}&\alpha_{2}&\ddots&\\\ &\ddots&\ddots&\beta_{i-1}\\\
&&\beta_{i-1}&\alpha_{i}\end{array}\right]$ (27)

where the matrix elements $\alpha_{i}$ and $\beta_{i}$ can essentially be
represented recursively in terms of Hamiltonian moments. Explicitly, it can be
shown

$\displaystyle\alpha_{1}$ $\displaystyle=\langle\hat{\mathcal{H}}\rangle,$
$\displaystyle\beta_{1}$
$\displaystyle=\sqrt{\langle\hat{\mathcal{H}}^{2}\rangle-\langle\hat{\mathcal{H}}\rangle^{2}},$
$\displaystyle\alpha_{2}$
$\displaystyle=\frac{\langle\hat{\mathcal{H}}^{3}\rangle-2\langle\hat{\mathcal{H}}^{2}\rangle\langle\hat{\mathcal{H}}\rangle+\langle\hat{\mathcal{H}}\rangle^{3}}{\langle\hat{\mathcal{H}}^{2}\rangle-\langle\hat{\mathcal{H}}\rangle^{2}},$
$\displaystyle\beta_{2}$
$\displaystyle=\frac{\langle\hat{\mathcal{H}}^{4}\rangle-(\alpha_{1}+\alpha_{2})^{2}\beta_{1}^{2}-(\alpha_{1}^{2}+\beta_{1}^{2})^{2}}{\beta_{1}^{2}},$
(28) $\displaystyle\vdots$

Now the strict upper bound of the ground state energy of the Hamiltonian can
be obtained by directly diagonalizing $\hat{\mathcal{T}}_{i}$. It has been
shown that the classical Lanczos scheme can directly be applied even employing
the noisy $\langle\hat{\mathcal{H}}^{n}\rangle$ Suchsland  _et al._ (2021),
which effectively provides a way to improve the energy estimate from real
quantum hardware without being subject to specific hardware and an explicit
description of the underlying noise. Preliminary tests on H2, H3, and four-
site tetrahedral Heisenberg model demonstrate that combining the Lanczos
scheme with VQE algorithm helps correct the noise and enhance the quality of
the ansatz.

Alternatively, the matrix elements $\alpha_{i}$ and $\beta_{i}$ can be
expressed in terms of the so-called connected moments (or cumulant) Horn and
Weinstein (1984). An $n$-th order connected moments is defined as

$\displaystyle\langle\hat{\mathcal{H}}^{n}\rangle_{c}\equiv\langle\phi|\hat{\mathcal{H}}^{n}|\phi\rangle-\sum_{m=0}^{n-1}\left(\begin{array}[]{c}n-1\\\
m-1\end{array}\right)\langle\hat{\mathcal{H}}^{m}\rangle_{c}\langle\phi|\hat{\mathcal{H}}^{n-m}|\phi\rangle,$
(31)

from which (28) can be generalized as $z$-expansions of $\alpha$ and $\beta$
Hollenberg (1993); Witte and Hollenberg (1994); Hollenberg and Witte (1996)

$\displaystyle\alpha(z)$
$\displaystyle=\langle\hat{\mathcal{H}}\rangle_{c}+z\left(\frac{\langle\hat{\mathcal{H}}^{3}\rangle_{c}}{\langle\hat{\mathcal{H}}^{2}\rangle_{c}}\right)+\mathcal{O}(z^{2}),$
$\displaystyle\beta^{2}(z)$
$\displaystyle=z\langle\hat{\mathcal{H}}^{2}\rangle_{c}+z^{2}\left(\frac{\langle\hat{\mathcal{H}}^{2}\rangle_{c}\langle\hat{\mathcal{H}}^{4}\rangle_{c}-\langle\hat{\mathcal{H}}^{3}\rangle_{c}^{2}}{2\langle\hat{\mathcal{H}}^{2}\rangle_{c}^{2}}\right)+\mathcal{O}(z^{3}),$
(32)

where $z\equiv i/V$ with $i$ the recursion index and $V$ the volume of the
system. From the bound analysis of orthogonal polynomials van Doorn (1987);
Ismail and Li (1992), it has been shown that the true ground state energy of
the Hamiltonian can be approximated through the greatest lower bound (i.e.
infimum) of the expansion (32) Hollenberg and Witte (1996),

$\displaystyle E_{0}\approx E_{0}^{\rm inf}=\inf_{z>0}[\alpha(z)-2\beta(z)].$
(33)

For example, the first order in $z$ gives Hollenberg and Witte (1994)

$\displaystyle E_{0}^{\rm
inf}=\langle\hat{\mathcal{H}}\rangle_{c}-\frac{\langle\hat{\mathcal{H}}^{2}\rangle_{c}^{2}}{\langle\hat{\mathcal{H}}^{3}\rangle_{c}^{2}-\langle\hat{\mathcal{H}}^{2}\rangle_{c}\langle\hat{\mathcal{H}}^{4}\rangle_{c}}\times$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\Bigg{(}\sqrt{3\langle\hat{\mathcal{H}}^{3}\rangle_{c}^{2}-2\langle\hat{\mathcal{H}}^{2}\rangle_{c}\langle\hat{\mathcal{H}}^{4}\rangle_{c}}-\langle\hat{\mathcal{H}}^{3}\rangle_{c}\Bigg{)}.$
(34)

Approach (34) has been employed to demonstrate the energy estimate of 2D
quantum magnetism model on lattices up to 25 qubits on IBM-Q quantum hardware
Vallury  _et al._ (2020), where the results show a consistent improvement of
the estimated energies in comparison to the VQE results for the same ansatz.
Another application of the Lanczos method is the continued fractions expansion
of the Green’s function. The Green’s function can be calculated in the Lehmann
representation using the $\alpha$ and $\beta$ parameters from the Lanczos
method. For full details, see e.g. Chapter 8 of Ref. 95.

It is worth mentioning that standard Lanczos is a variational approach while
the infimum approach (33) does not offer a strict upper bound of the true
ground state energy, which, as exemplified by (34), is essentially an
alternative size-extensive connected moment expansion of the ground state
energy. Other polynomial expansion approaches, in particular connected moment
expansions, will be discussed in details in Section III.3.

###  III.2 Real time evolution

Alternatively, one can resort to time evolution approaches to target the
energy of the Hamiltonian. Here the time evolution approaches include both
real time evolution (RTE) approach and its imaginary analogue. In RTE, the
dynamics of a quantum state $|\phi_{t}\rangle$ is governed by the time-
dependent Schrödinger equation

$\displaystyle i\frac{\partial|\phi_{t}\rangle}{\partial
t}=\hat{\mathcal{H}}|\phi_{t}\rangle.$ (35)

In the Schrödinger picture the solution of Eq. (35), $|\phi_{t}\rangle$, can
be expressed by acting a time-evolution operator $\hat{U}(t)=e^{-i\hat{H}t}$
(or its truncated Taylor expansion in practice) on an initial state
$|\phi_{0}\rangle$,

$\displaystyle|\phi_{t}\rangle$ $\displaystyle=\hat{U}(t)|\phi_{0}\rangle$
$\displaystyle=e^{-i\hat{\mathcal{H}}t}|\phi_{0}\rangle=\sum_{n=0}^{+\infty}\frac{(-it)^{n}}{n!}\hat{\mathcal{H}}^{n}|\phi_{0}\rangle.$
(36)

The introduction of $\hat{\mathcal{H}}^{n}$ in (36) effectively builds a non-
orthogonal Krylov subspace. Now consider an order-($r+1$) Krylov subspace
$\mathcal{K}_{r+1}=\\{\hat{\mathcal{H}}^{n}|\phi_{0}\rangle,n=0,1,\cdots,r\\}$,
the ground state can then be approximated as a linear combination of the
Krylov basis

$\displaystyle|\phi_{t}\rangle$
$\displaystyle=\sum_{n=0}^{r}c_{n}\hat{\mathcal{H}}^{n}|\phi_{0}\rangle$ (37)

under the constraint $\langle\phi_{t}|\phi_{t}\rangle=1$ with $c_{n}$’s being
the expansion coefficients to be determined. Plugging the expanded form (37)
to Eq. (35) and projecting the approximate evolution on to the subspace
$\mathcal{K}_{r+1}$, we end up with a set of time-dependent coupled equations
in terms of the first $2r+1$ Hamiltonian moments, which in matrix form can be
expressed as

$\displaystyle\dot{\mathbf{c}}=-i\mathbf{L}^{-1}\mathbf{R}\mathbf{c}$ (38)

where the vector $\mathbf{c}=(c_{0},c_{1},\cdots,c_{n})^{T}$ and
$\dot{\mathbf{c}}$ is its time-derivative. Hamiltonian moments matrices
$\mathbf{L}$ and $\mathbf{R}$ are defined as follows,

$\displaystyle\mathbf{L}$
$\displaystyle=\left(\begin{array}[]{cccc}1&\langle\hat{\mathcal{H}}\rangle&\cdots&\langle\hat{\mathcal{H}}^{r}\rangle\\\
\langle\hat{\mathcal{H}}\rangle&\langle\hat{\mathcal{H}}^{2}\rangle&\cdots&\langle\hat{\mathcal{H}}^{r+1}\rangle\\\
\vdots&\vdots&\ddots&\vdots\\\
\langle\hat{\mathcal{H}}^{r}\rangle&\langle\hat{\mathcal{H}}^{r+1}\rangle&\cdots&\langle\hat{\mathcal{H}}^{2r}\rangle\end{array}\right)$
(43) $\displaystyle\mathbf{R}$
$\displaystyle=\left(\begin{array}[]{cccc}\langle\hat{\mathcal{H}}\rangle&\langle\hat{\mathcal{H}}^{2}\rangle&\cdots&\langle\hat{\mathcal{H}}^{r+1}\rangle\\\
\langle\hat{\mathcal{H}}^{2}\rangle&\langle\hat{\mathcal{H}}^{3}\rangle&\cdots&\langle\hat{\mathcal{H}}^{r+2}\rangle\\\
\vdots&\vdots&\ddots&\vdots\\\
\langle\hat{\mathcal{H}}^{r+1}\rangle&\langle\hat{\mathcal{H}}^{r+2}\rangle&\cdots&\langle\hat{\mathcal{H}}^{2r+1}\rangle\end{array}\right).$
(48)

Eq. (38) is a first order, linear ordinary differential equation (ODE), whose
solution depends on the eigenvalues of the matrix $\mathbf{L}^{-1}\mathbf{R}$
that essentially approximate some of the eigenvalues of $\hat{\mathcal{H}}$.
In a recent study Guzman and Lacroix (2021), this approach is demonstrate to
be able to to target both the ground and excited states for pairing
Hamiltonian. Alternatively, as shown in Ref. 71, employing Rayleigh-Ritz
technique, $\mathbf{L}$ and $\mathbf{R}$ constitute general eigenvalue problem

$\displaystyle\mathbf{R}\mathbf{v}=\epsilon\mathbf{L}\mathbf{v}$ (49)

with the eigenpair ($\epsilon,\mathbf{v}$) approximate the energy and state
vector of the true ground state. To solve Eq. (49) for
($\epsilon,\mathbf{v}$), employ canonical orthogonalization to do the
following transformation

$\displaystyle\mathbf{U}^{\dagger}\mathbf{R}\mathbf{U}(\mathbf{U}^{-1}\mathbf{v})=\epsilon\mathbf{U}^{\dagger}\mathbf{R}\mathbf{v}=\epsilon(\mathbf{U}^{-1}\mathbf{v}),$
(50)

where $\mathbf{U}\mathbf{U}^{\dagger}=\mathbf{L}^{-1}$. Since $\mathbf{L}$ is
a hermitian that can be diagonalized by a unitary matrix $\mathbf{V}$ with
eigenvalues constituting a diagonal matrix $\mathbf{s}$, i.e.
$\mathbf{L}=\mathbf{V}\mathbf{s}\mathbf{V}^{\dagger}$, $\mathbf{U}$ can be
constructed through $\mathbf{U}=\mathbf{V}\mathbf{s}^{-1/2}$. Similar subspace
diagonalization schemes have been adopted in many other hybrid quantum-
classical algorithms such as quantum subspace expansion McClean  _et al._
(2017); Colless  _et al._ (2018), and some variants of VQE Parrish  _et al._
(2019); Nakanishi, Mitarai, and Fujii (2019); Huggins  _et al._ (2020).

###  III.3 Imaginary time evolution

The imaginary time evolution (ITE) approach has a long history of being a
robust computational approach to solve the ground state of a many-body quantum
system. To see how it works, let’s first assume the quantum state at time $t$,
$|\phi_{t}\rangle$, is now expanded in the eigen-space of $\hat{\mathcal{H}}$,
$\\{|\psi_{n}\rangle,n=0,1,\cdots\\}$, with $E_{n}$’s being the corresponding
eigenvalues, then as we show in the previous section $|\phi_{t}\rangle$ can be
expressed as

$\displaystyle|\phi_{t}\rangle=\sum_{n}c_{n}e^{-iE_{n}t}|\psi_{n}\rangle$ (51)

with $c_{n}$’s being the expansion coefficients at the initial time. Now, if
we consider a replacement $\tau=it$ ($\tau$ is often called imaginary time)
and confine $\tau\geq 0$, Eq. (51) becomes

$\displaystyle|\phi_{\tau}\rangle=\sum_{n}c_{n}e^{-E_{n}(\tau-\tau_{0})}|\psi_{n}\rangle.$
(52)

Comparing Eqs. (51) and (52), we see that the wavefunction is driven from an
“oscillating” superposition of the Hamiltonian eigenstates to an “exponential
decaying” superposition of the eigenstates with the decay rate proportional to
$E_{n}$. More importantly, in the limit of large $\tau$

$\displaystyle|\phi_{\tau\gg 1}\rangle\approx
c_{0}e^{-E_{0}\tau}|\psi_{0}\rangle,$ (53)

the ground state $\psi_{0}(x)$ is “screened out” because its exponential decay
rate is the smallest in the eigen spectrum of $\hat{\mathcal{H}}$. Briefly
speaking, given a trial wavefunction $|\phi\rangle$ that has non-zero overlap
with the true ground state wavefunction $\psi_{0}$,
$\langle\phi|\psi_{0}\rangle\neq 0$, the ITE approach, in comparison to the
RTE of the trial state, guarantees a monotonically decreasing energy
functional in $\tau$ that converges to the ground state energy at
$\tau\rightarrow\infty$, i.e.

$\displaystyle
E_{0}=\lim_{\tau\rightarrow+\infty}E(\tau)~{}~{}\text{with}~{}~{}E(\tau)=\langle\phi_{\tau}|\hat{\mathcal{H}}|\phi_{\tau}\rangle.$
(54)

Here the ITE of the trial wavefunction, $|\phi_{t}\rangle$, is defined as

$\displaystyle|\phi_{t}\rangle=\frac{e^{-\tau\hat{\mathcal{H}}/2}|\phi\rangle}{\left(\langle\phi|e^{-\tau\hat{\mathcal{H}}}|\phi\rangle\right)^{1/2}}.$
(55)

and the first derivative of $E(\tau)$ with respect to $\tau$ is non-positive

$\displaystyle\frac{\rm{d}E(\tau)}{\rm{d}\tau}=-\big{(}\langle\phi_{t}|\hat{\mathcal{H}}^{2}|\phi_{t}\rangle-\langle\phi_{t}|\hat{\mathcal{H}}|\phi_{t}\rangle^{2}\big{)}\leq
0,$ (56)

where the equality sign holds if and only if $|\phi_{t}\rangle$ is an exact
eigenstate of $\hat{\mathcal{H}}$ (e.g. $|\psi_{0}\rangle$).

The development and application of ITE approach targeting the ground state
wave function and energy dates back to 1970s, when the similar random-walk
imaginary-time technique were developed for diffusion Monte-Carlo methods
Davies  _et al._ (1980); Anderson (1975, 1979, 1980). Later on, in the pursuit
of non-perturbative analytical tool for a wide variety of Hamiltonian systems
that can be systematically improved, various polynomial expansions of Eq. (54)
have been studied extensively. For example, a straightforward way is to
utilize a finite number of Hamiltonian moments, $\hat{\mathcal{H}}^{n}$, to
constitute a truncated Taylor expansion of $E(\tau)$ which we will have a
detailed examination in Section IV.

t-expansion Beyond the conventional Taylor expansion, to reproduce asymptotic
behavior of $E(\tau)$ over a longer imaginary time, Horn and Weinstein
introduced a power series expansion in imaginary time $\tau$ for $E(\tau)$ in
the early 1980s Horn and Weinstein (1984),

$\displaystyle
E(\tau)=\sum_{n=0}^{\infty}\frac{(-\tau)^{n}}{n!}\langle\hat{\mathcal{H}}^{n+1}\rangle_{c},$
(57)

where $\langle\hat{\mathcal{H}}^{n+1}\rangle_{c}$ is the connected moments
defined in Eq. (31). The $\tau\rightarrow\infty$ behavior of $E(\tau)$ was
then reconstructed through Padé approximants to Eq. (57). Furthermore, it was
found that it might be more efficient and error-resilient to first construct
($L,L+M$)-Padé approximants ($M\geq 2$) to $\frac{\rm{d}E(\tau)}{\rm{d}\tau}$
rather than to $E(\tau)$, and then integrate the Padé approximant from 0 to
the critical $\tau$-value (where Eq. (56) becomes positive) to obtain a larger
set of approximations to $E(\tau)$. Early applications of the Padé
approximants of Eq. (57) on Heisenberg and Ising models demonstrated
remarkable improvement upon mean-field results.

Connected moment expansion Based on the Horn-Weinstein theorem, Cioslowski
further derived a more practical re-summation technique, the so-called
connected moment expansion (CMX), to Eq. (57) to address the algebraic
$\tau$-independent form of $E(\tau)$ Cioslowski (1987),

$\displaystyle\lim_{\tau\rightarrow+\infty}E(\tau)=\lim_{n\rightarrow\infty}E_{0}^{\rm
CMX(\it
n)}=\lim_{n\rightarrow\infty}\left(I_{1}-\frac{S_{2,1}^{2}}{S_{3,1}}\left(1+\frac{S_{2,2}^{2}}{S_{2,1}^{2}S_{3,2}}\left(1+\cdots\left(1+\frac{S_{2,n-1}^{2}}{S_{2,n-2}^{2}S_{3,n-1}}\right)\cdots\right)\right)\right),$
(58)

where the $S_{n,i}$’s are defined from the recursion

$\displaystyle S_{n,1}$
$\displaystyle=\langle\hat{\mathcal{H}}^{n}\rangle_{c},~{}~{}n=2,3,\ldots,$
$\displaystyle S_{n,i+1}$
$\displaystyle=S_{n,i}S_{n+2,i}-S_{n+1,i}^{2},~{}~{}i\geq 1.$

It is worth mentioning that further analysis of the CMX recursion had
suggested that Eq. (58) might be recast in a compact matrix form Knowles
(1987); Stubbins (1988)

$\displaystyle\lim_{n\rightarrow\infty}E^{{\rm CMX}(n)}$
$\displaystyle=\langle\hat{\mathcal{H}}\rangle_{c}-\lim_{n\rightarrow\infty}\left(\langle\hat{\mathcal{H}}^{2}\rangle_{c}~{}~{}\cdots~{}~{}\langle\hat{\mathcal{H}}^{n}\rangle_{c}\right)\mathbf{X}$
(59)

where vector $\mathbf{X}$ is the solution of the following linear system

$\displaystyle\left(\begin{array}[]{ccc}\langle\hat{\mathcal{H}}^{3}\rangle_{c}&\cdots&\langle\hat{\mathcal{H}}^{n+1}\rangle_{c}\\\
\vdots&\ddots&\vdots\\\
\langle\hat{\mathcal{H}}^{n+1}\rangle_{c}&\cdots&\langle\hat{\mathcal{H}}^{2n-1}\rangle_{c}\end{array}\right)\mathbf{X}=\left(\begin{array}[]{c}\langle\hat{\mathcal{H}}^{2}\rangle_{c}\\\
\vdots\\\ \langle\hat{\mathcal{H}}^{n}\rangle_{c}\end{array}\right).$ (66)

The analytical properties of CMX and its comparison with other methods (e.g.
Lanczos approach) have been extensively discussed in the literature Knowles
(1987); Prie  _et al._ (1994); Mancini, Zhou, and Meier (1994); Ullah (1995);
Mancini  _et al._ (1995); Fessatidis  _et al._ (2006, 2010). In particular,
similar to perturbational theories, CMX is conceptually simple and size-
extensive, and the accuracy can be easily tuned through the rank of the
connected moments included in the approximation and/or the quality of the
trial wave function. For the latter, there have been discussions on using
multi-configurational Cioslowski (1987) or even correlated wavefunction (e.g.
truncated CI or CC wave functions Noga, Szabados, and Surján (2002)) in the
CMX framework for accelerating the CMX convergence rate. Nevertheless, there
were two major problems associated with the CMX calculations. First, the
algebraic structure of the series expansion in the CMX could cause singularity
Mancini, Zhou, and Meier (1994). Early efforts tried to address this issue
with limited success by employing alternative moments expansion Mancini, Zhou,
and Meier (1994); Ullah (1995); Mancini  _et al._ (1995), generalized moments
expansion Fessatidis  _et al._ (2006), or generalized Padé expansion Knowles
(1987). Second, CMX results are not variational. In comparison with the
variational Lanczos methods Mancini, Prie, and Massano (1991); Prie  _et al._
(1994), it was found that CMX might be considered as a limiting case of the
strict Lanczos scheme, and their exact equivalence only holds in certain
regions of parameter space. A systematic analysis of low order Lanczos and CMX
ground state energy for model Hamiltonians such as harmonic oscilaltor,
anharmonic oscillator, and Kondo model concluded that the accuracy of both
approaches was dependent on the region of parameter space being studied Prie
_et al._ (1994).

Peeters-Devreese-Soldatov approach Regarding the variational expansion, there
had been another interesting yet less known moment approach proposed by
Peeters and Devreese Peeters and Devreese (1984), and further analyzed by
Soldatov Soldatov (1995) in the 1990s, which we will refer to as the PDS
approach. This approach originates in generalizing the Bogolubov inequality
Bogolubov (1947) and the Feynman inequality Feynman (1955) through the
analysis of their Laplace transforms to obtain the upper bounds for the free
energy. In particular, in the operator formalism, given a system Hamiltonian
operator $\hat{\mathcal{H}}$ with ground state $E_{0}$, and a complex scalar
$s\in\mathbb{C}$ with $\Re(s)>-E_{0}$, we can write the following inequality
Peeters and Devreese (1984); Devreese, Evrard, and Kartheuser (1975)

$\displaystyle\langle e^{-\hat{\mathcal{H}}}\rangle\geq
e^{-\langle\hat{\mathcal{H}}\rangle}\underset{\mathcal{L}^{-1}}{\stackrel{{\scriptstyle\mathcal{L}}}{{\Longleftrightarrow}}}\left\langle\frac{1}{s+\hat{\mathcal{H}}}\right\rangle\geq\frac{1}{s+\langle\hat{\mathcal{H}}\rangle}.$
(67)

with $\mathcal{L}$ and $\mathcal{L}^{-1}$ being the Laplace transform and
inverse Laplace transform operators, and $\langle\cdot\rangle$ the expectation
value with respect to a given trial wavefunction $|\phi\rangle$. The PDS
formalism is based on expanding and analyzing
$\left\langle\frac{1}{s+\hat{\mathcal{H}}}\right\rangle$ using a simple
identity (with introducing a parameter $a_{1}\in\mathbb{R}$, and
$a_{1}\neq-s$)

$\displaystyle\frac{1}{s+\hat{\mathcal{H}}}$
$\displaystyle=\frac{1}{s+a_{1}}-\frac{\hat{\mathcal{H}}-a_{1}}{(s+a_{1})^{2}}+\frac{1}{s+\hat{\mathcal{H}}}\cdot\frac{(\hat{\mathcal{H}}-a_{1})^{2}}{(s+a_{1})^{2}}.\;$
(68)

Remarkably, the expectation value of the third term on the right hand side of
identity (68), defined as a residual term $R_{1}(s,a_{1})$, is non-negative,

$\displaystyle
R_{1}(s,a_{1})=\left\langle\frac{1}{s+\hat{\mathcal{H}}}\cdot\frac{(\hat{\mathcal{H}}-a_{1})^{2}}{(s+a_{1})^{2}}\right\rangle\geq
0,~{}~{}(\Re(s)>-E_{0}).$ (69)

By analyzing its first and second derivatives with respect to the induced
parameter $a_{1}$, it is easy to show that when
$a_{1}=\langle\hat{\mathcal{H}}\rangle$, $R_{1}(s,a_{1})$ reaches its local
minima.

Now if we recursively expand $\frac{1}{s+\hat{\mathcal{H}}}$ in (68) one more
time, and introducing another parameter parameter $a_{2}\in\mathbb{R}$, and
$a_{2}\neq-s$, we have

$\displaystyle\frac{1}{s+\hat{\mathcal{H}}}$
$\displaystyle=\frac{1}{s+a_{1}}-\frac{\hat{\mathcal{H}}-a_{1}}{(s+a_{1})^{2}}+\left(\frac{1}{s+a_{2}}-\frac{\hat{\mathcal{H}}-a_{2}}{(s+a_{2})^{2}}\right)\times$
$\displaystyle~{}~{}~{}~{}\frac{(\hat{\mathcal{H}}-a_{1})^{2}}{(s+a_{1})^{2}}+R_{2}(s,a_{1},a_{2})$
(70)

with a refined residual term $R_{2}(s,a_{1},a_{2})$ also being non-negative
for $\Re(s)>-E_{0}$

$\displaystyle R_{2}(s,a_{1},a_{2})$
$\displaystyle=\left\langle\frac{1}{s+\hat{\mathcal{H}}}\cdot\frac{(\hat{\mathcal{H}}-a_{1})^{2}}{(s+a_{1})^{2}}\cdot\frac{(\hat{\mathcal{H}}-a_{2})^{2}}{(s+a_{2})^{2}}\right\rangle\geq
0.$ (71)

Examining the first and second derivatives of $R_{2}(s,a_{1},a_{2})$ with
respect to $a_{1}$, $a_{2}$ shows when

$\displaystyle\left\\{\begin{array}[]{l}\langle(\hat{\mathcal{H}}-a_{1})\cdot(\hat{\mathcal{H}}-a_{1})(\hat{\mathcal{H}}-a_{2})\rangle=0\\\
\langle(\hat{\mathcal{H}}-a_{2})\cdot(\hat{\mathcal{H}}-a_{1})(\hat{\mathcal{H}}-a_{2})\rangle=0\end{array}\right.,$
(74)

$R_{2}$ reaches its local minima.

Note that there is an important feature about $a_{i}$ ($i=1,2$), a little
transform of (74) shows

$\displaystyle
a_{i}=\frac{\displaystyle\left\langle\hat{\mathcal{H}}\prod_{\begin{subarray}{c}j=1\\\
j\neq
i\end{subarray}}^{2}(\hat{\mathcal{H}}-a_{j})^{2}\right\rangle}{\displaystyle\left\langle\prod_{\begin{subarray}{c}j=1\\\
j\neq
i\end{subarray}}^{2}(\hat{\mathcal{H}}-a_{j})^{2}\right\rangle},~{}~{}(i=1,2),$
(75)

from which it is straightforward to show

$\displaystyle\langle\hat{\mathcal{H}}\rangle\geq a_{i}\geq
E_{0},~{}~{}\text{and}~{}~{}s+a_{i}\geq 0,(i=1,2),$ (76)

therefore all the terms in (70) are non-negative, and

$\displaystyle R_{1}(s,a_{1})\geq R_{2}(s,a_{1},a_{2}).$ (77)

(76) and (77) tell us that by keeping recursively refining
$\frac{1}{s+\hat{\mathcal{H}}}$ using identity (68), we are able to get
tighter upper bounds to $E_{0}$. This then constitutes the basic idea of PDS
approach. For example, if we recursively expand (68) $K$ times
($K\in\mathbb{N},K\geq 2$) with each time introducing a new parameter $a_{i}$,
$i=1,\ldots,K$), we will get a more refined non-negative residual term
$R_{K}(s,a_{1},\ldots,a_{K})$

$\displaystyle R_{K}(s,a_{1},\ldots,a_{K})$
$\displaystyle=\left\langle\frac{1}{s+\hat{\mathcal{H}}}\prod_{i=1}^{K}\frac{(\hat{\mathcal{H}}-a_{i})^{2}}{(s+a_{i})^{2}}\right\rangle,$
(78)

and

$\displaystyle R_{K-1}(s,a_{1},\ldots,a_{K-1})\geq
R_{K}(s,a_{1},\ldots,a_{K}).$ (79)

The $K$-th order PDS formalism, PDS($K$), is then associated with determining
the introduced $K$ real parameters $(a_{1}^{(K)},\ldots,a_{K}^{(K)})$ that
minimize the value of $R_{K}$ through

$\frac{\partial R_{K}(s,a_{1}^{(K)},\ldots,a_{K}^{(K)})}{\partial
a_{i}^{(K)}}=0,~{}~{}(i=1,\ldots,K),$ (80)

which can be alternatively translated to

$\displaystyle\left\langle\prod_{j=1}^{K}(\hat{\mathcal{H}}-a_{j}^{(K)})\prod_{\begin{subarray}{c}m=1\\\
m\neq
i\end{subarray}}^{K}(\hat{\mathcal{H}}-a_{m}^{(K)})\right\rangle=0,~{}~{}(i=1,\ldots,K).$
(81)

with each $a_{i}^{(K)}$ ($i=1,\ldots,K$) being an strict upper bound to
$E_{0}$, and can be expressed as

$\displaystyle\langle\hat{\mathcal{H}}\rangle>\cdots>a_{i}^{(K-1)}>a_{i}^{(K)}=\frac{\displaystyle\left\langle\hat{\mathcal{H}}\prod_{\begin{subarray}{c}j=1\\\
j\neq
i\end{subarray}}^{K}(\hat{\mathcal{H}}-a_{j}^{(K)})^{2}\right\rangle}{\displaystyle\left\langle\prod_{\begin{subarray}{c}j=1\\\
j\neq
i\end{subarray}}^{K}(\hat{\mathcal{H}}-a_{j}^{(K)})^{2}\right\rangle}\geq
E_{0}.$ (82)

To numerically solve $a_{i}^{(K)}$’s in a PDS($K$) approach, take Eq. (74) as
an exmaple, we can first express Eq. (74) in a matrix form

$\displaystyle\mathbf{S}\cdot\mathbf{M}\cdot\mathbf{X}=\mathbf{0}$ (83)

with

$\displaystyle\mathbf{S}$ $\displaystyle=\left(\begin{array}[]{cc}1&-a_{1}\\\
1&-a_{2}\end{array}\right),~{}~{}\mathbf{M}=\left(\begin{array}[]{ccc}\langle\hat{\mathcal{H}}^{3}\rangle&\langle\hat{\mathcal{H}}^{2}\rangle&\langle\hat{\mathcal{H}}\rangle\\\
\langle\hat{\mathcal{H}}^{2}\rangle&\langle\hat{\mathcal{H}}\rangle&1\end{array}\right),$
(88) $\displaystyle\mathbf{X}$ $\displaystyle=\left(\begin{array}[]{c}1\\\
-a_{1}-a_{2}\\\
a_{1}a_{2}\end{array}\right),~{}~{}\mathbf{0}=\left(\begin{array}[]{c}0\\\
0\end{array}\right).$ (94)

Then, if we assume $a_{1}\neq a_{2}$, then $\det|\mathbf{S}|\neq 0$ and
$\mathbf{S}^{-1}$ exists, therefore we are able to act $\mathbf{S}^{-1}$ on
both sides of (83) to get

$\displaystyle\mathbf{M}\cdot\mathbf{X}=\mathbf{0}\Leftrightarrow\mathbf{\tilde{M}}\cdot\mathbf{\tilde{X}}=-\mathbf{\tilde{Y}}$
(95)

with

$\displaystyle\mathbf{\tilde{M}}$
$\displaystyle=\left(\begin{array}[]{cc}\langle\hat{\mathcal{H}}^{2}\rangle&\langle\hat{\mathcal{H}}\rangle\\\
\langle\hat{\mathcal{H}}\rangle&1\end{array}\right),~{}~{}\mathbf{\tilde{X}}=\left(\begin{array}[]{c}-a_{1}-a_{2}\\\
a_{1}a_{2}\end{array}\right),$ (100) $\displaystyle\mathbf{\tilde{Y}}$
$\displaystyle=\left(\begin{array}[]{c}\langle\hat{\mathcal{H}}^{3}\rangle\\\
\langle\hat{\mathcal{H}}^{2}\rangle\end{array}\right).$ (103)

It can be easily seen that given $\langle\hat{\mathcal{H}}^{n}\rangle$
($n=1,2,3$), $\mathbf{\tilde{X}}$ (or equivalently $\mathbf{X}$) can be solved
from (95), and $a_{i}$ ($i=1,2$) are then the solution of the polynomial

$\displaystyle\text{Poly}(a)=\sum_{i=0}^{2}X_{i}a^{2-i}.$ (104)

Now generalizing (83)-(104) we get the working equation for the PDS($K$)
approach, where one need to first solve (95) for an auxiliary vector
$\mathbf{\tilde{X}}=(X_{1},\cdots,X_{K})^{T}$ with matrix elements defined as

$\displaystyle\left\\{\begin{array}[]{cl}M_{ij}&=\langle H^{2K-i-j}\rangle\\\
Y_{i}&=\langle H^{2K-i}\rangle\end{array}\right.,~{}~{}(i,j=1,\cdots,K),$
(107)

and then solve the following polynomial

$\text{Poly}_{K}(a)=a^{K}+\sum_{i=1}^{K}\tilde{X}_{i}a^{K-i},$ (108)

for $(a_{1}^{(K)},\ldots,a_{K}^{(K)})$ that provide upper bounds to the exact
ground and excited state energies of the Hamiltonian characterized by either
discrete or continuous spectral resolutions, or both.

###  III.4 Variational simulation of real and imaginary time dynamics

In the RTE and ITE, the target quantum state is approximated through the
evolutions of the trial state as shown in Eqs. (36) and (55). In practice, the
trial state can usually be prepared by applying a sequence of parametrized
gates to the initial state $|0\rangle$. For example, we can write

$\displaystyle|\phi\rangle=|\phi(\vec{\theta})\rangle=U_{n}(\theta_{n})\cdots
U_{k}(\theta_{k})\cdots U_{1}(\theta)|0\rangle,$ (109)

where $U(\vec{\theta})=U_{n}(\theta_{n})\cdots U_{k}(\theta_{k})\cdots
U_{1}(\theta)$ with $\vec{\theta}=(\theta_{1},\cdots,\theta_{n})$ and
$U_{k}(\theta_{k})$ the $k$-th one-/two-qubit gate parametrized by rotation
$\theta_{k}$. In the conventional VQE, based on the measurement outcome of the
trial state, classical optimization routines are then employed to find optimal
$\vec{\theta}$ that minimizes the cost function defined as the measurement
outcome of the quantum state $|\phi(\vec{\theta})\rangle$,

$\displaystyle
E_{\min}=\min_{\vec{\theta}}E(\vec{\theta})=\min_{\vec{\theta}}\langle\phi(\vec{\theta})|\hat{\mathcal{H}}|\phi(\vec{\theta})\rangle,$
(110)

which gives a lowest upper bound to the true ground state state energy of the
many-body Hamiltonian $\hat{\mathcal{H}}$ in the parameter space. Similar
variational procedure can also be applied following the RTE and ITE approaches
by replacing the $|\phi\rangle$ in the VQE with the time-evolved state
approximated as a linear combination of Krylov bases as shown in Eqs. (36) and
(55). There are at least two advantages of combining variational optimization
with the time evolution of the quantum state. First of all, since the
variational optimization also contributes to bring down the cost function, the
combination would help reduce the number of Krylov bases employed to evolve
the trial state, which in turn helps reduce the measurement requirement.
Second, the Krylov bases employed to improve the fidelity of the trial state
helps restructure the potential energy surface in the same parameter space,
which in turn offers an effective approach to navigate the dynamics to be free
from getting trapped in the local minima of the old potential energy surface.

These advantages have been demonstrated by Peng and Kowalski in a recent study
Peng and Kowalski (2021), where Eq. (104) is effectively treated as an
alternative energy functional for the variational quantum solver whose
analytical energy derivative $\nabla E(\vec{\theta})$ can be exploted to drive
the optimization

$\displaystyle\vec{\theta}_{k+1}=\vec{\theta}_{k}-\eta\mathcal{R}^{-1}(\vec{\theta})\nabla
E(\vec{\theta}).$ (111)

Here $\eta$ is the learning rate for updating $\vec{\theta}$ at step $k$, and
$\mathcal{R}(\vec{\theta}_{k})$ is the Riemannian metric matrix at
$\vec{\theta}_{k}$ that is flexible to characterize the singular point of the
parameter space and is essentially related to the indistinguishability of
$E(\vec{\theta})$ Yamamoto (2019). Note that similar formulation is quite
often used in, for example, gradient descent approach Lemaéchal (2012) with
the metric being often replaced by the Hessian matrix. Different from the
conventional gradient descent approach, Eq. (111) has its origin in the
general nonlinear optimization framework featuring natural gradient for
targeting machine learning problems Amari (1998). Here, the natural gradient
accounts for the geometric structure of the parameter space, and is expressed
as the product of the Fisher information matrix and the gradient of the cost
function often with the attempt of circumventing the plateaus in the parameter
space McArdle  _et al._ (2019); Yamamoto (2019); Stokes  _et al._ (2020).
Analogously, in the quantum computing, the Fisher information matrix is
replaced by the quantum Fubini-Study metric to describe the curvature of the
ansatz class.

There have been discussions about deriving and applying quantum Fubini-Study
metric in the quantum simulations, and comparing its performance with
classical Fisher metric in some recent reports McArdle  _et al._ (2019);
Yamamoto (2019); Stokes  _et al._ (2020). For example, following Ref. 26 one
can show that the RTE or ITE approach can be translated to a variational
approach for a given ansatz
$|\phi_{\tau}\rangle=|\phi(\vec{\theta}(\tau))\rangle$,
$\vec{\theta}(\tau)\in\vec{\mathbb{R}}$ for targeting the target state
$|\phi\rangle$. Take ITE approach as an example, following McLachlan’s
variational principle McLachlan (1964)

$\displaystyle\delta\|\hat{\mathcal{Q}}|\phi_{\tau}\rangle\|=0,$ (112)

where

$\displaystyle\hat{\mathcal{Q}}\equiv\frac{\partial}{\partial\tau}+\hat{\mathcal{H}}-E(\tau)$
(113)

and

$\displaystyle\|\hat{\mathcal{Q}}|\phi_{\tau}\rangle\|^{2}$
$\displaystyle=\langle\phi_{\tau}|\hat{\mathcal{Q}}^{\dagger}\hat{\mathcal{Q}}|\phi_{\tau}\rangle$
$\displaystyle=\sum_{i,j}\frac{\partial\langle\phi_{\tau}|}{\partial\theta_{i}}\frac{\partial|\phi_{\tau}\rangle}{\partial\theta_{j}}\dot{\theta}_{i}\dot{\theta}_{j}+\sum_{i}\left(\frac{\partial\langle\phi_{\tau}|}{\partial\theta_{i}}(\hat{\mathcal{H}}-E(\tau))|\phi_{\tau}\rangle+\langle\phi_{\tau}|(\hat{\mathcal{H}}-E(\tau))|\frac{\partial|\phi_{\tau}\rangle}{\partial\theta_{i}}\right)\dot{\theta}_{i}+\langle\phi_{\tau}|(\hat{\mathcal{H}}-E(\tau))^{2}|\phi_{\tau}\rangle,$
(114)

with $\dot{\theta}_{i}\equiv\frac{\partial\theta_{i}}{\partial\tau}$, we have

$\displaystyle\frac{\partial\|\hat{\mathcal{Q}}|\phi_{\tau}\rangle\|^{2}}{\partial\dot{\theta}_{i}}$
$\displaystyle=\sum_{j}\left(\frac{\partial\langle\phi_{\tau}|}{\partial\theta_{i}}\frac{\partial|\phi_{\tau}\rangle}{\partial\theta_{j}}+\frac{\partial\langle\phi_{\tau}|}{\partial\theta_{j}}\frac{\partial|\phi_{\tau}\rangle}{\partial\theta_{i}}\right)\dot{\theta}_{j}+\frac{\partial\langle\phi_{\tau}|}{\partial\theta_{i}}(H-E_{\tau})|\phi_{\tau}\rangle+\langle\phi_{\tau}|(H-E{\tau})|\frac{\partial|\phi_{\tau}\rangle}{\partial\theta_{i}}=0.$
(115)

Note that since $|\phi_{\tau}\rangle$ is normalized,
$\langle\phi_{\tau}|\phi_{\tau}\rangle=1$, we also have

$\displaystyle
0=\frac{\partial\langle\phi_{\tau}|\phi_{\tau}\rangle}{\partial\theta_{i}}=\frac{\partial\langle\phi_{\tau}|}{\partial\theta_{i}}|\phi_{\tau}\rangle+\langle\phi_{\tau}|\frac{\partial|\phi_{\tau}\rangle}{\partial\theta_{i}},$
(116)

and Eq. (115) can be simplified as a linear system in matrix form

$\displaystyle\mathcal{R}\dot{\vec{\theta}}=-\nabla E(\vec{\theta})$ (117)

which is essentially the dynamics (111) with the matrix elements
$\mathcal{R}_{ij}$ given by

$\displaystyle\mathcal{R}_{ij}=\Re\left(\frac{\partial\langle\phi_{\tau}|}{\partial\theta_{i}}\frac{\partial|\phi_{\tau}\rangle}{\partial\theta_{j}}\right).$
(118)

Same dynamics and Riemannian metric might be derived from other variational
principles. As shown in Ref. 126, classically there are three variational
principles, namely the Dirac and Frenkel variational principle Dirac (1930);
Frenkel (1934), the McLachlan variational principle McLachlan (1964), and the
time-dependent variational principle Kramer and Saraceno (1981); Broeckhove
_et al._ (1988), that lead to the same evolution equation. However, when
parameters are confined to be real, the Dirac and Frenkel and the McLachlan
variational principles are equivalent, while the time-dependent variational
principle cannot lead to a nontrivial evolution of the parameters.

It is also worth mentioning that a caveat in the above derivation is that we
implicitly assume $|\phi_{\tau}\rangle=|\phi\rangle$. If there is time-
dependent global phase difference between $|\phi_{\tau}\rangle$ and
$|\phi\rangle$, their time-derivatives can be very different and the dynamics
(117) would be incorrect McArdle  _et al._ (2019). In essence, the
$\mathcal{R}$ defined in Eq. (118) is not gauge invariant thus not qualified
for measuring the quantum distance Cheng (2013). The problem can be fixed by
either explicitly introducing a time-dependent phase gate to the trial state,
or defining a gauge invariant quantum geometric tensor

$\displaystyle\mathcal{R}_{ij}=\Re\left(\frac{\partial\langle\phi_{\tau}|}{\partial\theta_{i}}\frac{\partial|\phi_{\tau}\rangle}{\partial\theta_{j}}\right)-\frac{\partial\langle\phi_{\tau}|}{\partial\theta_{i}}|\phi_{\tau}\rangle\langle\phi_{\tau}|\frac{\partial|\phi_{\tau}\rangle}{\partial\theta_{j}}$
(119)

which is essentially associated with quantum natural gradient Stokes  _et al._
(2020).

##  IV Quantum simulation of a model Hamiltonian

In the preceding sections we have given a brief review of hybrid quantum-
classical algorithms featuring quantum computed Hamiltonian moments and their
classical post-processing approaches that are in principle capable of
providing systematically improvable estimates for the target energy and
quantum state. In this section, we will go over a simple tutor to demonstrate
how one can use the algorithms of these kinds in practice for solving some
practical problems.

In this tutorial, we will study a model Hamiltonian

$\displaystyle\hat{\mathcal{H}}=J\sum_{\langle
i,j\rangle}\left(\sigma_{x}^{i}\sigma_{x}^{j}+\sigma_{y}^{i}\sigma_{y}^{j}\right)+U\sum_{\langle
i,j\rangle}\sigma_{z}^{i}\sigma_{z}^{j}+B\sum_{i}\sigma_{z}^{i}$ (120)

which, as suggested in Ref. 132, when combining with different choices of
scalars $J$, $U$, and $B$, abstracts a large range of physical topics such as
the integrability Essler and Fagotti (2016), quantum magnetism Vasiliev  _et
al._ (2018) and many-body localization Abanin and Papić (2017); Nandkishore
and Huse (2015). In the following we will show how to employ a Hamiltonian
moment based hybrid quantum-classical approach to compute the energy and
magnetism of this model Hamiltonian (120) governing four sites mapped to four
qubits.

###  IV.1 Computation preparation

To launch the quantum computation, we first need to prepare our trial state.
The following hardware efficient ansatz is used for the preparation of our
trial state

$\displaystyle|\phi_{0}\rangle=C_{z}^{0,1}C_{z}^{1,2}C_{z}^{2,3}R_{y}^{0}(\theta_{0})R_{y}^{2}(\theta_{1})R_{y}^{3}(\pi)|1100\rangle$
(121)

where $C_{z}^{p,q}$ is a controlled-Z operation with control qubit $p$ and
target qubit $q$, $R_{y}^{p}(\theta)$ is a rotation of $\theta$ on qubit $p$
around $y$-axis, and two rotations $(\theta_{0},\theta_{1})=(-2.0,1.0)$. We
employ the hybrid ITE approach for our demonstration where the Hamiltonian
moments are evaluated from IBM-Q NISQ quantum hardware, and the energy and
magnetism computations are performed using classical ITE approach on quantum
computed Hamiltonian moments. Specifically, the ITE propagator is approximated
using a $m$-order Taylor expansion,

$\displaystyle
e^{-\tau\hat{\mathcal{H}}/2}\simeq\sum_{n=0}^{m}\frac{(-\tau)^{n}}{n!}\hat{\mathcal{H}}^{n},$
(122)

and therefore $\langle
e^{-\tau\hat{\mathcal{H}}/2}\hat{\mathcal{H}}e^{-\tau\hat{\mathcal{H}}/2}\rangle$
and $\langle e^{-\tau\hat{\mathcal{H}}}\rangle$ can be approximated as linear
combinations of $\langle\hat{\mathcal{H}}^{n}\rangle$.

###  IV.2 Quantum measurement

As mentioned in preceding sections, since $\hat{\mathcal{H}}^{n}$ can be
generally expressed as a linear combination of Pauli strings, $P_{i}$’s,
$\langle\hat{\mathcal{H}}^{n}\rangle$ is then evaluated based on the
measurements of $\langle P_{i}\rangle$’s. In quantum computation, there are
generally two ways of measuring observables, direct and indirect measurement.
In the former, the measured state collapses to the measurement basis that one
can freely choose. In the latter, the measured state is not completely
destroyed. One simplest and important example is the Hadamard test Aharonov,
Jones, and Landau (2009). In the Hadamard test of $\langle\phi|U|\phi\rangle$
in the computational basis for a trial state $|\phi\rangle$ governed by a
unitary $U$, one is able to reuse the state $\frac{(I\pm
U)|\phi\rangle}{\sqrt{2}}$ after the measurement. This property can be
exploited for designing more efficient and less resource-demanding version of
some quantum algorithms (for example, iterative version of the quantum phase
estimation, see Refs. 138; 139). Here for the demonstration we first employ
the Hadamard test to measure each $\langle P_{i}\rangle$ that contributes
$\langle\hat{\mathcal{H}}^{n}\rangle$, and then to show how to group $P_{i}$’s
to enable simultaneous measurements.

As discussed in Section II.1, the Pauli cyclic rules (12) allow for some
reduction of the number of Pauli strings that need to be measured. Here, for
the four-site Heisenberg model, we found that 72 Pauli strings form a complete
basis for reformulating arbitrary $\mathcal{H}^{n}$. This finding suggests
that once $\langle P_{i}\rangle$’s are measured one can quickly (classically)
compute approximate $\langle e^{-\tau\mathcal{H}}\rangle$ and $\langle
e^{-\tau\mathcal{H}/2}\mathcal{H}e^{-\tau\mathcal{H}/2}\rangle$ using
arbitrary order expansion for any $\tau$. Moreover, the relatively few number
of required $P_{i}$’s allows the feasibility of measuring $\langle
P_{i}\rangle$’s on physical quantum computers. Here, four IBM-Q’s publicly
available devices, ibmq_quito, ibmq_santiago, ibmq_bogota, and ibmq_manila are
employed, whose topologies and noise metrics are shown in Fig. 3.

Figure 3: (a) Topologies of four IBM-Q NISQ devices employed in this work. All
the employed IBM-Q NISQ processors belong to Falcon family with basis gates
being CX, ID, RZ, SX, and X. The color of nodes implies frequency (GHz) of the
qubit or how fast a 1-qubit gate can be executed. The color of the connection
implies the gate time in nanoseconds for 2-qubit gate such as CX. (b) Average
CX error ($\overline{\epsilon}_{CX}$), average readout error
($\overline{\epsilon}_{Readout}$), average T1/T2 time
($\overline{T1}/\overline{T2}$) of the employed IBM-Q devices.

###  IV.3 Estimated ground state energy and magnetism

We estimated the ground state energies of the four-site Heisenberg model
Hamiltonian for varied $U$ and $J$ values with $B=1.0$ employing approximate
ITE approach, where $e^{-\tau\hat{\mathcal{H}}}$ is approximated by a 15th
order Taylor expansion. In all the calculations, we set $\tau=2.5$. A detailed
discussion with respect to different truncation of the Taylor expansion and
the choice of $\tau$ will be given in Section IV.4.3.

Fig. 4a exhibits the computed ground state energies of the model Hamiltonian
with varied $\\{U,J\\}\in(0,1)$, which are in an excellent agreement with the
exact solutions shown in Fig. 4b, and the average mean squared error between
the two is only 0.0025 a.u. Based on the measured $\langle P_{i}\rangle$’s we
are also able to evaluate the magnetization of the model Hamiltonian. The
magnetization operator is given by $\hat{\mathcal{M}}=\sum_{i}\sigma_{z}^{i}$.
The contour of the expectation values of magnetization with respect to the
estimated ground states under different $U$ and $J$ values is shown in Fig.
4c, which are again in great agreement with the exact solutions shown in Fig.
4d.

Figure 4: (a) Ground state energies of the model Hamiltonian with varied $U$
and $J$ values estimated by the hybrid ITE approach, and (b) the corresponding
exact ground state energies. (c) Magnetizations of the same model estimated by
the hybrid ITE approach, and (d) the corresponding exact magnetizations. In
the hybrid ITE approach, a 15-th order Taylor expansion is employed to
approximate $e^{-\tau\hat{\mathcal{H}}}$ with $\tau=2.5$, and totally 72
$\langle P_{i}\rangle$’s have been measured on the quantum computer
ibmq_quito. The measured $\langle P_{i}\rangle$’s are compared with their
analytical values in (e). The Mean Squared Error between the estimated and
exact ground state energies across over all the ($U$, $J$)’s is 0.0025 a.u.

###  IV.4 Discussions

####  IV.4.1 In comparison with VQE

Remarkably, the ground state energy and magnetization estimations in the
present study exhibit great improvement in comparison with the VQE solutions
of the same model on the real hardware reported in previous work Kandala  _et
al._ (2017). It is also worth mentioning that (i) the ansatz we used here is
much simpler and shallower, and thus less robust, than the one used in the
previous study, (ii) the measurement requirement for the model Hamiltonian is
trivial in the present study, and does not heuristically depends on the
quality of the ansatz and the number of iterations as it does in the
conventional VQE practice, and (iii) as shown in Fig. 4e the deviation between
measured and analytical $\langle P_{i}\rangle$’s could be as large as $\pm
0.3\sim 0.4$ a.u. given $\langle P_{i}\rangle\in(-1,1)$. Nevertheless, despite
of less robust ansatz employed and the noisy $\langle P_{i}\rangle$’s, the ITE
approach is able to give very accurate estimation based on trivial
measurements, thus features a robust error mitigation ability from the
algorithmic perspective. Similar feature has also been reported in the real
quantum application of other Krylov subspace methods, for example the hybrid
and quantum Lanczos approaches Suchsland  _et al._ (2021); Motta  _et al._
(2020).

####  IV.4.2 Grouping and readout error mitigation

As discussed in the preceding section, we can further reduce the number of
terms to be measured by grouping the Pauli strings that commute and performing
simultaneous measurements Yen, Verteletskyi, and Izmaylov (2020). Here, by
applying QWC grouping, we found that the number of terms to be measured can be
further reduced from 72 to 25 for the model Hamiltonian.

Another advantage of performing simultaneous measurement is that the readout
error of the simultaneous measurement could be largely mitigated through
calibration. To see how it works, suppose we perform $n$-qubit simultaneous
measurements for totally $2^{n}$ computational bases, the normalized resulting
counts from the measurement of each basis, when collected, would constitute a
matrix $\mathbf{J}$ of $2^{n}$ dimension. Under the ideal noise-free
condition, we know that $\mathbf{J}$ is an identity matrix of the same
dimension. In the real NISQ quantum computing, however, with the effect of
noise $\mathbf{J}$ deviates from the identity matrix, and the extent of
deviation implies how noisy the real measurement outcomes are. Based on this
relation, if $\mathbf{J}$ is known, we can algebraically build a connection
between any ideal counts vector $\mathbf{C}_{\rm ideal}$ and its noisy
analogue $\mathbf{C}_{\rm noisy}$ from the real measurement through

$\displaystyle\mathbf{C}_{\rm ideal}=\mathbf{J}^{-1}\mathbf{C}_{\rm noisy}.$
(123)

Figure 5: (a) and (b) Expectation values of 25 QWC bases measured from four
IBM-Q machines with five repeated executions on each machine without and with
measurement readout calibration respectively. (c) and (d) Expectation values,
as well as standard deviations, of 25 QWC bases averaged over all IBM-Q
machines and runs without and with measurement readout calibration
respectively. (e) The estimated ground state energies of the model Hamiltonian
with varied $U$’s and $J$’s based on the calibrated simultaneous measurements
of the 25 QWC bases.

Apparently, the construction of matrix $\mathbf{J}$ eats up the resource
quickly due to the exponential growth of the dimension. Nevertheless, for
small size problem like the four-site model Hamiltonian studied in the present
work, $\mathbf{J}$ is a $16\times 16$ matrix, and a direction construction of
$\mathbf{J}$ from simultaneously measuring 16 computation bases is still
feasible. Here, after grouping the commuting Pauli strings into 25 QWC bases,
we performed the simultaneous measurement in the computational basis for each
QWC basis five times on four IBM-Q machines, and then mitigated the
measurement outcomes employing (123). The raw and calibrated measurement
outcomes are exhibited in Fig. 5a-d. As can be seen, after calibration the
measurement outcomes are greatly improved showing great consistency with the
analytical ones and significantly suppressing the deviation. Based on the
measurement outcomes, the estimated ground state energies of the model
Hamiltonian are in excellent agreement with the exact solutions (see Fig. 5e),
and the MSE between the estimated and exact energies over all the ($U,J$)’s is
only 0.0045 a.u.

####  IV.4.3 Approximation in ITE

In the case of the model Hamiltonian studied here, as stated at the beginning
of this section we approximate the imaginary time operator
$e^{-\tau\mathcal{H}}$ employing truncated Taylor expansion. Generally
speaking, for $x\in[0,b]$, the truncated expansion of an exponential function
$e^{-x}$ is $\delta$-approximation when using $b+\log(\frac{1}{\delta})$ terms
in the expansion Sachdeva and Vishnoi (2014). Regarding the ploynomial
approximations to $e^{-\tau\hat{\mathcal{H}}}$, before figuring out the number
of terms in the expansion, the variable $\tau$ needs to be determined first as
it rescales the spectrum of $\hat{\mathcal{H}}$ and indirectly affects the
level of accuracy of the employed expansion. To do so, as exemplified in Fig.
6, we fix the number of terms in the expansion and evolve the energy as a
function of $\tau$ to find its optimal value that gives the lowest energy.
This is essentially a one dimensional optimization problem over $\tau$, and
can be solved classically by Golden-Section Search Kiefer (1953), or similar
methods. In the present study, we found $\tau_{\rm opt}=2.5$ providing well
converged ground state energies estimates for the given $U$ and $J$ ranges. It
worth mentioning that to improve the estimate of $e^{-\tau\hat{\mathcal{H}}}$
using polynomial expansions in the larger $\tau$ regime, one can alternatively
resort to other polynomial expansions, such as CMX as discussed in Section
III.3, or Padé expansion as originally proposed in Ref. 88 and recently
applied in Ref. 96.

Figure 6: (a) Energy evolutions and (b) corresponding error percentages of the
four-site Heisenberg model Hamiltonian with different $U$’s and $J$’s as
functions of imaginary time, $\tau$, estimated by the approximate ITE approach
employing a 15-th order Taylor expansion.

After fixing $\tau=\tau_{\rm opt}$, we can test the accuracy of the truncated
Taylor expansion using different numbers of terms. As shown in Fig. 7 the
level of accuracy improves as more terms used in the truncated Taylor
expansion, and it verifies that at least a 15th order expansion is needed in
order to give accurate energy estimates for the entire $U$ and $J$ ranges.

Figure 7: Ground state energies with varied $U$’s and $J$’s estimated by the
approximate ITE approaches employing (a) 9-th order, (b) 11-th order, (c)
15-th order, and (d) 17-th order Taylor expansion, respectively. For all the
approximations, we fixed $\tau=2.5$.

##  V Conclusion and outlook

In this tutorial review, we have gone over some recent developments of hybrid
quantum-classical algorithms based upon quantum computation of Hamiltonian
moments. In particular, we have given an overview about how to rely on quantum
resources to compute Hamiltonian moments, and what classical methods can be
used following the quantum computation of the Hamiltonian moments to give rise
to energy estimation for the target state. The levels of accuracy of these
hybrid quantum-classical approaches are usually systematically improvable,
even under the condition of working with noisy Hamiltonian moments, and thus
provide error mitigation at the algorithmic level regardless of detailed
description of the noise. Essentially, the algorithmic error mitigation can be
attributed to the fact that these hybrid approaches are exploring the Krylov
subspace of the entire Hilbert space, and if the ansatz is not totally
orthogonal the true state, the power iteration would in principle guarantee a
convergence given sufficiently high order Hamiltonian moments. On the other
hand, what comes with the systematically improvable accuracy that is based on
the cumulative construction of Krylov subspace is the price paid for obtaining
higher order Hamiltonian moments, which will significantly increase the
measurement requirement and bring hurdles for targeting larger quantum
systems. The ongoing efforts are focused on either introducing some more
efficient grouping rules to reduce the number of terms to be measured, or
polishing the quality of the ansatz based on some physical and chemical
inspiration such that lower order Hamiltonian moments will be sufficient for
approaching certain level of accuracy. The latter can even be combined with
variational quantum solver for the same purpose, which also helps the
variational solver circumvent the barren plateaus. Another active direction is
to combine these techniques with complete active space and/or Hamiltonian
downfolding techniques (for recent studies, see e.g. Refs. 141; 142) that can
target larger systems. All these combinations are worth exploring in the NISQ
era. Ultimately, from the accuracy and efficiency improvements brought by
these algorithms and/or their combinations with some other techniques, the
critical question becomes which problems are most benefited from quantum
computed Hamiltonian moments and how much quantum advantage can be eventually
gained in these contexts.

##  VI Acknowledgement

J. C. A., T. K., and B. P. acknowledge Dr. Karol Kowalski and Dr. Niri Govind
for their helpful discussions on some topics during the preparation of this
paper.

##  VII Funding Information

J. C. A. was supported in part by the U.S. Department of Energy, Office of
Science, Office of Workforce Development for Teachers and Scientists (WDTS)
under the Science Undergraduate Laboratory Internships Program (SULI). T. K.
is supported by the U. S. Department of Energy, Office of Science, Office of
Workforce Development for Teachers and Scientists, Office of Science Graduate
Student Research (SCGSR) program. The SCGSR program is administered by the Oak
Ridge Institute for Science and Education (ORISE) for the DOE. ORISE is
managed by ORAU under contract number DE-SC0014664. B. P. was supported by the
“Embedding QC into Many-body Frameworks for Strongly Correlated Molecular and
Materials Systems” project, which is funded by the U.S. Department of Energy,
Office of Science, Office of Basic Energy Sciences (BES), the Division of
Chemical Sciences, Geosciences, and Biosciences. This material is based upon
work supported by the U.S. Department of Energy, Office of Science, National
Quantum Information Science Research Centers.

##  VIII Data Availability Statement

The data and code that support the findings of this study are available at
https://github.com/jaulicino/HamiltonianMoment/

## References

  * Feynman (1982) R. P. Feynman, “Simulating physics with computers,” Int. J. Theor. Phys. 21, 467–488 (1982).
  * David Sherrill and Schaefer (1999) C. David Sherrill and H. F. Schaefer, “The configuration interaction method: Advances in highly correlated approaches,” (Academic Press, 1999) pp. 143–269.
  * Bartlett and Musiał (2007) R. J. Bartlett and M. Musiał, “Coupled-cluster theory in quantum chemistry,” Rev. Mod. Phys. 79, 291–352 (2007).
  * Sugisaki  _et al._ (2019) K. Sugisaki, S. Nakazawa, K. Toyota, K. Sato, D. Shiomi, and T. Takui, “Quantum chemistry on quantum computers: A method for preparation of multiconfigurational wave functions on quantum computers without performing post-hartree–fock calculations,” ACS Central Sci. 5, 167–175 (2019).
  * Aspuru-Guzik  _et al._ (2005) A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon, “Simulated quantum computation of molecular energies,” Science 309, 1704–1707 (2005).
  * O’Malley  _et al._ (2016) P. J. J. O’Malley, R. Babbush, I. D. Kivlichan, J. Romero, J. R. McClean, R. Barends, J. Kelly, P. Roushan, A. Tranter, N. Ding, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. G. Fowler, E. Jeffrey, E. Lucero, A. Megrant, J. Y. Mutus, M. Neeley, C. Neill, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, P. V. Coveney, P. J. Love, H. Neven, A. Aspuru-Guzik, and J. M. Martinis, “Scalable quantum simulation of molecular energies,” Phys. Rev. X 6, 031007 (2016).
  * Linke  _et al._ (2017) N. M. Linke, D. Maslov, M. Roetteler, S. Debnath, C. Figgatt, K. A. Landsman, K. Wright, and C. Monroe, “Experimental comparison of two quantum computing architectures,” Proc. Natl. Acad. Sci. U. S. A. 114, 3305–3310 (2017).
  * Kandala  _et al._ (2017) A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, “Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets,” Nature 549, 242–246 (2017).
  * Dumitrescu  _et al._ (2018) E. F. Dumitrescu, A. J. McCaskey, G. Hagen, G. R. Jansen, T. D. Morris, T. Papenbrock, R. C. Pooser, D. J. Dean, and P. Lougovski, “Cloud quantum computing of an atomic nucleus,” Phys. Rev. Lett. 120, 210501 (2018).
  * Klco  _et al._ (2018) N. Klco, E. F. Dumitrescu, A. J. McCaskey, T. D. Morris, R. C. Pooser, M. Sanz, E. Solano, P. Lougovski, and M. J. Savage, “Quantum-classical computation of schwinger model dynamics using quantum computers,” Phys. Rev. A 98, 032331 (2018).
  * Colless  _et al._ (2018) J. I. Colless, V. V. Ramasesh, D. Dahlen, M. S. Blok, M. E. Kimchi-Schwartz, J. R. McClean, J. Carter, W. A. de Jong, and I. Siddiqi, “Computation of molecular spectra on a quantum processor with an error-resilient algorithm,” Phys. Rev. X 8, 011021 (2018).
  * McCaskey  _et al._ (2019) A. J. McCaskey, Z. P. Parks, J. Jakowski, S. V. Moore, T. D. Morris, T. S. Humble, and R. C. Pooser, “Quantum chemistry as a benchmark for near-term quantum computers,” npj Quantum Inf. 5, 98 (2019).
  * Preskill (2018) J. Preskill, “Quantum Computing in the NISQ era and beyond,” Quantum 2, 79 (2018).
  * Peruzzo  _et al._ (2014) A. Peruzzo, J. R. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’brien, “A variational eigenvalue solver on a photonic quantum processor,” Nat. Commun. 5, 4213 (2014).
  * McClean  _et al._ (2016) J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, “The theory of variational hybrid quantum-classical algorithms,” New J. Phys. 18, 023023 (2016).
  * Romero  _et al._ (2018) J. Romero, R. Babbush, J. R. McClean, C. Hempel, P. J. Love, and A. Aspuru-Guzik, “Strategies for quantum computing molecular energies using the unitary coupled cluster ansatz,” Quantum Sci. Technol. 4, 014008 (2018).
  * Shen  _et al._ (2017) Y. Shen, X. Zhang, S. Zhang, J.-N. Zhang, M.-H. Yung, and K. Kim, “Quantum implementation of the unitary coupled cluster for simulating molecular electronic structure,” Phys. Rev. A 95, 020501 (2017).
  * Kandala  _et al._ (2019) A. Kandala, K. Temme, A. D. Corcoles, A. Mezzacapo, J. M. Chow, and J. M. Gambetta, “Error mitigation extends the computational reach of a noisy quantum processor,” Nature 567, 491–495 (2019).
  * Huggins  _et al._ (2020) W. J. Huggins, J. Lee, U. Baek, B. O’Gorman, and K. B. Whaley, “A non-orthogonal variational quantum eigensolver,” New J. Phys. 22, 073009 (2020).
  * Farhi, Goldstone, and Gutmann (2014) E. Farhi, J. Goldstone, and S. Gutmann, “A quantum approximate optimization algorithm,” arXiv preprint arXiv:1411.4028 (2014).
  * Bharti  _et al._ (2021) K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S.Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, W.-K. Mok, S. Sim, L.-C. Kwek, and A. Aspuru-Guzik, “Noisy intermediate-scale quantum (nisq) algorithms,” arXiv preprint arXiv:2101.08448 (2021).
  * Albash and Lidar (2018) T. Albash and D. A. Lidar, “Adiabatic quantum computation,” Rev. Mod. Phys. 90, 015002 (2018).
  * Aaronson and Arkhipov (2011) S. Aaronson and A. Arkhipov, “Proceedings of the forty-third annual acm symposium on theory of computing,” (2011) pp. 333–342.
  * Trabesinger (2012) A. Trabesinger, “Quantum simulation,” Nat. Phys. 8, 263 (2012).
  * Georgescu, Ashhab, and Nori (2014) I. M. Georgescu, S. Ashhab, and F. Nori, “Quantum simulation,” Rev. Mod. Phys. 86, 153–185 (2014).
  * McArdle  _et al._ (2019) S. McArdle, T. Jones, S. Endo, Y. Li, S. C. Benjamin, and X. Yuan, “Variational ansatz-based quantum simulation of imaginary time evolution,” npj Quantum Inf. 5, 1–6 (2019).
  * Motta  _et al._ (2020) M. Motta, C. Sun, A. T. K. Tan, M. J. O’Rourke, E. Ye, A. J. Minnich, F. G. S. L. Brandão, and G. K.-L. Chan, “Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution,” Nat. Phys. 16, 205–210 (2020).
  * Parrish and McMahon (2019) R. M. Parrish and P. L. McMahon, “Quantum filter diagonalization: Quantum eigendecomposition without full quantum phase estimation,” arXiv preprint arXiv:1909.08925 (2019).
  * Kyriienko (2020) O. Kyriienko, “Quantum inverse iteration algorithm for programmable quantum simulators,” npj Quantum Inf. 6, 1–8 (2020).
  * McClean  _et al._ (2018) J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, “Barren plateaus in quantum neural network training landscapes,” Nat. Commun. 9, 4812 (2018).
  * Cerezo  _et al._ (2021) M. Cerezo, A. Sone, T. Volkoff, L. Cincio, and P. J. Coles, “Cost function dependent barren plateaus in shallow parametrized quantum circuits,” Nat. Commun. 12, 1791 (2021).
  * Cerezo and Coles (2021) M. Cerezo and P. J. Coles, “Higher order derivatives of quantum neural networks with barren plateaus,” Quantum Sci. Technol. 6, 035006 (2021).
  * Herasymenko and O’Brien (2019) Y. Herasymenko and T. E. O’Brien, “A diagrammatic approach to variational quantum ansatz construction,” arXiv preprint arXiv:1907.08157 (2019).
  * Bonet-Monroig, Babbush, and O’Brien (2020) X. Bonet-Monroig, R. Babbush, and T. E. O’Brien, “Nearly optimal measurement scheduling for partial tomography of quantum states,” Phys. Rev. X 10, 031064 (2020).
  * Endo  _et al._ (2021) S. Endo, Z. Cai, S. C. Benjamin, and X. Yuan, “Hybrid quantum-classical algorithms and quantum error mitigation,” J. Phys. Soc. Jpn. 90, 032001 (2021).
  * Jordan and Wigner (1928) P. Jordan and E. Wigner, “Über das paulische Äquivalenzverbot,” Z. Physik 47, 631–651 (1928).
  * Bravyi and Kitaev (2002) S. B. Bravyi and A. Y. Kitaev, “Fermionic quantum computation,” Ann. Phys. 298, 210–226 (2002).
  * Seeley, Richard, and Love (2012) J. T. Seeley, M. J. Richard, and P. J. Love, “The bravyi-kitaev transformation for quantum computation of electronic structure,” J. Chem. Phys. 137, 224109 (2012).
  * Tranter  _et al._ (2015) A. Tranter, S. Sofia, J. Seeley, M. Kaicher, J. McClean, R. Babbush, P. V. Coveney, F. Mintert, F. Wilhelm, and P. J. Love, “The bravyi–kitaev transformation: Properties and applications,” Int. J. Quantum Chem. 115, 1431–1441 (2015).
  * Tyldesley (1975) J. R. Tyldesley,  _An Introduction to Tensor Analysis for Engineers and Applied Scientists_ (Longman, 1975).
  * Kowalski and Peng (2020) K. Kowalski and B. Peng, “Quantum simulations employing connected moments expansions,” J. Chem. Phys. 153, 201102 (2020).
  * Suchsland  _et al._ (2021) P. Suchsland, F. Tacchino, M. H. Fischer, T. Neupert, P. K. Barkoutsos, and I. Tavernelli, “Algorithmic Error Mitigation Scheme for Current Quantum Processors,” Quantum 5, 492 (2021).
  * Karp (2010) R. M. Karp, “Reducibility among combinatorial problems,” in  _50 Years of Integer Programming 1958-2008: From the Early Years to the State-of-the-Art_, edited by M. Jünger, T. M. Liebling, D. Naddef, G. L. Nemhauser, W. R. Pulleyblank, G. Reinelt, G. Rinaldi, and L. A. Wolsey (Springer Berlin Heidelberg, Berlin, Heidelberg, 2010) pp. 219–241.
  * Gokhale  _et al._ (2020) P. Gokhale, O. Angiuli, Y. Ding, K. Gui, T. Tomesh, M. Suchara, M. Martonosi, and F. T. Chong, “$o(n^{3})$ measurement cost for variational quantum eigensolver on molecular hamiltonians,” IEEE Trans. Qunatum Eng. 1, 1–24 (2020).
  * Yen, Verteletskyi, and Izmaylov (2020) T.-C. Yen, V. Verteletskyi, and A. F. Izmaylov, “Measuring all compatible operators in one series of single-qubit measurements using unitary transformations,” J. Chem. Theory Comput. 16, 2400–2409 (2020).
  * Vallury  _et al._ (2020) H. J. Vallury, M. A. Jones, C. D. Hill, and L. C. L. Hollenberg, “Quantum computed moments correction to variational estimates,” Quantum 4, 373 (2020).
  * Claudino  _et al._ (2021) D. Claudino, B. Peng, N. P. Bauman, K. Kowalski, and T. S. Humble, “Improving the accuracy and efficiency of quantum connected moments expansions,” Quantum Sci. Technol. 6, 034012 (2021).
  * Schwinger (1960) J. Schwinger, “Unitary operator bases,” Proc. Nat. Acad. Sci. USA 46, 570–579 (1960).
  * Klappenecker and Rötteler (2004) A. Klappenecker and M. Rötteler, “Constructions of mutually unbiased bases,” in  _Finite Fields and Applications_ , edited by G. L. Mullen, A. Poli, and H. Stichtenoth (Springer Berlin Heidelberg, Berlin, Heidelberg, 2004) pp. 137–144.
  * Altepeter, James, and Kwiat (2004) J. B. Altepeter, D. F. V. James, and P. G. Kwiat, “4 qubit quantum state tomography,” in  _Quantum State Estimation_, edited by M. Paris and J. Řeháček (Springer Berlin Heidelberg, Berlin, Heidelberg, 2004) pp. 113–145.
  * Lawrence, Brukner, and Zeilinger (2002) J. Lawrence, i. c. v. Brukner, and A. Zeilinger, “Mutually unbiased binary observable sets on n qubits,” Phys. Rev. A 65, 032320 (2002).
  * Izmaylov  _et al._ (2020) A. F. Izmaylov, T.-C. Yen, R. A. Lang, and V. Verteletskyi, “Unitary partitioning approach to the measurement problem in the variational quantum eigensolver method,” J. Chem. Theory Comput. 16, 190–195 (2020).
  * Huggins  _et al._ (2021) W. J. Huggins, J. R. McClean, N. C. Rubin, Z. Jiang, N. Wiebe, K. B. Whaley, and R. Babbush, “Efficient and noise resilient measurements for quantum chemistry on near-term quantum computers,” npj Quantum Inf. 7 (2021), 10.1038/s41534-020-00341-7.
  * Gonthier  _et al._ (2020) J. F. Gonthier, M. D. Radin, C. Buda, E. J. Doskocil, C. M. Abuan, and J. Romero, “Identifying challenges towards practical quantum advantage through resource estimation: the measurement roadblock in the variational quantum eigensolver,” (2020), arXiv:2012.04001 [quant-ph] .
  * Rubin, Babbush, and J. (2018) N. C. Rubin, R. Babbush, and M. J., “Application of fermionic marginal constraints to hybrid quantum algorithms,” New J. Phys. 20, 053020 (2018).
  * Wecker, Hastings, and Troyer (2015) D. Wecker, M. B. Hastings, and M. Troyer, “Progress towards practical quantum variational algorithms,” Phys. Rev. A 92, 042303 (2015).
  * Huang, Kueng, and Preskill (2020) H. Y. Huang, R. Kueng, and J. Preskill, “Predicting many properties of a quantum system from very few measurements,” Nat. Phys. 16, 1050–1057 (2020).
  * Huang, Kueng, and Preskill (2021) H.-Y. Huang, R. Kueng, and J. Preskill, “Efficient estimation of pauli observables by derandomization,” Phys. Rev. Lett. 127, 030503 (2021).
  * Aaronson (2018) S. Aaronson, “Shadow tomography of quantum states,” in  _Proceedings of the 50th Annual ACM SIGACT Symposium on Theory of Computing_, STOC 2018 (Association for Computing Machinery, New York, NY, USA, 2018) pp. 325–338.
  * Struchalin  _et al._ (2021) G. I. Struchalin, Y. A. Zagorovskii, E. V. Kovlakov, S. S. Straupe, and S. P. Kulik, “Experimental estimation of quantum state properties from classical shadows,” PRX Quantum 2, 010307 (2021).
  * Chen  _et al._ (2021) S. Chen, W. Yu, P. Zeng, and S. T. Flammia, “Robust shadow estimation,” arXiv preprint arXiv:2011.09636 (2021).
  * Zhao, Rubin, and Miyake (2021) A. Zhao, N. C. Rubin, and A. Miyake, “Fermionic partial tomography via classical shadows,” Phys. Rev. Lett. 127, 110504 (2021).
  * Acharya, Saha, and Sengupta (2021) A. Acharya, S. Saha, and A. M. Sengupta, “Informationally complete povm-based shadow tomography,” arXiv preprint arXiv:2105.05992 (2021).
  * Hadfield (2021) C. Hadfield, “Adaptive pauli shadows for energy estimation,” arXiv preprint arXiv:2105.12207 (2021).
  * Hillmich  _et al._ (2021) S. Hillmich, C. Hadfield, R. Raymond, A. Mezzacapo, and R. Wille, “Decision diagrams for quantum measurements with shallow circuits,” arXiv preprint arXiv:2105.06932 (2021).
  * Zhang  _et al._ (2021) T. Zhang, J. Sun, X.-X. Fang, X. Zhang, X. Yuan, and H. Lu, “Experimental quantum state measurement with classical shadows,” arXiv preprint arXiv:2106.10190 (2021).
  * Subramanian, Brierley, and Jozsa (2019) S. Subramanian, S. Brierley, and R. Jozsa, “Implementing smooth functions of a hermitian matrix on a quantum computer,” J. Phys. Commun. 3, 065002 (2019).
  * Venegas-Andraca (2012) S. Venegas-Andraca, “Quantum walks: a comprehensive review,” Quantum Inf. Process. 11, 1015–1106 (2012).
  * Kempe (2003) J. Kempe, “Quantum random walks: An introductory overview,” Contemp. Phys. 44, 307–327 (2003).
  * Patel and Priyadarsini (2018) A. Patel and A. Priyadarsini, “Efficient quantum algorithms for state measurement and linear algebra applications,” Int. J. Quantum Inf. 16, 1850048 (2018).
  * Seki and Yunoki (2021) K. Seki and S. Yunoki, “Quantum power method by a superposition of time-evolved states,” PRX Quantum 2, 010333 (2021).
  * Trotter (1959) H. F. Trotter, “On the product of semi-groups of operators,” Proc. Amer. Math. Soc. 10, 545–551 (1959).
  * Suzuki (1990) M. Suzuki, “Fractal decomposition of exponential operators with applications to many-body theories and monte carlo simulations,” Phys. Lett. A 146, 319–323 (1990).
  * Richardson and Gaunt (1927) L. F. Richardson and J. A. Gaunt, “Viii. the deferred approach to the limit,” Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character 226, 299–361 (1927).
  * Temme, Bravyi, and Gambetta (2017) K. Temme, S. Bravyi, and J. M. Gambetta, “Error mitigation for short-depth quantum circuits,” Phys. Rev. Lett. 119, 180509 (2017).
  * Sachdeva and Vishnoi (2014) S. Sachdeva and N. K. Vishnoi,  _Faster Algorithms via Approximation Theory_ (Now Foundations and Trends, 2014).
  * Childs, Kothari, and Somma (2017) A. M. Childs, R. Kothari, and R. D. Somma, “Quantum algorithm for systems of linear equations with exponentially improved dependence on precision,” SIAM Journal on Computing 46, 1920–1950 (2017).
  * Ekert  _et al._ (2002) A. K. Ekert, C. M. Alves, D. K. L. Oi, M. Horodecki, P. Horodecki, and L. C. Kwek, “Direct estimations of linear and nonlinear functionals of a quantum state,” Phys. Rev. Lett. 88, 217901 (2002).
  * Higgott, Wang, and Brierley (2019) O. Higgott, D. Wang, and S. Brierley, “Variational Quantum Computation of Excited States,” Quantum 3, 156 (2019).
  * Mitarai and Fujii (2019) K. Mitarai and K. Fujii, “Methodology for replacing indirect measurements with direct measurements,” Phys. Rev. Research 1, 013006 (2019).
  * Ge, Tura, and Cirac (2019) Y. Ge, J. Tura, and J. I. Cirac, “Faster ground state preparation and high-precision ground energy estimation with fewer qubits,” J. Math. Phys. 60, 1–25 (2019).
  * Brassard  _et al._ (2002) G. Brassard, P. Hoyer, M. Mosca, and A. Tapp, “Quantum amplitude amplification and estimation,” in  _Quantum Computation and Quantum Information_, Vol. 305, edited by J. S. J. Lomonaco (AMS Contemporary Mathematics, 2002) pp. 53–74.
  * Berry  _et al._ (2015) D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, “Simulating hamiltonian dynamics with a truncated taylor series,” Phys. Rev. Lett. 114, 090502 (2015).
  * Low and Chuang (2017) G. H. Low and I. L. Chuang, “Optimal hamiltonian simulation by quantum signal processing,” Phys. Rev. Lett. 118, 010501 (2017).
  * Gilyén  _et al._ (2019) A. Gilyén, Y. Su, G. H. Low, and N. Wiebe, “Quantum singular value transformation and beyond: Exponential improvements for quantum matrix arithmetics,” in  _Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing_, STOC 2019 (Association for Computing Machinery, New York, NY, USA, 2019) p. 193–204.
  * Childs and Wiebe (2012) A. M. Childs and N. Wiebe, “Hamiltonian simulation using linear combinations of unitary operations,” Quantum Inf. Comput. 12 (2012), 10.26421/qic12.11-12.
  * Lanczos (1950) C. Lanczos, “An iteration method for the solution of the eigenvalue problem of linear differential and integral operators,” J. Res. Natl. Bur. Stand. B 45, 255–283 (1950).
  * Horn and Weinstein (1984) D. Horn and M. Weinstein, “The $t$ expansion: A nonperturbative analytic tool for hamiltonian systems,” Phys. Rev. D 30, 1256–1270 (1984).
  * Hollenberg (1993) L. C. L. Hollenberg, “Plaquette expansion in lattice hamiltonian models,” Phys. Rev. D 47, 1640–1644 (1993).
  * Witte and Hollenberg (1994) N. S. Witte and L. C. L. Hollenberg, “Plaquette expansion proof and interpretation,” Z. Phys. B 95 (1994).
  * Hollenberg and Witte (1996) L. C. L. Hollenberg and N. S. Witte, “Analytic solution for the ground-state energy of the extensive many-body problem,” Phys. Rev. B 54, 16309–16312 (1996).
  * van Doorn (1987) E. A. van Doorn, “Representations and bounds for zeros of orthogonal polynomials and eigenvalues of sign-symmetric tri-diagonal matrices,” J. Approx. Theory 51, 254–266 (1987).
  * Ismail and Li (1992) M. E. H. Ismail and X. Li, “Bound on the extreme zeros of orthogonal polynomials,” Proc. Amer. Math. Soc. 115, 131–140 (1992).
  * Hollenberg and Witte (1994) L. C. L. Hollenberg and N. S. Witte, “General nonperturbative estimate of the energy density of lattice hamiltonians,” Phys. Rev. D 50, 3382–3386 (1994).
  * Pavarini  _et al._ (2011) E. Pavarini, E. Koch, A. Lichtenstein, and D. Vollhardt, eds.,  _The LDA+DMFT approach to strongly correlated materials_ , Schriften des Forschungszentrums Jülich. Reihe modeling and simulation, Vol. 1 (Forschungszenrum Jülich GmbH Zentralbibliothek, Verlag, Jülich, 2011).
  * Guzman and Lacroix (2021) E. A. R. Guzman and D. Lacroix, “Predicting ground state, excited states and long-time evolution of many-body systems from short-time evolution on a quantum computer,” arXiv preprint arXiv:2104.08181 (2021).
  * McClean  _et al._ (2017) J. R. McClean, M. E. Kimchi-Schwartz, J. Carter, and W. A. de Jong, “Hybrid quantum-classical hierarchy for mitigation of decoherence and determination of excited states,” Phys. Rev. A 95, 042308 (2017).
  * Parrish  _et al._ (2019) R. M. Parrish, E. G. Hohenstein, P. L. McMahon, and T. J. Martínez, “Quantum computation of electronic transitions using a variational quantum eigensolver,” Phys. Rev. Lett. 122, 230401 (2019).
  * Nakanishi, Mitarai, and Fujii (2019) K. M. Nakanishi, K. Mitarai, and K. Fujii, “Subspace-search variational quantum eigensolver for excited states,” Phys. Rev. Research 1, 033062 (2019).
  * Davies  _et al._ (1980) K. T. R. Davies, H. Flocard, S. Krieger, and M. S. Weiss, “Application of the imaginary time step method to the solution of the static hartree-fock problem,” Nuc. Phys. A 342, 111–123 (1980).
  * Anderson (1975) J. B. Anderson, “A random‐walk simulation of the schrödinger equation: H+3,” J. Chem. Phys. 63, 1499–1503 (1975).
  * Anderson (1979) J. B. Anderson, “Quantum chemistry by random walk: H4 square,” Int. J. Quantum Chem. 15, 109–120 (1979).
  * Anderson (1980) J. B. Anderson, “Quantum chemistry by random walk: Higher accuracy,” J. Chem. Phys. 73, 3897–3899 (1980).
  * Cioslowski (1987) J. Cioslowski, “Connected moments expansion: A new tool for quantum many-body theory,” Phys. Rev. Lett. 58, 83–85 (1987).
  * Knowles (1987) P. J. Knowles, “On the validity and applicability of the connected moments expansion,” Chem. Phys. Lett. 134, 512–518 (1987).
  * Stubbins (1988) C. Stubbins, “Methods of extrapolating the t-expansion series,” Phys. Rev. D 38, 1942 (1988).
  * Prie  _et al._ (1994) J. D. Prie, D. Schwall, J. D. Mancini, D. Kraus, and W. J. Massano, “On the relation between the connected-moments expansion and the lanczos variational scheme,” Nuov. Cim. D 16, 433–448 (1994).
  * Mancini, Zhou, and Meier (1994) J. D. Mancini, Y. Zhou, and P. F. Meier, “Analytic properties of connected moments expansions,” Int. J. Quantum Chem. 50, 101–107 (1994).
  * Ullah (1995) N. Ullah, “Removal of the singularity in the moment-expansion formalism,” Phys. Rev. A 51, 1808–1810 (1995).
  * Mancini  _et al._ (1995) J. D. Mancini, W. J. Massano, J. D. Prie, and Y. Zhuo, “Avoidance of singularities in moments expansions: a numerical study,” Phys. Lett. A 209, 107–112 (1995).
  * Fessatidis  _et al._ (2006) V. Fessatidis, J. D. Mancini, R. Murawski, and S. P. Bowen, “A generalized moments expansion,” Phys. Lett. A 349, 320–323 (2006).
  * Fessatidis  _et al._ (2010) V. Fessatidis, F. A. Corvino, J. D. Mancini, R. K. Murawski, and J. Mikalopas, “Analytic properties of moments matrices,” Phys. Lett. A 374, 2890–2893 (2010).
  * Noga, Szabados, and Surján (2002) J. Noga, A. Szabados, and P. Surján, “On the use of connected moments expansion with coupled cluster reference,” Int. J. Mol. Sci. 3, 508–521 (2002).
  * Mancini, Prie, and Massano (1991) J. D. Mancini, J. D. Prie, and W. J. Massano, “Approximations to the ground-state energy of the single-impurity kondo model using variational and connected-moments expansions,” Phys. Rev. A 43, 1777–1782 (1991).
  * Peeters and Devreese (1984) F. M. Peeters and J. T. Devreese, “Upper bounds for the free energy. a generalisation of the bogolubov inequality and the feynman inequality,” Journal of Physics A: Mathematical and General 17, 625 (1984).
  * Soldatov (1995) A. V. Soldatov, “Generalized variational principle in quantum mechanics,” Int. J. Mod. Phys. B 9, 2899–2936 (1995).
  * Bogolubov (1947) N. Bogolubov, “On the theory of superfluidity,” J. Phys. 11, 23–32 (1947).
  * Feynman (1955) R. P. Feynman, “Slow electrons in a polar crystal,” Phys. Rev. 97, 660–665 (1955).
  * Devreese, Evrard, and Kartheuser (1975) J. T. Devreese, R. Evrard, and E. Kartheuser, “Self-consistent equation-of-motion approach for polarons,” Phys. Rev. B 12, 3353–3367 (1975).
  * Peng and Kowalski (2021) B. Peng and K. Kowalski, “Variational quantum solver employing the pds energy functional,” Quantum 5, 473 (2021).
  * Yamamoto (2019) N. Yamamoto, “On the natural gradient for variational quantum eigensolver,” arXiv preprint arXiv:1909.05074 (2019).
  * Lemaéchal (2012) C. Lemaéchal, “Cauchy and the gradient method,” Doc. Math. Extra Optimization Stories, 251–254 (2012).
  * Amari (1998) S. Amari, “Natural gradient works efficiently in learning,” Neural Comput. 10, 251–276 (1998).
  * Stokes  _et al._ (2020) J. Stokes, J. Izaac, N. Killoran, and G. Carleo, “Quantum Natural Gradient,” Quantum 4, 269 (2020).
  * McLachlan (1964) A. D. McLachlan, “A variational solution of the time-dependent schrodinger equation,” Mol. Phys. 8, 39–44 (1964).
  * Yuan  _et al._ (2019) X. Yuan, S. Endo, Q. Zhao, Y. Li, and S. C. Benjamin, “Theory of variational quantum simulation,” Quantum 3, 191 (2019).
  * Dirac (1930) P. A. M. Dirac, “Note on exchange phenomena in the thomas atom,” Math. Proc. Camb. Philos. Soc. 26, 376–385 (1930).
  * Frenkel (1934) J. Frenkel,  _Wave Mechanics: Advanced General Theory_ (Clarendon Press Oxford, 1934).
  * Kramer and Saraceno (1981) P. Kramer and M. Saraceno, eds., “The time-dependent variational principle (tdvp),” in  _Geometry of the Time-Dependent Variational Principle in Quantum Mechanics_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 1981) pp. 3–14.
  * Broeckhove  _et al._ (1988) J. Broeckhove, L. Lathouwers, E. Kesteloot, and P. Van Leuven, “On the equivalence of time-dependent variational principles,” Chem. Phys. Lett. 149, 547–550 (1988).
  * Cheng (2013) R. Cheng, “Quantum geometric tensor (fubini-study metric) in simple quantum system: A pedagogical introduction,” arXiv preprint arXiv:1012.1337 (2013).
  * Smith  _et al._ (2019) A. Smith, M. Kim, F. Pollmann, and J. Knolle, “Simulating quantum many-body dynamics on a current digital quantum computer,” npj Quantum Inf. 5 (2019), 10.1038/s41534-019-0217-0.
  * Essler and Fagotti (2016) F. H. L. Essler and M. Fagotti, “Quench dynamics and relaxation in isolated integrable quantum spin chains,” J. Stat. Mech.: Theory Exp. 2016, 064002 (2016).
  * Vasiliev  _et al._ (2018) A. Vasiliev, O. Volkova, E. Zvereva, and M. Markina, “Milestones of low-d quantum magnetism,” npj Quant. Mater. 3, 18 (2018).
  * Abanin and Papić (2017) D. A. Abanin and Z. Papić, “Recent progress in many-body localization,” Ann. Phys. (Berl.) 529, 1700169 (2017).
  * Nandkishore and Huse (2015) R. Nandkishore and D. A. Huse, “Many-body localization and thermalization in quantum statistical mechanics,” Ann. Rev. Condens. Matter Phys. 6, 15–38 (2015).
  * Aharonov, Jones, and Landau (2009) D. Aharonov, V. Jones, and Z. Landau, “A polynomial quantum algorithm for approximating the jones polynomial,” Algorithmica 55, 395–421 (2009).
  * Knill, Ortiz, and Somma (2007) E. Knill, G. Ortiz, and R. D. Somma, “Optimal quantum measurements of expectation values of observables,” Phys. Rev. A 75, 012328 (2007).
  * Dobšíček  _et al._ (2007) M. Dobšíček, G. Johansson, V. Shumeiko, and G. Wendin, “Arbitrary accuracy iterative quantum phase estimation algorithm using a single ancillary qubit: A two-qubit benchmark,” Phys. Rev. A 76, 030306 (2007).
  * Kiefer (1953) J. Kiefer, “Sequential minimax search for a maximum,” Proc. Am. Math. Soc. 4, 502–506 (1953).
  * Bauman  _et al._ (2019) N. P. Bauman, E. J. Bylaska, S. Krishnamoorthy, G. H. Low, N. Wiebe, C. E. Granade, M. Roetteler, M. Troyer, and K. Kowalski, “Downfolding of many-body hamiltonians using active-space models: Extension of the sub-system embedding sub-algebras approach to unitary coupled cluster formalisms,” J. Chem. Phys. 151, 014107 (2019).
  * Mniszewski  _et al._ (2021) S. M. Mniszewski, P. A. Dub, S. Tretiak, P. M. Anisimov, Y. Zhang, and C. F. A. Negre, “Reduction of the molecular hamiltonian matrix using quantum community detection,” Sci. Rep. 11, 4099 (2021).