# Towards Efficient Quantum Computation of Molecular Ground State Energies
using Bayesian Optimization with Priors over Surface Topology

Farshud Sorourifar  Department of Chemical and Biomolecular Engineering, The
Ohio State University, Columbus, OH, USA  Mohamed Taha Rouabah  Constantine
Quantum Technologies, University Constantine 1 Fréres Mentouri, Constantine,
Algeria  Laboratoire de Physique Mathématique et Subatomique, University
Constantine 1 Fréres Mentouri, Constantine, Algeria  Nacer Eddine Belaloui
Constantine Quantum Technologies, University Constantine 1 Fréres Mentouri,
Constantine, Algeria  Laboratoire de Physique Mathématique et Subatomique,
University Constantine 1 Fréres Mentouri, Constantine, Algeria  Mohamed
Messaoud Louamri  Constantine Quantum Technologies, University Constantine 1
Fréres Mentouri, Constantine, Algeria  Theoretical Physics Laboratory,
University of Science and Technology Houari Boumediene, Algiers, Algeria
Diana Chamaki  Department of Chemistry, Columbia University, New York, NY, USA
Erik J. Gustafson  USRA Research Institute for Advanced Computer Science,
Mountain View, CA, USA  Norm M. Tubman  NASA Ames Research Center, Moffett
Field, CA, USA  Joel A. Paulson  Department of Chemical and Biomolecular
Engineering, The Ohio State University, Columbus, OH, USA  David E. Bernal
Neira  USRA Research Institute for Advanced Computer Science, Mountain View,
CA, USA  NASA Ames Research Center, Moffett Field, CA, USA  Davidson School of
Chemical Engineering, Purdue University, West Lafayette, IN, USA

###### Abstract

Variational Quantum Eigensolvers (VQEs) represent a promising approach to
computing molecular ground states and energies on modern quantum computers.
These approaches use a classical computer to optimize the parameters of a
trial wave function, while the quantum computer simulates the energy by
preparing and measuring a set of bitstring observations, referred to as shots,
over which an expected value is computed. Although more shots improve the
accuracy of the expected ground state, it also increases the simulation cost.
Hence, we propose modifications to the standard Bayesian optimization
algorithm to leverage few-shot circuit observations to solve VQEs with fewer
quantum resources. We demonstrate the effectiveness of our proposed approach,
Bayesian optimization with priors on surface topology (BOPT), by comparing
optimizers for molecular systems and demonstrate how current quantum hardware
can aid in finding ground state energies.

Keywords: Variational Quantum Eigensolver, Bayesian Optimization, Gaussian
Processes

##  1 Introduction

Solving electronic structure problems for molecular systems is a cornerstone
of modern chemical engineering research. These problems involve solving
Hamiltonian eigenvalue systems, which are used to define the ground-state
energy, the lowest possible configuration of a molecular system, which
directly influences the behavior of the system [1, 2]. These energies are used
to predict molecular properties [3], understand chemical reactions [4, 5],
control molecular systems [6, 7], as well as design pharmaceuticals [8],
polymers [9, 10], materials [11, 12, 13, 14, 15], and catalysts [16, 17, 18].
Consequently, these calculations affect the design and optimization of
industrial processes, with applications spanning healthcare, energy, and the
environment [19].

Calculating the ground state involves solving the Schrödinger equation, which
can only be solved exactly for a few selected cases. The interactions between
increasingly many electrons result in an exponential computational complexity,
making exact solutions practically intractable for large systems [20, 21].
Although approximate methods are used to compute the ground state [22, 23,
24], the scaling issue with many-body calculations persists, requiring further
approximations for larger systems [25]. Furthermore, these approximations are
known to erode the quality of the solutions to systems with strong electron-
electron interactions [26].

Solving many-body systems is a problem class where quantum computers are
expected to excel, making them a promising tool for efficiently computing
exact ground state energies [27, 28, 29, 30, 31], which also includes hardware
demonstrations [32, 29, 33, 34, 35, 36, 37, 38, 39]. An algorithm designed
explicitly for this task is known as the variational quantum eigensolver
(VQE), which has been used to great effect in many problems [40, 29, 41, 42,
43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54]. This approach leverages the
quantum computers’ capacity for efficiently simulating many-body systems,
which can be used to determine the ground state energy, $E_{0}$, from the
variational principle of quantum mechanics,

$E_{0}\leq\frac{\bra{\Psi(\theta)}\mathcal{H}\ket{\Psi(\theta)}}{\bra{\Psi(\theta)}\ket{\Psi(\theta)}}.$

Here, $\Psi(\theta)$ is a trial wave function of the system parameterized by
$\theta$, and $\mathcal{H}$ is the Hamiltonian acting on the system. Note that
in practice, the trial wave functions may be normalized so that
$\bra{\Psi(\theta)}\ket{\Psi(\theta)}=1$. The trial wave function represents
the state of the electronic configuration, whereas the Hamiltonian defines the
system’s energy. The VQE determines $E_{0}$ by formulating equation (1) as an
optimization problem,

$\min_{\theta\in\mathbb{R}^{d}}y(\theta)=\frac{1}{S}\sum^{S}\frac{\langle\Psi(\theta)|\mathcal{H}|\Psi(\theta)\rangle}{\langle\Psi(\theta)|\Psi(\theta)\rangle}$

which is solved on a non-quantum or _classical_ computer using optimization by
simulating $y(\theta)$ on the quantum computer for selected $\theta$. The
sample-based expectation value is calculated over $S$ repeated observations of
the final quantum state, referred to as shots. At the same time, this
theoretically allows for arbitrary precision, although a greater number of
measurements naturally incurs additional costs.

Recently, Bayesian optimization (BayesOpt) has been investigated as a means of
solving Equation (1) using VQE and has shown competitive performance against
other commonly used optimization routines [55, 56, 57, 58]. BayesOpt is a
general class of model-based derivative-free optimization (DFO) algorithms.
Model-based DFO methods use a surrogate model to determine a sampling
strategy, in contrast to model-free strategies (e.g., genetic algorithms),
which are less sample-efficient [59]. What differentiates BayesOpt from other
DFO methods is the use of a statistical surrogate to model a mean-based
prediction and the epistemic uncertainty (uncertainty caused by limited data)
to determine the most efficient sampling strategy. Gaussian processes (GPs)
are frequently used due to their excellent capacity for uncertainty
quantification, assuming the kernel function, which defines the GP model, has
been appropriately selected [60, 61, 62, 63, 64]. The mean and uncertainty
functions derived from the surrogate model are used to trade off between
exploring the parameter space (i.e., sampling in regions of greater
uncertainty) and exploiting the most promising regions (e.g., sampling near
the best-observed parameters). This exploration and exploitation trade-off
allows the global optimum to be found with limited evaluations.

BayesOpt for VQE problems has been the subject of growing research interest,
with most work focused on using standard implementations with limited
modifications [55, 56, 57, 58], or use BayesOpt as a component of a larger
optimization framework [65, 66]. Several trust-region-based BO methods have
also been investigated [67, 68], which provide convergence stability and allow
for parallel sampling. In [69], the authors propose a parallel BayesOpt
strategy for solving related problems, such as solving the VQE for an ensemble
of molecular geometries, by sharing sample data across the different problems,
and [70] builds on this work by investigating different information sharing
schema. The only work of which we are aware to propose modifications to the GP
model specialized for VQE is in [71] where the authors develop a novel
periodic kernel to construct the GP model, in addition to a novel acquisition
function for their trust-region based BaysOpt framework, inspired by the
Nakanishi-Fuji-Todo method [72].

In this work, we build on our proposed BayesOpt framework [73], which
constructs a prior on surface topology from noisy low-shot observations and
then optimizes for high-shot observations (i.e., the number of shots required
to achieve a prescribed accuracy) using a residual GP model. The prior on
surface topology uses a data-driven surrogate model as an initial belief of
$y(\theta)$ in lieu of the traditional zero-mean prior. Additionally, we find
that encoding information about the periodicity of the rotation gate
parameters translates into faster convergence. We modify the prior on surface
topology to accommodate larger datasets and consider the practical usage of
quantum computers. Specifically, we investigate the algorithm’s performance
for three noise models constructed from existing IBM quantum computers
(computational results) and on physical quantum hardware (experimental
results). While current quantum hardware tends to provide biased estimates of
the ground state, BayesOpt can still find optimal circuit parameters,
suggesting quantum computers may be useful for determining optimal circuit
parameters, i.e.,
$\theta^{\star}\in\operatorname*{argmin}_{\theta\in\mathbb{R}^{d}}y(\theta)$
in Equation (1), and computing the energy classically. The rest of the paper
is organized as follows: we first provide a detailed background on quantum
computing, variational quantum eigensolvers, and Bayesian optimization in
Section 2. We formally introduce the problem in Section 3 and present our
methods in Section 4. In Section 5, we present the VQE convergence results for
our method and several benchmarking methods for two molecular systems, with a
particular focus on how the algorithms perform on high-fidelity noise models
and when using noisy physical quantum computers. In addition, that section
discusses strategies for efficiently using currently available quantum
hardware to compute molecular ground-state energies relative to those of
purely classical methods. Lastly, we provide concluding remarks and directions
for future work in Section 6.

Figure 1: Bayesian Optimization with Priors on Topology (BOPT) Summary. Priors
on the quantum circuit parameter’s $2\pi$ periodicity and surface topology, in
the form of a periodic kernel and surrogate models, are used to initialize the
Bayesian Optimizations loop. At each iteration, the Gaussian process is
conditioned on the available circuit observations, which can be used to
construct the acquisition function. Optimizing the acquisition provides the
subsequent circuit parameters, which are used to query the circuit. With the
new observation, the loop is repeated until the sample budget is exhausted.

##  2 Background and Preliminaries

###  2.1 Quantum computing

Quantum information processing is a strategy to process information in
accordance with the principles of quantum mechanics. One area of quantum
information processing, quantum simulation, is used to solve problems
formulated as quantum systems. Many problem classes can be modified to fit
these criteria, such as factoring large numbers (for cryptographic
applications), simulating quantum systems, and solving optimization problems
[74]. Although it is important to note that quantum simulation is not
necessarily superior to classical methods, there are specific cases where it
is expected to excel, such as many-body problems [28].

Quantum computers, or quantum computing units (QPU), exploit quantum
mechanical phenomena to process information. This information processing is
performed physically using systems that exhibit multi/two-level quantum
behavior, such as superconducting qubits [75, 76] present in IBM quantum
devices. These two-qubit states are known as _qubits_. For quantum computers,
the computations are interpreted as simulations. These simulations are
represented by evolving a unit vector in the Hilbert space, which represents
the quantum state, based on the rules of quantum mechanics. Although
computationally expensive, these simulations can also be performed on
classical computers, where quantum effects are implemented by exponentially
growing states with the number of qubits. Conversely, leveraging quantum
mechanical phenomena directly in a quantum computer has the potential to scale
much more favorably since the computation is replaced by a controlled
simulation of the physical system in the computer.

Specifically, the quantum phenomena that are advantageous for quantum
simulations and improved information processing are superposition and
entanglement. Unlike classical bits, which only exist in one of two states (0
or 1), qubits, the fundamental unit of quantum information, can be in a
coherent superposition of both states simultaneously. The superposition of
states allows more information to be encoded in the base computing unit.
Entanglement is a phenomenon where the state of one qubit is intrinsically
correlated with the state of another qubit so that states cannot be factored
into products of individual states. In quantum computing, this allows
information that grows exponentially with the number of qubits to be processed
seamlessly. In addition, many-body systems often exhibit entanglement among
their constituent particles, making it a useful phenomenon for quantum
simulation. For a detailed discussion on the topic, we point the interested
reader to [77].

