††thanks: The research was partially supported by the Australian Research
Council’s Discovery Project DP200103015 and the Australian Centre of
Excellence for Mathematical and Statistical Frontiers (ACEMS).

# Quantum Speedup of Natural Gradient for Variational Bayes

Anna Lopatnikova<EMAIL_ADDRESS>Minh-Ngoc Tran  minh-
<EMAIL_ADDRESS>Discipline of Business Analytics, The University of
Sydney Business School, The University of Sydney, NSW 2006, Australia

###### Abstract

Variational Bayes (VB) is a critical method in machine learning and
statistics, underpinning the recent success of Bayesian deep learning. The
natural gradient is an essential component of efficient VB estimation, but it
is prohibitively computationally expensive in high dimensions. We propose a
computationally efficient regression-based method for natural gradient
estimation, with convergence guarantees under standard assumptions. The method
enables the use of quantum matrix inversion to further speed up VB. We
demonstrate that the problem setup fulfills the conditions required for
quantum matrix inversion to deliver computational efficiency. The method works
with a broad range of statistical models and does not require special-purpose
or simplified variational distributions.

Keywords. Quantum machine learning, quantum algorithm, Bayesian computation.

##  I Introduction

Our work demonstrates how a general quantum computer could help unlock the
next generation of deep learning applications. We focus on Variational Bayes
(VB), a critical method which has underpinned the recent surge in successful
industrial and research applications of Bayesian deep learning; see [1, 2, 3]
and references therein. Even though VB is efficient and scalable relative to
alternative methods, it remains too computationally intensive for many
practical applications, particularly when the number of variational parameters
$N$ is large (e.g. millions or more).

We propose a quantum-classical algorithm to speed up VB through efficient
computation of natural gradient. The use of natural gradient is seen as one of
the most promising ways to speed up VB [4, 5, 6, 7]. However, this method
requires the computation and inversion of the Fisher information matrix [8, 9,
10], which, for a model with a million parameters, would require a classical
computer to perform, at best, quintillions of operations – $\tilde{O}(N^{d})$
with $2<d\leq 3$ (Section II.1).

To make the problem of natural gradient estimation suitable for quantum
speedup, we reformulate it into a linear regression problem that exploits the
structure of the Fisher matrix (Section II.2). We show that the resulting
problem could be one of the best practical uses of the quantum linear systems
algorithms [11, 12, 13, 14, 15, 16, 17, 18]. These algorithms offer dramatic
speedups but have stringent requirements, rarely met by real-world problems
[13, 19]. We show that the problem of natural gradient estimation for VB meets
the critical requirements, enabling speedup to the classical limit of $O(N)$
needed to write down all the variational parameters (Section III).
Additionally, the linear regression formulation of natural gradient estimation
enables efficient classical computation with time complexity
$\tilde{O}(N^{b})$ with $1\leq b<2$ (Section II.3). We demonstrate the
method’s power in a classical simulation where the VB approximation is a
computationally intensive neural network (Section IV).

##  II A Regression-Based Approach to Natural Gradient Estimation for
Variational Bayes

In this section, we provide a review of VB that uses natural gradient
optimization, then demonstrate the regression-based approach and prove its
convergence.

###  II.1 Variational Bayes with Natural Gradient

One of the most promising classical advances in VB is the use of natural
gradient instead of Euclidean gradient in stochastic gradient descent [4, 5,
6, 7]. Being a geometric object that takes into account the information
geometry of the variational family, natural gradient often leads to faster and
more stable convergence than alternative methods [8]. In situations where
variational parameters lie on Riemannian manifolds, natural gradient is the
only reliable adaptive learning method that speeds up VB training [20].
However, in many cases, accurate computation of natural gradient is
infeasible, because it requires the analytic computation and inversion of the
Fisher information matrix. With computational complexity of
$\tilde{O}(N^{d})$, with $2<d\leq 3$ depending on various algorithms ($d=3$ in
practice), natural gradient estimation is prohibitively computationally
expensive in high-dimensional cases; current practice resorts to heuristic
workarounds that can affect the results of Bayesian inference [9]. We propose
a regression-based approach to natural gradient estimation for VB, which lends
itself to efficient computational methods, including quantum computing
algorithms (Section III).

Let $y$ be the data and $L(\theta)=p(y|\theta)$ the likelihood function based
on a statistical model, with $\theta$ the set of model parameters to be
estimated. Let $p(\theta)$ be the prior. Bayesian inference encodes all
information into the posterior distribution

$p(\theta|y)=\frac{p(\theta)L(\theta)}{p(y)}.$

The main task in Bayesian inference is to approximate the posterior
$p(\theta|y)$, and VB does so by approximating it by a probability
distribution with density $q_{\lambda}(\theta)$, $\lambda\in{\cal M}$ \- the
variational parameter space, belonging to some tractable family of
distributions such as Gaussians or neural networks. Denote by $N$ the size of
$\lambda$. The best $\lambda$ is found by maximizing the lower bound [3]

$\mathcal{L}(\lambda)=\int
q_{\lambda}(\theta)\log\frac{p(\theta)L(\theta)}{q_{\lambda}(\theta)}d\theta={\mathbb{E}}_{q_{\lambda}}\big{(}h_{\lambda}(\theta)\big{)},$
(1)

with

$h_{\lambda}(\theta)=\log\big(p(\theta)L(\theta)\big{missing})-\log
q_{\lambda}(\theta).$ (2)

The Euclidean gradient of the lower bound is [21]

$\nabla_{\lambda}\mathcal{L}(\lambda)={\mathbb{E}}_{q_{\lambda}}\big{(}\nabla_{\lambda}\log
q_{\lambda}(\theta)\times h_{\lambda}(\theta)\big{)}.$ (3)

For maximizing the lower bound objective function $\mathcal{L}(\lambda)$,
Amari [4] shows that the natural gradient, defined as

$\nabla_{\lambda}^{\text{nat}}\mathcal{L}(\lambda):=I^{-1}_{F}(\lambda)\nabla_{\lambda}\mathcal{L}(\lambda),$
(4)

where

$I_{F}(\lambda):=\text{Cov}_{q_{\lambda}}\big{(}\nabla_{\lambda}\log
q_{\lambda}(\theta),\nabla_{\lambda}\log q_{\lambda}(\theta)\big{)}$ (5)

is the Fisher information matrix of $q_{\lambda}$, works much more efficiently
than the Euclidean gradient. The reason is that the natural gradient takes
into account the geometry of the variational family
$\\{q_{\lambda},\lambda\in{\cal M}\\}$ in the optimization. Natural gradient
has been proven vital to the success of VB; see, e.g., Refs. [5, 1] among
others.

The classical VB algorithm with the use of natural gradient proceeds as
follows (Algorithm 1). Let
$\widehat{\nabla_{\lambda}^{\text{nat}}{\mathcal{L}}}(\lambda)$ be an estimate
of the natural gradient
${\nabla_{\lambda}^{\text{nat}}{\mathcal{L}}}(\lambda)$. At step $t$ of
stochastic natural gradient descent, the natural gradient estimate is used to
compute the momentum gradient $Y_{t+1}$, which in turn updates the variational
parameter $\lambda_{t+1}=\lambda_{t}+\alpha_{t}Y_{t+1}$.

Algorithm 1 VB Algorithm with Natural Gradient

  1. 1.

Initialize $\lambda_{0}$ and momentum gradient $Y_{0}$.

  2. 2.

For $t=0,1,...$ and until stopping criterion is met:

    1. a.

Estimate the natural gradient
$\widehat{\nabla_{\lambda}^{\text{nat}}{\mathcal{L}}}(\lambda_{t})$ and apply
gradient clipping.

    2. b.

Update the momentum gradient

$\displaystyle Y_{t+1}$ $\displaystyle=\omega
Y_{t}+(1-\omega)\widehat{\nabla_{\lambda}^{\text{nat}}{\mathcal{L}}}(\lambda_{t}).$
(6)

    3. c.

Update the variational parameter

$\displaystyle\lambda_{t+1}=\lambda_{t}+\alpha_{t}Y_{t+1}.$ (7)

The momentum gradient method, found highly useful in practice, smooths out the
natural gradient estimate as in (6) incorporating the latest estimate
$\widehat{\nabla_{\lambda}^{\text{nat}}{\mathcal{L}}}(\lambda_{t})$ with the
momentum weight $\omega\in(0,1)$. The learning rates $\alpha_{t}$ are required
to satisfy: $1\geq\alpha_{t}\downarrow 0$, $\sum\alpha_{t}=\infty$ and
$\sum\alpha_{t}^{2}<\infty$. We refer the interested reader to Ref. [21] for
more details.

The current practice of stochastic gradient descent for VB often includes
gradient clipping that normalizes the length of
$\widehat{\nabla_{\lambda}^{\text{nat}}{\mathcal{L}}}(\lambda_{t})$ [22].

The main computational bottleneck in Variational Bayes is computing the
gradient estimate
$\widehat{\nabla_{\lambda}^{\text{nat}}{\mathcal{L}}}(\lambda)$. We develop
the regression-based natural gradient estimation method to overcome this
problem, including with the use of a hybrid quantum-classical algorithm.

###  II.2 Estimating Natural Gradient with Regression

The obstacles to practical implementation of natural gradient descent are
well-known. It is difficult to compute the Fisher matrix let alone its
inverse. Even if one could compute the Fisher matrix analytically, solving the
equation in (4) is computationally infeasible in high dimensions. Our
regression-based approach to natural gradient estimation overcomes these
difficulties.

Denote $T(\theta)=\nabla_{\lambda}\log q_{\lambda}(\theta)$, a vector-valued
function often called the _score function_. By noting that
${\mathbb{E}}_{q_{\lambda}}(T(\theta))=0$, the Euclidean gradient of the lower
bound in (3) and the Fisher matrix in (5) can be rewritten as

$\displaystyle\nabla_{\lambda}\mathcal{L}(\lambda)=\text{\rm
Cov}_{q_{\lambda}}\big{(}T(\theta),h_{\lambda}(\theta)\big{)},\;\;\;$
$\displaystyle\text{and}\;\;\;I_{F}(\lambda)=\text{\rm
Cov}_{q_{\lambda}}\big{(}T(\theta),T(\theta)\big{)}.$ (8)

