# Quantum SDP-Solvers: Better upper and lower bounds

Joran van Apeldoorn QuSoft, CWI, the Netherlands. Supported by the Netherlands
Organization for Scientific Research, grant number 617.001.351.
<EMAIL_ADDRESS>András Gilyén QuSoft, CWI, the Netherlands. Supported by ERC
Consolidator Grant 615307-QPROGRESS<EMAIL_ADDRESS>Sander Gribling QuSoft,
CWI, the Netherlands. Supported by the Netherlands Organization for Scientific
Research, grant number 617.001.351<EMAIL_ADDRESS>Ronald de Wolf QuSoft, CWI
and University of Amsterdam, the Netherlands. Partially supported by ERC
Consolidator Grant 615307-QPROGRESS and QuantERA project QuantAlgo 680-91-034.
<EMAIL_ADDRESS>
###### Abstract

Brandão and Svore [BS16] recently gave quantum algorithms for approximately
solving semidefinite programs, which in some regimes are faster than the best-
possible classical algorithms in terms of the dimension $n$ of the problem and
the number $m$ of constraints, but worse in terms of various other parameters.
In this paper we improve their algorithms in several ways, getting better
dependence on those other parameters. To this end we develop new techniques
for quantum algorithms, for instance a general way to efficiently implement
smooth functions of sparse Hamiltonians, and a generalized minimum-finding
procedure.

We also show limits on this approach to quantum SDP-solvers, for instance for
combinatorial optimizations problems that have a lot of symmetry. Finally, we
prove some general lower bounds showing that in the worst case, the complexity
of every quantum LP-solver (and hence also SDP-solver) has to scale linearly
with $mn$ when $m\approx n$, which is the same as classical.

###### Contents

  1. 1 Introduction
    1. 1.1 Semidefinite programs
    2. 1.2 Classical solvers for LPs and SDPs
    3. 1.3 Quantum SDP-solvers: the Brandão-Svore algorithm
    4. 1.4 Our results
      1. 1.4.1 Improved quantum SDP-solver
      2. 1.4.2 Tools that may be of more general interest
      3. 1.4.3 Lower bounds
  2. 2 An improved quantum SDP-solver
    1. 2.1 The Arora-Kale framework for solving SDPs
    2. 2.2 Approximating the expectation value $\mbox{\rm Tr}\left(A\rho\right)$ using a quantum algorithm
      1. 2.2.1 General approach
      2. 2.2.2 The special case of diagonal matrices – for LP-solving
      3. 2.2.3 General case – for SDP-solving
    3. 2.3 An efficient 2-sparse oracle
    4. 2.4 Total runtime
  3. 3 Downside of this method: general oracles are restrictive
    1. 3.1 Sparse oracles are restrictive
    2. 3.2 General width-bounds are restrictive for certain SDPs
  4. 4 Lower bounds on the quantum query complexity
  5. 5 Conclusion
  6. A Classical estimation of the expectation value $\mbox{\rm Tr}\left(A\rho\right)$
  7. B Implementing smooth functions of Hamiltonians
    1. B.1 Implementation of smooth functions of Hamiltonians: general results
    2. B.2 Applications of smooth functions of Hamiltonians
  8. C Generalized minimum-finding algorithm
  9. D Sparse matrix summation
    1. D.1 A lower bound
    2. D.2 An upper bound
  10. E Equivalence of $R$, $r$, and $\varepsilon^{-1}$
  11. F Composing the adversary bound for multiple promise functions

##  1 Introduction

###  1.1 Semidefinite programs

In the last decades, particularly since the work of Grötschel, Lovász, and
Schrijver [GLS88], _semidefinite programs_ (SDPs) have become an important
tool for designing efficient optimization and approximation algorithms. SDPs
generalize and strengthen the better-known _linear_ programs (LPs), but (like
LPs) they are still efficiently solvable. The basic form of an SDP is the
following:

$\displaystyle\max\quad$ $\displaystyle\mbox{\rm Tr}(CX)$ (1) s.t.
$\displaystyle\mbox{\rm Tr}(A_{j}X)\leq b_{j}\quad\text{ for all }j\in[m],$
$\displaystyle X\succeq 0,$

where $[m]:=\\{1,\ldots,m\\}$. The input to the problem consists of Hermitian
$n\times n$ matrices $C,A_{1},\ldots,A_{m}$ and reals $b_{1},\ldots,b_{m}$.
For normalization purposes we assume $\left\lVert C\right\rVert,\left\lVert
A_{j}\right\rVert\leq 1$. The number of constraints is $m$ (we do not count
the standard $X\succeq 0$ constraint for this). The variable $X$ of this SDP
is an $n\times n$ positive semidefinite (psd) matrix. LPs correspond to the
case where all matrices are diagonal.

A famous example is the algorithm of Goemans and Williamson [GW95] for
approximating the size of a maximum cut in a graph $G=([n],E)$: the maximum,
over all subsets $S$ of vertices, of the number of edges between $S$ and its
complement $\bar{S}$. Computing MAXCUT$(G)$ exactly is NP-hard. It corresponds
to the following integer program

$\displaystyle\max\quad$ $\displaystyle\frac{1}{2}\sum_{\\{i,j\\}\in
E}(1-v_{i}v_{j})$ s.t. $\displaystyle v_{j}\in\\{+1,-1\\}\quad\text{ for all
}j\in[n],$

using the fact that $(1-v_{i}v_{j})/2=1$ if $v_{i}$ and $v_{j}$ are different
signs, and $(1-v_{i}v_{j})/2=0$ if they are the same. We can relax this
integer program by replacing the signs $v_{j}$ by unit vectors, and replacing
the product $v_{i}v_{j}$ in the objective function by the dot product
$v_{i}^{T}v_{j}$. We can implicitly optimize over such vectors (of unspecified
dimension) by explicitly optimizing over an $n\times n$ psd matrix $X$ whose
diagonal entries are 1. This $X$ is the Gram matrix of the vectors
$v_{1},\ldots,v_{n}$, so $X_{ij}=v_{i}^{T}v_{j}$. The resulting SDP is

$\displaystyle\max\quad$ $\displaystyle\frac{1}{2}\sum_{\\{i,j\\}\in
E}(1-X_{ij})$ s.t. $\displaystyle\mbox{\rm Tr}(E_{jj}X)=1\quad\text{ for all
}j\in[n],$ $\displaystyle X\succeq 0,$

where $E_{jj}$ is the $n\times n$ matrix that has a 1 at the $(j,j)$-entry,
and 0s elsewhere.

This SDP is a relaxation of a maximization problem, so it may overshoot the
correct value, but Goemans and Williamson showed that an optimal solution to
the SDP can be rounded to a cut in $G$ whose size is within a factor $\approx
0.878$ of MAXCUT$(G)$ (which is optimal under the Unique Games Conjecture
[KKMO07]). This SDP can be massaged into the form of (1) by replacing the
equality $\mbox{\rm Tr}(E_{jj}X)=1$ by inequality $\mbox{\rm Tr}(E_{jj}X)\leq
1$ (so $m=n$) and letting $C$ be a properly normalized version of the
Laplacian of $G$.

###  1.2 Classical solvers for LPs and SDPs

Ever since Dantzig’s development of the simplex algorithm for solving LPs in
the 1940s [Dan51], much work has gone into finding faster solvers, first for
LPs and then also for SDPs. The simplex algorithm for LPs (with some
reasonable pivot rule) is usually fast in practice, but has worst-case
exponential runtime. Ellipsoid methods and interior-point methods can solve
LPs and SDPs in polynomial time; they will typically _approximate_ the optimal
value to arbitrary precision. The best known general SDP-solvers [LSW15]
approximate the optimal value OPT of such an SDP up to additive error
$\varepsilon$, with complexity

${\mathcal{O}}(m(m^{2}+n^{\omega}+mns)\,\text{polylog}(m,n,R,1/\varepsilon)),$

where $\omega\in[2,2.373)$ is the (still unknown) optimal exponent for matrix
multiplication; $s$ is the _sparsity_ : the maximal number of non-zero entries
per row of the input matrices; and $R$ is an upper bound on the trace of an
optimal $X$.111See Lee, Sidford, and Wong [LSW15, Section 10.2 of arXiv
version 2], and note that our $m,n$ are their $n,m$, their $S$ is our $mns$,
and their $M$ is our $R$. The bounds for other SDP-solvers that we state later
also include another parameter $r$; it follows from the assumptions of [LSW15,
Theorem 45 of arXiv version 2] that in their setting $r\leq mR$, and hence $r$
is absorbed in the $\log^{O(1)}(mnR/\varepsilon)$ factor. The assumption here
is that the rows and columns of the matrices of SDP (1) can be accessed as
adjacency lists: we can query, say, the $\ell$th non-zero entry of the $k$th
row of matrix $A_{j}$ in constant time.

Arora and Kale [AK16] gave an alternative way to approximate OPT, using a
matrix version of the “multiplicative weights update” method.222See also
[AHK12] for a subsequent survey; the same algorithm was independently
discovered around the same time in the context of learning theory [TRW05,
WK12]. In Section 2.1 we will describe their framework in more detail, but in
order to describe our result we will start with an overly simplified sketch
here. The algorithm goes back and forth between candidate solutions to the
primal SDP and to the corresponding _dual_ SDP, whose variables are non-
negative reals $y_{1},\ldots,y_{m}$:

$\displaystyle\min\quad$ $\displaystyle b^{T}y$ (2) s.t.
$\displaystyle\sum_{j=1}^{m}y_{j}A_{j}-C\succeq 0,$ $\displaystyle y\geq 0.$

Under assumptions that will be satisfied everywhere in this paper, strong
duality applies: the primal SDP (1) and dual SDP (2) will have the same
optimal value OPT. The algorithm does a binary search for OPT by trying
different guesses $\alpha$ for it. Suppose we have fixed some $\alpha$, and
want to find out whether $\alpha$ is bigger or smaller than OPT. Start with
some candidate solution $X^{(1)}$ for the primal, for example a multiple of
the identity matrix ($X^{(1)}$ has to be psd but need not be a _feasible_
solution to the primal). This $X^{(1)}$ induces the following polytope:

$\displaystyle\mathcal{P}_{\varepsilon}(X^{(1)}):=\\{y\in\mathbb{R}^{m}:\ $
$\displaystyle b^{T}y\leq\alpha,$ $\displaystyle\mbox{\rm
Tr}\left(\Big{(}\sum_{j=1}^{m}y_{j}A_{j}-C\Big{)}X^{(1)}\right)\geq-\varepsilon,$
$\displaystyle y\geq 0\\}.$

This polytope can be thought of as a relaxation of the feasible region of the
dual SDP with the extra constraint that $\mbox{\rm OPT}\leq\alpha$: instead of
requiring that $\sum_{j}y_{j}A_{j}-C$ is psd, we merely require that its inner
product with the particular psd matrix $X^{(1)}$ is not too negative. The
algorithm then calls an “oracle” that provides a
$y^{(1)}\in\mathcal{P}_{\varepsilon}(X^{(1)})$, or outputs “fail” if
$\mathcal{P}_{0}(X^{(1)})$ is empty (how to efficiently implement such an
oracle depends on the application). In the “fail” case we know there is no
dual-feasible $y$ with objective value $\leq\alpha$, so we can increase our
guess $\alpha$ for OPT, and restart. In case the oracle produced a $y^{(1)}$,
this is used to define a Hermitian matrix $H^{(1)}$ and a new candidate
solution $X^{(2)}$ for the primal, which is proportional to $e^{-H^{(1)}}$.
Then the oracle for the polytope $\mathcal{P}_{\varepsilon}(X^{(2)})$ induced
by this $X^{(2)}$ is called to produce a candidate
$y^{(2)}\in\mathcal{P}_{\varepsilon}(X^{(2)})$ for the dual (or “fail”), this
is used to define $H^{(2)}$ and $X^{(3)}$ proportional to $e^{-H^{(2)}}$, and
so on.

Surprisingly, the average of the dual candidates $y^{(1)},y^{(2)},\ldots$
converges to a nearly-dual-feasible solution. Let $R$ be an upper bound on the
trace of an optimal $X$ of the primal, $r$ be an upper bound on the sum of
entries of an optimal $y$ for the dual, and $w^{*}$ be the “width” of the
oracle for a certain SDP: the maximum of
$\left\lVert\sum_{j=1}^{m}y_{j}A_{j}-C\right\rVert$ over all psd matrices $X$
and all vectors $y$ that the oracle may output for the corresponding polytope
$\mathcal{P}_{\varepsilon}(X)$. In general we will not know the width of an
oracle exactly, but only an upper bound $w\geq w^{*}$, that may depend on the
SDP; this is, however, enough for the Arora-Kale framework. In Section 2.1 we
will show that without loss of generality we can assume the oracle returns a
$y$ such that $\left\lVert y\right\rVert_{1}\leq r$. Because we assumed
$\left\lVert A_{j}\right\rVert,\left\lVert C\right\rVert\leq 1$, we have
$w^{*}\leq r+1$ as an easy width-bound. General properties of the
multiplicative weights update method guarantee that after
$T=\widetilde{\mathcal{O}}(w^{2}R^{2}/\varepsilon^{2})$ iterations333The
$\widetilde{\mathcal{O}}(\cdot)$ notation hides polylogarithmic factors in all
parameters., if no oracle call yielded “fail”, then the vector
$\frac{1}{T}\sum_{t=1}^{T}y^{(t)}$ is close to dual-feasible and satisfies
$b^{T}y\leq\alpha$. This vector can then be turned into a dual-feasible
solution by tweaking its first coordinate, certifying that $\mbox{\rm
OPT}\leq\alpha+\varepsilon$, and we can decrease our guess $\alpha$ for OPT
accordingly.

The framework of Arora and Kale is really a meta-algorithm, because it does
not specify how to implement the oracle. They themselves provide oracles that
are optimized for special cases, which allows them to give a very low width-
bound for these specific SDPs. For example for the MAXCUT SDP, they obtain a
solver with near-linear runtime in the number of edges of the graph. They also
observed that the algorithm can be made more efficient by not explicitly
calculating the matrix $X^{(t)}$ in each iteration: the algorithm can still be
made to work if instead of providing the oracle with $X^{(t)}$, we feed it
good estimates of $\mbox{\rm Tr}(A_{j}X^{(t)})$ and $\mbox{\rm Tr}(CX^{(t)})$.
Arora and Kale do not describe oracles for general SDPs, but as we show at the
end of Section 2.4 (using Appendix A to estimate $\mbox{\rm Tr}(A_{j}X^{(t)})$
and $\mbox{\rm Tr}(CX^{(t)})$), one can get a general classical SDP-solver in
their framework with complexity

$\widetilde{\mathcal{O}}\left(nms\left(\frac{Rr}{\varepsilon}\right)^{4}+ns\left(\frac{Rr}{\varepsilon}\right)^{7}\right).$
(3)

Compared to the complexity of the SDP-solver of [LSW15], this has much worse
dependence on $R$ and $\varepsilon$, but better dependence on $m$ and $n$.
Using the Arora-Kale framework is thus preferable over standard SDP-solvers
for the case where $Rr$ is small compared to $mn$, and a rough approximation
to OPT (say, small constant $\varepsilon$) is good enough. It should be noted
that for many specific cases, Arora and Kale get significantly better upper
bounds than (3) by designing oracles that are specifically optimized for those
cases.

###  1.3 Quantum SDP-solvers: the Brandão-Svore algorithm

Given the speed-ups that _quantum_ computers give over classical computers for
various problems [Sho97, Gro96, DHHM06, Amb07, HHL09], it is natural to ask
whether quantum computers can solve LPs and SDPs more efficiently as well.
Very little was known about this, until recently Brandão and Svore [BS16]
discovered quantum algorithms that significantly outperform classical SDP-
solvers in certain regimes. Because of the general importance of quickly
solving LPs and SDPs, and the limited number of quantum algorithms that have
been found so far, this is a very interesting development.

The key idea of the Brandão-Svore algorithm is to take the Arora-Kale approach
and to replace two of its steps by more efficient quantum subroutines. First,
given a vector $y^{(t-1)}$, it turns out one can use “Gibbs sampling” to
prepare the new primal candidate $X^{(t)}\propto e^{-H^{(t-1)}}$ _as a
$\log(n)$-qubit quantum state $\rho^{(t)}:=X^{(t)}/\mbox{\rm Tr}(X^{(t)})$_ in
much less time than needed to compute $X^{(t)}$ as an $n\times n$ matrix.
Second, one can efficiently implement the oracle for
$\mathcal{P}_{\varepsilon}(X^{(t)})$ based on a number of copies of
$\rho^{(t)}$, using those copies to estimate $\mbox{\rm Tr}(A_{j}\rho^{(t)})$
and $\mbox{\rm Tr}(A_{j}X^{(t)})$ when needed (note that $\mbox{\rm
Tr}(A\rho)$ is the expectation value of operator $A$ for the quantum state
$\rho$). This is based on something called “Jaynes’s principle.” The resulting
oracle is weaker than what is used classically, in the sense that it outputs a
sample $j\sim y_{j}/\left\lVert y\right\rVert_{1}$ rather than the whole
vector $y$. However, such sampling still suffices to make the algorithm work
(it also means we can assume the vector $y^{(t)}$ to be quite sparse).

Using these ideas, Brandão and Svore obtain a quantum SDP-solver of complexity

$\widetilde{\mathcal{O}}(\sqrt{mn}s^{2}R^{32}/\delta^{18}),$

with _multiplicative_ error $1\pm\delta$ for the special case where $b_{j}\geq
1$ for all $j\in[m]$, and $\mbox{\rm OPT}\geq 1$ (the latter assumption allows
them to convert additive error $\varepsilon$ to multiplicative error $\delta$)
[BS16, Corollary 5 in arXiv version 4]. They describe a reduction to transform
a general SDP of the form (1) to this special case, but that reduction
significantly worsens the dependence of the complexity on the parameters $R$,
$r$, and $\delta$.

Note that compared to the runtime (3) of our general instantiation of the
original Arora-Kale framework, there are quadratic improvements in both $m$
and $n$, corresponding to the two quantum modifications made to Arora-Kale.
However, the dependence on $R,r,s$ and $1/\varepsilon$ is much worse now than
in (3). This quantum algorithm thus provides a speed-up only in regimes where
$R,r,s,1/\varepsilon$ are fairly small compared to $mn$ (finding good examples
of SDPs in such regimes is an open problem).

###  1.4 Our results

In this paper we present two sets of results: improvements to the Brandão-
Svore algorithm, and better lower bounds for the complexity of quantum LP-
solvers (and hence for quantum SDP-solvers as well).

####  1.4.1 Improved quantum SDP-solver

Our quantum SDP-solver, like the Brandão-Svore algorithm, works by quantizing
some aspects of the Arora-Kale algorithm. However, the way we quantize is
different and faster than theirs.

First, we give a more efficient procedure to estimate the quantities
$\mbox{\rm Tr}(A_{j}\rho^{(t)})$ required by the oracle. Instead of first
preparing some copies of Gibbs state $\rho^{(t)}\propto e^{-H^{(t-1)}}$ as a
mixed state, we coherently prepare a purification of $\rho^{(t)}$, which can
then be used to estimate $\mbox{\rm Tr}(A_{j}\rho^{(t)})$ more efficiently
using amplitude-estimation techniques. Also, our purified Gibbs sampler has
logarithmic dependence on the error, which is exponentially better than the
Gibbs sampler of Poulin and Wocjan [PW09b] that Brandão and Svore invoke.
Chowdhury and Somma [CS16] also gave a Gibbs sampler with logarithmic error-
dependence, but assuming query access to the entries of $\sqrt{H}$ rather than
$H$ itself.