A well-defined initial state for each orbital in a molecular quantum system is
prepared to simulate it. These states are propagated through a quantum circuit
and are manipulated by quantum gates. The gates perform unitary
transformations, such as parameterized rotations and controlled gates, to
modify the quantum state. The set of gates or _circuit_ is designed to
represent the trial wave function, where the process of transforming the
initial state is known as state preparation. After the circuit execution, when
the final state of the circuit is obtained, it is _measured_. In this process,
the coherent state is projected or _collapsed_ onto a classical state
represented by a string of bits, e.g., 00110101. Each single execution is
known as a _shot_ , which is then mapped to an observable via an operator,
which in this case is the energy computed through the Hamiltonian. Lastly, a
sample-averaged expectation is taken over the collection of observed values.
On quantum hardware, running a shot is fast, and problems with 27 qubits using
thousands of shots can be simulated in seconds [78].

Although quantum computing is a promising and exciting technology for
simulating quantum systems, we emphasize that current computers are noisy
intermediate-scale quantum (NISQ) devices with open technical challenges [79].
The largest current quantum computers still have too few and noisy qubits
($\sim 500$) 111If not for current levels of noise, $\sim$ 100 qubits would
provide sufficient computing resources for the utility of quantum computing
[80]. However, NISQ error correction schemes often call for additional qubits
[81, 82]. , both limiting the size of molecular systems that can be simulated
and requiring the use of hybrid quantum-classical algorithms, like VQE, that
offload some of the computational burdens onto classical devices. The quantum
states of qubits are also fragile and prone to decoherence because of
environmental disturbances, such as thermal fluctuations. As a result, the
quantum states in NISQ devices can only be maintained for short periods before
decohering, limiting the complexity of circuits that can be implemented and of
problems that can be solved. Additionally, noise is introduced due to hardware
limitations, such as inaccuracies when implementing gate protocols. While
fully characterizing the quantum state is subject to the inherent uncertainty
of quantum mechanics, this can be addressed by using a large number of shots
to observe the circuit, which incurs additional costs. Thus, the practical
value of modern quantum computers depends on creative solutions for navigating
noisy observations and minimizing the computational load assigned to the
quantum computer.

Remark: While there is a myriad of quantum computing technologies that are
continuously being improved, we focus on a topical treatment of the subject.
The goal of this section is to provide high-level connections between
components in common hardware architectures and explain how they relate to the
mathematical concept of quantum simulation. While we strive to remain
technology agnostic, for brevity, we resort to gate-based models of quantum
computers to describe the process of simulating a trial wave function. We note
that the gate-based QPU architectures are one of the most widely studied
quantum computers and constitute the computational architecture used in our
hardware simulations.

###  2.2 Variational Quantum Eigensolvers

The VQE algorithm uses a hybrid quantum-classical approach to determine the
ground state of molecular systems by optimizing the parameters of a trial wave
function $\Psi(\theta)$, which is simulated in the quantum device (i.e.,
either a quantum computer or classical quantum simulator). First, a reference
state $\Psi_{\text{ref}}$ is selected to initialize the trial wave function,
such as the Hartree-Fock reference state. The principal challenge in VQE is
known as state preparation, which is the process of designing a trial wave
function to minimize the expected energy. This step can be separated into two
parts: the selection of a circuit structure known as _ansatz_ that represents
the trial wave function and the selection of parameters of the ansatz to drive
the system from the reference state to the trial state exhibiting the lowest
energy configuration.

Although the ansatz selection is still an open question, we present two that
are frequently encountered: those from the class of hardware efficient ansätze
(HEA) [32, 83, 84, 85, 86], and the unitary coupled cluster with singles and
doubles (UCCSD) ansatz [87, 88, 89]. In addition, the circuit depth $D$, which
is the number of repeating layers, is often treated as part of the ansatz.
Deeper circuits result in increased representational capabilities, ensuring
that $\Psi(\theta)$ can model the ground state but also increase the circuit
running time and the likelihood of decoherence.

UCCSD is a chemistry-inspired ansatz rooted in classical coupled cluster
theory, a well-established method for computing electronic structures of
molecules. The method is designed to capture the most important electronic
correlations in the molecular system by applying single- and double-excitation
operators to the reference state. Formally, the USSCD ansatz is defined as,

$\ket{\Psi({\theta})}=\big{[}\Pi_{i=1}^{D}e^{\hat{T}({\theta})-\hat{T}({\theta})^{{\dagger}}}\big{]}\ket{\Psi_{ref}},$
(1)

where $\theta$ represents the variational parameters that need to be
optimized, and
$\hat{T}({\theta})=\hat{T}_{1}({\theta})+\hat{T}_{2}(\vec{\theta})$ are the
aggregated single and double excitation operators.

The HEA, as its name suggests, is a more flexible class of ansätze designed to
use fewer gates and use NISQ-era QPUs more efficiently. As a result, they are
not limited to chemical systems. While these features sound appealing, it is
unclear whether tailored ansätze can potentially provide better utility with a
larger number of gates, e.g., by requiring fewer layers. A generic HEA can be
defined as

$\ket{\Psi(\vec{\theta})}=\big{[}\Pi_{i=1}^{D}\hat{\mathrm{U}}_{r}(\theta_{i})\times\hat{\mathrm{U}}_{\text{ent}}\big{]}\ket{\Psi_{\rm
ref}},$ (2)

where $\hat{\mathrm{U}}_{\text{ent}}$ are the entangling gates, and
$\hat{\mathrm{U}}_{r}(\theta_{i})$ are the single qubit gates. Here,
$\theta_{i}\in\mathbb{R}^{n_{g}}$ is the parameter vector for each layer,
where $n_{g}$ is the number of gates in each layer, such that $d=n_{g}D$. We
present a schematic overview of a quantum circuit using the HEA in Figure 2.
For more details on VQE algorithms and best practices, see [27].

Figure 2: (Left) Hartree-Fock initialization (shown in blue boxes with “X”)
HEA for $H_{2}$ (4 qubit system) using $D$ layers composed of $R_{y}$ single
qubit gates and CNOT entanglement gates acting on all possible qubit pairs.
The dashed box denotes the layer. (Right) $H_{2}$ molecular orbital diagram
showing the orbital qubit assignments.

###  2.3 Bayesian Optimization

As mentioned previously, BayesOpt is a model-based optimization algorithm
geared toward expensive black-box functions. It consists of two key
components: a statistical surrogate model and an acquisition function, the
latter of which defines how to use the model to identify subsequent sample
points. The algorithm proceeds sequentially. The GP is fit to available data,
and the predicted mean and variance are used to construct an acquisition
function. The acquisition function is maximized to determine the following
sample location. The black-box function is evaluated at the selected sample
point, which is then added to the data set, and the process is repeated until
a specified sample budget $B$ is exhausted. In practice, a fraction of the
initial budget $B_{\text{init}}$ is typically allocated to pseudo-random
sampling to generate an initial data set when little to no prior data are
available.

We begin by introducing the mathematical underpinning of GPs and their use in
BayesOpt; interested readers are referred to [60] for further details. GPs
generalize the notion of a multivariate Gaussian to function space,
representing an uncountable set of functions that fit the data. Specifically,
it constructs a posterior probability distribution over function space that
can be conditioned on the observed data, in accordance with Bayes’s rule,

$P\left(f(\theta)|\mathcal{D}^{\star}_{n}\right)\propto
P\left(\mathcal{D}^{\star}_{n}|f(\theta)\right)P\left(f(\theta)\right).$ (3)

Here, $P(f(\theta))$ denotes the prior probability of observing $f$ for a
given $\theta$, and $P(\mathcal{D}^{\star}_{n}|f(\theta))$ denotes the
conditional probability of observing the data
$\mathcal{D}^{\star}_{n}=\\{\left(\theta_{i},f(\theta_{i})\right)\\}_{i=0}^{n}$
(also known as the likelihood function). The key assumption in a GP, denoted
by $GP(\mu(\theta),k(\theta,\theta^{\prime}))$, is that the function values at
any subset of inputs are jointly Gaussian such that they are fully defined by
their mean function $\mu:\mathbb{R}^{d}\to\mathbb{R}$ and covariance or
_kernel_ function $k:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}_{\geq
0}$.

In general, $f(\theta)$ cannot be directly observed, and instead, the
measurement process introduces a source of noise that can be modeled as
$y=f(\theta)+\epsilon$, where $\epsilon\sim\mathcal{N}(0,\sigma^{2})$. Given a
noisy dataset $\mathcal{D}_{n}=\\{(\theta_{i},y_{i})\\}_{i=0}^{n}$, the
posterior probability of $f(\theta)$ can be modeled jointly with the noisy
observations $y_{1:n}$ as

$\left[\begin{array}[]{l}y_{1:n}\\\
f(\theta)\end{array}\right]\sim\mathcal{N}\left(\left[\begin{array}[]{c}{\mu_{0}}(\theta_{1:n})\\\
\mu_{0}(\theta)\end{array}\right],\left[\begin{array}[]{cc}\mathbf{K}_{n}+\sigma^{2}I_{n}&\mathbf{k}_{n}(\theta)\\\
\mathbf{k}_{n}(\theta)^{\top}&k_{0}(\theta,\theta)\end{array}\right]\right).$
(4)

Here, $\mathbf{k}_{n}(\theta)\in\mathbb{R}^{n}$ is a vector of covariance
terms between the observed inputs $\theta_{1:n}$ and the test point $\theta$,
i.e., $[\mathbf{k}_{n}(\theta)]_{i}=k(\theta_{i},\theta)$, and
$\mathbf{K}_{n}\in\mathbb{R}^{n\times n}$ is a symmetric covariance matrix
with elements $[\mathbf{K}_{n}]_{i,j}=k(\theta_{i},\theta_{j})$ for all
$i,j=1,\ldots,n$. Due to the properties of joint Gaussian random variables,
the posterior distribution $p(f(\theta)|\mathcal{D}_{n},\theta)$ given all
available observations, is normally distributed with the following posterior
mean and variance as a function of the test point

$\displaystyle\mu_{n}(\theta)$
$\displaystyle=\mu_{0}(\theta)+\mathbf{k}_{n}(\theta)^{\top}(\mathbf{K}_{n}+\sigma^{2}I_{n})(y_{1:n}-{\mu}_{0}(\theta_{1:n})),$
(5a) $\displaystyle\sigma^{2}_{n}(\theta)$
$\displaystyle=k_{0}(\theta,\theta)-\mathbf{k}_{n}(\theta)^{\top}(\mathbf{K}_{n}+\sigma^{2}I_{n})^{-1}\mathbf{k}_{n}(\theta).$
(5b)

Note that here, $\mu_{0}$ and $k_{0}$ refer to the prior mean and kernel
functions, respectively. The matrix inverse
$(\mathbf{K}_{n}+\sigma^{2}I_{n})^{-1}$ in equation (5) is often considered a
limiting factor for GP models - the computational complexity scales cubically
with the number of data points, making inference computationally challenging
for large data sets.

A zero-mean prior ($\mu_{0}(\theta)=0$) is generally a reasonable assumption
when normalizing the data to be zero-centered. The kernel function is
typically selected based on a prior belief about the degree of smoothness of
the underlying function. A common choice of the kernel class is that from the
Matérn class,