Hence, the natural gradient takes the form

$\displaystyle\nabla_{\lambda}^{\text{nat}}{\mathcal{L}}(\lambda)=I_{F}^{-1}(\lambda)\nabla_{\lambda}\mathcal{L}(\lambda)$
$\displaystyle=\big{[}\text{\rm
Cov}_{q_{\lambda}}(T(\theta),T(\theta))\big{]}^{-1}\text{\rm
Cov}_{q_{\lambda}}\big{(}T(\theta),h_{\lambda}(\theta)\big{)}.$ (9)

Equation (9) invites the interpretation of the natural gradient as the slope
vector $\beta$ of regression of the dependent variable $h_{\lambda}(\theta)$
on the vector of independent variables $T(\theta)$.

More concretely, let $\\{\theta_{i}\\}_{i=1}^{M}$ be $M$ samples from
$q_{\lambda}(\theta)$, and consider the linear regression

$h_{\lambda}(\theta_{i})=\beta_{0}+T(\theta_{i})^{\top}\beta+\epsilon_{i},$
(10)

where the $\epsilon_{i}$ are independently distributed with mean zero. We need
an intercept $\beta_{0}$ in the regression model (10) because
${\mathbb{E}}_{q_{\lambda}}(T(\theta))=0$. We note that
$\beta_{0}={\mathbb{E}}_{q_{\lambda}}(h_{\lambda}(\theta))=\mathcal{L}(\lambda)$;
hence, an estimate of $\beta_{0}$ provides an estimate of the lower bound - an
important quantity in the VB literature. Let $T$ be the design matrix whose
row vectors are $(1,T(\theta_{i})^{\top})$ and $h$ be the response vector of
$h_{\lambda}(\theta_{i})$. The minimum-norm ordinary least squares estimator
of $(\beta_{0},{\beta}^{\top})^{\top}$ is

$\displaystyle g=T^{+}h$ (11)

where $T^{+}$ denotes the pseudo-inverse of $T$ [23]. When $M>N$, $g$ is the
usual OLS estimator, $g=(T^{\top}T)^{-1}T^{\top}h$ and is an unbiased
estimator of $(\beta_{0},{\beta}^{\top})^{\top}$. When $M<N$,
$g=T^{\top}(TT^{\top})^{-1}h$ and is the minimum-norm solution to the
underconstrained linear equation $Tg=h$, and it can be shown that $g$ is a
consistent estimator (in $M$) of $(\beta_{0},{\beta}^{\top})^{\top}$.

In this paper, we focus mainly on the situation where the number of
$\\{\theta_{i}\\}$ samples $M$ is much smaller than the number of parameters
$N$, $M\ll N$. This is the situation observed in the vast majority of
practical applications of VB in deep learning [24, 25, 26].

Denote by $\widehat{\beta}$ the vector $g$ after deleting its very first
component. We use $\widehat{\beta}$ as an estimator of the natural gradient
$\nabla_{\lambda}^{\text{nat}}{\mathcal{L}}(\lambda)$. This regression-based
natural gradient formulation is inspired by the method of Salimans and Knowles
[27] and Malago _et al._ [28] who consider a regression-based method for
estimating the natural gradient when the variational distribution is limited
to the exponential family. Our method is general and can be easily applied to
situations with complicated variational distributions such as neural networks.
Hereafter, for notational simplicity, we will refer to $g$ in (11) as the
natural gradient estimator
$\widehat{\nabla_{\lambda}^{\text{nat}}{\mathcal{L}}}(\lambda)$, although in
practice one must remove the very first component of $g$ to obtain the natural
gradient estimate.

###  II.3 Computational Efficiency of Regression-Based Natural Gradient

The regression-based approach provides an estimate of the natural gradient
that avoids the requirement to compute the Fisher matrix $I_{F}$ analytically.
Pseudo-inversion of the short-and-wide $M\times N$, $M\ll N$ design matrix $T$
can be efficient. The time complexity of the popular Cholesky decomposition
method is $\tilde{O}(NM^{2})$, scaling favorably with $N$, particularly if
$M=O(N^{\alpha})$ with $0\leq\alpha<1/2$, where $\alpha=0$ denotes a
polylogarithmic dependence on $N$. The method delivers significant efficiency
over the naive methods to invert the $N\times N$ Fisher matrix, with time
complexity of $\tilde{O}(N^{d})$, where $2<d\leq 3$ (and usually $d=3$ in
practice). Note that the conjugate gradient method [87, 88], sometimes used as
a classical benchmark for quantum linear systems algorithms [13], does not
apply to the natural gradient problem. The conjugate gradient method requires
that the inverted matrix is symmetric and positive definite. The design matrix
$T$ is not symmetric and its symmetrized version $\begin{pmatrix}0&T\\\
T^{\dagger}&0\end{pmatrix}$ is not positive definite.

Theorem 1 in Section II.4 provides a convergence guarantee for VB optimization
with natural gradient (Algorithm 1) under standard conditions.

###  II.4 Convergence of VB with Regression-Based Natural Gradient

This section provides a convergence analysis of the VB algorithm using the
regression-based natural gradient described in the previous sections. To be
consistent with the conventional notation in the optimization literature, let
us denote $\tilde{\mathcal{L}}(\lambda):=-{\mathcal{L}}(\lambda)$ to be the
objective function, the VB optimization problem becomes

$\min_{\lambda}\tilde{\mathcal{L}}(\lambda).$

Write the estimate
$\widehat{\nabla_{\lambda}^{\text{nat}}\tilde{\mathcal{L}}}(\lambda_{t})$ of
the natural gradient as

$\widehat{\nabla_{\lambda}^{\text{nat}}\tilde{\mathcal{L}}}(\lambda_{t})=\nabla_{\lambda}^{\text{nat}}\tilde{\mathcal{L}}(\lambda_{t})+\epsilon_{t}$
(12)

with $\epsilon_{t}$ the estimate error. In Algorithm 1, by viewing that the
variational parameter $\lambda$ is updated first then the momentum gradient,
it can be equivalently rewritten as:

For $t=0,1,...$

  * •

$\lambda_{t+1}=\lambda_{t}-Y_{t}$.

  * •

$Y_{t+1}=\zeta_{t+1}Y_{t}+\gamma_{t+1}(\nabla_{\lambda}^{\text{nat}}\tilde{\mathcal{L}}(\lambda_{t+1})+\epsilon_{t+1})$.

Here, $\zeta_{t}=\alpha_{t}\omega$ and $\gamma_{t}=\alpha_{t}(1-\omega)$ with
$\omega$ the momentum weight and $\alpha_{t}$ the learning rates. For a
symmetric and positive definite matrix $\Sigma$, we denote by
$\sigma_{\max}(\Sigma)$, $\sigma_{\min}(\Sigma)$ the largest and smallest
eigenvalue of $\Sigma$, and by $\|\cdot\|$ the $\ell^{2}$-norm.

######  Theorem 1.

Assume that the learning rates satisfy the conditions that
$1\geq\alpha_{0}\geq\alpha_{1}\geq\cdots>0$, $\sum\alpha_{t}=\infty$ and
$\sum\alpha_{t}^{2}<\infty$. Furthermore, assume that the following regularity
conditions are satisfied in a neighborhood $\mathcal{N}$ in the
$\lambda$-space:

  * (A1)

The ordinary gradient is bounded, i.e.
$\|\nabla\tilde{\mathcal{L}}(\lambda)\|\leq b_{1}<\infty$.

  * (A2)

$0<b_{2}\leq\sigma_{\min}(I_{F}(\lambda))\leq\sigma_{\max}(I_{F}(\lambda))\leq
b_{3}<\infty$.

  * (A3)

The estimate errors are such that
${\mathbb{E}}(\|\epsilon_{t}\|)<b_{4}<\infty$ for all $t$, and
$\sum_{t}\alpha_{t}{\mathbb{E}}\|\epsilon_{t}\|<\infty$.

Let $\\{\lambda_{t},t=0,1,...\\}$ be the iterates generated from Algorithm 1.
Then, if $\lambda_{0}\in\mathcal{N}$,

$\min_{t=1,...,n}{\mathbb{E}}\|\nabla\tilde{\mathcal{L}}(\lambda_{t})\|\to
0,\;\;\;n\to\infty.$ (13)

Furthermore, if $\tilde{\mathcal{L}}(\lambda)$ is $L$-strongly convex on
$\mathcal{N}$ with the minimum point $\lambda^{*}$, i.e.,

$\tilde{\mathcal{L}}(\lambda)\geq\tilde{\mathcal{L}}(\eta)+\nabla\tilde{\mathcal{L}}(\eta)^{\top}(\lambda-\eta)+\frac{L}{2}\|\lambda-\eta\|^{2},\;\;\forall\lambda,\eta\in\mathcal{N},$

then

$\min_{t=1,...,n}{\mathbb{E}}\|\lambda_{t}-\lambda^{*}\|\to
0,\;\;\;n\to\infty.$ (14)

The practical implication of condition (A3) is that, as the natural gradient
estimate
$\widehat{\nabla_{\lambda}^{\text{nat}}\tilde{\mathcal{L}}}(\lambda_{t})$ is
no longer unbiased when $M<N$, it should be made more and more accurate along
with iterations $t$. That is, one might need to gradually increase the sample
size $M$ in the estimate (11). The result in (13) says that the algorithm does
stop at a stationary point of the target $\tilde{\mathcal{L}}(\lambda)$. The
result in (14) says that, under the strong convexity condition, the algorithm
converges to the minimum point $\lambda^{*}$. The proof can be found in the
Appendix.

Amari [4] demonstrated that (non-stochastic) natural gradient is Fisher-
efficient and, since Amari’s work was published, natural gradient has been
shown to require far fewer iterations than other methods for many models (see,
e.g., Ref. [8] and references therein). However, a rigorous treatment of the
convergence speed of (stochastic or non-stochastic) natural gradient is still
an active open area of research.

##  III Quantum Speedup of Variational Bayes