Second, we have a different implementation of the oracle, without using Gibbs
sampling or Jaynes’s principle (though, as mentioned above, we still use
purified Gibbs sampling for approximating the $\mbox{\rm Tr}(A\rho)$
quantities). We observe that the vector $y^{(t)}$ can be made very sparse:
_two_ non-zero entries suffice.444Independently of us, Ben David, Eldar, Garg,
Kothari, Natarajan, and Wright (at MIT), and separately Ambainis observed that
in the special case where all $b_{j}$ are at least 1, the oracle can even be
made 1-sparse, and the one entry can be found using one Grover search over $m$
points (in both cases personal communication 2017). The same happens
implicitly in our Section 2.3 in this case. However, in general 2 non-zero
entries are necessary in $y$. We then show how we can efficiently find such a
2-sparse vector (rather than merely sampling from it) using two applications
of a new generalization of the well-known quantum minimum-finding algorithm of
Dürr and Høyer [DH96], which is based on Grover search [Gro96].

These modifications both simplify and speed up the quantum SDP-solver,
resulting in complexity

$\widetilde{\mathcal{O}}(\sqrt{mn}s^{2}(Rr/\varepsilon)^{8}).$

The dependence on $m$, $n$, and $s$ is the same as in Brandão-Svore, but our
dependence on $R$, $r$, and $1/\varepsilon$ is substantially better. Note that
each of the three parameters $R$, $r$, and $1/\varepsilon$ now occurs with the
same 8th power in the complexity. This is no coincidence: as we show in
Appendix E, these three parameters can all be traded for one another, in the
sense that we can massage the SDP to make each one of them small at the
expense of making the others proportionally bigger. These trade-offs suggest
we should actually think of $Rr/\varepsilon$ as _one_ parameter of the primal-
dual pair of SDPs, not three separate parameters. For the special case of LPs,
we can improve the runtime to

$\widetilde{\mathcal{O}}(\sqrt{mn}(Rr/\varepsilon)^{5}).$

Like in Brandão-Svore, our quantum oracle produces very sparse vectors $y$, in
our case even of sparsity 2. This means that after $T$ iterations, the final
$\varepsilon$-optimal dual-feasible vector (which is a slightly tweaked
version of the average of the $T$ $y$-vectors produced in the $T$ iterations)
has only $\mathcal{O}(T)$ non-zero entries. Such sparse vectors have some
advantages, for example they take much less space to store than arbitrary
$y\in\mathbb{R}^{m}$. In fact, to get a sublinear running time in terms of
$m$, this is necessary. However, this sparsity of the algorithm’s output also
points to a weakness of these methods: if _every_ $\varepsilon$-optimal dual-
feasible vector $y$ has many non-zero entries, then the number of iterations
needs to be large. For example, if every $\varepsilon$-optimal dual-feasible
vector $y$ has $\Omega(m)$ non-zero entries, then these methods require
$T=\Omega(m)$ iterations before they can reach an $\varepsilon$-optimal dual-
feasible vector. Since
$T={\mathcal{O}}\left(\frac{R^{2}r^{2}}{\varepsilon^{2}}\ln(n)\right)$ this
would imply that $\frac{Rr}{\varepsilon}=\Omega(\sqrt{m/\ln(n)})$, and hence
many classical SDP-solvers would have a better complexity than our quantum
SDP-solver. As we show in Section 3, this will actually be the case for
families of SDPs that have a lot of symmetry.

####  1.4.2 Tools that may be of more general interest

Along the way to our improved SDP-solver, we developed some new techniques
that may be of independent interest. These are mostly tucked away in
appendices, but here we will highlight two.

##### Implementing smooth functions of a given Hamiltonian.

In Appendix B we describe a general technique to apply a function $f(H)$ of a
sparse Hamiltonian $H$ to a given state $|\phi\rangle$. Roughly speaking, what
this means is that we want a unitary circuit that maps $|0\rangle|\phi\rangle$
to $|0\rangle f(H)|\phi\rangle+|1\rangle|*\rangle$. If need be, we can then
combine this with amplitude amplification to boost the $|0\rangle
f(H)|\phi\rangle$ part of the state. If the function
$f:\mathbb{R}\to\mathbb{C}$ can be approximated well by a low-degree Fourier
series, then our preparation will be efficient in the sense of using few
queries to $H$ and few other gates. The novelty of our approach is that we
construct a good Fourier series from the polynomial that approximates $f$ (for
example a truncated Taylor series for $f$). Our Theorem 40 can be easily
applied to various smooth functions without using involved integral
approximations, unlike previous works building on similar techniques. Our most
general result Corollary 42 only requires that the function $f$ can be nicely
approximated locally around each possible eigenvalues of $H$, improving on
Theorem 40.

In this paper we mostly care about the function $f(x)=e^{-x}$, which is what
we want for generating a purification of the Gibbs state corresponding to $H$;
and the function $f(x)=\sqrt{x}$, which is what we use for estimating
quantities like $\mbox{\rm Tr}(A\rho)$. However, our techniques apply much
more generally than these two functions. For example, they also simplify the
analysis of the improved linear-systems solver of Childs et al. [CKS15], where
the relevant function is $f(x)=1/x$. As in their work, the Linear Combination
of Unitaries technique of Childs et al. [CW12, BCC+15, BCK15] is a crucial
tool for us.

##### A generalized minimum-finding algorithm.

Dürr and Høyer [DH96] showed how to find the minimal value of a function
$f:[N]\to\mathbb{R}$ using ${\mathcal{O}}(\sqrt{N})$ queries to $f$, by
repeatedly using Grover search to find smaller and smaller elements of the
range of $f$. In Appendix C we describe a more general minimum-finding
procedure. Suppose we have a unitary $U$ which prepares a quantum state
$U|0\rangle=\sum_{k=1}^{N}|\psi_{k}\rangle|x_{k}\rangle$, where the
$|\psi_{k}\rangle$ are unnormalized states. Our procedure can find the minimum
value $x_{k^{*}}$ among the $x_{k}$’s that have support in the second
register, using roughly ${\mathcal{O}}(1/\left\lVert\psi_{k^{*}}\right\rVert)$
applications of $U$ and $U^{-1}$. Also, upon finding the minimal value $k^{*}$
the procedure actually outputs the normalized state proportional to
$|\psi_{k^{*}}\rangle|x_{k^{*}}\rangle$. This immediately gives the Dürr-Høyer
result as a special case, if we take $U$ to produce
$U|0\rangle=\frac{1}{\sqrt{N}}\sum_{k=1}^{N}|k\rangle|f(k)\rangle$ using one
query to $f$. Unlike Dürr-Høyer, we need not assume direct query access to the
individual values $f(k)$.

More interestingly for us, for a given $n$-dimensional Hamiltonian $H$, if we
combine our minimum-finder with phase estimation using unitary $U=e^{iH}$ on
one half of a maximally entangled state, then we obtain an algorithm for
estimating the smallest eigenvalue of $H$ (and preparing its ground state)
using roughly ${\mathcal{O}}(\sqrt{n})$ applications of phase estimation with
$U$. A similar result on approximating the smallest eigenvalue of a
Hamiltonian was already shown by Poulin and Wocjan [PW09a], but we improve on
their analysis to be able to apply it as a subroutine in our procedure to
estimate $\mbox{\rm Tr}(A_{j}\rho)$.

####  1.4.3 Lower bounds

What about lower bounds for quantum SDP-solvers? Brandão and Svore already
proved that a quantum SDP-solver has to make $\Omega(\sqrt{n}+\sqrt{m})$
queries to the input matrices, for some SDPs. Their lower bound is for a
family of SDPs where $s,R,r,1/\varepsilon$ are all constant, and is by
reduction from a search problem.

In this paper we prove lower bounds that are quantitatively stronger in $m$
and $n$, but for SDPs with non-constant $R$ and $r$. The key idea is to
consider a Boolean function $F$ on $N=abc$ input bits that is the composition
of an $a$-bit majority function with a $b$-bit OR function with a $c$-bit
majority function. The known quantum query complexities of majority and OR,
combined with composition properties of the adversary lower bound, imply that
every quantum algorithm that computes this function requires
$\Omega(a\sqrt{b}c)$ queries. We define a family of LPs, with constant
$1/\varepsilon$ but non-constant $r$ and $R$ (we could massage this to make
$R$ or $r$ constant using the results of Appendix E, but not
$Rr/\varepsilon$), such that constant-error approximation of OPT computes $F$.
Choosing $a$, $b$, and $c$ appropriately, this implies a lower bound of

$\Omega\left(\sqrt{\max\\{n,m\\}}\left(\min\\{n,m\\}\right)^{3/2}\right)$

queries to the entries of the input matrices for quantum LP-solvers. Since LPs
are SDPs with sparsity $s=1$, we get the same lower bound for quantum SDP-
solvers. If $m$ and $n$ are of the same order, this lower bound is
$\Omega(mn)$, the same scaling with $mn$ as the classical general
instantiation of Arora-Kale (3). In particular, this shows that we cannot have
an $O(\sqrt{mn})$ upper bound without simultaneously having polynomial
dependence on $Rr/\varepsilon$. The construction of our lower bound implies
that for the case $m\approx n$, this polynomial dependence has to be at least
$(Rr/\varepsilon)^{1/4}$.

##### Subsequent work.

Following the first version of our paper, improvements in the running time
were obtained in [BKL+17, AG18], the latter providing a runtime of
$\widetilde{\mathcal{O}}\left((\sqrt{m}+\sqrt{n}\frac{Rr}{\epsilon})s\left(\frac{Rr}{\epsilon}\right)^{4}\right)$.
More recently, a quantum interior point method for solving SDPs and LPs was
obtained by Kerenidis and Prakash [KP18]. It is hard to compare the latter
algorithm to the other SDP-solvers for two reasons. First, the output of their
algorithm consists only of almost-feasible solutions to the primal and dual
(their algorithm has a polynomial dependence on the distance to feasibility).
It is therefore not clear what their output means for the optimal value of the
SDPs. Secondly, the runtime of their algorithm depends polynomially on the
condition number of the matrices that the interior point method encounters,
and no explicit bounds for these condition numbers are given.

Our results on implementing smooth functions of a given Hamiltonian have been
extended to more general input models (block-encodings) in [GSLW18]. This
recent paper builds on some of our techniques, but achieves slightly improved
complexities by directly implementing the transformations without using
Hamiltonian simulation as a subroutine.

Very recently van Apeldoorn et al. [AGGW18] and Chakrabarti et al. [CCLW18]
developed quantum algorithms for general black-box convex optimization, where
one optimizes over a general convex set $K$, and the access to $K$ is via
membership and/or separation oracles. Since we work in a model where we are
given access directly to the constraints defining the problem, our results are
incomparable to theirs.

##### Organization.

The paper is organized as follows. In Section 2 we start with a description of
the Arora-Kale framework for SDP-solvers, and then we describe how to quantize
different aspects of it to obtain a quantum SDP-solver with better dependence
on $R$, $r$, and $1/\varepsilon$ (or rather, on $Rr/\varepsilon$) than Brandão
and Svore got. In Section 3 we describe the limitations of primal-dual SDP-
solvers using general oracles (not optimized for specific SDPs) that produce
sparse dual solutions $y$: if good solutions are dense, this puts a lower
bound on the number of iterations needed. In Section 4 we give our lower
bounds. A number of the proofs are relegated to the appendices: how to
classically approximate $\mbox{\rm Tr}(A_{j}\rho)$ (Appendix A), how to
efficiently implement smooth functions of Hamiltonians on a quantum computer
(Appendix B), our generalized method for minimum-finding (Appendix C), upper
and lower bounds on how efficiently we can query entries of sums of sparse
matrices (Appendix D), how to trade off $R$, $r$, and $1/\varepsilon$ against
each other (Appendix E), and the composition property of the adversary method
that we need for our lower bounds (Appendix F).

##  2 An improved quantum SDP-solver

Here we describe our quantum SDP-solver. In Section 2.1 we describe the
framework designed by Arora and Kale for solving semidefinite programs. As in
the recent work by Brandão and Svore, we use this framework to design an
efficient quantum algorithm for solving SDPs. In particular, we show that the
key subroutine needed in the Arora-Kale framework can be implemented
efficiently on a quantum computer. Our implementation uses different
techniques than the quantum algorithm of Brandão and Svore, allowing us to
obtain a faster algorithm. The techniques required for this subroutine are
developed in Sections 2.2 and 2.3. In Section 2.4 we put everything together
to prove the main theorem of this section (the notation is explained below):

######  Theorem 1.

Instantiating Meta-Algorithm 1 using the trace calculation algorithm from
Section 2.2 and the oracle from Section 2.3 (with width-bound $w:=r+1$), and
using this to do a binary search for $\mbox{\rm OPT}\in[-R,R]$ (using
different guesses $\alpha$ for OPT), gives a quantum algorithm for solving
SDPs of the form (1), which (with high probability) produces a feasible
solution $y$ to the dual program which is optimal up to an additive error
$\varepsilon$, and uses

$\widetilde{\mathcal{O}}\left(\sqrt{nm}s^{2}\left(\frac{Rr}{\varepsilon}\right)^{\\!\\!8}\right)$

queries to the input matrices and the same order of other gates.

##### Notation/Assumptions.

We use $\log$ to denote the logarithm in base $2$. We denote the all-zero
matrix and vector by $0$. Throughout we assume each element of the input
matrices can be represented by a bitstring of size $\text{poly}(\log n,\log
m)$. We use $s$ as the sparsity of the input matrices, that is, the maximum
number of non-zero entries in a row (or column) of any of the matrices
$C,A_{1},\ldots,A_{m}$ is $s$. Recall that for normalization purposes we
assume $\left\lVert A_{1}\right\rVert,\ldots,\left\lVert
A_{m}\right\rVert,\left\lVert C\right\rVert\leq 1$. We furthermore assume that
$A_{1}=I$ and $b_{1}=R$, that is, the trace of primal-feasible solutions is
bounded by $R$ (and hence also the trace of primal-optimal solutions is
bounded by $R$). The analogous quantity for the dual SDP (2), an upper bound
on $\sum_{j=1}^{m}y_{j}$ for an optimal dual solution $y$, will be denoted by
$r$. However, we do not add the constraint $\sum_{j=1}^{m}y_{j}\leq r$ to the
dual. We will assume $r\geq 1$. For $r$ to be well-defined we have to make the
explicit assumption that the optimal solution in the dual is attained. In
Section 3 it will be necessary to work with the best possible upper bounds: we
let $R^{*}$ be the smallest trace of an optimal solution to SDP (1), and we
let $r^{*}$ be the smallest $\ell_{1}$-norm of an optimal solution to the
dual. These quantities are well-defined; indeed, both the primal and dual
optimum are attained: the dual optimum is attained by assumption, and due to
the assumption $A_{1}=I$, the dual SDP is strictly feasible, which means that
the optimum in (1) is attained.

Unless specified otherwise, we always consider _additive_ error. In
particular, an $\varepsilon$-optimal solution to an SDP will be a feasible
solution whose objective value is within additive error $\varepsilon$ of the
optimum.

##### Input oracles:

We assume sparse black-box access to the elements of the matrices
$C,A_{1},\ldots,A_{m}$ defined in the following way: for input
$(j,k,\ell)\in(\\{0\\}\cup[m])\times[n]\times[s]$ we can query the location
and value of the $\ell$th non-zero entry in the $k$th row of the matrix
$A_{j}$ (where $j=0$ would indicate the $C$ matrix).

Specifically in the quantum case, as described in [BCK15], we assume access to
an oracle $O_{I}$ which serves the purpose of sparse access. $O_{I}$
calculates the $\text{index}_{A_{j}}:[n]\times[s]\to[n]$ function, which for
input $(k,\ell)$ gives the column index of the $\ell$th non-zero element in
the $k$th row of $A_{j}$. We assume this oracle computes the index “in place”:

$O_{I}|j,k,\ell\rangle=|j,k,\text{index}_{A_{j}}(k,\ell)\rangle.$ (4)

(In the degenerate case where the $k$th row has fewer than $\ell$ non-zero
entries, $\text{index}_{A_{j}}(k,\ell)$ is defined to be $\ell$ together with
some special symbol.) We also assume we can apply the inverse of $O_{I}$.

We also need another oracle $O_{M}$, returning a bitstring representation of
$(A_{j})_{ki}$ for any $j\in\\{0\\}\cup[m]$ and $k,i\in[n]$:

$O_{M}|j,k,i,z\rangle=|j,k,i,z\oplus{(A_{j})_{ki}}\rangle.$ (5)

The slightly unusual “in place” definition of oracle $O_{I}$ is not too
demanding. In particular, if instead we had an oracle that computed the non-
zero entries of a row in order, then we could implement both $O_{I}$ and its
inverse using $\log(s)$ queries (we can compute $\ell$ from $k$ and
$\text{index}_{A_{j}}(k,\ell)$ using binary search) [BCK15].

##### Computational model:

As our computational model, we assume a slight relaxation of the usual quantum
circuit model: a classical control system that can run quantum subroutines. We
limit the classical control system so that its number of operations is at most
a polylogarithmic factor bigger than the gate complexity of the quantum
subroutines, i.e., if the quantum subroutines use $C$ gates, then the
classical control system may not use more than
${\mathcal{O}}(C\,\text{polylog}(C))$ operations.

When we talk about gate complexity, we count the number of two-qubit quantum
gates needed for implementation of the quantum subroutines. Additionally, we
assume for simplicity that there exists a unit-cost QRAM gate that allows us
to store and retrieve qubits in a memory, by means of a swap of two registers
indexed by another register:

$QRAM:|i,x,r_{1},\dots,r_{K}\rangle\mapsto|i,r_{i},r_{1},\dots,r_{i-1},x,r_{i+1},\ldots,r_{K}\rangle,$

where the registers $r_{1},\dots,r_{K}$ are only accessible through this gate.
The QRAM gate can be seen as a quantum analogue of pointers in classical
computing. The only place where we need QRAM is in Appendix D, for a data
structure that allows efficient access to the non-zero entries of a sum of
sparse matrices; for the special case of LP-solving it is not needed.

###  2.1 The Arora-Kale framework for solving SDPs

In this section we give a short introduction to the Arora-Kale framework for
solving semidefinite programs. We refer to [AK16, AHK12] for a more detailed
description and omitted proofs.

The key building block is the Matrix Multiplicative Weights (MMW) algorithm
introduced by Arora and Kale in [AK16]. The MMW algorithm can be seen as a
strategy for you in a game between you and an adversary. We first introduce
the game. There is a number of rounds $T$. In each round you present a density
matrix $\rho$ to an adversary, the adversary replies with a loss matrix $M$
satisfying $-I\preceq M\preceq I$. After each round you have to pay $\mbox{\rm
Tr}\left(M\rho\right)$. Your objective is to pay as little as possible. The
MMW algorithm is a strategy for you that allows you to lose not too much, in a
sense that is made precise below. In Algorithm 1 we state the MMW algorithm,
the following theorem shows the key property of the output of the algorithm.

