Low-depth gradient measurements can improve convergence in variational hybrid quantum-classical algorithms

Introduction

where $A_{j}$ is the Hermitian operator which generates pulse $j$. This is the form of variational state most commonly encountered in the literature on variational algorithms, and is also motivated theoretically [MRBAG16, YRS+17, BJ18]. Note that in the presence of noise, the parameterized family of states which may be prepared will actually consist of mixed states. We only consider the noiseless case in this paper. It may also be the case that there are more pulses than independent parameters. For instance, one may impose a constraint like $\bm{\uptheta}_{i}=\bm{\uptheta}_{j}$. We assume for simplicity that the parameters are independent, but comment on how our results could be easily extended to this case.

It is assumed that the quantum device is controlled by a classical “outer loop”, and the quantum device is used only for preparing simple quantum states and making simple measurements. The classical outer loop uses this measurement information to perform a classical optimization of some function $f(\bm{\uptheta})$ over the feasible set $\mathcal{X}$, where the objective function $f(\bm{\uptheta})$ is induced by some Hermitian objective observable $H$, via the relation $f(\bm{\uptheta}):=\expectationvalue{H}{\bm{\uptheta}}$. Algorithms of this family have been proposed in the context of quantum simulation (e.g. variational quantum eigensolvers [PMS+14, WHT15]), combinatorial optimization (e.g. QAOA [FGG14]), and machine learning (e.g. quantum classifiers [FN18, MNKF18, SK18, SBSW18, HCT+18]).

As a simple example, in a simulation context, $H$ could be some physical Hamiltonian for which we want to approximately obtain the ground state energy. If the true ground state (or a state close to the true ground state) belongs to the parameterized family $\{\ket{\bm{\uptheta}}\}_{\bm{\uptheta}}$ for $\bm{\uptheta}\in\mathcal{X}$, approximately minimizing $f(\bm{\uptheta})$ yields an approximation to the ground state energy. The typical way a variational algorithm would obtain information about $f(\bm{\uptheta})$ with little quantum resources is by expanding $f(\bm{\uptheta})=\expectationvalue{H}{\bm{\uptheta}}=\sum_{i=1}^{m}\alpha_{i}\expectationvalue{P_{i}}{\bm{\uptheta}}$ where $P_{i}$ are tensor products of Pauli operators (recall that any operator may be expanded in this way), or are some other observables which may be easily measured. In this paper, we assume the $P_{i}$ are products of Pauli operators. Assuming it is possible to easily measure these operators, one can estimate $f(\bm{\uptheta})$ by estimating each $\expectationvalue{P_{i}}{\bm{\uptheta}}$ separately and combining the results according to the coefficients $\alpha_{i}$. Of course, due to the randomness of the measurement outcomes, many preparations of $\ket{\bm{\uptheta}}$ and measurements may be required to obtain a good estimate of $f(\bm{\uptheta})$.

This effect occurs generally in variational hybrid algorithms: the randomness of the quantum measurement outcomes translates into the classical outer loop only having stochastic access to the objective function $f$. That is, at point $\bm{\uptheta}$ in parameter space, it cannot directly observe the function value $f(\bm{\uptheta})$, but rather some random variable whose expectation value is $f(\bm{\uptheta})$. Numerical simulations that overlook this fact may be misleading. For instance, letting $\epsilon$ denote the optimization error, some problems that admit $\log(1/\epsilon)$ convergence rates given noiseless access to the objective function values admit $\text{poly}(1/\epsilon)$ convergence rates in the stochastic setting [Bub15].

But even worse for prospects of optimization, the resulting classical stochastic optimization problem will generally be complicated and nonconvex, and hence may be intractable. However, one can hope for heuristics which find a reasonable approximate solution. Furthermore, if the algorithm is in a convex vicinity of an optimum, or local optimum, algorithms like stochastic gradient descent which are known to converge for convex problems may converge to the local optimum despite the problem being globally nonconvex. For some proposed applications of variational algorithms one desires a very precise solution, so it is likely that much time will be spent converging in the vicinity of an optimum and this situation may be especially relevant. In this paper, we focus on this latter scenario of convergence within a convex region containing a local optimum, either assuming or proving that this case applies.

In the usual formulation of variational algorithms, the classical outer loop is assumed to take some number of quantum measurements at a point $\bm{\uptheta}$ in parameter space in order to approximate the objective function value $f(\bm{\uptheta})$ at that point, and then perform an optimization based on these values. However, one can imagine more complicated algorithms which, instead of taking measurements to estimate $f(\bm{\uptheta})$ at point $\bm{\uptheta}$, take measurements corresponding to some other property of the optimization problem. Indeed, a natural alternative choice is to take measurements corresponding to $\nabla f(\bm{\uptheta})$, the gradient of the objective function at point $\bm{\uptheta}$. Such a strategy was proposed in the context of combinatorial optimization [GS17], quantum chemistry [RBM+18], and machine learning [MNKF18, SBSW18, FN18]. Similar ideas were also proposed in the context of implementing Hamiltonian evolution in a low-depth variational setting [LB17]. Low-depth procedures for directly measurement gradients in variational algorithms typically require only marginally greater quantum circuit depth than that required for measuring the objective function. Some such methods are reviewed and extended in [BIS+18, SBG+18].

However, it was not clear that such gradient-measurement based strategies could confer any advantage over objective-measurement based strategies. For one, note that it is possible to obtain an estimate of $\nabla f(\bm{\uptheta})$ using only estimates of the objective function $f(\bm{\uptheta})$. To see how, note that for small $\epsilon$,

2 Summary of results

In order to rigorously prove a lower bound on the number of quantum measurements required for zeroth-order variational optimization, it is convenient to introduce a black-box formalism. To motivate the introduction of a black-box formalism, note that any variational algorithm can be simulated by a purely classical algorithm that makes zero quantum measurements. Of course, the time complexity of the classical simulation may be exponential in the problem size.

One encounters a similar problem in the classical setting of trying to quantify the complexity of optimization. The objective function to be minimized could be extremely complicated and difficult to study analytically (for instance, it could correspond to the output of some complicated algorithm) or otherwise inaccessible. A black-box model was therefore developed for the study of convex optimization [NY83] and remains popular in current research in the field. In this setting, the function to be optimized is encoded in an oracle, and we define the query complexity of an algorithm for optimizing the function to be the number of calls made to the oracle. The algorithm may be promised that the objective function has certain properties, but is not given an exact description of the function. This black-box formalism provides a natural and general setting for proving bounds in convex optimization.

It should be noted that, when variational algorithms are used in practice, commuting terms in the Pauli decomposition of $H$ can be measured on a single trial state. To simplify the analysis, our black-box model does not take advantage of this possible speedup. Hence, the query cost in our model corresponds to the number of products of Pauli operators measured, rather than the number of state preparations required, which is a related quantity but could be lower. (We comment on how taking this into account would affect our bounds for the toy model we study in Table 2.)

We consider a restriction of the variational problem to a convex region of parameter space $\mathcal{X}$ on which the objective function $f(\bm{\uptheta})$ is assumed to be convex. We import known results from the stochastic optimization literature to obtain convergence rates for various optimization strategies and for various assumptions about $f(\bm{\uptheta})$, such as strong convexity. We derive upper bounds on the query cost for strategies that use analytic gradient measurements in conjunction with stochastic gradient descent or stochastic mirror descent with an $l_{1}$ setup. The upper bounds are functions of the dimension of parameter space $p$, precision $\epsilon$, geometry of the feasible set, and coefficients in the Pauli expansions of the objective observable and pulse generators.

