# Variational quantum solver employing the PDS energy functional

Bo Peng<EMAIL_ADDRESS>Physical and Computational Science Division, Pacific
Northwest National Laboratory, Richland, Washington 99354, United States of
America  Karol Kowalski<EMAIL_ADDRESS>Physical and Computational
Science Division, Pacific Northwest National Laboratory, Richland, Washington
99354, United States of America

###### Abstract

Recently a new class of quantum algorithms that are based on the quantum
computation of the connected moment expansion has been reported to find the
ground and excited state energies. In particular, the Peeters-Devreese-
Soldatov (PDS) formulation is found variational and bearing the potential for
further combining with the existing variational quantum infrastructure. Here
we find that the PDS formulation can be considered as a new energy functional
of which the PDS energy gradient can be employed in a conventional variational
quantum solver. In comparison with the usual variational quantum eigensolver
(VQE) and the original static PDS approach, this new variational quantum
solver offers an effective approach to navigate the dynamics to be free from
getting trapped in the local minima that refer to different states, and
achieve high accuracy at finding the ground state and its energy through the
rotation of the trial wave function of modest quality, thus improves the
accuracy and efficiency of the quantum simulation. We demonstrate the
performance of the proposed variational quantum solver for toy models, H2
molecule, and strongly correlated planar H4 system in some challenging
situations. In all the case studies, the proposed variational quantum approach
outperforms the usual VQE and static PDS calculations even at the lowest
order. We also discuss the limitations of the proposed approach and its
preliminary execution for model Hamiltonian on the NISQ device.

##  1 Introduction

Quantum computing (QC) techniques attract much attention in many mathematics,
physics, and chemistry areas by providing means to address insurmountable
computational barriers for simulating quantum systems on classical
computers.[47, 61, 53, 3, 2, 43] One of the focus areas for quantum computing
is quantum chemistry, where Hamiltonians can be effectively mapped into qubit
registers. In this area, several quantum computing algorithms, including
quantum phase estimator (QPE) [38, 15, 5, 12, 58, 69, 26, 52] and variational
quantum eigensolver (VQE), [50, 44, 55, 60, 31, 32, 24, 16, 29] have been
extensively tested on benchmark systems corresponding to the description of
chemical reactions involving bond-forming and breaking processes, excited
states, and strongly correlated molecular systems. In more recent
applications, several groups reported quantum algorithms for imaginary time
evolution,[42, 46] quantum filter diagonalization,[48] quantum inverse
iteration algorithms,[36] and quantum power/moments methods. [59, 67] The main
thrust that drives this field is related to the efficient encoding of the
electron correlation effects that are needed to describe molecular systems.

Basic methodological questions related to an efficient way of incorporating
large degrees of freedom required to capture a subtle balance between static
and dynamical correlations effects still need to be appropriately addressed. A
typical way of addressing these challenges in VQE approaches is by
incorporating more and more parameters (usually corresponding to excitation
amplitudes in a broad class of unitary coupled-cluster methods [27, 4, 64, 35,
37, 17]). Unfortunately, this brute force approach is quickly stumbling into
insurmountable problems associated with the resulting quantum circuit
complexity and problems with numerical optimization procedures performed on
classical machines (the so-called barren plateau problem reported in Refs.[45,
11, 68, 10, 51, 41, 66]).

In this paper, we propose a new solution to these problems. Instead of adding
more parameters to the trial wave function, we choose to optimize a new class
of energy functionals (or quasi-functionals, where the energy is calculated as
a simple equation solution) that already encompasses information about high-
order static and dynamical correlation effects. An ideal choice for such high-
level functional is based on the Peeters, Devreese, and Soldatov (PDS)
formalism,[49, 62] where variational energy is obtained as a solution of
simple equations expressed in terms of the Hamiltonian’s moments or
expectations values of the powers of the Hamiltonians operator defined for the
trial wave function. In Ref. [34] we demonstrated that in such calculations
high-level of accuracy can be achieved even with very simple parametrization
of the trial wave functions (capturing only essential correlation effects) and
low-rank moments. We believe that merging the PDS formalism with the quantum
gradient based variational approach would be considered as a more interesting
alternative for by-passing main problems associated with the excessive number
of amplitudes that need to be included to reach the so-called chemical
accuracy.

In the following sections we will briefly introduce the PDS formalism and
describe how the PDS energy functional can be incorporated with the
minimization procedures that are based on the quantum gradient approach [25,
57, 45, 42, 71] to produce a new class of variational quantum solver (which is
called PDS($K$)-VQS for short in the rest of the paper) to target the ground
state and its energy in a quantum-classical hybrid manner. Furthermore, we
will test its performance, in particular the performance of the more
affordable lower order PDS($K$)-VQS ($K=2,3,4$) approaches combining with the
trial wave function expressed in low-depth quantum circuits, at finding the
ground state and its energy for the Hamiltonians describing toy models and H2
molecular system, as well as the strongly correlated planar H4 system, in some
challenging situations where the barren plateau problem precludes the
effective utilization of the standard VQE approach.

##  2 Method

###  2.1 PDS formalism

In this section we will give a brief description of the PDS formalism. The
detailed discussion of the PDS methodology and highly relevant connected
moment expansion (CMX) formalisms have been given in the original work[49, 62]
as well as our recent work[34, 14] and many earlier literatures (see for
example Refs. [33, 54, 39, 65, 40, 18, 19]). The many-body techniques used in
the derivation of PDS expansions originate in the effort to provide upper
bounds for the free energies, and to provide alternative re-derivation of the
Bogolubov’s [8] and Feynman’s [20] inequalities. Since the Gibbs-Bogolubov
inequality reduces to the Rayleight-Ritz variational principle in zero
temperature limit, these formulations can be directly applied to quantum
chemistry. Here we only provide an overview of basic steps involved in the
derivations of the PDS formulation.

A starting point of the studies of upper bounds for the exact ground-state
energy $E_{0}$ is the analysis of function $\Gamma(t)$ (defined for trial wave
function $|\phi\rangle$ having non-zero overlap with the ground-state wave
function)

$\Gamma(t)=\langle\phi|e^{-tH}|\phi\rangle\;,$ (1)

and its Laplace transform $f(s)$

$f(s)=\int_{0}^{+\infty}e^{-st}\Gamma(t)dt\;.$ (2)

It can be proved that, for a complex scalar $s$, Eq. (2) exists if the real
part of $s$ $\Re(s)>-E_{0}$. Under this condition, for Hamiltonian $H$ defined
by discrete energy levels $E_{i}$ and corresponding eigenvectors
$|\Psi_{i}\rangle$ ($i=0,1,\ldots,M$)

$H=\sum_{i=0}^{M}E_{i}|\Psi_{i}\rangle\langle\Psi_{i}|\;,$ (3)

$f(s)$ takes the form

$f(s)=\sum_{i=0}^{M}\frac{\omega(E_{i})}{s+E_{i}}$ (4)

where $\omega(E_{i})=|\langle\Psi_{i}|\phi\rangle|^{2}$. The PDS formalism is
based on introducing parameters into expansion (4) using a simple identity
(with a real parameter $a$)

$\displaystyle\frac{1}{s+E_{n}}$ $\displaystyle=$
$\displaystyle\frac{1}{s+a}-\frac{E_{n}-a}{(s+a)^{2}}+\frac{(E_{n}-a)^{2}}{(s+E_{n})(s+a)^{2}}\;.$
(5)

When the above identity is applied for the first time to Eq. (4) (introducing
the first parameter $a_{1}$) one gets the following expression for the $f(s)$
function

$\displaystyle f(s)$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{M}\omega(E_{i})\left[\frac{1}{s+a_{1}}-\frac{E_{i}-a_{1}}{(s+a_{1})^{2}}+\frac{(E_{i}-a_{1})^{2}}{(s+E_{i})(s+a_{1})^{2}}\right]$
(6)

The transformation (6) can be repeated $K$ times (with each time introducing a
new parameter $a_{i}$, $i=1,\ldots,K$) to reformulate the $f(s)$ function as

$f(s)=R_{K}(s,a_{1},\ldots,a_{K})+W_{K}(s,a_{1},\ldots,a_{K})\;,$ (7)

where