Input

Parameter $\eta\leq 1$, number of rounds $T$.

Rules

In each round player $1$ (you) presents a density matrix $\rho$, player $2$
(the adversary) replies with a matrix $M$ satisfying $-I\preceq M\preceq I$.

Output

A sequence of symmetric $n\times n$ matrices $M^{(1)},\ldots,M^{(T)}$
satisfying $-I\preceq M^{(t)}\preceq I$, for $t\in[T]$ and a sequence of
$n\times n$ psd matrices $\rho^{(1)},\ldots,\rho^{(T)}$ satisfying $\mbox{\rm
Tr}\left(\rho^{(t)}\right)=1$ for $t\in[T]$.

Strategy of player $1$:

Take $\rho^{(1)}:=I/n$

In round $t$:

  1. 1.

Show the density matrix $\rho^{(t)}$ to the adversary.

  2. 2.

Obtain the loss matrix $M^{(t)}$ from the adversary.

  3. 3.

Update the density matrix as follows:

$\rho^{(t+1)}:=\left.\exp\left(-\eta\sum_{\tau=1}^{t}M^{(\tau)}\right)\right/\mbox{\rm
Tr}\left(\exp\left(-\eta\sum_{\tau=1}^{t}M^{(\tau)}\right)\right)$

Algorithm 1 Matrix Multiplicative Weights (MMW) Algorithm

######  Theorem 2 ([AK16, Theorem 3.1]).

For every adversary, the sequence $\rho^{(1)},\ldots,\rho^{(T)}$ of density
matrices constructed using the Matrix Multiplicative Weights Algorithm (1)
satisfies

$\sum_{t=1}^{T}\mbox{\rm
Tr}\left(M^{(t)}\rho^{(t)}\right)\leq\lambda_{\min}\left(\sum_{t=1}^{T}M^{(t)}\right)+\eta\sum_{t=1}^{T}\mbox{\rm
Tr}\left((M^{(t)})^{2}\rho^{(t)}\right)+\frac{\ln(n)}{\eta}.$

Arora and Kale use the MMW algorithm to construct an SDP-solver. For that,
they construct an adversary who promises to satisfy an additional condition:
in each round $t$, the adversary returns a matrix $M^{(t)}$ whose trace inner
product with the density matrix $\rho^{(t)}$ is non-negative. The above
theorem shows that then, after $T$ rounds, the average of the adversary’s
responses satisfies the stronger condition that its smallest eigenvalue is not
too negative:
$\lambda_{\min}\left(\frac{1}{T}\sum_{t=1}^{T}M^{(t)}\right)\geq-\eta-\frac{\ln(n)}{\eta
T}$. More explicitly, the MMW algorithm is used to build a vector $y\geq 0$
such that

$\frac{1}{T}\sum_{t=1}^{T}M^{(t)}\propto\sum_{j=1}^{m}y_{j}A_{j}-C$

and $b^{T}y\leq\alpha$. That is, the smallest eigenvalue of the matrix
$\sum_{j=1}^{m}y_{j}A_{j}-C$ is only slightly below zero and $y$’s objective
value is at most $\alpha$. Since $A_{1}=I$, increasing the first coordinate of
$y$ makes the smallest eigenvalue of $\sum_{j}y_{j}A_{j}-C$ bigger, so that
this matrix becomes psd and hence dual-feasible. By the above we know how much
the minimum eigenvalue has to be shifted, and with the right choice of
parameters it can be shown that this gives a dual-feasible vector
$\overline{y}$ that satisfies $b^{T}\overline{y}\leq\alpha+\varepsilon$. In
order to present the algorithm formally, we require some definitions.

Given a candidate solution $X\succeq 0$ for the primal problem (1) and a
parameter $\varepsilon\geq 0$, define the polytope

$\displaystyle\mathcal{P}_{\varepsilon}(X):=\\{y\in\mathbb{R}^{m}:\ $
$\displaystyle b^{T}y\leq\alpha,$ $\displaystyle\mbox{\rm
Tr}\left(\Big{(}\sum_{j=1}^{m}y_{j}A_{j}-C\Big{)}X\right)\geq-\varepsilon,$
$\displaystyle y\geq 0\\}.$

One can verify the following:

######  Lemma 3 ([AK16, Lemma 4.2]).

If for a given candidate solution $X\succeq 0$ the polytope
$\mathcal{P}_{0}(X)$ is empty, then a scaled version of $X$ is primal-feasible
and of objective value at least $\alpha$.

The Arora-Kale framework for solving SDPs uses the MMW algorithm where the
role of the adversary is taken by an $\varepsilon$-approximate oracle:

Input

An $n\times n$ psd matrix $X$, a parameter $\varepsilon$, and the input
matrices and reals of (2).

Output

Either the Oracleε returns a vector $y$ from the polytope
$\mathcal{P}_{\varepsilon}(X)$ or it outputs “fail”. It may only output fail
if $\mathcal{P}_{0}(X)=\emptyset$.

Algorithm 2 $\varepsilon$-approximate Oracleε for maximization SDPs

As we will see later, the runtime of the Arora-Kale framework depends on a
property of the oracle called the _width_ :

######  Definition 4 (_Width_ of Oracleε).

The _width_ of Oracleε for an SDP is the smallest $w^{*}\geq 0$ such that for
every primal candidate $X\succeq 0$, the vector $y$ returned by Oracleε
satisfies $\left\lVert\sum_{j=1}^{m}y_{j}A_{j}-C\right\rVert\leq w^{*}$.

In practice, the width of an oracle is not always known. However, it suffices
to work with an upper bound $w\geq w^{*}$: as we can see in Meta-Algorithm 1,
the purpose of the width is to rescale the matrix $M^{(t)}$ in such a way that
it forms a valid response for the adversary in the MMW algorithm.

Input

The input matrices and reals of SDP (1) and trace bound $R$. The current guess
$\alpha$ of the optimal value of the dual (2). An additive error tolerance
$\varepsilon>0$. An $\frac{\varepsilon}{3}$-approximate oracle Oracleε/3 as in
Algorithm 2 with width-bound $w$.

Output

Either “Lower” and a vector $\overline{y}\in\mathbb{R}^{m}_{+}$ feasible for
(2) with $b^{T}\overline{y}\leq\alpha+\varepsilon$
or “Higher” and a symmetric $n\times n$ matrix $X$ that, when scaled suitably,
is primal-feasible with objective value at least $\alpha$.

$T:=\left\lceil\frac{9w^{2}R^{2}\ln(n)}{\varepsilon^{2}}\right\rceil$.

$\eta:=\sqrt{\frac{\ln(n)}{T}}$.

$\rho^{(1)}:=I/n$

for $t=1,\dots,T$ do

Run Oracleε/3 with $X^{(t)}=R\rho^{(t)}$.

if Oracleε/3 outputs “fail” then

return “Higher” and a description of $X^{(t)}$.

end if

Let $y^{(t)}$ be the vector generated by Oracleε/3.

Set $M^{(t)}=\frac{1}{w}\left(\sum_{j=1}^{m}y_{j}^{(t)}A_{j}-C\right)$.

Define $H^{(t)}=\sum_{\tau=1}^{t}M^{(\tau)}$.

Update the state matrix as follows: $\rho^{(t+1)}:=\exp\left(-\eta
H^{(t)}\right)/\mbox{\rm Tr}\left(\exp\left(-\eta H^{(t)}\right)\right)$.

end for

If Oracleε/3 does not output “fail” in any of the $T$ rounds, then output the
dual solution
$\overline{y}=\frac{\varepsilon}{R}e_{1}+\frac{1}{T}\sum_{t=1}^{T}y^{(t)}$
where $e_{1}=(1,0,\ldots,0)\in\mathbb{R}^{m}$.

Meta-Algorithm 1 Primal-Dual Algorithm for solving SDPs

The following theorem shows the correctness of the Arora-Kale primal-dual
meta-algorithm for solving SDPs, stated in Meta-Algorithm 1:

######  Theorem 5 ([AK16, Theorem 4.7]).

Given an SDP of the form (1) with input matrices $A_{1}=I,A_{2},\ldots,A_{m}$
and $C$ having operator norm at most $1$, and input reals
$b_{1}=R,b_{2},\ldots,b_{m}$. Assume Meta-Algorithm 1 does not output “fail”
in any of the rounds, then the returned vector $\overline{y}$ is feasible for
the dual (2) with objective value at most $\alpha+\varepsilon$. If Oracleε/3
outputs “fail” in the $t$-th round then a suitably scaled version of $X^{(t)}$
is primal-feasible with objective value at least $\alpha$.

The SDP-solver uses
$T=\left\lceil\frac{9w^{2}R^{2}\ln(n)}{\varepsilon^{2}}\right\rceil$
iterations. In each iteration several steps have to be taken. The most
expensive two steps are computing the matrix exponential of the matrix $-\eta
H^{(t)}$ and the application of the oracle. Note that the only purpose of
computing the matrix exponential is to allow the oracle to compute the values
$\mbox{\rm Tr}\left(A_{j}X\right)$ for all $j$ and $\mbox{\rm
Tr}\left(CX\right)$, since the polytope depends on $X$ only through those
values. To obtain faster algorithms it is important to note, as was done
already by Arora and Kale, that the primal-dual algorithm also works if we
provide a (more accurate) oracle with approximations of $\mbox{\rm
Tr}\left(A_{j}X\right)$. Let $a_{j}:=\mbox{\rm
Tr}\left(A_{j}\rho\right)=\mbox{\rm Tr}\left(A_{j}X\right)/\mbox{\rm
Tr}\left(X\right)$ and $c:=\mbox{\rm Tr}\left(C\rho\right)=\mbox{\rm
Tr}\left(CX\right)/\mbox{\rm Tr}\left(X\right)$. Then, given a list of reals
$\tilde{a}_{1},\ldots,\tilde{a}_{m},\tilde{c}$ and a parameter $\theta\geq 0$,
such that $|\tilde{a}_{j}-a_{j}|\leq\theta$ for all $j$, and
$|\tilde{c}-c|\leq\theta$, we define the polytope

$\displaystyle\tilde{\mathcal{P}}(\tilde{a}_{1},\ldots,\tilde{a}_{m},\tilde{c}-(r+1)\theta):=\\{y\in\mathbb{R}^{m}:\
$ $\displaystyle b^{T}y\leq\alpha,$ $\displaystyle\sum_{j=1}^{m}y_{j}\leq r,$
$\displaystyle\sum_{j=1}^{m}\tilde{a}_{j}y_{j}\geq\tilde{c}-(r+1)\theta$
$\displaystyle y\geq 0\\}.$

For convenience we will denote
$\tilde{a}=(\tilde{a}_{1},\ldots,\tilde{a}_{m})$ and
$c^{\prime}:=\tilde{c}-(r+1)\theta$. Notice that $\tilde{\mathcal{P}}$ also
contains a new type of constraint: $\sum_{j}y_{j}\leq r$. Recall that $r$ is
defined as a positive real such that there exists an optimal solution $y$ to
SDP (2) with $\left\lVert y\right\rVert_{1}\leq r$. Hence, using that
$\mathcal{P}_{0}(X)$ is a _relaxation_ of the feasible region of the dual
(with bound $\alpha$ on the objective value), we may restrict our oracle to
return only such $y$:

$\mathcal{P}_{0}(X)\neq\emptyset\Rightarrow\mathcal{P}_{0}(X)\cap\\{y\in\mathbb{R}^{m}:\sum_{j=1}^{m}y_{j}\leq
r\\}\neq\emptyset.$

The benefit of this restriction is that an oracle that always returns a vector
with bounded $\ell_{1}$-norm automatically has a width $w^{*}\leq r+1$, due to
the assumptions on the norms of the input matrices. The downside of this
restriction is that the analogue of Lemma 3 does not hold for
$\mathcal{P}_{0}(X)\cap\\{y\in\mathbb{R}^{m}:\sum_{j}y_{j}\leq r\\}$.555Using
several transformations of the SDP, from Appendix E and Lemma 2 of [BS16], one
can show that there is a way to remove the need for this restriction. Hence,
after these modifications, if for a given candidate solution $X\succeq 0$ the
oracle outputs that the set $\mathcal{P}_{0}(X)$ is empty, then a scaled
version of $X$ is primal feasible for this new SDP, with objective value at
least $\alpha$. This scaled version of $X$ can be modified to a near-feasible
solution to the original SDP (it will be psd, but it might violate the linear
constraints a little bit) with nearly the same objective value.

The following shows that an oracle that always returns a vector
$y\in\tilde{\mathcal{P}}(\tilde{a},c^{\prime})$ if one exists, is a
$4Rr\theta$-approximate oracle as defined in Algorithm 2.

######  Lemma 6.

Let $\tilde{a}_{1},\ldots,\tilde{a}_{m}$ and $\tilde{c}$ be
$\theta$-approximations of $\mbox{\rm
Tr}\left(A_{1}\rho\right),\ldots,\mbox{\rm Tr}\left(A_{m}\rho\right)$ and
$\mbox{\rm Tr}\left(C\rho\right)$, respectively, where $X=R\rho$. Then the
following holds:

$\mathcal{P}_{0}(X)\cap\\{y\in\mathbb{R}^{m}:\sum_{j=1}^{m}y_{j}\leq
r\\}\subseteq\tilde{\mathcal{P}}(\tilde{a},c^{\prime})\subseteq\mathcal{P}_{4Rr\theta}(X).$

###### Proof.

First, suppose
$y\in\mathcal{P}_{0}(X)\cap\\{y\in\mathbb{R}^{m}:\sum_{j}y_{j}\leq r\\}$. We
now have

$\sum_{j=1}^{m}\tilde{a}_{j}y_{j}-c\geq-\sum_{j=1}^{m}|\tilde{a}_{j}-\mbox{\rm
Tr}\left(A_{j}\rho\right)|y_{j}-|\tilde{c}-\mbox{\rm
Tr}\left(C\rho\right)|\geq-\theta\left\lVert
y\right\rVert_{1}-\theta\geq-(r+1)\theta$

which shows that $y\in\tilde{\mathcal{P}}(\tilde{a},c^{\prime})$.

Next, suppose $y\in\tilde{\mathcal{P}}(\tilde{a},c^{\prime})$. We show that
$y\in\mathcal{P}_{4Rr\theta}(X)$. Indeed, since $|\mbox{\rm
Tr}\left(A_{j}\rho\right)-\tilde{a}_{j}|\leq\theta$ we have

$\mbox{\rm
Tr}\left(\\!\\!\left(\sum_{j=1}^{m}y_{j}A_{j}-C\right)\rho\\!\right)\geq\left(\sum_{j=1}^{m}\tilde{a}_{j}y_{j}+\tilde{c}\right)-(r+1)\theta\geq-(2+r+\left\lVert
y\right\rVert_{1})\theta\geq\\!-4r\theta$

where the last inequality used our assumptions $r\geq 1$ and $\left\lVert
y\right\rVert_{1}\leq r$. Hence

$\mbox{\rm
Tr}\left(\left(\sum_{j=1}^{m}y_{j}A_{j}-C\right)X\right)\geq-4r\mbox{\rm
Tr}\left(X\right)\theta=-4Rr\theta.$

For the latter inequality we use $\mbox{\rm Tr}\left(X\right)=R$. ∎

We have now seen the Arora-Kale framework for solving SDPs. To obtain a
quantum SDP-solver it remains to provide a quantum oracle subroutine. By the
above discussion it suffices to set $\theta=\varepsilon/(12Rr)$ and to use an
oracle that is based on $\theta$-approximations of $\mbox{\rm
Tr}\left(A\rho\right)$ (for $A\in\\{A_{1},\ldots,A_{m},C\\}$), since with that
choice of $\theta$ we have
$\mathcal{P}_{4Rr\theta}(X)=\mathcal{P}_{\varepsilon/3}(X)$. In the section
below we first give a quantum algorithm for approximating $\mbox{\rm
Tr}\left(A\rho\right)$ efficiently (see also Appendix A for a classical
algorithm). Then, in Section 2.3, we provide an oracle using those estimates.
The oracle will be based on a simple geometric idea and can be implemented
both on a quantum computer and on a classical computer (of course, resulting
in different runtimes). In Section 2.4 we conclude with an overview of the
runtime of our quantum SDP-solver. We want to stress that our solver is meant
to work for any SDP. In particular, our oracle does not use the structure of a
specific SDP. As we will show in Section 3, any oracle that works for all SDPs
necessarily has a large width-bound. To obtain quantum speedups for a
_specific_ class of SDPs it will be necessary to develop oracles tuned to that
problem, we view this as an important direction for future work. Recall from
the introduction that Arora and Kale also obtain fast classical algorithms for
problems such as MAXCUT by developing specialized oracles.

###  2.2 Approximating the expectation value $\mbox{\rm Tr}\left(A\rho\right)$
using a quantum algorithm

In this section we give an efficient quantum algorithm to approximate
quantities of the form $\mbox{\rm Tr}\left(A\rho\right)$. We are going to work
with Hermitian matrices $A,H\in\mathbb{C}^{n\times n}$, such that $\rho$ is
the Gibbs state $e^{-H}/\mbox{\rm Tr}\left(e^{-H}\right)$. Note the analogy
with quantum physics: in physics terminology $\mbox{\rm Tr}\left(A\rho\right)$
is simply called the “expectation value” of $A$ for a quantum system in a
thermal state corresponding to $H$.

The general approach is to separately estimate $\mbox{\rm
Tr}\left(Ae^{-H}\right)$ and $\mbox{\rm Tr}\left(e^{-H}\right)$, and then to
use the ratio of these estimates as an approximation of $\mbox{\rm
Tr}\left(A\rho\right)=\mbox{\rm Tr}\left(Ae^{-H}\right)/\mbox{\rm
Tr}\left(e^{-H}\right)$. Both estimations are done using state preparation to
prepare a pure state with a flag, such that the probability that the flag is 0
is proportional to the quantity we want to estimate, and then to use amplitude
estimation to estimate that probability. Below in Section 2.2.1 we first
describe the general approach. In Section 2.2.2 we then instantiate this for
the special case where all matrices are diagonal, which is the relevant case
for LP-solving. In Section 2.2.3 we handle the general case of arbitrary
matrices (needed for SDP-solving); the state-preparation part will be
substantially more involved there, because in the general case we need not
know the diagonalizing bases for $A$ and $H$, and $A$ and $H$ may not be
simultaneously diagonalizable.

####  2.2.1 General approach

To start, consider the following lemma about the multiplicative approximation
error of a ratio of two real numbers that are given by multiplicative
approximations:

######  Lemma 7.

Let $0\leq\theta\leq 1$ and let $\alpha,\tilde{\alpha},Z,\tilde{Z}$ be
positive real numbers such that $|\alpha-\tilde{\alpha}|\leq\alpha\theta/3$
and $|Z-\tilde{Z}|\leq Z\theta/3$. Then

$\left|\frac{\alpha}{Z}-\frac{\tilde{\alpha}}{\tilde{Z}}\right|\leq\theta\frac{\alpha}{Z}$

###### Proof.

Observe $|\tilde{Z}|\geq|Z|2/3$, thus