Stochastic mirror decent (see Appendix A.2, or e.g. [NJLS09, JN11, Bub15] for more detailed reviews), which we will abbreviate as SMD, can be thought of as a generalization of SGD to non-Euclidean spaces. As motivation for why SMD might be relevant in our setting, note that the 1-norm of a parameter vector has a very natural interpretation. Namely, since $e^{-iA_{j}\bm{\uptheta}_{j}/2}$ is essentially the unitary evolution generated by $A_{j}/2$ for time $\bm{\uptheta}_{j}$, $\|\bm{\uptheta}\|_{1}:=\sum_{i}|\bm{\uptheta}_{i}|$ may be interpreted as the total amount of time that the starting state $\ket{\Psi}$ is evolved for to reach the trial state $\ket{\bm{\uptheta}}:=e^{-iA_{p}\bm{\uptheta}_{p}/2}\cdots e^{-iA_{1}\bm{\uptheta}_{1}/2}\ket{\Psi}$, and $\|\bm{\uptheta}_{i}-\bm{\uptheta}_{j}\|_{1}$ is associated with the amount of time for which the associated pulse sequences differ. On the other hand, SGD is appropriate for Euclidean geometriesa. Taking our norm to be $\|\cdot\|_{1}$ instead of $\|\cdot\|_{2}$ and using a suitable version of SMD yields nearly quadratically better scaling with respect to the dimension of parameter space $p$ in some settings, as compared with SGD. In other settings, SGD outperforms SMD. We record upper bounds based on both strategies. It is an open problem to understand what sort of values the parameters in the upper bounds typically take in practice, and relatedly, whether a Euclidean geometry or some other geometry is more appropriate.

For comparison, we also record a rigorous upper bound for zeroth-order strategies by applying two of the best known zeroth-order upper bounds [FKM05, AFH+11] from the stochastic optimization literature. Our results on general upper bounds are displayed in Table 1. It should be noted that these zeroth-order upper bounds are the best rigorous zeroth-order upper bounds that we are aware of, but it is likely that other derivative-free algorithms would significantly outperform these bounds in practice in many instances. For instance, methods based on trust regions and surrogate models [CSV09] often perform very well in practice, despite not necessarily having strong theoretical guarantees.

This point is an example of a more general limitation that applies to all of the upper bounds we report, relating to the difference between theoretical and empirical results. For one, these are upper bounds for convergence that apply when the algorithm has a trusted region which it knows contains a local optimum, and which the objective function is convex over. Furthermore, while strong convexity of the objective function in the domain can greatly improve the convergence, the rigorous upper bounds that exploit this property apply when the algorithm has a good estimate of this strong convexity parameter. In other language, the upper bounds apply in a promise setting, where the algorithm is promised that a certain convex subset of parameter space contains the optimum, the objective function is convex in this domain, and any other relevant parameters take certain values. However, this may very well not be the case in practice. Furthermore, these are worst-case upper bounds, and may be outperformed in practice.

All of these challenges are well-known in machine learning tasks such as deep learning, in which one may want to converge to a local optimum of some complicated nonconvex optimization problem. In practice, one can do hyperparameter optimization, which would involve yet another classical outer loop varying over the parameters used to define the optimization algorithm itself. Parameters could also be set adaptively (e.g. [KB14]). Another possibility is that the algorithm could try to construct a surrogate model for the objective function which is valid in some region, and estimate relevant parameters from the surrogate model. Unfortunately, such methods are often not backed up by strong theoretical results, even when the practical performance is very good. While there may be a sizable gap between theory and practice, we hope that the theoretical convergence upper bounds will nonetheless be useful for guiding practical implementations and expectations.

Finally, we point out that these upper bounds do not explicitly depend on the number of terms in the Pauli expansion of the objective observable or pulse generators, but rather on sums of coefficients in the expansion. This is a desirable feature for applications such as quantum chemistry, where it is often the case that the electronic structure Hamiltonian is written as a sum of a very large number of terms, many of which have very small norm.

After establishing an oracular setting and recording some upper bounds, for a given parameter $\epsilon>0$ we define a certain class $\mathcal{H}_{n}^{\epsilon}$ of simple 1-local objective observables on $n$ qubits which we use to demonstrate a separation in query complexity between variational algorithms which only make zeroth-order queries to the oracle and those which make first-order queries. The optima of the observables in $\mathcal{H}^{\epsilon}_{n}$ are $O(\epsilon)$-close to each other in objective function value, in the sense that for all $H,H^{\prime}\in\mathcal{H}^{\epsilon}_{n}$, if $\ket{\psi}$ is an optimum (i.e. ground state) of $H$, then $\expectationvalue{H^{\prime}}{\psi}-\lambda_{\min}(H^{\prime})\leq O(\epsilon)$, where $\lambda_{\min}(H^{\prime})$ is the smallest eigenvalue of $H^{\prime}$.

We show that for any precision parameter $0<\epsilon<\Theta(n)$, any zeroth-order variational algorithm which optimizes any observable $H\in\mathcal{H}^{\epsilon}_{n}$ to precision $\epsilon$ and only queries the oracle with states in the vicinity of the optimum of $H$ must make at least $\Omega(n^{3}/\epsilon^{2})$ queries. Here, the “vicinity of the optimum” is essentially the set of states that are $O(\epsilon)$-close to optimal in objective function value. On the other hand, we show that after making a good choice of variational ansatz, an SGD algorithm that performs analytic gradient measurements to get gradient estimates optimizes any objective observable in the family to expected precision $\epsilon$ with only $O(n^{2}/\epsilon)$ queries of states in the vicinity of the optimum.

Suppose $\mathcal{A}$ is a variational algorithm that only makes zeroth-order queries to $\mathcal{O}_{H}$, only queries states in the vicinity of the optima of $\mathcal{H}^{\epsilon}_{n}$, and for any $H\in\mathcal{H}^{\epsilon}_{n}$ that is realized, outputs a description of a state whose expected objective function value is $\epsilon$-close to $\lambda_{\min}(H)$. Then $\mathcal{A}$ must make $\Omega\quantity(\frac{n^{3}}{\epsilon^{2}})$ queries.

There exists a variational algorithm that only makes first-order queries to $\mathcal{O}_{H}$, only queries states in the vicinity of the optima of $\mathcal{H}^{\epsilon}_{n}$, and for any $H\in\mathcal{H}^{\epsilon}_{n}$, makes $O\quantity(\frac{n^{2}}{\epsilon})$ queries and outputs a description of a state whose expected objective function value is $\epsilon$-close to $\lambda_{\min}(H)$. An algorithm that achieves this rate is a simple stochastic gradient descent strategy.

Suppose $\mathcal{A}$ is a variational algorithm that may make queries of any order to $\mathcal{O}_{H}$, may query the oracle with any state, and for any $H\in\mathcal{H}^{\epsilon}_{n}$ outputs a description of a state whose expected objective function value is $\epsilon$-close to the optimum. Then $\mathcal{A}$ makes at least $\Omega\quantity(\frac{n^{2}}{\epsilon})$ queries.

We summarize the above results in Table 2, where for comparison we also include a rigorous zeroth-order upper bound based on the results of [FKM05] and [AFH+11].

3 Related work

In this paper, we are primarily interested in the low-depth setting. The methods we consider for measuring the gradient yield an unbiased, but possibly very noisy estimate of the gradient in low depth. An alternative approach for measuring the gradient in variational algorithms was recently proposed in [GAW17], which builds on Jordan’s gradient measurement algorithm [Jor05]. Their algorithm offers significantly better performance for obtaining precise estimates of the gradient, but also requires significantly more quantum resources, with coherence time requirements increasing with the desired precision.

A lower bound for a class of derivative-free stochastic convex optimization problems was shown in [JNR12]. Our separation result in Section 5 is similar in spirit to their result, and the proof strategies share similarities. However, our setting of variational hybrid algorithms is very different from theirs, preventing their result from being ported to variational quantum algorithms. In particular, we are interested specifically in stochastic optimization problems induced by a quantum observable and variational ansatz. Furthermore, in our setting the classical outer loop’s optimization problem is not fixed; different variational ansätze induce different optimization problems. Our lower bound for zeroth-order algorithms takes this extra freedom into account, applying for any choice of variational ansatz (and also allowing the algorithm to change ansatz over the course of the optimization). Our proof strategy for the zeroth-order lower bound also borrows some techniques from [AWBR09], which showed lower bounds for classes of first-order stochastic optimization problems. In turn, these techniques are inspired by methods in statistical minimax and learning theory.

4 Organization

In Section 2, we state the conventions we adhere to and record known results from stochastic convex optimization that we later use. In Section 3, we introduce a black-box setting for variational algorithms and define the oracle $\mathcal{O}_{H}$ encoding an objective observable $H$. In Section 4, we present general upper bounds on the query cost of optimization for different algorithms and different assumptions on the objective function. In Section 5, we define a parameterized class of variational optimization problems on $n$ qubits $\mathcal{H}_{n}^{\epsilon}$ and use this class of problems to prove a query complexity separation between zeroth-order and first-order optimization strategies. We conclude and mention some open questions in Section 6. Appendix A contains some relevant background on first-order stochastic convex optimization algorithms.