$\displaystyle R_{K}(s,a_{1},\ldots,a_{K})$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{M}\left[\frac{\omega(E_{i})}{s+E_{i}}\prod_{j=1}^{K}\frac{(E_{i}-a_{j})^{2}}{(s+a_{j})^{2}}\right]\geq-
E_{0}\leavevmode\nobreak\ \leavevmode\nobreak\ (\text{if}\leavevmode\nobreak\
\leavevmode\nobreak\ \Re(s)>-E_{0}),$ (8) $\displaystyle
W_{K}(s,a_{1},\ldots,a_{K})$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{M}\left\\{\omega(E_{i})\sum_{j=1}^{K}\left[\Big{(}\frac{1}{s+a_{j}}-\frac{E_{i}-a_{j}}{(s+a_{j})^{2}}\Big{)}\prod_{n=1}^{j-1}\frac{(E_{i}-a_{n})^{2}}{(s+a_{n})^{2}}\right]\right\\}\;.$
(9)

The $K$-th order PDS formalism (PDS($K$) for short henceforth) is then
associated with defining the introduced $K$ real parameters
$(a_{1},\ldots,a_{K})$ that minimize the value of
$R_{K}(s,a_{1},\ldots,a_{K})$. In this minimization process the necessary
extreme conditions are given by the system of equations

$\frac{\partial R_{K}(s,a_{1},\ldots,a_{K})}{\partial
a_{i}}=0,\leavevmode\nobreak\ \leavevmode\nobreak\ (i=1,\ldots,K),$ (10)

which can be alternatively represented by the matrix system of equations for
an auxiliary vector $\mathbf{X}=(X_{1},\cdots,X_{K})^{T}$

$\mathbf{MX}=-\mathbf{Y}.$ (11)

Here, the matrix elements of $\bf M$ and vector $\bf Y$ are defined as the
expectation values of Hamiltonian powers (i.e. moments),
$M_{ij}=\langle\phi|H^{2K-i-j}|\phi\rangle$,
$Y_{i}=\langle\phi|H^{2K-i}|\phi\rangle$ ($i,j=1,\cdots,K$) (for simplicity,
we will use the notation $\langle
H^{n}\rangle\equiv\langle\phi|H^{n}|\phi\rangle$). It can be shown that the
optimal parameters in the PDS($K$) formalism,
$(a_{1}^{(K)},\ldots,a_{K}^{(K)})$, are the roots of the polynomial
$P_{K}(\mathcal{E})$,

$P_{K}(\mathcal{E})=\mathcal{E}^{K}+\sum_{i=1}^{K}X_{i}\mathcal{E}^{K-i},$
(12)

and these roots provide upper bounds for the exact ground and excited state
energies, e.g., for the ground state energy we have

$E_{0}\leq{\rm
min}(a_{1}^{(K)},\ldots,a_{K}^{(K)})\leq\langle\phi|H|\phi\rangle\;.$ (13)

Note that, as shown in Refs.[49, 62] the PDS formalism also applies to the
Hamiltonian characterized by discrete and continuous spectral resolutions
together.

###  2.2 PDS($K$)-VQS formalism

In the variational method, we approximate the quantum state using parametrized
trial state $|\Psi\rangle\approx|\phi\rangle$. Using a quantum circuit, the
trial state can be prepared by applying a sequence of parametrized unitary
gates on the initial state $|0\rangle$,

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

($\vec{\theta}=\\{\theta_{1},\cdots,\theta_{n}\\}$). Here $U_{k}(\theta_{k})$
is the $k$-th unitary single- or two-qubit gate that is controlled by
parameter $\theta_{k}$. The goal is to approach the ground-state energy of a
many-body Hamiltonian, $H$, by finding the values of these parameters,
$\vec{\theta}$, that minimize the expectation value of the Hamiltonian

$\displaystyle
E_{\min}=\min_{\vec{\theta}}\langle\phi(\vec{\theta})|H|\phi(\vec{\theta})\rangle.$
(15)

To do this, the conventional VQE starts by constructing the ansatz
$|\phi(\vec{\theta})\rangle$ and measuring the corresponding expectation value
of the Hamiltonian using a quantum computer, and then relies on a classical
optimization routine to obtain new $\vec{\theta}$. During the parameter
optimization (or dynamics), the set of parameters that is updated at the
$k$-th step ($k>1$) can be written as

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

where
$\nabla\mathcal{E}(\vec{\theta})=\partial\mathcal{E}/\partial\vec{\theta}$ is
the energy gradient vector, and $\eta$ is the step size (or learning rate).
$\mathcal{R}(\vec{\theta})$ is the Riemannian metric matrix at $\vec{\theta}$
that is flexible to characterize the singular point in the parameter space and
is essentially related to the indistinguishability of
$\mathcal{E}(\vec{\theta})$.[71] It is worth mentioning that Eq. (5)
originates from natural gradient learning method in the general nonlinear
optimization framework especially targeting machine learning problems.[1]
Here, the natural gradient is the optimizer that accounts for the geometric
structure of the parameter space. For the curved (or nonorthonormal) parameter
manifold that exhibits the Riemannian character (e.g. in large neural
networks), natural gradient learning method is often employed to avoid the
plateaus in the parameter space.[42, 71, 63]

Note that when the parameter space is a Euclidean space with orthonormal
coordinate system the Riemannian metric tensor will reduce to the unity matrix
(see Tab. 1). In VQE setting, one can define the Riemannian metric as the
quantum Fubini-Study metric, which is the quantum analog of the Fisher
information matrix in the classical natural gradient,[1] to measure the
distance in the space of pure quantum states. The quantum Fubini-Study metric
describes the curvature of the ansatz class rather than the learning
landscape, but often performs as well as Hessian based methods (e.g. BFGS
optimizer that approximates the Hessian of the cost function using first-order
gradient, see Ref. [70] for a recent detailed discussion). There are also some
other options for the Riemannian metric including imaginary-time evolution
(ITE) or even classical Fisher metric that have been discussed in some recent
reports.[42, 71, 63] In Tab. 1, three commonly used flavors of the Riemannian
metric matrix $\mathcal{R}(\vec{\theta})$ are listed and will be used in the
following case studies. Remarkably, as pointed out in Refs. [73, 42], the
difference between natural gradient descent (NGD) and ITE accounts for the
global phase, and if introducing a time-dependent phase gate to the trial
state, the Riemannian metric employing NGD will be equivalent to the metric
employing ITE.

Riemannian metric $\mathcal{R}_{ij}(\vec{\theta})$ GD $\delta_{ij}$ NGD
$\Re\Big{(}\frac{\partial\langle\phi(\vec{\theta})|}{\partial\theta_{i}}\frac{\partial|\phi(\vec{\theta})\rangle}{\partial\theta_{j}}\Big{)}-\frac{\partial\langle\phi(\vec{\theta})|}{\partial\theta_{i}}|\phi(\vec{\theta})\rangle\langle\phi(\vec{\theta})|\frac{\partial|\phi(\vec{\theta})\rangle}{\partial\theta_{j}}$
ITE
$\Re\Big{(}\frac{\partial\langle\phi(\vec{\theta})|}{\partial\theta_{i}}\frac{\partial|\phi(\vec{\theta})\rangle}{\partial\theta_{j}}\Big{)}$

Table 1: Three Rimannian metric forms, ordinary gradient descent (GD), natural
gradient descent (NGD), and imaginary time evolution (ITE), exploited in the
present study. Figure 1: The workflow of variational quantum solver employing
the PDS energy functional.

To get the energy gradient in the PDS framework, take the derivative w.r.t.
$\theta_{i}$ on both sides of Eq. (12), and after reorganizing the terms we
can express the energy derivative as

$\displaystyle\frac{\partial\mathcal{E}}{\partial\theta_{i}}=\frac{-1}{K\mathcal{E}^{K-1}+\displaystyle\sum_{i=1}^{K-1}(K-i)X_{i}\mathcal{E}^{K-i-1}}\left(\begin{array}[]{c}\mathcal{E}^{K-1}\\\
\vdots\\\
1\end{array}\right)^{T}\frac{\partial\mathbf{X}}{\partial\theta_{i}},$ (20)

where $\frac{\partial\mathbf{X}}{\partial\theta_{i}}$ is associated with the
$\theta_{i}$-derivative of Eq. (11),

$\displaystyle\mathbf{M}\frac{\partial\mathbf{X}}{\partial\theta_{i}}=-\frac{\partial\mathbf{Y}}{\partial\theta_{i}}-\frac{\partial\mathbf{M}}{\partial\theta_{i}}\mathbf{X},$
(21)

and can be obtained by solving Eq. (21) as a linear equation with $\partial
Y_{i}/\partial\theta_{k}=\partial\langle H^{2K-i}\rangle/\partial\theta_{k}$
and $\partial M_{ij}/\partial\theta_{k}=\partial\langle
H^{2K-i-j}\rangle/\partial\theta_{k}$. Fig. 1 summarizes the workflow of
PDS($K$)-VQS, where on the classical side the PDS($K$) module includes two
steps, (i) solving two consecutive linear problems to get $\mathbf{X}$ and
$\partial\mathbf{X}/\partial\theta_{i}$, and (ii) solving for roots of
polynomial (12) and computing Eq. (20). On the quantum side, in comparison
with the conventional VQE, the present PDS($K$)-VQS infrastructure relies on
quantum circuits to measure $\langle H^{n}\rangle$ and their
$\vec{\theta}$-derivatives.