$\displaystyle\left|\frac{\alpha}{Z}-\frac{\tilde{\alpha}}{\tilde{Z}}\right|$
$\displaystyle=\left|\frac{\alpha\tilde{Z}-\tilde{\alpha}Z}{Z\tilde{Z}}\right|=\left|\frac{\alpha\tilde{Z}-\alpha
Z+\alpha Z-\tilde{\alpha}Z}{Z\tilde{Z}}\right|$
$\displaystyle\leq\left|\frac{\alpha\tilde{Z}-\alpha
Z}{Z\tilde{Z}}\right|+\left|\frac{\alpha
Z-\tilde{\alpha}Z}{Z\tilde{Z}}\right|\leq\frac{\alpha}{Z}\left|\frac{\tilde{Z}-Z}{\tilde{Z}}\right|+\left|\frac{\theta\alpha}{3\tilde{Z}}\right|$
$\displaystyle\leq\frac{3}{2}\left(\frac{\alpha}{Z}\left|\frac{\tilde{Z}-Z}{Z}\right|+\frac{\theta\alpha}{3Z}\right)\leq\frac{3}{2}\left(\frac{2\theta}{3}\frac{\alpha}{Z}\right)=\theta\frac{\alpha}{Z}.$

∎

######  Corollary 8.

Let $A$ be such that $\left\lVert A\right\rVert\leq 1$. A multiplicative
$\theta/9$-approximation of both $\mbox{\rm
Tr}\left(\frac{I+A/2}{4}e^{-H}\right)$ and $\mbox{\rm
Tr}\left(\frac{I}{4}e^{-H}\right)$ suffices to get an additive
$\theta$-approximation of $\frac{\mbox{\rm Tr}\left(Ae^{-H}\right)}{\mbox{\rm
Tr}\left(e^{-H}\right)}$.

###### Proof.

According to Lemma 7 by dividing the two multiplicative approximations we get

$\frac{\theta}{3}\frac{\mbox{\rm
Tr}\left(\frac{I+A/2}{4}e^{-H}\right)}{\mbox{\rm
Tr}\left(\frac{I}{4}e^{-H}\right)}=\frac{\theta}{3}\left(1+\frac{\mbox{\rm
Tr}\left(\frac{A}{2}e^{-H}\right)}{\mbox{\rm
Tr}\left(e^{-H}\right)}\right)\leq\frac{\theta}{3}\left(1+\frac{\left\lVert
A\right\rVert}{2}\right)\leq\theta/2,$

i.e., an additive $\theta/2$-approximation of

$1+\frac{\mbox{\rm Tr}\left(\frac{A}{2}e^{-H}\right)}{\mbox{\rm
Tr}\left(e^{-H}\right)},$

which yields an additive $\theta$-approximation to $\mbox{\rm
Tr}\left(A\rho\right)$. ∎

It thus suffices to approximate both quantities from the corollary separately.
Notice that both are of the form $\mbox{\rm
Tr}\left(\frac{I+A/2}{4}e^{-H}\right)$, the first with the actual $A$, the
second with $A=0$. Furthermore, a multiplicative $\theta/9$-approximation to
both can be achieved by approximating both up to an additive error
$\theta\mbox{\rm Tr}\left(e^{-H}\right)/72$, since $\mbox{\rm
Tr}\left(\frac{I}{8}e^{-H}\right)\leq\mbox{\rm
Tr}\left(\frac{I+A/2}{4}e^{-H}\right)$.

For now, let us assume we can construct a unitary $U_{A,H}$ such that if we
apply it to the state $|0\ldots 0\rangle$ then we get a probability
$\frac{\mbox{\rm Tr}\left((I+A/2)e^{-H}\right)}{4n}$ of outcome 0 when
measuring the first qubit. That is:

$\left\lVert(\langle 0|\otimes I)U_{A,H}|0\ldots
0\rangle\right\rVert^{2}=\frac{\mbox{\rm Tr}\left((I+A/2)e^{-H}\right)}{4n}.$

In practice we will not be able to construct such a $U_{A,H}$ exactly, instead
we will construct a $\tilde{U}_{A,H}$ that yields a sufficiently close
approximation of the correct probability. Since we have access to such a
unitary, the following lemma allows us to use amplitude estimation to estimate
the probability and hence $\mbox{\rm Tr}\left(\frac{I+A/2}{4}e^{-H}\right)$ up
to the desired error.

######  Lemma 9.

Suppose we have a unitary $U$ acting on $q$ qubits such that $U|0\ldots
0\rangle=|0\rangle|\psi\rangle+|\Phi\rangle$ with $(\langle 0|\otimes
I)|\Phi\rangle=0$ and $\left\lVert\psi\right\rVert^{2}=p\geq p_{\min}$ for
some known bound $p_{\min}$. Let $\mu\in(0,1]$ be the allowed multiplicative
error in our estimation of $p$. Then with
${\mathcal{O}}\left(\frac{1}{\mu\sqrt{p_{\min}}}\right)$ uses of $U$ and
$U^{-1}$ and using ${\mathcal{O}}\left(\frac{q}{\mu\sqrt{p_{\min}}}\right)$
gates on the $q$ qubits we obtain a $\tilde{p}$ such that
$|p-\tilde{p}|\leq\mu p$ with probability at least $4/5$.

###### Proof.

We use the amplitude-estimation algorithm of [BHMT02, Theorem 12] with $M$
applications of $U$ and $U^{-1}$. This provides an estimate $\tilde{p}$ of
$p$, that with probability at least $8/\pi^{2}>4/5$ satisfies

$|p-\tilde{p}|\leq
2\pi\frac{\sqrt{p(1-p)}}{M}+\frac{\pi^{2}}{M^{2}}\leq\frac{\pi}{M}\left(2\sqrt{p}+\frac{\pi}{M}\right).$

Choosing $M$ the smallest power of $2$ such that $M\geq
3\pi/(\mu\sqrt{p_{\min}})$, with probability at least $4/5$ we get

$|p-\tilde{p}|\leq\mu\frac{\sqrt{p_{\min}}}{3}\left(2\sqrt{p}+\mu\frac{\sqrt{p_{\min}}}{3}\right)\leq\mu\frac{\sqrt{p}}{3}\left(3\sqrt{p}\right)\leq\mu
p.$

The $q$ factor in the gate complexity comes from the implementation of the
amplitude amplification steps needed in amplitude-estimation. The gate
complexity of the whole amplitude-estimation procedure is dominated by this
contribution, proving the final gate complexity. ∎

######  Corollary 10.

Suppose we are given the positive numbers $z\leq\mbox{\rm
Tr}\left(e^{-H}\right)$, $\theta\in(0,1]$, and unitary circuits
$\tilde{U}_{A^{\prime},H}$ for $A^{\prime}=0$ and $A^{\prime}=A$ with
$\left\lVert A\right\rVert\leq 1$, each acting on at most $q$ qubits such that

$\left|\left\lVert(\langle 0|\otimes I)\tilde{U}_{A^{\prime},H}|0\ldots
0\rangle\right\rVert^{2}-\frac{\mbox{\rm
Tr}\left((I+A^{\prime}/2)e^{-H}\right)}{4n}\right|\leq\frac{\theta z}{144n}.$

(To clarify the notation: if $|\psi\rangle$ is a 2-register state, then
$(\langle 0|\otimes I)|\psi\rangle$ is the (unnormalized) state in the 2nd
register that results from projecting on $|0\rangle$ in the 1st register.)

Applying the procedure of Lemma 9 to $\tilde{U}_{A^{\prime},H}$ (both for
$A^{\prime}=0$ and for $A^{\prime}=A$) with $p_{\min}=\frac{z}{9n}$ and
$\mu=\theta/19$, and combining the results using Corollary 8 yields an
additive $\theta$-approximation of $\mbox{\rm Tr}\left(A\rho\right)$ with high
probability. The procedure uses

${\mathcal{O}}\left(\frac{1}{\theta}\sqrt{\frac{n}{z}}\right)$

applications of $\tilde{U}_{A,H}$, $\tilde{U}_{0,H}$ and their inverses, and
${\mathcal{O}}\left(\frac{q}{\theta}\sqrt{\frac{n}{z}}\right)$ additional
gates.

###### Proof.

First note that since $I+A^{\prime}/2\succeq I/2$ we have

$p:=\frac{\mbox{\rm
Tr}\left((I+A^{\prime}/2)e^{-H}\right)}{4n}\geq\frac{\mbox{\rm
Tr}\left(e^{-H}\right)}{8n}$

and thus

$\left|\left\lVert(\langle 0|\otimes I)\tilde{U}_{A^{\prime},H}|0\ldots
0\rangle\right\rVert^{2}-p\right|\leq\frac{\theta
z}{144n}\leq\frac{\theta}{18}\cdot\frac{\mbox{\rm
Tr}\left(e^{-H}\right)}{8n}\leq\frac{\theta}{18}p\leq\frac{p}{18}.$ (6)

Therefore

$\left\lVert(\langle 0|\otimes I)\tilde{U}_{A^{\prime},H}|0\ldots
0\rangle\right\rVert^{2}\geq\left(1-\frac{1}{18}\right)p\geq\left(1-\frac{1}{18}\right)\frac{\mbox{\rm
Tr}\left(e^{-H}\right)}{8n}>\frac{\mbox{\rm
Tr}\left(e^{-H}\right)}{9n}\geq\frac{z}{9n}=p_{\min}.$

Also by (6) we have

$\left\lVert(\langle 0|\otimes I)\tilde{U}_{A^{\prime},H}|0\ldots
0\rangle\right\rVert^{2}\leq\left(1+\frac{\theta}{18}\right)p\leq\frac{19}{18}p.$

Therefore using Lemma 9, with $\mu=\theta/19$, with high probability we get a
$\tilde{p}$ satisfying

$\left|\tilde{p}-\left\lVert(\langle 0|\otimes
I)\tilde{U}_{A^{\prime},H}|0\ldots
0\rangle\right\rVert^{2}\right|\leq\frac{\theta}{19}\cdot\left\lVert(\langle
0|\otimes I)\tilde{U}_{A^{\prime},H}|0\ldots
0\rangle\right\rVert^{2}\leq\frac{\theta}{18}p.$ (7)

By combining (6)-(7) using the triangle inequality we get

$\left|p-\tilde{p}\right|\leq\frac{\theta}{9}p,$

so that Corollary 8 can indeed be applied. The complexity statement follows
from Lemma 9 and our choices of $p_{\min}$ and $\mu$. ∎

Notice the $1/\sqrt{z}\geq 1/\sqrt{\mbox{\rm Tr}\left(e^{-H}\right)}$ factor
in the complexity statement of the last corollary. To make sure this factor is
not too large, we would like to ensure $\mbox{\rm
Tr}\left(e^{-H}\right)=\Omega(1)$. This can be achieved by substituting
$H_{+}=H-\lambda_{\min}I$, where $\lambda_{\min}$ is the smallest eigenvalue
of $H$. It is easy to verify that this will not change the value $\mbox{\rm
Tr}\left(Ae^{-H}/\mbox{\rm Tr}\left(e^{-H}\right)\right)$.

It remains to show how to compute $\lambda_{\min}$ and how to apply
$\tilde{U}_{A,H}$. Both of these steps are considerably easier in the case
where all matrices are diagonal, so we will consider this case first.

####  2.2.2 The special case of diagonal matrices – for LP-solving

In this section we consider diagonal matrices, assuming oracle access to $H$
of the following form:

$O_{H}|i\rangle|z\rangle=|i\rangle|z\oplus H_{ii}\rangle$

and similarly for $A$. Notice that this kind of oracle can easily be
constructed from the general sparse matrix oracle (5) that we assume access
to.

######  Lemma 11.

Let $A,H\in\mathbb{R}^{n\times n}$ be diagonal matrices such that $\left\lVert
A\right\rVert\leq 1$ and $H\succeq 0$, and let $\mu>0$ be an error parameter.
Then there exists a unitary $\tilde{U}_{A,H}$ such that

$\left|\left\lVert(\langle 0|\otimes I)\tilde{U}_{A,H}|0\ldots
0\rangle\right\rVert^{2}-\mbox{\rm
Tr}\left(\frac{I+A/2}{4n}e^{-H}\right)\right|\leq\mu,$

which uses $1$ quantum query to $A$ and $H$ and
${\mathcal{O}}(\log^{{\mathcal{O}}(1)}(1/\mu)+\log(n))$ other gates.

###### Proof.

For simplicity assume $n$ is a power of two. This restriction is not
necessary, but makes the proof a bit simpler to state.

In the first step we prepare the state $\sum_{i=1}^{n}|i\rangle/\sqrt{n}$
using $\log(n)$ Hadamard gates on $|0\rangle^{\otimes\log(n)}$. Then we query
the diagonal values of $H$ and $A$ to get the state
$\sum_{i=1}^{n}|i\rangle|H_{ii}\rangle|A_{ii}\rangle/\sqrt{n}$. Using these
binary values we apply a finite-precision arithmetic circuit to prepare

$\frac{1}{\sqrt{n}}\sum_{i=1}^{n}|i\rangle|H_{ii}\rangle|A_{ii}\rangle|\beta_{i}\rangle\text{,
where
}\beta_{i}:=\arcsin\left(\sqrt{\frac{1+A_{ii}/2}{4}e^{-H_{ii}}+\delta_{i}}\right)/\pi\text{,
and }|\delta_{i}|\leq\mu.$

Note that the error $\delta_{i}$ comes from writing down only a finite number
of bits $b_{1}.b_{2}b_{3}\dots b_{\log(8/\mu)}$. Due to our choice of $A$ and
$H$, we know that $\beta_{i}$ lies in $[0,1]$. We proceed by first adding an
ancilla qubit initialized to $|1\rangle$ in front of the state, then we apply
$\log(8/\mu)$ controlled rotations to this qubit: for each $b_{j}=1$ we apply
a rotation by angle $\pi 2^{-j}$. In other words, if $b_{1}=1$, then we rotate
$|1\rangle$ fully to $|0\rangle$. If $b_{2}=1$, then we rotate halfway, and we
proceed further by halving the angle for each subsequent bit. We will end up
with the state:

$\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left(\sqrt{\frac{1+A_{ii}/2}{4}e^{-H_{ii}}+\delta_{i}}|0\rangle+\sqrt{1-\frac{1+A_{ii}/2}{4}e^{-H_{ii}}-\delta_{i}}|1\rangle\right)|i\rangle|A_{ii}\rangle|H_{ii}\rangle|\beta_{i}\rangle.$

It is now easy to see that the squared norm of the $|0\rangle$-part of this
state is as required:

$\left\lVert\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\sqrt{\frac{1+A_{ii}/2}{4}e^{-H_{ii}}+\delta_{i}}|i\rangle\right\rVert^{2}=\frac{1}{n}\sum_{i=1}^{n}\left(\frac{1+A_{ii}/2}{4}e^{-H_{ii}}+\delta_{i}\right)=\frac{\mbox{\rm
Tr}\left((I+A/2)e^{-H}\right)}{4n}+\sum_{i=1}^{n}\frac{\delta_{i}}{n},$

which is an additive $\mu$-approximation since
$\left|\sum_{i=1}^{n}\frac{\delta_{i}}{n}\right|\leq\mu$. ∎

######  Corollary 12.

Let $A,H\in\mathbb{R}^{n\times n}$ be diagonal matrices, with $\left\lVert
A\right\rVert\leq 1$. An additive $\theta$-approximation of

$\mbox{\rm Tr}\left(A\rho\right)=\frac{\mbox{\rm
Tr}\left(Ae^{-H}\right)}{\mbox{\rm Tr}\left(e^{-H}\right)}$

can be computed using ${\mathcal{O}}\left(\frac{\sqrt{n}}{\theta}\right)$
queries to $A$ and $H$ and
$\widetilde{\mathcal{O}}\left(\frac{\sqrt{n}}{\theta}\right)$ other gates.

###### Proof.

Since $H$ is a diagonal matrix, its eigenvalues are exactly its diagonal
entries. Using the quantum minimum-finding algorithm of Dürr and Høyer [DH96]
one can find (with high success probability) the minimum $\lambda_{\min}$ of
the diagonal entries using ${\mathcal{O}}(\sqrt{n})$ queries to the matrix
elements. Applying Lemma 11 and Corollary 10 to $H_{+}=H-\lambda_{\min}I$,
with $z=1$, gives the stated bound. ∎

####  2.2.3 General case – for SDP-solving

In this section we will extend the ideas from the last section to non-diagonal
matrices. There are a few complications that arise in this more general case.
These mostly follow from the fact that we now do not know the eigenvectors of
$H$ and $A$, which were the basis states before, and that these eigenvectors
might not be the same for both matrices. For example, to find the minimal
eigenvalue of $H$, we can no longer simply minimize over its diagonal entries.
To solve this, in Appendix C we develop new techniques that generalize
minimum-finding.

Furthermore, the unitary $\tilde{U}_{A,H}$ in the LP case could be seen as
applying the operator

$\sqrt{\frac{I+A/2}{4}e^{-H}}$

to a superposition of its eigenvectors. This is also more complicated in the
general setting, due to the fact that the eigenvectors are no longer the basis
states. In Appendix B we develop general techniques to apply smooth functions
of Hamiltonians to a state. Among other things, this will be used to create an
efficient purified Gibbs sampler.

Our Gibbs sampler uses similar methods to the work of Chowdhury and Somma
[CS16] for achieving logarithmic dependence on the precision. However, the
result of [CS16] cannot be applied to our setting, because it implicitly
assumes access to an oracle for $\sqrt{H}$ instead of $H$. Although their
paper describes a way to construct such an oracle, it comes with a large
overhead: they construct an oracle for $\sqrt{H^{\prime}}=\sqrt{H+\nu I}$,
where $\nu\in\mathbb{R}_{+}$ is some possibly large positive number. This
shifting can have a huge effect on $Z^{\prime}=\mbox{\rm
Tr}\left(e^{-H^{\prime}}\right)=e^{-\nu}\mbox{\rm Tr}\left(e^{-H}\right)$,
which can be prohibitive due to the $\sqrt{1/Z^{\prime}}$ factor in the
runtime, which blows up exponentially in $\nu$.

In the following lemma we show how to implement $\tilde{U}_{A,H}$ using the
techniques we developed in Appendix B.

######  Lemma 13.

Let $A,H\in\mathbb{C}^{n\times n}$ be Hermitian matrices such that
$\left\lVert A\right\rVert\leq 1$ and $I\preceq H\preceq KI$ for a known
$K\in\mathbb{R}_{+}$. Assume $A$ is $s$-sparse and $H$ is $d$-sparse with
$s\leq d$. Let $\mu>0$ be an error parameter. Then there exists a unitary
$\tilde{U}_{A,H}$ such that

$\left|\left\lVert(\langle 0|\otimes I)\tilde{U}_{A,H}|0\ldots
0\rangle\right\rVert^{2}-\mbox{\rm
Tr}\left(\frac{I+A/2}{4n}e^{-H}\right)\right|\leq\mu$

that uses $\widetilde{\mathcal{O}}\left(Kd\right)$ queries to $A$ and $H$, and
the same order of other gates.

###### Proof.

The basic idea is that we first prepare a maximally entangled state
$\sum_{i=1}^{n}|i\rangle|i\rangle/\sqrt{n}$, and then apply the (norm-
decreasing) maps $e^{-H/2}$ and $\sqrt{\frac{I+A/2}{4}}$ to the first
register. Note that we can assume without loss of generality that $\mu\leq 1$,
otherwise the statement is trivial.