One of the major benefits of the linear regression formulation of natural
gradient estimation is that it enables quantum speedup of VB. Multiple quantum
algorithms for efficient estimation of linear regression exist [30, 31, 32,
33, 34, 35], under certain conditions delivering material speedup over the
best classical solutions. The quantum linear regression algorithms are based
on quantum linear systems algorithms [11, 12, 13, 14, 15, 16, 17, 18], the
most efficient of which, when provided data input in the form of quantum
states, deliver exponential speedup over the best classical algorithms.
Accounting for quantum state preparation and readout, quantum computers can
estimate quantum natural gradient polynomially faster than the best classical
counterparts, even though the regression formulation supports efficient
classical methods. In this section, we demonstrate that, for certain commonly
met conditions, where the sample size $M$ and the condition number $\kappa$ of
design matrix $T$ (the ratio of the highest and lowest singular values of $T$,
$\kappa\equiv\tau_{max}/\tau_{min}$) scale as a fractional power of the model
size $N$, quantum algorithms (known at the time of writing) can deliver
further polynomial computational advantage.

###  III.1 Quantum Natural Gradient Estimation

The expression for the natural gradient estimate $g=T^{+}h$ in (11) is a
solution to the (underconstrained) linear system $Tg=h$. Quantum linear
systems algorithms, first proposed by Harrow, Hassidim, and Lloyd (HHL) [11],
solve the quantum analogue of the linear systems problem. The algorithms take
as an input the quantum state $\ket{h}$, which encodes the vector $h$ in
amplitude encoding [36, 37], and output the quantum state
$\ket{g}\equiv\ket{T^{+}h}$, which encodes the natural gradient estimate $g$
in a quantum state in amplitude encoding.

The runtime and resource requirements of quantum linear systems algorithms
depend on the dimensions $M$ and $N$ of the matrix $T$, its condition number
$\kappa$ and the required precision $\epsilon$. The runtime of HHL quantum
linear systems algorithm with updated Hamiltonian simulation [38] scales as
$O(M\kappa^{2}/\epsilon\;\text{poly}\,\log(N\kappa/\epsilon))$ [16]. When
$M=O(\text{poly}\;\log(N))$ and $\kappa=O(\text{poly}\;\log(N))$, the HHL
algorithm outputs $\ket{g}$ exponentially faster than the best classical
algorithm outputs $g$. An alternative quantum linear systems algorithm by
Wossnig, Zhao, and Prakash (WZP) [15] has time complexity of
$O(\kappa^{2}\sqrt{N}\text{poly}\;\log N/\epsilon)$ for a bounded spectral
norm of $T$; this method is attractive relative to HHL when $M>O(\sqrt{N})$.
The algorithms assume sparse access to a black-box Hamiltonian oracle and an
address oracle [39] to make the entries of the design matrix $T$ available to
the algorithm. Note that the need to estimate the non-zero entries of $T$
limits the efficiency of any algorithm – classical or quantum. The computation
of $T$ requires $O(NM)$ operations and, unless an efficient way of estimating
$T$ is found (e.g. for specific variational distributions
$q_{\lambda}(\theta)$), $O(NM)$ bounds from below the computational complexity
of natural gradient estimation by classical, quantum, or the so-called
quantum-inspired methods.

When $\kappa=O(\text{poly}\,\log(N))$ and $M=O(\text{poly}\,\log(N))$ and the
required precision does not have to be high, i.e.
$\epsilon^{-1}=O(\text{poly}\,\log(N))$, quantum-inspired algorithms [40] can
be most efficient. Quantum-inspired algorithms are randomized classical
algorithms designed to emulate some computational advantages of quantum
algorithms using sampling that mimics quantum measurements [41]. The method’s
time complexity does not (directly) depend on the size of the problem $N$, but
scales as $\tilde{O}(\kappa^{16}k^{6}\|T\|_{F}^{6}/\epsilon^{6})$, where
$\tilde{O}(f(x))\equiv\tilde{O}(f(x)\text{poly}\;\log(g(x)))$, $k$ is the rank
of $T$ (we can assume $k=M)$, and $\|\cdot\|_{F}$ denotes the Frobenius norm.

###  III.2 Hybrid Quantum-Classical VB Algorithm

The quantum linear system algorithm provides an efficient way to estimate the
natural gradient $g$, but its output is a quantum state $\ket{g}$. To use
$\ket{g}$ in VB, we face a design choice: either use $\ket{g}$ in a quantum
algorithm – i.e. run the entire VB gradient descent algorithm on a quantum
computer – or read out the information encoded in $\ket{g}$ for use on a
classical computer. Even though readout can be a resource-intensive step, we
propose to use a hybrid quantum-classical approach, illustrated in Figure 1,
where the quantum computer calculates $\ket{g}$, and the classical computer
uses the information to update parameters $\lambda$, to sample from the
variational distribution $q_{\lambda}(\theta)$, and to estimate the likelihood
vector $h$ and the design matrix $T$.

Figure 1: A schematic representation of a single iteration of the hybrid
quantum-classical VB algorithm. The update of the variational parameter vector
$\lambda_{t}$ starts in the upper left corner and proceeds clockwise. The
dashed line highlights the steps performed by a quantum computer, where
$\mathcal{A}$ stands for the quantum algorithm. The output of the quantum
algorithm is an consistent estimate of the natural gradient $\hat{g}_{t}$.
This estimate feeds into the momentum gradient $Y_{t+1}$, defined in Eq. (6),
used to update the variational parameter vector and obtain its next iteration
$\lambda_{t+1}$. If the stopping criterion is reached, the algorithm outputs
$\lambda_{t+1}$ as the desired estimate of the optimal variational parameter
vector; if the stopping criterion is not satisfied, the algorithm proceeds
with the update of $\lambda_{t+1}$.

The hybrid quantum-classical approach is more suitable for VB gradient descent
because iterative quantum algorithms are, in general, resource intensive.
Quantum transformations are linear; to create effectively non-linear
transformations repeated quantum registers are used, some of which are
destroyed using projective quantum measurements [see, e.g. 18, 42]. If the
projective measurement yields the desired result with probability $p<1$
(usually $p\leq 1/2)$, the measurement destroys, on average, $1/p-1$ quantum
registers. For an iterative quantum algorithm to yield the desired result
after $T$ iterations, $O(p^{-T})$ qubit registers are required. Therefore,
end-to-end quantum algorithms work for iterative computations only if a few
iterative steps are sufficient. Bausch [89] proposes a recurrent method with
post-selection in the context of quantum Recurrent Neural Networks and argues
that, for approximate post-selection it is possible to reduce this overhead.
Another possible approach is to run the iterative algorithm without post-
selection for $T$ steps and then amplify the desired end state using Quantum
Amplitude Amplification (QAA) [50]. Even with these improvements, the
exponential dependence on the number of iterations $T/z$ (where $z=2$ for QAA
method) remains. For example, Rebentrost et al. [42] propose an end-to-end
quantum algorithm to find the minimum of a homogeneous polynomial using
Newton’s method (which leverages the density matrix exponentiation method of
Lloyd _et al._ [90]), suitable for situations where only a few iterative steps
are required. If an optimization requires more than a few iterative steps, a
hybrid quantum-classical algorithm, such as the one proposed in this paper, is
a better choice, provided it can handle the classical-quantum-classical
handoff efficiently.

The hybrid quantum-classical VB algorithm with quantum natural descent is
summarized in Algorithm 2. It extends the classical Algorithm 1 to embed a
quantum subroutine to estimate the natural gradient. We consider two versions
of the quantum natural gradient estimator: The first, Algorithm 3, includes
full readout of the output state $\ket{g}$. The second, Algorithm 4, includes
partial readout of $\ket{g}$ focused on the most significant coordinates of
$g$. We describe the subroutines in the next two subsections.

###  III.3 Quantum Natural Gradient with Full Readout

The quantum natural gradient subroutine with full readout (Algorithm 3) starts
with the initialization of the quantum state $\ket{h}$ and making available
the entries of the design matrix $T$. Note that because $h$ is a vector of
length $M$ and $M\ll N$, preparation of the input state $\ket{h}$ requires
only $O(M/\log M)$ operations, negligible relative to the complexity of other
parts of the algorithm. As discussed above, the computation of the $M\times N$
entries of matrix $T$ bounds from below the time complexity of both classical
and quantum algorithms for natural gradient estimation. Once matrix $T$ is
stored on a classical computer, time-efficient methods have been proposed to
make it available for quantum computation either in a quantum random access
memory (QRAM) [45, 46, 47] or a classical read-only memory for quantum access
(QROM) [14]. With these inputs, the quantum computer executes a quantum linear
systems algorithm (e.g. HHL or WPZ) to perform singular value inversion and
the computation of the quantum state $\ket{g}$. The last step of the quantum
subroutine is the readout of state $\ket{g}$ using a series of quantum
measurements.

Algorithm 2 VB with Quantum Natural Gradient

  1. 1.

On a classical computer, initialize $\lambda_{0}$ and momentum gradient
$Y_{0}$.

  2. 2.

For $t=0,1,...$and until the stopping criterion is met:

    1. a.

Generate $\theta_{j}\sim q_{\lambda^{(t)}}(\theta)$, $j=1,...,M$, and form the
matrix $T$ and vector $h$.

    2. b.

Use a classical, quantum (Algorithm 3 or 4), or quantum-inspired algorithm to
obtain the natural gradient estimate $\widehat{g}_{t}$.

    3. c.

Update the momentum gradient

$Y_{t+1}=\omega Y_{t}+(1-\omega)\widehat{g_{t}}.$

    4. d.

Update the variational parameter

$\displaystyle\lambda_{t+1}=\lambda_{t}+\alpha_{t}Y_{t+1}.$ (15)

Since quantum measurement of state $\ket{g}$ in the basis
$\\{\ket{i}\\}_{i=1}^{N}$ yields the state $\ket{i}$ with probability
$|g_{i}|^{2}/\|g\|^{2}$, repeated measurements allow us to estimate
$|g_{i}|/\|g\|$ – the normalized absolute values of each gradient coordinate.
To extract the sign of $g_{i}$, an additional step is required. At the
beginning of the quantum computation, we pull in an additional auxiliary qubit
and apply a Hadamard gate [see, e.g. 48, 37] to create two computational
branches. In the computational branch of the auxiliary qubit in state
$\ket{0}$ we run the HHL or WPZ algorithm to create the state $\ket{g}$; in
the $\ket{1}$-branch we create the uniform state
$\ket{\gamma}\equiv\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\ket{i}$. The resulting
state is