Preliminaries

We will assume throughout that variational states are parameterized according to an ansatz of the form

Given a parameterization $\Theta$ and a feasible set $\mathcal{X}$, the classical objective function to be minimized is induced by some Hermitian operator $H$, which we refer to as the objective observable. In particular, the classical objective function is given by

In the context of variational algorithms, it is assumed that the quantum device is capable of measuring some subset of quantum observables. We assume that the set of observables which may be measured is the set of all Pauli operators.

We collect notation and parameters in Table 3.

2 Requisite results about stochastic convex optimization

We will obtain upper bounds for variational algorithms in convex regions by combining well known classical convergence results with sampling strategies for estimating the gradient. Here, we record the classical optimization results we will need. Background on stochastic gradient descent and stochastic mirror descent may be found in Appendix A (see e.g. [NJLS09, JN11, Bub15] for more thorough reviews). First, we define strong convexity.

For $\lambda>0$, the real-valued function $f$ is $\lambda$-strongly convex with respect to norm $\|\cdot\|$ on some convex domain $\mathcal{X}$ if $\forall\mathbf{x},\mathbf{y}\in\mathcal{X}$,

Note that a twice-differentiable function is $\lambda$-strongly convex with respect to the $2$-norm if all of the eigenvalues of the Hessian matrix at each point in the domain are at least $\lambda$. More generally, $f$ is $\lambda$-strongly convex at $\mathbf{\mathbf{x}}$ w.r.t. an arbitrary norm $\|\cdot\|$ if $\mathbf{h}^{\top}\nabla^{2}f(\mathbf{x})\mathbf{h}\geq\lambda\|\mathbf{h}\|^{2}$ for all $\mathbf{h}$, where $\nabla^{2}f(\mathbf{x})$ is the Hessian of $f$ at $\mathbf{x}$. In contrast, $f$ is convex if the Hessians are merely positive semidefinite. Intuitively, if $f$ is strongly convex, then it is lower bounded by a quadratic function. Strong convexity can often be used to accelerate optimization [Bub15].

In this section, we record known upper bounds for optimizing convex functions given access to noisy, unbiased gradient information. For the first two results below, we follow the presentation of the review on algorithms for convex optimization [Bub15].

Assume $\mathcal{X}$ is contained in a Euclidean ball of radius $R_{2}$ and $f$ is convex on $\mathcal{X}$. Then projected SGD with fixed step size $\eta=\frac{R_{2}}{G_{2}}\sqrt{\frac{2}{T}}$ satisfies

where $\mathbf{x}_{1}$ is the starting point, and the algorithm visits points $\mathbf{x}_{1},\dots,\mathbf{x}_{T}$.

Assume $f$ is $\lambda_{2}$-strongly convex on $\mathcal{X}$ with respect to $\|\cdot\|_{2}$. Then SGD with step size $\eta_{s}=\frac{2}{\lambda_{2}(s+1)}$ at iteration $s$ satisfies

where $\mathbf{x}_{1}$ is the starting point, and the algorithm visits points $\mathbf{x}_{1},\dots,\mathbf{x}_{T}$.

Assume $\mathcal{X}$ is contained in a 1-ball of radius $R_{1}$, and $f$ is convex on $\mathcal{X}$. Then stochastic mirror descent with an appropriate $l_{1}$ setup and and step size $\eta=\frac{R_{1}}{G_{\infty}}\sqrt{\frac{2}{T}}$, satisfies

where $\mathbf{x}_{1}$ is the starting point, and the algorithm queries points $\mathbf{x}_{1},\dots,\mathbf{x}_{T}$.

Assume $f$ is $\lambda_{1}$-strongly convex on $\mathcal{X}$ with respect to norm $\|\cdot\|_{1}$. Then a certain SMD-like algorithm, running for $T$ iterations, outputs a (random) vector $\bar{\mathbf{x}}$ such that

2.2 Upper bounds for stochastic zeroth-order (derivative-free) optimization

Compared to stochastic first-order optimization, less is known about rigorous upper bounds for stochastic zeroth-order optimization. The few rigorous upper bounds for zeroth-order optimization that are known [FKM05, AD10, AFH+11, JNR12, Sha13] are usually weaker than their first-order counterparts. The upper bounds we use in this paper are from [FKM05] and [AFH+11], in which the authors prove rigorous upper bounds for stochastic zeroth-order convex optimization in which the expected error in objective function value converges to zero like $\sqrt{\frac{p^{2}}{T}}$ and $\sqrt{\frac{p^{32}}{T}}$, respectively, where $T$ is the number of iterations. We record their results below, adapted for our purposes. Note that in these papers, the authors state the results in terms of online optimization. However, it is straightforward to convert these into results for stochastic optimization. Similar adaptations of these same results are also noted in [JNR12] and [Sha13].

There also exists an algorithm [AFH+11] that makes $T$ queries and outputs $\bar{\mathbf{x}}$ such that

Note that this implies that $\min\quantity(O\quantity(\frac{p^{2}E^{2}R_{2}^{2}(L+E/r_{2})^{2}}{\epsilon^{4}}),\widetilde{O}(\frac{p^{32}E^{2}}{\epsilon^{2}}))$ queries are needed to optimize to expected precision $\epsilon$. Other rigorous upper bounds for derivative-free stochastic convex optimization exist [AD10, JNR12, Sha13], which apply in settings which are less applicable to this paper, or to specific families of functions.

Black-box formulation

One of our main goals is to prove a rigorous lower bound on the number of quantum measurements required to minimize an objective function $f(\bm{\uptheta})=\expectationvalue{H}{\bm{\uptheta}}$ to within some precision $\epsilon$. However, examining the formulation of this question already reveals a subtlety. Namely, a classical computer can simulate the quantum part of a variational algorithm in (in general) exponential time, which implies that any variational algorithm can be simulated with a variational algorithm which makes zero quantum measurements.

Of course, in practice it would generally be intractable for a classical computer to simulate a variational algorithm. In fact, it is known that the existence of an efficient classical algorithm for sampling from the output distribution of the commonly-considered variational algorithm QAOA would imply a collapse of the polynomial hierarchy [FH16]. We therefore seek a formulation of the problem which better captures the behavior of realistic algorithms. One way to do this is to strip away the classical outer loop’s knowledge of the specific objective observable $H$ that it is trying to optimize by encoding the observable in the black box, and only giving the outer loop the black box $\mathcal{O}_{H}$ and a promise that $H$ belongs to some particular family $\mathcal{H}$. In practice, this essentially means that the black-box picture is applicable when the classical optimization algorithm that the outer loop runs does not depend on the details of the objective observable $H$ itself (but could depend on the family $\mathcal{H}$), but is rather some more general-purpose algorithm like gradient descent, Nelder-Mead, SPSA, BFGS, etc. We are not aware of any proposed variational algorithm that is not in this class. Note that, while the classical optimization component of a variational algorithm does not exploit detailed structure of the objective observable in current proposals, the variational ansatz sometimes does (e.g. [PMS+14, FGG14, WHT15]). For example, in QAOA, the ansatz involves pulses of the form $e^{iH\bm{\uptheta}_{j}}$ where $H$ is the objective observable. For such cases, our query upper bounds are still fully applicable. Our lower bound for zeroth-order variational optimization is ansatz-independent in the sense that for any ansatz the algorithm chooses, the lower bound still holds. Hence, the zeroth-order lower bound is still fully applicable in the setting where the ansatz may be a function of $H$. However, there is a subtlety that if the algorithm is promised that the ansatz is a certain function of $H$, it could use this information to learn the ground state of $H$ with fewer queries, and the lower bound may no longer apply. Essentially, this means that the zeroth-order lower bound applies for zeroth-order algorithms which do not cleverly exploit the dependence of the variational ansatz on the objective observable. This includes all “general-purpose” zeroth-order methods such as SPSA, Nelder-Mead, and gradient descent with gradients estimated via finite-differences, to name a few.