$k_{\text{Mat\'{e}rn}}\left(\theta,\theta^{\prime}\right)=\sigma_{k}^{2}\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\sqrt{2\nu}{\sum_{i=1}^{d}\rho_{i}\left(\theta_{i}-\theta_{i}^{\prime}\right)}\right)^{\nu}K_{\nu}\left(\sqrt{2\nu}{\sum_{i=1}^{d}\rho_{i}\left(\theta_{i}-\theta_{i}^{\prime}\right)}\right),$
(6)

where $\Gamma$ is the gamma function and $K_{\nu}$ is the modified Bessel
function of the second kind. The parameters $\rho_{i}$ for all $i=1,\ldots,d$
are inverse squared length scale parameters, which models the response
sensitivity for each dimension of $\theta\in\mathbb{R}^{d}$, and
$\sigma_{k}^{2}$ is a scale parameter meant to capture expected variation in
$f$. The parameter $\nu$ controls the smoothness of the GP model, with an
integer value $\nu$ resulting in a model that is $\nu-1$ time differentiable.
In the limit $\nu\rightarrow\infty$, the well-known squared exponential (or
radial basis) function is recovered.

Another much less commonly used kernel is the periodic kernel. Considering
that the parameters used in quantum circuits are periodic, and this kernel is
a natural choice. The periodic kernel is defined as

$k_{\text{Periodic}}\left(\mathbf{\theta},\mathbf{\theta}^{\prime}\right)=\sigma_{k}^{2}\exp\left(-2\sum_{i}{\sin^{2}\left(\frac{\pi}{p}\left(\theta_{i}-\theta_{i}^{\prime}\right)\right)\rho_{i}}\right).$
(7)

While $p$ is typically an unknown parameter, in VQE, setting $p=2\pi$ for
$\theta_{i}\in[0,2\pi]$ can be used to enforce periodicity. Furthermore, we
note that this kernel can be considered another topological prior, as it
represents the relationship between the cost function and the circuit
parameters, which is a sine curve with period $2\pi$ [72].

Part of selecting the kernel also requires choosing the hyperparameters
$\gamma=\left\\{\rho_{1},\ldots,\rho_{D},\sigma_{k}^{2}\right\\}$. This
selection is made after each new observation by minimizing the negative
marginal log-likelihood,

$\gamma^{\star}=\operatorname*{argmin}_{\gamma}\underbrace{\frac{1}{2}(\mathbf{y}_{n}-\boldsymbol{\mu}_{n})^{\top}\Sigma_{\mathbf{y}_{n}}^{-1}(\mathbf{y}_{n}-\boldsymbol{\mu}_{n})}_{\text{Data
fit term
}}+\underbrace{\frac{1}{2}\log\left(\left|\Sigma_{\mathbf{y}_{n}}\right|\right)}_{\text{Complexity
penalty }}+\frac{n_{y}}{2}\log(2\pi).$ (8)

where $\mathbf{y}_{n}=y_{1:n}$, $\boldsymbol{\mu}_{n}=\mu(\theta_{1:n})$, and
$\Sigma_{\mathbf{y}_{n}}=\mathbf{K}_{n}+\hat{\sigma}_{f}^{2}I_{n}$. The
negative marginal log-likelihood is a loss function that provides some
resilience to overfitting. The data-fit term encourages the model to match the
available data, while the complexity term provides a regularization effect
that encourages a more simple model.

For a given dataset, we can use equation (8) to fit a GP to a dataset, which
provides a mean and covariance function defined by equations (5). Now, we turn
to the acquisition function, which is used to select subsequent samples in
search of optimal circuit parameters. One clear example of how (5) can be used
to select $\theta_{n+1}$ to balance competing objectives (exploration versus
exploitation) is the lower confidence bound (LCB) acquisition function

$\alpha_{\text{LCB},n}(\theta)=\underbrace{\mu_{n}(\theta)}_{\text{Exploitation}}-\underbrace{\sqrt{\beta}\sigma_{n}(\theta)}_{\text{Exploration}},$
(9)

which determines the new sample point by minimizing
$\alpha_{\text{LCB},n}(\theta)$ over some compact candidate set
$\Theta\subseteq\mathbb{R}^{d}$

$\theta_{n+1}=\operatorname*{argmin}_{\theta\in\Theta}\alpha_{\text{LCB},n}(\theta),$
(10)

which can be solved with standard gradient-based optimization. Here, the mean
function promotes the exploitation of regions that are known to be good, while
the covariance function, scaled by the confidence interval parameter
$\sqrt{\beta}$, promotes exploitation. Together, the acquisition function
assigns optimistic values to each candidate point.

Furthermore, we introduce the more commonly used expected improvement
acquisition function, which measures the amount by which the current objective
is expected to improve over the current best or _incumbent_ solution, $\eta$.
Formally, the expected improvement criterion for minimizing an objective
function is defined as

$\alpha_{\text{EI},n}(\theta)=\mathbb{E}_{n}\big{[}\max\left(\eta-f(\theta),0\right)\big{]}.$
(11)

The EI acquisition function is commonly used since an analytic solution can be
derived for GP models [62]

$\displaystyle\alpha_{\text{EI},n}(\theta)=\sigma_{n}(\theta)\left(z_{n}(\theta)\Phi(z_{n}(\theta)\right)+\phi\left(z_{n}(\theta)\right),$
(12)

where $z_{n}(\theta)=(\eta-\mu_{n}(\theta))/\sigma_{n}(\theta)$, $\Phi$ is the
standard Gaussian cumulative density, and $\phi$ is the probability density.
While $\eta$ is often set to the best-observed value
$\eta=\min_{i\in\\{1,\ldots,n\\}}y_{i}$, it has been discussed in [90] that
this choice causes issues in a noisy setting. Since noise is an important
characteristic for the objective functions considered in this work, we instead
use the best posterior mean value as our incumbent, i.e.,
$\eta=\min_{\theta\in\Theta}\mu_{n}(\theta)$.

##  3 Problem Formulation

Before introducing the proposed approach, we formally present the VQE problem.
First, let $f(\theta)$ represent a theoretical circuit evaluated at a given
vector of parameters $\theta$. The circuit $f(\theta)$ is theoretical in the
sense that it represents a noise-free simulation with an exact expectation
when observing the state (i.e., in the infinite shot limit). We assume that
querying the circuit produces observations that are corrupted by the
underlying generative process

$\displaystyle
y_{\theta}=f(\theta)+\epsilon_{S}+\epsilon_{\text{QPU}}(\theta),$ (13)

where $\epsilon_{S}\sim\mathcal{N}\left(0,(\sigma_{s}/\sqrt{S})^{2}\right)$ is
an independent and identically distributed Gaussian noise term, with a
standard deviation proportional to $1/\sqrt{S}$, where $S$ is the number of
shots chosen to evaluate the circuit (i.e., the standard error in the mean
estimate), and the term
$\epsilon_{\text{QPU}},\sim\mathcal{N}\left(\mu_{\text{QPU}}(\theta),\sigma_{\text{QPU}}^{2}(\theta)\right)$,
represents the noise/errors due to hardware. Classical quantum simulators can
generally be modeled with $\epsilon_{\text{QPU}}=0$. In contrast, quantum
hardware will have an unknown and hardware-specific structure for
$\mu_{\text{QPU}}(\theta)$ and $\sigma_{\text{QPU}}^{2}(\theta)$. However,
data-driven models exist to simulate QPU performance in classical simulations.
Formally, the VQE we aim to solve is

$\theta^{\star}=\operatorname*{argmin}_{\theta\in\Theta}f(\theta)$ (14)

where $\Theta\in[0,2\pi]^{d}$ is the parameter space. When the number of shots
$S$ is sufficiently large, $f(\theta)+\epsilon_{S}\approx f(\theta)$, hence we
assume that $S$ is selected to ensure that classically simulated systems are
arbitrarily accurate, i.e., $y_{\theta}\approx f(\theta)$. This formulation
allows for solving equation (14) using observations $y_{\theta}$.

Although hardware noise mitigation is an active research area, a true noise-
resilient QPU at a reasonable scale is probably years away. As a result,
directly observing simulated energies from a QPU does not have much practical
value currently, as the errors are often greater than the accuracy required to
determine the ground state of these molecular systems, frequently called
_chemical accuracy_. Solving VQEs with classical simulation is computationally
prohibitive for systems exceeding $\sim 40$ qubits [91] ($\sim 64$ qubits and
more are possible if using approximate simulators [92, 93]) due to the state
preparation. However, a single simulation with a provided set of circuit
parameters is tractable for small molecular systems. Hence, a possible
strategy for the use of near-term quantum computers is to solve equation (14)
on a QPU and validate $\theta^{\star}$ by evaluating (1) with a classical
quantum simulation.

The form we present in equation (13), using additive noise components, is a
useful model; however, the actual noise would enter the system at a variety of
intermediate points in the simulation. For example, gate errors would take the
form of $f\left(\theta+\epsilon_{\text{gates}}(\theta)\right)$, where
$\epsilon_{\text{gates}}$ is a model of the QPU specific gate implementation
error. The errors introduced from having too few shots are due to a poor
sample-averaged approximation of the expectation taken in equation (1).
However, this form provides a reasonable approximation of the circuit that
facilitates more efficient search strategies under the BayesOpt paradigm.

##  4 Methods

###  4.1 Surface Topology Prior

The proposed framework constructs a prior on the surface topology by fitting a
surrogate model to low-shot observations and deploying BayesOpt to learn a
residual model between a low-shot surrogate and high-shot observations
$y_{\theta}$. This framework, which we refer to as a Bayesian optimization
with priors on surface topology (BOPT), is inspired by the reference models
used in [94].

First, we assume that we can generate noisy low-shot observations,

$\displaystyle\tilde{y}_{\theta}=f(\theta)+\epsilon_{s}+\epsilon_{\text{QPU}}(\theta),$
(15)

where $s<S$ is the shot count used to generate noisy observations. The costs
of querying $y_{\theta}$ and $\tilde{y}_{\theta}$ are defined as $1$ and $s/S$
respectively, i.e., the costs are linearly dependent on the number of shots,
consistent with commercial QPU cost structures [95, 96]. Let
$\mathcal{D}^{\text{low-
shot}}_{\text{init}}=\\{\theta_{i},\tilde{y}_{\theta_{i}}\\}^{m}$ be the set
of $m$ low-shot observations obtained from spending the initialization budget
$B_{\text{init}}$ uniformly over the parameter space. The prior on surface
topology is defined in terms of a general surrogate model
$\tilde{\mu}_{y}(\theta)$ for predicting $f(\theta)$ constructed from
$\mathcal{D}^{\text{low-shot}}_{\text{init}}$. Our proposed Gaussian process
with priors on surface topology (GP-TP) can then be expressed as

$\displaystyle f(\theta)\sim
GP(\tilde{\mu}_{y}(\theta),k_{\text{Periodic}}(\theta,\theta^{\prime})).$ (16)

As opposed to using a constant-mean prior function in (5) (the default
assumption), we use low-shot observations to learn a more effective prior that
shifts the high-shot observations in a nonlinear way. Note that we can
actually interpret this process as training a standard zero-mean GP on
residual data
$\\{\theta_{i},y_{\theta_{i}}-\tilde{\mu}_{y}(\theta_{i})\\}_{i=1}^{n}$. The
key innovation in BOPT is thus the use of this new type of prior, which we
show is capable of speeding up convergence in VQE applications.

We present an illustration of the motivation behind BOPT in Figure 3, where we
use additive normal noise terms, with 1 high-shot observation having the same
cost as 16 low-shot observations (i.e.,
$\epsilon_{S}\sim\mathcal{N}(0,0.1^{2})$, and
$\epsilon_{s}\sim\mathcal{N}(0,0.4^{2})$). In the leftmost plot, a standard GP
is trained on 5 high-shot samples; the model is poorly fit due to the
insufficient volume of data. On the other hand, the center-left plot shows
that 16 noisy observations can be used to learn an accurate model of the
surface topology $\tilde{\mu}_{y}(\theta)$, providing a better fit than the
standard GP model with just 20% of the budget expended. The center-right plot
shows the true residual $f(\theta)-\tilde{\mu}_{y}(\theta)$ and its
corresponding GP model using 4 high-shot observations $y_{\theta}$. The
rightmost plot represents the complete GP-TP model of the true objective
function $f(\theta)$, which uses the same budget spent on constructing the
standard GP shown.

Figure 3: BOPT illustrations: (Left) Standard GP modeling approach, fit
directly to 5 high-shot observations. (Center left) the topological prior was
constructed from 16 noisy observations. (Center right) GP model of the
residual fit using only 4 high-shot observations. (Right) the GP with
topological prior of the objective function using the same budget spent on
constructing the standard GP shown in the leftmost plot.

BOPT can take advantage of any standard BO acquisition function by simply
replacing the standard GP model with the proposed GP-TP model, whose
predictive posterior distribution is denoted by
$f(\theta)|\mathcal{D}_{n},\text{TP}\sim\mathcal{N}(\mu_{n,\text{TP}}(\theta),\sigma^{2}_{n,\text{TP}}(\theta))$.
For example, LCB can represented as follows

$\alpha_{\text{BOPT-
LCB},n}(\theta)=\mu_{n,\text{TP}}(\theta)-\sqrt{\beta}\sigma_{n,\text{TP}}(\theta),$
(17)

To extend EI to the BOPT setting, we replace the analytical solution in (12)
with the so-called noisy EI [61] approach that can more effectively handle
noise in the observations. Noisy EI performs an expectation over (11) using
values drawn from the predictive posterior distribution. In this work, we
focus predominantly on LCB. We summarize BOPT for the LCB acquisition function
in Algorithm 1.

Algorithm 1 Bayesian Optimization with Priors on Surface Topology (BOPT)
Algorithm

1:Input: Circuit $f$; high-shot count $S$; low-shot count $s$; low-shot data
$\mathcal{D}^{\text{low-shot}}_{\text{init}}$ constructed by spending
$B_{\text{init}}$ on $\tilde{y}_{\theta}$; and total budget $B$.

2:Initialize: Domain $\Theta=[0,2\pi]^{d}$; High-shot data
$\mathcal{D}^{\text{high-shot}}\leftarrow\emptyset$; and $k=B_{\text{init}}$.

3:Construct surrogate $\tilde{\mu}_{y}$ using $\mathcal{D}^{\text{low-
shot}}_{\text{init}}$

4:Initialize GP-TP model via (16)

5:while $k\leq B$ do

6: Minimize $\alpha_{\text{BOPT-LCB},k}(\theta)$ to find $\theta_{k+1}$

7: Evaluate quantum circuit $y(\theta_{k+1},s=S)$

8: Update high-shot dataset $\mathcal{D}^{\text{high-
shot}}\leftarrow\\{\theta_{i},y(\theta_{i},S)-\tilde{\mu}_{y}(\theta_{i})\\}_{i=0}^{k+1}$

9: Update posterior GP-TP model using $\mathcal{D}^{\text{high-shot}}$

10: Increment spent budget $k\leftarrow k+1$

11:end while

###  4.2 Sparse Gaussian Processes

In the BOPT framework, the algorithm is initialized by taking low-shot samples
$\mathcal{D}^{\text{low-shot}}_{\text{init}}$, which contains a relatively
large number of samples, to construct $\tilde{\mu}_{y}$ appearing in the
topological prior. While this surrogate can be constructed using any data-
driven model, it must be capable of handling large data sets containing
significant amounts of noise. In this work, we are interested in using GPs due
to their simple and effective uncertainty quantification properties,
particularly when dealing with shot noise in quantum circuits [58]. However,
the covariance matrix inverse in equation (5b) results in quadratic memory and
cubic computational complexity, making the construction and inference with a
GP-TP from large datasets a limiting factor. To improve the limit on the
number of low-shot observations that can be used, we use sparse GP models [97]
to construct $\tilde{\mu}_{y}(\theta)$.

Assume that a given dataset $\mathcal{D}^{\text{low-shot}}_{\text{init}}$
contains $m$ data points and that computational or memory requirements limit
us to constructing a GP with at most $\ell$ data points, where $\ell<m$. To
fit the GP model, a set of _inducing points_ , $\hat{\mathcal{D}}_{\ell}$ must
be learned to minimize the statistical distance between the true GP
conditioned on $\mathcal{D}^{\text{low-shot}}_{\text{init}}$ and an
approximate GP conditioned on $\hat{\mathcal{D}}_{\ell}$. A common approach is
to learn an approximate predictive posterior as a function of a latent
distribution over the inducing points $q_{\xi}(\hat{\Theta})$ with
hyperparameters $\xi$, where
$\hat{\Theta}=\\{\hat{\theta}_{i}\\}_{1}^{\ell}\in\hat{D}_{\ell}$ [97, 98, 99,
100], known as a sparse variational 222The term _variational_ in machine
learning parlance denotes a framework to approximate an intractable posterior
distribution over the latent variables, with a more tractable distribution.
This is a different usage of _variational_ in quantum algorithms, which
denotes a class of methods used to approximate the ground state of a quantum
system. GP (SVGP). The SVGP learns $q_{\xi}(\hat{\Theta})$ by maximizing the
evidence lower bound objective (ELBO) function,

$\operatorname*{argmax}_{\gamma,\xi,\hat{\Theta}}\operatorname{ELBO}(\gamma,\xi,\hat{\Theta})=\sum_{i=1}^{m}q_{\xi}\left(\hat{\Theta}\right)\ln
p_{\gamma}\left(\tilde{y}_{\theta_{i}}\mid\hat{\Theta},\theta_{i}\right)-\operatorname{KL}\left(q_{\xi}\left(\hat{\Theta}\right)\|p_{\gamma}\left(\hat{\Theta}\right)\right).$
(18)

Here, $\operatorname{KL}(\cdot)$ denotes the Kullback-Leibler divergence
between the approximate model and the true posterior
$p_{\gamma}(\hat{\Theta})$. The ELBO function does not require the functional
form of the true posterior to optimize the approximate model, as the
intractable posterior is computed by decomposing Bayes’ rule. The summation in
the ELBO is based on the assumption that the training data consists of i.i.d.
samples (which is generally the case for VQE problems).

For this work, we use the stochastic variational inference (SVI) strategy
[98], which provides better scaling by using an explicit representation of the
inducing variables. GP regression can be performed using natural gradients
[101] of the ELBO with respect to the variational parameters to learn the
inducing variables through mini-batched stochastic gradient descent. The SVI
method iteratively proceeds by selecting a mini-batch to compute the natural
gradients, which can be used to individually update the latent variable
distributions
$\bigl{\\{}q_{\xi}\left(\hat{\theta_{i}}\right)\bigr{\\}}_{i=1}^{m}$ using
inference over $q_{\xi}\left(\Theta\right)$ where
$\Theta=\\{\theta_{i}\\}_{i=1}^{m}$. By maximizing the ELBO loss function to
find the hyperparameters $\xi$ and inducing points $\hat{\Theta}$, we can
construct a sparse GP model to approximate the posterior distribution to scale
more favorably in terms of the amount of training data. This SVI approach to
constructing the sparse GP provides the mean function that serves as the
surrogate-based prior in step 3 of algorithm (1).

##  5 Results and Discussion

In this paper, we investigate the performance of our proposed BOPT framework
in solving VQE problems. We consider two molecular systems, $H_{2}$ and
$H_{3}^{+}$, and separate the case studies based on the method of quantum
simulation. Specifically, we consider a state vector simulator (SVS), a fake
backend, and a physical quantum computer. The SVS is a classical quantum
simulation strategy that allows the simulation of the system without hardware
noise by fully determining the state vector $|\Psi\rangle,$ executed in the
circuit. This method would represent the ideal quantum computer, but the cost
of simulating such a circuit grows exponentially with the number of qubits,
limiting its usage to small systems. The fake backend simulations perform the
SVS simulations using a data-driven noise model that resembles specific
quantum computers, allowing for detailed performance analysis on noisy
simulators. These noise models allow the simulation of various transpilation
strategies, which map the circuit to a specific quantum processor topology.
The fake backends used in this work are based on IBM’s quantum computers,
Cairo, Kolkata, and Mumbai. All of them implement the Falcon processor layout
and use an unoptimized transpilation. These fake backends more closely
resemble the behavior of existing quantum hardware, yet making it classically
simulatable. We aimed to corroborate our classical simulations with hardware
runs on the same quantum hardware. Unfortunately, because of the quick
turnaround and development of QPUs, IBM made all classically simulated noise
models unavailable while drafting this manuscript. Finally, we execute our
algorithm on actual quantum hardware using IBM’s Torino quantum computer. This
computer uses the Heron processor and is considerably larger than the
previously mentioned devices (133 vs. 27 qubits), but a noise model is not yet
available for it.

In both SVS and noisy simulators, we use the Hartree-Fock initialization with
a HEA using $R_{y}$ gates and entangling CNOT gates. We consider ansatz depths
of 4 layers, and the chemistry Hamiltonian is constructed on the STO-3G basis.
The $H_{2}$ problem is represented with 4 qubits (16 parameters), and the
$H_{3}^{+}$ problem is represented with 6 qubits (24 parameters). Optimization
runs use a total budget $B=150$, with the BayesOpt methods allocating
$B_{\text{init}}=30$ for initialization. The circuits are sampled using a
high-shot count $\bar{s}=100,000$, except when initializing BOPT methods,
which use a low-shot count of $\underline{s}=1,000$.

In all results, we present our method compared to the Powell optimization
algorithm [102], which has been identified as a well-balanced optimization
algorithm for VQE problems [103]. Our proposed algorithm, uses the GP-TP and
the LCB acquisition function, denoted by _BOPT_ , while using the EI
acquisition function is denoted as _BOPT (EI)_. In addition, we used several
variants of the BayesOpt algorithm to ablate the proposed modifications. The
proposed method using a Matérn 2.5 is denoted _BOPT (Matérn)_ (or _BOPT
(Matérn, EI)_). The standard GP models with a Matérn 2.5 is denoted by the
choice of acquisition function (_EI_ or _LCB_), while methods that use the
periodic kernel are denoted by _EI-p_ or _LCB-p_.

###  5.1 State Vector Simulators

We present the convergence results for the SVS-based VQE in Figure 4. The
upper row shows the results for $H_{2}$, and the lower row for $H_{3}^{+}$.
The leftmost column compares Powell, BOPT and BOPT (Matérn), and the standard
BayesOpt for both EI and LCB acquisition functions. The plot shows that even
standard BayesOpt converges notably faster than Powell, with the fastest
convergence achieved by the BOPT methods. Furthermore, BOPT (Matérn) and
LCB-p/EI-p methods show that both the periodic kernel and the topological
prior have a significant effect on the convergence rate. However, the kernel
effects without the topological prior are less pronounced on the $H_{2}$
problem relative to the $H_{3}^{+}$ problem. A possible explanation for this
behavior is that $H_{2}$ has fewer parameters but the same number of samples,
resulting in a more dense set of samples relative to $H_{3}^{+}$. If there are
sufficient samples, the periodicity of the parameter space can be learned from
observations without the explicit prior, and thus, it would not be expected to
accelerate convergence. In contrast, the less dense set of samples over the
$H_{3}^{+}$ parameter space may be insufficient to identify periodic behavior,
thus benefiting from the explicit inclusion of this prior.

Figure 4: (top) $H_{2}$, (bottom) $H_{3}^{+}$. (left) Ablation of the periodic
kernel and topological prior benchmarked against Powell. (center) Comparison
of BOPT methods with EI and LCB acquisition function. (right) Individual runs
for the BOPT with periodic kernel relative to Powell. All convergence curves
in this figure are generated from 10 runs and show the standard error in the
shaded region.

In the center column of Figure 4, the two BOPT methods (with and without the
periodic kernel) are presented for each acquisition function. The results show
that both acquisition functions perform comparably, with the EI finding
slightly higher energy on $H_{2}$ with slightly faster convergence. These
results, coupled with those from the first column, suggest that the choice of
acquisition function does not impact performance as much as choices affecting
the model. On the other hand, the BOPT framework provides a clear advantage
over the standard GP model, and the periodic kernel provides an advantage
(particularly when used in combination with BOPT) over a standard Matérn 2.5
kernel.

In the third column of Figure 4, the individual runs by Powell and BOPT are
shown. These results show that Powell, which uses a sample-based approximation
of the gradients, requires many samples before finding a direction of
improvement. Conversely, BOPT rapidly converges in the early iterations,
likely due to exploiting the topological prior, but appears to have a tendency
to get “stuck,” preventing further improvement. The exploration-exploitation
trade-off allows BayesOpt to find near-optimal solutions quickly. However, in
some cases, it may over-explore once it finds a near-optimal solution, which
results in behavior similar to that observed in the BOPT method. This
pathology has been observed in other VQE problems and addressed by [65], where
the authors propose using a local optimization strategy initialized from the
best solution found by BayesOpt after expending a given sample budget.

###  5.2 Noisy Simulators

Next, we present the optimization results on using noisy quantum simulators
based on IBM’s quantum computers Cairo, Kolkata, and Mumbai for the $H_{2}$
and $H_{3}^{+}$. First, we present the convergence results in Figure 5 for
Powell, LCB, BOPT (Matérn), and BOPT. The four methods show similar
performance across all three simulators, with the highest average performance
achieved using Mumbai. Additionally, we note that the best-observed energies
are worse for all methods relative to the noise-free case. In contrast, the
relative performance across the 4 methods is similar to the SVS, with Powel
performing the worst, followed by LCB, BOPT (Matérn), and BOPT performing the
best.

Figure 5: (top) $H_{2}$, (bottom) $H_{3}$; convergence results for several
different optimization frameworks on fake backend (left) Kolkata, (center)
Mumbai, and (right) Cairo.

Figure 6: (left) $H_{2}$, (right) $H_{3}$; Comparison of the optimization
framework on the different QPU models (solid lines), and the energies observed
from evaluating the same parameter sequence from optimizing the QPU models on
an SVS (dotted lines). The $H_{2}$ convergence data is generated from 10 runs,
and for $H_{3}$, it is generated from 5 runs. The shaded regions show the
standard error.

On average, the energies observed from the noisy simulators were higher than
those observed from the SVS. However, Figure 6 shows that the sample trace
trajectories (the sequence of parameters corresponding to the best-observed
energy up to the current iteration), which were generated from optimizing the
fake backends, are reevaluated on the SVS. For all three cases, we see that
the actual energies computed from the selected parameters are much lower than
what the fake backends report. Furthermore, we see that the parameters
selected by the BOPT framework on both problems are comparable to those
generated by optimizing the SVS. We also note that the behavior of the trace
trajectories exhibits similar iteration dynamics, such as an increasing bias
for lower energies and similar changes in the standard error. Addressing these
limitations in accurately measuring the observables from a noisy quantum
circuit, yet obtaining the correct convergence toward promising parameters has
been observed in other studies [104] and motivates research in error-
mitigation techniques in NISQ devices toward useful application of these
devices [80].

The errors between the fake backends and the validation results from the
simulations presented in Figure 6 are presented in Table 1 for both problems,
using parity plots and histograms. The histograms show that Kolkata has the
smallest noise, while Mumbai and Cairo have similar noise values. The parity
plots show that a bias is present and appears inversely proportional to the
magnitude of the observed energy. We also note that the zero-bias occurs
around $-0.2$, which is the average initial value for the optimization methods
in Figures 4 and 6. The parity plots also show a small spread in the energies
observed on fake backends at a particular value of the SVS.

Table 1: Errors between measurements on an SVS compared to the fake backends, using parity plots and histograms of errors between the two models, from data generated from the optimization and validation runs in Figure 5.  $H_{2}$ | $H_{3}^{+}$
---|---
|
Kolkata
|
Mumbai
|
Cairo

These results suggest that the probability distribution representing the
hardware noise may have a mean function linear with respect to $\theta$, with
a small variance. These properties of the hardware noise support the use of
BayesOpt methods for efficient state preparation on quantum hardware. The
linear bias does not affect the location of the minima, while the small
variance allows the GP model to model the response surface accurately using
fewer samples. However, the bias in the hardware noise prevents the usage of
current QPUs for determining the ground state energies and requires classical
validation, as shown in Figure 6.

###  5.3 Quantum Hardware

For the last case study, we compare the proposed method with the Powell method
to find the ground state of the $H_{2}$ molecule on physical quantum hardware
and validate it on a classical computer. For the $H_{2}$ problem, we use an
ansatz that restricts the state to real amplitudes, with an initial state of 4
qubits in an equal superposition, $\ket{\Psi_{\text{ref}}}=H^{\otimes
4}\ket{0000}$, where $H^{\otimes 4}$ denotes $4$ Hadamard gates applied on
each of the $4$ qubits of the circuit. The ansatz, shown in Figure 7 (top), is
parameterized by $8$ variables and contains a single entangling layer of $3$
$\operatorname{CNOT}$ gates in a linear layout. The final circuit that runs on
the IBM Torino quantum processor after transpilation is also shown in Figure 7
(bottom). Each optimization run is repeated three times, using a total budget
$B=180$, with $B_{\text{init}}=30$ allocated for initializing the BOPT. The
circuits are sampled using a high-shot count $\bar{s}=8000$ (the maximum
number of shots on IBM Torino being $8192$) and a low-shot count of
$\underline{s}=1000$ for BOPT initialization.

The IBM Quantum job run time for each iteration of the VQEs, comprising $5$
circuits and measurements for the $H_{2}$ Hamiltonian, was $13$ seconds for
both methods. Function calls on the QPU also incur additional communication
time overhead, which includes data transfer, queue, and data retrieval times,
as well as pre-processing and post-processing time. Data transfer time is the
time it takes for job data to be sent from the user’s local machine to the
server and for the job to be created on the server, while data retrieval is
the time taken for the results to be sent back to the user’s local machine
[105]. Note that data pre- and post-processing steps might occur concurrently
with the communication phases. Classical pre-processing and IBM server
communication time could take up to $45$ seconds per iteration for BOPT and up
to $4$ seconds per iteration for Powell, while classical SVS evaluations and
optimization took an average of $25$ seconds per iteration for BOPT with
$B=$180 and $0.1$ seconds per iteration for Powell, using 8 threads of a
4-core i5-9300H CPU - 2.40GHz processor. The three BOPT VQEs took more than
$6$ hours to run within a single IBM Quantum session, while the three Powell
VQEs took almost $5$ hours. Note that the BOPT run time per iteration is small
for the few first iterations, similar to the time taken by Powell, and starts
growing as the number of iterations increases because of inference scaling
cubically with the number of observations.

Figure 7: The original ansatz circuit (top), and the transpiled circuit
(bottom): to run the circuit on IBM Torino, a transpilation step is necessary
to express the ansatz only in terms of gates which are natively supported by
the QPU (up to a global phase of $3\pi/2$), in this case; $\operatorname{RZ}$,
$\operatorname{\sqrt{X}}$, and $\operatorname{CZ}$. It is worth noting that
only $\operatorname{\sqrt{X}}$ and $\operatorname{CZ}$ introduce noise, as
$\operatorname{RZ}$ is implemented as a noiseless virtual gate on IBM QPUs.
Figure 8:  IBM Torino QPU layout. The circuit used for $H_{2}$ ground energy
estimation is mapped to the highlighted qubits within the red box.

The results of the computation on quantum hardware are presented in Figure 9
(left). A dotted line indicates the energies observed when optimizing on the
QPU, and a dash-dotted line shows the energies from evaluating the QPU trace
trajectories on an SVS. The solid lines show the energies observed when
optimizing on the SVS directly. The QPU results show a bias similar to those
observed on the fake backends, which is evident by similar changes in the
trajectories of both option methods. Powell and BOPT achieve similar
performance on the QPU (when validated on the SVS) as they do on the SVS,
suggesting that the variance of the noise does not significantly affect
convergence. Powell’s resilience to hardware noise suggests it would be a
suitable local optimization strategy to deploy from the best solution found by
BOPT. A parity plot between the QPU simulations and SVS is presented in Figure
9 (right). Relative to the fake backends, the QPU noise exhibits higher
variance and appears to have a nonlinear bias.

Figure 9:  (Left) Comparison of BOPT to Powell on quantum hardware simulations
for $H_{2}$. Dotted lines are QPU optimization results on IBM’s Torino, dashed
lines are the SVS validation of the QPU trace trajectories, and solid lines
denote the optimization on the noise-free SVS. (Right) Parity plot of the
evaluated circuit parameters on IBM Torino versus the SVS.

##  6 Conclusions and Future Work

In this work, we present a shot-efficient Bayesian optimization framework to
address state preparation in the variational quantum eigensolver, referred to
as Bayesian optimization with priors on topology (BOPT). The framework uses a
sparse Gaussian process model and low-shot circuit observations to construct a
topological prior, where a residual is learned while performing Bayesian
Optimization. Statistical optimization results generated from state vector
simulators (SVS), SVS with data-driven noise models, and a physical quantum
computer show that this strategy significantly improves the standard Bayesian
optimization algorithm and substantially outperforms the solver Powell,
particularly when using periodic kernels in the Gaussian process models.

Quantum computers are still too noisy for accurate ground state estimation;
however, quantum computers may provide a computationally efficient means for
state preparation using algorithms like the one presented here, which can be
validated by other means, such as classical simulation. Future work to be
considered are noise mitigation strategies. Several strategies are available,
typically provided through the quantum computing interface, to mitigate
hardware noise. These methods typically incur a significant additional
computational cost and may thus be reserved for validating prepared trial
states, eliminating the need for any classical simulation.

Practical near-term quantum computing for chemistry and material science still
requires further algorithmic developments. Future work investigating
modifications that allow efficient Bayesian optimization in high-dimensional
parameter spaces is crucial to solving larger molecules of interest and
simulating systems with more complex ansätze. Furthermore, including geometry
relaxation while solving VQE may increase the dimension of the optimization
variables, depending on how the problem is formulated. Local optimization
methods deployed from the best solutions found with BOPT may make it possible
to find solutions within chemical accuracy. Although these modifications may
make practical problems tractable on quantum computers, extensive benchmarking
against classical methods will be needed to quantify the value of quantum
computation of molecular ground state energies.

##  7 Acknowledgements

F.S. acknowledges funding support from the NSF Graduate Research Fellowship.
M.T.R., N.E.B., and M.M.L. acknowledge the financial support of the European
Union (Grant no. DCIPANAF/2020/420-028) through the African Research
Initiative for Scientific Excellence (ARISE) pilot program. J.A.P.
acknowledges funding support from NSF Grant 2237616. D.B.N. is supported by
the startup grant of the Davidson School of Chemical Engineering at Purdue
University. N.M.T. acknowledges support from the U.S. Department of Energy,
Office of Science, National Quantum Information Science Research Centers, Co-
Design Center for Quantum Advantage under Contract No. DE-SC0012704 (C2QA).

The authors are grateful for support from NASA Ames Research Center. This
research used resources of the National Energy Research Scientific Computing
Center, a DOE Office of Science User Facility supported by the Office of
Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231
using NERSC award ASCR-ERCAP0024469. We acknowledge the use of IBM Quantum
services through the Quantum Collaborative for this work. The views expressed
are those of the authors and do not reflect the official policy or position of
IBM or the IBM Quantum team. The authors thank the National Academies of
Science, Engineering, and Medicine for their support through the Arab-American
Frontiers Fellowship that enabled M.T.R.’s visit to D.B.N.

## References

  * [1] Norm M. Tubman, C. Daniel Freeman, Daniel S. Levine, Diptarka Hait, Martin Head-Gordon, and K. Birgitta Whaley.  Modern Approaches to Exact Diagonalization and Selected Configuration Interaction with the Adaptive Sampling CI Method.  Journal of Chemical Theory and Computation, 16(4):2139–2159, 2020.
  * [2] Attila Szabo and Neil S Ostlund.  Modern quantum chemistry: introduction to advanced electronic structure theory.  Courier Corporation, 2012.
  * [3] Felix A. Faber, Luke Hutchison, Bing Huang, Justin Gilmer, Samuel S. Schoenholz, George E. Dahl, Oriol Vinyals, Steven Kearnes, Patrick F. Riley, and O. Anatole von Lilienfeld.  Prediction Errors of Molecular Machine Learning Models Lower than Hybrid DFT Error.  Journal of Chemical Theory and Computation, 13(11):5255–5264, 2017.
  * [4] Zhengyuan Li, Peng Wang, Xiang Lyu, Vamsi Krishna Reddy Kondapalli, Shuting Xiang, Juan D. Jimenez, Lu Ma, Takeshi Ito, Tianyu Zhang, Jithu Raj, Yanbo Fang, Yaocai Bai, Jianlin Li, Alexey Serov, Vesselin Shanov, Anatoly I. Frenkel, Sanjaya D. Senanayake, Shize Yang, Thomas P. Senftle, and Jingjie Wu.  Directing CO2 electroreduction pathways for selective C2 product formation using single-site doped copper catalysts.  Nature Chemical Engineering, 1(2):159–169, February 2024.
  * [5] Johannes C. B. Dietschreit, Dennis J. Diestler, and Rafael Gómez-Bombarelli.  Entropy and Energy Profiles of Chemical Reactions.  Journal of Chemical Theory and Computation, 19(16):5369–5379, 2023.  PMID: 37535443.
  * [6] Iman Nodozi, Charlie Yan, Mira Khare, Abhishek Halder, and Ali Mesbah.  Neural Schrödinger Bridge With Sinkhorn Losses: Application to Data-Driven Minimum Effort Control of Colloidal Self-Assembly.  IEEE Transactions on Control Systems Technology, 32(3):960–973, 2024.
  * [7] Rolf Findeisen, Martha A. Grover, Christian Wagner, Michael Maiworm, Ruslan Temirov, F. Stefan Tautz, Murti V. Salapaka, Srinivasa Salapaka, Richard D. Braatz, and S.O. Reza Moheimani.  Control on a molecular scale: A perspective.  In 2016 American Control Conference (ACC), pages 3069–3082, 2016.
  * [8] K. Boone, F. Abedin, M.R. Anwar, and K.V. Camarda.  Chapter 8 - Molecular Design in the Pharmaceutical Industries.  In Mariano Martín, Mario R. Eden, and Nishanth G. Chemmangattuvalappil, editors, Tools For Chemical Product Design, volume 39 of Computer Aided Chemical Engineering, pages 221–238. Elsevier, 2016.
  * [9] Akash Arora, David C. Morse, Frank S. Bates, and Kevin D. Dorfman.  Accelerating self-consistent field theory of block polymers in a variable unit cell.  The Journal of Chemical Physics, 146(24):244902, 06 2017.
  * [10] Pengyu Chen and Kevin D. Dorfman.  A soft crystalline packing with no metallic analogue.  Nature Materials, 23(4):455–456, April 2024.
  * [11] Justin C. Bui, Eric W. Lees, Daniela H. Marin, T. Nathan Stovall, Lihaokun Chen, Ahmet Kusoglu, Adam C. Nielander, Thomas F. Jaramillo, Shannon W. Boettcher, Alexis T. Bell, and Adam Z. Weber.  Multi-scale physics of bipolar membranes in electrochemical processes.  Nature Chemical Engineering, 1(1):45–60, January 2024.
  * [12] Rafael Gómez-Bombarelli, Jorge Aguilera-Iparraguirre, Timothy D. Hirzel, Dong-Gwang Ha, Markus Einzinger, Tony Wu, Marc A. Baldo, and Alán Aspuru-Guzik.  Turbocharged molecular discovery of OLED emitters: from high-throughput quantum simulation to highly efficient TADF devices.  In Franky So, Chihaya Adachi, and Jang-Joo Kim, editors, Organic Light Emitting Materials and Devices XX, volume 9941, page 99410A. International Society for Optics and Photonics, SPIE, 2016.
  * [13] Bridgette J. Befort, Ryan S. DeFever, Edward J. Maginn, and Alexander W. Dowling.  Machine Learning-Enabled Optimization of Force Fields for Hydrofluorocarbons.  In Yoshiyuki Yamashita and Manabu Kano, editors, 14th International Symposium on Process Systems Engineering, volume 49 of Computer Aided Chemical Engineering, pages 1249–1254. Elsevier, 2022.
  * [14] Yuri Alexeev, Maximilian Amsler, Marco Antonio Barroca, Sanzio Bassini, Torey Battelle, Daan Camps, David Casanova, Young Jay Choi, Frederic T. Chong, Charles Chung, Christopher Codella, Antonio D. Córcoles, James Cruise, Alberto Di Meglio, Ivan Duran, Thomas Eckl, Sophia Economou, Stephan Eidenbenz, Bruce Elmegreen, Clyde Fare, Ismael Faro, Cristina Sanz Fernández, Rodrigo Neumann Barros Ferreira, Keisuke Fuji, Bryce Fuller, Laura Gagliardi, Giulia Galli, Jennifer R. Glick, Isacco Gobbi, Pranav Gokhale, Salvador de la Puente Gonzalez, Johannes Greiner, Bill Gropp, Michele Grossi, Emanuel Gull, Burns Healy, Matthew R. Hermes, Benchen Huang, Travis S. Humble, Nobuyasu Ito, Artur F. Izmaylov, Ali Javadi-Abhari, Douglas Jennewein, Shantenu Jha, Liang Jiang, Barbara Jones, Wibe Albert de Jong, Petar Jurcevic, William Kirby, Stefan Kister, Masahiro Kitagawa, Joel Klassen, Katherine Klymko, Kwangwon Koh, Masaaki Kondo, Doga Murat Kurkcuoglu, Krzysztof Kurowski, Teodoro Laino, Ryan Landfield, Matt Leininger, Vicente Leyton-Ortega, Ang Li, Meifeng Lin, Junyu Liu, Nicolas Lorente, Andre Luckow, Simon Martiel, Francisco Martin-Fernandez, Margaret Martonosi, Claire Marvinney, Arcesio Castaneda Medina, Dirk Merten, Antonio Mezzacapo, Kristel Michielsen, Abhishek Mitra, Tushar Mittal, Kyungsun Moon, Joel Moore, Sarah Mostame, Mario Motta, Young-Hye Na, Yunseong Nam, Prineha Narang, Yu ya Ohnishi, Daniele Ottaviani, Matthew Otten, Scott Pakin, Vincent R. Pascuzzi, Edwin Pednault, Tomasz Piontek, Jed Pitera, Patrick Rall, Gokul Subramanian Ravi, Niall Robertson, Matteo A.C. Rossi, Piotr Rydlichowski, Hoon Ryu, Georgy Samsonidze, Mitsuhisa Sato, Nishant Saurabh, Vidushi Sharma, Kunal Sharma, Soyoung Shin, George Slessman, Mathias Steiner, Iskandar Sitdikov, In-Saeng Suh, Eric D. Switzer, Wei Tang, Joel Thompson, Synge Todo, Minh C. Tran, Dimitar Trenev, Christian Trott, Huan-Hsin Tseng, Norm M. Tubman, Esin Tureci, David García Valiñas, Sofia Vallecorsa, Christopher Wever, Konrad Wojciechowski, Xiaodi Wu, Shinjae Yoo, Nobuyuki Yoshioka, Victor Wen zhe Yu, Seiji Yunoki, Sergiy Zhuk, and Dmitry Zubarev.  Quantum-centric supercomputing for materials science: A perspective on challenges and future directions.  Future Generation Computer Systems, 160:666–710, 2024.
  * [15] Lindsay Bassman Oftelie, Miroslav Urbanek, Mekena Metcalf, Jonathan Carter, Alexander F Kemper, and Wibe A de Jong.  Simulating quantum materials with digital quantum computers.  Quantum Science and Technology, 6(4):043002, sep 2021.
  * [16] Soonho Kwon, Kelsey A. Stoerzinger, Reshma Rao, Liang Qiao, William A. III Goddard, and Yang Shao-Horn.  Facet-Dependent Oxygen Evolution Reaction Activity of IrO2 from Quantum Mechanics and Experiments.  Journal of the American Chemical Society, 146(17):11719–11725, 2024.  PMID: 38636103.
  * [17] Bryan R. Goldsmith, Jacques Esterhuizen, Jin-Xun Liu, Christopher J. Bartel, and Christopher Sutton.  Machine learning for heterogeneous catalyst design and discovery.  AIChE Journal, 64(7):2311–2323, 2018.
  * [18] Nima Nikbin, Phuong T. Do, Stavros Caratzoulas, Raul F. Lobo, Paul J. Dauenhauer, and Dionisios G. Vlachos.  A DFT study of the acid-catalyzed conversion of 2,5-dimethylfuran and ethylene to p-xylene.  Journal of Catalysis, 297:35–43, 2013.
  * [19] Claire S. Adjiman, Nikolaos V. Sahinidis, Dionisios G. Vlachos, Bhavik Bakshi, Christos T. Maravelias, and Christos Georgakis.  Process Systems Engineering Perspective on the Design of Materials and Molecules.  Industrial & Engineering Chemistry Research, 60(14):5194–5206, 2021.
  * [20] Jan Derezinski.  Asymptotic Completeness of Long-Range N-Body Quantum Systems.  Annals of Mathematics, 138(2):427–476, 1993.
  * [21] Aspuru-Guzik Alan, Anthony D. Dutoi, Peter J. Love, and Martin Head-Gordon.  Simulated Quantum Computation of Molecular Energies.  Science, September 2005.
  * [22] Janus J Eriksen, Tyler A Anderson, J Emiliano Deustua, Khaldoon Ghanem, Diptarka Hait, Mark R Hoffmann, Seunghoon Lee, Daniel S Levine, Ilias Magoulas, Jun Shen, et al.  The ground state electronic energy of benzene.  The journal of physical chemistry letters, 11(20):8922–8929, 2020.
  * [23] Jeongnim Kim, Andrew D Baczewski, Todd D Beaudet, Anouar Benali, M Chandler Bennett, Mark A Berrill, Nick S Blunt, Edgar Josué Landinez Borda, Michele Casula, David M Ceperley, Simone Chiesa, Bryan K Clark, Raymond C Clay, Kris T Delaney, Mark Dewing, Kenneth P Esler, Hongxia Hao, Olle Heinonen, Paul R C Kent, Jaron T Krogel, Ilkka Kylänpää, Ying Wai Li, M Graham Lopez, Ye Luo, Fionn D Malone, Richard M Martin, Amrita Mathuriya, Jeremy McMinis, Cody A Melton, Lubos Mitas, Miguel A Morales, Eric Neuscamman, William D Parker, Sergio D Pineda Flores, Nichols A Romero, Brenda M Rubenstein, Jacqueline A R Shea, Hyeondeok Shin, Luke Shulenburger, Andreas F Tillack, Joshua P Townsend, Norm M Tubman, Brett Van Der Goetz, Jordan E Vincent, D ChangMo Yang, Yubo Yang, Shuai Zhang, and Luning Zhao.  QMCPACK: an open source ab initio quantum Monte Carlo package for the electronic structure of atoms, molecules and solids.  Journal of Physics: Condensed Matter, 30(19):195901, April 2018.
  * [24] Norm M. Tubman, Joonho Lee, Tyler Y. Takeshita, Martin Head-Gordon, and K. Birgitta Whaley.  A deterministic alternative to the full configuration interaction quantum Monte Carlo method.  The Journal of Chemical Physics, 145(4):044112, 2016.
  * [25] Mihail Bogojeski, Leslie Vogt-Maranto, Mark E. Tuckerman, Klaus-Robert Müller, and Kieron Burke.  Quantum chemical accuracy from density functional approximations via machine learning.  Nature Communications, 11(1):5223, October 2020.
  * [26] Yi Shi, Yuming Shi, and Adam Wasserman.  Strong electron correlation from partition density functional theory.  The Journal of Chemical Physics, 159(22):224108, December 2023.
  * [27] Jules Tilly, Hongxiang Chen, Shuxiang Cao, Dario Picozzi, Kanav Setia, Ying Li, Edward Grant, Leonard Wossnig, Ivan Rungger, George H. Booth, and Jonathan Tennyson.  The Variational Quantum Eigensolver: A review of methods and best practices.  Physics Reports, 986:1–128, 2022.  The Variational Quantum Eigensolver: a review of methods and best practices.
  * [28] Benedikt Fauseweh.  Quantum many-body simulations on digital quantum computers: State-of-the-art and future challenges.  Nature Communications, 15(1):2123, March 2024.
  * [29] Erik J Gustafson, Juha Tiihonen, Diana Chamaki, Farshud Sorourifar, J Wayne Mullinax, Andy CY Li, Filip B Maciejewski, Nicolas PD Sawaya, Jaron T Krogel, David E Bernal Neira, and Norman Tubman.  Surrogate optimization of variational quantum circuits, 2024.
  * [30] Harper R. Grimsley, George S. Barron, Edwin Barnes, Sophia E. Economou, and Nicholas J. Mayhall.  Adaptive, problem-tailored variational quantum eigensolver mitigates rough parameter landscapes and barren plateaus.  npj Quantum Information, 9(1), March 2023.
  * [31] Eric R. Anschuetz and Bobak T. Kiani.  Quantum variational algorithms are swamped with traps.  Nature Communications, 13(1), December 2022.
  * [32] 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(7671):242–246, September 2017.
  * [33] Katherine Klymko, Carlos Mejuto-Zaera, Stephen J. Cotton, Filip Wudarski, Miroslav Urbanek, Diptarka Hait, Martin Head-Gordon, K. Birgitta Whaley, Jonathan Moussa, Nathan Wiebe, Wibe A. de Jong, and Norm M. Tubman.  Real-time evolution for ultracompact hamiltonian eigenstates on quantum hardware.  PRX Quantum, 3:020323, May 2022.
  * [34] Andy C. Y. Li, M. Sohaib Alam, Thomas Iadecola, Ammar Jahin, Joshua Job, Doga Murat Kurkcuoglu, Richard Li, Peter P. Orth, A. Barı ş Özgüler, Gabriel N. Perdue, and Norm M. Tubman.  Benchmarking variational quantum eigensolvers for the square-octagon-lattice kitaev model.  Phys. Rev. Res., 5:033071, Aug 2023.
  * [35] Siva Darbha, Milan Kornjača, Fangli Liu, Jan Balewski, Mark R. Hirsbrunner, Pedro Lopes, Sheng-Tao Wang, Roel Van Beeumen, Daan Camps, and Katherine Klymko.  False vacuum decay and nucleation dynamics in neutral atom systems, 2024.
  * [36] Siva Darbha, Milan Kornjača, Fangli Liu, Jan Balewski, Mark R. Hirsbrunner, Pedro Lopes, Sheng-Tao Wang, Roel Van Beeumen, Katherine Klymko, and Daan Camps.  Long-lived oscillations of false and true vacuum states in neutral atom systems, 2024.
  * [37] Javier Robledo-Moreno, Mario Motta, Holger Haas, Ali Javadi-Abhari, Petar Jurcevic, William Kirby, Simon Martiel, Kunal Sharma, Sandeep Sharma, Tomonori Shirakawa, Iskandar Sitdikov, Rong-Yang Sun, Kevin J. Sung, Maika Takita, Minh C. Tran, Seiji Yunoki, and Antonio Mezzacapo.  Chemistry Beyond Exact Solutions on a Quantum-Centric Supercomputer.  arXiv e-prints, page arXiv:2405.05068, May 2024.
  * [38] Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph C. Bardin, Rami Barends, Andreas Bengtsson, Sergio Boixo, Michael Broughton, Bob B. Buckley, David A. Buell, Brian Burkett, Nicholas Bushnell, Yu Chen, Zijun Chen, Yu-An Chen, Ben Chiaro, Roberto Collins, Stephen J. Cotton, William Courtney, Sean Demura, Alan Derk, Andrew Dunsworth, Daniel Eppens, Thomas Eckl, Catherine Erickson, Edward Farhi, Austin Fowler, Brooks Foxen, Craig Gidney, Marissa Giustina, Rob Graff, Jonathan A. Gross, Steve Habegger, Matthew P. Harrigan, Alan Ho, Sabrina Hong, Trent Huang, William Huggins, Lev B. Ioffe, Sergei V. Isakov, Evan Jeffrey, Zhang Jiang, Cody Jones, Dvir Kafri, Kostyantyn Kechedzhi, Julian Kelly, Seon Kim, Paul V. Klimov, Alexander N. Korotkov, Fedor Kostritsa, David Landhuis, Pavel Laptev, Mike Lindmark, Erik Lucero, Michael Marthaler, Orion Martin, John M. Martinis, Anika Marusczyk, Sam McArdle, Jarrod R. McClean, Trevor McCourt, Matt McEwen, Anthony Megrant, Carlos Mejuto-Zaera, Xiao Mi, Masoud Mohseni, Wojciech Mruczkiewicz, Josh Mutus, Ofer Naaman, Matthew Neeley, Charles Neill, Hartmut Neven, Michael Newman, Murphy Yuezhen Niu, Thomas E. O’Brien, Eric Ostby, Bálint Pató, Andre Petukhov, Harald Putterman, Chris Quintana, Jan-Michael Reiner, Pedram Roushan, Nicholas C. Rubin, Daniel Sank, Kevin J. Satzinger, Vadim Smelyanskiy, Doug Strain, Kevin J. Sung, Peter Schmitteckert, Marco Szalay, Norm M. Tubman, Amit Vainsencher, Theodore White, Nicolas Vogt, Z. Jamie Yao, Ping Yeh, Adam Zalcman, and Sebastian Zanker.  Observation of separated dynamics of charge and spin in the fermi-hubbard model, 2020.
  * [39] Pauline J Ollitrault, Matthias Loipersberger, Robert M Parrish, Alexander Erhard, Christine Maier, Christian Sommer, Juris Ulmanis, Thomas Monz, Christian Gogolin, Christofer S Tautermann, et al.  Estimation of electrostatic interaction energies on a trapped-ion quantum computer.  ACS Central Science, 10(4):882–889, 2024.
  * [40] Panagiotis G. Anastasiou, Yanzhu Chen, Nicholas J. Mayhall, Edwin Barnes, and Sophia E. Economou.  TETRIS-ADAPT-VQE: An adaptive algorithm that yields shallower, denser circuit ansätze.  Physical Review Research, 6(1):013254, 2024.
  * [41] Jarrod R. McClean, Jonathan Romero, Ryan Babbush, and Alán Aspuru-Guzik.  The theory of variational hybrid quantum-classical algorithms.  New Journal of Physics, 18(2):023023, February 2016.
  * [42] Lennart Bittel, Jens Watty, and Martin Kliesch.  Fast gradient estimation for variational quantum algorithms.  arXiv preprint arXiv:2210.06484, 2022.
  * [43] K. Mitarai, M. Negoro, M. Kitagawa, and K. Fujii.  Quantum circuit learning.  Phys. Rev. A, 98:032309, September 2018.
  * [44] Aram W. Harrow and John C. Napp.  Low-Depth Gradient Measurements Can Improve Convergence in Variational Hybrid Quantum-Classical Algorithms.  Phys. Rev. Lett., 126:140502, April 2021.
  * [45] Tyson Jones and Julien Gacon.  Efficient calculation of gradients in classical simulations of variational quantum algorithms.  arXiv preprint arXiv:2009.02823, 2020.
  * [46] Erik Berglund, Sarit Khirirat, and Xiaoyu Wang.  Zeroth-Order Randomized Subspace Newton Methods.  In ICASSP 2022 - 2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 6002–6006, 2022.
  * [47] Zhewei Yao, Amir Gholami, Sheng Shen, Mustafa Mustafa, Kurt Keutzer, and Michael W. Mahoney.  ADAHESSIAN: An Adaptive Second Order Optimizer for Machine Learning.  arXiv e-prints, page arXiv:2006.00719, June 2020.
  * [48] Ryan Shaffer, Lucas Kocia, and Mohan Sarovar.  Surrogate-based optimization for variational quantum algorithms.  Phys. Rev. A, 107:032415, March 2023.
  * [49] Mark R. Hirsbrunner, J. Wayne Mullinax, Yizhi Shen, David B. Williams-Young, Katherine Klymko, Roel Van Beeumen, and Norm M. Tubman.  Diagnosing local minima and accelerating convergence of variational quantum eigensolvers with quantum subspace techniques, 2024.
  * [50] Oskar Leimkuhler and K. Birgitta Whaley.  A quantum eigenvalue solver based on tensor networks.  arXiv e-prints, page arXiv:2404.10223, April 2024.
  * [51] Pauline J Ollitrault, Cristian L Cortes, Jerome F Gonthier, Robert M Parrish, Dario Rocca, Gian-Luca Anselmetti, Matthias Degroote, Nikolaj Moll, Raffaele Santagati, and Michael Streif.  Enhancing initial state overlap through orbital optimization for faster molecular electronic ground-state energy estimation.  arXiv preprint arXiv:2404.08565, 2024.
  * [52] Tianyi Hao, Zichang He, Ruslan Shaydulin, Marco Pistoia, and Swamit Tannu.  Variational Quantum Algorithm Landscape Reconstruction by Low-Rank Tensor Completion.  arXiv e-prints, page arXiv:2405.10941, May 2024.
  * [53] Jernej Rudi Finžgar, Martin J. A. Schuetz, J. Kyle Brubaker, Hidetoshi Nishimori, and Helmut G. Katzgraber.  Designing quantum annealing schedules using bayesian optimization.  Physical Review Research, 6(2), April 2024.
  * [54] Lixue Cheng, Yu-Qin Chen, Shi-Xin Zhang, and Shengyu Zhang.  Quantum approximate optimization via learning-based adaptive optimization.  Communications Physics, 7(1), March 2024.
  * [55] J. S. Otterbach, R. Manenti, N. Alidoust, A. Bestwick, M. Block, B. Bloom, S. Caldwell, N. Didier, E. Schuyler Fried, S. Hong, P. Karalekas, C. B. Osborn, A. Papageorge, E. C. Peterson, G. Prawiroatmodjo, N. Rubin, Colm A. Ryan, D. Scarabelli, M. Scheer, E. A. Sete, P. Sivarajah, Robert S. Smith, A. Staley, N. Tezak, W. J. Zeng, A. Hudson, Blake R. Johnson, M. Reagor, M. P. da Silva, and C. Rigetti.  Unsupervised Machine Learning on a Hybrid Quantum Computer.  arXiv preprint arXiv:1712.05771, 2017.
  * [56] Simone Tibaldi, Davide Vodola, Edoardo Tignone, and Elisa Ercolessi.  Bayesian Optimization for QAOA.  IEEE Transactions on Quantum Engineering, 2023.
  * [57] Anthony N. Ciavarella and Ivan A. Chernyshev.  Preparation of the SU(3) lattice Yang-Mills vacuum with variational quantum methods.  Phys. Rev. D, 105:074504, April 2022.
  * [58] Giovanni Iannelli and Karl Jansen.  Noisy bayesian optimization for variational quantum eigensolvers.  arXiv preprint arxiv:2112.00426, 2021.
  * [59] David Eriksson, Michael Pearce, Jacob Gardner, Ryan D Turner, and Matthias Poloczek.  Scalable Global Optimization via Local Bayesian Optimization.  In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
  * [60] C. E. Rasmussen and C. K. I. Williams.  Gaussian Processes for Machine Learning. MIT Press, 2006.
  * [61] Benjamin Letham, Brian Karrer, Guilherme Ottoni, and Eytan Bakshy.  Constrained Bayesian Optimization with Noisy Experiments.  Bayesian Analysis, 14(2):495 – 519, 2019.
  * [62] Donald R Jones, Matthias Schonlau, and William J Welch.  Efficient global optimization of expensive black-box functions.  Journal of Global Optimization, 13:455–492, 1998.
  * [63] Maximilian Balandat, Brian Karrer, Daniel Jiang, Samuel Daulton, Ben Letham, Andrew G Wilson, and Eytan Bakshy.  BoTorch: A Framework for Efficient Monte-Carlo Bayesian Optimization.  Advances in Neural Information Processing Systems, 33:21524–21538, 2020.
  * [64] P. I. Frazier.  A Tutorial on Bayesian Optimization.  arXiv preprint arXiv:1807.02811, 2018.
  * [65] J. Muller, W. Lavrijsen, C. Iancu, and W. de Jong.  Accelerating Noisy VQE Optimization with Gaussian Processes.  In 2022 IEEE International Conference on Quantum Computing and Engineering (QCE), pages 215–225. IEEE Computer Society, September 2022.
  * [66] Shiro Tamiya and Hayata Yamasaki.  Stochastic gradient line Bayesian optimization for efficient noise-robust optimization of parameterized quantum circuits.  npj Quantum Information, 8(1):90, July 2022.
  * [67] Ryan Shaffer, Lucas Kocia, and Mohan Sarovar.  Surrogate-based optimization for variational quantum algorithms.  Phys. Rev. A, 107:032415, March 2023.
  * [68] Lixue Cheng, Yu-Qin Chen, Shi-Xin Zhang, and Shengyu Zhang.  Quantum approximate optimization via learning-based adaptive optimization.  Communications Physics, 7(1):83, March 2024.
  * [69] Chris N. Self, Kiran E. Khosla, Alistair W. R. Smith, Frédéric Sauvage, Peter D. Haynes, Johannes Knolle, Florian Mintert, and M. S. Kim.  Variational quantum algorithm with information sharing.  npj Quantum Information, 7(1):116, July 2021.
  * [70] Milena Röhrs, Alexey Bochkarev, and Arcesio C Medina.  Bayesian optimisation with improved information sharing for the variational quantum eigensolver.  arXiv preprint arXiv:2405.14353, 2024.
  * [71] Kim Nicoli, Christopher J. Anders, Lena Funcke, Tobias Hartung, Karl Jansen, Stefan Kühn, Klaus-Robert Müller, Paolo Stornati, Pan Kessel, and Shinichi Nakajima.  Physics-Informed Bayesian Optimization of Variational Quantum Circuits.  In A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, editors, Advances in Neural Information Processing Systems, volume 36, pages 18341–18376. Curran Associates, Inc., 2023.
  * [72] Ken M. Nakanishi, Keisuke Fujii, and Synge Todo.  Sequential minimal optimization for quantum-classical hybrid algorithms.  Phys. Rev. Res., 2:043158, October 2020.
  * [73] F Sorourifar, D Chamaki, J Paulson, N Tubman, and D.E. Bernal Neira.  Bayesian Optimization Priors for Efficient Variational Quantum Algorithms.  In Proceedings of the 34th European Symposium on Computer Aided Process Engineering / 15th International Symposium on Process Systems Engineering (ESCAPE34/PSE2024), 2024.  To be published.
  * [74] David E. Bernal, Akshay Ajagekar, Stuart M. Harwood, Spencer T. Stober, Dimitar Trenev, and Fengqi You.  Perspectives of quantum computing for chemical engineering.  AIChE Journal, 68(6):e17651, 2022.
  * [75] Davide Castelvecchi.  Quantum computers ready to leap out of the lab in 2017.  Nature, 541(7635):9–10, January 2017.
  * [76] Norbert M. Linke, Dmitri Maslov, Martin Roetteler, Shantanu Debnath, Caroline Figgatt, Kevin A. Landsman, Kenneth Wright, and Christopher Monroe.  Experimental comparison of two quantum computing architectures.  Proceedings of the National Academy of Sciences, 114(13):3305–3310, 2017.
  * [77] Michael A Nielsen and Isaac L Chuang.  Quantum computation and quantum information.  Cambridge university press, 2010.
  * [78] Megan L. Dahlhauser and Travis S. Humble.  Benchmarking characterization methods for noisy quantum circuits.  Phys. Rev. A, 109:042620, April 2024.
  * [79] John Preskill.  Quantum Computing in the NISQ era and beyond.  Quantum, 2:79, August 2018.
  * [80] Youngseok Kim, Andrew Eddins, Sajant Anand, Ken Xuan Wei, Ewout Van Den Berg, Sami Rosenblatt, Hasan Nayfeh, Yantao Wu, Michael Zaletel, Kristan Temme, et al.  Evidence for the utility of quantum computing before fault tolerance.  Nature, 618(7965):500–505, 2023.
  * [81] Peter W. Shor.  Scheme for reducing decoherence in quantum computer memory.  Phys. Rev. A, 52:R2493–R2496, October 1995.
  * [82] Abhinav Kandala, Kristan Temme, Antonio D. Córcoles, Antonio Mezzacapo, Jerry M. Chow, and Jay M. Gambetta.  Error mitigation extends the computational reach of a noisy quantum processor.  Nature, 567(7749):491–495, March 2019.
  * [83] Anbang Wu, Gushu Li, Yuke Wang, Boyuan Feng, Yufei Ding, and Yuan Xie.  Towards Efficient Ansatz Architecture for Variational Quantum Algorithms.  arXiv preprint arXiv:2111.13730, 2021.
  * [84] Xiaoxiao Xiao, Hewang Zhao, Jiajun Ren, Wei-Hai Fang, and Zhendong Li.  Physics-Constrained Hardware-Efficient Ansatz on Quantum Computers that is Universal, Systematically Improvable, and Size-consistent.  Journal of Chemical Theory and Computation, 2024.
  * [85] 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.  Nature communications, 5(1):4213, 2014.
  * [86] Panagiotis G Anastasiou, Nicholas J Mayhall, Edwin Barnes, and Sophia E Economou.  How to really measure operator gradients in ADAPT-VQE.  arXiv preprint arXiv:2306.03227, 2023.
  * [87] 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 Science and Technology, 4(1):014008, 2018.
  * [88] Yudong Cao, Jonathan Romero, Jonathan P. Olson, Matthias Degroote, Peter D. Johnson, Má ria Kieferová, Ian D. Kivlichan, Tim Menke, Borja Peropadre, Nicolas P. D. Sawaya, Sukin Sim, Libor Veis, and Alán Aspuru-Guzik.  Quantum Chemistry in the Age of Quantum Computing.  Chemical Reviews, 119(19):10856–10915, August 2019.
  * [89] Abhinav Anand, Philipp Schleich, Sumner Alperin-Lea, Phillip W. K. Jensen, Sukin Sim, Manuel Díaz-Tinoco, Jakob S. Kottmann, Matthias Degroote, Artur F. Izmaylov, and Alán Aspuru-Guzik.  A quantum computing view on unitary coupled cluster theory.  Chemical Society Reviews, 51(5):1659–1684, 2022.
  * [90] Z. Wang and N. de Freitas.  Theoretical analysis of Bayesian optimisation with unknown Gaussian process hyperparameters.  arXiv preprint arXiv:1406.7758, 2014.
  * [91] Thomas Häner and Damian S. Steiger.  0.5 Petabyte Simulation of a 45-Qubit Quantum Circuit.  arXiv e-prints, page arXiv:1704.01127, April 2017.
  * [92] J Wayne Mullinax and Norm M Tubman.  Large-scale sparse wavefunction circuit simulator for applications with the variational quantum eigensolver.  arXiv preprint arXiv:2301.05726, 2023.
  * [93] Abid Khan, Bryan K. Clark, and Norm M. Tubman.  Pre-optimizing variational quantum eigensolvers with tensor networks.  arXiv e-prints, page arXiv:2310.12965, October 2023.
  * [94] Qiugang Lu, Leonardo D González, Ranjeet Kumar, and Victor M Zavala.  Bayesian Optimization with Reference Models: A Case Study in MPC for HVAC Central Plants.  Computers & Chemical Engineering, 154:107491, 2021.
  * [95] Amazon.  Amazon Web Services, Amazon Braket Pricing, 2024.
  * [96] Microsoft.  Microsoft Build, Azure Quantum Pricing, 2024.
  * [97] Michalis Titsias.  Variational Learning of Inducing Variables in Sparse Gaussian Processes.  In David van Dyk and Max Welling, editors, Proceedings of the Twelth International Conference on Artificial Intelligence and Statistics, volume 5 of Proceedings of Machine Learning Research, pages 567–574. PMLR, 16–18 Apr 2009.
  * [98] James Hensman, Nicoló Fusi, and Neil D. Lawrence.  Gaussian Processes for Big Data.  CoRR, abs/1309.6835, 2013.
  * [99] David R. Burt, Carl Edward Rasmussen, and Mark van der Wilk.  Convergence of Sparse Variational Inference in Gaussian Processes Regression.  Journal of Machine Learning Research, 21(131):1–63, 2020.
  * [100] Henry B Moss, Sebastian W Ober, and Victor Picheny.  Inducing point allocation for sparse Gaussian processes in high-throughput Bayesian optimisation.  In International Conference on Artificial Intelligence and Statistics, pages 5213–5230. PMLR, 2023.
  * [101] Shun-ichi Amari.  Natural Gradient Works Efficiently in Learning.  Neural Computation, 10(2):251–276, 1998.
  * [102] Michael JD Powell.  An efficient method for finding the minimum of a function of several variables without calculating derivatives.  The computer journal, 7(2):155–162, 1964.
  * [103] Harshdeep Singh, Sonjoy Majumder, and Sabyashachi Mishra.  Benchmarking of different optimizers in the variational quantum algorithms for applications in quantum chemistry.  The Journal of Chemical Physics, 159(4):044117, 07 2023.
  * [104] Stefan H Sack and Daniel J Egger.  Large-scale quantum approximate optimization on nonplanar graphs with machine learning noise mitigation.  Physical Review Research, 6(1):013223, 2024.
  * [105] IBM.  Estimate job run time - ibm quantum documentation.