In the present work, due to the relatively small system size, we directly
exploit the Hadamard test to compute the real part of $\langle H^{n}\rangle$
for the Hamiltonians that are represented as a sum of Pauli strings. It is
worth mentioning that for typical molecular systems that can be represented by
$N$ qubits, the number of $\langle H^{n}\rangle$ measurement scales as
$\mathcal{O}(N^{4n})$, which nevertheless can be reduced once the Pauli
strings are multiplied and their expectation values are re-used as the
contributions to the higher order moments.[34] For example, as we will show
later for the H4 system that comprises 184 Pauli strings in the Hamiltonian,
the effective number of Pauli strings required for arbitrary $\langle
H^{n}\rangle$ ($n=2,3,4$) measurements can be dropped from 1842, 1843, and
1844 to 1774, 3702, and 4223, respectively, after the Pauli reduction, and the
4223 strings will not be changed for more complex $\langle H^{n}\rangle$’s
($n>4$). Similar findings have also been reported in Ref. [67], where by
grouping the Pauli strings into tensor-product basis sets the authors examined
the operator counts for $\langle H^{4}\rangle$ of Heisenberg model defined on
different lattice geometries for the number of qubits ranging from 2 up to 36,
and found that the effective number of Pauli strings to be measured drops by
several orders of magnitude with sub-linear scaling in a given number of
qubits. For larger systems, the number of measurements can be further reduced
by introducing active space and local approximation. Alternatively, one can
approximate $\langle H^{n}\rangle$ by a linear combination of the time-
evolution operators as introduced in some recent reports.[59, 7] For the
estimation of $\partial\langle H^{n}\rangle/\partial\theta_{k}$, in the
present work we limit $U_{k}(\theta_{k})$ exploited in the state preparation
to be only one-qubit rotations. Then, following Ref. [57], $\partial\langle
H^{n}\rangle/\partial\theta_{k}$ can be obtained by measuring $\langle
H^{n}\rangle$ twice using the same circuit but shifting $\theta_{k}$ by
$\pm\frac{\pi}{2}$ separately, i.e.

$\displaystyle\frac{\partial\langle
H^{n}\rangle_{(\cdots,\theta_{k},\cdots)}}{\partial\theta_{k}}=\frac{1}{2}\Big{(}\langle
H^{n}\rangle_{(\cdots,\theta_{k}+\frac{\pi}{2},\cdots)}-\langle
H^{n}\rangle_{(\cdots,\theta_{k}-\frac{\pi}{2},\cdots)}\Big{)}.$ (22)

If $\theta_{k}$ parametrizes more than one one-qubit rotations in the circuit,
then based on the product rule $\partial\langle
H^{n}\rangle/\partial\theta_{k}$ will have contributions from all one-qubit
$\theta_{k}$ rotations, each of which will be obtained by applying (22) on the
corresponding rotation.

##  3 Numerical examples

In this section, with several examples, we will demonstrate how the
PDS($K$)-VQS performs in some challenging situations, and its difference in
comparison to the conventional VQE and static PDS($K$) expansions.

###  3.1 Toy Hamoltonians

We first test the PDS($K$)-VQS on two toy Hamiltonians

$\displaystyle H_{A}$ $\displaystyle=1.5I_{4\times 4}+0.5(I_{2\times
2}\otimes\sigma_{z}-2\sigma_{z}\otimes\sigma_{z})$
$\displaystyle=\left(\begin{array}[]{cccc}1&0&0&0\\\ 0&2&0&0\\\ 0&0&3&0\\\
0&0&0&0\end{array}\right),$ (27) $\displaystyle H_{B}$
$\displaystyle=I_{4\times 4}+0.5(I_{2\times
2}\otimes\sigma_{z}-\sigma_{z}\otimes\sigma_{z})$
$\displaystyle=\left(\begin{array}[]{cccc}1&0&0&0\\\ 0&1&0&0\\\ 0&0&2&0\\\
0&0&0&0\end{array}\right),$ (32)

with ansatze

$\displaystyle|\phi_{A}(\theta_{1},\theta_{2})\rangle=\tilde{R}_{Y}^{0,1}(\theta_{2})R_{X}^{0}(\theta_{1})|00\rangle,$
$\displaystyle|\phi_{B}(\theta_{1},\theta_{2})\rangle=\tilde{R}_{Y}^{0,1}(\theta_{2})R_{X}^{0}(\theta_{1})R_{X}^{1}(\theta_{1})|01\rangle,$

that have been exploited by McArdle et al. [42] to demonstrate the performance
of different Riemannian metrics in the conventional VQE approach for finding
the ground-state energy of the same Hamiltoinans. Here,
$\tilde{R}_{Y}^{p,q}(\theta)$ is a controlled $Y$ rotation of $\theta$ with
control qubit $p$ and target qubit $q$, and $R_{X}^{p}(\theta)$ is a rotation
of $\theta$ on qubit $p$ around the $x$-axis. The rotation about the $j$-axis
is defined as $R_{\sigma_{j}}(\theta)=e^{-\text{i}\theta\sigma_{j}/2}$ with
$\sigma_{j}$ being one of the Pauli spin matrices.

Figure 2: Variational trajectories on the PDS(2) energy surface (left panels)
and original potential energy surface (right panels) discovering the ground
state energy of Hamiltonian, $H_{A}$, explored by gradient descent (top
panels) and natural gradient descent/imaginary time evolution (bottom panels).
On the background energy surfaces, the dark blue and white colors correspond
to the global maximum and minimum energies, respectively. The arrows indicate
the trajectories of the dymanics, and are colored green if the trajectory
converges to the ground state energy, and red otherwise. The step size
$\eta=0.05$ in all the calculations. Figure 3: Variational trajectories on the
PDS(2) energy surface (left panels) and original potential energy surfaces
(right panels) discovering the ground state energy of Hamiltonian, $H_{B}$,
explored by gradient descent (top panels), natural gradient descent (middle
panels), and variational imaginary time (bottom panels). On the background
energy surfaces, the dark blue and white colors correspond to the global
maximum and minimum energies, respectively. The arrows indicate the
trajectories of the methods, and are colored green if the trajectory converges
to the true ground state energy, and red otherwise. The step size $\eta=0.05$
in all the calculationss.

Figs. 2 and 3 show the performances of the proposed PDS($K$)-VQS ($K=2$, i.e.
PDS(2)-VQS) and the conventional VQE approaches for finding the ground state
energy of the toy Hamiltonians. As can be seen, the ability of VQE navigation
to avoid the local minima on the conventional PES depends on the Riemannian
metric exploited. For system A, in comparison to GD, the NGD (or equivalently
ITE in this case) is able to avoid the local minimum at
$(\theta_{1},\theta_{2})=(0,0)$. This is because the Riemannian metric,

$\displaystyle\mathcal{R}=\left(\begin{array}[]{cc}\frac{1}{4}&0\\\
0&\frac{1}{4}\sin^{2}(\frac{\theta_{1}}{2})\end{array}\right),$ (35)

used in the NGD/ITE correctly characterizes any rotation pair with
$\theta_{1}=0$ as a singular point (i.e. $\det|R|=0$) such that
$\mathcal{R}^{-1}$ will numerically navigate the dynamics (e.g. via singular
value decomposition) to avoid collapsing in this local minimum once the
trajectory is getting close. Therefore, if the metric is unable to
characterize the local minima as singular points, the VQE would still get
trapped. This can be observed from the VQE performance for system B, where
both NGD and ITE fail to escape the local minima,
$(\theta_{1},\theta_{2})\sim(\pm\frac{3\pi}{8},0)$, in the dynamics due to the
fact that the local minima are not the singular points of $R$ in either NGD or
ITE.