In the remainder of this section, we define our black-box formulation of variational algorithms. Specifically, we define an oracle $\mathcal{O}_{H}$ encoding the objective observable $H$. This definition allows us to talk about the query complexity of optimizing some family of objective observables $\mathcal{H}$. The classical outer loop is given as input an oracle $\mathcal{O}_{H}$ under the promise that $H\in\mathcal{H}$, and attempts to find an approximate optimum using as few queries as possible.

If the parameterization $\Theta$ is clear from context, we may not explicitly note that $\Theta$ is provided as input to the oracle. If we speak of querying the oracle with a state $\ket{\psi}$, we mean querying the oracle with a parameterization and parameter vector that describe the state $\ket{\psi}$. In the remainder of this section, we first define how the oracle behaves for zeroth-order queries. We then briefly review how to measure gradients (and higher-order derivatives) in low depth in variational algorithms, and then define the behavior of the oracle upon first- and higher-order queries.

Let $H$ be some objective observable. Decompose $H$ into a linear combination of $m$ products of Pauli operators as $H=\sum_{i=1}^{m}\alpha_{i}P_{i}$ where $\alpha_{i}>0$. (The coefficients may all be assumed to be positive by absorbing the phase into the operator.) Now, defining the normalization factor $E:=\sum_{i}\alpha_{i}$ and the probability distribution $p_{i}:=\alpha_{i}/E$, we may write

From this expression, it is clear that by sampling index $i$ with probability $p_{i}$, measuring $P_{i}$ with respect to the state $\ket{\bm{\uptheta}}$, and then multiplying the outcome by $E$ we obtain an unbiased estimator for $f(\bm{\uptheta})=\expectationvalue{H}{\bm{\uptheta}}$. Furthermore, since the measurement outcome of $P_{i}$ is either $+1$ or $-1$, the output of this estimator is $\pm E$-valued. It is also clear that $|f(\bm{\uptheta})|\leq E$ for all $\bm{\uptheta}$. We define the behavior of the sampling oracle $\mathcal{O}_{H}$ for zeroth-order sampling to be essentially the above process.

Let $H=E\sum_{i=1}^{m}p_{i}P_{i}$ be a decomposition of an objective observable as above, where $E>0$ and $p_{i}$ is a probability distribution. Given as input a parameterization $\Theta$, a parameter $\bm{\uptheta}$, and an empty coordinate multiset $S=\varnothing$, the oracle behaves as follows. It internally prepares $\ket{\bm{\uptheta}}$ and measures the observable $P_{i}$ with probability proportional to $p_{i}$. It then multiplies the outcome by $E$ and outputs the resulting $\pm E$-valued estimator.

2 Analytic gradient measurements

Variational algorithms typically aim to estimate the objective function $f(\bm{\uptheta}):=\expectationvalue{H}{\bm{\uptheta}}$ at some point $\bm{\uptheta}$ in parameter space, and use this information along with previous estimates of $f$ to propose a new point $\bm{\uptheta}^{\prime}$ in parameter space. In this case, the classical outer loop essentially has a stochastic zeroth-order oracle for the objective function.

Some recent works [GS17, MNKF18, RBM+18, SBSW18, SBG+18] have instead suggested a different optimization strategy, in which one directly extracts information about the gradient of the objective function $f(\bm{\uptheta})$ by measuring corresponding quantum observables. Quantum measurements of this type that correspond to estimates of the gradient of the objective function are often referred to as analytic gradient measurements. In this section, we review these strategies.

For notational convenience, we define $U_{i}:=e^{-iA_{i}\bm{\uptheta}_{i}/2}$ to be the unitary corresponding to pulse $i$, and for $i\leq j$ we define $U_{i:j}:=e^{-iA_{j}\bm{\uptheta}_{j}/2}\cdots e^{-iA_{i}\bm{\uptheta}_{i}/2}$ to be the sequence of pulses from $i$ through $j$, inclusive. Note that in using this notation we are hiding the dependence on $\bm{\uptheta}$ for visual clarity. Recall that the objective function corresponding to objective observable $H$ is given by

It is straightforward to calculate the following relation via the chain rule applied to the above expression:

We now describe how the above quantity could be measured in a variational algorithm. Denote the Pauli decomposition of $A_{j}$ as $A_{j}=\sum_{k=1}^{n_{j}}\beta^{(j)}_{k}Q^{(j)}_{k}$ where $Q^{(j)}_{k}$ are products of Pauli operators. As in the previous sections, denote the Pauli decomposition of $H$ as $H=\sum_{i=1}^{m}\alpha_{i}P_{i}$. Then by linearity we can rewrite the above derivative as

Now, we can obtain an unbiased estimator for $\imaginary\expectationvalue{U_{1:j}^{\dagger}Q^{(j)}_{k}U_{(j+1):p}^{\dagger}P_{l}U_{1:p}}{\Psi}$ via a (generalized) Hadamard test. In particular, the following procedure may be used for estimating $\imaginary\expectationvalue{U_{1:j}^{\dagger}Q^{(j)}_{k}U_{(j+1):p}^{\dagger}P_{l}U_{1:p}}{\Psi}$.

Hadamard test for estimating $-\imaginary\expectationvalue{U_{1:j}^{\dagger}Q^{(j)}_{k}U_{(j+1):p}^{\dagger}P_{l}U_{1:p}}{\Psi}$ 1. Initialize Register $A$ in the qubit state $\ket{+}_{A}$. Initialize Register $B$ in the state $\ket{\Psi}_{B}$. 2. Apply $U_{1:j}$ to Register $B$. 3. Apply a Controlled-$Q^{(j)}_{k}$ gate to Register $B$, controlled on Register $A$. 4. Apply $U_{(j+1):p}$ to Register $B$. 5. Apply a Controlled-$P_{l}$ gate to Register $B$, controlled on Register $A$. 6. Measure the Pauli $Y$ operator on Register $A$. The above procedure yields a $\pm 1$-valued unbiased estimator for $-\imaginary\expectationvalue{U_{1:j}^{\dagger}Q^{(j)}_{k}U_{(j+1):p}^{\dagger}P_{l}U_{1:p}}{\Psi}$, requiring one quantum measurement.

Algorithm 1: Generalized Hadamard test [EAO+02, LB17, GS17, RBM+18]

Hence, one may estimate $\nabla f(\bm{\uptheta})$ by expanding the derivatives as above, and then estimating each term of the expansion using Algorithm 1. Alternative methods of analyticly measuring derivatives are described in [MNKF18, SBG+18], which require similar quantum resources to the scheme we just described (but do not necessarily require controlled-Pauli gates). We note that throughout this paper, one could estimate gradients using a strategy based on these methods instead, and the results would be essentially unchanged.

We now describe an equivalent way of understanding analytic gradients. Observe that

Hence, if we define the Hermitian operators

and we define $\vec{G}:=(G_{1},\dots,G_{p})^{\top}$, then we may write $\nabla f(\bm{\uptheta})=\expectationvalue{\vec{G}}{\bm{\uptheta}}$. An alternative commutator expression for the derivatives was noted in [MBS+18].

There are cases in which one may want to impose a constraint that some parameters are always equal. For example, this situation occurs for the “Hamiltonian variational” ansatz proposed in [WHT15]. Hence, one could have a $p$-pulse ansatz but a smaller number of independent variational parameters. We note that this situation is easily addressed within the framework of this paper. For example, consider the case in which $\bm{\uptheta}_{i}$ is constrained to always equal $\bm{\uptheta}_{j}$, i.e. $\bm{\uptheta}_{i}=\bm{\uptheta}_{j}:=\xi$. It is straightforward to show by linearity that $\partialderivative{f}{\xi}(\bm{\uptheta})=\expectationvalue{(G_{i}+G_{j})}{\bm{\uptheta}}$, where $G_{i}$ and $G_{j}$ are defined as above. Hence, $\partialderivative{f}{\xi}(\bm{\uptheta})$ may be estimated via Algorithm 1 just as in the unconstrained case. For simplicity, we assume that there are no such constraints on the parameters. However, all results in this paper can be easily generalized to work with such constraints via this observation.

3 First-order sampling

