# Unleashed from Constrained Optimization: Quantum Computing for Quantum
Chemistry Employing Generator Coordinate Method

Muqing Zheng  Lehigh University, Bethlehem, Pennsylvania, 18015, USA  Pacific
Northwest National Laboratory, Richland, Washington, 99354, USA  Bo Peng
<EMAIL_ADDRESS>Pacific Northwest National Laboratory, Richland, Washington,
99354, USA  Ang Li  Pacific Northwest National Laboratory, Richland,
Washington, 99354, USA  Xiu Yang  Lehigh University, Bethlehem, Pennsylvania,
18015, USA  Karol Kowalski  Pacific Northwest National Laboratory, Richland,
Washington, 99354, USA

(December 2023)

Highlights

  * •

We categorize the challenges in hybrid quantum-classical approaches, such as
the barren plateau and the exactness of the ansätze, as constrained
optimization problems.

  * •

We demonstrate the interconnection between constrained optimization and the
generalized eigenvalue problem in quantum chemistry computations. This is
achieved using a distinctive class of non-orthogonal and overcomplete basis
sets generated by Givens rotation-type canonical transformations on a
reference state.

  * •

Using the generator coordinate approach, we represent the wave function with
these basis sets. The resulting generalized eigenvalue problem establishes
rigorous lower bounds on energy, surpassing the performance of the
conventional variational quantum eigensolver (VQE) that utilizes the same
canonical transformations in its ansätz.

  * •

For practical applications, we introduce an adaptive scheme to select these
transformations. This strategy features the linear expansion of the non-
orthogonal basis set and bypasses the optimization, ensuring a balanced blend
of accuracy and efficiency in hybrid quantum-classical simulations.

  * •

Our methodology paves the way for broader applications in quantum chemistry,
which encompasses simultaneous ground and excited state simulations and the
generation of tailored projected/effective Hamiltonians for large scale
simulations.

Abstract Hybrid quantum-classical approaches offer potential solutions for
quantum chemistry problems, but they also introduce challenges such as the
barren plateau and the exactness of the ansätze. These challenges often
manifest as constrained optimization problems without a guarantee of
identifying global minima. In this work, we highlight the interconnection
between constrained optimization and generalized eigenvalue problems, using a
unique class of non-orthogonal and overcomplete basis sets generated by Givens
rotation-type canonical transformations on a reference state. Employing the
generator coordinate approach, we represent the wave function in terms of
these basis sets. The ensuing generalized eigenvalue problem yields rigorous
lower bounds on energy, outperforming the conventional variational quantum
eigensolver (VQE) that employs the same canonical transformations in its
ansätze. Our approach effectively tackles the barren plateau issue and the
heuristic nature of numerical minimizers in the standard VQE, making it ideal
for intricate quantum chemical challenges. For real-world applications, we
propose an adaptive scheme for selecting these transformations, emphasizing
the linear expansion of the non-orthogonal basis sets. This ensures a
harmonious balance between accuracy and efficiency in hybrid quantum-classical
simulations. Our analysis and suggested methodology further broaden the
applications of quantum computing in quantum chemistry. Notably, they pave the
way for alternative strategies in excited state computation and Hamiltonian
downfolding, laying the groundwork for sophisticated quantum simulations in
chemistry.

##  I Introduction

Accurately obtaining ground and excited state energies, along with the
corresponding many-body wave functions, is pivotal in comprehending diverse
physical phenomena in molecules and materials. This ranges from high-
temperature superconductivity in materials like cuprates Imada  _et al._
(1998) and bond-breaking chemical reactions to complex electronic processes in
biological and synthetic catalysts with transition metals Witzke  _et al._
(2020) or $f$-block atoms Gould  _et al._ (2022). The associated spin,
electronic properties, and dynamics are crucial for deciphering the structure-
property-function correlation in various fields, including catalysis, sensors,
and quantum materials. A prime example is the oxygen-evolving complex in
photosystem II of green plants, essential for water molecule oxidation and
photosynthesis Askerka  _et al._ (2017); Raymond and Blankenship (2008). This
complex, involving manganese atoms in various oxidation states, facilitates
electron-transfer reactions for water oxidation. However, this task becomes
exceptionally challenging in the presence of non-trivial quantum effects, such
as strong electron correlation, which influence the evolution of nuclei,
electrons, and spins under external stimuli. Conventional wave function
methodologies, like configuration interaction, coupled cluster, and many-body
perturbation theory, are tailored for diverse electron correlation scenarios
Shavitt and Bartlett (2009). Still, they often fall short in handling complex
cases or exhibit prohibitive scaling with increasing system size.
Consequently, this has become a vigorous area of computational research,
encompassing both classical and burgeoning quantum computing studies. The
primary objective is to strike an optimal balance between accuracy and
computational scalability.

In the realm of quantum computing, significant strides have been made towards
promising near-term hybrid quantum-classical strategies, which includes
variational quantum algorithm 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); Cao  _et
al._ (2019); Grimsley  _et al._ (2019a, b); Verteletskyi  _et al._ (2020);
McArdle  _et al._ (2019, 2020); Tilly  _et al._ (2021), quantum approximate
optimization algorithm Farhi  _et al._ (2014); Stein  _et al._ (2022), quantum
annealing Bharti  _et al._ (2021); Albash and Lidar (2018), Gaussian boson
sampling Aaronson and Arkhipov (2011), analog quantum simulation Trabesinger
(2012); Georgescu  _et al._ (2014), iterative quantum assisted eigensolver
McArdle  _et al._ (2019); Motta  _et al._ (2020); Parrish and McMahon (2019);
Kyriienko (2020), and many others. These approaches typically delegate certain
computational tasks to classical computers, thereby conserving quantum
resources in contrast to exclusively to quantum methods. Within this
framework, the variational quantum eigensolver (VQE) and its adaptive
derivatives are seen as the frontrunners in leveraging near-future quantum
advantages Peruzzo  _et al._ (2014); Grimsley  _et al._ (2019a). However,
these anticipations are also potentially impeded by the heuristic nature
inherent in the critical optimization processes. Issues such as the rigor of
ansätze exactness, the challenge in navigating potential energy surfaces
replete with numerous local minima, and the numerical nuances in minimization
techniques, remain nebulous. These queries are further convoluted when
considering the scalability of these methods concerning the number of
operators or the depth of quantum circuits involved.

As an alternative to VQE, which employs a highly nonlinear parametrization of
the wave function, other near-term strategies aim to construct and traverse a
subspace within the Hilbert space to closely approximate the desired state and
energy. Typical examples include the quantum subspace expansion Colless  _et
al._ (2018); McClean  _et al._ (2017); Takeshita  _et al._ (2020); McClean
_et al._ (2020), the hybrid and quantum Lanczos approaches Suchsland  _et al._
(2021); Motta  _et al._ (2020); Tkachenko  _et al._ (2023), quantum computed
moment approaches Seki and Yunoki (2021); Kyriienko (2020); Kowalski and Peng
(2020); Vallury  _et al._ (2020); Aulicino  _et al._ (2022), and quantum
equation-of-motion approach Ollitrault  _et al._ (2020). These strategies
often draw inspiration from truncated configuration interaction approaches or
explore the Krylov subspace.However, the comparison and interrelation of these
methods, especially when contrasting subspace expansion with nonlinear
optimization, frequently remain unclear. This uncertainty is exacerbated by
the choice of ansätze, which rarely guarantees exactness.

Recently, we employed the generator coordinate method (GCM) Hill and Wheeler
(1953); Griffin and Wheeler (1957); Rodríguez-Guzmán  _et al._ (2002); Bender
_et al._ (2003); Ring and Schuck (2004); Yao  _et al._ (2010); Egido (2016);
Hizawa  _et al._ (2021) and nonorthogonal quantum eigensolver as an
alternative near-term approach Zheng  _et al._ (2023), potentially addressing
the limitations identified in the VQE. Similar to its contemporaries, the GCM
employs low-depth quantum circuits and efficiently utilizes existing ansätze
to explore a subspace, targeting specific states and energy levels.
Nonetheless, the original GCM demands a priori knowledge of the system to
meticulously select both the ansätze and generator, thereby circumventing
heuristic approaches but potentially affecting its scalability and efficiency.

In this study, we discerned that the Unitary Coupled Cluster (UCC) excitation
generator, primarily implemented through Givens rotations, provides a
framework for delineating the differences and connections between the VQE and
GCM. This understanding illuminated why, given a consistent ansätze, solving a
generalized eigenvalue problem within the GCM framework could establish a
lower boundary compared to the customarily constrained optimization problem
encountered in the VQE method. Equipped with this knowledge, we implemented a
disentangled UCC ansätze and an adaptive generator selection scheme, aiming to
overcome the challenges present in the original GCM. Our results suggest that
the new ADAPT-GCM, guided by the subspace expansion using the employed UCC
disentangled ans”atze, outperforms the ADAPT-VQE method. This superiority is
evidenced by numerical results across various benchmark systems, particularly
in instances of strong electron correlation.

##  II GCM: A Brief Overview

Originating from studies of nuclear structure and nuclear reactions, including
those occurring in heavy-ion collisions or astrophysical scenarios, the
Generator Coordinate Method (GCM), proposed by Griffin-Hill-Wheeler (GHW) Hill
and Wheeler (1953); Griffin and Wheeler (1957), approximates the many-body
wave function $|\Psi\rangle$ of the Hamiltonian $H$ by describing collective
motion through the variation of one or more collective generator coordinates,
denoted as ${\bm{\alpha}}=\\{\alpha_{i}\\}$, where $i=1,\cdots,N_{\alpha}$ and
$N_{\alpha}>1$. This variation leads to the representation of the many-body
wave function of a quantum system in terms of a set of basis states (also
called generating functions) $\\{|\phi(\bm{\alpha})\rangle\\}$:

$\displaystyle|\Psi\rangle\approx|\psi_{\rm GCM}\rangle=\int{\rm
d}{\bm{\alpha}}|\phi({\bm{\alpha}})\rangle f({\bm{\alpha}}).$ (1)

Here, $f({\bm{\alpha}})$ is an unknown ${\bm{\alpha}}$-dependent weight
function. The choice of these generator coordinates ${\bm{\alpha}}$ depends on
the specific system under study. The fundamental concept of the GCM is to
treat these generator coordinates as variational parameters, which are then
used to approximate the target wave function as a superposition of
$\\{|\phi({\bm{\alpha}})\rangle\\}$. The GCM’s ability to adapt to changing
collective coordinates and its variational nature make it well-suited for
describing complex states that break symmetries, allowing for a flexible and
efficient description of properties such as ground states, excited states, and
transitions between them.