$\displaystyle\frac{1}{\sqrt{2}}\ket{0}\ket{g}+\frac{1}{\sqrt{2}}\ket{1}\ket{\gamma}.$
(16)

which is equivalent to

$\displaystyle\frac{1}{2}\sum_{i=1}^{N}(g_{i}+\frac{1}{\sqrt{N}})\ket{+}\ket{i}+\frac{1}{2}\sum_{i=1}^{N}(g_{i}-\frac{1}{\sqrt{N}})\ket{-}\ket{i},$

where $\ket{\pm}\equiv\frac{1}{\sqrt{2}}(\ket{0}\pm\ket{1})$ represent
Hadamard basis states [48, 37] of the auxiliary qubit. To streamline notation
and without loss of generality, we assumed $\|g\|=1$ (we discuss the
implications of measuring $g/\|g\|$ rather than $g$ at the end of this
subsection).

Algorithm 3 Quantum Natural Gradient with Full Readout

  1. 1.

Initialize registers and inputs:

    1. a.

On a classical computer, initialize two counter vectors – length $N$ integer
vectors $n^{+}$ and $n^{-}$ to $\vec{0}$.

    2. b.

On a quantum computer, initialize a data register of $\text{ceil}(\log N)$
qubits and an auxiliary qubit.

    3. c.

On a quantum or classical computer, initialize appropriate data structure for
efficient quantum access to elements of $T$

  2. 2.

Obtain natural gradient estimate $\widehat{g}$:

    1. a.

For $j=0,1,\ldots$ and until stopping criterion is met

       * •

Put the auxiliary qubit in the state $\frac{\ket{0}+\ket{1}}{\sqrt{2}}$.

       * •

In the data register, in the $\ket{1}$-branch of the auxiliary qubit, create a
uniform state $\ket{\gamma}$.

       * •

In the data register, in the $\ket{0}$-branch of the auxiliary qubit:

         * –

Initialize the state $\ket{h}$.

         * –

Perform quantum matrix inversion using HHL or WZP algorithm. (Additional
auxiliary registers and qubits are appended as needed in this step.) Obtain
quantum state $\ket{g}$.

       * •

Perform simultaneous quantum measurement of the data register in the
computational basis and the auxiliary qubit in the Hadamard basis. Receive
paired outcome $[i;p]$, were $i=0,\ldots,N-1$ and $p\in\\{+,-\\}$.

       * •

Update counter vector $n^{p}_{i}=n^{p}_{i}+1$.

    2. b.

Construct the unbiased estimate $\widehat{g}$ as in (18).

Simultaneous measurements of the register holding the superposition of
$\ket{g}$ and $\ket{\gamma}$, which we call the data register, and the
auxiliary register yield results that comprise an unbiased estimator of
${g}_{i}$. We measure the data register in the computational basis and the
auxiliary qubit in the Hadamard basis. The probability of measuring the
$\ket{\pm}$ state in the auxiliary qubit and $\ket{i}$ in the data register is
$|g_{i}\pm\frac{1}{\sqrt{N}}|^{2}/4$. We perform repeated measurements and
count the outcomes. Let $n^{\pm}_{i}$ be the number of measurements of the
data register that yield the state $\ket{i}$ when the measurement of the
auxiliary qubit yields $\ket{\pm}$, and
$n_{T}=\sum_{i=1}^{N}(n^{+}_{i}+n^{-}_{i})$ be the total number of
measurements taken. Let $\mathbb{E}_{m}$ indicate expectation over measurement
results. Then,

$\displaystyle\mathbb{E}_{m}(n^{+}_{i}-n^{-}_{i})=$
$\displaystyle=n_{T}\Big{(}\frac{|g_{i}+1/\sqrt{N}|^{2}}{4}-\frac{|g_{i}-1/\sqrt{N}|^{2}}{4}\Big{)}=n_{T}\frac{g_{i}}{\sqrt{N}},$
(17)

which implies that measurement results yield an unbiased estimator of $g_{i}$:

$\displaystyle\widehat{g}_{i}=\sqrt{N}\frac{n^{+}_{i}-n^{-}_{i}}{n_{T}},\;\;\;\mathbb{E}_{m}(\widehat{g}_{i})$
$\displaystyle=g_{i}.$ (18)

We denote by $\widehat{g}=(\widehat{g}_{1},...,\widehat{g}_{N})^{\top}$ the
quantum natural gradient which is a classical vector.

Technically, this algorithm results in a unbiased estimator of the normalized
natural gradient, i.e. $g/\|g\|$, because of the unit norm requirement of
quantum states in quantum computation. This makes quantum natural gradient
suit naturally to gradient clipping (Section II.1) - a method used in
stochastic gradient descent that clips off the length of the gradient vector
to help stabilize the optimization [22].

###  III.4 Quantum Natural Gradient with Gauss-Southwell Rule

Full readout of $g/\|g\|$ requires $O(N/\epsilon^{2})$ quantum measurements to
estimate the vector $g$ within precision $\epsilon$. Full readout, however,
may not be required when estimation of top gradient coordinates (by absolute
value) is sufficient for VB convergence. We propose a new method we call the
quantum Gauss-Southwell rule developed specifically for quantum natural
gradient estimation (Algorithm 4). The method extracts one of the highest-
value directions of the natural gradient estimate vector $g$ within precision
$\epsilon$ using $O(\sqrt{N}/\epsilon)$ measurements.

The quantum Gauss-Southwell rule is a quantum analogue of the Gauss-Southwell
rule for efficient classical coordinate descent [49]. At each step of
coordinate descent, the classical Gauss-Southwell rule selects the top
gradient coordinate (by absolute value). Classical Gauss-Southwell rule can
deliver efficient linear convergence even for high-dimensional models; its
primary shortcoming is its computational intensity relative to randomized
coordinate descent, because it requires the computation of the full gradient
vector. The quantum approach estimates the full gradient vector highly
efficiently in quantum parallel, making quantum Gauss-Southwell rule efficient
and able to handle _natural gradient_ Gauss-Southwell descent, infeasible on a
classical computer.

Algorithm 4 Quantum Natural Gradient with Gauss-Southwell Rule

  1. 1.

Initialize registers and inputs:

    1. a.

On a quantum computer, initialize a data register of $\text{ceil}(\log N)$
qubits and an auxiliary qubit.

    2. b.

On a quantum or classical computer, initialize appropriate data structure for
efficient quantum access to elements of $T$

    3. c.

In the data register, initialize the state $\ket{h}$.

  2. 2.

Determine high-value gradient coordinate $i$:

    1. a.

In the data register, perform quantum matrix inversion using HHL or WZP
algorithm. (Additional auxiliary registers and qubits are appended as needed
in this step.)

    2. b.

Perform quantum measurement of the data register in the computational basis.
Receive paired outcome $i$, were $i=0,\ldots,N-1$.

  3. 3.

Obtain natural gradient coordinate estimate $\widehat{g_{i}}$:

    1. a.

Uncompute the data register

    2. b.

Put the auxiliary qubit in the state $\frac{\ket{0}+\ket{1}}{\sqrt{2}}$.

    3. c.

In the data register, in the $\ket{1}$-branch of the auxiliary qubit, create a
uniform state $\ket{\gamma}$.

    4. d.

In the data register, in the $\ket{0}$-branch of the auxiliary qubit:

       * •

Initialize the state $\ket{h}$.

       * •

Perform quantum matrix inversion using HHL or WZP algorithm. (Additional
auxiliary registers and qubits are appended as needed in this step.) Obtain
quantum state $\ket{g}$.

    5. e.

Perform Quantum Amplitude Estimation with $\ket{i}\ket{+}$ as the reference
state. Receive estimate of
$\frac{1}{2}\big{|}\frac{g_{i}}{\|g\|}+\frac{1}{\sqrt{N}}\big{|}$.

    6. f.

Repeat Steps 3a through 3d

    7. g.

Perform Quantum Amplitude Estimation with $\ket{i}\ket{-}$ as the reference
state. Receive estimate of
$\frac{1}{2}\big{|}\frac{g_{i}}{\|g\|}-\frac{1}{\sqrt{N}}\big{|}$.

    8. h.

Construct the unbiased estimate $\widehat{g_{i}}$ as in (19).

Let $\ket{g}$ be the output state of the quantum natural gradient estimation
algorithm. We measure $\ket{g}$ in the computational basis. The measurement
yields a computational basis state $\ket{i}$ with probability $|g_{i}|^{2}$.
The quantum measurement naturally yields the index of one of the highest
absolute-value gradient directions with high probability, analogously to the
classical Gauss-Southwell coordinate descent which selects the highest
absolute-value gradient direction for the parameter update. Having determined
the index $i$ of a significant gradient coordinate, we use Quantum Amplitude
Estimation [50] with target state $\ket{i}$ to extract the value of $g_{i}$ in
$O(\sqrt{N}/\epsilon)$ operations. Applying Quantum Amplitude Estimation
directly to the state $\ket{g}$ yields the absolute value $|g_{i}|$ rather
than the signed value $g_{i}$. As in the case of full readout (Section III.1),
in order to extract the sign of $g_{i}$, we append an auxiliary qubit in order
to be able to estimate
$\frac{1}{2}\big{|}\frac{g_{i}}{\|g\|}\pm\frac{1}{\sqrt{N}}\big{|}$ as
described in Algorithm 4:

$\displaystyle\frac{1}{4}\Big{|}\frac{g_{i}}{\|g\|}+\frac{1}{\sqrt{N}}\Big{|}^{2}-\frac{1}{4}\Big{|}\frac{g_{i}}{\|g\|}-\frac{1}{\sqrt{N}}\Big{|}^{2}=\frac{1}{\sqrt{N}}\frac{g_{i}}{\|g\|}$
(19)

Formal proof of convergence of quantum VB with Gauss-Southwell rule is outside
the scope of this work and is a topic of future research.

###  III.5 Complexity of Quantum Natural Gradient Estimation Algorithms