In the previous section, we described how information about the derivatives of the objective function can be extracted in low depth using a generalized Hadamard test. In this section, we describe a specific estimator of a derivative of the objective function which requires one Pauli measurement. We will use this estimator to define the behavior of the oracle $\mathcal{O}_{H}$ upon a first-order query.

As in the previous section, denote the Pauli expansion of $H$ as $H=\sum_{i=1}^{m}\alpha_{i}P_{i}$ and the Pauli expansion of $A_{j}$ as $A_{j}=\sum_{k=1}^{n_{j}}\beta^{(j)}_{k}Q^{(j)}_{k}$, where all $\alpha$ and $\beta$ coefficients are positive real numbers. Then we may write $\partialderivative{f}{\bm{\uptheta}_{j}}(\bm{\uptheta})$ as the following expansion:

We now rewrite this expansion as a certain expectation value, similarly to what we did in the definition of zeroth-order sampling. First, we observe that some of the commutators in the expansion may trivially be zero, if the operators $U_{(j+1):p}Q^{(j)}_{k}U_{(j+1):p}^{\dagger}$ and $P_{l}$ act nontrivially on disjoint sets of qubits. This will often be the case in the toy model we analyze in Section 5. Removing terms that are trivially zero will improve convergence in our optimization algorithms. To this end, we define a new set of coefficients:

where $\text{qubits}(U_{(j+1):p}Q^{(j)}_{k}U_{(j+1):p}^{\dagger})$ denotes the set of qubits on which $U_{(j+1):p}Q^{(j)}_{k}U_{(j+1):p}^{\dagger}$ acts nontrivially, after removing pulses which trivially commute through $Q^{(j)}_{k}$ and cancel the corresponding inverse pulse. We define the associated normalization factors

and probability distributions $q^{(j)}_{kl}:=\frac{1}{\Gamma_{j}}\gamma^{(j)}_{kl}$ over the indices $k$ and $l$, where $j$ is considered fixed. Note that we have the bound $\Gamma_{j}\leq EB_{j}$ where $B_{j}:=\sum_{k=1}^{n_{j}}\beta^{(j)}_{k}$. Equipped with these definitions, we may write

It is straightforward to see that $|\partialderivative{f}{\bm{\uptheta}_{j}}(\bm{\uptheta})|\leq\Gamma_{j}$ for all $\bm{\uptheta}$. Given the above representation of $\partialderivative{f}{\bm{\uptheta}_{j}}(\bm{\uptheta})$, it is clear that the following procedure provides an unbiased estimator for $\partialderivative{f}{\bm{\uptheta}_{j}}(\bm{\uptheta})$ which requires a single measurement.

An unbiased one-measurement estimator for $\partialderivative{f}{\bm{\uptheta}_{j}}(\bm{\uptheta})$. 1. Sample $(K,L)$ from the distribution $q^{(j)}_{KL}$ as defined above. 2. Use a Hadamard test (Algorithm 1) to obtain a one-measurement unbiased estimate of $\expectationvalue{\frac{i}{2}\commutator{U_{(j+1):p}Q^{(j)}_{K}U_{(j+1):p}^{\dagger}}{P_{L}}}{\bm{\uptheta}}=-\imaginary\expectationvalue{U_{1:j}^{\dagger}Q^{(j)}_{K}U^{\dagger}_{(j+1):p}P_{L}U_{1:p}}{\Psi}$. 3. Multiply the resulting number by $\Gamma_{j}$. The estimator for $\partialderivative{f}{\bm{\uptheta}_{j}}(\bm{\uptheta})$ described above is $\pm\Gamma_{j}$-valued.

Algorithm 2: unbiased, one-measurement estimator for $\partialderivative{f}{\bm{\uptheta}_{j}}(\bm{\uptheta})$.

Motivated by these derivative-estimating procedures, we now define the first-order behavior of the oracle $\mathcal{O}_{H}$.

Let $H$ denote an objective observable. Upon input of parameterization $\Theta$, parameter $\bm{\uptheta}$, and a coordinate multiset $S=\{j\}$ for some $j\in[p]$, the oracle internally prepares the state $\ket{\bm{\uptheta}}$ and runs Algorithm 2 above. It outputs the resulting $\pm\Gamma_{j}$-valued estimator for $\partialderivative{f}{\bm{\uptheta}_{j}}(\bm{\uptheta})$.

4 Higher-order sampling

The sampling procedure we have described above for obtaining unbiased estimates of derivatives in low depth generalizes to higher-order derivatives. In this section, we outline how the procedure would work. Start by recalling the derivative operators we derived above:

where, since $G_{j}$ is independent of $\bm{\uptheta}_{k}$, this result follows from arguments identical to those we used to derive the expression for $G_{j}$. Also, note that from the original definition $\expectationvalue{H}{\bm{\uptheta}}:=\expectationvalue{e^{iA_{1}\bm{\uptheta}_{1}/2}\cdots e^{iA_{p}\bm{\uptheta}_{p}/2}He^{-iA_{p}\bm{\uptheta}_{p}/2}\cdots e^{-iA_{1}\bm{\uptheta}_{1}/2}}{\Psi}$, it is clear that $\partialderivative{f}{\bm{\uptheta}_{k}}{\bm{\uptheta}_{j}}=\partialderivative{f}{\bm{\uptheta}_{j}}{\bm{\uptheta}_{k}}$. We therefore have, for $k\leq j$,

To see how to estimate this in low depth, note that we have

From the above expression, we see how to generalize the first-order sampling procedure to higher orders. To obtain an unbiased estimate of $\partialderivative{f}{\bm{\uptheta}_{k}}{\bm{\uptheta}_{j}}$ with a single measurement, first expand $A_{k}$, $A_{j}$, and $H$ as linear combinations of products of Paulis. In turn, this yields an expansion of $\partialderivative{f}{\bm{\uptheta}_{k}}{\bm{\uptheta}_{j}}$ as a linear combination of real parts of inner products of states that are acted on with pulses and Paulis. For a one-measurement estimator, randomly choose one of these inner products with probability proportional to the magnitude its coefficient, and then get an unbiased estimate of the inner product by performing a Hadamard test, similarly to what we described for the first-order case. Note that, in the second-order case (or more generally for the even-order case), the Hadamard test will involve an $X$-basis measurement instead of $Y$-basis measurement, since a real part is being estimated.

5 Query complexity in the black-box formalism

Having defined the sampling oracle $\mathcal{O}_{H}$, we may now quantify the cost of an algorithm by the number of queries it makes to the oracle. The general setup for a variational optimization problem in the black-box setting is that the classical “outer loop” is promised that the objective observable $H$ to be minimized belongs to a family $\mathcal{H}$ of observables, and is given access to the sampling oracle $\mathcal{O}_{H}$. Note that, from the perspective of the outer loop, the problem of minimizing the objective function is a purely classical black-box optimization problem since it gives classical input to the oracle and receives classical output. We now formalize the notion of the “error” associated with some variational algorithm $\mathcal{A}$ for optimizing a family $\mathcal{H}$ of objective observables.

Let $\mathcal{H}$ denote a set of objective observables, and $\mathcal{A}$ be a (possibly randomized) classical algorithm which has access to a sampling oracle $\mathcal{O}_{H}$ for some $H\in\mathcal{H}$ and outputs a description of a quantum state $\ket{\psi}$. Then the optimization error of $\mathcal{A}$ with respect to $\mathcal{H}$, $\operatorname*{Err}(\mathcal{A},\mathcal{H})$, is defined to be

where the expectation is over the possible randomness of the output state $\ket{\psi}$.

In other words, $\operatorname*{Err}(\mathcal{A},\mathcal{H})$ is the worst-case expected error in objective function value that $\mathcal{A}$ makes over all objective observables in the set $\mathcal{H}$. We now make a few more definitions that will be convenient later.

Define the $\delta$-optimum of an observable $H$ to be the set of all states $\ket{\psi}$ such that $\expectationvalue{H}{\psi}-\lambda_{\min}(H)\leq\delta$. Define the $\delta$-optimum of a set of observables $\mathcal{H}$ to be the union of the $\delta$-optima of each observable in the set. We say a black-box algorithm $\mathcal{A}$ is a $\delta$-vicinity algorithm for $\mathcal{H}$ if it only queries the black box with descriptions of states that are in the $\delta$-optimum of $\mathcal{H}$.