In contrast, the PDS(2)-VQS robustly converge to the true ground state for
both systems regardless of the employed Riemannian metric. The success of
PDS($K$)-VQS in these toy examples can be essentially attributed to the fact
that, in comparison to the original PES where the local minima corresponding
to a non-ground state, the entire PDS($K$) energy surface, except the singular
areas (see the infinitesimal white strips on the left panels of Fig. 2 at
$\theta_{1}=0$) where the fidelity of the trial wave function w.r.t the target
state is strictly zero, provides an upper bound energy surface for the same
(ground) state. This state-specific nature makes the PDS($K$)-VQS essentially
explore a lower upper bound of the ground state at a given PDS order, and
therefore the dynamics will not be trapped at a location that is associated
with a different state. It worth mentioning that, a lower bound of the ground
state energy can also be obtained from a static, and more costly, higher order
PDS($K$) standalone calculation as demonstrated in our previous work.[34] From
this perspective, the PDS($K$)-VQS approach provides an effective way to
explore the possibility of pushing the low order PDS($K$) results towards high
accuracy that would otherwise require higher-order and more expensive PDS($K$)
calculations. Besides, since generalized variational principle applies in the
PDS framework,[49, 62] if other roots of Eq. (12) are concerned, the
PDS($K$)-VQS will also be able to navigate the dynamics to give lower upper
bounds for excited states as long as the fidelity of the trial wave function
with respect to the target state is non-zero.

###  3.2 $H_{2}$ and $H_{4}$ systems

We further employ the proposed PDS($K$)-VQS approach to find the ground state
energy of H2 and H4 molecular systems. For H2 molecule, we exploit an
effective Hamiltonian and an ansatz exploited by Yamamoto [71] and Bravyi et
al.[9],

$\displaystyle H=0.4(\sigma_{z}\otimes
I+I\otimes\sigma_{z})+0.2\sigma_{x}\otimes\sigma_{x}$
$\displaystyle|\phi(\vec{\theta})\rangle=R_{Y}^{0}(2\theta_{3})R_{Y}^{1}(2\theta_{4})\tilde{U}_{N}^{0,1}R_{Y}^{0}(2\theta_{1})R_{Y}^{1}(2\theta_{2})|00\rangle$

where $\tilde{U}_{N}^{p,q}$ denotes the CNOT gate with control qubit $p$ and
target qubit $q$.

Figure 4: The computed ground state energy (top panels), energy deviation
w.r.t. exact energy (middle panels), and fidelity of the trial state (bottom
panels) of the H2 molecule iterate in the conventional VQE and PDS($K$)-VQS
($K=2,3,4$) infrastructures employing gradient descent (left panels) and
natural gradient descent/imaginary time evolution (right panels). The initial
rotation is given by $\vec{\theta}=(7\pi/32,\pi/2,0,0)$. The step size
$\eta=0.05$ in all the calculations.

Fig. 4 compares the VQE and PDS($K$)-VQS performances exploiting the above-
mentioned ansatz to find the ground state energy of the H2 Hamiltonian. As can
be seen, starting from the given initial rotation, the VQE is unable to
converge to the ground state energy within 100 iterations, but rather drops to
an excited state energy ($-0.2$ a.u. in this case). Actually, it has been
shown that,[71] starting from the same initial rotation, the VQE needs to go
through a “plateau” that resides at this energy value and spreads over
$\sim$400 iterations before hitting the ground state energy ($\sim$$-0.8$ a.u.
in this case) regardless of the employed Riemannian metric.

To achieve a higher level of accuracy (e.g. chemical accuracy
$\|E(\vec{\theta})-E_{exact})\|<1.5\times 10^{-3}$ a.u.), low order
PDS($K$)-VQS typically needs more iterations than high order PDS($K$)-VQS. As
shown in the middle left panel of Fig. 4, by employing GD in the dynamics, it
takes the PDS(4)-VQS $<$10 iterations to converge to the ground state energy
with energy deviation being $<10^{-14}$ a.u. regardless of the employed
Riemannian metric, while it takes the PDS(2)/PDS(3)-VQS almost 100 steps to
bound the deviation to be $<10^{-3}$ a.u. Remarkably, the performance can be
improved when GD is replaced by NGD/ITE in the PDS(2)/PDS(3)-VQS dynamics. In
particular, within 80 iterations the PDS(3)-VQS employing NGD/ITE can converge
to the accuracy level that is almost same as that of PDS(4)-VQS.

On the other hand, the quality of the trial wave function is more
significantly improved in the low order PDS($K$)-VQS dynamics than in the high
order PDS($K$)-VQS dynamics. For example, the fidelity of the trial wave
function w.r.t the exact ground state gradually increases from almost zero to
$\sim$0.35 within 100 iterations using PDS(2)-VQS regardless of the employed
Riemannian metric, and this change is significantly steeper than the almost
flat curves of PDS(3)/PDS(4)-VQS as shown at the bottom of Fig. 4. However, in
comparison to GD, employing NGD/ITE in the PDS(3)/PDS(4)-VQS quickly improves
the fidelity of the trial wave function from $<$0.02 to 0.2$\sim$0.3 within 10
iterations. It is worth mentioning that since the fidelity of the trial wave
function at the initial rotation is almost zero, both VQE and the static
PDS($K$) ($K=2,3,4$) calculations alone cannot help identify the ground state
energy in this case, which makes PDS($K$)-VQS a necessary and effective
approach to target ground state energy and improve the trial wave function.

Remarkably, the improvement of the trial wave function employing PDS($K$)-VQS
approach might be limited. This can be seen from the flat fidelity curves of
the trial state driven in the PDS(3/4)-VQS dynamics after first several
iterations as shown at the bottom right of Fig. 4. This is due to the fact
that the PDS($K$) formalism does not require the ansatz to sufficiently
approximate the target state, while is still able to provide systematically
improvable upper bounds of the expectation value of the target state by
exploring the Krylov subspace. The benefit is the great simplification of the
state preparation. The limitation is also obvious in that it sometimes would
be challenging to further improve the quality of the trial state within the
PDS($K$)-VQS framework if the energies were well converged already, which
would then compromise the accuracy of the property calculations that usually
requires a sufficiently accurate description of the target wave function.

Figure 5: The performance of VQE and PDS($K$)-VQS ($K=2,3,4$) employing
ordinary gradient descent (GD) to compute the ground state energy of a planar
H4 system with $R_{H-H}=2.0$ a.u. (Top left) The circuit used to generate the
ansatz with 16 rotation parameters that is inspired by the basis rotation
ansatz for a linear hydrogen chain in Ref. [2]. Here, we consider the planar
H4 system in 3-21G basis. The generated Hamiltonian acts on eight qubits, and
considers an active space of four electrons in eight spin-orbitals. With this
Hamiltonian, the ground state of the planar H4 system is a triplet state with
the exact energy $E_{\text{exact}}=-2.00591266$ a.u. (Bottom left) the
deviations of the VQE and PDS($K$)-VQS energies and (bottom right) the
fidelity change of the trial wave function w.r.t. true ground state during the
PDS($K$)-VQS calculations. The initial values of all the rotations are set to
0.001. The step size $\eta=1.0$ in all the calculations.

We also test the proposed PDS($K$)-VQS approach for a slightly larger system,
the planar H4 system, where a 8-qubit circuit with 16 rotation parameters
shown at the top of Fig. 5 is employed to prepare the trial wave function for
finding the ground state energy. The state preparation circuit is inspired by
the similar circuit that has been reported being successfully applied for
preparing the trial state for the Hartree Fock state of the  linear hydrogen
chain systems.[2] For the planar H4 system whose ground state is a triplet,
the circuit with close-to-zero initial rotations would generate a trial state
that is almost singlet, which makes the conventional VQE and the static
PDS($K$) ($K=2,3,4$) simply fail. On the other hand, as shown at the bottom of
Fig. 5, the PDS($K$)-VQS ($K=2,3,4$) are capable of dealing with such a tough
situation and again outperform. As can be seen, within 200 iterations,
PDS($K$)-VQS ($K=2,3,4$) are able to converge to the ground state energy well
below chemical accuracy and improving the fidelity of the trial wave function
to be $>$0.96. It is worth noting that even though the converged rotations
obtained from the PDS($K$)-VQS calculations generate a high fidelity state,
the expectation value of the generated state is still $\sim 0.02$ a.u. above
the exact energy, and it then becomes challenging to further improve the
fidelity employing the same circuit infrastructure through varying the
rotations. Therefore, the circuit used here might not be sufficient for
preparing true ground state in practice if higher fidelity is desired. We here
intentionally employ the circuit to artificially generate an extreme
challenging case to show the performance difference between conventional VQE
and PDS($K$)-VQS approaches.

##  4 Discussion

From Section III, it has been seen that the PDS($K$)-VQS approach bears the
potential of speeding up the iterations in comparison with the conventional
VQE approach. However, it is worth noting that the measurement effort of
evaluating $\langle H^{n}\rangle$’s ($n>1$) and their derivatives are usually
more expensive than that of $\langle H\rangle$ and its derivative, and the
actual cost saving will therefore be compromised.

To have a close look at the measurement of the $\langle H^{n}\rangle$ (and its
impact on the total cost), we employ the following metric to give an estimate
for the number of measurements, $M$,[23, 56, 69]