Technically, to determine $f({\bm{\alpha}})$, substituting (1) into the
Schrödinger equation, and then the variational process leads to the so-called
GHW integral equation:

$\displaystyle\int{\rm
d}{\bm{\alpha}}^{\prime}\bigg{(}\langle\phi(\bm{\alpha})|H|\phi(\bm{\alpha}^{\prime})\rangle-E\langle\phi(\bm{\alpha})|\phi(\bm{\alpha}^{\prime})\rangle\bigg{)}f({\bm{\alpha}}^{\prime})=0.$
(2)

This integral equation can be discretized in a finite-dimensional subspace of
the many-body Hilbert space, leading to a general eigenvalue problem:

$\displaystyle\sum_{j=1}^{N_{\alpha}}\bigg{(}\mathbf{H}_{ij}-E\mathbf{S}_{ij}\bigg{)}f(\alpha_{j})=0,$
(3)

where

$\displaystyle\mathbf{H}_{ij}$
$\displaystyle=\langle\phi(\alpha_{i})|H|\phi(\alpha_{j})\rangle,$
$\displaystyle\mathbf{S}_{ij}$
$\displaystyle=\langle\phi(\alpha_{i})|\phi(\alpha_{j})\rangle.$

While the GCM is a powerful method, it can be computationally demanding,
especially when dealing with a large number of basis states or when symmetry-
breaking effects are intricate. For strongly correlated molecular systems,
these requirements are often associated with incorporating sufficient electron
correlation effects into GCM calculations, which typically involves using more
sophisticated electronic structure methods than mean-field theory. Techniques
such as configuration interaction (CI), coupled-cluster (CC) methods, and
many-body perturbation theory (MBPT) can be employed to capture electron
correlation effects accurately Shavitt and Bartlett (2009). The basis
functions generated from these methods should be capable of representing the
electron correlation effects relevant to the quantum phenomena under
investigation. Various levels of electron correlation can be included in GCM
calculations, depending on the chosen basis functions and computational
resources available. It is interesting to observe the GCM as a generalization
of the non-orthogonal determinants (e.g., non-orthogonal configuration
interaction). Even though non-entangled single-particle states were employed
in the GCM for simplifying the implementation in our previous report Zheng
_et al._ (2023), the entangled many-body states can also be employed in the
framework to involve complex correlations among particles, making it suitable
for describing systems with complex shapes and collective motion.

##  III GCM vs. VQE: general eigenvalue process vs. constrained optimization

Given the unique properties of GCM, we believe it holds the potential to
address the prevalent issues identified in VQE practices. One primary
advantage is that in GCM, variation occurs on the generating function, not
directly on the scalar generator coordinate. This distinction allows the
variational problem to be tackled through a straightforward eigenvalue
process. This approach stands in contrast to the method of solving a
constrained optimization problem, the benefits of which can even be
illustrated with a simple toy model.

As depicted in Fig. 1, the toy model—comprising two electrons and four spin
orbitals—employs a VQE approach with a Trotterized UCC single-type ansätze:

$\displaystyle|\psi_{\rm
VQE}(\vec{\theta})\rangle=G_{2,4}(\theta_{2})G_{1,3}(\theta_{1})|\phi_{0}\rangle,$
(4)

where the single excitation Givens rotation here is defined as

$\displaystyle G_{p_{i},q_{i}}(\theta_{i})$
$\displaystyle=\exp[\theta_{i}(A_{p_{i},q_{i}}-A_{p_{i},q_{i}}^{\dagger})].$
(5)

with $A_{p,q}=a_{p}^{\dagger}a_{q}$. Essentially, this VQE approach explores a
two-dimensional parameter subspace. The consecutive actions of $G_{2,4}$ and
$G_{1,3}$ generate a total of four distinct configurations, including the
reference. However, since the number of free parameters in this case is fewer
than the number of distinct configurations required, the VQE approach is
constrained by this limitation to be unable to fully explore the subspace to
guarantee the most optimal solution within the subspace (which happened to be
the exact target state), irrespective of the numerical minimizer used.

In contrast, GCM employs the same Givens rotations, separately or in the
product form as discussed by Fukutome Fukutome (1981), to create four
generating functions. These functions correspond to a subspace consisting of
four non-orthogonal superpositioned states, providing a scope that is
sufficiently expansive to encapsulate the target state. Here, each Givens
rotation is essentially equivalent to a UCC single circuit, simplifying the
quantum implementation compared to its double or higher excitation analogs.
It’s important to note that

  * •

Even though the exactness of the disentangled UCC ansätz has been proven
Evangelista  _et al._ (2019), the ansätze employed in GCM, or in the
generalized eigenvalue process, need not be exact. The only prerequisite for
the ansätze used in this process is their ability to generate sufficient
superpositioned states so that a better approximate to the target state can be
found in the corresponding subspace.

  * •

The Givens rotations used in both VQE and GCM may originate from the same set,
ensuring a comparable level of complexity in the quantum circuits involved (to
achieve the same level of exactness, VQE would require a deeper circuit by
incorporating additional Givens rotations). While GCM necessitates more
measurements than a single VQE iteration, a direct comparison of the total
quantum resources utilized between the two methods is not straightforward; it
hinges on both the number of generating functions required in GCM and the
number of iterations needed in VQE.

  * •

In the context of non-orthogonal quantum eigensolver including GCM approaches,
to cover a wide range of correlation, entangled basis sets are generated by
applying UCC operators to the reference single determinant (e.g., the
k-UpCCGSD ansätz Lee  _et al._ (2019), UCC-derived ansatz and a cluster-
Jastrow modification Baek  _et al._ (2023)). These methods often require
domain-specific knowledge available in classical quantum chemistry to
construct high-quality wave functions, and rely on the modified Hadamard test
of Ref. 14 for evaluating the off-diagonal elements at a polynomial cost.
Nevertheless, how to automate the sophisticated ansätz construction process
such that the deployment of the non-orthgonal quantum solver can be easier
even for the non-experts is still a challenge.

These observations then lead us to the pivotal question in advancing a more
efficient GCM-based algorithm for both classical and quantum computation: “How
can we guide the selection of generating functions, or the Givens rotations,
to construct a subspace in the Hilbert space that encompasses the target
state?”

Figure 1: (a) A toy model consists of two-electron (one alpha electron and one
beta electron) in four spin-orbitals (two alpha spin-orbitals and two beta
spin-orbitals), where the spin-flip transition is assumed forbidden. (b) The
projection of the exact wave function on each configuration (yellow shadown).
A constrained optimization of the free parameters will put limits on the
projection (dashed line). (c) Comparative demonstration between GCM and VQE.
The GCM generates a set of non-orthogonal bases, called generating functions.
Then the GCM explores the projection of the system on these generating
functions and solves a corresponding generalized eigenvalue problem for target
state and its energy. The conventional VQE essentially explores the parameter
subspace for a given wave function ansätz through numerical optimization that
can be usually constrained by many factors ranging from ansätz inexactness to
barren plateaus and else. For given Givens rotations that generated excited
state configurations, the lowest eigenvalue obtained from GCM guarantees a
lower bound of the most optimal solution from the VQE.

##  IV The Scaling of Number of Generating Functions: From Exponential To
Linear

Figure 2: Number of generating functions as functions of the number of Givens
rotations (denoted as $m$) included in the GCM wave function ansätz at
different truncation levels.

In constructing a non-trivial wave function ansätz for a molecular system that
comprises $n_{e}$ electrons in $N$ spin orbitals, one might consider applying
a sequence of $m\leq n_{e}$ Givens rotations to a reference state
$|\phi_{0}\rangle$ to generate all the possible configurations. Each Givens
rotation defined in Eq. (5) corresponds to a single excitation, and generates
a superposition of no more than two states. The ansätz can then be expressed
as:

$\displaystyle|\psi(\vec{\theta})\rangle=\prod_{i=1}^{m}G_{p_{i},q_{i}}(\theta_{i})|\phi_{0}\rangle,$
(6)

Consequently, the total number of configurations, $n_{c}$ contained within the
superpositioned ansätz (6), where $\vec{\theta}=\\{\theta_{i}|i=1,\cdots,m\\}$
varies, cannot exceed $2^{m}$. This figure represents the maximum number of
generating functions utilized to achieve the most optimal solution within the
corresponding subspace, exhibiting an exponential scaling relative to the
count of Givens rotations applied.

A strategy to potentially generate the full set of $2^{m}$ generating
functions entails applying $1,2,\cdots,N$ Givens rotations separately to the
reference, as shown below:

  * •

$|\phi_{0}\rangle,$

  * •

$G_{p_{i},q_{i}}(\theta_{i})|\phi_{0}\rangle,\leavevmode\nobreak\
\leavevmode\nobreak\ 1\leq i\leq m,$

  * •

$G_{p_{i},q_{i}}(\theta_{i})G_{p_{j},q_{j}}(\theta_{j})|\phi_{0}\rangle,\leavevmode\nobreak\
\leavevmode\nobreak\ 1\leq i<j\leq m,$

  * •

$G_{p_{i},q_{i}}(\theta_{i})G_{p_{j},q_{j}}(\theta_{j})G_{p_{k},q_{k}}(\theta_{k})|\phi_{0}\rangle,\leavevmode\nobreak\
\leavevmode\nobreak\ 1\leq i<j<k\leq m,$

  * •

$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \vdots$

  * •

$\prod_{i=1}^{m}G_{p_{i},q_{i}}(\theta_{i})|\phi_{0}\rangle$

Here, the number of generating functions that employ $k$ Givens rotations is
equivalent to the number of $k$-combinations drawn from the set of $m$ Givens
rotations, aligning with the combinatorial identity:

$\displaystyle\sum_{k=0}^{m}\bigg{(}\begin{array}[]{c}m\\\
k\end{array}\bigg{)}=2^{m}.$ (9)

Importantly, when $m=n_{e}$, this methodology for constructing generating
functions can encapsulate the complete configuration space of the molecular
system, thereby guaranteeing exactness in this limit. This characteristic
bears notable similarity to full configuration interaction treatments, with
the distinctive aspect that this approach constructs a non-orthogonal many-
body basis set.