General upper bounds for variational algorithms in a convex region

In this section, we give general upper bounds on the query cost of variational algorithms in a region where the objective function is convex. This amounts to applying the known upper bounds for stochastic convex optimization from Section 2.2 to the setting in which estimates of the objective function, or derivatives of the objective function, come from the oracle specified in Section 3 (which is easy to implement in low depth in practice). Note that the oracle returns estimates of partial derivatives w.r.t. specific components. However, there are multiple ways of using these derivative estimates to construct a gradient estimator. We describe two such estimators. The first is designed to be used with SGD, and the second is designed to be used with SMD with an $l_{1}$ setup.

For the remainder of this section, fix some objective observable with Pauli expansion $H=\sum_{i=1}^{m}\alpha_{i}P_{i}$ and some parameterization $\Theta$ whose pulse generators $A_{j}$ have Pauli expansions $A_{j}=\sum_{k=1}^{n_{j}}\beta^{(j)}_{k}Q^{(j)}_{k}$. As in Section 3, define $E:=\sum_{i=1}^{m}\alpha_{i}$, $B_{j}=\sum_{k=1}^{n_{j}}\beta^{(j)}_{k}$. Define $\Gamma_{j}$ to be the normalization factor associated with coordinate $j$ as defined in Section 3 (see also Table 3). We collect these $\Gamma_{j}$ into a vector as

First, we specify some unbiased estimators for the gradient that we will use. The estimator of Algorithm 3 is based on $l_{1}$ sampling and designed with the goal in mind of achieving a smaller $2$-norm of the estimator and will be used in conjunction with SGD. The estimator of Algorithm 4 is based on $l_{2}$ sampling and designed with the goal of achieving a smaller $\infty$-norm of the estimator and will be used in conjunction with SMD. The latter estimator also requires a mild assumption on the $\infty$-norm of the objective function. The estimators also differ in their number of samples: the former estimator uses a single sample while the latter could be called a “mini-batch” estimator which uses an asymptotically growing number of samples. We first define the two estimators, and then prove their correctness and bound them in the subsequent lemmas.

Algorithm 3: $l_{1}$-sampling estimator for $\nabla f(\bm{\uptheta})$.

Algorithm 4: $l_{2}$-sampling estimator for $\nabla f(\bm{\uptheta})$.

for $t\geq 0$. Using this bound with $t=\frac{\|\vec{\Gamma}\|_{2}}{\sqrt{2p}}$ and recalling $N_{j}=\left\lceil p\frac{\Gamma_{j}^{2}}{\|\vec{\Gamma}\|_{2}^{2}}\ln\quantity(4p^{2}\frac{\|\vec{\Gamma}\|_{\infty}^{2}}{\|\vec{\Gamma}\|_{2}^{2}})\right\rceil$ yields

By the union bound, the probability that $|\hat{G}_{j}-(\nabla f)_{j}|\geq\frac{\|\vec{\Gamma}\|_{2}}{\sqrt{2p}}$ for some $j$ is upper bounded by $\frac{\|\vec{\Gamma}\|_{2}^{2}}{2p\|\vec{\Gamma}\|_{\infty}^{2}}$. If this event occurs, then we only have the trivial upper bound $\|\hat{\mathbf{g}}\|_{\infty}^{2}\leq\|\vec{\Gamma}\|_{\infty}^{2}$. Conditioned on this “bad” event not occurring, we have the bound $\|\hat{\mathbf{g}}\|_{\infty}^{2}\leq 2\frac{\|\vec{\Gamma}\|_{2}^{2}}{p}$, where we used the assumption that $|(\nabla f)_{j}|\leq\frac{\|\vec{\Gamma}\|_{2}}{\sqrt{2p}}$ for all $j$. It follows that

2 Upper bounds

Fix an objective observable $H$, parameterization $\Theta$, and a closed, convex feasible set $\mathcal{X}$. Define $f(\bm{\uptheta}):=\expectationvalue{H}{\bm{\uptheta}}$, and define $\vec{\Gamma}$ as above (see also Table 3). Let $\bm{\uptheta}^{*}$ denote a minimizer of $f(\bm{\uptheta})$ on $\mathcal{X}$.

Use the 1-query estimator of Algorithm 3 for $\nabla f(\bm{\uptheta})$ in conjunction with Theorem 2.1. ∎

Use the 1-query estimator of Algorithm 3 for $\nabla f(\bm{\uptheta})$ in conjunction with Theorem 2.2. ∎

Use the estimator of Algorithm 4 for $\nabla f(\bm{\uptheta})$ in conjunction with Theorem 2.4. ∎

For comparison, we also present an upper bound for the case in which we only make zeroth-order queries to $\mathcal{O}_{H}$.

Note that the outputs of zeroth-order queries to $\mathcal{O}_{H}$ have magnitude $E$, and apply Theorem 2.5. ∎

3 When is SMD superior to SGD?

In Section 1.2, we gave intuition for why we might hope that using the 1-norm instead of 2-norm and using SMD with an $l_{1}$ setup instead of SGD might be beneficial in some cases. In particular, we noted that the 1-norm of a parameter vector $\bm{\uptheta}$ has a natural interpretation as the duration of evolution from the starting state $\ket{\Psi}$ to the trial state associated with $\bm{\uptheta}$, $\ket{\bm{\uptheta}}$. The $l_{1}$-distance between two parameter vectors may be interpreted as the amount of time for while the two associated pulse sequences differ.

Comparing the upper bounds from the previous section, we see that where SGD has a factor of $\|\vec{\Gamma}\|_{1}^{2}$, SMD with an $l_{1}$ setup has instead a factor of $\|\vec{\Gamma}\|_{2}^{2}$. Note that $\|\vec{\Gamma}\|_{2}^{2}$ is never larger than $\|\vec{\Gamma}\|_{1}^{2}$, and in fact can a factor of $p$ smaller. Consider for example the case in which $\Gamma_{1}\approx\Gamma_{2}\approx\cdots\approx\Gamma_{p}$, which may be a realistic scenario in practice. In this case, we have $\|\vec{\Gamma}\|_{2}^{2}\approx p\Gamma_{1}^{2}$ for the SMD bound, whereas we have $\|\vec{\Gamma}\|_{1}^{2}\approx p^{2}\Gamma_{1}^{2}$ for the SGD bound, which is quadratically worse in the dimension of parameter space.

On the other hand, where the SGD bounds involve a factor of $R_{2}^{2}$, the SMD bounds involve a factor of $R_{1}^{2}$. $R_{2}^{2}$ is never larger than $R_{1}^{2}$, and can be significantly smaller. This could be the case when, for example, the feasible set $\mathcal{X}$ is a Euclidean ball. On the other hand, if $\mathcal{X}$ is a 1-ball, then $R_{1}=R_{2}$, and SMD could potentially achieve substantially better performance than SGD due to the $\|\vec{\Gamma}\|_{2}^{2}$ versus $\|\vec{\Gamma}\|_{1}^{2}$ discrepancy.

Another consideration is the issue of strong convexity. As is evident from the above bounds, the presence of strong convexity can substantially accelerate the optimization. SGD can take advantage of strong convexity w.r.t. the 2-norm, but SMD in the $l_{1}$ setup measures strong convexity w.r.t. the 1-norm, and in fact it is straightforward to show that the strong convexity parameters are related by $\lambda_{1}\leq\lambda_{2}\leq p\lambda_{1}$. In the toy problem we analyze in Section 5, we have $\lambda_{2}=\Theta(1)$ but $\lambda_{1}=\Theta(1/n)$ where $n$ is the number of qubits. However, $\|\vec{\Gamma}\|_{1}^{2}=\Theta(n^{2})$ while $\|\vec{\Gamma}\|_{2}^{2}=\Theta(n)$, so up to log factors and constants, SGD and SMD achieve the same asymptotic convergence rate for this toy model.

In conclusion, it is not clear from the upper bounds in the previous section or from the toy model we study in Section 5 whether SGD or SMD with an $l_{1}$ setup would typically achieve better upper bounds in practice. It is an interesting problem for future work to understand whether an $l_{2}$ (Euclidean) setup or an $l_{1}$ setup is usually more natural for variational algorithms.

Oracle separation between zeroth-order and first-order optimization strategies for variational algorithms