The hybrid quantum-classical VB algorithm provides a flexible framework that
can leverage one of a number of quantum, classical, or even classical quantum-
inspired [51] algorithms for natural gradient estimation. The key factors in
determining which natural gradient estimation algorithm is to be used are the
number of samples $M$ and the condition number $\kappa$ of the design matrix
$T$. The desired precision $\epsilon$ is also a driver of computational
complexity, but it is not the limiting variable in the the stochastic
approximation to the natural gradient. We compare two quantum algorithms – HHL
and WPZ, with two measurement schemes (full readout and Gauss-Southwell rule)
to read off $\hat{g}$ – with a classical singular value decomposition (SVD)
algorithm and a quantum-inspired algorithm [40].

Table 1 summarizes the time complexity of algorithms for natural gradient
estimation. The time complexity of quantum algorithms takes into account
quantum state preparation and readout. Each of the algorithms we consider is
the most efficient algorithm in some regime defined by the scaling of $M$ and
$\kappa$ with the size of the model $N$. Quantum algorithms are most efficient
when $M=\tilde{O}(N^{\alpha})$, where $0<\alpha<1$, and $\kappa<C_{\kappa}M$,
where $C_{\kappa}$ is a constant. When $\alpha=0$ and
$\kappa=\tilde{O}(\text{poly}\,\log(N))$, i.e. the matrix $T$ is sparse and
well-conditioned, quantum-inspired algorithms are most efficient. Table 2
summarizes the most efficient algorithms given $M$ and $\kappa$. The summary
assumes that the scaling of all algorithms with problem size is bounded from
below by the need to estimate the design matrix $T$; should more efficient
ways to estimate non-zero entries of $T$ become available, quantum algorithms
will become the most efficient choice across most combinations of $M$ and
$\kappa$.

Table 1: Time complexity of pseudoinverse algorithms as a function of parameter vector size $N$, number of samples $M$, condition number $\kappa$ of the design matrix $T$ (the ratio of the highest and lowest singular values of $T$, $\kappa\equiv\tau_{max}/\tau_{min}$), the required precision $\epsilon$, and the Frobenius norm $\|T\|_{F}$. Quantum algorithms by Harrow, Hassidim, and Lloyd [11, 38] (_HHL_) and Wossnig, Zhao, and Prakash [15] (_WZP_) with full readout (_FR_) or Gauss-Southwell rule (_GR_) are compared with classical SVD algorithms (_CL_) and quantum-inspired classical algorithm by Gilyen, Lloyd, and Tang [40] (_QI_).  Algorithm | Complexity
---|---
_Quantum:_
| HHL-FR | $\max(O(\frac{NM\kappa^{2}}{\epsilon^{3}}\text{poly}\log(\frac{N\kappa}{\epsilon})),O(NM))$
| HHL-GS | $\max(O(\frac{\sqrt{N}M\kappa^{2}}{\epsilon^{2}}\text{poly}\log(\frac{N\kappa}{\epsilon})),O(NM))$
| WZP-FR | $\max(O(\frac{N^{3/2}\kappa^{2}}{\epsilon^{3}}\text{poly}\log(N)),O(NM))$
| WZP-GS | $\max(O(\frac{N\kappa^{2}}{\epsilon^{2}}\text{poly}\log(N)),O(NM))$
_Classical:_
| CL | $\tilde{O}(NM^{2})$
_Quantum-Inspired:_
| QI | $\tilde{O}(M^{6}\kappa^{16}\|T\|_{F}^{6}/\epsilon^{6})$
Table 2: Relative efficiency of pseudoinverse algorithms across design matrix sparsity/rank and conditioning regimes. Table lists most efficient algorithm(s) for each regime. The number of samples $M$ determines the sparsity or rank of the design matrix $T$. The conditioning number $\kappa$ is the ratio of the maximum and minimum singular values of $T$, $\kappa=\tau_{max}/\tau_{min}$. Scaling of $M$ and $\kappa$ with model size $N$, $M=\tilde{O}(N^{\alpha})$ and $\kappa=\tilde{O}(N^{\delta})$, defines regimes that favor one or a few of selected algorithms. The scaling exponents $\alpha$ and $\delta$ equal to $0$ denote polylogarithmic dependence on $N$. The algorithms include: Classical SVD algorithms (_CL_); Quantum-inspired – Gilyen, Lloyd, and Tang [40] (_QI_); Quantum – Harrow, Hassidim, and Lloyd [11, 38] (_HHL_) and Wossnig, Zhao, and Prakash [15] (_WZP_) with full readout (_FR_) or Gauss-Southwell rule (_GR_). (Quantum algorithm with no readout scheme specified indicates that the algorithm is favored independently of readout method.) |  | $M=\tilde{O}(N^{\alpha})$ |
---|---|---|---
|  |  $\alpha=0$abcdefg |  $0<\alpha\leq 1/4$abcd |  $1/4<\alpha\leq 1/2$ab |  $1/2<\alpha<1$ab |
$\kappa=\tilde{O}(N^{\delta})$a | $\delta=0$ | QI | HHL, WZP-GS a | HHL, WZP a | HHL, WZP |
$0<\delta\leq 1/4$ | CL, HHL-GS a | HHL-GS | HHL-GS, WZP-GS | HHL-GS, WZP |
$1/4<\delta\leq 1/2$a | CL | CL, HHL-GS | CL, HHL-GS | WZP-GS |
|  |  | WZP-GS a |  |
| $1/2<\delta<1$ | CL | CL | CL | CL |
|  |  |  |  | WZP-GS |

##  IV A Numerical Example: Stiefel Neural Network VB

This section provides a numerical example that highlights the need for quantum
speed up to unlock further applications in deep learning. We demonstrate,
using classical simulation, how quantum natural gradient can enable inference
on a high-dimensional dataset using an expressive, but computationally
demanding deep neural network VB framework.

The recent VB literature has called for the use of flexible and expressive
variational distributions $q_{\lambda}$ that are able to approximate
sufficiently well a wide range of posterior distributions. We consider such a
flexible variational distribution by the following construction

$Z_{0}\sim N(0,I),\;\;Z_{k}=\zeta_{k}(W_{k}Z_{k-1}+b_{k}),\;\;k=0,...,K,$ (20)

where $N(0,I)$ denotes the multivariate standard normal distribution, the
$\zeta_{k}$ are activation functions such as tanh or sigmoid, $W_{k}$ and
$b_{k}$ are coefficient matrices and vectors respectively. The distribution of
$Z_{K}$ gives us the variational distribution $q_{\lambda}$ with
$\lambda=(W_{1},b_{1},...,W_{K},b_{K})$. The construction in (20) is a neural
network with $K$ layers and can be considered as an example of normalizing
flows - a class of methods for constructing expressive probability
distributions via a composition of simple bijective transformations [52].
There are two desirable properties of a normalizing flow: the transformation
from $Z_{k-1}$ to $Z_{k}$ must be invertible and, for computational
efficiency, it must be easy to compute the Jacobian determinant. To this end,
we impose the following constraint on the $W_{k}$:

$W_{k}^{\top}W_{k}=I,k=1,...,K,$

that is, $W_{k}$ belongs to the Stiefel manifold. This constraint makes it
easy to compute the inverse
$Z_{k-1}=W_{k}^{\top}(\zeta_{k}^{-1}(Z_{k})-b_{k})$, and the Jacobian is a
diagonal matrix. We refer to the construction of probability distributions in
(20) as the Stiefel normalizing flow or the Stiefel neural network. As
$\lambda$ has a rich geometric structure, adaptive learning methods such as
Adam and AdaGrad fail to work, leaving the natural gradient method the only
option for training $\lambda$.

Being a deep neural network, the Stiefel normalizing flow is expressive and
can approximate a wide range of probability distributions. However, it is
extremely computationally expensive to use this normalizing flow in VB as the
number of parameters in $\lambda$ is quadratic in the number of model
parameters, which poses a real challenge for the current computational
technologies.

This section simulates and demonstrates the performance of quantum natural
gradient for the Stiefel neural network VB in a small dimensional setting. We
make it clear upfront that this numerical example is not run on a real quantum
computer; despite many hardware advances made recently, the availability of
such a computer for general use is still years away. This example uses a
quantum-like surrogate to test the performance of the quantum measurement step
for reading off the quantum state $\ket{g}$ into a classical vector
$\widehat{g}$ as in (18). Given the manifold constraints of the $W_{k}$, we
use the VB method on manifolds of Tran _et al._ [20].

We run the simulated algorithm on a standalone basis instead of contrasting
its efficiency against that of alternative methodologies, such as stochastic
gradient descent or adaptive learning methods (e.g. Adam or AdaGrad), because
these methodologies fail to work for the Stiefel neural network.

We use the gene expression dataset, Colon, from
http://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary.html. This
dataset has 62 observations on a binary response variable and 2000 covariates
(we only use the first 100 covariates, plus an intercept, in this example). We
consider fitting a logistic regression model to this dataset, and approximate
its posterior by a Stiefel neural network with $K=2$, hence, the size of
$\lambda$ is $N=20,604$. We deliberately select $K=2$ to keep the size of the
problem numerically manageable on a classical computer; and focus on
simulating the less-explored quantum measurement step, which is the
computational bottleneck of the HHL algorithm. More precisely, we compute $g$
classically as in (11), normalize it, then estimate the $\widehat{g}_{i}$ as
in (18). That is, the frequency numbers $n_{i}^{+}$ and $n_{i}^{-}$ are
obtained from a multinomial distribution with probabilities
$|g_{i}\pm\frac{1}{\sqrt{N}}|^{2}/4$. The number of measurements used in this
example is $n_{T}=500$, with $M=1000$ Monte Carlo samples. Note that the
present example does not really satisfy the situation $M\ll N$ because we keep
$N$ relatively small due to the complexity of the Stiefel neural network and
the limitations of classical computers. The condition $M\ll N$ applies to
large-scale deep learning applications where the number of parameters $N$ is
hundreds of thousands or higher and $M\approx 1000$ remains sufficient.