The strategy outlined above can also be employed to establish a hierarchy of
approximations to the GCM wave function, aiming towards a complete
configuration subspace where the ansätz (6) can be fully projected out. Within
this framework, the approximation at level $k$ ($0\leq k<m$) is defined by a
working subspace that incorporates both the ansätze (6) and some auxiliary
generating functions produced by acting a product of at most $k$ Givens
rotations on the reference. Notably, at level $k=0$, the approximation narrows
to the subspace containing only the ansätz (6) and $|\phi_{0}\rangle$. In this
context, the lowest eigenvalue of the Hamiltonian $H$ can still be less than
or equal to its variational expectation value with respect to the ansätz (6)
obtained by optimizing the rotations in the ansätz. When $k=m$, independent of
the rotations in the generating functions, the configuration subspace has
expanded sufficiently for the ansätz (6) to be fully projected out. For other
approximation levels where $0<k<m$, a variational procedure can be undertaken
before the configurational subspace expansion, aiming to deduce a lower bound
for the optimal expectation value of $H$ with respect to the ansätz (6). In
this fashion, the GCM can act as a subsequent step to further refine the
energy, especially when the VQE and its ADAPT version, imposing constraints
due to the inexactness of the ansätz and local minima.

It’s important to recognize that the generating functions, including the
ansätz, are not generally orthogonal. This non-orthogonality necessitates
extra resources for evaluating the associated overlap matrix. Additionally, if
there are higher-order Givens rotations involved, the way of the configuration
subspace expansion is influenced by the ordering of the Givens rotations in
the ansätz. However, this ordering limitation can be partially circumvented by
permitting all permutations of the Givens rotations within the auxiliary
generating functions. Such a step increases the measurement pre-factors
without affecting the scaling behavior. For instance, complete permutations
can be permitted when incorporating the auxiliary generating functions
produced by the product of two Givens rotations respectively corresponding to
one single excitation and one double excitation.

##  V Integrating GCM with ADAPT Approach

As mentioned in the previous section, the GCM can be implemented as a
subsequent step either after each ADAPT-VQE iteration or following a complete
ADAPT-VQE calculation. In the ensuing discussion, we designate these
approaches as ADAPT-VQE-GCM and ADAPT-VQE-GCM(1), respectively (refer to
Appendix C). The numeral ‘1’ in ADAPT-VQE-GCM(1) signifies a ‘one-shot’ GCM
calculation performed at the conclusion of an ADAPT-VQE calculation to enhance
the results.

To facilitate a more robust selection of the auxiliary generating functions,
we consider the fluctuation of the GCM energy when a new generating function
is included in the basis set. In the VQE context, similar considerations have
led to its ADAPT version. In the GCM framework, several routines can be
proposed for a gradient-based ADAPT-GCM approach. One straightforward method
is to directly compute the GCM energy gradient with respect to the scalar
rotation (see Appendix B). The energy gradient computed in this manner depends
on the corresponding eigenvector, which in turn requires solving the GCM
general eigenvalue problem for every new basis. However, when solving a
generalized eigenvalue problem–the eigenvalue (i.e. the energy) is more
sensitive to the working subspace than to the scalar rotation in the employed
generating functions. This is exemplified in the Lemma 1 of Theorem 1 in
Appendix A, where, for a $2\times 2$ generalized eigenvalue problem, the
choice of the scalar rotations in the employed generating functions is
relatively flexible.

Figure 3: (a) Schematic depiction of the ADAPT-GCM algorithm. The working
subspace and its corresponding configuration subspace are characterized by a
surrogate state $|\psi\rangle_{s}$. (b) In a two-configuration subspace, a
Givens rotation can couple two configurations, effectively generating a
surrogate state, a superposition of two configurations, that can be employed
as a metric to characterize this subspace. The quality of the working
subspace, comprising one configuration and the surrogate state, is not
sensitive to the specific choice of scalar rotation.

Our proposed gradient-based ADAPT-GCM algorithm in this paper is summarized in
Fig. 3a, in which we consider a surrogate product state $|\psi\rangle_{s}$
that, irrespective of its detailed non-singular scalar rotations, primarily
serves as a metric to (1) abstract the change in the GCM working subspace and
(2) approximate the GCM energy gradient calculation. Through this metric, we
can then compute the gradient in a straightforward way similar to that in the
conventional ADAPT-VQE. But significantly different from ADAPT-VQE, we have to
emphasize that this ADAPT-GCM approach skips the optimization of the phases
and its energy is not from the measurement of the expectation value with
respect to the surrogate state, but obtained from solving a generalized
eigenvalue problem in the working subspace.

In the proposed ADAPT-VQE-GCM(1), the number of bases in the working subspace
equals to the number of Givens rotations in the ansätz generated at the last
ADAPT-VQE iteration, incremented by one. This increment accounts for the
product of these Givens rotations corresponding to the ansätz itself.
Conversely, in the ADAPT-VQE-GCM and ADAPT-GCM approaches, the number of bases
is twice the number of iterations. To illustrate this, consider the $n$-th
iteration: when the Givens rotation $G_{k}$ with the greatest gradient is
selected, two associated basis vectors,

  * •

$G_{k}(\theta_{k})|\phi_{\rm HF}\rangle$ and

  * •

$G_{k}(\theta_{k})\prod_{i=n-1}^{1}G_{i}(\theta_{i})|\phi_{\rm HF}\rangle$

are added to the working subspace. So the number of bases is increased by two
in every iteration. The key difference is that in ADAPT-GCM, the phase is set
to a constant scalar (e.g., $\pi/4$) for any Givens rotation in the ansätz,
while in ADAPT-VQE-GCM, the phases $\theta_{1}$ to $\theta_{n-1}$ are
optimized during the $(n-1)$-th iteration.

##  VI Numerical experiments

We have performed numerical experiments on four molecular systems in six
geometries for testing convergence performance of the proposed adapt versions
of GCM approaches under different electronic configurations. For example, for
$H_{4}$ and $H_{6}$ molecular systems, we consider two geometries for each,
i.e. linear and square geometries for $H_{4}$, and two $H_{6}$ geometries with
different $H-H$ length. The more pronounced quasidegeneracy in the square
$H_{4}$ and the stretched $H_{6}$ exhibit strong static correlations that are
challenging for traditional single-reference methods Jankowski and Paldus
(1980). The operator pool employed in the ADAPT approaches consists of all
possible singly and doubly spin-adapted excitation operators. In all the
ADAPT-VQE simulations, the convergence is achieved when the norm of the
gradient vector is less than $10^{-4}$. Here, we use the more strict criteria
than in Ref. 56 to show the full picture of the performance of ADAPT-VQE.
Other criteria, such as the magnitude of the change of the lowest eigenvalue,
can be used for earlier convergence with the assistance of ADAPT-VQE-GCM. More
details are listed in Appendix E.

The performance results are shown in Fig. 4, where the conventional ADAPT-VQE
results are also included for comparison. Generally speaking, the GCM-related
approaches provide several recipes for locating a lower bound to the ADAPT-VQE
energy. As can be seen, in all the test cases, the ADAPT-VQE-GCM and ADAPT-GCM
approaches exhibit a steady monotonic converging curve until close to the
machine precision. The “one-shot” GCM, ADAPT-VQE-GCM(1), also slightly
improves the ADAPT-VQE energy when the latter converges. Notably, the ADAPT-
VQE-GCM gives a lower energy than the ADAPT-VQE _at every single iteration_ ,
and thus converges faster. As a trade-off of skipping the phase optimization,
the ADAPT-GCM generally requires more iterations to converge than ADAPT-VQE
and ADAPT-VQE-GCM, but both ADAPT-GCM and ADAPT-VQE-GCM are capable to achieve
exact numerical solutions for all the testing cases as indicated by the
corresponding energy difference from the exact are $\leq 10^{-13}$ a.u. at the
convergence.

One issue in the GCM simulations is the numerical instability of the
generalized eigenvalue problem. It is possible that the overlap matrix
$\mathbf{S}$ becomes numerically “singular ” at certain iterations. In
classical computing, the issue of singularity can be perfectly avoided by
orthogonalizing the generating function basis set. However, in quantum
computing, the quantum version of the orthogonalization algorithm like the
quantum Gram-Schmidt process usually requires extra quantum resources, such as
QRAM Zhang  _et al._ (2021). Assuming the orthogonalization is enabled, the
number of iterations in the employed ADAPT procedures can be greatly reduced
(see Appendix E). In the present study, when a singular overlap matrix is
encountered, we choose to reconstruct an approximate but more stable overlap
matrix by imposing a small disturbance ($\leq 10^{-14}$) on all diagonal
entries of the overlap matrix.

Figure 4: Convergence performance of ADAPT-VQE-GCM, ADAPT-VQE-GCM(1), and
ADAPT-GCM in computing the ground states of six molecules. For comparison,
results from ADAPT-VQE are also included. The shaded region is when error is
smaller than the chemical accuracy. ADAPT-VQE and ADAPT-VQE-GCM are terminated
when the sum of the magnitudes of all gradients falls below $10^{-4}$. ADAPT-
GCM is set to converge if the energy differences are under $10^{-6}$ a.u. for
a certain number of consecutive iterations, where this number depends on the
size of the operator pool.

##  VII Quantum resource comparison

In comparison to ADAPT-VQE, the ADAPT-VQE-GCM and ADAPT-VQE-GCM(1) require
additional quantum resource to improve the ADAPT-VQE energy. However, the
resource comparison between ADAPT-VQE and ADAPT-GCM is not straightforward.
Therefore, in this section, we focus on the differences in quantum resources,
specifically the number of CNOT gates and measurements required between ADAPT-
VQE and ADAPT-GCM.

The number of CNOT gates in these two ADAPT approaches depends on how the
basis sets are prepared at each iteration and the number of iterations needed
to achieve a certain accuracy. In Fig. 5, we compare the CNOT counts for these
two ADAPT approaches at each iteration, both before and after quantum circuit
transpilation. For relatively small test cases, the CNOT counts in ADAPT-VQE
state preparation are typically smaller than those in ADAPT-GCM, which
includes only the single and product bases. For larger systems with stronger
correlation, higher accuracy in ADAPT-VQE may require more CNOT gates than in
ADAPT-GCM. It is important to note that this comparison focuses only on state
or basis preparation in the two ADAPT approaches and does not account for the
number of optimizations governed by the minimizer in ADAPT-VQE, which could
potentially increase the CNOT counts. The actual numbers of classical
optimizations used in each ADAPT iteration in numerical experiments are
recorded in Fig. 6. Additionally, the circuit used in obtaining the
$\mathbf{H}$ and $\mathbf{S}$ matrices in ADAPT-GCM differs from the Hadamard
circuit used in ADAPT-VQE. Note that obtaining an off-diagonal matrix element
requires more CNOT gates than the diagonal ones since two bases need to be
prepared simultaneously in the former.