$\displaystyle M=\Bigg{(}\frac{\sum_{G}\sqrt{\sum_{i,j,\in
G}h_{i}h_{j}\text{cov}\big{(}P_{i},P_{j}\big{)}}}{\epsilon}\Bigg{)}^{2},$ (36)

where $\epsilon$ is the desired precision and $h_{i}$’s and $P_{i}$’s are the
coefficients and Pauli strings representing a moment (i.e.
$H^{n}=\sum_{i}h_{i}P_{i}$) and having been partitioned into certain groups,
$G$’s, in which simultaneous measurement can be performed.
$\text{cov}\big{(}P_{i},P_{j}\big{)}$ is the covariance between two Pauli
strings bounded by

$\displaystyle\text{cov}\big{(}P_{i},P_{j}\big{)}\leq|\sqrt{\text{var}(P_{i})\cdot\text{var}(P_{j})}|$
(37)

with the variance being computed from $\text{var}(P_{i})=1-\langle
P_{i}\rangle^{2}$. Here, we assume the covariances between different Pauli
strings to be zero for the brevity of the discussion. We can apply the above
metric to, for example, estimate the number of measurements of $H^{n}$
($n=1,2,3$) required by the PDS(2)-VQS calculation for the complete active
space (4 electrons, 4 spin-orbitals) of the planar H4 system. Given
$\epsilon\sim 0.5$mHartree, since $H^{n}$ ($n=1,2,3$) can be generated from at
most $\sim 3700$ Pauli strings, the estimated number of measurements needs to
be done is $\sim 4.8\times 10^{9}$, which is one order of magnitude higher
than that for $\langle H\rangle$ ($\sim 1.2\times 10^{8}$). Thus, given the
same trial state in this H4 case, if the number of conventional VQE iterations
is no more than one order of magnitude larger than that of the PDS(2)-VQS
iterations, VQE would outperform PDS(2)-VQS in terms of total number of
measurements, and PDS(2)-VQS outperforms otherwise. It is worth mentioning
that, during the PDS($K$)-VQS process for the ground state and energy, the
excited state energies can also be estimated directly from the higher roots of
the polynomial (12) without any additional measurement (although accurate
excited state energies would require higher order PDS($K$)-VQS calculations).
In contrast, the conventional VQE would need distinct trial states, and thus
different measurements, for targeting different states.

Generally speaking, as long as the relatively large number of measurements of
the Pauli strings becomes manageable, the PDS($K$)-VQS approach can be
potentially applied for targeting the exact solutions for the system sizes
that are not classically tractable, in particular for the systems whose true
ground and excited states we have little knowledge of, or are challenging to
obtain classically. To reduce the measurement demand, typical strategy is to
partition the Pauli strings (that contribute to the moments) into commuting
subsets that follow a certain rule, e.g. qubit-wise commutativity (QWC)[31,
44], general commutativity[22, 72], unitary partitioning[30], and/or Fermionic
basis rotation grouping[28] to name a few. The applications of these commuting
rules to the single Hamiltonian have shown that, at a cost of introducing
additional one-/multi-qubit unitary transformation before the measurement, the
total number of required measurements can be significantly reduced from
$\mathcal{O}(N^{4})$ to $\mathcal{O}(N^{2\sim 3})$, or for simpler cases even
$\mathcal{O}(N)$. For higher order moments, as we mentioned in the method
section, early study of applying QWC bases to Heisenberg models represented by
up to 36 qubits exhibits a sub-linear scaling of the number of measurements in
the number of the qubits (Ref. [67]), which then leads us to expect similar
scaling behaviors of the number of required measurements for evaluating
moments for molecular systems. Beside exploring the commutativity of Pauli
strings, other approaches including the linear combinations of unitary
operators (LCU) technique[13], direct block-encoding[21, 6], and quantum power
methods[59] might also be worth studying for reducing the number of
measurements at the cost of circuit depth. In the light of that, we plan to
perform a comprehensive benchmark as a follow-up work.

Since the PDS($K$)-VQS formalism involves solving linear system of equations
and polynomial, there is a concern of numerical instability when applying the
PDS($K$)-VQS approach in optimization. Theoretically, the numerical
instability of the PDS($K$)-VQS approach might come from two sources, (a) the
singularity and ill-conditioning of the matrix $M$ in Eq. (11) that might
consist of high order moments, and (b) the singularity of the Riemannian
metric ($\mathcal{R}$) used in the dynamics (16). In particular, the
singularity of matrix $M$ can be easily observed if the trial vector becomes
very close to the exact wave function ($\rm det|M|=0$ if we replace the trial
vector with exact vector). Numerically, the singularity problem can be avoided
by adding a small positive number (e.g. $10^{-6}$) to the eigenvalue of the
matrix $M$ or $\mathcal{R}$ via singular value decomposition (SVD). However,
it is worthing noting that adding small perturbation to $M$ might violate the
variationality of the PDS approach, and would not be recommended to use if the
strict upper-bounds to the true energy are concerned. The ill-conditioning of
matrix $M$ could occur in the high order PDS calculations, where high order
moments could make the condition number of matrix $M$ very large. Thus, from
the practical point of view, due to the potentially larger number of
measurements and ill-conditioning arising from high Hamiltonian powers, lower
oder PDS($K$)-VQS approaches are usually more feasible.

Figure 6: The performance of VQE and PDS(2)-VQS employing ordinary gradient
descent (GD) to compute the ground state energy of the four-site 2D Heisenberg
model. (Top right) The circuit employed to generate the trial vector, where
only the first rotation in RY gate is treated as a variational parameter
$\theta$, and other three rotations are fixed to $(0,3,3)$. Two initial
rotations $\theta_{0}=-2.0$ and $\theta_{0}=-3.0$ are chosen for performance
comparison. The exact ground state energy of the 2D Heisenberg model is
$E_{\text{exact}}=-3.6$ a.u. (Center) The VQE and PDS(2)-VQS energies and
(bottom) the corresponding fidelity changes of the trial vectors w.r.t. true
ground state in the first ten iterations in the conventional VQE and
PDS(2)-VQS noise-free calculations. The step size $\eta=1.0/\text{Iteration}$
in all the calculations. Figure 7: The computed ground state energy (top left)
and magnetization (top right) of the four-site 2D Heisenberg model and the
corresponding changes of the fidelity (bottom left) and variational parameter
$\theta$ (bottom right) in the first ten PDS($K$)-VQS ($K=2,3,4$) iterations
running on IBM Toronto quantum hardware. The physical setup, error sources,
and computed expectation values of Hamiltonian moments (up to $\langle
H^{7}\rangle$) and the associated standard deviations are shown in Fig. 8. In
all the calculations ordinary gradient descent (GD) is employed. The initial
rotation $\vec{\theta}_{0}=-3.0$. The exact ground state energy and
magnetization of the 2D Heisenberg model are $E_{\text{exact}}=-3.6$ a.u. and
$\sum_{i}\langle\sigma_{z_{i}}\rangle=-4.0$ a.u., respectively. The step size
$\eta=1.0/\text{Iteration}$ in all the calculations.

Ultimately, one would be concerned about how the PDS($K$)-VQS applies to
general models and how it performs on the real quantum hardware subject to the
device noise. To address these concerns and explore the potential of the
PDS($K$)-VQS approach, we have started to launch the PDS($K$)-VQS calculations
for more general Hamiltonians on both simulator and the real quantum hardware.
Figs. 6 and 7 exhibit some preliminary results for a four-site 2D Heisenberg
model with external magnetic field, $H=J\sum_{\langle
ij\rangle}\big{(}X_{i}X_{j}+Y_{i}Y_{j}+Z_{i}Z_{j}\big{)}+B\sum_{i}Z_{i}$ with
$J/B=0.1$. The simple circuit employed for the state preparation in both VQE
and PDS($K$)-VQS simulations is shown in Fig. 6, where, for the brevity of our
discussion, we only treat one rotation in the state preparation as the
variational parameter, and fix all other three rotations. As can be seen from
the noise-free simulations in Fig. 6, the PDS(2)-VQS results quickly converge
within five iterations achieving $\sim 0.99$ fidelity, while the performance
of VQE exhibits strong dependence on the initial rotation (for
$\vec{\theta}_{0}=-2.0$, the conventional VQE is able to converge in 10
iterations with $\Delta E<0.05$ a.u. and Fidelity $\sim 0.97$). When running
the PDS($K$)-VQS simulations for the same model on the IBM Toronto quantum
hardware, as shown in Fig. 7, in comparison to the ideal curves, the
PDS(2/3)-VQS optimization curves on the real hardware significantly slows
down, and deviate from the exact solutions due to the error from the real
machine. However, if we increase the PDS order to perform PDS(4)-VQS
calculations, the accuracy of the results systematically improves. For
example, in the PDS(4)-VQS approach both the computed ground state energy and
the trial state (and thus the magnetization) converge within 10 iterations
being very close to the exact solutions.