Figure 2: Lower bound over iterations for Stiefel neural network VB with a
classical simulation of quantum natural gradient. The smoothly-increasing
lower bound indicates a well-behaved training VB procedure. The underlying
dataset has 62 observations on a binary response variable and 2000 covariates.
The present example uses the first 100 covariates, with the likelihood based
on logistic regression. Stiefel neural network with $K=2$ ($N=20,604$)
provides a variational approximation of the posterior. The classical
simulation focuses on the quantum measurement step – the computational
bottleneck of the algorithm.

Figure 2 plots the lower bound estimates over the iterations. The smooth
increasing of the lower bound objective function indicates that the VB
algorithm with quantum natural gradient is converging and well behaved.

##  V Conclusion

We have proposed a computationally efficient regression-based approach to
estimating the natural gradient for VB, a step towards enabling a greater
range of Bayesian inference applications. We proved that VB with the
regression-based natural gradient is guaranteed to converge.

The regression-based approach supports quantum speed up of natural gradient
estimation. Quantum computing is the next-generation computing paradigm that
will transform computationally-intensive spheres of academic research and
industry. But quantum computers’ power can only apply to problems with
specific features, such as those where the handoff of data from a classical to
a quantum computer and back can be made economical, or those requiring
primarily linear transformation. We demonstrate that the regression-based
approach can make nearly arbitrary VB inference problems suitable for quantum
speedup.

One of the algorithmic tools we propose is the natural gradient readout with
quantum Gauss-Southwell rule, a quantum analogue of classical Gauss-Southwell
rule. An intriguing aspect of the quantum Gauss-Southwell method for natural
gradient is that it is an example of useful computation infeasible on a
classical computer but matter-of-course on a quantum computer. It provides a
taste for novel statistical and data science methods that will co-evolve with
the development of quantum computing [37].

Our work contributes to quantum machine learning: We propose a concrete
practical application of quantum linear systems algorithm using a quantum
regression method [53, 31, 54, 32, 55, 33, 56]. Our method works for a broad
class of variational models; it does not require a special-purpose learning
model, as in, e.g. [57, 58, 59, 60, 61, 62, 63, 64, 65, 66] and references
therein; quantum speedup for VB with special classes of models include [67,
68]. Our algorithm may help speed up quantum variational algorithms in in a
similar way it helps classical VB. A broader class of problems could benefit
from linearization technique, making them candidates for quantum speedup.

We also contribute to general advances in efficient quantum gradient descent
[69, 70, 42, 71, 72, 73, 60, 74, 75, 64, 76]. For example, Schuld _et al._
[70] propose the evaluation of the gradient for special-form quantum
variational distributions. We use an alternative approach to avoid encoding a
distribution function directly into a quantum state, which can be
prohibitively expensive, requiring at least $O(N)$ qubits [42]111Consider a
discretization with $d_{1}$ points for each dimension. It would require an
$O(h^{N})$ dimensional Hilbert space and $O(N\log d_{1})$ qubits..

The proposed algorithm adds to the growing list of quantum computing solutions
to real-life problems of high practical importance,including, for example,
modelling quantum chemistry for development of new materials, medicine,
agriculture, and energy ([78, 79, 80, 81, 82, 83, 84] and references therein);
modeling biological systems [85]; and option pricing [60, 86].

## Appendix

###### Proof of Theorem 1.

We first need the following lemma:

######  Lemma 1.

For every $t\geq 0$, ${\mathbb{E}}\|Y_{t}\|\leq\gamma_{t}C$ with $C>0$ some
finite constant.

###### Proof of Lemma 1.

We have that

$\displaystyle\|Y_{0}\|$ $\displaystyle=$
$\displaystyle\|\widehat{\nabla_{\lambda}^{\text{nat}}\tilde{\mathcal{L}}}(\lambda_{0})\|\leq\|\nabla_{\lambda}^{\text{nat}}\tilde{\mathcal{L}}(\lambda_{0})\|+\|\epsilon_{0}\|$
$\displaystyle\leq$
$\displaystyle\|I_{F}^{-1}(\lambda_{0})\|\|\nabla\tilde{\mathcal{L}}(\lambda_{0})\|+\|\epsilon_{0}\|$
$\displaystyle=$
$\displaystyle\sigma_{\max}(I_{F}^{-1}(\lambda_{0}))\|\nabla\tilde{\mathcal{L}}(\lambda_{0})\|+\|\epsilon_{0}\|.$

Hence,

$\displaystyle{\mathbb{E}}\|Y_{0}\|\leq\frac{b_{1}}{b_{2}}+b_{4}=C_{1}\leq\gamma_{0}C,$

with $C\geq C_{1}/\gamma_{0}$. Assume that
${\mathbb{E}}\|Y_{t}\|\leq\gamma_{t}C$ for some $t\geq 1$. Then,

$\displaystyle{\mathbb{E}}\|Y_{t+1}\|$ $\displaystyle\leq$
$\displaystyle\zeta_{t+1}{\mathbb{E}}\|Y_{t}\|$
$\displaystyle+\gamma_{t+1}{\mathbb{E}}(\|I_{F}^{-1}(\lambda_{t+1})\|\|\nabla\tilde{\mathcal{L}}(\lambda_{t+1})\|+\|\epsilon_{t+1}\|)$
$\displaystyle\leq$ $\displaystyle\zeta_{t+1}C\gamma_{t}+\gamma_{t+1}C_{1}$
$\displaystyle=$ $\displaystyle C\gamma_{t+1}(\omega\alpha_{t}+C_{1}/C).$

As $\omega\alpha_{t}<\omega<1$, with $C\geq C_{1}/(1-\omega)$ we have
$\omega\alpha_{t}+C_{1}/C\leq 1$ for all $t$, and hence,

${\mathbb{E}}\|Y_{t+1}\|\leq C\gamma_{t+1}.$

The proof is concluded by reduction with
$C=\max\\{C_{1}/\gamma_{0},C_{1}/(1-\omega)\\}$. ∎

We now prove Theorem 1. First,

$\displaystyle\langle Y_{t+1},\nabla\tilde{\mathcal{L}}(\lambda_{t+1})\rangle$
$\displaystyle=\zeta_{t+1}\langle
Y_{t},\nabla\tilde{\mathcal{L}}(\lambda_{t+1})\rangle$
$\displaystyle+\gamma_{t+1}\langle\nabla_{\lambda}^{\text{nat}}\tilde{\mathcal{L}}(\lambda_{t+1}),\nabla_{\lambda}\tilde{\mathcal{L}}(\lambda_{t+1})\rangle$
$\displaystyle+\gamma_{t+1}\langle\epsilon_{t+1},\nabla\tilde{\mathcal{L}}(\lambda_{t+1})\rangle$
(21)

By Lemma 1,

${\mathbb{E}}|\langle
Y_{t},\nabla\tilde{\mathcal{L}}(\lambda_{t+1})\rangle|\leq{\mathbb{E}}\|Y_{t}\|\|\nabla\tilde{\mathcal{L}}(\lambda_{t+1})\|\leq
b_{1}C\gamma_{t}.$

We have that,

${\mathbb{E}}\langle\epsilon_{t+1},\nabla\tilde{\mathcal{L}}(\lambda_{t+1})\rangle\leq
b_{1}{\mathbb{E}}\|\epsilon_{t+1}\|.$

Taking expectation of (Proof of Theorem 1.) results in

$\displaystyle{\mathbb{E}}\langle
Y_{t+1},\nabla\tilde{\mathcal{L}}(\lambda_{t+1})\rangle$
$\displaystyle=O(\alpha_{t}^{2})+\gamma_{t+1}{\mathbb{E}}\|\nabla_{\lambda}\tilde{\mathcal{L}}(\lambda_{t+1})\|^{2}_{I_{F}}$
$\displaystyle+O\big{(}\gamma_{t+1}{\mathbb{E}}\|\epsilon_{t+1}\|\big{)}.$
(22)

In the above, $\|v\|_{I_{F}}$ denotes the Fisher-Rao metric of vector $v$.
Using Taylor’s expansion,

$\displaystyle\tilde{\mathcal{L}}(\lambda_{t+1})$ $\displaystyle=$
$\displaystyle\tilde{\mathcal{L}}(\lambda_{t})-\langle
Y_{t},\nabla\tilde{\mathcal{L}}(\lambda_{t})\rangle+O(\|Y_{t}\|^{2}).$

By Lemma 1, taking the expectation of both sides and summing over
$t=0,1,...,n$,

$\sum_{t=0}^{n}{\mathbb{E}}\big{(}\langle
Y_{t},\nabla\tilde{\mathcal{L}}(\lambda_{t})\rangle\big{)}={\mathbb{E}}\tilde{\mathcal{L}}(\lambda_{0})-{\mathbb{E}}\tilde{\mathcal{L}}(\lambda_{n+1})+\sum_{t=0}^{n}O(\alpha_{t}^{2}).$
(23)

Comparing (Proof of Theorem 1.) with (23), and noting condition (A3) and that
$\sum\alpha_{t}^{2}<\infty$, we have that

$\sum_{t=1}^{\infty}\gamma_{t}{\mathbb{E}}\|\nabla_{\lambda}\tilde{\mathcal{L}}(\lambda_{t})\|^{2}_{I_{F}}<\infty.$

As

$\Big{(}\min_{t=1,...,n}{\mathbb{E}}\|\nabla_{\lambda}\tilde{\mathcal{L}}(\lambda_{t})\|^{2}_{I_{F}}\Big{)}\sum_{t=1}^{n}\gamma_{t}\leq\sum_{t=1}^{n}\gamma_{t}{\mathbb{E}}\|\nabla_{\lambda}\tilde{\mathcal{L}}(\lambda_{t})\|^{2}_{I_{F}},$

taking $n\to\infty$ and because $\sum_{t=1}^{n}\gamma_{t}\to\infty$, one must
have that

$\min_{t=1,...,n}{\mathbb{E}}\|\nabla_{\lambda}\tilde{\mathcal{L}}(\lambda_{t})\|^{2}_{I_{F}}\to
0,\;\;\;n\to\infty.$

From Assumption (A2),
$\|\nabla_{\lambda}\tilde{\mathcal{L}}(\lambda_{t})\|^{2}_{I_{F}}\geq\|\nabla_{\lambda}\tilde{\mathcal{L}}(\lambda_{t})\|^{2}/b_{3}$,
hence