Let $\tilde{W}_{0}=\left(\langle 0|\otimes
I\right)\tilde{W}\left(|0\rangle\otimes I\right)$ be a $\mu/5$-approximation
of the map $e^{-H/2}$ (in operator norm) implemented by using Theorem 43, and
let $\tilde{V}_{0}=\left(\langle 0|\otimes
I\right)\tilde{V}\left(|0\rangle\otimes I\right)$ be a $\mu/5$-approximation
of the map $\sqrt{\frac{I+A/2}{4}}$ implemented by using Theorem 40. We define
$\tilde{U}_{A,H}:=\tilde{V}\tilde{W}$, noting that there is a hidden $\otimes
I$ factor in both $\tilde{V}$ and $\tilde{W}$ corresponding to each other’s
ancilla qubit. As in the linear programming case, we are interested in $p$,
the probability of measuring a $00$ in the first register (i.e., the two
“flag” qubits) after applying $\tilde{U}_{A,H}$. We will analyze this in terms
of these operators below. We will make the final approximation step precise in
the next paragraph.

$\displaystyle p^{\prime}$ $\displaystyle:=\left\lVert\left(\langle 00|\otimes
I\right)\tilde{U}_{A,H}\left(|00\rangle\otimes
I\right)\sum_{i=1}^{n}\frac{|i\rangle|i\rangle}{\sqrt{n}}\right\rVert^{2}$
$\displaystyle=\left\lVert\tilde{V}_{0}\tilde{W}_{0}\sum_{i=1}^{n}\frac{|i\rangle|i\rangle}{\sqrt{n}}\right\rVert^{2}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\langle
i|\tilde{W}_{0}^{\dagger}\tilde{V}_{0}^{\dagger}\tilde{V}_{0}\tilde{W}_{0}|i\rangle$
$\displaystyle=\frac{1}{n}\mbox{\rm
Tr}\left(\tilde{W}_{0}^{\dagger}\tilde{V}_{0}^{\dagger}\tilde{V}_{0}\tilde{W}_{0}\right)$
$\displaystyle=\frac{1}{n}\mbox{\rm
Tr}\left(\tilde{V}_{0}^{\dagger}\tilde{V}_{0}\tilde{W}_{0}\tilde{W}_{0}^{\dagger}\right)$
(8) $\displaystyle\approx\frac{1}{n}\mbox{\rm
Tr}\left(\frac{I+A/2}{4}e^{-H}\right).$ (9)

Note that for all matrices $B,\tilde{B}$ with $\left\lVert B\right\rVert\leq
1$, we have

$\displaystyle\left\lVert
B^{\dagger}B-\tilde{B}^{\dagger}\tilde{B}\right\rVert$
$\displaystyle=\left\lVert(B^{\dagger}-\tilde{B}^{\dagger})B+B^{\dagger}(B-\tilde{B})-(B^{\dagger}-\tilde{B}^{\dagger})(B-\tilde{B})\right\rVert$
$\displaystyle\leq\left\lVert(B^{\dagger}-\tilde{B}^{\dagger})B\right\rVert+\left\lVert
B^{\dagger}(B-\tilde{B})\right\rVert+\left\lVert(B^{\dagger}-\tilde{B}^{\dagger})(B-\tilde{B})\right\rVert$
$\displaystyle\leq\left\lVert
B^{\dagger}-\tilde{B}^{\dagger}\right\rVert\left\lVert
B\right\rVert+\left\lVert B^{\dagger}\right\rVert\left\lVert
B-\tilde{B}\right\rVert+\left\lVert
B^{\dagger}-\tilde{B}^{\dagger}\right\rVert\left\lVert
B-\tilde{B}\right\rVert$ $\displaystyle\leq 2\left\lVert
B-\tilde{B}\right\rVert+\left\lVert B-\tilde{B}\right\rVert^{2}.$

Since $\mu\leq 1$, and hence $2\mu/5+(\mu/5)^{2}\leq\mu/2$, this implies (with
$B=e^{-H/2}$ and $\tilde{B}=\tilde{W}_{0}^{\dagger}$) that $\left\lVert
e^{-H}-\tilde{W}_{0}\tilde{W}_{0}^{\dagger}\right\rVert\leq\mu/2$, and also
(with $B=\sqrt{(I+A/2)/4}$ and $\tilde{B}=\tilde{V}_{0}$)
$\left\lVert(I+A/2)/4-\tilde{V}_{0}^{\dagger}\tilde{V}_{0}\right\rVert\leq\mu/2$.
Let $\left\lVert\cdot\right\rVert_{1}$ denote the trace norm (a.k.a. Schatten
1-norm). Note that for all $C,D,\tilde{C},\tilde{D}$:

$\displaystyle\left|\mbox{\rm Tr}\left(CD\right)-\mbox{\rm
Tr}\left(\tilde{C}\tilde{D}\right)\right|$ $\displaystyle\leq\left\lVert
CD-\tilde{C}\tilde{D}\right\rVert_{1}$
$\displaystyle=\left\lVert(C-\tilde{C})D+C(D-\tilde{D})-(C-\tilde{C})(D-\tilde{D})\right\rVert_{1}$
$\displaystyle\leq\left\lVert(C-\tilde{C})D\right\rVert_{1}+\left\lVert
C(D-\tilde{D})\right\rVert_{1}+\left\lVert(C-\tilde{C})(D-\tilde{D})\right\rVert_{1}$
$\displaystyle\leq\left\lVert C-\tilde{C}\right\rVert\left\lVert
D\right\rVert_{1}+\left\lVert D-\tilde{D}\right\rVert\left(\left\lVert
C\right\rVert_{1}+\left\lVert C-\tilde{C}\right\rVert_{1}\right).$

Which, in our case (setting $C=(I+A/2)/4$, $D=e^{-H}$,
$\tilde{C}=\tilde{V}_{0}^{\dagger}\tilde{V}_{0}$, and
$\tilde{D}=\tilde{W}_{0}\tilde{W}_{0}^{\dagger}$) implies that

$\left|\mbox{\rm Tr}\left(\left(I+A/2\right)e^{-H}/4\right)-\mbox{\rm
Tr}\left(\tilde{V}_{0}^{\dagger}\tilde{V}_{0}\tilde{W}_{0}\tilde{W}_{0}^{\dagger}\right)\right|\leq(\mu/2)\mbox{\rm
Tr}\left(e^{-H}\right)+(\mu/2)(1/2+\mu/2)n.$

Dividing both sides by $n$ and using equation (8) then implies

$\displaystyle\left|\mbox{\rm
Tr}\left(\left(I+A/2\right)e^{-H}\right)/(4n)-p^{\prime}\right|$
$\displaystyle\leq\frac{\mu}{2}\frac{\mbox{\rm
Tr}\left(e^{-H}\right)}{n}+\frac{\mu}{2}\left(\frac{1}{2}+\frac{\mu}{2}\right)$
$\displaystyle\leq\frac{\mu}{2}+\frac{\mu}{2}$ $\displaystyle=\mu.$

This proves the correctness of $\tilde{U}_{A,H}$. It remains to show that the
complexity statement holds. To show this we only need to specify how to
implement the map $\sqrt{\frac{I+A/2}{4}}$ using Theorem 40, since the map
$e^{H/2}$ is already dealt with in Theorem 43. To use Theorem 40, we choose
$x_{0}:=0$, $K:=1$ and $r:=1$, since $\left\lVert A\right\rVert\leq 1$.
Observe that
$\sqrt{\frac{1+x/2}{4}}=\frac{1}{2}\sum_{k=0}^{\infty}\binom{1/2}{k}\left(\frac{x}{2}\right)^{k}$
whenever $|x|\leq 1$. Also let $\delta=1/2$, so $r+\delta=\frac{3}{2}$ and
$\frac{1}{2}\sum_{k=0}^{\infty}\left|\binom{1/2}{k}\right|\left(\frac{3}{4}\right)^{k}\leq
1=:B$. Recall that $\tilde{V}$ denotes the unitary that Theorem 40 constructs.
Since we choose the precision parameter to be $\mu/5=\Theta(\mu)$, Theorem 40
shows $\tilde{V}$ can be implemented using
${\mathcal{O}}\left(d\log^{2}\left(1/\mu\right)\right)$ queries and
${\mathcal{O}}\left(d\log^{2}\left(1/\mu\right)\left[\log(n)+\log^{2.5}\left(1/\mu\right)\right]\right)$
gates. This cost is negligible compared to our implementation cost of
$e^{-H/2}$ with $\mu/5$ precision: Theorem 43 uses
${\mathcal{O}}\left(Kd\log^{2}\left(K/\mu\right)\right)$ queries and
${\mathcal{O}}\left(Kd\log^{2}\left(Kd/\mu\right)\left[\log(n)+\log^{2.5}\left(Kd/\mu\right)\right]\right)$
gates to implement $\tilde{W}$. ∎

######  Corollary 14.

Let $A,H\in\mathbb{C}^{n\times n}$ be Hermitian matrices such that
$\left\lVert A\right\rVert\leq 1$ and $\left\lVert H\right\rVert\leq K$ for a
known bound $K\in\mathbb{R}_{+}$. Assume $A$ is $s$-sparse and $H$ is
$d$-sparse with $s\leq d$. An additive $\theta$-approximation of

$\mbox{\rm Tr}\left(A\rho\right)=\frac{Ae^{-H}}{\mbox{\rm
Tr}\left(e^{-H}\right)}$

can be computed using
$\widetilde{\mathcal{O}}\left(\frac{\sqrt{n}dK}{\theta}\right)$ queries to $A$
and $H$, while using the same order of other gates.

###### Proof.

Start by computing an estimate $\tilde{\lambda}_{\min}$ of
$\lambda_{\min}(H)$, the minimum eigenvalue of $H$, up to additive error
$\varepsilon=1/2$ using Lemma 50. We define
$H_{+}:=H-(\tilde{\lambda}_{\min}-3/2)I$, so that $I\preceq H_{+}$ but
$2I\nprec H_{+}$.

Applying Lemma 13 and Corollary 10 to $H_{+}$ with $z=e^{-2}$ gives the stated
bound. ∎

###  2.3 An efficient 2-sparse oracle

Recall from the end of Section 2.1 that $\tilde{a}_{j}$ is an additive
$\theta$-approximation to $\mbox{\rm Tr}\left(A_{j}\rho\right)$, $\tilde{c}$
is a $\theta$-approximation to $\mbox{\rm Tr}\left(C\rho\right)$ and
$c^{\prime}=\tilde{c}-r\theta-\theta$. We first describe our quantum 2-sparse
oracle assuming access to a unitary which acts as
$|j\rangle|0\rangle|0\rangle\mapsto|j\rangle|\tilde{a}_{j}\rangle|\psi_{j}\rangle$,
where $|\psi_{j}\rangle$ is some workspace state depending on $j$. We then
briefly discuss how to modify the analysis when we are given an oracle which
acts as
$|j\rangle|0\rangle|0\rangle\mapsto|j\rangle\sum_{i}\beta_{j}^{i}|\tilde{a}^{i}_{j}\rangle|\psi^{i}_{j}\rangle$
(where each $\tilde{a}^{i}_{j}$ is an additive $\theta$-approximation to
$\mbox{\rm Tr}\left(A_{j}\rho\right)$), since this is the output of the trace-
estimation procedure of the previous section.

Our goal is to find a $y\in\tilde{\mathcal{P}}(\tilde{a},c^{\prime})$, i.e., a
$y$ such that

$\displaystyle\left\lVert y\right\rVert_{1}$ $\displaystyle\leq r$
$\displaystyle b^{T}y$ $\displaystyle\leq\alpha$ $\displaystyle\tilde{a}^{T}y$
$\displaystyle\geq c^{\prime}$ $\displaystyle y$ $\displaystyle\geq 0$

If $\alpha\geq 0$ and $c^{\prime}\leq 0$, then $y=0$ is a solution and our
oracle can return it. If not, then we may write $y=Nq$ with $N=\left\lVert
y\right\rVert_{1}>0$ and hence $\left\lVert q\right\rVert_{1}=1$. So we are
looking for an $N$ and a $q$ such that

$\displaystyle b^{T}q$ $\displaystyle\leq\alpha/N$ (10)
$\displaystyle\tilde{a}^{T}q$ $\displaystyle\geq c^{\prime}/N$
$\displaystyle\left\lVert q\right\rVert_{1}$ $\displaystyle=1$ $\displaystyle
q$ $\displaystyle\geq 0$ $\displaystyle 0$ $\displaystyle<N\leq r$

We can now view $q\in\mathbb{R}_{+}^{m}$ as the coefficients of a convex
combination of the points $p_{i}=(b_{i},\tilde{a}_{i})$ in the plane. We want
such a combination that lies to the upper left of
$g_{N}=(\alpha/N,c^{\prime}/N)$ for some $0<N\leq r$. Let $\mathcal{G}_{N}$
denote the upper-left quadrant of the plane starting at $g_{N}$.

######  Lemma 15.

If there is a $y\in\tilde{\mathcal{P}}(\tilde{a},c^{\prime})$, then there is a
$2$-sparse $y^{\prime}\in\tilde{\mathcal{P}}(\tilde{a},c^{\prime})$ such that
$\left\lVert y\right\rVert_{1}=\left\lVert y^{\prime}\right\rVert_{1}$.

###### Proof.

Consider $p_{i}=(b_{i},\tilde{a}_{i})$ and $g=(\alpha/N,c^{\prime}/N)$ as
before, and write $y=Nq$ where $\sum_{j=1}^{m}q_{j}=1$, $q\geq 0$. The vector
$q$ certifies that a convex combination of the points $p_{i}$ lies in
$\mathcal{G}_{N}$. But then there exist $j,k\in[m]$ such that the line segment
$\overline{p_{j}p_{k}}$ intersects $\mathcal{G}_{N}$. All points on this line
segment are convex combinations of $p_{j}$ and $p_{k}$, hence there is a
convex combination of $p_{j}$ and $p_{k}$ that lies in $\mathcal{G}_{N}$. This
gives a $2$-sparse $q^{\prime}$, and
$y^{\prime}=Nq^{\prime}\in\tilde{\mathcal{P}}(\tilde{a},c^{\prime})$. ∎

We can now restrict our search to $2$-sparse $y$. Let
$\mathcal{G}=\bigcup_{N\in(0,r]}\mathcal{G}_{N}$. Then we want to find two
points $p_{j},p_{k}$ that have a convex combination in $\mathcal{G}$, since
this implies that a scaled version of their convex combination gives a
$y\in\tilde{\mathcal{P}}(\tilde{a},c^{\prime})$ and $\left\lVert
y\right\rVert_{1}\leq r$.

######  Lemma 16.

There is an oracle that returns a 2-sparse vector
$y\in\tilde{\mathcal{P}}(\tilde{a},c^{\prime})$, if one exists, using one
search and two minimizations over the $m$ points
$p_{j}=(b_{j},\tilde{a}_{j})$. This gives a classical algorithm that uses
${\mathcal{O}}(m)$ calls to the subroutine that gives the entries of
$\tilde{a}$ and ${\mathcal{O}}(m)$ other operations, and a quantum algorithm
that (in order to solve the problems with high probability) uses
${\mathcal{O}}(\sqrt{m})$ calls to the subroutine that gives the entries of
$\tilde{a}$ and $\widetilde{\mathcal{O}}\left(\sqrt{m}\right)$ other gates.

###### Proof.

The algorithm can be summarized as follows:

  1. 1.

Check if $\alpha\geq 0$ and $c^{\prime}\leq 0$. If so, output $y=0$.

  2. 2.

Check if there is a $p_{i}\in\mathcal{G}$. If so, let $q=e_{i}$ and
$N=\frac{c^{\prime}}{\tilde{a}_{i}}$.

  3. 3.

Find $p_{j},p_{k}$ so that the line segment $\overline{p_{j}p_{k}}$ goes
through $\mathcal{G}$. This gives coefficients $q$ of a convex combination
that can be scaled by $N=\frac{c^{\prime}}{\tilde{a}^{T}q}$ to give $y$. The
main realization is that we can search separately for $p_{j}$ and $p_{k}$.

First we will need a better understanding of the shape of $\mathcal{G}$ (see
Figure 1 for illustration). This depends on the sign of $\alpha$ and
$c^{\prime}$. If we define $\mbox{\rm sign}(0)=1$:

  * (a)

If $\mbox{\rm sign}(\alpha)=-1,\mbox{\rm sign}(c^{\prime})=-1$. The corner
point of $\mathcal{G}$ is $(\alpha/r,c^{\prime}/r)$. One edge goes up
vertically and an other follows the line segment
$\lambda\cdot(\alpha,c^{\prime})$ for $\lambda\in[1/r,\infty)$ starting at the
corner.

  * (b)

If $\mbox{\rm sign}(\alpha)=-1,\mbox{\rm sign}(c^{\prime})=1$. Here
$\mathcal{G}_{N}\subseteq\mathcal{G}_{r}$ for $N\leq r$. So
$\mathcal{G}=\mathcal{G}_{r}$. The corner point is again
$(\alpha/r,c^{\prime}/r)$, but now one edge goes up vertically and one goes to
the left horizontally.

  * (c)

If $\mbox{\rm sign}(\alpha)=1,\mbox{\rm sign}(c^{\prime})=-1$. This is the
case where $y=0$ is a solution, $\mathcal{G}$ is the whole plane and has no
corner.

  * (d)

If $\mbox{\rm sign}(\alpha)=1,\mbox{\rm sign}(c^{\prime})=1$. The corner point
of $\mathcal{G}$ is again $(\alpha/r,c^{\prime}/r)$. From there one edge goes
to the left horizontally and one edge follows the line segment
$\lambda\cdot(\alpha,c^{\prime})$ for $\lambda\in[1/r,\infty)$.

(a) $\mbox{\rm sign}(\alpha)=-1,\mbox{\rm sign}(c^{\prime})=-1$

(b) $\mbox{\rm sign}(\alpha)=-1,\mbox{\rm sign}(c^{\prime})=1$

(c) $\mbox{\rm sign}(\alpha)=1,\mbox{\rm sign}(c^{\prime})=-1$

(d) $\mbox{\rm sign}(\alpha)=1,\mbox{\rm sign}(c^{\prime})=1$

Figure 1: The region $\mathcal{G}$ in light blue. The borders of two quadrants
$\mathcal{G}_{N}$ have been drawn by thick dashed blue lines. The red dot at
the beginning of the arrow is the point $(\alpha/r,c^{\prime}/r)$.

Since $\mathcal{G}$ is always an intersection of at most $2$ halfspaces, steps
1-2 of the algorithm are easy to perform. In step 1 we handle case (c) by
simply returning $y=0$. For the other cases $(\alpha/r,c^{\prime}/r)$ is the
corner point of $\mathcal{G}$ and the two edges are simple lines. Hence in
step 2 we can easily search through all the points to find out if there is one
lying in $\mathcal{G}$; since $\mathcal{G}$ is a very simple region, this only
amounts to checking on which side of the two lines a point lies.