##  5 Conclusion

In summary, we propose a new variational quantum solver that employs the PDS
energy gradient. In comparison with the usual VQE, the PDS($K$)-VQS helps
identify an upper bound energy surface for the ground state, and thus frees
the dynamics from being trapped at local minima that refer to non-ground
states. In comparison with the static PDS($K$) expansions, the PDS($K$)-VQS
guides the rotation of the trial wave function of modest quality, and is able
to achieve high accuracy at the expense of low order PDS($K$) expansions. We
have demonstrated the capability of the PDS($K$)-VQS approach at finding the
ground state and its energy for toy models, H2 molecule, and strongly
correlated planar H4 system in some challenging situations. In all the case
studies, the PDS($K$)-VQS outperforms the standalone VQE and static PDS($K$)
calculations in terms of efficiency even at the lowest order. We also
discussed the limitations of the PDS($K$)-VQS approach at the current stage.
In particular, the PDS($K$)-VQS approach may suffer from large amount of
measurements for large systems, which can nevertheless be reduced at the cost
of circuit depth by working together with some measurement reduction methods.
Finally, we have started to launch PDS($K$)-VQS simulations for more general
Hamiltonians on IBM quantum hardware. Preliminary results for Heisenberg model
indicate that higher order PDS($K$)-VQS approach exhibits better noise-
resistance than the lower order ones. The discussed approach can be extended
to any variational formulation based on the utilization of $\langle
H^{n}\rangle$ moments (e.g. Krylov subspace algorithms).

Figure 8: (a) Quantum processor device map for ibmq_toronto showing the four
qubits (Qn, $n=0-3$) used in the present computation. (b) Average CNOT error,
1-qubit readout assignment error, and thermal relaxation time constant (T1)
and dephasing time constant (T2) in the four qubits used in the present
computation. (c) The expectation values of the Hamiltonian moments, $\langle
H^{n}\rangle$ ($n=1-7$), assembled from the measurements of the expectation
values of 21 QWC bases for four-site 2D Heisenberg model $H=J\sum_{\langle
ij\rangle}\big{(}X_{i}X_{j}+Y_{i}Y_{j}+Z_{i}Z_{j}\big{)}+B\sum_{i}Z_{i}$ with
$J/B=0.1$. The data points correspond to mean value from the calculations on
IBM Quantum processor ibmq_toronto with statistical error bars corresponding
to $5\times 8192$ shots (per point). The trial state is constructed using the
circuit given in Fig. 6 with initial rotation $\theta_{0}=-3.0$.

##  6 Acknowledgement

B. P. and K. K. were 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. B. P. and K. K. acknowledge the use of the IBMQ for this work.
The views expressed are those of the authors and do not reflect the official
policy or position of IBM or the IBMQ team.