$\min_{t=1,...,n}{\mathbb{E}}\|\nabla_{\lambda}\tilde{\mathcal{L}}(\lambda_{t})\|^{2}\to
0,\;\;\;n\to\infty,$

which proves (13). To proof (14), note that

$\displaystyle\tilde{\mathcal{L}}(\lambda^{*})$ $\displaystyle\geq$
$\displaystyle\tilde{\mathcal{L}}(\lambda_{t})+\nabla\tilde{\mathcal{L}}(\lambda_{t})^{\top}(\lambda^{*}-\lambda_{t})+\frac{L}{2}\|\lambda^{*}-\lambda_{t}\|^{2}$
$\displaystyle\geq$
$\displaystyle\tilde{\mathcal{L}}(\lambda_{t})-\|\nabla\tilde{\mathcal{L}}(\lambda_{t})\|\|\lambda^{*}-\lambda_{t}\|+\frac{L}{2}\|\lambda^{*}-\lambda_{t}\|^{2},$

where we have used the Cauchy-Schwarz inequality in the second inequality.
Then,

$\|\nabla\tilde{\mathcal{L}}(\lambda_{t})\|\|\lambda^{*}-\lambda_{t}\|-\frac{L}{2}\|\lambda^{*}-\lambda_{t}\|^{2}\geq\tilde{\mathcal{L}}(\lambda_{t})-\tilde{\mathcal{L}}(\lambda^{*})\geq
0$

which implies

$\|\lambda^{*}-\lambda_{t}\|\leq\frac{2}{L}\|\nabla\tilde{\mathcal{L}}(\lambda_{t})\|.$

Taking expectation, we have

$\min_{t=1,...,n}{\mathbb{E}}\|\lambda^{*}-\lambda_{t}\|\leq\frac{2}{L}\min_{t=1,...,n}{\mathbb{E}}\|\nabla\tilde{\mathcal{L}}(\lambda_{t})\|\to
0,\;\;n\to\infty.$

∎