In this section, we prove a separation between algorithms which make only zeroth-order queries to the sampling oracle, and those which make first-order queries to the sampling oracle, within the vicinity of the global optimum. This separation is proven with respect to a certain simple parameterized family $\mathcal{H}_{n}^{\epsilon}$ of objective observables on $n$ qubits. The optima of the observables in $\mathcal{H}^{\epsilon}_{n}$ are $O(\epsilon)$ close to each other, in the sense that for any $H,H^{\prime}\in\mathcal{H}^{\epsilon}_{n}$, the ground state of $H$ is an $O(\epsilon)$ optimum of $H^{\prime}$. Precisely, we will prove the following.

For any $n\geq 15$ and $\epsilon\leq 0.01n$, let $\mathcal{A}$ be any zeroth-order, $100\epsilon$-vicinity algorithm for the family $\mathcal{H}^{\epsilon}_{n}$ that makes $T$ queries to the oracle. Then, if $\operatorname*{Err}(\mathcal{A},\mathcal{H}^{\epsilon}_{n})\leq\epsilon$, it must hold that $T\geq\Omega\quantity(\frac{n^{3}}{\epsilon^{2}})$ where the implicit factor is some fixed constant.

On the other hand, we prove that this same class of variational problems can be optimized substantially faster if the algorithm makes first-order queries to the oracle, as quantified in the following theorem. In fact, the algorithm that achieves this convergence rate is a simple stochastic gradient descent strategy. Hence, not only is the query complexity much better in this case, but the classical algorithm achieving this query complexity can be implemented efficiently. For comparison, we also obtain a zeroth-order upper bound for this class of problems using the algorithms of [FKM05] and [AFH+11]. Our first-order upper bound is given in the following theorem.

For any $\epsilon\leq 0.01n$, there exists a first-order, $100\epsilon$-vicinity algorithm $\mathcal{A}$ for the family $\mathcal{H}^{\epsilon}_{n}$ that makes $O\quantity(\frac{n^{2}}{\epsilon})$ queries and achieves an error $\operatorname*{Err}(\mathcal{A},\mathcal{H}^{\epsilon}_{n})\leq\epsilon$. Moreover, $\mathcal{A}$ is a simple stochastic gradient descent algorithm.

For any $n\geq 15$ and $\epsilon\leq 0.01n$, suppose $\mathcal{A}$ is an algorithm that makes $T$ queries and satisfies $\operatorname*{Err}(\mathcal{A},\mathcal{H}^{\epsilon}_{n})\leq\epsilon$. Then $T\geq\Omega\quantity(\frac{n^{2}}{\epsilon})$.

Since this lower bound is achieved (up to a possible constant factor) by the upper bound of SGD, we see that SGD is essentially optimal among all black-box strategies for optimizing $\mathcal{H}_{n}^{\epsilon}$.

The subset of objective observables we consider are perturbed around a very simple 1-local Hamiltonian.

Intuitively, for a fixed small parameter $\delta$, the set of $2^{n}$ observables $\{H^{\delta}_{v}\}_{v}$ are perturbed around $H^{0}=-\frac{1}{\sqrt{2}}\sum_{i=1}^{n}\quantity(X_{i}+Z_{i})$. The parameter $\delta$ characterizes the strength of the perturbation, and the binary vector $v$ encodes the direction of the perturbation. It is straightforward to see that the ground state of $H^{0}$ is $\ket{\pi/4}^{\otimes n}$, where we have defined $\ket{\pi/4}:=\cos(\pi/8)\ket{0}+\sin(\pi/8)\ket{1}$. Geometrically, the state $\ket{\pi/4}$ corresponds to the pure qubit state with polarization $\frac{1}{\sqrt{2}}(\hat{x}+\hat{z})$. In the remainder of this section, we record some facts about these Hamiltonians, and define some quantities.

First, note that we may write $H^{\delta}_{v}=-\sum_{i=1}^{n}\hat{n}^{v_{i}\delta}\cdot\vec{\sigma}_{i}$ where $\hat{n}^{v_{i}\delta}=\quantity(\sin(\frac{\pi}{4}+v_{i}\delta),0,\cos(\frac{\pi}{4}+v_{i}\delta))$ and $\vec{\sigma}_{i}$ is the vector of Pauli operators acting on qubit $i$. We may now read off $\lambda_{\min}(H^{\delta}_{v})=-n$, and the associated eigenvector is

Next, we calculate the expectation value of $H^{\delta}_{v}$ with respect to any quantum state on $n$ qubits.

Suppose $\rho$ is a quantum state such that the polarization of $\rho_{i}$, the reduced state of $\rho$ on qubit $i$, is $\vec{r}_{i}$. Then $\tr\quantity[H_{v}^{\delta}\rho]=-\sum_{i=1}^{n}\vec{r}_{i}\cdot\hat{n}^{\delta v_{i}}$.

Finally, we define the set $\mathcal{H}_{n}^{\epsilon}$ which we will prove the separation with respect to. To do so, we first define a bias parameter $\delta(\epsilon)$ associated with the precision parameter $\epsilon$.

For a given “precision parameter” $\epsilon$, define the associated “bias parameter”

Now, we define $\mathcal{H}^{\epsilon}_{n}$ to be the set of such observables with bias parameter $\delta(\epsilon)$.

$\mathcal{H}_{n}^{\epsilon}:=\{H^{\delta(\epsilon)}_{v}\,:\,\forall v\in\{-1,1\}^{n}\}$.

For the remainder of the paper, we often hide the dependence of $\delta$ on $\epsilon$ for notational simplicity, and simply write $\delta$ where we implicitly mean $\delta(\epsilon)$. Note that our constraint $\epsilon\leq 0.01n$ implies $\delta<0.7$.

In this section, we prove Theorem 5.1. Our proof strategy for the lower bound is to reduce a statistical learning problem to the optimization problem, and then lower bound the number of oracle calls required to solve the learning problem. Precisely, we will take an appropriate subset $\mathcal{M}^{\epsilon}_{n}\subset\mathcal{H}^{\epsilon}_{n}$, parameterized by some subset $\mathcal{V}$ of the $n$-dimensional hypercube $\{-1,+1\}^{n}$. That is, we will have $\mathcal{M}^{\epsilon}_{n}=\{H^{\delta(\epsilon)}_{v}\,:\,v\in\mathcal{V}\}$ where $\mathcal{V}\subset\{-1,1\}^{n}$ will be strategically chosen. We prove that, if there exists an algorithm $\mathcal{A}$ that satisfies $\operatorname*{Err}(\mathcal{A},\mathcal{M}^{\epsilon}_{n})\leq\epsilon$, then the same algorithm could be used to identify the hidden parameter $v\in\mathcal{V}$ associated with the objective observable $H^{\delta}_{v}\in\mathcal{M}^{\epsilon}_{n}$. By employing information theoretic methods, we will lower bound the number of oracle calls required to identify the parameter $v$, which in turn lower bounds the number of calls required to optimize to precision $\epsilon$.

Our proof in some parts adapts techniques from [AWBR09] and [JNR12], which lower bound the query cost of certain convex first-order and derivative-free optimization problems. These results in turn draw on methods from statistical minimax and learning theory.

We begin by defining, for fixed $\epsilon$, a subset $\mathcal{M}^{\epsilon}_{n}\subset\mathcal{H}_{n}^{\epsilon}$ of objective observables that are well-separated, in the sense that if a state is close to the optimal of $H_{v}^{\delta}\in\mathcal{M}^{\epsilon}_{n}$, then it must be far from the optimal of $H_{v^{\prime}}^{\delta}\in\mathcal{M}^{\epsilon}_{n}$ for any other parameter $v^{\prime}$. We make this precise below.

We make use of the following classical fact about packings of the hypercube (see for example [Gun11] for a simple proof).

There exists a subset $\mathcal{V}$ of the $n$-dimensional hypercube $\{-1,1\}^{n}$ of size $|\mathcal{V}|\geq e^{n/8}$ such that, if $\Delta(v,v^{\prime})$ denotes the Hamming distance between $v$ and $v^{\prime}$,

for all $v\neq v^{\prime}$ with $v,v^{\prime}\in\mathcal{V}$.