$L_{2}$$L_{1}$$p_{j}$$p_{k}$$\angle L_{2}\ell_{k}$$\angle
L_{1}L_{2}$$\angle\ell_{j}L_{1}$ $(\alpha/r,c^{\prime}/r)$ Figure 2:
Illustration of $\mathcal{G}$ with the points $p_{j},p_{k}$ and the angles
$\angle\ell_{j}L_{1},\angle L_{1}L_{2},\angle L_{2}\ell_{k}$ drawn in. Clearly
the line $\overline{p_{j}p_{k}}$ only crosses $\mathcal{G}$ when the total
angle is less than $\pi$.

Now, if we cannot find a single point in $\mathcal{G}$, we need a combination
of two points in step 3. Let $L_{1},L_{2}$ be the edges of $\mathcal{G}$ and
let $\ell_{j}$ and $\ell_{k}$ be the line segments from
$(\alpha/r,c^{\prime}/r)$ to $p_{j}$ and $p_{k}$, respectively. Then, as can
be seen in Figure 2, the line segment $\overline{p_{j}p_{k}}$ goes through
$\mathcal{G}$ if and only if (up to relabeling $p_{j}$ and $p_{k}$)
$\angle\ell_{j}L_{1}+\angle L_{1}L_{2}+\angle L_{2}\ell_{k}\leq\pi$. Since
$\angle L_{1}L_{2}$ is fixed, we can simply look for a $j$ such that
$\angle\ell_{j}L_{1}$ is minimized and a $k$ such that $\angle L_{2}\ell_{k}$
is minimized. If $\overline{p_{j}p_{k}}$ does not pass through $\mathcal{G}$
for this pair of points, then it does not for any of the pairs of points.

Notice that these minimizations can be done separately and hence can be done
in the stated complexity. Given the minimizing points $p_{j}$ and $p_{k}$, it
is easy to check if they give a solution by calculating the angle between
$\ell_{j}$ and $\ell_{k}$. The coefficients of the convex combination $q$ are
then easy to compute. It remains to compute the scaling $N$. This is done by
rewriting the constraints of (10):

$\frac{c^{\prime}}{q^{T}\tilde{a}}\leq N\leq\frac{\alpha}{q^{T}b}$

So we can pick any value in this range for $N$. ∎

The analysis above applies if we are given an oracle which acts as
$|j\rangle|0\rangle|0\rangle\mapsto|j\rangle|\tilde{a}_{j}\rangle|\psi_{j}\rangle$.
However, the trace estimation procedure of Corollary 14 acts as
$|j\rangle|0\rangle|0\rangle\mapsto|j\rangle\sum_{i}\beta_{j}^{i}|\tilde{a}^{i}_{j}\rangle|\psi^{i}_{j}\rangle$
where each $|\tilde{a}^{i}_{j}\rangle$ is an approximation of $a_{j}$ and the
amplitudes $\beta_{j}^{i}$ are such that measuring the second register with
high probability returns an $\tilde{a}^{i}_{j}$ which is $\theta$-close to
$a_{j}$. Since we can exponentially reduce the probability that we obtain an
$\tilde{a}^{i}_{j}$ which is further than $\theta$ away from $a_{j}$, we will
for simplicity assume that for all $i,j$ we have
$|\tilde{a}^{i}_{j}-a_{j}|\leq\theta$; the neglected exponentially small
probabilities will only affect the analysis in negligible ways. Note that
while we do not allow our oracle enough time to obtain classical descriptions
of all $\tilde{a}_{j}$s (we aim for a runtime of
$\widetilde{\mathcal{O}}\left(\sqrt{m}\right)$), we do have enough time to
compute $\tilde{c}$ once initially. Knowing $\tilde{c}$, we can compute the
angles defined by the points $(b_{j},\tilde{a}_{j}^{i})$ with respect to the
corner point of $(\alpha/r,(\tilde{c}-\theta)/r-\theta)$ and the lines
$L_{1},L_{2}$ (see Figure 2). We now apply our generalized minimum-finding
algorithm with runtime $\widetilde{\mathcal{O}}\left(\sqrt{m}\right)$ (see
Theorem 49) starting with a uniform superposition over the $j$s to find
$\ell,k\in[m]$ with points $(\tilde{a}_{\ell},b_{\ell})$ and
$(\tilde{a}_{k},b_{k})$ approximately minimizing the respective angles to
lines $L_{1},L_{2}$. It follows from Lemma 6 and Lemma 15 that if
$\mathcal{P}_{0}(X)\cap\\{y\in\mathbb{R}^{m}:\sum_{j}y_{j}\leq r\\}$ is non-
empty, then some convex combination of $(\tilde{a}_{\ell},b_{\ell})$ and
$(\tilde{a}_{k},b_{k})$ lies in $\mathcal{G}$. On the other hand, if
$\mathcal{P}_{4Rr\theta}(X)\cap\\{y\in\mathbb{R}^{m}:\sum_{j}y_{j}\leq r\\}$
is empty, then by the same lemmas the respective angles will be such that we
correctly conclude that
$\mathcal{P}_{0}(X)\cap\\{y\in\mathbb{R}^{m}:\sum_{j}y_{j}\leq r\\}$ is empty.

###  2.4 Total runtime

We are now ready to add our quantum implementations of the trace calculations
and the oracle to the classical Arora-Kale framework. See 1

###### Proof.

Using our implementations of the different building blocks, it remains to
calculate what the total complexity will be when they are used together.

Cost of the oracle for $H^{(t)}$.

The first problem in each iteration is to obtain access to an oracle for
$H^{(t)}$. In each iteration the oracle will produce a $y^{(t)}$ that is at
most $2$-sparse, and hence in the $(t+1)$th iteration, $H^{(t)}$ is a linear
combination of $2t$ of the $A_{j}$ matrices and the $C$ matrix.

We can write down a sparse representation of the coefficients of the linear
combination that gives $H^{(t)}$ in each iteration by adding the new terms
coming from $y^{(t)}$. This will clearly not take longer than
$\widetilde{\mathcal{O}}\left(T\right)$, since there are only a constant
number of terms to add for our oracle. As we will see, this term will not
dominate the complexity of the full algorithm.

Using such a sparse representation of the coefficients, one query to a sparse
representation of $H^{(t)}$ will cost $\widetilde{\mathcal{O}}\left(st\right)$
queries to the input matrices and $\widetilde{\mathcal{O}}\left(st\right)$
other gates. For a detailed explanation and a matching lower bound, see
Appendix D.

Cost of the oracle for $\mbox{\rm Tr}\left(A_{j}\rho\right)$.

In each iteration $M^{(t)}$ is made to have operator norm at most $1$. This
means that

$\left\lVert-\eta H^{(t)}\right\rVert\leq\eta\sum_{\tau=1}^{t}\left\lVert
M^{(\tau)}\right\rVert\leq\eta t.$

Furthermore we know that $H^{(t)}$ is at most $d:=s(2t+1)$-sparse. Calculating
$\mbox{\rm Tr}\left(A_{j}\rho\right)$ for one index $j$ up to an additive
error of $\theta:=\varepsilon/(12Rr)$ can be done using the algorithm from
Corollary 14. This will take

$\widetilde{\mathcal{O}}\left(\sqrt{n}\frac{\left\lVert H\right\rVert
d}{\theta}\right)=\widetilde{\mathcal{O}}\left(\sqrt{n}s\frac{\eta
t^{2}Rr}{\varepsilon}\right)$

queries to the oracle for $H^{(t)}$ and the same order of other gates. Since
each query to $H^{(t)}$ takes $\widetilde{\mathcal{O}}\left(st\right)$ queries
to the input matrices, this means that

$\widetilde{\mathcal{O}}\left(\sqrt{n}s^{2}\frac{\eta
t^{3}Rr}{\varepsilon}\right)$

queries to the input matrices will be made, and the same order of other gates,
for each approximation of a $\mbox{\rm Tr}\left(A_{j}\rho\right)$ (and
similarly for approximating $\mbox{\rm Tr}\left(C\rho\right)$).

Total cost of one iteration.

Lemma 16 tells us that we will use
$\widetilde{\mathcal{O}}\left(\sqrt{m}\right)$ calculations of $\mbox{\rm
Tr}\left(A_{j}\rho\right)$, and the same order of other gates, to calculate a
classical description of a $2$-sparse $y^{(t)}$. This brings the total cost of
one iteration to

$\widetilde{\mathcal{O}}\left(\sqrt{nm}s^{2}\frac{\eta
t^{3}Rr}{\varepsilon}\right)$

queries to the input matrices, and the same order of other gates.

Total quantum runtime for SDPs.

Since $w\leq r+1$ we can set
$T=\widetilde{\mathcal{O}}\left(\frac{R^{2}r^{2}}{\varepsilon^{2}}\right)$.
With $\eta=\sqrt{\frac{\ln(n)}{T}}$, summing over all iterations in one run of
the algorithm gives a total cost of

$\displaystyle\widetilde{\mathcal{O}}\left(\sum_{t=1}^{T}\sqrt{nm}s^{2}\frac{\eta
t^{3}Rr}{\varepsilon}\right)$
$\displaystyle=\widetilde{\mathcal{O}}\left(\sqrt{nm}s^{2}\frac{\eta
T^{4}Rr}{\varepsilon}\right)$
$\displaystyle=\widetilde{\mathcal{O}}\left(\sqrt{nm}s^{2}\left(\frac{Rr}{\varepsilon}\right)^{\\!\\!8}\right)$

queries to the input matrices and the same order of other gates.

∎

##### Total quantum runtime for LPs.

The final complexity of our algorithm contains a factor
$\widetilde{\mathcal{O}}\left(sT\right)$ that comes from the sparsity of the
$H^{(t)}$ matrix. This assumes that when we add the input matrices together,
the rows become less sparse. This need not happen for certain SDPs. For
example, in the SDP relaxation of MAXCUT, the $H^{(t)}$ will always be
$d$-sparse, where $d$ is the degree of the graph. A more important class of
examples is that of linear programs: since LPs have diagonal $A_{j}$ and $C$,
their sparsity is $s=1$, and even the sparsity of the $H^{(t)}$ is always $1$.
This, plus the fact that the traces can be computed without a factor
$\left\lVert H\right\rVert$ in the complexity (as shown in Corollary 12 in
Section 2.2.2), means that our algorithm solves LPs with

$\widetilde{\mathcal{O}}\left(\sqrt{nm}\left(\frac{Rr}{\varepsilon}\right)^{\\!\\!5}\right)$

queries to the input matrices and the same order of other gates.

##### Total classical runtime.

Using the classical techniques for trace estimation from Appendix A, and the
classical version of our oracle (Lemma 16), we are also able to give a general
classical instantiation of the Arora-Kale framework. The final complexity will
then be

$\widetilde{\mathcal{O}}\left(nms\left(\frac{Rr}{\varepsilon}\right)^{\\!\\!4}+ns\left(\frac{Rr}{\varepsilon}\right)^{\\!\\!7}\right).$

The better dependence on $Rr/\varepsilon$ and $s$, compared to our quantum
algorithm, comes from the fact that we now have the time to write down
intermediate results explicitly. For example, we do not need to recalculate
parts of $H^{(t)}$ for every new query to it, instead we can just calculate it
once at the start of the iteration by adding $M^{(t)}$ to $H^{(t-1)}$ and
writing down the result.

##### Further remarks.

We want to stress again that our solver is meant to work for _all_ SDPs. In
particular, it does not use the structure of a specific SDP. As we show in the
next section, every oracle that works for all SDPs must have large width. To
obtain quantum speedups for a _specific_ class of SDPs, it will be necessary
to develop oracles tuned to that problem. We view this as an important
direction for future work. Recall from the introduction that Arora and Kale
also obtain fast classical algorithms for problems such as MAXCUT by doing
exactly that: they develop specialized oracles for those problems.

##  3 Downside of this method: general oracles are restrictive

In this section we show some of the limitations of a method that uses sparse
or general oracles, i.e., ones that are not optimized for the properties of
specific SDPs. We will start by discussing sparse oracles in the next section.
We will use a counting argument to show that sparse solutions cannot hold too
much information about a problem’s solution. In Section 3.2 we will show that
width-bounds that do not depend on the specific structure of an SDP are for
many problems not efficient. As in the rest of the paper, we will assume the
notation of Section 2, in particular of Meta-Algorithm 1.

###  3.1 Sparse oracles are restrictive

######  Lemma 17.

If, for some specific SDP, every $\varepsilon$-optimal dual-feasible vector
has at least $\ell$ non-zero elements, then the width $w$ of any $k$-sparse
oracle for this SDP is such that
$\frac{Rw}{\varepsilon}=\Omega\left(\sqrt{\frac{\ell}{k\ln(n)}}\right)$.

###### Proof.

The vector $\bar{y}$ returned by Meta-Algorithm 1 is, by construction, the
average of $T$ vectors $y^{(t)}$ that are all $k$-sparse, plus one extra
$1$-sparse term of $\frac{\varepsilon}{R}e_{1}$, and hence $\ell\leq kT+1$.
The stated bound on $\frac{Rw}{\varepsilon}$ then follows directly by
combining this inequality with
$T={\mathcal{O}}(\frac{R^{2}w^{2}}{\varepsilon^{2}}\ln(n))$. ∎

The oracle presented in Section 2.3 always provides a $2$-sparse vector $y$.
This implies that if an SDP requires an $\ell$-sparse dual solution, we must
have $\frac{Rw}{\varepsilon}=\Omega(\sqrt{\ell/\ln(n)})$. This in turn means
that the upper bound on the runtime of our algorithm will be of order
$\ell^{4}\sqrt{nm}s^{2}$. This is clearly bad if $\ell$ is of the order $n$ or
$m$.

Of course it could be the case that almost every (useful) SDP has a sparse
approximate dual solution (or can easily be rewritten so that it does), and
hence sparseness might be not a restriction at all. However, this does not
seem to be the case. We will prove that for certain kinds of SDPs, no useful
dual solution can be very sparse. Let us first define what we mean by
_useful_.

######  Definition 18.

A problem is defined by a function $f$ that, for every element $p$ of the
problem domain $\mathcal{D}$, gives a subset of the solution space
$\mathcal{S}$, consisting of the solutions that are considered correct. We say
a family of SDPs, $\\{SDP^{(p)}\\}_{p\in\mathcal{D}}$, solves the problem via
the dual if there is an $\varepsilon\geq 0$ and a function $g$ such that for
every $p\in\mathcal{D}$ and every $\varepsilon$-optimal dual-feasible vector
$y^{(p)}$ to $SDP^{(p)}$:

$g(y^{(p)})\in f(p).$

In other words, an $\varepsilon$-optimal dual solution can be converted into a
correct solution of the original problem without more knowledge of $p$.

For these kinds of SDP families we will prove a lower bound on the sparsity of
the dual solutions. The idea for this bound is as follows. If you have a lot
of different instances that require different solutions, but the SDPs are
equivalent up to permuting the constraints and the coordinates of
$\mathbb{R}^{n}$, then a dual solution should have a lot of unique
permutations and hence cannot be too sparse.

######  Theorem 19.

Consider a problem and a family of SDPs as in Definition 18. Let
$\mathcal{T}\subseteq\mathcal{D}$ be such that for all $p,q\in\mathcal{T}$:

  * •

$f(p)\cap f(q)=\emptyset$. That is, a solution to $p$ is not a solution to $q$
and vice-versa.

  * •

Let $A_{j}^{(p)}$ be the constraints of $SDP^{(p)}$ and $A_{j}^{(q)}$ those
from $SDP^{(q)}$ (and define $C^{(p)}$, $C^{(q)}$, $b_{j}^{(p)}$, and
$b_{j}^{(q)}$ in the same manner). Then there exist $\sigma\in S_{n}$, $\pi\in
S_{m}$ s.t. $\sigma^{-1}A^{(p)}_{\pi(j)}\sigma=A^{(q)}_{j}$ (and
$\sigma^{-1}C^{(p)}\sigma=C^{(q)}$). That is, the SDPs are the same up to
permutations of the labels of the constraints and permutations of the
coordinates of $\mathbb{R}^{n}$.

If $y^{(p)}$ is an $\varepsilon$-optimal dual-feasible vector to $SDP^{(p)}$
for some $p\in\mathcal{T}$, then $y^{(p)}$ is at least
$\frac{\log(|\mathcal{T}|)}{\log m}$-dense (i.e., has at least that many non-
zero entries).

###### Proof.

We first observe that, with $SDP^{(p)}$ and $SDP^{(q)}$ as in the lemma, if
$y^{(p)}$ is an $\varepsilon$-optimal dual-feasible vector of $SDP^{(p)}$,
then $y^{(q)}$ defined by

$y^{(q)}_{j}:=y^{(p)}_{\pi(j)}=\pi(y^{(p)})_{j}$

is an $\varepsilon$-optimal dual vector for $SDP^{(q)}$. Here we use the fact
that a permutation of the $n$ coordinates in the primal does not affect the
dual solutions. Since $f(p)\cap f(q)=\emptyset$ we know that $g(y^{(p)})\neq
g(y^{(q)})$ and so $y^{(p)}\neq y^{(q)}$. Since this is true for every $q$ in
$\mathcal{T}$, there should be at least $|\mathcal{T}|$ different vectors
$y^{(q)}=\pi(y^{(p)})$.

A $k$-sparse vector can have $k$ different non-zero entries and hence the
number of possible unique permutations of that vector is at most

$\binom{m}{k}k!=\frac{m!}{(m-k)!}=\prod_{t=m-k+1}^{m}t\leq m^{k}$

so

$\frac{\log|\mathcal{T}|}{\log m}\leq k.$

∎

##### Example.

Consider the $(s,t)$-mincut problem, i.e., the dual of the $(s,t)$-maxflow.
Specifically, consider a simple instance of this problem: the union of two
complete graphs of size $z+1$, where $s$ is in one subgraph and $t$ in the
other. Let the other vertices be labeled by $\\{1,2,\dots,2z\\}$. Every
assignment of the labels over the two halves gives a unique mincut, in terms
of which labels fall on which side of the cut. There is exactly one partition
of the vertices in two sets that cuts no edges (namely the partition consists
of the two complete graphs), and every other partition cuts at least $z$
edges. Hence a $z/2$-approximate cut is a mincut. This means that there are
$\binom{2z}{z}$ problems that require a different output. So for every family
of SDPs that is symmetric under permutation of the vertices and for which a
$z/2$-approximate dual solution gives an $(s,t)$-mincut, the sparsity of a
$z/2$-approximate dual solution is at least666Here $m$ is the number of
constraints, not the number of edges in the graph.

$\frac{\log{\binom{2z}{z}}}{\log m}\geq\frac{z}{\log{m}},$

where we used that $\binom{2z}{z}\geq\frac{2^{2z}}{2\sqrt{z}}$.

###  3.2 General width-bounds are restrictive for certain SDPs

In this section we will show that width-bounds can be restrictive when they do
not consider the specific structure of an SDP.

######  Definition 20.

A function $w(n,m,s,r,R,\varepsilon)$ is called a _general width-bound_ for an
oracle if $w(n,m,s,r,R,\varepsilon)$ is a correct width-bound for that oracle,
for every SDP with parameters $n,m,s,r,R,\varepsilon$. In particular, the
function $w$ may not depend on the structure of the $A_{j}$, $C$, and $b$.

We will show that general width-bounds need to scale with $r^{*}$ (recall that
$r^{*}$ denotes the smallest $\ell_{1}$-norm of an optimal solution to the
dual). We then go on to show that if two SDPs in a class can be combined to
get another element of that class in a natural manner, then, under some mild
conditions, $r^{*}$ will be of the order $n$ and $m$ for some instances of the
class.