## References

  * Hoffman  _et al._ [2013] M. D. Hoffman, D. M. Blei, C. Wang, and J. Paisley, Stochastic variational inference, Journal of Machine Learning Research 14, 1303 (2013).
  * Kingma and Welling [2014] D. Kingma and M. Welling, Auto-encoding Variational Bayes, in  _2nd International Conference on Learning Representations_ (2014).
  * Zhang  _et al._ [2018] C. Zhang, J. Bütepage, H. Kjellström, and S. Mandt, Advances in variational inference, IEEE transactions on pattern analysis and machine intelligence 41, 2008 (2018).
  * Amari [1998] S.-I. Amari, Natural gradient works efficiently in learning, Neural Computation 10, 251 (1998).
  * Sato [2001] M. Sato, Online model selection based on the variational Bayes, Neural Computation 13, 1649 (2001).
  * Tran  _et al._ [2017] M. Tran, D. Nott, and R. Kohn, Variational Bayes with intractable likelihood, Journal of Computational and Graphical Statistics 26, 873 (2017).
  * Mishkin  _et al._ [2018] A. Mishkin, F. Kunstner, D. Nielsen, M. Schmidt, and M. E. Khan, Slang: Fast structured covariance approximations for Bayesian deep learning with natural gradient, Neural Information Processing Systems (NIPS) (2018).
  * Martens [2020] J. Martens, New insights and perspectives on the natural gradient method, Journal of Machine Learning Research 21, 1 (2020).
  * Martens and Grosse [2015] J. Martens and R. Grosse, Optimizing neural networks with Kronecker-factored approximate curvature, in  _Proceedings of the 32nd International Conference on Machine Learning_ (2015) pp. 2408–2417.
  * Tran  _et al._ [2020] M.-N. Tran, N. Nguyen, D. Nott, and R. Kohn, Bayesian deep net glm and glmm, Journal of Computational and Graphical Statistics 29, 97 (2020).
  * Harrow  _et al._ [2009] A. W. Harrow, A. Hassidim, and S. Lloyd, Quantum algorithm for linear systems of equations, Physical Review Letters 103, 150502 (2009).
  * Ambainis [2012] A. Ambainis, Variable time amplitude amplification and quantum algorithms for linear algebra problems, in  _STACS’12 (29th Symposium on Theoretical Aspects of Computer Science)_ , Vol. 14 (LIPIcs, 2012) pp. 636–647.
  * Clader  _et al._ [2013] B. D. Clader, B. C. Jacobs, and C. R. Sprouse, Preconditioned quantum linear system algorithm, Physical Review Letters 110, 250504 (2013).
  * Kerenidis and Prakash [2016] I. Kerenidis and A. Prakash, Quantum recommendation systems, arXiv preprint arXiv:1603.08675 (2016).
  * Wossnig  _et al._ [2018] L. Wossnig, Z. Zhao, and A. Prakash, Quantum linear system algorithm for dense matrices, Physical Review Letters 120, 050502 (2018).
  * Childs  _et al._ [2017] A. M. Childs, R. Kothari, and R. D. Somma, Quantum algorithm for systems of linear equations with exponentially improved dependence on precision, SIAM Journal on Computing 46, 1920 (2017).
  * Subaşı  _et al._ [2019] Y. Subaşı, R. D. Somma, and D. Orsucci, Quantum algorithms for systems of linear equations inspired by adiabatic quantum computing, Physical Review Letters 122, 060504 (2019).
  * Gilyén  _et al._ [2019] A. Gilyén, Y. Su, G. H. Low, and N. Wiebe, Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics, in  _Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing_ (2019) pp. 193–204.
  * Aaronson [2015] S. Aaronson, Read the fine print, Nature Physics 11, 291 (2015).
  * Tran  _et al._ [2021a] M.-N. Tran, D. H. Nguyen, and D. Nguyen, Variational Bayes on manifolds, Statistics and Computing 31 (2021a).
  * Tran  _et al._ [2021b] M.-N. Tran, T.-N. Nguyen, and V.-H. Dao, A practical tutorial on variational bayes, arXiv:2103.01327 (2021b).
  * Goodfellow  _et al._ [2016] I. Goodfellow, Y. Bengio, and A. Courville,  _Deep Learning_ (MIT Press, 2016) http://www.deeplearningbook.org.
  * Kobak  _et al._ [2020] D. Kobak, J. Lomond, and B. Sanchez, The optimal ridge penalty for real-world high-dimensional data can be zero or negative due to the implicit ridge regularization, Journal of Machine Learning Research 21, 1 (2020).
  * Liang  _et al._ [2019] T. Liang, T. Poggio, A. Rakhlin, and J. Stokes, Fisher-rao metric, geometry, and complexity of neural networks, in  _The 22nd International Conference on Artificial Intelligence and Statistics_ (PMLR, 2019) pp. 888–896.
  * Bartlett  _et al._ [2020] P. L. Bartlett, P. M. Long, G. Lugosi, and A. Tsigler, Benign overfitting in linear regression, Proceedings of the National Academy of Sciences 117, 30063 (2020).
  * Dar  _et al._ [2021] Y. Dar, V. Muthukumar, and R. G. Baraniuk, A farewell to the bias-variance tradeoff? an overview of the theory of overparameterized machine learning, arXiv preprint arXiv:2109.02355 (2021).
  * Salimans  _et al._ [2013] T. Salimans, D. A. Knowles, _et al._ , Fixed-form variational posterior approximation through stochastic linear regression, Bayesian Analysis 8, 837 (2013).
  * Malago  _et al._ [2013] L. Malago, M. Matteucci, and G. Pistone, Natural gradient, fitness modelling and model selection: A unifying perspective, in  _2013 IEEE Congress on Evolutionary Computation_ (IEEE, 2013) pp. 486–493.
  * Note [1] Note that the conjugate gradient method [87, 88], sometimes used as a classical benchmark for quantum linear systems algorithms [13], does not apply to the natural gradient problem. The conjugate gradient method requires that the inverted matrix is symmetric and positive definite. The design matrix $T$ is not symmetric and its symmetrized version $\begin{pmatrix})0&T\\\ T^{\dagger}&0\end{pmatrix}$ is not positive definite.
  * Wiebe  _et al._ [2014] N. Wiebe, A. Kapoor, and K. M. Svore, Quantum deep learning, arXiv preprint arXiv:1412.3489 (2014).
  * Schuld  _et al._ [2016] M. Schuld, I. Sinayskiy, and F. Petruccione, Prediction by linear regression on a quantum computer, Physical Review A 94, 022342 (2016).
  * Wang [2017] G. Wang, Quantum algorithm for linear regression, Physical Review A 96, 012335 (2017).
  * Liu  _et al._ [2019] H.-L. Liu, C.-H. Yu, Y.-S. Wu, S.-J. Pan, S.-J. Qin, F. Gao, and Q.-Y. Wen, Quantum algorithm for logistic regression, arXiv preprint arXiv:1906.03834 (2019).
  * Yu  _et al._ [2021] C.-H. Yu, F. Gao, and Q. Wen, An improved quantum algorithm for ridge regression, IEEE Transactions on Knowledge and Data Engineering 33, 858 (2021).
  * Potok  _et al._ [2021] T. Potok _et al._ , Adiabatic quantum linear regression, Scientific Reports 11, 1 (2021).
  * Schuld and Petruccione [2018] M. Schuld and F. Petruccione,  _Supervised learning with quantum computers_ (Springer, 2018).
  * Lopatnikova  _et al._ [2021] A. Lopatnikova, M.-N. Tran, and S. A. Sisson, An introduction to quantum computing for statisticians and data scientists, arXiv preprint arXiv:2112.06587 (2021).
  * Berry  _et al._ [2015] D. W. Berry, A. M. Childs, and R. Kothari, Hamiltonian simulation with nearly optimal dependence on all parameters, in  _2015 IEEE 56th Annual Symposium on Foundations of Computer Science_ (IEEE, 2015) pp. 792–809.
  * Aharonov and Ta-Shma [2003] D. Aharonov and A. Ta-Shma, Adiabatic quantum state generation and statistical zero knowledge, in  _Proceedings of the 35th Annual ACM Symposium on Theory of Computing_ (2003) pp. 20–29.
  * Gilyén  _et al._ [2018] A. Gilyén, S. Lloyd, and E. Tang, Quantum-inspired low-rank stochastic regression with logarithmic dependence on the dimension, arXiv preprint arXiv:1811.04909 (2018).
  * Frieze  _et al._ [2004] A. Frieze, R. Kannan, and S. Vempala, Fast monte-carlo algorithms for finding low-rank approximations, Journal of the ACM (JACM) 51, 1025 (2004).
  * Rebentrost  _et al._ [2019] P. Rebentrost, M. Schuld, L. Wossnig, F. Petruccione, and S. Lloyd, Quantum gradient descent and Newton’s method for constrained polynomial optimization, New Journal of Physics 21, 073023 (2019).
  * Note [2] Bausch [89] proposes a recurrent method with post-selection in the context of quantum Recurrent Neural Networks and argues that, for approximate post-selection it is possible to reduce this overhead. Another possible approach is to run the iterative algorithm without post-selection for $T$ steps and then amplify the desired end state using Quantum Amplitude Amplification (QAA) [50]. Even with these improvements, the exponential dependence on the number of iterations $T/z$ (where $z=2$ for QAA method) remains.
  * Note [3] The algorithm leverages the density matrix exponentiation method of Lloyd _et al._ [90].
  * Giovannetti  _et al._ [2008a] V. Giovannetti, S. Lloyd, and L. Maccone, Quantum random access memory, Physical Review Letters 100, 160501 (2008a).
  * Giovannetti  _et al._ [2008b] V. Giovannetti, S. Lloyd, and L. Maccone, Architectures for a quantum random access memory, Physical Review A 78, 052310 (2008b).
  * Hann  _et al._ [2021] C. T. Hann, G. Lee, S. Girvin, and L. Jiang, Resilience of quantum random access memory to generic noise, PRX Quantum 2, 020311 (2021).
  * Nielsen and Chuang [2002] M. A. Nielsen and I. Chuang, Quantum computation and quantum information (2002).
  * Nutini  _et al._ [2015] J. Nutini, M. Schmidt, I. Laradji, M. Friedlander, and H. Koepke, Coordinate descent converges faster with the gauss-southwell rule than random selection, in  _International Conference on Machine Learning_ (PMLR, 2015) pp. 1632–1641.
  * Brassard  _et al._ [2002] G. Brassard, P. Hoyer, M. Mosca, and A. Tapp, Quantum amplitude amplification and estimation, Contemporary Mathematics 305, 53 (2002).
  * Tang [2018] E. Tang, Quantum-inspired classical algorithms for principal component analysis and supervised clustering, arXiv preprint arXiv:1811.00414 (2018).
  * Papamakarios  _et al._ [2021] G. Papamakarios, E. Nalisnick, D. Rezende, S. Mohamed, and B. Lakshminarayanan, Normalizing flows for probabilistic modeling and inference, Journal of Machine Learning Research 22, 1 (2021).
  * Wiebe  _et al._ [2012] N. Wiebe, D. Braun, and S. Lloyd, Quantum algorithm for data fitting, Physical Review Letters 109, 050505 (2012).
  * Duan  _et al._ [2017] B. Duan, J. Yuan, Y. Liu, and D. Li, Quantum algorithm for support matrix machines, Physical Review A 96, 032301 (2017).
  * Wang  _et al._ [2019] D. Wang, O. Higgott, and S. Brierley, Accelerated variational quantum eigensolver, Physical Review Letters 122, 140504 (2019).
  * Kerenidis and Prakash [2020] I. Kerenidis and A. Prakash, Quantum gradient descent for linear systems and least squares, Physical Review A 101, 022316 (2020).
  * Low  _et al._ [2014] G. H. Low, T. J. Yoder, and I. L. Chuang, Quantum inference on Bayesian networks, Physical Review A 89, 062315 (2014).
  * Rebentrost  _et al._ [2014] P. Rebentrost, M. Mohseni, and S. Lloyd, Quantum support vector machine for big data classification, Physical Review Letters 113, 130503 (2014).
  * McClean  _et al._ [2016] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, The theory of variational hybrid quantum-classical algorithms, New Journal of Physics 18, 023023 (2016).
  * Rebentrost  _et al._ [2018] P. Rebentrost, T. R. Bromley, C. Weedbrook, and S. Lloyd, Quantum Hopfield neural network, Physical Review A 98, 042308 (2018).
  * Harrow and Napp [2019] A. Harrow and J. Napp, Low-depth gradient measurements can improve convergence in variational hybrid quantum-classical algorithms, arXiv preprint arXiv:1901.05374 (2019).
  * McArdle  _et al._ [2019] S. McArdle, T. Jones, S. Endo, Y. Li, S. C. Benjamin, and X. Yuan, Variational ansatz-based quantum simulation of imaginary time evolution, npj Quantum Information 5, 1 (2019).
  * Schuld and Killoran [2019] M. Schuld and N. Killoran, Quantum machine learning in feature Hilbert spaces, Physical Review Letters 122, 040504 (2019).
  * Schuld  _et al._ [2020] M. Schuld, A. Bocharov, K. M. Svore, and N. Wiebe, Circuit-centric quantum classifiers, Physical Review A 101, 032308 (2020).
  * Abbas  _et al._ [2020] A. Abbas, D. Sutter, C. Zoufal, A. Lucchi, A. Figalli, and S. Woerner, The power of quantum neural networks, arXiv preprint arXiv:2011.00027 (2020).
  * Park and Kastoryano [2020] C.-Y. Park and M. J. Kastoryano, Geometry of learning neural quantum states, Physical Review Research 2, 023232 (2020).
  * Miyahara and Sughiyama [2018] H. Miyahara and Y. Sughiyama, Quantum extension of variational Bayes inference, Physical Review A 98, 022330 (2018).
  * Zhao  _et al._ [2019] Z. Zhao, A. Pozas-Kerstjens, P. Rebentrost, and P. Wittek, Bayesian deep learning on a quantum computer, Quantum Machine Intelligence 1, 41 (2019).
  * Jordan [2005] S. P. Jordan, Fast quantum algorithm for numerical gradient estimation, Physical Review Letters 95, 050501 (2005).
  * Schuld  _et al._ [2019] M. Schuld, V. Bergholm, C. Gogolin, J. Izaac, and N. Killoran, Evaluating analytic gradients on quantum hardware, Physical Review A 99, 032331 (2019).
  * Stokes  _et al._ [2020] J. Stokes, J. Izaac, N. Killoran, and G. Carleo, Quantum natural gradient, Quantum 4, 269 (2020).
  * Sweke  _et al._ [2020] R. Sweke, F. Wilde, J. J. Meyer, M. Schuld, P. K. Fährmann, B. Meynard-Piganeau, and J. Eisert, Stochastic gradient descent for hybrid quantum-classical optimization, Quantum 4, 314 (2020).
  * Gacon  _et al._ [2021] J. Gacon, C. Zoufal, G. Carleo, and S. Woerner, Simultaneous perturbation stochastic approximation of the quantum fisher information, arXiv preprint arXiv:2103.09232 (2021).
  * Mitarai  _et al._ [2018] K. Mitarai, M. Negoro, M. Kitagawa, and K. Fujii, Quantum circuit learning, Physical Review A 98, 032309 (2018).
  * Zoufal  _et al._ [2019] C. Zoufal, A. Lucchi, and S. Woerner, Quantum generative adversarial networks for learning and loading random distributions, npj Quantum Information 5, 1 (2019).
  * Zoufal  _et al._ [2020] C. Zoufal, A. Lucchi, and S. Woerner, Variational quantum Boltzmann machines, arXiv preprint arXiv:2006.06004 (2020).
  * Note [4] Consider a discretization with $d_{1}$ points for each dimension. It would require an $O(h^{N})$ dimensional Hilbert space and $O(N\log d_{1})$ qubits.
  * Reiher  _et al._ [2017] M. Reiher, N. Wiebe, K. M. Svore, D. Wecker, and M. Troyer, Elucidating reaction mechanisms on quantum computers, Proceedings of the National Academy of Sciences 114, 7555 (2017).
  * Aspuru-Guzik  _et al._ [2018] A. Aspuru-Guzik, R. Lindh, and M. Reiher, The matter simulation (r) evolution, ACS Central Science 4, 144 (2018).
  * Cao  _et al._ [2018] Y. Cao, J. Romero, and A. Aspuru-Guzik, Potential of quantum computing for drug discovery, IBM Journal of Research and Development 62, 6 (2018).
  * McArdle  _et al._ [2020] S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Quantum computational chemistry, Reviews of Modern Physics 92, 015003 (2020).
  * Babbush  _et al._ [2018] R. Babbush, N. Wiebe, J. McClean, J. McClain, H. Neven, and G. K.-L. Chan, Low-depth quantum simulation of materials, Physical Review X 8, 011044 (2018).
  * Macridin  _et al._ [2018] A. Macridin, P. Spentzouris, J. Amundson, and R. Harnik, Electron-phonon systems on a universal quantum computer, Physical Review Letters 121, 110504 (2018).
  * Ajagekar and You [2019] A. Ajagekar and F. You, Quantum computing for energy systems optimization: Challenges and opportunities, Energy 179, 76 (2019).
  * Robert  _et al._ [2021] A. Robert, P. K. Barkoutsos, S. Woerner, and I. Tavernelli, Resource-efficient quantum algorithm for protein folding, npj Quantum Information 7, 1 (2021).
  * Stamatopoulos  _et al._ [2020] N. Stamatopoulos, D. J. Egger, Y. Sun, C. Zoufal, R. Iten, N. Shen, and S. Woerner, Option pricing using quantum computers, Quantum 4, 291 (2020).
  * Shewchuk  _et al._ [1994] J. R. Shewchuk _et al._ , An introduction to the conjugate gradient method without the agonizing pain (1994).
  * Saad [2003] Y. Saad,  _Iterative methods for sparse linear systems_ (SIAM, 2003).
  * Bausch [2020] J. Bausch, Recurrent quantum neural networks, Advances in Neural Information Processing Systems 33 (2020).
  * Lloyd  _et al._ [2014] S. Lloyd, M. Mohseni, and P. Rebentrost, Quantum principal component analysis, Nature Physics 10, 631 (2014).