##  7 Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Amari [1998] S. Amari.  Natural gradient works efficiently in learning.  _Neural Computation_ , 10(2):251–276, 1998.  doi: 10.1162/089976698300017746.
  * Arute et al. [2020] Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph C. Bardin, Rami Barends, Sergio Boixo, Michael Broughton, Bob B. Buckley, David A. Buell, Brian Burkett, Nicholas Bushnell, Yu Chen, Zijun Chen, Benjamin Chiaro, Roberto Collins, William Courtney, Sean Demura, Andrew Dunsworth, Edward Farhi, Austin Fowler, Brooks Foxen, Craig Gidney, Marissa Giustina, Rob Graff, Steve Habegger, Matthew P. Harrigan, Alan Ho, Sabrina Hong, Trent Huang, William J. Huggins, Lev Ioffe, Sergei V. Isakov, Evan Jeffrey, Zhang Jiang, Cody Jones, Dvir Kafri, Kostyantyn Kechedzhi, Julian Kelly, Seon Kim, Paul V. Klimov, Alexander Korotkov, Fedor Kostritsa, David Landhuis, Pavel Laptev, Mike Lindmark, Erik Lucero, Orion Martin, John M. Martinis, Jarrod R. McClean, Matt McEwen, Anthony Megrant, Xiao Mi, Masoud Mohseni, Wojciech Mruczkiewicz, Josh Mutus, Ofer Naaman, Matthew Neeley, Charles Neill, Hartmut Neven, Murphy Yuezhen Niu, Thomas E. O’Brien, Eric Ostby, Andre Petukhov, Harald Putterman, Chris Quintana, Pedram Roushan, Nicholas C. Rubin, Daniel Sank, Kevin J. Satzinger, Vadim Smelyanskiy, Doug Strain, Kevin J. Sung, Marco Szalay, Tyler Y. Takeshita, Amit Vainsencher, Theodore White, Nathan Wiebe, Z. Jamie Yao, Ping Yeh, and Adam Zalcman.  Hartree-fock on a superconducting qubit quantum computer.  _Science_ , 369(6507):1084–1089, 2020.  doi: 10.1126/science.abb9811.
  * Babbush et al. [2018] Ryan Babbush, Nathan Wiebe, Jarrod McClean, James McClain, Hartmut Neven, and Garnet Kin-Lic Chan.  Low-depth quantum simulation of materials.  _Phys. Rev. X_ , 8:011044, 2018.  doi: 10.1103/PhysRevX.8.011044.
  * Bartlett et al. [1989] Rodney J. Bartlett, Stanisław A. Kucharski, and Jozef Noga.  Alternative coupled-cluster ansätze II. the unitary coupled-cluster method.  _Chem. Phys. Lett._ , 155(1):133–140, 1989.  doi: 10.1016/s0009-2614(89)87372-5.
  * Berry et al. [2007] Dominic W Berry, Graeme Ahokas, Richard Cleve, and Barry C Sanders.  Efficient quantum algorithms for simulating sparse hamiltonians.  _Comm. Math. Phys._ , 270(2):359–371, 2007.  doi: 10.1007/s00220-006-0150-x.
  * Berry et al. [2015] Dominic W. Berry, Andrew M. Childs, Richard Cleve, Robin Kothari, and Rolando D. Somma.  Simulating hamiltonian dynamics with a truncated taylor series.  _Phys. Rev. Lett._ , 114:090502, 2015.  doi: 10.1103/PhysRevLett.114.090502.
  * Bespalova and Kyriienko [2020] Tatiana A. Bespalova and Oleksandr Kyriienko.  Hamiltonian operator approximation for energy measurement and ground state preparation.  _preprint_ , arXiv:2009.03351, 2020.  URL https://arxiv.org/abs/2009.03351.
  * Bogolubov [1947] N. Bogolubov.  On the theory of superfluidity.  _J. Phys._ , 11:23–32, 1947.  doi: 10.1142/9789814612524_0001.
  * Bravyi et al. [2017] Sergey Bravyi, Jay M. Gambetta, Antonio Mezzacapo, and Kristan Temme.  Tapering off qubits to simulate fermionic hamiltonians.  _preprint_ , arXiv:1701.08213, 2017.  URL https://arxiv.org/abs/1701.08213.
  * Cerezo and Coles [2021] M Cerezo and Patrick J Coles.  Higher order derivatives of quantum neural networks with barren plateaus.  _Quantum Sci. Technol._ , 6(3):035006, jun 2021\.  doi: 10.1088/2058-9565/abf51a.
  * Cerezo et al. [2021] M. Cerezo, Akira Sone, Tyler Volkoff, Lukasz Cincio, and Patrick J. Coles.  Cost function dependent barren plateaus in shallow parametrized quantum circuits.  _Nat. Commun._ , 12(1), 2021.  doi: 10.1038/s41467-021-21728-w.
  * Childs [2010] Andrew M Childs.  On the relationship between continuous-and discrete-time quantum walk.  _Comm. Math. Phys._ , 294(2):581–603, 2010.  doi: 10.1007/s00220-009-0930-1.
  * Childs and Wiebe [2012] Andrew M. Childs and Nathan Wiebe.  Hamiltonian simulation using linear combinations of unitary operations.  _Quantum Inf. Comput._ , 12(11-12):0901–0924, 2012.  doi: 10.26421/qic12.11-12.
  * Claudino et al. [2021] Daniel Claudino, Bo Peng, Nicholas P. Bauman, Karol Kowalski, and Travis S. Humble.  Improving the accuracy and efficiency of quantum connected moments expansions.  _Quantum Sci. Technol._ , accepted, 2021.  doi: 10.1088/2058-9565/ac0292.
  * Cleve et al. [1998] Richard Cleve, Artur Ekert, Chiara Macchiavello, and Michele Mosca.  On quantum algorithms.  _Proc. R. Soc. Lond. A_ , 454(1969):339–354, 1998\.  doi: 10.1002/(SICI)1099-0526(199809/10)4:1<33::AID-CPLX10>3.0.CO;2-U.
  * 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.  doi: 10.1103/PhysRevX.8.011021.
  * Evangelista et al. [2019] Francesco A Evangelista, Garnet Kin-Lic Chan, and Gustavo E Scuseria.  Exact parameterization of fermionic wave functions via unitary coupled cluster theory.  _J. Chem. Phys._ , 151(24):244112, 2019.  doi: 10.1063/1.5133059.
  * Fessatidis et al. [2006] Vassilios Fessatidis, Jay D Mancini, Robert Murawski, and Samuel P Bowen.  A generalized moments expansion.  _Phys. Lett. A_ , 349(5):320–323, 2006.  doi: 10.1016/j.physleta.2005.09.039.
  * Fessatidis et al. [2010] Vassilios Fessatidis, Frank A Corvino, Jay D Mancini, Robert K Murawski, and John Mikalopas.  Analytic properties of moments matrices.  _Phys. Lett. A_ , 374(28):2890–2893, 2010.  doi: 10.1016/j.physleta.2010.05.010.
  * Feynman [1955] R. P. Feynman.  Slow electrons in a polar crystal.  _Phys. Rev._ , 97:660–665, 1955.  doi: 10.1103/PhysRev.97.660.
  * Gilyén et al. [2019] András Gilyén, Yuan Su, Guang Hao Low, and Nathan 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, page 193–204, New York, NY, USA, 2019. Association for Computing Machinery.  ISBN 9781450367059.  doi: 10.1145/3313276.3316366.
  * Gokhale et al. [2020] Pranav Gokhale, Olivia Angiuli, Yongshan Ding, Kaiwen Gui, Teague Tomesh, Martin Suchara, Margaret Martonosi, and Frederic T. Chong.  $o(n^{3})$ measurement cost for variational quantum eigensolver on molecular hamiltonians.  _IEEE Trans. Qunatum Eng._ , 1:1–24, 2020.  doi: 10.1109/TQE.2020.3035814.
  * Gonthier et al. [2020] Jérôme F. Gonthier, Maxwell D. Radin, Corneliu Buda, Eric J. Doskocil, Clena M. Abuan, and Jhonathan Romero.  Identifying challenges towards practical quantum advantage through resource estimation: the measurement roadblock in the variational quantum eigensolver.  _preprint_ , arXiv:2012.04001, 2020.  URL https://arxiv.org/abs/2012.04001.
  * Grimsley et al. [2019] Harper R Grimsley, Sophia E Economou, Edwin Barnes, and Nicholas J Mayhall.  An adaptive variational algorithm for exact molecular simulations on a quantum computer.  _Nat. Commun._ , 10(1):1–9, 2019.  doi: 10.1038/s41467-019-10988-2.
  * Guerreschi and Smelyanskiy [2017] Gian Giacomo Guerreschi and Mikhail Smelyanskiy.  Practical optimization for hybrid quantum-classical algorithms.  _preprint_ , arXiv:1701.01450, 2017.  URL https://arxiv.org/abs/1701.01450.
  * Häner et al. [2016] T. Häner, D. S. Steiger, M. Smelyanskiy, and M. Troyer.  High performance emulation of quantum circuits.  In _SC ’16: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_ , pages 866–874, Nov 2016\.  doi: 10.1109/SC.2016.73.
  * Hoffmann and Simons [1988] Mark R Hoffmann and Jack Simons.  A unitary multiconfigurational coupled-cluster method: Theory and applications.  _J. Chem. Phys._ , 88(2):993–1002, 1988.  doi: 10.1063/1.454125.
  * Huggins et al. [2021] William J. Huggins, Jarrod R. McClean, Nicholas C. Rubin, Zhang Jiang, Nathan Wiebe, K. Birgitta Whaley, and Ryan Babbush.  Efficient and noise resilient measurements for quantum chemistry on near-term quantum computers.  _npj Quantum Inf._ , 7(1):23, 2021.  doi: 10.1038/s41534-020-00341-7.
  * Huggins et al. [2020] William James Huggins, Joonho Lee, Unpil Baek, Bryan O’Gorman, and K Birgitta Whaley.  A non-orthogonal variational quantum eigensolver.  _New J. Phys._ , 22:073009, 2020.  doi: 10.1088/1367-2630/ab867b.
  * Izmaylov et al. [2020] Artur F. Izmaylov, Tzu-Ching Yen, Robert A. Lang, and Vladyslav Verteletskyi.  Unitary partitioning approach to the measurement problem in the variational quantum eigensolver method.  _J. Chem. Theory Comput._ , 16(1):190–195, 2020\.  doi: 10.1021/acs.jctc.9b00791.
  * Kandala et al. [2017] Abhinav Kandala, Antonio Mezzacapo, Kristan Temme, Maika Takita, Markus Brink, Jerry M. Chow, and Jay M. Gambetta.  Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets.  _Nature_ , 549:242–246, 2017.  doi: 10.1038/nature23879.
  * Kandala et al. [2019] Abhinav Kandala, Kristan Temme, Antonio D Corcoles, Antonio Mezzacapo, Jerry M Chow, and Jay M Gambetta.  Error mitigation extends the computational reach of a noisy quantum processor.  _Nature_ , 567:491–495, 2019.  doi: 10.1038/s41586-019-1040-7.
  * Knowles [1987] Peter J Knowles.  On the validity and applicability of the connected moments expansion.  _Chem. Phys. Lett._ , 134(6):512–518, 1987.  doi: 10.1016/0009-2614(87)87184-1.
  * Kowalski and Peng [2020] Karol Kowalski and Bo Peng.  Quantum simulations employing connected moments expansions.  _J. Chem. Phys._ , 153(20):201102, 2020.  doi: 10.1063/5.0030688.
  * Kutzelnigg [1991] Werner Kutzelnigg.  Error analysis and improvements of coupled-cluster theory.  _Theor. Chim. Acta_ , 80(4-5):349–386, 1991.  doi: 10.1007/BF01117418.
  * Kyriienko [2020] Oleksandr Kyriienko.  Quantum inverse iteration algorithm for programmable quantum simulators.  _npj Quantum Inf._ , 6(1):1–8, 2020.  doi: 10.1038/s41534-019-0239-7.
  * Lee et al. [2019] Joonho Lee, William J. Huggins, Martin Head-Gordon, and K. Birgitta Whaley.  Generalized unitary coupled cluster wave functions for quantum computation.  _J. Chem. Theory Comput._ , 15(1):311–324, 2019\.  doi: 10.1021/acs.jctc.8b01004.
  * Luis and Peřina [1996] A Luis and J Peřina.  Optimum phase-shift estimation and the quantum description of the phase difference.  _Phys. Rev. A_ , 54(5):4564, 1996.  doi: 10.1103/PhysRevA.54.4564.
  * Mancini et al. [1994] Jay D Mancini, Yu Zhou, and Peter F Meier.  Analytic properties of connected moments expansions.  _Int. J. Quantum Chem._ , 50(2):101–107, 1994\.  doi: 10.1002/qua.560500203.
  * Mancini et al. [1995] Jay D Mancini, William J Massano, Janice D Prie, and Yu Zhuo.  Avoidance of singularities in moments expansions: a numerical study.  _Phys. Lett. A_ , 209(1-2):107–112, 1995.  doi: 10.1016/0375-9601(95)00757-2.
  * Marrero et al. [2020] Carlos Ortiz Marrero, Mária Kieferová, and Nathan Wiebe.  Entanglement induced barren plateaus.  _preprint_ , arXiv:2010.15968, 2020.  URL https://arxiv.org/abs/2010.15968.
  * McArdle et al. [2019] Sam McArdle, Tyson Jones, Suguru Endo, Ying Li, Simon C Benjamin, and Xiao Yuan.  Variational ansatz-based quantum simulation of imaginary time evolution.  _npj Quantum Inf._ , 5(1):1–6, 2019.  doi: 10.1038/s41534-019-0187-2.
  * McArdle et al. [2020] Sam McArdle, Suguru Endo, Alan Aspuru-Guzik, Simon C Benjamin, and Xiao Yuan.  Quantum computational chemistry.  _Rev. Mod. Phys._ , 92(1):015003, 2020.  doi: 10.1103/RevModPhys.92.015003.
  * McClean et al. [2016] Jarrod R McClean, Jonathan Romero, Ryan Babbush, and Alán Aspuru-Guzik.  The theory of variational hybrid quantum-classical algorithms.  _New J. Phys._ , 18(2):023023, 2016.  doi: 10.1088/1367-2630/18/2/023023.
  * McClean et al. [2018] Jarrod R McClean, Sergio Boixo, Vadim N Smelyanskiy, Ryan Babbush, and Hartmut Neven.  Barren plateaus in quantum neural network training landscapes.  _Nat. Commun._ , 9(1):1–6, 2018.  doi: 10.1038/s41467-018-07090-4.
  * Motta et al. [2020] Mario Motta, Chong Sun, Adrian TK Tan, Matthew J O’Rourke, Erika Ye, Austin J Minnich, Fernando GSL Brandão, and Garnet Kin-Lic Chan.  Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution.  _Nat. Phys._ , 16(2):205–210, 2020.  doi: 10.1038/s41567-019-0704-4.
  * Nielsen and Chuang [2011] Michael A. Nielsen and Isaac L. Chuang.  _Quantum Computation and Quantum Information: 10th Anniversary Edition_.  Cambridge University Press, New York, NY, USA, 10th edition, 2011.  ISBN 1107002176, 9781107002173.  doi: 10.1017/CBO9780511976667.
  * Parrish and McMahon [2019] Robert M Parrish and Peter L McMahon.  Quantum filter diagonalization: Quantum eigendecomposition without full quantum phase estimation.  _preprint_ , arXiv:1909.08925, 2019.  URL https://arxiv.org/abs/1909.08925.
  * Peeters and Devreese [1984] François M Peeters and Jozef T Devreese.  Upper bounds for the free energy. a generalisation of the bogolubov inequality and the feynman inequality.  _J. Phys. A: Math. Gen._ , 17(3):625, 1984.  doi: 10.1088/0305-4470/17/3/024.
  * Peruzzo et al. [2014] Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J Love, Alán Aspuru-Guzik, and Jeremy L O’brien.  A variational eigenvalue solver on a photonic quantum processor.  _Nat. Commun._ , 5:4213, 2014.  doi: 10.1038/ncomms5213.
  * Pesah et al. [2020] Arthur Pesah, M Cerezo, Samson Wang, Tyler Volkoff, Andrew T Sornborger, and Patrick J Coles.  Absence of barren plateaus in quantum convolutional neural networks.  _preprint_ , arXiv:2011.02966, 2020.  URL https://arxiv.org/abs/2011.02966.
  * Poulin et al. [2018] David Poulin, Alexei Kitaev, Damian S. Steiger, Matthew B. Hastings, and Matthias Troyer.  Quantum algorithm for spectral measurement with a lower gate count.  _Phys. Rev. Lett._ , 121:010501, 2018.  doi: 10.1103/PhysRevLett.121.010501.
  * Preskill [2018] John Preskill.  Quantum computing in the nisq era and beyond.  _Quantum_ , 2:79, 2018.  doi: 10.22331/q-2018-08-06-79.
  * Prie et al. [1994] Janice D Prie, D Schwall, Jay D Mancini, D Kraus, and William J Massano.  On the relation between the connected-moments expansion and the lanczos variational scheme.  _Nuov. Cim. D_ , 16(5):433–448, 1994.  doi: 10.1007/BF02463732.
  * Romero et al. [2018] Jonathan Romero, Ryan Babbush, Jarrod R McClean, Cornelius Hempel, Peter J Love, and Alán Aspuru-Guzik.  Strategies for quantum computing molecular energies using the unitary coupled cluster ansatz.  _Quantum Sci. Technol._ , 4(1):014008, 2018.  doi: 10.1088/2058-9565/aad3e4.
  * Rubin et al. [2018] Nicholas C Rubin, Ryan Babbush, and McClean Jarrod.  Application of fermionic marginal constraints to hybrid quantum algorithms.  _New J. Phys._ , 20:053020, 2018.  doi: 10.1088/1367-2630/aab919.
  * Schuld et al. [2019] Maria Schuld, Ville Bergholm, Christian Gogolin, Josh Izaac, and Nathan Killoran.  Evaluating analytic gradients on quantum hardware.  _Phys. Rev. A_ , 99:032331, 2019.  doi: 10.1103/PhysRevA.99.032331.
  * Seeley et al. [2012] Jacob T. Seeley, Martin J. Richard, and Peter J. Love.  The bravyi-kitaev transformation for quantum computation of electronic structure.  _J. Chem. Phys._ , 137(22):224109, 2012.  doi: 10.1063/1.4768229.
  * Seki and Yunoki [2021] Kazuhiro Seki and Seiji Yunoki.  Quantum power method by a superposition of time-evolved states.  _PRX Quantum_ , 2:010333, 2021.  doi: 10.1103/PRXQuantum.2.010333.
  * Shen et al. [2017] Yangchao Shen, Xiang Zhang, Shuaining Zhang, Jing-Ning Zhang, Man-Hong Yung, and Kihwan Kim.  Quantum implementation of the unitary coupled cluster for simulating molecular electronic structure.  _Phys. Rev. A_ , 95:020501, 2017.  doi: 10.1103/PhysRevA.95.020501.
  * Shor [1999] Peter W Shor.  Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer.  _SIAM Rev._ , 41(2):303–332, 1999.  doi: 10.1137/S0036144598347011.
  * Soldatov [1995] Andrey V Soldatov.  Generalized variational principle in quantum mechanics.  _Int. J. Mod. Phys. B_ , 9(22):2899–2936, 1995\.  doi: 10.1142/S0217979295001087.
  * Stokes et al. [2020] James Stokes, Josh Izaac, Nathan Killoran, and Giuseppe Carleo.  Quantum Natural Gradient.  _Quantum_ , 4:269, 2020.  doi: 10.22331/q-2020-05-25-269.
  * Taube and Bartlett [2006] Andrew G Taube and Rodney J Bartlett.  New perspectives on unitary coupled-cluster theory.  _Int. J. Quantum Chem._ , 106(15):3393–3401, 2006\.  doi: 10.1002/qua.21198.
  * Ullah [1995] Nazakat Ullah.  Removal of the singularity in the moment-expansion formalism.  _Phys. Rev. A_ , 51(3):1808, 1995.  doi: 10.1103/PhysRevA.51.1808.
  * Uvarov and Biamonte [2021] Alexey Uvarov and Jacob Biamonte.  On barren plateaus and cost function locality in variational quantum algorithms.  _J. Phys. A: Math. and Theo._ , 54:245301, 2021.  doi: 10.1088/1751-8121/abfac7.
  * Vallury et al. [2020] Harish J. Vallury, Michael A. Jones, Charles D. Hill, and Lloyd C. L. Hollenberg.  Quantum computed moments correction to variational estimates.  _Quantum_ , 4:373, 2020.  doi: 10.22331/q-2020-12-15-373.
  * Wang et al. [2020] Samson Wang, Enrico Fontana, Marco Cerezo, Kunal Sharma, Akira Sone, Lukasz Cincio, and Patrick J Coles.  Noise-induced barren plateaus in variational quantum algorithms.  _preprint_ , arXiv:2007.14384, 2020.  URL https://arxiv.org/abs/2007.14384.
  * Wecker et al. [2015] Dave Wecker, Matthew B. Hastings, and Matthias Troyer.  Progress towards practical quantum variational algorithms.  _Phys. Rev. A_ , 92:042303, 2015.  doi: 10.1103/PhysRevA.92.042303.
  * Wierichs et al. [2020] David Wierichs, Christian Gogolin, and Michael Kastoryano.  Avoiding local minima in variational quantum eigensolvers with the natural gradient optimizer.  _Phys. Rev. Research_ , 2:043246, 2020.  doi: 10.1103/PhysRevResearch.2.043246.
  * Yamamoto [2019] Naoki Yamamoto.  On the natural gradient for variational quantum eigensolver.  _preprint_ , arXiv:1909.05074, 2019.  URL https://arxiv.org/abs/1909.05074.
  * Yen et al. [2020] Tzu-Ching Yen, Vladyslav Verteletskyi, and Artur F. Izmaylov.  Measuring all compatible operators in one series of single-qubit measurements using unitary transformations.  _J. Chem. Theory Comput._ , 16(4):2400–2409, 2020\.  doi: 10.1021/acs.jctc.0c00008.
  * Yuan et al. [2019] Xiao Yuan, Suguru Endo, Qi Zhao, Ying Li, and Simon C. Benjamin.  Theory of variational quantum simulation.  _Quantum_ , 3:191, 2019.  doi: 10.22331/q-2019-10-07-191.