We start by showing, for specifically constructed LPs, a lower bound on the
width of any oracle. Although these LPs will not solve any useful problem,
every general width-bound should also apply to these LPs. This gives a lower
bound on general width-bounds.

######  Lemma 21.

For every $n\geq 4$, $m\geq 4$, $s\geq 1$, $R^{*}>0$, $r^{*}>0$, and
$\varepsilon\leq 1/2$ there is an SDP with these parameters for which every
oracle has width at least $\frac{1}{2}r^{*}$.

###### Proof.

We will construct an LP for $n=m=3$. This is enough to prove the lemma since
LPs are a subclass of SDPs and we can increase $n$, $m$, and $s$ by adding
more dimensions and $s$-dense SDP constraints that do not influence the
analysis below. For some $k>0$, consider the following LP

$\displaystyle\max\ \ \ $ $\displaystyle(1,0,0)x$ s.t.
$\displaystyle\begin{bmatrix}1&1&1\\\ 1/k&1&0\\\
-1&0&-1\end{bmatrix}x\leq\begin{bmatrix}R\\\ 0\\\ -R\end{bmatrix}$
$\displaystyle x\geq 0$

where the first row is the primal trace constraint. Notice that
$x_{1}=x_{2}=0$ due to the second constraint. This implies that $\mbox{\rm
OPT}=0$ and, due to the last constraint, that $x_{3}\geq R$. Notice that
$(0,0,R)$ is actually an optimal solution, so $R^{*}=R$.

To calculate $r^{*}$, look at the dual of the LP:

$\displaystyle\min\ \ \ $ $\displaystyle(R,0,-R)y$ s.t.
$\displaystyle\begin{bmatrix}1&1/k&-1\\\ 1&1&0\\\
1&0&-1\end{bmatrix}y\geq\begin{bmatrix}1\\\ 0\\\ 0\end{bmatrix}$
$\displaystyle y\geq 0,$

due to strong duality its optimal value is $0$ as well. This implies
$y_{1}=y_{3}$, so the first constraint becomes $y_{2}\geq k$. This in turn
implies $r^{*}\geq k$, which is actually attained (by $y=(0,k,0)$) so
$r^{*}=k$.

Since the oracle and width-bound should work for every
$x\in\mathbb{R}^{3}_{+}$ and every $\alpha$, they should in particular work
for $x=(R,0,0)$ and $\alpha=0$. In this case the polytope for the oracle
becomes

$\displaystyle\mathcal{P}_{\varepsilon}(x):=\\{y\in\mathbb{R}^{m}:\ $
$\displaystyle y_{1}-y_{3}\leq 0,$ $\displaystyle y_{1}-y_{3}+y_{2}/k\geq
1-\varepsilon,$ $\displaystyle y\geq 0\\}.$

since $b^{T}y=y_{1}-y_{3}$, $c^{T}x=1$, $a_{1}^{T}x=1$, $a_{2}^{T}x=1/k$ and
$a_{3}^{T}x=-1$. This implies that for every
$y\in\mathcal{P}_{\varepsilon}(x)$, we have $y_{2}\geq k/(1-\varepsilon)\geq
k/2\geq r^{*}/2$.

Notice that the term

$\left\lVert\sum^{m}_{j=1}y_{j}A_{j}-C\right\rVert$

in the definition of width for an SDP becomes

$\left\lVert A^{T}y-c\right\rVert_{\infty}$

in the case of an LP. In our case, due to the second constraint in the dual,
we know that

$\left\lVert A^{T}y-c\right\rVert_{\infty}\geq y_{1}+y_{2}\geq\frac{r^{*}}{2}$

for every vector $y$ from $\mathcal{P}_{\varepsilon}(x)$. This shows that any
oracle has width at least $r^{*}/2$ for this LP. ∎

######  Corollary 22.

For every general width-bound $w(n,m,s,r,R,\varepsilon)$, if $n,m\geq 3$,
$s\geq 1$, $r>0$, $R>0$, and $\varepsilon\leq 1/2$, then

$w(n,m,s,r,R,\varepsilon)\geq\frac{r}{2}.$

Note that this bound applies to both our algorithm and the one given by
Brandão and Svore. It turns out that for many natural classes of SDPs,
$r^{*},R^{*}$, $\varepsilon$, $n$ and $m$ can grow linearly for some
instances. In particular, this is the case if SDPs in a class combine in a
natural manner. Take for example two SDP relaxations for the MAXCUT problem on
two graphs $G^{(1)}$ and $G^{(2)}$ (on $n^{(1)}$ and $n^{(2)}$ vertices,
respectively):

$\displaystyle\max\quad$ $\displaystyle\mbox{\rm
Tr}\left(L(G^{(1)})X^{(1)}\right)$ s.t. $\displaystyle\mbox{\rm
Tr}\left(X^{(1)}\right)\leq n^{(1)}$ $\displaystyle\mbox{\rm
Tr}\left(E_{jj}X^{(1)}\right)\leq 1\text{ for }j=1,\dots,n^{(1)}$
$\displaystyle X^{(1)}\succeq 0$

$\displaystyle\max\quad$ $\displaystyle\mbox{\rm
Tr}\left(L(G^{(2)})X^{(2)}\right)$ s.t. $\displaystyle\mbox{\rm
Tr}\left(X^{(2)}\right)\leq n^{(2)}$ $\displaystyle\mbox{\rm
Tr}\left(E_{jj}X^{(2)}\right)\leq 1\text{ for }j=1,\dots,n^{(2)}$
$\displaystyle X^{(2)}\succeq 0$

Where $L(G)$ is the Laplacian of a graph. Note that this is not normalized to
operator norm $\leq 1$, but for simplicity we ignore this here. If we denote
the direct sum of two matrices by $\oplus$, that is

$A\oplus B=\begin{bmatrix}A&0\\\ 0&B\end{bmatrix},$

then, for the disjoint union of the two graphs, we have

$L(G^{(1)}\cup G^{(2)})=L(G^{(1)})\oplus L(G^{(2)}).$

This, plus the fact that the trace distributes over direct sums of matrices,
means that the SDP relaxation for MAXCUT on $G^{(1)}\cup G^{(2)}$ is the same
as a natural combination of the two separate maximizations:

$\displaystyle\max\quad$ $\displaystyle\mbox{\rm
Tr}\left(L(G^{(1)})X^{(1)}\right)+\mbox{\rm Tr}\left(L(G^{(2)})X^{(2)}\right)$
s.t. $\displaystyle\mbox{\rm Tr}\left(X^{(1)}\right)+\mbox{\rm
Tr}\left(X^{(2)}\right)\leq n^{(1)}+n^{(2)}$ $\displaystyle\mbox{\rm
Tr}\left(E_{jj}X^{(1)}\right)\leq 1\text{ for }j=1,\dots,n^{(1)}$
$\displaystyle\mbox{\rm Tr}\left(E_{jj}X^{(2)}\right)\leq 1\text{ for
}j=1,\dots,n^{(2)}$ $\displaystyle X^{(1)},X^{(2)}\succeq 0.$

It is easy to see that the new value of $n$ is $n^{(1)}+n^{(2)}$, the new
value of $m$ is $m^{(1)}+m^{(2)}-1$ and the new value of $R^{*}$ is
$n^{(1)}+n^{(2)}=R^{*(1)}+R^{*(2)}$. Since it is natural for the MAXCUT
relaxation that the additive errors also add, it remains to see what happens
to $r^{*}$, and so, for general width-bounds, what happens to $w$. As we will
see later in this section, under some mild conditions, these kind of
combinations imply that there are MAXCUT-relaxation SDPs for which $r^{*}$
also increases linearly, but this requires a bit more work.

######  Definition 23.

We say a class of SDPs (with associated allowed approximation errors) is
_combinable_ if there is a $k\geq 0$ so that for every two elements in this
class, $(SDP^{(a)},\varepsilon^{(a)})$ and $(SDP^{(b)},\varepsilon^{(b)})$,
there is an instance in the class, $(SDP^{(c)},\varepsilon^{(c)})$, that is a
combination of the two in the following sense:

  * •

$C^{(c)}=C^{(a)}\oplus C^{(b)}$.

  * •

$A^{(c)}_{j}=A^{(a)}_{j}\oplus A^{(b)}_{j}$ and
$b_{j}^{(c)}=b_{j}^{(a)}+b_{j}^{(b)}$ for $j\in[k]$.

  * •

$A^{(c)}_{j}=A^{(a)}_{j}\oplus\mathbf{0}$ and $b^{(c)}_{j}=b^{(a)}_{j}$ for
$j=k+1,\dots,m^{(a})$.

  * •

$A^{(c)}_{m^{(a)}+j-k}=\mathbf{0}\oplus A^{(b)}_{j}$ and
$b^{(c)}_{m^{(a)}+j-k}=b^{(b)}_{j}$ for $j=k+1,\dots,m^{(b)}$.

  * •

$\varepsilon^{(c)}\leq\varepsilon^{(a)}+\varepsilon^{(b)}$.

In other words, some fixed set of constraints are summed pairwise, and the
remaining constraints get added separately.

Note that this is a natural generalization of the combining property of the
MAXCUT relaxations (in that case $k=1$ to account for the trace bound).

######  Theorem 24.

If a class of SDPs is combinable and there is an element $SDP^{(1)}$ for which
every optimal dual solution has the property that

$\sum_{j=k+1}^{m}y_{m}\geq\delta$

for some $\delta>0$, then there is a sequence $(SDP^{(t)})_{t\in\mathbb{N}}$
in the class such that $\frac{R^{*(t)}r^{*(t)}}{\varepsilon^{(t)}}$ increases
linearly in $n^{(t)}$, $m^{(t)}$ and $t$.

###### Proof.

The sequence we will consider is the $t$-fold combination of $SDP^{(1)}$ with
itself. If $SDP^{(1)}$ is  $\displaystyle\max\quad$ $\displaystyle\mbox{\rm
Tr}(CX)$ s.t. $\displaystyle\mbox{\rm Tr}(A_{j}X)\leq b_{j}\quad\text{ for
}j\in[m^{(1)}],$ $\displaystyle X\succeq 0$ $\displaystyle\min\quad$
$\displaystyle\sum_{j=1}^{m^{(1)}}b_{j}y_{j}$ s.t.
$\displaystyle\sum_{j=1}^{m^{(1)}}y_{j}A_{j}-C\succeq 0,$ $\displaystyle y\geq
0$ then $SDP^{(t)}$ is

$\displaystyle\max\quad$ $\displaystyle\ \sum_{i=1}^{t}\mbox{\rm
Tr}\left(CX_{i}\right)$ s.t. $\displaystyle\sum_{i=1}^{t}\mbox{\rm
Tr}(A_{j}X_{i})\leq tb_{j}$ $\displaystyle\text{ for }j\in[k],$
$\displaystyle\mbox{\rm Tr}(A_{j}X_{i})\leq b_{j}$ $\displaystyle\text{ for
}j=k+1,\dots,m^{(1)}\text{ and }i=1,\dots,t$ $\displaystyle X_{i}\succeq 0$
$\displaystyle\text{ for all }i=1,\dots,t$

with dual

$\displaystyle\min\quad$
$\displaystyle\sum_{j=1}^{k}tb_{j}y_{j}+\sum_{i=1}^{t}\sum_{j=k+1}^{m^{(1)}}b_{j}y_{j}^{i}$
s.t.
$\displaystyle\sum_{j=1}^{k}y_{j}A_{j}+\sum_{j=k+1}^{m^{(1)}}y_{j}^{i}A_{j}\succeq
C\text{ for }i=1,\dots,t$ $\displaystyle y,y^{i}\geq 0.$

First, let us consider the value of $\mbox{\rm OPT}^{(t)}$. Let $X^{(1)}$ be
an optimal solution to $SDP^{(1)}$ and for all $i\in[t]$ let $X_{i}=X^{(1)}$.
Since these $X_{i}$ form a feasible solution to $SDP^{(t)}$, this shows that
$\mbox{\rm OPT}^{(t)}\geq t\cdot\mbox{\rm OPT}^{(1)}$. Furthermore, let
$y^{(1)}$ be an optimal dual solution of $SDP^{(1)}$, then
$(y^{(1)}_{1},\dots,y^{(1)}_{k})\oplus\left(y^{(1)}_{k+1},\cdots,y^{(1)}_{m^{(1)}}\right)^{\oplus
t}$ is a feasible dual solution for $SDP^{(t)}$ with objective value
$t\cdot\mbox{\rm OPT}^{(1)}$, so $\mbox{\rm OPT}^{(t)}=t\cdot\mbox{\rm
OPT}^{(1)}$.

Next, let us consider the value of $r^{*(t)}$. Let $\tilde{y}\oplus
y^{1}\oplus\dots\oplus y^{t}$ be an optimal dual solution for $SDP^{(t)}$,
split into the parts of $y$ that correspond to different parts of the
combination. Then $\tilde{y}\oplus y^{i}$ is a feasible dual solution for
$SDP^{(1)}$ and hence $b^{T}(\tilde{y}\oplus y^{i})\geq\mbox{\rm OPT}^{(1)}$.
On the other hand we have

$t\cdot\mbox{\rm OPT}^{(1)}=\mbox{\rm
OPT}^{(t)}=\sum_{i=1}^{t}b^{T}(\tilde{y}\oplus y^{i}),$

this implies that each term in the sum is actually equal to $\mbox{\rm
OPT}^{(1)}$. But if $(\tilde{y}\oplus y^{i})$ is an optimal dual solution of
$SDP^{(1)}$ then $\left\lVert(\tilde{y}\oplus y^{i})\right\rVert_{1}\geq
r^{*(1)}$ by definition and $\left\lVert y^{i}\right\rVert_{1}\geq\delta$. We
conclude that $r^{*(t)}\geq r^{*(1)}-\delta+t\delta$.

Now we know the behavior of $r^{*}$ under combinations, let us look at the
primal to find a similar statement for $R^{*(t)}$. Define a new SDP,
$\widehat{SDP}^{(t)}$, in which all the constraints are summed when combining,
that is, in Definition 23 we take $k=n^{(1)}$, however, contrary to that
definition, we even sum the psd constraints:

$\displaystyle\max\quad$ $\displaystyle\ \sum_{i=1}^{t}\mbox{\rm
Tr}\left(CX_{i}\right)$ s.t. $\displaystyle\sum_{i=1}^{t}\mbox{\rm
Tr}(A_{j}X_{i})\leq tb_{j}\quad\text{ for }j\in[m^{(1)}],$
$\displaystyle\sum_{i=1}^{t}X_{i}\succeq 0.$

This SDP has the same objective function as $SDP^{(t)}$ but a larger feasible
region: every feasible $X_{1},\dots,X_{t}$ for $SDP^{(t)}$ is also feasible
for $\widehat{SDP}^{(t)}$. However, by a change of variables,
$X:=\sum_{i=1}^{t}X_{i}$, it is easy to see that $\widehat{SDP}^{(t)}$ is
simply a scaled version of $SDP^{(1)}$. So, $\widehat{SDP}^{(t)}$ has optimal
value $t\cdot\mbox{\rm OPT}^{(1)}$. Since optimal solutions to
$\widehat{SDP}^{(t)}$ are scaled optimal solutions to $SDP^{(1)}$, we have
$\hat{R}^{*(t)}=t\cdot R^{*(1)}$. Combining the above, it follows that every
optimal solution to $SDP^{(t)}$ is optimal to $\widehat{SDP}^{(t)}$ as well,
and hence has trace at least $t\cdot R^{*(1)}$, so $R^{*(t)}\geq t\cdot
R^{*(1)}$.

We conclude that

$\frac{R^{*(t)}r^{*(t)}}{\varepsilon^{(t)}}\geq\frac{tR^{*(1)}(r^{*(1)}+(t-1)\delta)}{t\varepsilon^{(1)}}=\Omega\left(t\right)$

and $n^{(t)}=tn^{(1)}$, $m^{(t)}=t(m^{(1)}-k)+k$. ∎

This shows that for many natural SDP formulations for combinatorial problems,
such as the MAXCUT relaxation or LPs that have to do with resource management,
$R^{*}r^{*}/\varepsilon$ increases linearly in $n$ and $m$ for some instances.
Hence, using $R^{*}\leq R$ and Lemma 21, $Rw/\varepsilon$ grows at least
linearly when a general width-bound is used.

##  4 Lower bounds on the quantum query complexity

In this section we will show that every LP-solver (and hence every SDP-solver)
that can distinguish two optimal values with high probability needs
$\Omega\left(\sqrt{\max\\{n,m\\}}\left(\min\\{n,m\\}\right)^{3/2}\right)$
quantum queries in the worst case.

For the lower bound on LP-solving we will give a reduction from a composition
of Majority and OR functions.

######  Definition 25.

Given input bits $Z_{ij\ell}\in\\{0,1\\}^{a\times b\times c}$ the problem of
calculating

$\displaystyle MAJ_{a}($ $\displaystyle
OR_{b}(MAJ_{c}(Z_{111},\dots,Z_{11c}),\dots,MAJ_{c}(Z_{1b1},\dots,Z_{1bc})),$
$\displaystyle\dots,$ $\displaystyle
OR_{b}(MAJ_{c}(Z_{a11},\dots,Z_{a1c}),\dots,MAJ_{c}(Z_{ab1},\dots,Z_{abc}))$
$\displaystyle)$

with the promise that

  * •

Each inner $\mbox{\rm MAJ}_{c}$ is a boundary case, in other words
$\sum_{\ell=1}^{c}Z_{ij\ell}\in\\{c/2,c/2+1\\}$ for all $i,j$.

  * •

The outer $\mbox{\rm MAJ}_{a}$ is a boundary case, in other words, if
$\tilde{Z}\in\\{0,1\\}^{a}$ is the bitstring that results from all the OR
calculations, then $|\tilde{Z}|\in\\{a/2,a/2+1\\}$.

is called the promise $\mbox{\rm MAJ}_{a}\mathchar 45\relax\mbox{\rm
OR}_{b}\mathchar 45\relax\mbox{\rm MAJ}_{c}$ problem.

######  Lemma 26.

It takes at least $\Omega(a\sqrt{b}c)$ queries to the input to solve the
promise $\mbox{\rm MAJ}_{a}\mathchar 45\relax\mbox{\rm OR}_{b}\mathchar
45\relax\mbox{\rm MAJ}_{c}$ problem.

###### Proof.

The promise version of $\mbox{\rm MAJ}_{k}$ is known to require $\Omega(k)$
quantum queries. Likewise, it is known that the ORk function requires
$\Omega(\sqrt{k})$ queries. Furthermore, the adversary bound tells us that
query complexity is multiplicative under composition of functions (even
promise functions). For completeness we include a proof of this composition
property as Theorem 61 in Appendix F. The proof is a straightforward
modification of a proof of Lee et al. [LMR+11, Section 5] for the case where
the outer function is a total Boolean function; Kimmel [Kim13] already gave a
similar modification of the proof of Lee et al. for the case where a promise
function is composed with itself. ∎

######  Lemma 27.

Determining the value

$\sum_{i=1}^{a}\max_{j\in[b]}\sum_{\ell=1}^{c}Z_{ij\ell}$