Fix $\mathcal{V}$ to be such a subset of $\{-1,1\}^{n}$, and define $\mathcal{M}^{\epsilon}_{n}:=\{H^{\delta}_{v}\,:\,v\in\mathcal{V}\}$. The Hamming distance provides a natural distance measure between points of the hypercube. We now define a notion of distance $d$ between objective observables $H_{v}^{\delta}$ and $H_{v^{\prime}}^{\delta}$. Intuitively, if $d(v,v^{\prime})$ is large, then a state that is close to the optimal of $H_{v}^{\delta}$ cannot be close to the optimal of $H_{v^{\prime}}^{\delta}$.

For $v,v^{\prime}\in\{-1,1\}^{n}$, we define the semimetric

where the minimization is over all normalized pure states on $n$ qubits.

Note that $\lambda_{\min}(H^{\delta}_{v})$ is simply $-n$, but we oftentimes write $\lambda_{\min}(H^{\delta}_{v})$ for clarity. We now define a packing parameter $\beta$ which quantifies how packed the subset $\mathcal{V}$ is, with respect to the semimetric $d$.

The packing parameter $\beta$ corresponding to the above subset $\mathcal{V}\subset\{-1,1\}^{n}$ and semimetric $d$ on the hypercube is defined to be

Suppose that for some state $\ket{\psi}$ and parameter $v\in\mathcal{V}$, $\expectationvalue{H_{v}^{\delta}}{\psi}-\lambda_{\min}(H^{\delta}_{v})\leq\beta/3$. Then for all $v^{\prime}\neq v$ with $v^{\prime}\in\mathcal{V}$, $\expectationvalue{H_{v^{\prime}}^{\delta}}{\psi}-\lambda_{\min}(H^{\delta}_{v^{\prime}})>\beta/3$.

Suppose there exists some parameter $v^{\prime}\in\mathcal{V}$, $v^{\prime}\neq v$ for which $\expectationvalue{H_{v^{\prime}}^{\delta}}{\psi}-\lambda_{\min}(H^{\delta}_{v^{\prime}})\leq\beta/3$. From Definition 5.4, this implies that $d(v,v^{\prime})\leq 2\beta/3$, which contradicts the assumption that $\beta$ is the packing parameter. ∎

We now show that any algorithm which optimizes the observables in the set $\mathcal{M}^{\epsilon}_{n}$ with error $\epsilon$ can be used to identify the parameter $v$ with high probability.

Suppose that $\mathcal{A}$ is an algorithm such that $\operatorname*{Err}(\mathcal{A},\mathcal{M}^{\epsilon}_{n})\leq\beta/9$. Then, one may use the output of $\mathcal{A}$ to construct an estimator $\hat{v}$ such that, if the objective observable is $H_{v}^{\delta}$ for $v\in\mathcal{V}$, then $\Pr[\hat{v}=v]\geq 2/3$.

By assumption, if the observable that is realized is $H^{\delta}_{v}$ for $v\in\mathcal{V}$, $\mathcal{A}$ outputs a description $\psi$ of a quantum state $\ket{\psi}$ such that

Define the estimator $\hat{v}(\psi):=\operatorname*{argmin}_{v^{\prime}\in\mathcal{V}}\expectationvalue{H_{v^{\prime}}^{\delta}}{\psi}-\lambda_{\min}(H^{\delta}_{v^{\prime}})=\operatorname*{argmin}_{v^{\prime}\in\mathcal{V}}\expectationvalue{H_{v^{\prime}}^{\delta}}{\psi}$. Lemma 5.3 implies that, if $\expectationvalue{H_{v}^{\delta}}{\psi}-\lambda_{\min}(H^{\delta}_{v})\leq\beta/3$, this estimator returns $\hat{v}=v$ with probability one. Since this event occurs with probability at least $2/3$, the estimator returns $\hat{v}=v$ with probability at least $2/3$. ∎

We have shown that the ability to optimize $\mathcal{M}^{\epsilon}_{n}$ well implies the ability to identify the hidden parameter $v\in\mathcal{V}$ with high probability. We now compute the packing parameter $\beta$ for the family $\mathcal{M}^{\epsilon}_{n}$.

For the subset $\mathcal{V}$, semimetric $d$, and packing parameter $\beta$ as defined above,

Recall that for all $v,v^{\prime}\in\{-1,1\}^{n}$,

where the minimization is over all normalized pure states on $n$ qubits. Therefore, to compute $d(v,v^{\prime})$, it suffices to compute the smallest eigenvalue of $H_{v}^{\delta}+H_{v^{\prime}}^{\delta}$.

where we used the trigonometric identities $\sqrt{2}\cos(\delta)=\cos(\pi/4+\delta)+\cos(\pi/4-\delta)=\sin(\pi/4+\delta)+\sin(\pi/4-\delta)$. From this expression, it is clear that the smallest eigenvalue of $H_{v}^{\delta}+H_{v^{\prime}}^{\delta}$ is $-2(n-\Delta(v,v^{\prime}))-2\cos(\delta)\Delta(v,v^{\prime})$, from which it follows that $d(v,v^{\prime})=2\Delta(v,v^{\prime})\quantity(1-\cos(\delta))$. By construction, for all $v\neq v^{\prime}$ with $v,v^{\prime}\in\mathcal{V}$, we have $\Delta(v,v^{\prime})\geq n/4$. It follows that $\beta\geq\frac{n}{2}(1-\cos(\delta))$.

The final inequality follows from the fact that $\cos(\delta)\leq 1-\frac{2\delta^{2}}{5}$ for $\delta\leq 0.7$. ∎

Any algorithm $\mathcal{A}$ for which $\operatorname*{Err}(\mathcal{A},\mathcal{M}^{\epsilon}_{n})\leq\epsilon$ can be used to construct an estimator $\hat{v}$ which correctly identifies the parameter $v$ of the realized observable $H^{\delta}_{v}\in\mathcal{M}^{\epsilon}_{n}$ with probability at least $2/3$.

By Lemma 5.5, the packing parameter is at least $\frac{n\delta^{2}}{5}$. Then by Lemma 5.4, if we can optimize observables in the set $\mathcal{M}^{\epsilon}_{n}$ with expected error at most $\frac{1}{9}\frac{n\delta^{2}}{5}=\frac{n\delta^{2}}{45}=\epsilon$, we can identify $v$ with probability at least $2/3$. ∎

Our proof will proceed as follows. We restrict to the subset $\mathcal{M}^{\epsilon}_{n}\subset\mathcal{H}_{n}^{\epsilon}$ and prove a lower bound on the number of zeroth-order, $100\epsilon$-vicinity queries one must make in order to identify the hidden parameter $v$ associated with the realized objective observable $H^{\delta}_{v}\in\mathcal{M}^{\epsilon}_{n}$. By Lemma 5.6, this number also lower bounds the number of such queries an algorithm $\mathcal{A}$ must make to satisfy $\operatorname*{Err}(\mathcal{A},\mathcal{M}^{\epsilon}_{n})\leq\epsilon$. Since $\mathcal{M}^{\epsilon}_{n}$ is a subset of $\mathcal{H}^{\epsilon}_{n}$, optimizing $\mathcal{M}^{\epsilon}_{n}$ is no harder than optimizing $\mathcal{H}^{\epsilon}_{n}$, and so this number also lower bounds the number of such queries needed to optimize $\mathcal{H}^{\epsilon}_{n}$ to precision $\epsilon$.

We next prove two simple lemmas we will need.

Suppose $\ket{\phi}$ is in the $\mu$-optimum of $H^{\delta}_{v}$, i.e. $\expectationvalue{H^{\delta}_{v}}{\phi}-\lambda_{\min}(H^{\delta}_{v})\leq\mu$. Let $\vec{r}_{i}$ be the polarization of the reduced state of $\ket{\phi}$ on qubit $i$, and let $\alpha_{i}\in[0,\pi]$ be the (unoriented) angle between the vector $\vec{r}_{i}$ and the unit vector $\hat{n}^{\delta v_{i}}$. Then $\frac{1}{n}\sum_{i=1}^{n}\alpha_{i}^{2}\leq\frac{10\mu}{n}$ and $\frac{1}{n}\sum_{i=1}^{n}\alpha_{i}\leq\sqrt{\frac{10\mu}{n}}$.