Regarding the number of measurements in the fault tolerant quantum
computation, ADAPT-VQE requires measurements for two parts, VQE and gradient,
at each iteration. The number of the measurements in VQE is not deterministic
and depends on the average number of optimizations in the minimizer,
$\tilde{N}_{\rm opt}$, at each iteration. Assuming there are $N_{\rm iter}$
iterations, the number of measurements in this part is $\tilde{N}_{\rm
opt}\times N_{\rm iter}$. Meanwhile, the number of gradient measurements
roughly equals the size of the pool, which is the same as $N_{\rm iter}$,
multiplied by the number of terms in the Hamiltonian, $N_{\rm term}$. The
latter number can be further reduced by applying grouping techniques
extensively discussed in the previous report Zheng  _et al._ (2023).
Therefore, the total number of measurements in ADAPT-VQE scales as
$\mathcal{O}(\tilde{N}_{\rm opt}\times N_{\rm iter}+N_{\rm term}\times N_{\rm
iter})$.

In ADAPT-GCM, there are also two parts of measurements: the gradient
measurements and the measurements for constructing the projected Hamiltonian
$\mathbf{H}$ and overlap matrices $\mathbf{S}$ at each iteration. The former
is the same as in ADAPT-VQE. The number of measurements for constructing the
entire $\mathbf{H}$ and $\mathbf{S}$ matrices scales as $\mathcal{O}(N_{\rm
GF}^{2})$, with $N_{\rm GF}$ being the number of generating functions.
However, at each iteration, since only new matrix elements associated with new
generating functions are required, the actual number of measurements only
scales as $\mathcal{O}(N_{\rm GF})$ if the number of the generating functions
grows linearly with the number of iterations, as studied in this paper. Since
$N_{\rm GF}\sim\mathcal{O}(N_{\rm iter})$, the measurement comparison between
the two ADAPT approaches then boils down to $\mathcal{O}(\tilde{N}_{\rm
opt}\times N_{\rm iter(VQE)})$ vs. $\mathcal{O}(N_{\rm iter(GCM)}^{2})$.

The above discussion is within the fault-tolerance context. In real quantum
measurements, noise from various sources can be significant. For solving the
generalized eigenvalue problem, the uncertainty in the finite shots of
measurements is often a concern. In our previous study Zheng  _et al._ (2023),
assuming the measurement data to be normally distributed, we’ve shown that
stronger correlation requires orders of magnitude more measurements than
relatively weak correlation. Additionally, it has been shown that finite
sampling error in the overlap $\mathbf{S}$ matrix seems more influential than
in the $\mathbf{H}$ matrix. A more detailed discussion on this topic in the
ADAPT-GCM context is given in Appendix F.

Figure 5: Energy difference as a function of the number of CNOTs in the
corresponding circuits of newly selected ansätze. Note that the number of
CNOTs for ADAPT-GCM at each iteration is the sum of the selected single basis
and a product basis. All circuits are generated from the order-1 trottered
Fermionic operators by Qiskit Qiskit contributors (2023). The reduced counts
are obtained using qiskit.transpile with optimization level 3, which includes
canceling back-to-back CNOTs, “commutative cancellation,” and unitary
synthesis. The implementation of Givens rotation is discussed in Appendix D.
Figure 6: The performance of classical minimizer (in terms of number of
optimization) at each ADAPT-VQE iteration.

##  VIII Conclusion and outlook

In this study, we have demonstrated a new hybrid quantum chemical approach for
accurate and efficient quantum simulations. This approach represents another
subspace expansion method, employing non-orthogonal Givens rotations as bases
for subspace expansion. This differs from typical constrained optimization
methods such as VQE and its ADAPT variant. Compared with other subspace
expansion approaches, our approach automates the basis selection. Compared
with the ADAPT-VQE, our Givens rotation-based method uses the same ansätze
pool but solves a generalized eigenvalue problem within a working subspace.
Here, the global minimum ensures a lower bound to the ADAPT-VQE solutions. Our
results indicate that the numerical performance of the new ADAPT-GCM surpasses
that of the ADAPT-VQE method across several benchmark systems, particularly in
instances of strong electron correlation.

Unlike other subspace expansion methods and those that feature solving
generalized eigenvalue problems, the GCM-based approach is size-extensive.
This is due to the fact that the bases are essentially UCC ansätze. However,
since the bases are generally non-orthogonal, sampling in the overlap matrix
may have a greater impact than in the projected Hamiltonian matrix. Despite
this, the GCM-based approach can be seamlessly applied directly to qubit
space. This is achieved by using the exponential of antihermitian Pauli
strings as generating functions, thereby improving scalability for larger
systems. A similar approach has been reported in the development of qubit-
ADAPT-VQE as a scalable version of ADAPT-VQE Tang  _et al._ (2021).

The simplicity of the GCM Ansatz renders it an appealing solver in constant-
circuit-depth quantum algorithms, such as the QFlow algorithm Kowalski and
Bauman (2023). However, utilizing the GCM solver in QFlow necessitates robust
procedures for approximating GCM by UCC-type ansätz, which can be achieved
through the maximization of the overlap $|\langle\psi_{\rm GCM}|\psi_{\rm
UCC}\rangle|$.

##  IX Data and code availability

The data and the code for the related numerical simulation is available in
GitHub repository Zheng  _et al._ (2022).

##  X Acknowledgements

M. Z., 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 under FWP 72689. A. L. was supported by the U.S.
Department of Energy, Office of Science, National Quantum Information Science
Research Centers, Quantum Science Center. X. Y. was supported by National
Science Foundation CAREER DMS-2143915. M. Z. and X. Y. both also were
supported by Defense Advanced Research Projects Agency as part of the project
W911NF2010022: The Quantum Computing Revolution and Optimization: Challenges
and Opportunities. This research used resources of the Oak Ridge Leadership
Computing Facility, which is a DOE Office of Science User Facility supported
under Contract DE-AC05-00OR22725. The Pacific Northwest National Laboratory is
operated by Battelle for the U.S. Department of Energy under Contract DE-
AC05-76RL01830.

##  Appendix A Mapping optimization to a configuration subspace

Theorem 1. For a Hamiltonian $H$ and a parametrized normalized ansätz
$|\Phi(\bm{\theta})\rangle=U(\bm{\theta})|\phi\rangle$, where $|\phi\rangle$
is a fixed reference state and $U(\bm{\theta})$ is a unitary operator
dependent on the parameter
$\bm{\theta}=\\{\theta_{1},\theta_{2},\cdots\theta_{P}\\}$, the minimal
expectation value

$\displaystyle\epsilon=\min_{\bm{\theta}}\langle\phi|U^{\dagger}(\bm{\theta})|H|U(\bm{\theta})|\phi\rangle$
(10)

is the upper bound of the minimal eigenvalue of $H$ in an arbitrary linearly
independent $M$-configuration subspace
$\\{\psi_{1},\psi_{2},\cdots,\psi_{M}\\}$ ($M\geq P$) where
$U(\bm{\theta})|\phi\rangle$ can be completely projected out, i.e.,

$\displaystyle
U(\bm{\theta})|\phi\rangle=\sum_{i=1}^{M}c_{i}(\bm{\theta})|\psi_{i}\rangle.,$
(11)

where

$\displaystyle
c_{i}(\bm{\theta})=\sum_{j=1}^{M}(\mathbf{S}^{-1})_{ij}\langle\psi_{j}|U(\bm{\theta})|\phi\rangle$
(12)

and $\mathbf{S}$ is the the overlap matrix with

$\displaystyle\mathbf{S}_{ij}=\langle\psi_{i}|\psi_{j}\rangle.$ (13)

Proof : Given (11), $U(\bm{\theta})|\phi\rangle$ is one state in the linearly
independent $M$-configuration subspace, thus the minmal expectation value from
minimizing (10) is constrained by the form of $U(\bm{\theta})$, and does not
necessarily give the minimal expectation value of $H$ in the subspace.
Instead, searching the minimal expectation value of $H$ in the subspace can be
re-formulated to minimizing a generalized Rayleigh quotient for the projected
Hamiltonian matrix $\mathbf{H}$ and the overlap matrix $\mathbf{S}$

$\displaystyle\min_{f}\leavevmode\nobreak\
R(f)=\frac{f^{\ast}\mathbf{H}f}{f^{\ast}\mathbf{S}f},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \text{subjec
to}\leavevmode\nobreak\ \leavevmode\nobreak\ f^{\ast}\mathbf{S}f=1,$ (14)

with matrix elements defined as

$\displaystyle\mathbf{H}_{ij}$
$\displaystyle=\langle\psi_{i}|H|\psi_{j}\rangle,$ (15)
$\displaystyle\mathbf{S}_{ij}$
$\displaystyle=\langle\psi_{i}|\psi_{j}\rangle.$ (16)

By introducing a Lagrange multiplier and forming the Lagrangian,

$\displaystyle
L(f,\epsilon)=f^{\ast}\mathbf{H}f-\epsilon(f^{\ast}\mathbf{S}f-1),$ (17)

the minimization corresponds to taking the derivative of $L(f,\epsilon)$ with
respect to $f$ and setting it to zero, leading to a generalized eigenvalue
equation:

$\displaystyle\mathbf{H}f=\epsilon\mathbf{S}f,$ (18)

where the lowest eigenvalue corresponds to the global minimum in the
$M$-configuration subspace.

Lemma 1 Consider a unitary operation $U$ parametrized by a single scalar
$\theta$ in a variational ansätze. If, for any value of $\theta$, the action
of $U$ on $|\phi\rangle$ produces at most two distinct configurations, then
the task of minimizing
$\langle\phi|U^{\dagger}(\theta)|H|U(\theta)|\phi\rangle$ can be equivalently
approached by solving the following $2\times 2$ generalized eigenvalue
problem:

$\displaystyle\left(\begin{array}[]{cc}\langle\phi|U^{\dagger}(\theta_{1})HU(\theta_{1})|\phi\rangle&\langle\phi|U^{\dagger}(\theta_{1})HU(\theta_{2})|\phi\rangle\\\
\langle\phi|U^{\dagger}(\theta_{2})HU(\theta_{1})|\phi\rangle&\langle\phi|U^{\dagger}(\theta_{2})HU(\theta_{2})|\phi\rangle\end{array}\right)f$
(21) $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
=\epsilon\left(\begin{array}[]{cc}1&\langle\phi|U^{\dagger}(\theta_{1})U(\theta_{2})|\phi\rangle\\\
\langle\phi|U^{\dagger}(\theta_{2})U(\theta_{1})|\phi\rangle&1\end{array}\right)f.$
(24)

Here $\theta_{1}$ and $\theta_{2}$ are chosen such that
$\|\langle\phi|U^{\dagger}(\theta_{1})U(\theta_{2})|\phi\rangle\|^{2}\neq 1$.
The lowest eigenvalue of the $2\times 2$ problem yields the minimal energy
within the defined subspace.

Proof : For any $\theta$, let’s express

$\displaystyle
U(\theta)|\phi\rangle=a(\theta)|\phi_{1}\rangle+b(\theta)|\phi_{2}\rangle$
(25)

where $|\phi_{1}\rangle$ and $|\phi_{2}\rangle$ are orthogonal, i.e.,
$\langle\phi_{1}|\phi_{2}\rangle=0$, and $a(\theta)$ and $b(\theta)$ are
scalar functions of $\theta$ with $|a(\theta)|^{2}+|b(\theta)|^{2}=1$

According to Theorem 1, the minimization of
$\langle\phi|U^{\dagger}(\theta)|H|U(\theta)|\phi\rangle$ over all $\theta$
can be recast into a $2\times 2$ eigenvalue problem with the Hamiltonian
matrix

$\displaystyle\mathbf{H}=\left(\begin{array}[]{cc}\langle\phi_{1}|H|\phi_{1}\rangle&\langle\phi_{1}|H|\phi_{2}\rangle\\\
\langle\phi_{2}|H|\phi_{1}\rangle&\langle\phi_{2}|H|\phi_{2}\rangle\end{array}\right).$
(28)

Alternatively, the expansion of $U(\theta)|\phi\rangle$ can also be expanded
by non-orthogonal basis, e.g.

$\displaystyle
U(\theta)|\phi\rangle=c_{1}U(\theta_{1})|\phi\rangle+c_{2}U(\theta_{2})|\phi\rangle,$
(29)

where $\theta_{1}$ and $\theta_{2}$ are two scalars satisfying
$\|\langle\phi|U^{\dagger}(\theta_{1})U(\theta_{2})|\phi\rangle\|^{2}\neq 1$,
and $c_{1}$ and $c_{2}$ depends on the choice of
$\theta,\theta_{1},\theta_{2}$ and can be obtained through

$\displaystyle(a(\theta_{1})b(\theta_{2})-b(\theta_{1})a(\theta_{2}))c_{1}$
$\displaystyle=a(\theta)b(\theta_{2})-b(\theta)a(\theta_{2}),$ (30)
$\displaystyle(a(\theta_{2})b(\theta_{1})-b(\theta_{2})a(\theta_{1}))c_{2}$
$\displaystyle=a(\theta)b(\theta_{1})-b(\theta)a(\theta_{1}).$ (31)

If the expansion is done in such a non-orthogonal subspace, the minimization
of $\langle\phi|U^{\dagger}(\theta)|H|U(\theta)|\phi\rangle$ transforms into a
generalized eigenvalue problem (18) with matrices defined as

$\displaystyle\mathbf{H}=\left(\begin{array}[]{cc}\langle\phi|U^{\dagger}(\theta_{1})HU(\theta_{1})|\phi\rangle&\langle\phi|U^{\dagger}(\theta_{1})HU(\theta_{2})|\phi\rangle\\\
\langle\phi|U^{\dagger}(\theta_{2})HU(\theta_{1})|\phi\rangle&\langle\phi|U^{\dagger}(\theta_{2})HU(\theta_{2})|\phi\rangle\end{array}\right).$
(34)

and

$\displaystyle\mathbf{S}=\left(\begin{array}[]{cc}1&\langle\phi|U^{\dagger}(\theta_{1})U(\theta_{2})|\phi\rangle\\\
\langle\phi|U^{\dagger}(\theta_{2})U(\theta_{1})|\phi\rangle&1\end{array}\right).$
(37)

Solving this generalized eigenvalue problem allows us to determine the
combination of $U(\theta_{1})|\phi\rangle$ and $U(\theta_{2})|\phi\rangle$
that corresponds to the lowest expectation value for $H$ in this two-
configuration subspace.

##  Appendix B Energy gradient in ADAPT-GCM framework

In the GCM framework, the $k$-th eigenvalue of the projected $H$,
$\epsilon_{k}$, can be represented as

$\epsilon_{k}=\frac{f^{T}_{k}\mathbf{H}f_{k}}{f^{T}_{k}\mathbf{S}f_{k}},$ (38)

where $f_{k}$ is the corresponding normalized eigenvectors,
$f^{T}_{k}f_{k}=1$. Assume the working subspace comprises $M$ generating
functions, $\\{|\psi_{i}\rangle\\}_{1\leq i\leq M}$, in which $M-1$ are
generated by acting one Givens rotation on the reference, and only one is
generated by acting the product form of Givens rotation on the reference, i.e.