for a $Z$ from the promise $\mbox{\rm MAJ}_{a}\mathchar 45\relax\mbox{\rm
OR}_{b}\mathchar 45\relax\mbox{\rm MAJ}_{c}$ problem up to additive error
$\varepsilon=1/3$, solves the promise $\mbox{\rm MAJ}_{a}\mathchar
45\relax\mbox{\rm OR}_{b}\mathchar 45\relax\mbox{\rm MAJ}_{c}$ problem.

###### Proof.

Notice that due to the first promise,
$\sum_{\ell=1}^{c}Z_{ij\ell}\in\\{c/2,c/2+1\\}$ for all $i\in[a],j\in[b]$.
This implies that

  * •

If the $i$th OR is $0$, then all of its inner MAJ functions are $0$ and hence

$\max_{j\in[b]}\sum_{\ell=1}^{c}Z_{ij\ell}=\frac{c}{2}$

  * •

If the $i$th OR is $1$, then at least one of its inner MAJ functions is $1$
and hence

$\max_{j\in[b]}\sum_{\ell=1}^{c}Z_{ij\ell}=\frac{c}{2}+1$

Now, if we denote the string of outcomes of the OR functions by
$\tilde{Z}\in\\{0,1\\}^{a}$, then

$\sum_{i=1}^{a}\max_{j\in[b]}\sum_{\ell=1}^{c}Z_{ij\ell}=a\frac{c}{2}+|\tilde{Z}|$

Hence determining the left-hand side will determine $|\tilde{Z}|$; this
Hamming weight is either $\frac{a}{2}$ if the full function evaluates to $0$,
or $\frac{a}{2}+1$ if it evaluates to $1$. ∎

######  Lemma 28.

For an input $Z\in\\{0,1\\}^{a\times b\times c}$ there is an LP with $m=c+a$
and $n=c+ab$ for which the optimal value is

$\sum_{i=1}^{a}\max_{j\in[b]}\sum_{\ell=1}^{c}Z_{ij\ell}$

Furthermore, a query to an entry of the input matrix or vector costs at most
$1$ query to $Z$.

###### Proof.

Let $Z^{(i)}$ be the matrix one gets by fixing the first index of $Z$ and
putting the entries in a $c\times b$ matrix, so $Z^{(i)}_{\ell j}=Z_{ij\ell}$.
We define the following LP:

$\displaystyle\mbox{\rm OPT}=\text{max}\ \ \ $
$\displaystyle\sum_{k=1}^{c}w_{k}$ s.t.
$\displaystyle\begin{bmatrix}I&-Z^{1}&\cdots&-Z^{a}\\\ 0&\mathbf{1}^{T}&&\\\
0&&\ddots&\\\ 0&&&\mathbf{1}^{T}\end{bmatrix}\begin{bmatrix}w\\\ v^{(1)}\\\
\vdots\\\ v^{(a)}\\\ \end{bmatrix}\leq\begin{bmatrix}0\\\ \mathbf{1}\\\
\vdots\\\ \mathbf{1}\end{bmatrix}$ $\displaystyle
v^{1},\dots,v^{a}\in\mathbb{R}_{+}^{b},w\in\mathbb{R}_{+}^{c}$

Notice every $Z^{(i)}$ is of size $c\times b$, so that indeed $m=c+a$ and
$n=c+ab$.

For every $i\in[a]$ there is a constraint that says

$\sum_{j=1}^{b}v^{(i)}_{j}\leq 1.$

The constraints involving $w$ say that for every $k\in[c]$

$w_{k}\leq\sum_{i=1}^{a}\sum_{j=1}^{b}v^{(i)}_{j}Z^{(i)}_{kj}=\sum_{i=1}^{a}(Z^{(i)}v^{(i)})_{k}$

where $(Z^{(i)}v^{(i)})_{k}$ is the $k$th entry of the matrix-vector product
$Z^{(i)}v^{(i)}$. Clearly, for an optimal solution these constraints will be
satisfied with equality, since in the objective function $w_{k}$ has a
positive weight. Summing over $k$ on both sides, we get the equality

OPT $\displaystyle=\sum_{k=1}^{c}w_{k}$
$\displaystyle=\sum_{k=1}^{c}\sum_{i=1}^{a}(Z^{(i)}v^{(i)})_{k}$
$\displaystyle=\sum_{i=1}^{a}\sum_{k=1}^{c}(Z^{(i)}v^{(i)})_{k}$
$\displaystyle=\sum_{i=1}^{a}\left\lVert Z^{(i)}v^{(i)}\right\rVert_{1}$

so in the optimum $\left\lVert Z^{(i)}v^{(i)}\right\rVert_{1}$ will be
maximized. Note that we can use the $\ell_{1}$-norm as a shorthand for the sum
over vector elements since all elements are positive. In particular, the value
of $\left\lVert Z^{(i)}v^{(i)}\right\rVert_{1}$ is given by

max $\displaystyle\left\lVert Z^{(i)}v^{(i)}\right\rVert_{1}$ s.t.
$\displaystyle\left\lVert v^{(i)}\right\rVert_{1}\leq 1$ $\displaystyle
v^{(i)}\geq 0$

Now $\left\lVert\smash{Z^{(i)}v^{(i)}}\right\rVert_{1}$ will be maximized by
putting all weight in $v^{(i)}$ on the index that corresponds to the column of
$Z^{(i)}$ that has the highest Hamming weight. In particular in the optimum
$\left\lVert\smash{Z^{(i)}v^{(i)}}\right\rVert_{1}=\max_{j\in[b]}\sum_{\ell=1}^{c}Z^{(i)}_{\ell
j}$. Putting everything together gives:

$\mbox{\rm OPT}=\sum_{i=1}^{a}\left\lVert
Z^{(i)}v^{(i)}\right\rVert_{1}=\sum_{i=1}^{a}\max_{j\in[b]}\sum_{\ell=1}^{c}Z^{(i)}_{\ell
j}=\sum_{i=1}^{a}\max_{j\in[b]}\sum_{\ell=1}^{c}Z_{ij\ell}$

∎

######  Theorem 29.

There is a family of LPs, with $m\leq n$ and two possible integer optimal
values, that require at least $\Omega(\sqrt{n}m^{3/2})$ quantum queries to the
input to distinguish those two values.

###### Proof.

Let $a=c=m/2$ and $b=\frac{n-c}{a}=\frac{2n}{m}-1$, so that $n=c+ab$ and
$m=c+a$. By Lemma 28 there exists an LP with $n=c+ab$ and $m=c+a$ that
calculates

$\sum_{i=1}^{a}\max_{j\in[b]}\sum_{\ell=1}^{c}Z_{ij\ell}$

for an input $Z$ to the promise $\mbox{\rm MAJ}_{a}\mathchar 45\relax\mbox{\rm
OR}_{b}\mathchar 45\relax\mbox{\rm MAJ}_{c}$ problem. By Lemma 27, calculating
this value will solve the promise $\mbox{\rm MAJ}_{a}\mathchar
45\relax\mbox{\rm OR}_{b}\mathchar 45\relax\mbox{\rm MAJ}_{c}$ problem. By
Lemma 26 the promise $\mbox{\rm MAJ}_{a}\mathchar 45\relax\mbox{\rm
OR}_{b}\mathchar 45\relax\mbox{\rm MAJ}_{c}$ problem takes
$\Omega(a\sqrt{b}c)$ quantum queries in the worst case. This implies a lower
bound of

$\Omega\left(m^{2}\sqrt{\frac{n}{m}}\right)=\Omega(m^{3/2}\sqrt{n})$

quantum queries on solving these LPs. ∎

######  Corollary 30.

Distinguishing two optimal values of an LP (and hence also of an SDP) with
additive error $\varepsilon<1/2$ requires

$\Omega\left(\sqrt{\max\\{n,m\\}}\left(\min\\{n,m\\}\right)^{3/2}\right)$

quantum queries to the input matrices in the worst case.

###### Proof.

This follows from Theorem 29 and LP duality. ∎

It is important to note that the parameters $R$ and $r$ from the Arora-Kale
algorithm are not constant in this family of LPs
($R,r=\Theta(\min\\{n,m\\}^{2})$ here), and hence this lower bound does not
contradict the scaling with $\sqrt{mn}$ of the complexity of our SDP-solver or
Brandão and Svore’s. Since we show in the appendix that one can always rewrite
the LP (or SDP) so that $2$ of the parameters $R,r,\varepsilon$ are constant,
the lower bound implies that any algorithm with a sub-linear dependence on $m$
or $n$ has to depend at least polynomially on $Rr/\varepsilon$. For example,
the above family of LPs shows that an algorithm with a $\sqrt{mn}$ dependence
has to have an $(Rr/\varepsilon)^{\kappa}$ factor in its complexity with
$\kappa\geq 1/4$. It remains an open question whether a lower bound of
$\Omega(\sqrt{mn})$ can be proven for a family of LPs/SDPs where
$\varepsilon$, $R$ and $r$ all constant.

##  5 Conclusion

In this paper we gave better algorithms and lower bounds for quantum SDP-
solvers, improving upon recent work of Brandão and Svore [BS16]. Here are a
few directions for future work:

  * •

Better upper bounds. The runtime of our algorithm, like the earlier algorithm
of Brandão and Svore, has better dependence on $m$ and $n$ than the best
classical SDP-solvers, but worse dependence on $s$ and on $Rr/\varepsilon$. It
may be possible to improve the dependence on $s$ to linear and/or the
dependence on $Rr/\varepsilon$ to less than our current 8th power.

  * •

Applications of our algorithm. As mentioned, both our and Brandão-Svore’s
quantum SDP-solvers only improve upon the best classical algorithms for a
specific regime of parameters, namely where $mn\gg Rr/\varepsilon$.
Unfortunately, we don’t know particularly interesting problems in
combinatorial optimization in this regime. As shown in Section 3, many natural
SDP formulations will not fall into this regime. However, it would be
interesting to find useful SDPs for which our algorithm gives a significant
speed-up.

  * •

New algorithms. As in the work by Arora and Kale, it might be more promising
to look at oracles (now quantum) that are designed for specific SDPs. Such
oracles could build on the techniques developed here, or develop totally new
techniques. It might also be possible to speed up other classical SDP solvers,
for example those based on interior-point methods.

  * •

Better lower bounds. Our lower bounds are probably not optimal, particularly
for the case where $m$ and $n$ are not of the same order. The most interesting
case would be to get lower bounds that are simultaneously tight in the
parameters $m$, $n$, $s$, and $Rr/\varepsilon$.

##### Acknowledgments.

We thank Fernando Brandão for sending us several drafts of [BS16] and for
answering our many questions about their algorithms, Stacey Jeffery for
pointing us to [Kim13], and Andris Ambainis and Robin Kothari for useful
discussions and comments. We also thank the anonymous FOCS’17 referees for
helpful comments that improved the presentation.

## References

  * [AG18] Joran van Apeldoorn and András Gilyén.  Improvements in quantum SDP-solving with applications.  arXiv: 1804.05058, 2018.
  * [AGGW18] Joran van Apeldoorn, András Gilyén, Sander Gribling, and Ronald de Wolf.  Convex optimization using quantum oracles.  arXiv: 1809.00643, 2018.
  * [AHK12] S. Arora, E. Hazan, and S. Kale.  The multiplicative weights update method: a meta-algorithm and applications.  Theory of Computing, 8(6):121--164, 2012.
  * [AK16] S. Arora and S. Kale.  A combinatorial, primal-dual approach to semidefinite programs.  Journal of the ACM, 63(2):12, 2016.  Earlier version in STOC’07.
  * [Amb07] A. Ambainis.  Quantum walk algorithm for element distinctness.  SIAM Journal on Computing, 37(1):210--239, 2007.  Earlier version in FOCS’04. quant-ph/0311001.
  * [AS08] N. Alon and J. H. Spencer.  The Probabilistic Method.  Wiley-Interscience, third edition, 2008.
  * [BBHT98] M. Boyer, G. Brassard, P. Høyer, and A. Tapp.  Tight bounds on quantum searching.  Fortschritte der Physik, 46(4--5):493--505, 1998.  Earlier version in Physcomp’96. quant-ph/9605034.
  * [BCC+15] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma.  Simulating Hamiltonian dynamics with a truncated Taylor series.  Physical Review Letters, 114:090502, 2015.  arXiv:1412.4687.
  * [BCK15] D. W. Berry, A. M. Childs, and R. Kothari.  Hamiltonian simulation with nearly optimal dependence on all parameters.  Proceedings of 56th IEEE FOCS, pages 792--809, 2015.  arxiv:1501.01715.
  * [Bel13] A. Belovs.  Applications of the Adversary Method in Quantum Query Algorithms.  PhD thesis, University of Latvia, 2013.  quant-ph/0703264.
  * [BHMT02] G. Brassard, P. Høyer, M. Mosca, and A. Tapp.  Quantum amplitude amplification and estimation.  In Quantum Computation and Quantum Information: A Millennium Volume, volume 305 of AMS Contemporary Mathematics Series, pages 53--74. 2002.  quant-ph/0005055.
  * [BKL+17] Fernando G. S. L. Brandão, Amir Kalev, Tongyang Li, Cedric Yen-Yu Lin, Krysta M. Svore, and Xiaodi Wu.  Quantum SDP solvers: Large speed-ups, optimality, and applications to quantum learning.  arXiv: 1710.02581, 2017.
  * [BS16] F. Brandão and K. Svore.  Quantum speed-ups for semidefinite programming.  To appear in Proceedings of FOCS’17. arXiv:1609.05537v3, 2016.
  * [CCLW18] Shouvanik Chakrabarti, Andrew M. Childs, Tongyang Li, and Xiaodi Wu.  Quantum algorithms and lower bounds for convex optimization.  arXiv: 1809.01731, 2018.
  * [CEMM98] R. Cleve, A. Ekert, C. Macchiavello, and M. Mosca.  Quantum algorithms revisited.  In Proceedings of the Royal Society of London, volume A454, pages 339--354, 1998.  quant-ph/9708016.
  * [CKS15] A. M. Childs, R. Kothari, and R. D. Somma.  Quantum linear systems algorithm with exponentially improved dependence on precision.  arxiv:1511.02306, 7 November 2015.
  * [CS16] A. N. Chowdhury and R. D. Somma.  Quantum algorithms for Gibbs sampling and hitting-time estimation.  arXiv:1603.02940, 9 Mar 2016.
  * [CW12] A. M. Childs and N. Wiebe.  Hamiltonian simulation using linear combinations of unitary operations.  Quantum Information and Computation, 12(11-12):901--924, 2012.  arXiv:1202.5822.
  * [Dan51] G. B. Dantzig.  Maximization of a linear function of variables subject to linear inequalities.  In Activity Analysis of Production and Allocation, Cowles Commission Monograph No. 13, pages 339--347. John Wiley & Sons Inc., New York, N. Y., 1951.
  * [DH96] C. Dürr and P. Høyer.  A quantum algorithm for finding the minimum.  quant-ph/9607014, 18 Jul 1996.
  * [DHHM06] C. Dürr, M. Heiligman, P. Høyer, and M. Mhalla.  Quantum query complexity of some graph problems.  SIAM Journal on Computing, 35(6):1310--1328, 2006.  Earlier version in ICALP’04.
  * [GLS88] M. Grötschel, L. Lovász, and A. Schrijver.  Geometric Algorithms and Combinatorial Optimization.  Springer, 1988.
  * [GR02] L. Grover and T. Rudolph.  Creating superpositions that correspond to efficiently integrable probability distributions.  quant-ph/0208112, 15 Aug 2002.
  * [Gro96] L. K. Grover.  A fast quantum mechanical algorithm for database search.  In Proceedings of 28th ACM STOC, pages 212--219, 1996.  quant-ph/9605043.
  * [GSLW18] András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe.  Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics.  arXiv: 1806.01838, 2018.
  * [GW95] M. X. Goemans and D. P. Williamson.  Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming.  Journal of the ACM, 42(6):1115--1145, 1995.  Earlier version in STOC’94.
  * [HHL09] A. Harrow, A. Hassidim, and S. Lloyd.  Quantum algorithm for solving linear systems of equations.  Physical Review Letters, 103(15):150502, 2009.  arXiv:0811.3171.
  * [HLŠ07] P. Høyer, T. Lee, and R. Špalek.  Negative weights make adversaries stronger.  In Proceedings of 39th ACM STOC, pages 526--535, 2007.  quant-ph/0611054.
  * [Kim13] S. Kimmel.  Quantum adversary (upper) bound.  Chicago Journal of Theoretical Computer Science, 2013.
  * [KKMO07] S. Khot, G. Kindler, E. Mossel, and R. O’Donnell.  Optimal inapproximability results for MAX-CUT and other 2-variable CSPs?  SIAM Journal on Computing, 37(1):319--357, 2007.  Earlier version in FOCS’04.
  * [KP18] Iordanis Kerenidis and Anupam Prakash.  A quantum interior point method for LPs and SDPs.  arXiv: 1808.09266, 2018.
  * [LC16] G. H. Low and I. L. Chuang.  Hamiltonian simulation by qubitization.  arXiv:1610.06546, 20 Oct 2016.
  * [LC17] G. H. Low and I. L. Chuang.  Optimal Hamiltonian simulation by quantum signal processing.  Physical Review Letters, 118:010501, 2017.  arXiv:1606.02685.
  * [LMR+11] T. Lee, R. Mittal, B. Reichardt, R. Špalek, and M. Szegedy.  Quantum query complexity of state conversion.  In Proceedings of 52nd IEEE FOCS, pages 344--353, 2011.  arXiv:1011.3020.
  * [LSW15] Y. T. Lee, A. Sidford, and S. C. Wong.  A faster cutting plane method and its implications for combinatorial and convex optimization.  In Proceedings of 56th IEEE FOCS, pages 1049--1065, 2015.  arXiv:1508.04874.
  * [NC00] M. A. Nielsen and I. L. Chuang.  Quantum Computation and Quantum Information.  Cambridge University Press, 2000.
  * [PW09a] D. Poulin and P. Wocjan.  Preparing ground states of quantum many-body systems on a quantum computer.  Physical Review Letters, 102:130503, 2009.  arXiv:0809.2705.
  * [PW09b] D. Poulin and P. Wocjan.  Sampling from the thermal quantum Gibbs state and evaluating partition functions with a quantum computer.  Physical Review Letters, 103:220502, 2009.  arXiv:0905.2199.
  * [Rei09] B. Reichardt.  Span programs and quantum query complexity: The general adversary bound is nearly tight for every Boolean function.  In Proceedings of 50th IEEE FOCS, pages 544--551, 2009.
  * [Sho97] P. W. Shor.  Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer.  SIAM Journal on Computing, 26(5):1484--1509, 1997.  Earlier version in FOCS’94. quant-ph/9508027.
  * [TRW05] K. Tsuda, G. Rätsch, and M. K. Warmuth.  Matrix exponentiated gradient updates for on-line learning and Bregman projection.  Journal of Machine Learning Research, 6:995--1018, 2005.  Earlier version in NIPS’04.
  * [WK12] M. K. Warmuth and D. Kuzmin.  Online variance minimization.  Machine Learning, 87(1):1--32, 2012.  Earlier version in COLT’06.

##  Appendix A Classical estimation of the expectation value $\mbox{\rm
Tr}\left(A\rho\right)$