$\displaystyle|\psi_{i}\rangle=\left\\{\begin{array}[]{rc}G_{p_{i},q_{i}}(\theta_{i})|\phi_{0}\rangle,&1\leq
i\leq M-1\\\
\prod_{j=1}^{M-1}G_{p_{j},q_{j}}(\theta_{j})|\phi_{0}\rangle,&i=M\end{array}\right.\leavevmode\nobreak\
\leavevmode\nobreak\ .$ (41)

Then ${\bf H}$ and ${\bf S}$ are $M\times M$ matrices projected in this
working subspace with their matrix elements being defined as

$\displaystyle{\bf H}_{ij}$ $\displaystyle=\langle\psi_{i}|H|\psi_{j}\rangle,$
(42) $\displaystyle{\bf S}_{ij}$
$\displaystyle=\langle\psi_{i}|\psi_{j}\rangle,$ (43)

and the derivative of the eigenvalue $\epsilon_{k}$ with respect to the
rotation $\theta_{s}$ can be computed through

$\displaystyle\frac{\partial\epsilon_{k}}{\partial\theta_{s}}=\frac{\langle\frac{\partial\mathbf{H}}{\partial\theta_{s}}\rangle_{k}\cdot\langle\mathbf{S}\rangle_{k}-\langle\mathbf{H}\rangle_{k}\cdot\langle\frac{\partial\mathbf{S}}{\partial\theta_{s}}\rangle_{k}}{\langle\mathbf{S}\rangle_{k}^{2}},$
(44)

with $\langle\mathbf{O}\rangle_{k}=f^{T}_{k}\mathbf{O}f_{k}$ for matrix
$\mathbf{O}$. The matrix elements of
$\frac{\partial\mathbf{S}}{\partial\theta_{s}}$ and
$\frac{\partial\mathbf{H}}{\partial\theta_{s}}$ are the following,

  * •

if neither $|\psi_{i}\rangle$ nor $|\psi_{j}\rangle$ is parametrized by
$\theta_{s}$

$\displaystyle\frac{\partial(\mathbf{H}_{ij})}{\partial\theta_{s}}=\frac{\partial(\mathbf{S}_{ij})}{\partial\theta_{s}}=0,$
(45)

  * •

if only $|\psi_{i}\rangle$ is parametrized by $\theta_{s}$ and $i\neq j$

$\displaystyle\frac{\partial(\mathbf{H}_{ij})}{\partial\theta_{s}}$
$\displaystyle=\langle\psi_{i}|(a_{q_{s}}^{\dagger}a_{p_{s}}-a_{p_{s}}^{\dagger}a_{q_{s}})H|\psi_{j}\rangle,$
(46) $\displaystyle\frac{\partial(\mathbf{S}_{ij})}{\partial\theta_{s}}$
$\displaystyle=\langle\psi_{i}|(a_{q_{s}}^{\dagger}a_{p_{s}}-a_{p_{s}}^{\dagger}a_{q_{s}})|\psi_{j}\rangle,$
(47)

  * •

if only $|\psi_{j}\rangle$ is parametrized by $\theta_{s}$ and $i\neq j$

$\displaystyle\frac{\partial(\mathbf{H}_{ij})}{\partial\theta_{s}}$
$\displaystyle=\langle\psi_{i}|H(a_{p_{s}}^{\dagger}a_{q_{s}}-a_{q_{s}}^{\dagger}a_{p_{s}})|\psi_{j}\rangle,$
(48) $\displaystyle\frac{\partial(\mathbf{S}_{ij})}{\partial\theta_{s}}$
$\displaystyle=\langle\psi_{i}|(a_{p_{s}}^{\dagger}a_{q_{s}}-a_{q_{s}}^{\dagger}a_{p_{s}})|\psi_{j}\rangle,$
(49)

  * •

if both $|\psi_{j}\rangle$ are parametrized by $\theta_{s}$

$\displaystyle\frac{\partial(\mathbf{H}_{ij})}{\partial\theta_{s}}$
$\displaystyle=\langle\psi_{i}|[H,(a_{p_{s}}^{\dagger}a_{q_{s}}-a_{q_{s}}^{\dagger}a_{p_{s}})]|\psi_{j}\rangle,$
(50) $\displaystyle\frac{\partial(\mathbf{S}_{ij})}{\partial\theta_{s}}$
$\displaystyle=0.$ (51)

As can be seen matrices $\frac{\partial\mathbf{S}}{\partial\theta_{s}}$ and
$\frac{\partial\mathbf{H}}{\partial\theta_{s}}$ are very sparse with the non-
zero elements only residing in the column and row whose index is equal to $p$,
which can greatly reduce the computational complexity. The gradient employed
in the proposed ADAPT-GCM in Fig. 3 can be considered as an approximation to
the above gradient calculations, where only (50) and (51) are taken into
consideration.

##  Appendix C ADAPT-VQE-GCM and ADAPT-VQE-GCM(1) algorithms

The main structure of ADAPT-VQE-GCM in Alg. 1 is built upon the a complete
ADAPT-VQE process from Ref. 56. In fact, the essential use of ADAPT-VQE is to
provide the unitary operators and phase values (except the phase for the
latest selected ansätze) to form the bases for GCM. ADAPT-VQE-GCM(1) only
includes the selected single Givens rotations from all previous ADAPT
iterations and the final product of the Givens rotations when ADAPT-VQE
converges in the basis set.

Algorithm 1 ADAPT-VQE-GCM

1:Operator pool $\\{G_{l}\\}_{(l)\in L}$ with $L$ the orbital index pool,
initial state $\ket{\psi_{(0)}^{VQE}}:=\ket{\phi_{HF}}$, VQE optimal parameter
set $\Theta^{*}_{(0)}=\emptyset$, GCM basis set $\mathbf{B}=\emptyset$, and
Hamiltonian $H$

2:$n\leftarrow 0$

3:while convergence criteria is not met do

4: $n\leftarrow n+1$

5: Select operator: $G_{(n)}\leftarrow\max_{l\in
L}\left\\{\frac{\partial\epsilon}{\partial\theta_{l}}\big{|}_{\theta_{l}=0}=\bra{\psi^{VQE}_{(n-1)}}\left[H,A_{l}\right]\ket{\psi^{VQE}_{(n-1)}}\right\\}$;

6: Prepare VQE ansätz and optimize free parameters:
$\ket{\tilde{\psi}_{(n)}^{VQE}(\Theta_{(n)})}=\prod_{i=n}^{1}G_{(i)}(\theta_{i})\ket{\phi_{HF}}$;
$\Theta^{*}_{(n)}={\rm
argmin}_{\Theta_{(n)}}\langle\tilde{\psi}_{(n)}^{VQE}(\Theta_{(n)})|H\ket{\tilde{\psi}_{(n)}^{VQE}(\Theta_{(n)})}$;
$\ket{\psi_{(n)}^{VQE}}\leftarrow\ket{\tilde{\psi}_{(n)}^{VQE}(\Theta^{\ast}_{(n)})}$;

7: Expand the basis set $B$: $\mathbf{B}\leftarrow
G_{(n)}(\theta^{\ast}_{n})|\phi_{HF}\rangle$;
$\mathbf{B}\leftarrow\ket{\psi_{(n)}^{VQE}}$;

8: Compute matrices and solve for eigenvalues:
$\mathbf{H}=\mathbf{B}^{T}H\mathbf{B}$; $\mathbf{S}=\mathbf{B}^{T}\mathbf{B}$;
$\mathbf{H}f=\epsilon\mathbf{S}f$;

9:end while

10:return eigenvalues $\epsilon$ after convergence

##  Appendix D Implementing the Givens rotations

For skew-Hermitian operators

$\displaystyle A_{p_{i},q_{i}}$
$\displaystyle=a_{p_{i}}^{\dagger}a_{q_{i}}-a_{q_{i}}^{\dagger}a_{p_{i}},\leavevmode\nobreak\
\leavevmode\nobreak\ p_{i}<q_{i}$ (52) $\displaystyle
A_{p_{i},r_{i},q_{i},s_{i}}$
$\displaystyle=a_{p_{i}}^{\dagger}a_{r_{i}}^{\dagger}a_{q_{i}}a_{s_{i}}-a_{s_{i}}^{\dagger}a_{q_{i}}^{\dagger}a_{r_{i}}a_{p_{i}},\leavevmode\nobreak\
\leavevmode\nobreak\ q_{i}<s_{i}<p_{i}<r_{i}$ (53)

the number of CNOT gates in their corresponding Givens rotation circuits can
be estimated as follows. For single excitations, Eq. (52), its corresponding
Givens rotation circuit for two adjacent indices is given in Fig. 1, i.e., two
CNOT gates. If the two indices are not adjacent, $q_{i}-p_{i}-1$ FSWAP gates
will be employed. If each FSWAP contains three CNOT gates, Hashim  _et al._
(2021) the total number of CNOT gates for implementing a single excitation
Givens rotation would be $3(q_{i}-p_{i})-1$. For double excitations, Eq. (53),
the typical implementation involves the Trotterization of the exponential of
the Pauli representation of the second quantized expression. For example,
through Jordan-Wigner transformation,

$\displaystyle a_{p_{i}}=\frac{1}{2}(X_{p_{i}}+{\rm
i}Y_{p_{i}})\prod_{k=0}^{p_{i}-1}Z_{k},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
a_{p_{i}}^{\dagger}=\frac{1}{2}(X_{p_{i}}-{\rm
i}Y_{p_{i}})\prod_{k=0}^{p_{i}-1}Z_{k},$ (54)

the exponential of Eq. (53) can be expressed as

$\displaystyle\exp(\theta A_{p_{i},r_{i},q_{i},s_{i}})=\exp\bigg{(}-{\rm
i}\frac{\theta}{8}(X_{q_{i}}Y_{s_{i}}X_{p_{i}}X_{r_{i}}+Y_{q_{i}}X_{s_{i}}X_{p_{i}}X_{r_{i}}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+Y_{q_{i}}Y_{s_{i}}Y_{p_{i}}X_{r_{i}}+Y_{q_{i}}Y_{s_{i}}X_{p_{i}}Y_{r_{i}}-X_{q_{i}}X_{s_{i}}Y_{p_{i}}X_{r_{i}}-X_{q_{i}}X_{s_{i}}X_{p_{i}}Y_{r_{i}}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
-Y_{q_{i}}X_{s_{i}}Y_{p_{i}}Y_{r_{i}}-X_{q_{i}}Y_{s_{i}}Y_{p_{i}}Y_{r_{i}})\prod_{k=p_{i}+1}^{r_{i}-1}Z_{k}\prod_{l=q_{i}+1}^{s_{i}-1}Z_{l}\bigg{)},$
(55)

whose first Trotter approximation can be implemented using 48 CNOT gate if
$p_{i},r_{i},q_{i},s_{i}$ are adjacent to each other. Yordanov  _et al._
(2020) If not, longer CNOT ladders will be employed to connect mapped qubits.
For all the singlet operators in the spin-adapted pool in Section VI, under
the Jordan-Wigner transformation, the average number of CNOTs gates of an
individual operator is

$\displaystyle\bar{N}_{Pauli\;per\;string}\times\bar{N}_{string}\times 2$ (56)

where $\bar{N}_{Pauli\;per\;string}$ is the average number of non-identity
Pauli gates in all strings, $\bar{N}_{Paulistring}$ is the average number of
Pauli strings in an operator after mapping, and the multiplier $2$ comes from
the symmetries of the CNOT ladders.

##  Appendix E Other computational details

The implementation of ADAPT-VQE is the original code in Ref. 56 from the
corresponding repository Grimsley  _et al._ (2022), and its computation is
based on the SciPy linear algebra package Virtanen  _et al._ (2020). To have
the same foundation for comparison, calculations in Section VI for three GCM
related algorithms follow the same manner. However, the quantum resources
estimation, shown in Section VII, utilizes Qiskit Qiskit contributors (2023).
All the numerical experiments are carried on a MacBook Pro with 2.6 GHz 6-Core
Intel Core i7 CPU and Intel UHD Graphics 630 GPU.

As mentioned in Section VI, some small disturbances are used to mitigate
numerical instability in the generalized eigenvalue problems of ADAPT-GCM.
This problem can be resolved with orthgonalization of basis vectors, while the
current quantum implementation requires special hardware like QRAM Zhang  _et
al._ (2021), which is beyond the scope of this paper. To illustrate the
effects, numerical simulations with classical orthgonalization on both $H_{6}$
examples is conducted and demonstrated in Fig. 7. When both ADAPT-VQE-GCM and
ADAPT-GCM converge in either example, the number of orthogonal bases are both
at around 100.

Figure 7: Number of bases used in each iteration before and after
orthgonalization in ADAPT-VQE-GCM and ADAPT-VQE in both $H_{6}$ examples. The
“orth” in the labels mean the calculation was done with the “orthogonalized
bases.”

##  Appendix F Randomness due to finite number of shots

Recall that in both ADAPT-GCM and GCM part of the ADAPT-VQE-GCM, we compute
project Hamiltonian, $\mathbf{H}$, and overlap matrix, $\mathbf{S}$, entry by
entry

$\displaystyle\mathbf{H}_{ij}$
$\displaystyle=\langle\psi_{i}|H|\psi_{j}\rangle=\langle\psi_{HF}|B_{i}^{\dagger}HB_{j}|\psi_{HF}\rangle$
(57) $\displaystyle\mathbf{S}_{ij}$
$\displaystyle=\langle\psi_{i}|\psi_{j}\rangle=\langle\psi_{HF}|B_{i}^{\dagger}B_{j}|\psi_{HF}\rangle,$
(58)

where $B_{i}$ and $B_{j}$ are short-hand notations for some products of Givens
rotation matrices for simplicity. These expectations can be computed by
converting $B_{i}^{\dagger}HB_{j}$ and $B_{i}^{\dagger}B_{j}$ to the sums of
Pauli strings, respectively, then grouping the Pauli strings to do
simultaneous measurements. The quantum and classical complexities and effects
of shot noise have been explicitly discussed in Ref. 50.

Instead, in this section, we illustrate the shot noise calculation when
employing Hadamard test for the expectation computation. Due to the randomness
of the collapse of qubits, we will transfer $\mathbf{H}_{ij}$ to a random
variable where the number of shots controls its variance.

Assume the Pauli decomposition of Hamiltonian
$H=\sum_{k=1}^{N_{term}}c_{k}P_{k}$, $c_{k}\in\mathbb{R}$ for all $k$, has
been provided, we have

$\displaystyle\mathbf{H}_{ij}=\sum_{k=1}^{N_{term}}c_{k}\langle\psi_{HF}|B_{i}^{\dagger}P_{k}B_{j}|\psi_{HF}\rangle.$
(59)

Without the loss of the generality, we only focus on $\mathbf{H}_{ij}$ since
the calculation of $\mathbf{S}_{ij}$ can be regarded as a special case for
(59) when $N_{term}=1$, $c_{k}=1$, and $P_{k}$ is an identity matrix. It is
known that, Hadamard test gives

$\displaystyle\text{Pr}\left(\text{ancillary is }\ket{0}\right)$
$\displaystyle=\frac{1}{2}+\frac{1}{2}\text{Re}\langle\psi_{HF}|B_{i}^{\dagger}P_{k}B_{j}|\psi_{HF}\rangle$
(60) $\displaystyle\text{Pr}\left(\text{ancillary is }\ket{1}\right)$
$\displaystyle=\frac{1}{2}-\frac{1}{2}\text{Re}\langle\psi_{HF}|B_{i}^{\dagger}P_{k}B_{j}|\psi_{HF}\rangle,$
(61)

after measuring the ancillary qubit in the circuit. Thus, a binary random
variable can be defined by

$\displaystyle\footnotesize\Gamma^{(k)}_{l}=\begin{cases}1&\text{if ancillary
is $\ket{0}$ in $l^{th}$ shot when compute $\langle
B_{i}^{\dagger}P_{k}B_{j}\rangle$}\\\ -1&\text{if ancillary is $\ket{1}$ in
$l^{th}$ shot when compute $\langle
B_{i}^{\dagger}P_{k}B_{j}\rangle$}\end{cases}$ (62)

with

$\displaystyle p_{k}$
$\displaystyle:=\mathbb{E}\left[\Gamma^{(k)}_{l}\right]=\text{Re}\langle\psi_{HF}|B_{i}^{\dagger}P_{k}B_{j}|\psi_{HF}\rangle,$
(63) $\displaystyle\text{Var}\left[\Gamma^{(k)}_{l}\right]$
$\displaystyle=\mathbb{E}\left[\left(\Gamma^{(k)}_{l}\right)^{2}\right]-\left(\mathbb{E}\left[\Gamma^{(k)}_{l}\right]\right)^{2}=1-p_{k}^{2},$
(64)

for all $l$ by assuming $\Gamma^{(k)}_{l}$’s are mutually independent and
identically distributed. For simplicity, we omit the indices $i$ and $j$ for
now. In other words, each measurement on the ancillary qubit in the Hadamard
test circuit is modelled by sampling the random variable $\Gamma^{(k)}_{l}$.

We further construct another random variable, omitting indices $i$ and $j$
again,

$\displaystyle\Lambda^{(k)}=\sum_{i=l}^{N^{(k)}_{shot}}\Gamma^{(k)}_{l}$ (65)

where $N^{(k)}_{shot}$ is the number of shots for the Hadamard test circuit of
$\langle B_{i}^{\dagger}P_{k}B_{j}\rangle$, and it is considered as a user-
defined constant. Then, it follows that

$\displaystyle\frac{1}{2}\left(\Lambda^{(k)}+N^{(k)}_{shot}\right)\sim\text{Bin}\left(N^{(k)}_{shot},\frac{1}{2}+\frac{1}{2}p_{k}\right),$
(66)

by realizing $\frac{1}{2}+\frac{1}{2}\Gamma^{(k)}_{i}$ is a standard Bernoulli
random variable. Further combining (59) and (66), another discrete random
variable

$\displaystyle\Xi^{(ij)}=\sum_{k=1}^{N_{term}}\frac{c_{k}}{N^{(k)}_{shot}}\Lambda^{(k)}=\sum_{k=1}^{N_{term}}\frac{c_{k}}{N^{(k)}_{shot}}\sum_{l=1}^{N^{(k)}_{shot}}\Gamma^{(k)}_{l}$
(67)

has

$\displaystyle\mathbb{E}\left[\Xi^{(ij)}\right]$
$\displaystyle=\sum_{k=1}^{N_{term}}\frac{c_{k}}{N^{(k)}_{shot}}N^{(k)}_{shot}p_{k}=\mathbf{H}_{ij}$
(68) $\displaystyle\text{Var}\left[\Xi^{(ij)}\right]$
$\displaystyle=\sum_{k=1}^{N_{term}}\left(\frac{c_{k}}{N^{(k)}_{shot}}\right)^{2}N^{(k)}_{shot}\text{Var}\left[\Gamma^{(k)}_{l}\right]$
$\displaystyle=\sum_{k=1}^{N_{term}}\frac{c^{2}_{k}}{N^{(k)}_{shot}}\left(1-p_{k}^{2}\right).$
(69)

In this case, $\mathbf{H}_{ij}$ has a corresponding random variable
$\Xi^{(ij)}$ where the randomness is caused by the uncertainty from the
collapse of the qubits during the measurement. The role of the number of
shots, as (69) implies, is to govern the scale of the variance. The variance
(69) goes to 0 when every $N_{shot}^{(k)}$ goes to infinity.

Chebyshev’s inequality shows that, for every $a>0$,

$\displaystyle\text{Pr}\left(\left|\Xi^{(ij)}-\mathbf{H}_{ij}\right|\geq
a\right)\leq\frac{\text{Var}\left[\Xi^{(ij)}\right]}{a^{2}}=:\eta_{ij}.$ (70)

Setting $N^{(k)}_{shot}$ to $N_{shot}$ for all $i,j,k$, Table 1 provides the
relation between $N_{shot}$ and $\max_{i,j}\,\eta_{ij}$ when $a$ is fixed at
$10^{-4}$. However, the tightness of the bound is still unclear for this
discrete variable, and this is not sufficient to illustrate the effects of the
disturbance of the matrix entries on the eigenvalues. This leads to numerical
simulations.

Table 1: $N_{shot}$ for the various choices of $\max_{i,j}\,\eta_{ij}$ when $a=10^{-4}$. Molecule | $\max\,\eta_{ij}=5\%$ | $\max\,\eta_{ij}=1\%$
---|---|---
$H4$ (square) | $6.32\times 10^{10}$ | $3.16\times 10^{11}$
$LiH$ | $1.34\times 10^{11}$ | $6.68\times 10^{11}$
$H_{6}$ ($2.0\,\,a.u.)$ | $1.26\times 10^{11}$ | $6.29\times 10^{11}$
$H_{6}$ ($3.5\,\,a.u.)$ | $6.17\times 10^{10}$ | $3.08\times 10^{11}$

Since it is very common to have huge $N^{(k)}_{shot}$ in practice, we can
assume $\Lambda^{(k)}$ is closely approximated by

$\displaystyle\Lambda^{(k)}_{gauss}\sim\text{N}\left(N^{(k)}_{shot}p_{k},\,\,N^{(k)}_{shot}\left(1-p_{k}^{2}\right)\right)$
(71)

and mutually independent for different $k$. Thus, $\Xi^{(ij)}$ is also
approximated by a Gaussian random variable $\Xi^{(ij)}_{gauss}$ with the same
mean and variance as in (68) and (69), respectively.

The idea of importance sampling (IS) can be applied here to significantly
reduce the number of shots. By setting

$\displaystyle N^{(k)}_{shot}:=|c_{k}|\tau$ (72)

for some large constant $\tau$, the random variable

$\displaystyle\Xi^{(ij)}_{IS}=\frac{1}{\tau}\sum_{k=1}^{N_{term}}\text{sign}(c_{k})\Lambda^{(k)}_{gauss}$
(73)

is also Gaussian and has the mean $\mathbf{H}_{ij}$ and the variance

$\displaystyle\text{Var}\left[\Xi^{(ij)}_{IS}\right]=\frac{1}{\tau}\sum_{k=1}^{N_{term}}|c_{k}|\left(1-p_{k}^{2}\right).$
(74)

It is suggested to use larger $N^{k}_{shot}$ and $\tau$ for the overlap matrix
$\mathbf{S}$ than $\mathbf{H}$ to improve the numerical stability. In the
following simulations, $\mathbf{S}$ will always have one hundred times larger
$N^{k}_{shot}$ and $\tau$ than $\mathbf{H}$.

The Monte Carlo simulations for squared $H_{4}$, $LiH$, and both $H_{6}$ are
conducted by treating each matrix entry as a normal random variable follows
(68) and (69). We select a pair of $\mathbf{H}$ and $\mathbf{S}$ from a
specific iteration of ADAPT-GCM of each molecule such that the error of
ground-energy is around $10^{-11}$ to $10^{-14}$ a.u. without the finite-
sampling noise.

However, the initial simulations show that the shot noise severely aggravates
the numerical instability in the eigenvalue problems. This makes us impose
another strategy when solve the eigenvalue problems. Besides adding $10^{-10}$
on the diagonal entries of every overlap matrix, we also discard eigenvalues
that are smaller than an lower bound based on the HF energy, $\left\lfloor
E_{\text{HF}}\right\rfloor-1$. In practice, this bound can be modified from
the energy in the previous iteration.

As summarized in Fig. 8, without IS, $H_{4}$ and $LiH$ require
$N^{(k)}_{shot}=\tau_{simu}$ to be about $10^{7}$ to reach the chemical
accuracy stably while both $H_{6}$ instances need the number to be around
$10^{11.5}$ and $10^{12.5}$, respectively. When IS is applied, setting
$\tau=\tau_{simu}$, $H_{6}$ ($3.5\,\,a.u.$) gives the largest accuracy drop,
but as Table 2 shows, the reduction by IS on the total number of shots for
$\mathbf{H}$ when $N^{k}_{shot}=\tau=\tau_{simu}$ is around $96.2\%$ for
$H_{6}$ ($3.5\,\,a.u.$). So importance sampling is still effectively reduce
the number of shots in all four molecules.

Table 2: Reduction by IS on the total number of shots for $\mathbf{H}$ when $N^{(k)}_{shot}=\tau=\tau_{simu}$ regardless the change of accuracy. Molecule | $N_{term}$ | $\sum_{k}|c_{k}|$ |  $\%$ reduction by IS
---|---|---|---
$H4$ (square) | 529 | 43.4 | $91.8\%$
$LiH$ | 1861 | 65.8 | $96.5\%$
$H_{6}$ ($2.0\,\,a.u.)$ | 2629 | 115.0 | $95.6\%$
$H_{6}$ ($3.5\,\,a.u.)$ | 2629 | 99.2 | $96.2\%$
Figure 8: Monte Carlo simulations for the effects of the shot noise on the
accuracy of the selected iterations in ADAPT-GCM on different molecules. The
error under infinite number of shots are in $10^{-11}$ to $10^{-14}$. The
black and red-shaded area is the $95\%$ confidence interval estimated from 100
runs for each given $\tau_{simu}$. When not using importance sampling,
$N^{k}_{shot}=\tau_{simu}$ for a projected Hamiltonian $\mathbf{H}$ and
$N^{k}_{shot}=100\tau_{simu}$ for an overlap matrix $\mathbf{S}$. With
importance sampling, $N^{k}_{shot}=|c_{k}|\tau_{simu}$ for $\mathbf{H}$ and
$N^{k}_{shot}=100|c_{k}|\tau_{simu}$ for $\mathbf{S}$. Note that the total
number of shots for estimating a matrix entry is different when
$N^{(k)}_{shot}=\tau=\tau_{simu}$.

## References

  * Imada  _et al._ (1998) M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. 70, 1039 (1998).
  * Witzke  _et al._ (2020) R. J. Witzke, D. Hait, K. Chakarawet, M. Head-Gordon, and T. D. Tilley, ACS Catalysis 10, 7800 (2020).
  * Gould  _et al._ (2022) C. A. Gould, K. R. McClain, D. Reta, J. G. C. Kragskow, D. A. Marchiori, E. Lachman, E.-S. Choi, J. G. Analytis, R. D. Britt, N. F. Chilton, B. G. Harvey, and J. R. Long, Science 375, 198 (2022).
  * Askerka  _et al._ (2017) M. Askerka, G. W. Brudvig, and V. S. Batista, Accounts of Chemical Research 50, 41 (2017).
  * Raymond and Blankenship (2008) J. Raymond and R. E. Blankenship, Coordination Chemistry Reviews 252, 377 (2008), the Role of Manganese in Photosystem II.
  * Shavitt and Bartlett (2009) I. Shavitt and R. Bartlett,  _Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory_, Cambridge Molecular Science (Cambridge University Press, 2009).
  * 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, Nat. Commun. 5, 4213 (2014).
  * McClean  _et al._ (2016) J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, 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, 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, Phys. Rev. A 95, 020501 (2017).
  * Kandala  _et al._ (2017) A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Nature 549, 242 (2017).
  * Kandala  _et al._ (2019) A. Kandala, K. Temme, A. D. Corcoles, A. Mezzacapo, J. M. Chow, and J. M. Gambetta, Nature 567, 491 (2019).
  * 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, Phys. Rev. X 8, 011021 (2018).
  * Huggins  _et al._ (2020) W. J. Huggins, J. Lee, U. Baek, B. O’Gorman, and K. B. Whaley, New J. Phys. 22, 073009 (2020).
  * Cao  _et al._ (2019) Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. Sawaya, _et al._ , Chemical reviews 119, 10856 (2019).
  * Grimsley  _et al._ (2019a) H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, Nature communications 10, 1 (2019a).
  * Grimsley  _et al._ (2019b) H. R. Grimsley, D. Claudino, S. E. Economou, E. Barnes, and N. J. Mayhall, Journal of chemical theory and computation 16, 1 (2019b).
  * Verteletskyi  _et al._ (2020) V. Verteletskyi, T.-C. Yen, and A. F. Izmaylov, The Journal of chemical physics 152, 124114 (2020).
  * McArdle  _et al._ (2019) S. McArdle, T. Jones, S. Endo, Y. Li, S. C. Benjamin, and X. Yuan, npj Quantum Inf. 5, 1 (2019).
  * McArdle  _et al._ (2020) S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Reviews of Modern Physics 92, 015003 (2020).
  * Tilly  _et al._ (2021) J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth, _et al._ , arXiv preprint arXiv:2111.05176 (2021).
  * Farhi  _et al._ (2014) E. Farhi, J. Goldstone, and S. Gutmann, arXiv preprint arXiv:1411.4028 (2014).
  * Stein  _et al._ (2022) S. Stein, N. Wiebe, Y. Ding, P. Bo, K. Kowalski, N. Baker, J. Ang, and A. Li, in  _Proceedings of the 49th Annual International Symposium on Computer Architecture_ (2022) pp. 59–71.
  * 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, arXiv preprint arXiv:2101.08448 (2021).
  * Albash and Lidar (2018) T. Albash and D. A. Lidar, Rev. Mod. Phys. 90, 015002 (2018).
  * Aaronson and Arkhipov (2011) S. Aaronson and A. Arkhipov, in  _Proceedings of the Forty-Third Annual ACM Symposium on Theory of Computing_, STOC ’11 (Association for Computing Machinery, New York, NY, USA, 2011) p. 333–342.
  * Trabesinger (2012) A. Trabesinger, Nat. Phys. 8, 263 (2012).
  * Georgescu  _et al._ (2014) I. M. Georgescu, S. Ashhab, and F. Nori, Rev. Mod. Phys. 86, 153 (2014).
  * 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, Nat. Phys. 16, 205 (2020).
  * Parrish and McMahon (2019) R. M. Parrish and P. L. McMahon, arXiv preprint arXiv:1909.08925 (2019).
  * Kyriienko (2020) O. Kyriienko, npj Quantum Inf. 6, 1 (2020).
  * McClean  _et al._ (2017) J. R. McClean, M. E. Kimchi-Schwartz, J. Carter, and W. A. de Jong, Phys. Rev. A 95, 042308 (2017).
  * Takeshita  _et al._ (2020) T. Takeshita, N. C. Rubin, Z. Jiang, E. Lee, R. Babbush, and J. R. McClean, Phys. Rev. X 10, 011004 (2020).
  * McClean  _et al._ (2020) J. McClean, Z. Jiang, N. Rubin, R. Babbush, and H. Neven, Nat. Commun. 11, 636 (2020).
  * Suchsland  _et al._ (2021) P. Suchsland, F. Tacchino, M. H. Fischer, T. Neupert, P. K. Barkoutsos, and I. Tavernelli, Quantum 5, 492 (2021).
  * Tkachenko  _et al._ (2023) N. V. Tkachenko, L. Cincio, A. I. Boldyrev, S. Tretiak, P. A. Dub, and Y. Zhang, arXiv , 2204.10741 (2023).
  * Seki and Yunoki (2021) K. Seki and S. Yunoki, PRX Quantum 2, 010333 (2021).
  * Kowalski and Peng (2020) K. Kowalski and B. Peng, J. Chem. Phys. 153, 201102 (2020).
  * Vallury  _et al._ (2020) H. J. Vallury, M. A. Jones, C. D. Hill, and L. C. L. Hollenberg, Quantum 4, 373 (2020).
  * Aulicino  _et al._ (2022) J. C. Aulicino, T. Keen, and B. Peng, International Journal of Quantum Chemistry 122, e26853 (2022).
  * Ollitrault  _et al._ (2020) P. J. Ollitrault, A. Kandala, C.-F. Chen, P. K. Barkoutsos, A. Mezzacapo, M. Pistoia, S. Sheldon, S. Woerner, J. M. Gambetta, and I. Tavernelli, Phys. Rev. Res. 2, 043140 (2020).
  * Hill and Wheeler (1953) D. L. Hill and J. A. Wheeler, Physical Review 89, 1102 (1953).
  * Griffin and Wheeler (1957) J. J. Griffin and J. A. Wheeler, Phys. Rev. 108, 311 (1957).
  * Rodríguez-Guzmán  _et al._ (2002) R. Rodríguez-Guzmán, J. Egido, and L. Robledo, Nuclear Physics A 709, 201 (2002).
  * Bender  _et al._ (2003) M. Bender, P.-H. Heenen, and P.-G. Reinhard, Reviews of Modern Physics 75, 121 (2003).
  * Ring and Schuck (2004) P. Ring and P. Schuck,  _The nuclear many-body problem_ (Springer Science & Business Media, 2004).
  * Yao  _et al._ (2010) J. Yao, J. Meng, P. Ring, and D. Vretenar, Physical Review C 81, 044311 (2010).
  * Egido (2016) J. L. Egido, Physica Scripta 91, 073003 (2016).
  * Hizawa  _et al._ (2021) N. Hizawa, K. Hagino, and K. Yoshida, Physical Review C 103, 034313 (2021).
  * Zheng  _et al._ (2023) M. Zheng, B. Peng, N. Wiebe, A. Li, X. Yang, and K. Kowalski, Phys. Rev. Res. 5, 023200 (2023).
  * Fukutome (1981) H. Fukutome, Progress of Theoretical Physics 65, 809 (1981).
  * Evangelista  _et al._ (2019) F. A. Evangelista, G. K.-L. Chan, and G. E. Scuseria, J. Chem. Phys. 151, 244112 (2019).
  * Lee  _et al._ (2019) J. Lee, W. J. Huggins, M. Head-Gordon, and K. B. Whaley, Journal of Chemical Theory and Computation 15, 311 (2019), https://doi.org/10.1021/acs.jctc.8b01004 .
  * Baek  _et al._ (2023) U. Baek, D. Hait, J. Shee, O. Leimkuhler, W. J. Huggins, T. F. Stetina, M. Head-Gordon, and K. B. Whaley, PRX Quantum 4, 030307 (2023).
  * Jankowski and Paldus (1980) K. Jankowski and J. Paldus, International Journal of Quantum Chemistry 18, 1243 (1980).
  * Grimsley  _et al._ (2019c) H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, Nature Communications 10 (2019c), 10.1038/s41467-019-10988-2.
  * Zhang  _et al._ (2021) K. Zhang, M.-H. Hsieh, L. Liu, and D. Tao, Physical Review Research 3 (2021), 10.1103/physrevresearch.3.043095.
  * Qiskit contributors (2023) Qiskit contributors, “Qiskit: An open-source framework for quantum computing,”  (2023), qiskit version 0.45.0.
  * Tang  _et al._ (2021) H. L. Tang, V. Shkolnikov, G. S. Barron, H. R. Grimsley, N. J. Mayhall, E. Barnes, and S. E. Economou, PRX Quantum 2, 020310 (2021).
  * Kowalski and Bauman (2023) K. Kowalski and N. P. Bauman, Phys. Rev. Lett. 131, 200601 (2023).
  * Zheng  _et al._ (2022) M. Zheng, B. Peng, A. Li, X. Yang, and K. Kowalski, “QuGCM,” https://github.com/pnnl/QuGCM (2022).
  * Hashim  _et al._ (2021) A. Hashim, R. Rines, V. Omole, R. K. Naik, J. M. Kreikebaum, D. I. Santiago, F. T. Chong, I. Siddiqi, and P. Gokhale, arXiv preprint (2021).
  * Yordanov  _et al._ (2020) Y. S. Yordanov, D. R. M. Arvidsson-Shukur, and C. H. W. Barnes, Phys. Rev. A 102, 062612 (2020).
  * Grimsley  _et al._ (2022) H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, “adapt-vqe,” https://github.com/mayhallgroup/adapt-vqe (2022).
  * Virtanen  _et al._ (2020) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, Nature Methods 17, 261 (2020), scipy version 1.11.3.