Quantum SDP-Solvers: Better upper and lower bounds

TL;DR

This paper improves quantum algorithms for semidefinite programming by enhancing their efficiency and establishing fundamental lower bounds, advancing the understanding of quantum optimization methods.

Contribution

The paper introduces new techniques for quantum algorithms, improves existing SDP-solvers, and establishes lower bounds on their complexity, especially for symmetric problems.

Findings

Enhanced quantum SDP algorithms with better parameter dependence

Development of methods to implement smooth functions of sparse Hamiltonians

Proven lower bounds indicating linear scaling of quantum SDP-solvers in worst-case scenarios

Abstract

Brand\~ao and Svore very 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.

Equations

$$\max$$

$$\text{mbox}\mathrm{Tr}(A_{j}X)\le b_j\text{ for all }j\in [m],$$

$$X⪰0,$$

$$\max$$

$$v_j\in \{+1,-1\}\text{ for all }j\in [n],$$

$$\max$$

$$\text{mbox}\mathrm{Tr}(E_{jj}X)=1\text{ for all }j\in [n],$$

$$X⪰0,$$

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

$$\min$$

$$\sum _{j=1}^m{y}_j{A}_j-C⪰0,$$

$$y\ge 0.$$

$$\displaystyle\mathcal{P}_{\varepsilon}(X^{(1)}):=\{y\in\mathbb{R}^{m}:\$$

$$\displaystyle\mbox{\rm Tr}\left(\Big{(}\sum_{j=1}^{m}y_{j}A_{j}-C\Big{)}X^{(1)}\right)\geq-\varepsilon,$$

$$y\ge 0\}.$$

$$\tilde{\mathcal{O}}\left(nms{\left(\frac{Rr}{\epsilon }\right)}^4+ns{\left(\frac{Rr}{\epsilon }\right)}^7\right).$$

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

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

$$\tilde{\mathcal{O}}(\sqrt{mn}(Rr/\epsilon {)}^5).$$

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

$$\tilde{\mathcal{O}}\left(\sqrt{nm}{s}^2{\left(\frac{Rr}{\epsilon }\right)}^8\right)$$

$$O_I|j,k,\ell ⟩=|j,k,{\text{index}}_{A_j}(k,\ell )⟩.$$

$$O_M|j,k,i,z⟩=|j,k,i,z\oplus (A_j{)}_{ki}⟩.$$

$$QRAM:|i,x,r_1,\dots ,r_K⟩\mapsto |i,r_i,r_1,\dots ,r_{i-1},x,r_{i+1},\dots ,r_K⟩,$$

$${\rho }^{(t+1)}:=\exp \left(-\eta \sum _{\tau =1}^t{M}^{(\tau )}\right)/\text{mbox}\mathrm{Tr}\left(\exp \left(-\eta \sum _{\tau =1}^t{M}^{(\tau )}\right)\right)$$

$$\sum _{t=1}^T\text{mbox}\mathrm{Tr}\left(M^{(t)}{\rho }^{(t)}\right)\le {\lambda }_{\min }\left(\sum _{t=1}^T{M}^{(t)}\right)+\eta \sum _{t=1}^T\text{mbox}\mathrm{Tr}\left((M^{(t)}{)}^2{\rho }^{(t)}\right)+\frac{\ln (n)}{\eta }.$$

$$\frac{1}{T}\sum _{t=1}^T{M}^{(t)}\propto \sum _{j=1}^m{y}_j{A}_j-C$$

$$\displaystyle\mathcal{P}_{\varepsilon}(X):=\{y\in\mathbb{R}^{m}:\$$

$$\displaystyle\mbox{\rm Tr}\left(\Big{(}\sum_{j=1}^{m}y_{j}A_{j}-C\Big{)}X\right)\geq-\varepsilon,$$

$$y\ge 0\}.$$

$$\displaystyle\tilde{\mathcal{P}}(\tilde{a}_{1},\ldots,\tilde{a}_{m},\tilde{c}-(r+1)\theta):=\{y\in\mathbb{R}^{m}:\$$

$$\sum _{j=1}^m{y}_j\le r,$$

$$\sum _{j=1}^m{\tilde{a}}_j{y}_j\ge \tilde{c}-(r+1)\theta$$

$$y\ge 0\}.$$

$${\mathcal{P}}_0(X)\ne \varnothing \Rightarrow {\mathcal{P}}_0(X)\cap \{y\in {\mathbb{R}}^m:\sum _{j=1}^m{y}_j\le r\}\ne \varnothing .$$

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

$$\sum _{j=1}^m{\tilde{a}}_j{y}_j-\tilde{c}\ge \sum _{j=1}^m({\tilde{a}}_j-\text{mbox}\mathrm{Tr}\left(A_j\rho \right))y_j-(\tilde{c}-\text{mbox}\mathrm{Tr}\left(C\rho \right))\ge -\theta {\Vert y\Vert }_1-\theta \ge -(r+1)\theta ,$$

Full text

We also show limits on this approach to quantum SDP-solvers, for instance for combinatorial optimization 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.

1. 1 Introduction

2. 1.1 Semidefinite programs

3. 1.2 Classical solvers for LPs and SDPs

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

5. 1.4 Our results

6. 1.4.1 Improved quantum SDP-solver

7. 1.4.2 Tools that may be of more general interest

8. 1.4.3 Lower bounds

9. 2 An improved quantum SDP-solver

10. 2.1 The Arora-Kale framework for solving SDPs

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

12. 2.2.1 General approach

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

14. 2.2.3 General case – for SDP-solving

15. 2.3 An efficient 2-sparse oracle

16. 2.4 Total runtime

17. 3 Downside of this method: general oracles are restrictive

18. 3.1 Sparse oracles are restrictive

19. 3.2 General width-bounds are restrictive for certain SDPs

20. 4 Lower bounds on the quantum query complexity

21. 5 Conclusion

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

23. B Implementing smooth functions of Hamiltonians

24. B.1 Implementation of smooth functions of Hamiltonians: general results

25. B.2 Applications of smooth functions of Hamiltonians

26. C Generalized minimum-finding algorithm

27. D Sparse matrix summation

28. D.1 A lower bound

29. D.2 An upper bound

30. E Equivalence of $R$, $r$, and $\varepsilon^{-1}$

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 the better-known linear programs (LPs), but (like LPs) they are still efficiently solvable. Thanks to their stronger expressive power, SDPs can sometimes give better approximation algorithms than LPs.

The basic form of an SDP is the following:

[TABLE]

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

[TABLE]

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

[TABLE]

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. The ellipsoid method and interior-point methods can solve LPs and SDPs in polynomial time; they will typically approximate the optimal value to arbitrary precision [GLS81, NN94]. The best known general SDP-solvers [LSW15] approximate the optimal value OPT of such an SDP up to additive error $\varepsilon$, with complexity

[TABLE]

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 the optimization community, first-order methods for semidefinite programming have been considered for instance in [Ren16, Ren19]. 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}$:

[TABLE]

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:

[TABLE]

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

[TABLE]

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 when Brandão and Svore [BS17] 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

[TABLE]

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$) [BS17, 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 [CS17] 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

[TABLE]

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

[TABLE]

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 runtime 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.

Let $f:\mathbb{R}\to\mathbb{C}$ be a function. 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$.555Here a univariate function $f:\mathbb{R}\to\mathbb{C}$ is applied to a Hermitian matrix $H$ in the standard way, by acting on the eigenvalues of $H$: if $H$ has diagonalization $H=U^{-1}DU$, with $D$ the diagonal matrix of $H$’s eigenvalues, then $f(H)$ is the matrix $U^{-1}f(D)U$, where the diagonal matrix $f(D)$ is obtained from $D$ by applying $f$ to its diagonal entries. Roughly speaking, what applying $f(H)$ to $|\phi\rangle$ 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$ 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 41 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 43, only requires that the function $f$ can be nicely approximated locally around each possible eigenvalue of $H$, improving on Theorem 41.

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 to these two functions. For example, they also simplify the analysis of the improved linear-systems solver of Childs et al. [CKS17], 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 non-constant $R$ and $r$ such that a constant-error approximation of OPT computes $F$. Choosing $a$, $b$, and $c$ appropriately, this implies a lower bound of

[TABLE]

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 runtime of quantum SDP-solvers were obtained in [BKL*+*19, vAG19a], 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)$. For the special case of LP-solving, where $s=1$, [vAG19b] further improved the runtime to $\widetilde{\mathcal{O}}\left((\sqrt{m}+\sqrt{n})\left(\frac{Rr}{\epsilon}\right)^{3}\right)$.

In a different algorithmic direction, Kerenidis and Prakash [KP18] recently obtained a quantum interior-point method for solving SDPs and LPs. 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 [GSLW19]. 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.

Recently van Apeldoorn et al. [vAGGdW20] and Chakrabarti et al. [CCLW20] 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), and how to trade off the parameters $R$, $r$, and $1/\varepsilon$ against each other (Appendix E).

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

[TABLE]

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 [math] and the all-one vector by $\mathbf{1}$. 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 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”:

[TABLE]

(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]$:

[TABLE]

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:

[TABLE]

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.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

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 $X\succeq 0$ and $\alpha\in[-R,R]$, 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.

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]).**

Suppose we are 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

[TABLE]

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$:

[TABLE]

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\}$.666Using several transformations of the SDP, from Appendix E and Lemma 2 of [BS17], 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 $\rho=X/\mathrm{Tr}(X)$ where $\mathrm{Tr}(X)\leq R$. 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. Then the following holds:

[TABLE]

Proof.

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

[TABLE]

where in the first inequality we subtracted $\sum_{j=1}^{m}\mbox{\rm Tr}\left(A_{j}\rho\right)y_{j}-\mbox{\rm Tr}\left(C\rho\right)\geq 0$.

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

[TABLE]

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

[TABLE]

For the latter inequality we used $\mbox{\rm Tr}\left(X\right)\leq 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

[TABLE]

Proof.

The inequality can be proven as follows

[TABLE]

where the last step used $\tilde{Z}\geq\frac{2}{3}Z$. ∎

Corollary 8**.**

Let $A$ be such that $\left\lVert A\right\rVert\leq 1$. A multiplicative $\frac{\theta}{9}$-approximation of both $\mbox{\rm Tr}\left(\frac{I}{4}e^{-H}\right)$ and $\mbox{\rm Tr}\left(\frac{I+A/2}{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

[TABLE]

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

[TABLE]

which yields an additive $\theta$-approximation to $\mbox{\rm Tr}\left(Ae^{-H}\right)/\mbox{\rm Tr}\left(e^{-H}\right)$. ∎

It thus suffices to separately approximate both quantities from the corollary. 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:

[TABLE]

(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.)

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. When 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 and some ancilla 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

[TABLE]

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

[TABLE]

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

[TABLE]

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 probability at least $4/5$. The procedure uses

[TABLE]

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

[TABLE]

and thus

[TABLE]

Therefore

[TABLE]

Also by (6) we have

[TABLE]

Therefore using Lemma 9 and setting $\mu=\theta/19$, with probability at least $4/5$ we get a $\tilde{p}$ satisfying

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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

Proof.

First we prepare the state $\sum_{i=1}^{n}|i\rangle/\sqrt{n}$ with ${\mathcal{O}}\left(\log(n)\right)$ one- and two-qubit gates. If $n$ is a power of 2 we do this by applying $\log_{2}(n)$ Hadamard gates on $|0\rangle^{\otimes\log_{2}(n)}$; in the general case it is still possible to prepare the state $\sum_{i=1}^{n}|i\rangle/\sqrt{n}$ with ${\mathcal{O}}\left(\log(n)\right)$ two-qubit gates, for example by preparing the state $\sum_{i=1}^{k}|i\rangle/\sqrt{k}$ for $k=2^{\lceil\log_{2}(n)\rceil}$ and then using (exact) amplitude amplification in order to remove the $i>n$ from the superposition.

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

[TABLE]

The error $\delta_{i}$ is because we only write down 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:

[TABLE]

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

[TABLE]

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

Corollary 12**.**

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

[TABLE]

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

[TABLE]

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 [CS17] for achieving logarithmic dependence on the precision. However, the result of [CS17] 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

[TABLE]

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 44, 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 41. 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 their respective ancilla qubit. As in the linear programming case, we are interested in the probability $p$ of measuring outcome $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.

[TABLE]

Now we show that the above quantity is a good approximation of

[TABLE]

For this we show that $\tilde{V}_{0}^{\dagger}\tilde{V}_{0}\approx(I+A/2)/4$ and $\tilde{W}_{0}\tilde{W}_{0}^{\dagger}\approx e^{-H}$. To see this, first note that for all matrices $B,\tilde{B}$ with $\left\lVert B\right\rVert\leq 1$, we have

[TABLE]

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}$:

[TABLE]

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

[TABLE]

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

[TABLE]

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 41 (see Appendix B), since the map $e^{H/2}$ is already dealt with in Theorem 44. To use Theorem 41, 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 41 constructs. Since we choose the precision parameter to be $\mu/5=\Theta(\mu)$, Theorem 41 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 the cost of our implementation of $e^{-H/2}$ with $\mu/5$ precision: Theorem 44 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

[TABLE]

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 51. 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

To motivate the problem below, 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 an additive $\theta$-approximation to $\mbox{\rm Tr}\left(C\rho\right)$ and $c^{\prime}=\tilde{c}-r\theta-\theta$. We first describe a simplified version of our quantum 2-sparse oracle (Lemma 16) that assumes access to a unitary acting 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$.

At the end of this section we then discuss how to modify the analysis when we are given an oracle that 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}$. We do so since the output of the trace-estimation procedure of the previous section is of this more complicated form.

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

[TABLE]

Our first observation is that the polytope $\tilde{\mathcal{P}}(\tilde{a},c^{\prime})$ is extremely simple: it has only three non-trivial constraints and, if it is non-empty, then it contains a point with at most $2$ non-zero coordinates. The latter follows from general properties of linear programs: any feasible LP with at most $k$ constraints has a $k$-sparse optimal solution (see, e.g., [Sch86, Ch. 7]). Note that non-emptiness of $\tilde{\mathcal{P}}(\tilde{a},c^{\prime})$ is equivalent to the optimal value of

[TABLE]

being at most $r$ (the latter being an LP with only $2$ non-trivial constraints). We will give an alternative, constructive proof that we can obtain a $2$-sparse solution in Lemma 15. This will also illustrate the intuition behind the definition of our $2$-sparse oracle.

Before doing so, let us give a first naive approach to find a $y\in\tilde{\mathcal{P}}(\tilde{a},c^{\prime})$ which will not be sufficiently efficient for our purposes. Using the formulation in Equation (10), we could attempt to find a $y\in\tilde{\mathcal{P}}(\tilde{a},c^{\prime})$ by solving ${\mathcal{O}}(m^{2})$ linear programs with only $2$ variables and $2$ constraints each (these LPs are obtained from (10) by setting all but 2 variables equal to zero) and searching for an optimal solution whose value is at most $r$. Here each linear program is determined by the values $\tilde{a}_{j}$, $b_{j}$ and $c^{\prime}$, and thus we can decide if $\tilde{\mathcal{P}}(\tilde{a},c^{\prime})$ is non-empty using ${\mathcal{O}}(m^{2})$ classical time (given these values).

We use a more sophisticated approach to show that ${\mathcal{O}}(m)$ classical operations (and queries) suffice. Our approach is amenable to a quantum speedup: it also implies that only $\widetilde{\mathcal{O}}\left(\sqrt{m}\right)$ quantum queries suffice. In particular, we now show how to reduce the problem of finding a $y\in\tilde{\mathcal{P}}(\tilde{a},c^{\prime})$ to finding a convex combination of points $(b_{j},\tilde{a}_{j})$ that lies within a certain region of the plane.

First observe that if $\alpha\geq 0$ and $c^{\prime}\leq 0$, then $y=0$ is a solution and our oracle can return it. From now on we will assume that $\alpha<0$ or $c^{\prime}>0$. Then for a feasible point $y$ 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

[TABLE]

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.

Let $y\in\tilde{\mathcal{P}}(\tilde{a},c^{\prime})$, and $N=\left\lVert y\right\rVert_{1}$. Consider $p_{i}=(b_{i},\tilde{a}_{i})$ and $g=(\alpha/N,c^{\prime}/N)$ as before, and write $q=y/N$ such that $\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}$, see Figure 1 for the shape of $\mathcal{G}$. 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})$ with $\left\lVert y\right\rVert_{1}\leq r$ (this scaling can be computed efficiently given $p_{j}$ and $p_{k}$).

Furthermore, regarding the possible (non-)emptiness of $\mathcal{G}$ we know the following by Lemma 6 and Lemma 15:

• •

If $\mathcal{P}_{0}(X)\cap\{y\in\mathbb{R}^{m}\colon\sum_{j}y_{j}\leq r\}$ is non-empty, then some convex combination of two of the $p_{j}$’s lies in $\mathcal{G}$.

• •

If $\mathcal{P}_{4Rr\theta}(X)\cap\{y\in\mathbb{R}^{m}\colon\sum_{j}y_{j}\leq r\}$ is empty, then no convex combination of the $p_{j}$’s lies in $\mathcal{G}$.

We first prove a simplified version of the main result of this section. The analysis below applies if there are $m$ points $p_{j}=(b_{j},\tilde{a}_{j})$, where $j\in[m]$, and we are given a unitary which acts as $|j\rangle|0\rangle|0\rangle\mapsto|j\rangle|\tilde{a}_{j}\rangle|\psi_{j}\rangle$. We later prove the result for more general oracles, which may give superpositions of approximations instead of just the one value $\tilde{a}_{j}$.

Lemma 16** (Simple version).**

There is an algorithm that returns a 2-sparse vector $q$ (with $q\geq 0$ and $\|q\|_{1}=1$) such that $\sum_{j=1}^{m}q_{j}p_{j}\in\mathcal{G}$, 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 problem with high probability) uses ${\mathcal{O}}(\sqrt{m})$ calls to an (exact quantum) 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:

Check if $\alpha\geq 0$ and $c^{\prime}\leq 0$. If so, then return $q=0$. 2. 2.

Check if there is a $p_{i}\in\mathcal{G}$. If so, then return $q=e_{i}$ 3. 3.

Find $p_{j},p_{k}$ so that the line segment $\overline{p_{j}p_{k}}$ goes through $\mathcal{G}$ and return the corresponding $q$. 4. 4.

If the first three steps did not return a vector $q$, then output ‘Fail’.

The main realization is that in step 3 we can search separately for $p_{j}$ and $p_{k}$. We explain this in more detail below, but first we will need a better understanding of the shape of $\mathcal{G}$ (see Figure 1 for illustration). The shape of $\mathcal{G}$ depends on the sign of $\alpha$ and $c^{\prime}$.

• (a)

If $\alpha<0$ and $c^{\prime}<0$. 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 $\alpha<0$ and $c^{\prime}\geq 0$. 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 $\alpha\geq 0$ and $c^{\prime}\leq 0$. This is the case where $y=0$ is a solution, $\mathcal{G}$ is the whole plane and has no corner.

• (d)

If $\alpha\geq 0$ and $c^{\prime}>0$. 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)$.

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.

Now, if we cannot find a single point in $\mathcal{G}$ in step 2, then 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. ∎

We now consider the more general case where we are given access to a unitary which for each $j$ provides a superposition over different values $\tilde{a}_{j}$. We do so because the trace estimation procedure of Corollary 14 provides an oracle of this form.

Lemma 17** (General version).**

Assume that we are given an oracle that acts as

[TABLE]

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}$.

There is a quantum algorithm that uses $\widetilde{\mathcal{O}}\left(\sqrt{m}\right)$ calls to the oracle described above, and the same order of two-qubit gates, and (with high probability) has the following guarantees.

• •

If $\mathcal{P}_{0}(X)\cap\{y\in\mathbb{R}^{m}\colon\sum_{j}y_{j}\leq r\}$ is non-empty, then the algorithm returns a $2$-sparse vector in $\mathcal{P}_{4Rr\theta}(X)\cap\{y\in\mathbb{R}^{m}\colon\sum_{j}y_{j}\leq r\}$.

• •

If $\mathcal{P}_{4Rr\theta}(X)\cap\{y\in\mathbb{R}^{m}\colon\sum_{j}y_{j}\leq r\}$ is empty, then the algorithm correctly concludes that $\mathcal{P}_{0}(X)\cap\{y\in\mathbb{R}^{m}\colon\sum_{j}y_{j}\leq r\}$ is empty.

Proof.

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 quantum algorithm 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 (after this measurement, $\mathcal{G}$ is well-defined). Knowing $\tilde{c}$, we can compute the angles defined by the points $p_{j}^{i}=(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 50) starting with a uniform superposition over the $j$s to find $k,\ell\in[m]$ and points $p_{k}^{i}$ and $p_{\ell}^{i^{\prime}}$ approximately minimizing the respective angles to lines $L_{1},L_{2}$. Here ‘approximately minimizing’ means that there is no $j\in[m]$ such that for all $i^{\prime\prime}$ the angle of $p_{j}^{i^{\prime\prime}}=(b_{j},\tilde{a}_{j}^{i^{\prime\prime}})$ with $L_{1}$ is smaller than that of $p_{k}^{i}$ with $L_{1}$ (and similar for $\ell$ and $L_{2}$). From this point on we can simply consider the model in the simple version of this lemma, since by the analysis above there exists an approximation $\tilde{a}\in\mathbb{R}^{m}$ with $\tilde{a}_{k}=\tilde{a}_{k}^{i}$ and $\tilde{a}_{\ell}=\tilde{a}_{\ell}^{i^{\prime}}$ and where $k$ and $\ell$ are the correct minimizers.

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 the same lemmas guarantee 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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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 18**.**

If, for some specific SDP of the form (1), every $\varepsilon$-optimal dual-feasible vector has at least $\ell$ non-zero elements, then the width $w$ of any $k$-sparse Oracleε/3 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 SDP of interest 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, as we will see below, this is not the case. We will prove that for certain kinds of SDPs, no “useful” dual solution can be very sparse. Intuitively, a dual solution to an SDP is “useful” if it can be turned into a solution of the problem that the SDP is trying to solve. We make this more precise in the definition below.

Definition 19**.**

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)}$:

[TABLE]

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 20**.**

Consider a problem and a family of SDPs as in Definition 19. 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.

• •

The number of constraints $m$ and the primal variable size $n$ are the same for $SDP^{(p)}$ and $SDP^{(q)}$.

• •

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)}$, $b^{(p)}_{\pi(j)}=b^{(q)}_{j}$). 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

[TABLE]

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

[TABLE]

so

[TABLE]

∎

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 least777Here $m$ is the number of constraints, not the number of edges in the graph.

[TABLE]

where we used that $\binom{2z}{z}\geq\frac{2^{2z}}{2\sqrt{z}}$. In particular this holds for the standard linear programming formulation of the $(s,t)$-maxflow/$(s,t)$-mincut problem.

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 21**.**

An algorithm $O$ is called a general oracle if, when provided with an error parameter $\varepsilon$, it correctly implements an Oracleε (as in Algorithm 2) for all inputs. We use $O_{\varepsilon}$ to denote the algorithm provided by $O$ with error parameter $\varepsilon$ fixed. A function $w(n,m,s,r,R,\varepsilon)$ is called a general width-bound for a general oracle if, for every $0<\varepsilon<1/2$, the value $w(n,m,s,r,R,\varepsilon)$ is a correct width-bound (see Definition 4) for $O_{\varepsilon}$ for every SDP with parameters $n,m,s,r$, and $R$. In particular, the function $w$ may not depend on the structure of the input $A_{1},\ldots,A_{m}$, $C$, $b$ or on the value of $\alpha$.

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 22**.**

For every $n\geq 3$, $m\geq 3$, $s\geq 1$, $R^{*}\geq 1$, $r^{*}>0$, there is an SDP with these parameters such that for any $0\leq\varepsilon\leq 1/2$ any Oracleε for this SDP has width at least $r^{*}/2$.

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

[TABLE]

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:

[TABLE]

due to strong duality its optimal value is [math] 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

[TABLE]

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/R)\geq k/2=r^{*}/2$, where the second inequality follows from the assumptions that $\varepsilon\leq 1/2$ and $R\geq 1$.

Notice that the term

[TABLE]

in the definition of width for an SDP becomes

[TABLE]

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

[TABLE]

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

Corollary 23**.**

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

[TABLE]

Proof.

Consider the LP given by Lemma 22 with $r^{*}=r$. Using a general oracle with general width-bound $w$ for this LP implies the corollary. ∎

Note that this bound applies to both our algorithm and the one given by Brandão and Svore. It turns out that for many classes of SDPs it is natural to assume that $m$, $r^{*}$, and $R^{*}$ grow linearly with $n$, and that the “logical” choice of $\varepsilon$ also scales linearly with $n$. We now argue that this is for instance the case when SDPs in a class combine in a natural manner. As an example, consider the MAXCUT problem. For, e.g., $d$-regular graphs the MAXCUT value grows linearly with the number of vertices $n$. Therefore, one is generally interested in finding a constant multiplicative approximation to the optimal value. For $d$-regular graphs this would thus translate to an additive error which scales linearly with the number of vertices. We argue below that for the SDP-relaxation it is also natural to assume that $r$ and $R$ grow linearly in $n$. 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):

[TABLE]

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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)}$. Likewise, it is natural to assume that the desired approximation error for the combined SDP is the sum of the desired errors for the individual SDPs: starting with feasible solutions $X^{(i)}$ that are $\varepsilon^{(i)}$-approximate solutions to the two SDP-relaxations ($i=1,2$), the matrix $X^{(1)}\oplus X^{(2)}$ is an $(\varepsilon^{(1)}+\varepsilon^{(2)})$-approximate solution to the combined SDP. 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 24**.**

We say that a class of SDPs (each with an associated allowed approximation error) 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)}=\varepsilon^{(a)}+\varepsilon^{(b)}$.

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

The motivation for the above definition reflects the following: if $X^{(a)}$ and $X^{(b)}$ are feasible solutions to $SDP^{(a)}$ and $SDP^{(b)}$ that are $\varepsilon^{(a)}$-approximate and $\varepsilon^{(b)}$-approximate solutions respectively, then $X^{(a)}\oplus X^{(b)}$ is an $(\varepsilon^{(a)}+\varepsilon^{(b)})$-approximate solution to $SDP^{(c)}$.

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

Theorem 25**.**

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

[TABLE]

for some $\delta>0$, then there is a sequence $(SDP^{(t)},\varepsilon^{(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$, and $SDP^{(t)}$ is the $t$-fold combination of $SDP^{(1)}$ with itself.

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

[TABLE]

with dual

[TABLE]

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

[TABLE]

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 24 we take $k=n^{(1)}$, however, contrary to that definition, we even sum the psd constraints:

[TABLE]

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

[TABLE]

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 22, $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 26**.**

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

[TABLE]

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 27**.**

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; Kimmel [Kim13, Lemma A.3 (Lemma 6 in the arXiv version)] showed that this even holds for promise functions. ∎

Lemma 28**.**

Determining the value

[TABLE]

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 [math], then all of its inner MAJ functions are [math] and hence

[TABLE]

• •

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

[TABLE]

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

[TABLE]

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

Lemma 29**.**

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

∎

Theorem 30**.**

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 29 there exists an LP with $n=c+ab$ and $m=c+a$ that calculates

[TABLE]

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 28, 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 27 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

[TABLE]

quantum queries on solving these LPs. ∎

Corollary 31**.**

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

[TABLE]

quantum queries to the input matrices in the worst case.

Proof.

This follows from Theorem 30 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 [BS17]. 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$. In subsequent work (see the introduction), these dependencies have been improved, but there is room for further improvement.

• •

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.

An earlier version of this paper appeared in the proceedings of the FOCS’17 conference [vAGGdW17].

We thank Fernando Brandão for sending us several drafts of [BS17] and for answering our many questions about their algorithms. We thank Stacey Jeffery for pointing us to [Kim13], and Andris Ambainis and Robin Kothari for useful discussions and comments. Thanks to Monique Laurent for finding a bug in an earlier version, which we fixed. We also thank the anonymous referees of FOCS’17 and Quantum for very helpful comments that improved the presentation.

JvA and SG are supported by the Netherlands Organization for Scientific Research (NWO), grant number 617.001.351. AG and RdW are supported by ERC Consolidator Grant 615307-QPROGRESS. RdW is also partially supported by NWO through Gravitation-grant Quantum Software Consortium - 024.003.037, and through QuantERA project QuantAlgo 680-91-034.

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

To provide contrast to Section 2.2, here we describe a classical procedure to efficiently estimate $\mathrm{Tr}(A\rho)$ where $A$ is a Hermitian matrix such that $\left\lVert A\right\rVert\leq 1$, and $\rho=\exp(-H)/\mbox{\rm Tr}\left(\exp(-H)\right)$ for some Hermitian matrix $H$. The results in this section can be seen as a generalization of [AK16, Section 7]. The key observation is that if we are given a Hermitian matrix $B\succeq 0$, and if we take a random vector $u=(u_{1},\ldots,u_{n})$ where $u_{i}\in\{\pm 1\}$ is uniformly distributed, then, using $\mathbb{E}[u_{i}]=0$, $\mathbb{E}[u_{i}^{2}]=1$, we have

[TABLE]

We now show that $u^{T}\sqrt{B}A\sqrt{B}u$ is highly concentrated around its mean by Chebyshev’s inequality.

Lemma 32**.**

Given a Hermitian matrix $A$, with $||A||\leq 1$, a psd matrix $B$, and a parameter $0<\theta\leq 1$. With probability $1-1/16$, the average of $k={\mathcal{O}}\left(1/\theta^{2}\right)$ independent samples from the distribution $u^{T}\sqrt{B}A\sqrt{B}u$ is at most $\theta\mbox{\rm Tr}\left(B\right)$ away from $\mbox{\rm Tr}\left(AB\right)$. Here $u=(u_{i})$ and each $u_{i}\in\{\pm 1\}$ is i.i.d. uniformly distributed.

Proof.

We let $F_{k}$ be the random variable $\frac{1}{k}\sum_{i=1}^{k}(u^{(i)})^{T}\sqrt{B}A\sqrt{B}u^{(i)}$, where each of the vectors $u^{(i)}\in\{\pm 1\}^{n}$ is sampled from the distribution described above. By the above it is clear that $\mathbb{E}[F_{k}]=\mathrm{Tr}(AB)$.

We will use Chebyshev’s inequality which, in our setting, states that for every $t>0$

[TABLE]

here $\sigma_{k}^{2}$ is the variance of $F_{k}$. We will now upper bound the variance of $F_{k}$. First note that $\mathrm{var}(F_{k})=\frac{1}{k}\mathrm{var}(u^{T}\sqrt{B}A\sqrt{B}u)$. It therefore suffices to upper bound the variance $\sigma^{2}$ of $u^{T}\sqrt{B}A\sqrt{B}u$. We first write

[TABLE]

We then calculate $(*)$ using $\mathbb{E}[u_{i}]=0$, $\mathbb{E}[u_{i}^{2}]=1$, and the independence of the $u_{i}$’s:

[TABLE]

Therefore, using Cauchy-Schwarz, we have

[TABLE]

where on the last line we use $\left\lVert A\right\rVert\leq 1$ and $\mbox{\rm Tr}\left(AY\right)\leq\left\lVert A\right\rVert\mbox{\rm Tr}\left(Y\right)$ for any $Y\succeq 0$, in particular for $BA^{2}B$ and $B^{2}$.

It follows that $\sigma_{k}^{2}\leq 2\mbox{\rm Tr}\left(B\right)^{2}/k$. Chebyshev’s inequality (12) therefore shows that for $k=\lceil 32/\theta^{2}\rceil$ and $t=4$,

[TABLE]

∎

A simple computation shows that the success probability in the above lemma can be boosted to $1-\delta$ by picking the median of ${\mathcal{O}}(\log(1/\delta))$ repetitions. To show this, let $K=\lceil\log(1/\delta)\rceil$ and for each $i\in[K]$ let $(F_{k})_{i}$ be the average of $k$ samples of $u^{T}\sqrt{B}A\sqrt{B}u$. Let $z_{K}$ denote the median of those $K$ numbers. We have

[TABLE]

We upper bound $(*)$:

[TABLE]

Analogously, one can show that \mathrm{Pr}\big{(}z_{K}\leq-\theta\mbox{\rm Tr}\left(B\right)+\mbox{\rm Tr}\left(AB\right)\big{)}\leq\frac{1}{2}\delta. Hence

[TABLE]

This proves the following lemma:

Lemma 33**.**

Given a Hermitian matrix $A$, with $||A||\leq 1$, a psd matrix $B$, and parameters $0<\delta\leq 1/2$ and $0<\theta\leq 1$. Using $k={\mathcal{O}}\left(\log(\frac{1}{\delta})/\theta^{2}\right)$ samples from the distribution $u^{T}\sqrt{B}A\sqrt{B}u$, one can find an estimate of $\mbox{\rm Tr}\left(AB\right)$ that, with probability $1-\delta$, has additive error at most $\theta\mbox{\rm Tr}\left(B\right)$. Here $u=(u_{i})$ and the $u_{i}\in\{\pm 1\}$ are i.i.d. uniformly distributed.

Looking back at Meta-Algorithm 1, we would like to apply the above lemma to $B=\exp(-H)$.888For ease of notation, we write $H$ for $\eta H$. Since it is expensive to compute the exponent of a matrix, it is of interest to consider samples from $v^{T}Av$, where $v$ is an approximation of $\sqrt{B}u$.

Say $\left\lVert\sqrt{B}u-v\right\rVert\leq\kappa$ and, as always, $\left\lVert A\right\rVert\leq 1$. Then

[TABLE]

Now observe that we are interested in

[TABLE]

Suppose an upper bound $K$ on $\left\lVert H\right\rVert$ is known, then we can consider $H^{\prime}=H-KI$ which satisfies $H^{\prime}\preceq 0$. It follows that $\left\lVert\exp(-H^{\prime})\right\rVert\geq 1$ and, with $B=\exp(-H^{\prime})$, therefore $\left\lVert\smash{\sqrt{B}}\right\rVert\leq\left\lVert B\right\rVert\leq\mbox{\rm Tr}\left(B\right)$. Hence, taking $\kappa\leq\min\{\theta/\left\lVert u\right\rVert,\left\lVert\smash{\sqrt{B}}u\right\rVert\}=\theta/\left\lVert u\right\rVert$,999Here we assume that $\theta\leq 1$. Then, since $\lambda_{\min}(B)\geq 1$, we trivially have $\theta/\left\lVert u\right\rVert\leq 1/\sqrt{n}\leq\sqrt{n}\leq\left\lVert\smash{\sqrt{B}}u\right\rVert$. we find

[TABLE]

This shows that the additional error incurred by sampling from $v^{T}Av$ is proportional to $\theta\mbox{\rm Tr}\left(B\right)$. Finally, a $\kappa$-approximation of $\smash{\sqrt{B}}u$, with $\kappa=\theta/\left\lVert u\right\rVert$ can be obtained by using the truncated Taylor series of $\exp(-H^{\prime}/2)$ of degree $p=\max\{2e\left\lVert H^{\prime}\right\rVert,\log\left(\frac{\sqrt{n}}{\theta}\right)\}$:

[TABLE]

Lemma 34**.**

Given a Hermitian $s$-sparse matrix $A$, with $||A||\leq 1$, a psd matrix $B=\exp(-H)$ with $H\preceq 0$, for a $d$-sparse $H$, and parameters $0<\delta\leq 1/2$ and $0<\theta\leq 1$. With probability $1-\delta$, using $k={\mathcal{O}}\left(\log(\frac{1}{\delta})/\theta^{2}\right)$ samples from the distribution $v^{T}Av$, one can find an estimate that is at most $\theta\mbox{\rm Tr}\left(B\right)$ away from $\mbox{\rm Tr}\left(AB\right)$. Here

[TABLE]

where $p={\mathcal{O}}(\max\{\left\lVert H\right\rVert,\log\left(\frac{\sqrt{n}}{\theta}\right)\})$, and $u=(u_{j})$ where the $u_{j}\in\{\pm 1\}$ are i.i.d. uniformly distributed.

Lemma 35**.**

Given $m$ Hermitian $s$-sparse $n\times n$ matrices $A_{1}=I,A_{2},\ldots,A_{m}$, with $\left\lVert A_{j}\right\rVert\leq 1$ for all $j$, a Hermitian $d$-sparse $n\times n$ matrix $H$ with $\left\lVert H\right\rVert\leq K$, and parameters $0<\delta\leq 1/2$ and $0<\theta\leq 1$. With probability $1-\delta$, we can compute $\theta$-approximations $a_{1},\ldots,a_{m}$ of $\mbox{\rm Tr}\left(A_{1}B\right)/\mbox{\rm Tr}\left(B\right),\ldots,\mbox{\rm Tr}\left(A_{m}B\right)/\mbox{\rm Tr}\left(B\right)$ where $B=\exp(-H)$, using

[TABLE]

queries to the entries of $A_{1},\ldots,A_{m}$, $H$ and arithmetic operations.

Proof.

As observed above, for every matrix $A$,

[TABLE]

Lemma 34 states that for $B^{\prime}=\exp(-H+KI)$ and $A\in\{A_{1},\ldots,A_{m}\}$ using $k={\mathcal{O}}\left(\log(\frac{m}{\delta})/\theta^{2}\right)$ samples from the distribution $v^{T}Av$, one can find an estimate that is at most $\theta\mbox{\rm Tr}\left(B\right)$ away from $\mbox{\rm Tr}\left(AB\right)$ with probability $1-\delta/m$. Here

[TABLE]

where $p={\mathcal{O}}(\max\{K,\log\left(\frac{\sqrt{n}}{\theta}\right)\})$, and $u=(u_{j})$ with the $u_{j}\in\{\pm 1\}$ i.i.d. uniformly distributed.

Observe that the $k$ samples from $v^{T}Av$ are really obtained from $k$ samples of vectors $u=(u_{j})$ combined with some post-processing, namely obtaining $v=\sum_{i=0}^{p}\frac{((H-KI)/2)^{i}}{i!}u$ and two more sparse matrix vector products.

We can therefore obtain $k$ samples from each of $v^{T}A_{1}v,\ldots,v^{T}A_{m}v$ by once calculating $k$ vectors $v=\sum_{i=0}^{p}\frac{((H-KI)/2)^{i}}{i!}u$, and then, for each of the $m$ matrices $A_{j}$ computing the $k$ products $v^{T}A_{j}v$. The $k$ vectors $v$ can be constructed using

[TABLE]

queries to the entries of $H$ and arithmetic operations. The $mk$ matrix vector products can be computed using

[TABLE]

arithmetic operations and queries to the entries of $A_{1},\ldots,A_{m}$ and $H$. This leads to total complexity

[TABLE]

for computing $k$ samples from each of $v^{T}A_{1}v,\ldots,v^{T}A_{m}v$.

The results of Lemma 34 say that for each $j$, using those $k$ samples of $v^{T}A_{j}v$ we can construct a $\theta\mbox{\rm Tr}\left(B^{\prime}\right)/4$-approximation $a_{j}^{\prime}$ of $\mbox{\rm Tr}\left(A_{j}B^{\prime}\right)$, with probability $1-\delta/(2m)$. Therefore, by a union bound, with probability $1-\delta/2$ we can construct $\theta\mbox{\rm Tr}\left(B^{\prime}\right)/4$-approximations $a_{1}^{\prime},\ldots,a_{m}^{\prime}$ of $\mbox{\rm Tr}\left(A_{1}B^{\prime}\right),\ldots,\mbox{\rm Tr}\left(A_{m}B^{\prime}\right)$. Therefore, for each $j$, with probability at least $1-\delta$, by Lemma 7 we have that $a_{j}=a_{j}^{\prime}/a_{1}^{\prime}$ is a $\theta$-approximation of $\mbox{\rm Tr}\left(A_{j}B^{\prime}\right)/\mbox{\rm Tr}\left(B^{\prime}\right)$, and hence it is a $\theta$-approximation of $\mbox{\rm Tr}\left(A_{j}B\right)/\mbox{\rm Tr}\left(B\right)$. ∎

Appendix B Implementing smooth functions of Hamiltonians

In this appendix we show how to efficiently implement smooth functions of a given Hamiltonian. First we explain what we mean by a function of a Hamiltonian $H\in\mathbb{C}^{n\times n}$, i.e., a Hermitian matrix. Since Hermitian matrices are diagonalizable using a unitary matrix, we can write $H=U^{\dagger}\text{diag}(\lambda)U$, where $\lambda\in\mathbb{R}^{n}$ is the vector of eigenvalues. Then for a function $f:\mathbb{R}\rightarrow\mathbb{C}$ we define $f(H):=U^{\dagger}\text{diag}(f(\lambda))U$ with a slight abuse of notation, where we apply $f$ to the eigenvalues in $\lambda$ one-by-one. Note that if we approximate $f$ by $\tilde{f}$, then $\left\lVert\tilde{f}(H)-f(H)\right\rVert=\left\lVert\text{diag}(\tilde{f}(\lambda))-\text{diag}(f(\lambda))\right\rVert$. Suppose $D\subseteq\mathbb{R}$ is such that $\lambda\in D^{n}$, then we can upper bound this norm by the maximum of $|\tilde{f}(x)-f(x)|$ over $x\in D$. Finally we note that $D=[-\left\lVert H\right\rVert,\left\lVert H\right\rVert]$ is always a valid choice.

The main idea of the method presented below, is to implement a map $\tilde{f}(H)$, where $\tilde{f}$ is a good (finite) Fourier approximation of $f$ for all $x\in[-\left\lVert H\right\rVert,\left\lVert H\right\rVert]$. The novelty in our approach is that we construct a Fourier approximation based on some polynomial approximation. In the special case, when $f$ is analytic and $\left\lVert H\right\rVert$ is less than the radius of convergence of the Taylor series, we can obtain good polynomial approximation functions simply by truncating the Taylor series, with logarithmic dependence on the precision parameter. Finally we implement the Fourier series using Hamiltonian simulation and the Linear Combination of Unitaries (LCU) trick [CW12, BCC*+*15, BCK15].

This approach was already used in several earlier papers, particularly in [CKS17, CS17]. There the main technical difficulty was to obtain a good truncated Fourier series. This is a non-trivial task, since on top of the approximation error, one needs to optimize two other parameters of the Fourier approximation that determine the complexity of implementation, namely:

$\bullet$ the largest time parameter $t$ that appears in some Fourier term $e^{-itH}$, and

$\bullet$ the total weight of the coefficients, by which we mean the $1$-norm of the vector of coefficients.

Earlier works used clever integral approximations and involved calculus to construct a good Fourier approximation for a specific function $f$. We are not aware of a general result.

In contrast, our Theorem 41 and Corollary 43 avoids the usage of any integration. It obtains a low-weight Fourier approximation function using the Taylor series. The described method is completely general, and has the nice property that the maximal time parameter $t$ depends logarithmically on the desired approximation precision. Since it uses the Taylor series, it is easy to apply to a wide range of smooth functions.

The circuit we describe for the implementation of the linear operator $f(H):\mathbb{C}^{n}\rightarrow\mathbb{C}^{n}$ is going to depend on the specific function $f$, but not on $H$; the $H$-dependence is only coming from Hamiltonian simulation. Since the circuit for a specific $f$ can be constructed in advance, we do not need to worry about the (polynomial) cost of constructing the circuit, making the analysis simpler. When we describe gate complexity, we count the number of two-qubit gates needed for a quantum circuit implementation, just as in Section 2.

Since this appendix presents stand-alone results, here we will deviate slightly from the notation used throughout the rest of the paper, to conform to the standard notation used in the literature (for example, $\varepsilon$, $r$, $\theta$ and $a$ have a different meaning in this appendix). For simplicity we also assume, that the Hamiltonian $H$ acts on $\mathbb{C}^{n}$, where $n$ is a power of 2. Whenever we write $\log(\text{formula})$ in some complexity statement we actually mean $\log_{2}(2+\text{formula})$ in order to avoid incorrect near-[math] or even negative expressions in complexity bounds that would appear for small values of the formula.

Hamiltonian simulation.

We implement each term in a Fourier series using a Hamiltonian simulation algorithm, and combine the terms using the LCU Lemma. Specifically we use [BCK15], but in fact our techniques would work with any kind of Hamiltonian simulation algorithm.101010For example there is a more recent method for Hamiltonian simulation [LC19, LC17] that could possibly improve on some of the log factors we get from [BCK15], but one could even consider completely different input models allowing different simulation methods. The following definition describes what we mean by controlled Hamiltonian simulation.

Definition 36**.**

Let $M=2^{J}$ for some $J\in\mathbb{N}$, $\gamma\in\mathbb{R}$ and $\epsilon\geq 0$. We say that the unitary

[TABLE]

implements controlled $(M,\gamma)$-simulation of the Hamiltonian $H$, where $|m\rangle$ denotes a (signed) bitstring $|b_{J}b_{J-1}\ldots b_{0}\rangle$ such that $m=-b_{J}2^{J}+\sum_{j=0}^{J-1}b_{j}2^{j}$. The unitary $\tilde{W}$ implements controlled $(M,\gamma,\varepsilon)$-simulation of the Hamiltonian $H$, if

[TABLE]

Note that in this definition we assume that both positive and negative powers of $e^{iH}$ are simulated. This is necessary for our Fourier series, but sometimes we use only positive powers, e.g., for phase estimation; in that case we can simply ignore the negative powers.

The following lemma is inspired by the techniques of [CKS17]. It calculates the cost of such controlled Hamiltonian simulation in terms of queries to the input oracles (4)-(5) as described in Section 2.

Lemma 37**.**

Let $H\in\mathbb{C}^{n\times n}$ be a $d$-sparse Hamiltonian. Suppose we know an upper bound $K\in\mathbb{R_{+}}$ on the norm of $H$, i.e., $\left\lVert H\right\rVert\leq K$, and let $\tau:=M\gamma K$. If $\varepsilon>0$ and $\gamma=\Omega(1/(Kd))$, then a controlled $(M,\gamma,\varepsilon)$-simulation of $H$ can be implemented using ${\mathcal{O}}(\tau d\log(\tau/\varepsilon)/\log\log(\tau/\varepsilon))$ queries and ${\mathcal{O}}\left(\tau d\log(\tau/\varepsilon)/\log\log(\tau/\varepsilon)\left[\log(n)+\log^{\frac{5}{2}}(\tau/\varepsilon)\right]\right)$ gates.

Proof.

We use the results of [BCK15, Lemma 9-10], which tell us that a $d$-sparse Hamiltonian $H$ can be simulated for time $t$ with $\varepsilon$ precision in the operator norm using

[TABLE]

queries and gate complexity

[TABLE]

Now we use a standard trick to remove log factors from the implementation cost, and write the given unitary $W$ as the product of some increasingly precisely implemented controlled Hamiltonian simulation unitaries. For $b\in\{0,1\}$ let us introduce the projector $|b\rangle\!\langle b|_{j}:=I_{2^{j}}\otimes|b\rangle\!\langle b|\otimes I_{2^{J-j}}$, where $J=\log(M)$. Observe that

[TABLE]

The $j$-th operator $e^{\pm i2^{j}\gamma H}$ in the product (15) can be implemented with $2^{j-J-1}\epsilon$ precision using ${\mathcal{O}}\left(2^{j}\gamma Kd\log\left(\frac{2^{j}\gamma K}{\varepsilon 2^{j-J-1}}\right)/\log\log\left(\frac{2^{j}\gamma K}{\varepsilon 2^{j-J-1}}\right)\right)={\mathcal{O}}(2^{j}\gamma Kd\log(\tau/\epsilon)/\log\log(\tau/\epsilon))$ queries by (13) and using ${\mathcal{O}}(2^{j}d\log(\tau/\epsilon)/\log\log(\tau/\epsilon)[\log(n)+\log^{\frac{5}{2}}(\tau/\epsilon)])$ gates by (14). Let us denote by $\tilde{W}$ the concatenation of all these controlled Hamiltonian simulation unitaries. Adding up the costs we see that our implementation of $\tilde{W}$ uses ${\mathcal{O}}(\tau d\log(\tau/\varepsilon)/\log\log(\tau/\varepsilon))$ queries and has gate complexity ${\mathcal{O}}\left(\tau d\log(\tau/\varepsilon)/\log\log(\tau/\varepsilon)\left[\log(n)+\log^{\frac{5}{2}}(\tau/\varepsilon)\right]\right)$. Using the triangle inequality repeatedly, it is easy to see that $\left\lVert W\!-\!\tilde{W}\right\rVert\leq\sum_{j=0}^{J}2^{j-J-1}\varepsilon\leq\varepsilon$. ∎

B.1 Implementation of smooth functions of Hamiltonians: general results

The first lemma we prove provides the basis for our approach. It shows how to turn a polynomial approximation of a function $f$ on the interval $[-1,1]$ into a nice Fourier series in an efficient way, while not increasing the weight of coefficients. This is useful, because we can implement a function given by a Fourier series using the LCU Lemma, but only after scaling it down with the weight of the coefficients.

Lemma 38**.**

Let $\delta,\varepsilon\in\!(0,1)$ and $f:\mathbb{R}\rightarrow\mathbb{C}$ s.t. $\left|f(x)\!-\!\sum_{k=0}^{K}a_{k}x^{k}\right|\leq\varepsilon/4$ for all $x\in\![-1+\delta,1-\delta]$. Then $\exists\,c\in\mathbb{C}^{2M+1}$ such that

[TABLE]

for all $x\in\![-1+\delta,1-\delta]$, where $M=\max\left(2\left\lceil\ln\left(\frac{4\left\lVert a\right\rVert_{1}}{\varepsilon}\right)\frac{1}{\delta}\right\rceil,0\right)$ and $\left\lVert c\right\rVert_{1}\leq\left\lVert a\right\rVert_{1}$. Moreover $c$ can be efficiently calculated on a classical computer in time $\text{poly}(K,M,\log(1/\varepsilon))$.

Proof.

Let us introduce the notation $\left\lVert f\right\rVert_{\infty}=\sup\{|f(x)|:x\in\![-1+\delta,1-\delta]\}$. First we consider the case when $\left\lVert a\right\rVert_{1}<\varepsilon/2$. Then $\left\lVert f\right\rVert_{\infty}\leq\left\lVert f(x)-\sum_{k=0}^{K}a_{k}x^{k}\right\rVert_{\infty}+\left\lVert\sum_{k=0}^{K}a_{k}x^{k}\right\rVert_{\infty}<\varepsilon/4+\varepsilon/2<\varepsilon$. So in this case the statement holds with $M=0$ and $c=0$, i.e., even with an empty sum.

From now on we assume $\left\lVert a\right\rVert_{1}\geq\varepsilon/2$. We are going to build up our approximation gradually. Our first approximate function $\tilde{f}_{1}(x):=\sum_{k=0}^{K}a_{k}x^{k}$ satisfies $\left\lVert f-\tilde{f}_{1}\right\rVert_{\infty}\!\!\leq\varepsilon/4$ by assumption. In order to construct a Fourier series, we will work towards a linear combination of sines. To that end, note that $\forall x\in[-1,1]$: $\tilde{f}_{1}(x)\!=\!\sum_{k=0}^{K}a_{k}\left(\frac{\arcsin\left(\sin(x\pi/2)\right)}{\pi/2}\right)^{\!k}\!$. Let $b^{(k)}$ denote the series of coefficients such that $\left(\frac{\arcsin(y)}{\pi/2}\right)^{k}=\sum_{\ell=0}^{\infty}b_{\ell}^{(k)}y^{\ell}$ for all $y\in[-1,1]$. For $k=1$ the coefficients are just $\frac{2}{\pi}$ times the coefficients of the Taylor series of $\arcsin$ so we know that $b^{(1)}_{2\ell}=0$ while $b^{(1)}_{2\ell+1}=\binom{2\ell}{\ell}\frac{2^{-2\ell}}{2\ell+1}\frac{2}{\pi}$. Since $\left(\frac{\arcsin(y)}{\pi/2}\right)^{\!k+1}\!\!=\left(\frac{\arcsin(y)}{\pi/2}\right)^{\!k}\left(\sum_{\ell=0}^{\infty}b_{\ell}^{(1)}y^{\ell}\right)$, we obtain the formula $b^{(k+1)}_{\ell}=\sum_{\ell^{\prime}=0}^{\ell}b^{(k)}_{\ell^{\prime}}b^{(1)}_{\ell-\ell^{\prime}}$, so one can recursively calculate each $b^{(k)}$. As $b^{(1)}\geq 0$ one can use the above identity inductively to show that $b^{(k)}\geq 0$. Therefore $\left\lVert b^{(k)}\right\rVert_{1}=\sum_{\ell=0}^{\infty}b_{\ell}^{(k)}1^{\ell}=\left(\frac{\arcsin(1)}{\pi/2}\right)^{k}=1$. Using the above definitions and observations we can rewrite

[TABLE]

To obtain the second approximation function, we want to truncate the summation over $\ell$ at $L=\ln\left(\frac{4\left\lVert a\right\rVert_{1}}{\varepsilon}\right)\frac{1}{\delta^{2}}$ in the above formula. We first estimate the tail of the sum. We are going to use that for all $\delta\in[0,1]$: $\sin((1-\delta)\pi/2)\leq 1-\delta^{2}$. For all $k\in\mathbb{N}$ and $x\in\![-1+\delta,1-\delta]$ we have:

[TABLE]

Thus we have $\left\lVert\tilde{f}_{1}-\tilde{f}_{2}\right\rVert_{\infty}\leq\varepsilon/4$ for

[TABLE]

To obtain our third approximation function, we will approximate $\sin^{\ell}(x\pi/2)$. First observe that

[TABLE]

which, as we will show (for $M^{\prime}$ much larger than $\sqrt{\ell}$) is very well approximated by

[TABLE]

Truncating the summation in (16) based on this approximation reduces the maximal time evolution parameter (i.e., the maximal value of the parameter $t$ in the $\exp(izt)$ terms) quadratically. To make this approximation precise, we use Chernoff’s inequality [AS08, A.1.7] for the binomial distribution, or more precisely its corollary for sums of binomial coefficients, stating

[TABLE]

Let $M^{\prime}=\left\lceil\ln\left(\frac{4\left\lVert a\right\rVert_{1}}{\varepsilon}\right)\frac{1}{\delta}\right\rceil$ and suppose $\ell\leq L$, then this bound implies that

[TABLE]

where for the last inequality we use the assumption $\varepsilon\leq 2\left\lVert a\right\rVert_{1}$. By combining (16) and (17) we get that for all $\ell\leq L$

[TABLE]

Substituting $z=x\pi/2$ into this bound we can see that $\left\lVert\tilde{f}_{2}-\tilde{f}_{3}\right\rVert_{\infty}\leq\varepsilon/2$, for

[TABLE]

using $\sum_{k=0}^{K}|a_{k}|\sum_{\ell=0}^{\lfloor L\rfloor}\left|b^{(k)}_{\ell}\right|\leq\sum_{k=0}^{K}|a_{k}|=\left\lVert a\right\rVert_{1}$. Therefore we can conclude that $\tilde{f}_{3}$ is an $\varepsilon$-approximation to $f$:

[TABLE]

Observe that in (18) the largest value of $|m-\ell|$ in the exponent is upper bounded by $2M^{\prime}=M$. So by rearranging the terms in $\tilde{f}_{3}$ we can write $\tilde{f}_{3}(x)=\sum_{m=-M}^{M}c_{m}e^{\frac{i\pi m}{2}x}$. Now let us fix a value $k$ in the first summation of (18). Observe that after taking the absolute value of each term, the last two summations still yield a value $\leq 1$, since $\left\lVert b^{(k)}\right\rVert_{1}=1$ and $\sum_{m=0}^{\ell}\binom{\ell}{m}=2^{\ell}$. It follows that $\left\lVert c\right\rVert_{1}\leq\left\lVert a\right\rVert_{1}$. From the construction of the proof, it is easy to see that (an $\varepsilon$-approximation of) $c$ can be calculated in time $\text{poly}(K,M,\log(1/\varepsilon))$. ∎

Now we present the Linear Combination of Unitaries (LCU) Lemma [CW12, BCC*+*15, BCK15], which we will use for combining the Fourier terms in our quantum circuit. Since we intend to use LCU for implementing non-unitary operations, we describe a version without the final amplitude amplification step. We provide a short proof for completeness.

Lemma 39** (LCU Lemma [CW12, BCC*+*15, BCK15]).**

Let $U_{1},U_{2},\ldots,U_{m}$ be unitaries on a Hilbert space $\mathcal{H}$, and $L=\sum_{i=1}^{m}a_{i}U_{i}$, where $a\in\mathbb{R}_{+}^{m}\setminus\{0\}$. Let $V=\sum_{i=1}^{m}|i\rangle\!\langle i|\otimes U_{i}$ and $A\in\mathbb{C}^{m\times m}$ be a unitary such that $A|0\rangle=\sum_{i=1}^{m}\sqrt{\frac{a_{i}}{\left\lVert a\right\rVert_{1}}}|i\rangle$. Then $\frac{L}{\left\lVert a\right\rVert_{1}}=\left(\langle 0|\otimes I\right)\left(A^{\dagger}\otimes I\right)V\left(A\otimes I\right)\left(|0\rangle\otimes I\right)$, i.e., for every $|\psi\rangle\in\mathcal{H}$ we have $\left(A^{\dagger}\otimes I\right)V\left(A\otimes I\right)|0\rangle|\psi\rangle=|0\rangle\frac{L}{\left\lVert a\right\rVert_{1}}|\psi\rangle+|\Phi^{\perp}\rangle$, where the vector $|\Phi^{\perp}\rangle$ satisfies $\left(|0\rangle\!\langle 0|\otimes I\right)|\Phi^{\perp}\rangle=0$.

Proof.

[TABLE]

∎

The next result summarizes how to efficiently implement a Fourier series of a Hamiltonian.

Lemma 40**.**

Suppose $f(x)=\sum_{m=-M}^{M-1}c_{m}e^{im\gamma x}$, for a given $c\in\mathbb{C}^{2M}\setminus\{0\}$. We can construct a unitary $\tilde{U}$ which implements the operator $\frac{f(H)}{\left\lVert c\right\rVert_{1}}=\sum_{m=-M}^{M-1}\frac{c_{m}}{\left\lVert c\right\rVert_{1}}e^{im\gamma H}$ with $\varepsilon$ precision, i.e., such that

[TABLE]

using ${\mathcal{O}}(M(\log(M)+1))$ two-qubit gates and a single use of a circuit implementing controlled $(M,\gamma,\varepsilon)$-simulation of $H$.

Proof.

This is a direct corollary of Lemma 39. To work out the details, note that we can always extend $c$ with some [math] values, so we can assume without loss of generality that $M$ is a power of $2$. This is useful, because then we can represent each $m\in[-M,M-1]$ as a $(J+1)$-bit signed integer for $J=\log(M)$.

The implementation of the operator $A$ in Lemma 39 does not need any queries and it can be constructed exactly using ${\mathcal{O}}(M(\log(M)+1))$ two-qubit gates, e.g., by the techniques of [GR02]. We sketch the basic idea which is based on induction. For $J=1$ the operator $A$ is just a two-qubit unitary. Suppose we proved the claim for bitstrings of length $J$ and want to prove the claim for length $J+1$. Let $a\in\mathbb{R}_{+}^{2^{J+1}}$ be such that $\left\lVert a\right\rVert_{1}=1$ and define $\tilde{a}\in\mathbb{R}_{+}^{2^{J}}$ such that $\tilde{a}_{b}=a_{b,0}+a_{b,1}$ for all bitstrings $b\in\{0,1\}^{J}$. Then we have a circuit $\tilde{A}$ that uses ${\mathcal{O}}(2^{J}(J+1))$ gates and satisfies $\sqrt{\tilde{a}_{b}}=\langle b|\tilde{A}|0\ldots 0\rangle$ for all $b\in\{0,1\}^{J}$. We can add an extra $|0\rangle$-qubit and implement a controlled rotation gate $R_{b}$ on it for each $b\in\{0,1\}^{J}$. Let $R_{b}$ have rotation angle $\arccos\left(\sqrt{a_{b,0}/\tilde{a}_{b}}\right)$ and be controlled by $b$. It is easy to see that the new unitary $A$ satisfies $\sqrt{a_{b^{\prime}}}=\langle b^{\prime}|A|0\ldots 0\rangle$ for each $b^{\prime}\in\{0,1\}^{J}$. Each $R_{b}$ can be implemented using ${\mathcal{O}}(J)$ two-qubit gates and ancilla qubits, justifying the gate complexity and concluding the induction.

What remains is to implement the operator $V\kern-0.56905pt=\sum_{m=-M}^{M-1}|m\rangle\!\langle m|\otimes\frac{c_{m}}{|c_{m}|}e^{im\gamma H}$ from Lemma 39. We implement $V=PW$ in two steps, where $P=\sum_{m=-M}^{M-1}|m\rangle\!\langle m|\otimes\frac{c_{m}}{|c_{m}|}I$. This $P$ can be implemented exactly using ${\mathcal{O}}(M(\log(M)+1))$ gates simply by building a controlled gate that adds the right phase for each individual bitstring. Since the bitstring on which we want to do a controlled operation has length $\log(M)+1$, each controlled operation can be constructed using ${\mathcal{O}}(\log(M)+1)$ gates and ancilla qubits resulting in the claimed gate complexity. We use a circuit implementing controlled $(M,\gamma,\varepsilon)$-simulation of $H$, denoted by $\tilde{W}$, which is an $\varepsilon$-approximation of $W$ by definition.

Finally $\tilde{U}:=(A^{\dagger}\otimes I)P\tilde{W}\left(A\otimes I\right)$. This yields an $\varepsilon$-precise implementation, since

[TABLE]

∎

Now we can state the main result of this appendix, which tells us how to efficiently turn a function (provided with its Taylor series) of a Hamiltonian $H$, into a quantum circuit by using controlled Hamiltonian simulation.

In the following theorem we assume that the eigenvalues of $H$ lie in a radius-$r$ ball around $x_{0}$. The main idea is that if even $r+\delta$ is less than the radius of convergence of the Taylor series, then we can obtain an $\varepsilon$-approximation of $f$ by truncating the series at logarithmically high powers. $B$ will be an upper bound on the absolute value of the function within the $r+\delta$ ball around $x_{0}$, in particular $\left\lVert f(H)/B\right\rVert\leq 1$. Therefore we can implement $f(H)/B$ as a block of some larger unitary. It turns out that apart from the norm and sparsity of $H$ and precision parameters, the complexity depends on the ratio of $\delta$ and $r$.

Theorem 41** (Implementing a smooth function of a Hamiltonian).**

Let $x_{0}\in\mathbb{R}$ and $r>0$ be such that $f(x_{0}+x)=\sum_{\ell=0}^{\infty}a_{\ell}x^{\ell}$ for all $x\in\![-r,r]$. Suppose $B>0$ and $\delta\in(0,r]$ are such that $\sum_{\ell=0}^{\infty}(r+\delta)^{\ell}|a_{\ell}|\leq B$. If $\left\lVert H-x_{0}I\right\rVert\leq r$ and $\varepsilon\in\!\left(0,\frac{1}{2}\right]$, then we can implement a unitary $\tilde{U}$ such that $\left\lVert(\langle 0|\otimes I)\tilde{U}(|0\rangle\otimes I)-\frac{f(H)}{B}\right\rVert\leq\varepsilon$, using ${\mathcal{O}}\left(r/\delta\log\left(r/(\delta\varepsilon)\right)\log\left(1/\varepsilon\right)\right)$ gates and a single use of a circuit for controlled $\left({\mathcal{O}}(r\log(1/\varepsilon)/\delta),{\mathcal{O}}(1/r),\varepsilon/2\right)$-simulation of $H$.

Suppose we are given $K$ such that $\left\lVert H\right\rVert\leq K$ and $r={\mathcal{O}}(K)$. If, furthermore, $H$ is $d$-sparse and is accessed via oracles (4)-(5), then the whole circuit can be implemented using

[TABLE]

Proof.

The basic idea is to combine Lemma 38 and Lemma 40 and apply them to a transformed version of the function. First we define $\delta^{\prime}:=\delta/(r+\delta)$, which is at most $1/2$ by assumption. Then, for all $\ell\in\mathbb{N}$ let $b_{\ell}:=a_{\ell}(r+\delta)^{\ell}$ and define the function $g:[-1+\delta^{\prime},1-\delta^{\prime}]\rightarrow\mathbb{R}$ by $g(y):=\sum_{\ell=0}^{\infty}b_{\ell}y^{\ell}$ so that

[TABLE]

Now we set $L:=\left\lceil\frac{1}{\delta^{\prime}}\log\left(\frac{8}{\varepsilon}\right)\right\rceil$. Then for all $y\in\![-1+\delta^{\prime},1-\delta^{\prime}]$

[TABLE]

We would now like to obtain a Fourier-approximation of $g$ for all $y\in[-1+\delta^{\prime},1-\delta^{\prime}]$, with precision $\varepsilon^{\prime}=\frac{\varepsilon B}{2}$. Let $b^{\prime}:=(b_{0},b_{1},\ldots,b_{L-1})$ and observe that $\left\lVert b^{\prime}\right\rVert_{1}\leq\left\lVert b\right\rVert_{1}\leq B$. We apply Lemma 38 to the function $g$, using the polynomial approximation corresponding to the truncation to the first $L$ terms, i.e., using the coefficients in $b^{\prime}$. Then we obtain a Fourier $\varepsilon^{\prime}$-approximation $\tilde{g}(y):=\sum_{m=-M}^{M}\tilde{c}_{m}e^{\frac{i\pi m}{2}y}$ of $g$, with

[TABLE]

such that the vector of coefficients $\tilde{c}\in\mathbb{C}^{2M+1}$ satisfies $\left\lVert\tilde{c}\right\rVert_{1}\leq\left\lVert b^{\prime}\right\rVert_{1}\leq\left\lVert b\right\rVert_{1}\leq B$. Let

[TABLE]

by (19) we see that $\tilde{f}$ is an $\varepsilon^{\prime}$-precise Fourier approximation of $f$ on the interval $[x_{0}-r,x_{0}+r]$. To transform this Fourier series to its final form, we note that $\tilde{f}(z)=\sum_{m=-M}^{M}\tilde{c}_{m}e^{\frac{i\pi m}{2(r+\delta)}(z-x_{0})}$, so by defining $c_{m}:=\tilde{c}_{m}e^{-\frac{i\pi m}{2(r+\delta)}x_{0}}$ we get a Fourier series in $z$, while preserving $\left\lVert c\right\rVert_{1}=\left\lVert\tilde{c}\right\rVert_{1}\leq B$.

In the trivial case, when $c=0$, we choose a unitary $\tilde{U}$, such that it maps the $|0\rangle$ ancilla state to $|1\rangle$, then clearly $(\langle 0|\otimes I)\tilde{U}(|0\rangle\otimes I)=0=\tilde{f}(H)$. Clearly such a $\tilde{U}$ can be implemented using ${\mathcal{O}}(1)$ gates and [math] queries. Otherwise we can apply Lemma 40 to this modified Fourier series to construct a unitary circuit $\tilde{V}$ implementing an $\frac{\varepsilon}{2}$-approximation of $\tilde{f}(H)/\left\lVert c\right\rVert_{1}$. We can further scale down the amplitude of the $|0\rangle$-part of the output by a factor of $\left\lVert c\right\rVert_{1}/B\leq 1$, to obtain an approximation of $\tilde{f}(H)/B$ as follows. We simply add an additional ancilla qubit initialized to $|0\rangle$ on which we act with the one-qubit unitary

[TABLE]

Finally we define $\tilde{U}:=Rot\otimes\tilde{V}$, and define $|0\rangle|0\rangle$ as the new success indicator, where the first qubit is the new ancilla. We show that $\tilde{U}$ implements $f(H)/B$ with $\varepsilon$ precision: (if $c=0$, let us use the definition $\tilde{f}(H)/\left\lVert c\right\rVert_{1}:=0$)

[TABLE]

Lemma 40 uses ${\mathcal{O}}\left(M\log(M+1)\right)={\mathcal{O}}\left(r/\delta\log(1/\varepsilon)\log\left(r/(\delta\varepsilon)\right)\right)$ gates and a single use of a controlled $\left(M,\gamma=\pi/(2r+2\delta),\varepsilon/2\right)$-simulation of $H$. If $\left\lVert H\right\rVert={\mathcal{O}}(K)$, we can use Lemma 37 to conclude ${\mathcal{O}}\left(M\gamma Kd\log\left(\frac{1}{\varepsilon}\right)\log\left(\frac{M\gamma K}{\varepsilon}\log\left(\frac{1}{\varepsilon}\right)\right)\right)={\mathcal{O}}\left(\frac{Kd}{\delta}\log\left(\frac{K}{\delta\varepsilon}\right)\log\left(\frac{1}{\varepsilon}\right)\right)$ query and

[TABLE]

gate complexity. Finally note that the polynomial cost of calculating $c$ that is required by Lemma 38 does not affect the query complexity or the circuit size, it only affects the description of the circuit. ∎

Remark 42**.**

Note that in the above theorem we can relax the criterion $\left\lVert H-x_{0}I\right\rVert\leq r$. Suppose we have an orthogonal projector $\Pi$, which projects to eigenvectors with eigenvalues in $[x_{0}-r,x_{0}+r]$, i.e., $[H,\Pi]=0$ and $\left\lVert\Pi\left(H-x_{0}I\right)\Pi\right\rVert\leq r$. Then the circuit $\tilde{U}$ constructed in Theorem 41 satisfies

[TABLE]

The following corollary shows how to implement functions piecewise in small “patches” using Remark 42. The main idea is to first estimate the eigenvalues of $H$ up to $\theta$ precision, and then implement the function using the Taylor series centered around a point close to the eigenvalue.

This approach has multiple advantages. First, the function may not have a Taylor series that is convergent over the whole domain of possible eigenvalues of $H$. Even if there is such a series, it can have very poor convergence properties, making $B$ large and therefore requiring a lot of amplitude amplification. Nevertheless, for small enough neighborhoods the Taylor series always converges quickly, overcoming this difficulty.

Corollary 43**.**

Suppose $(x_{\ell})\in\mathbb{R}^{L}$ and $r,\theta\in\mathbb{R_{+}}$ are such that the spectrum of $H$ lies in the domain $\bigcup_{\ell=1}^{L}[x_{\ell}-(r-2\theta),x_{\ell}+(r+2\theta)]$.111111This way, even if we make $\theta$ error during the estimation of an eigenvalue $\lambda_{i}$, the closest $x_{\ell}$ will still contain $\lambda_{i}$ in its radius-$r$ neighborhood. Suppose there exist coefficients $a_{k}^{(\ell)}\in\mathbb{R}$ such that for all $\ell\in[L]$ and $x\in\![-r,r]$ we have $f(x_{\ell}+x)=\sum_{k=0}^{\infty}a_{k}^{(\ell)}x^{k}$, and $\sum_{k=0}^{\infty}(r+\delta)^{k}|a_{k}^{(\ell)}|\leq B$ for some fixed $\delta\in[0,r]$ and $B>0$. If $\left\lVert H\right\rVert\leq K$ and $\varepsilon\in\!\left(0,\frac{1}{2}\right]$, then we can implement a unitary $\tilde{U}$ such that

[TABLE]

using ${\mathcal{O}}\left(Lr/\delta\log\left(r/(\delta\varepsilon)\right)\log\left(1/\varepsilon\right)+\log(K/\theta)\log\log(K/(\theta\varepsilon))\right)$ gates, and with ${\mathcal{O}}(\log(1/\varepsilon))$ uses of an $({\mathcal{O}}(1/\theta),\pi/K,\Omega(\varepsilon^{2}/\log(1/\varepsilon)))$-simulation of $H$ and a single use of a circuit for controlled $\left({\mathcal{O}}(r\log(1/\varepsilon)/\delta),{\mathcal{O}}(1/r),\varepsilon/2\right)$-simulation of $H$. If $r={\mathcal{O}}(K)$, $\theta\leq r/4$, $\theta=\Omega(\delta)$, $\left\lVert H\right\rVert\leq K$, $H$ is $d$-sparse and is accessed via oracles (4)-(5), then the circuit can be implemented using

[TABLE]

Sketch of the proof..

We start by performing phase estimation on $e^{iH}$ with $\approx\theta$ resolution in phase. We boost the success probability by taking the median outcome of ${\mathcal{O}}(\log(1/\varepsilon))$ parallel repetitions, so that we get a worse-than-$\theta$ estimation with probability at most ${\mathcal{O}}(\varepsilon^{2})$. This way the boosted phase estimation circuit is ${\mathcal{O}}(\varepsilon)$-close in operator norm to an “idealized” phase estimation unitary that never makes approximation error greater than $\theta$ (for more details on this type of argument, see the proof of Lemma 51). Phase estimation uses controlled $({\mathcal{O}}(K/\theta),\pi/K,\Omega(\varepsilon^{2}/\log(1/\varepsilon)))$-simulation of $H$ and a Fourier transform on ${\mathcal{O}}(\log(K/\theta))$-bit numbers which can be implemented using ${\mathcal{O}}(\log(K/\theta)\log\log(K/(\theta\varepsilon)))$ gates. The probability-boosting uses ${\mathcal{O}}(\log(1/\varepsilon))$ repetitions. Controlled on the phase estimate $\tilde{\lambda}$ that we obtained, we implement a $1/B$-scaled version of the corresponding function “patch” $\left.f(x)\right|_{[x_{\ell}-r,x_{\ell}+r]}$ centered around $\arg\min|x_{\ell}-\tilde{\lambda}|$ using Theorem 41 and Remark 42. The additional gate complexities of the “patches” add up to ${\mathcal{O}}\left(Lr/\delta\log\left(r/(\delta\varepsilon)\right)\log\left(1/\varepsilon\right)\right)$, but since each “patch” uses the same controlled $\left({\mathcal{O}}(r\log(1/\varepsilon)/\delta),{\mathcal{O}}(1/r),\varepsilon/2\right)$-simulation of $H$, we only need to implement that once. Finally we uncompute phase estimation.121212Note that phase estimation on some eigenvector of $e^{iH}$ can produce a superposition of different estimates of the phase. If some intervals $[x_{\ell}-r,x_{\ell}+r]$ overlap for $\ell$ and $\ell^{\prime}$, those estimates could lead to different implementations of $f(H)$ (one based on the coefficients $a_{k}^{(\ell)}$ and one based on $a_{k}^{(\ell^{\prime})}$). However, this causes no difficulty; since we used the same normalization $1/B$ for all implementations, both implementations lead to essentially the same state after postselecting on the $|0\rangle$ ancilla state. For the final complexity, note that we can assume without loss of generality that $L(r-2\theta)={\mathcal{O}}(K)$, since otherwise we can just remove some redundant intervals from the domain. Hence $Lr={\mathcal{O}}(K)$ and Lemma 37 implies the stated complexities. ∎

This corollary is essentially as general and efficient as we can hope for. Let $D$ denote the domain of possible eigenvalues of $H$. If we want to implement a reasonably smooth function $f$, then it probably satisfies the following: there is some $r=\Omega(1)$, such that for each $x\in D$, the Taylor series in the radius-$r$ neighborhood of $x$ converges quickly, more precisely the Taylor coefficients $a_{k}^{(x)}$ for the $x$-centered series satisfy $\sum_{k=0}^{\infty}|a_{k}^{(x)}|r^{k}={\mathcal{O}}(\left\lVert f\right\rVert_{\infty})$, where we define $\left\lVert f\right\rVert_{\infty}:=\sup_{x\in D}|f(x)|$. If this is the case, covering $D$ with radius-${\mathcal{O}}(r)$ intervals, choosing $\theta=\Theta(r)$ and $\delta=\Theta(r)$, Corollary 43 provides an $\widetilde{\mathcal{O}}\left(\left\lVert H\right\rVert d\right)$ query and gate complexity implementation of $f(H)/B$, where $B={\mathcal{O}}(\left\lVert f\right\rVert_{\infty})$. The value of $B$ is optimal up to constant factors, since $f(H)/B$ must have norm at most $1$. Also the $\left\lVert H\right\rVert d$ factor in the complexity is very reasonable, and we achieve the logarithmic error dependence which is the primary motivation of the related techniques. An application along the lines of this discussion can be found in Lemma 47.

Also note that in the above corollary we added up the gate complexities of the different “patches.” Since these gates prepare the Fourier coefficients of the function corresponding to the different Taylor series at different points, one could use this structure to implement all coefficients with a single circuit. This can potentially result in much smaller circuit sizes, which could be beneficial when the input model allows more efficient Hamiltonian simulation (which then would no longer be the bottleneck in the complexity).

B.2 Applications of smooth functions of Hamiltonians

In this subsection we use the input model for the $d$-sparse matrix $H$ as described at the start of Section 2. We calculate the implementation cost in terms of queries to the input oracles (4)-(5), but it is easy to convert the results to more general statements as in the previous subsection.

The following theorem shows how to efficiently implement the function $e^{-H}$ for some $H\succeq I$. We use this result in the proof of Lemma 13 to estimate expectation values of the quantum state $\rho=e^{-H}/\mbox{\rm Tr}\left(e^{-H}\right)$ (for the application we ensure that $H\succeq I$ by adding some multiple of $I$).

Theorem 44**.**

Suppose that $I\preceq H$ and we are given $K\in\mathbb{R}_{+}$ such that $\left\lVert H\right\rVert\leq 2K$. If $\varepsilon\in(0,1/3)$, then we can implement a unitary $\tilde{U}$ such that $\left\lVert(\langle 0|\otimes I)\tilde{U}(|0\rangle\otimes I)-e^{-H}\right\rVert\leq\varepsilon$ using

[TABLE]

Proof.

In order to use Theorem 41 we set $x_{0}:=K+1/2$ so that $\left\lVert H-x_{0}I\right\rVert\leq K=:r$, and use the function

[TABLE]

We choose $\delta:=1/2$ so that $e^{-x_{0}}\sum_{\ell=0}^{\infty}\frac{(r+\delta)^{\ell}}{\ell!}=e^{-x_{0}}\sum_{\ell=0}^{\infty}\frac{x_{0}^{\ell}}{\ell!}=1$, therefore we set $B:=1$. Theorem 41 tells us that we can implement a unitary $\tilde{U}$, such that $\tilde{f}(H):=(\langle 0|\otimes I)\tilde{U}(|0\rangle\otimes I)$ is an $\varepsilon$-approximation of $f(H)/B=e^{-H}$, using

[TABLE]

∎

To conclude this appendix, we now sketch the proofs of a few interesting consequences of the earlier results in this appendix. These will, however, not be used in the body of the paper.

First, we show how to use the above subroutine together with amplitude amplification to prepare a Gibbs state with cost depending logarithmically on the precision parameter, as shown by the following lemma. To our knowledge this is the first Gibbs sampler that achieves logarithmic dependence on the precision parameter without assuming access to the entries of $\sqrt{H}$ as in [CS17]. This can mean a significant reduction in complexity; for more details see the introduction of Section 2.2.3.

Lemma 45**.**

We can probabilistically prepare a purified Gibbs state $|\tilde{\gamma}\rangle_{AB}$ such that with high probability $\left\lVert\mbox{\rm Tr}_{B}\left(|\tilde{\gamma}\rangle\!\langle\tilde{\gamma}|_{AB}\right)\kern-0.56905pt-\kern-0.56905pte^{-H}\!/\mbox{\rm Tr}\left(e^{-H}\right)\kern-0.85358pt\right\rVert_{1}\leq\varepsilon$ holds, using an expected cost $\widetilde{\mathcal{O}}\left(\!\sqrt{n/\mbox{\rm Tr}\left(e^{-H}\right)}\right)$ times the complexity of Theorem 44. If we are given a number $z\leq\mbox{\rm Tr}\left(e^{-H}\right)$, then we can also prepare $|\tilde{\gamma}\rangle_{AB}$ in a unitary fashion with cost $\widetilde{\mathcal{O}}\left(\sqrt{n/z}\right)$ times the complexity of Theorem 44.

Sketch of proof.

First we show how to prepare a purified sub-normalized Gibbs state. Then we use the exponential search algorithm of Boyer et al. [BBHT98] (with exponentially decreasing guesses for the norm $a$ of the subnormalized Gibbs state, and hence exponentially increasing number of amplitude amplification steps) to postselect on this sub-normalized state in a similar fashion as in Algorithm 3. There is a possible caveat here: if we postselect on a state with norm $a$, then it gets rescaled by $1/a$ and its preparation error is rescaled by $1/a$ as well. Therefore during the rounds of the search algorithm we always increase the precision of implementation to compensate for the increased (error) amplification. Since the success of postselection in a round is upper bounded by the square of ${\mathcal{O}}\left(a\cdot\#\{\text{amplification steps in the round}\}\right)$, the probability for the postselection to succeed in any of the early rounds is small.

Now we describe how to prepare a purified sub-normalized Gibbs state. We use the decomposition $H\!=\!\sum_{j=1}^{n}E_{j}|\phi_{j}\rangle\!\langle\phi_{j}|$, where $\{|\phi_{j}\rangle:j\in[n]\}$ is an orthonormal eigenbasis of $H$. Due to the invariance of maximally entangled states under transformations of the form $W\otimes W^{*}$ for unitary $W$, we have

[TABLE]

Suppose we can implement a unitary $U$ such that $(\langle 0|\otimes I)U(|0\rangle\otimes I)=e^{-H/2}$. If we apply $U$ to the A-register of the state (20), then we get a state $|\gamma\rangle$ such that $\mbox{\rm Tr}_{B}\left((\langle 0|\otimes I)|\gamma\rangle\!\langle\gamma|(|0\rangle\otimes I)\right)=e^{-H}/n$.

If we implement $e^{-H/2}$ with sufficient precision using Theorem 44 in the exponential search algorithm, then after ${\mathcal{O}}\left(\sqrt{n/\mbox{\rm Tr}\left(e^{-H}\right)}\right)$ rounds of amplitude amplification, with high probability we obtain a Gibbs state using the claimed expected runtime.

If we also know a lower bound $z\leq\mbox{\rm Tr}\left(e^{-H}\right)$, then we have an upper bound on the expected runtime, therefore we can turn the procedure into a unitary circuit using standard techniques. ∎

We can also recover the Gibbs sampler of Chowdhury and Somma [CS17]: if we apply our Corollary 43 to the function $e^{-x^{2}}$ assuming access to $\sqrt{H}$ for some psd matrix $H$, then we get a Gibbs sampler for the state $e^{-H}/\mbox{\rm Tr}(e^{-H})$, similar to [CS17]. The advantage of the presented approach is that it avoids the usage of involved integral transformations, and can be presented without writing down a single integral sign, also due to our general results the proof is significantly shorter. Before we prove the precise statement in Lemma 47, we need some preparation for the application of Corollary 43:

Lemma 46**.**

For all $k\in\mathbb{N}$ we have

[TABLE]

and for all $x\in\mathbb{R}$

[TABLE]

Therefore

[TABLE]

Proof.

We prove both claims by induction. $\partial_{x}^{0}e^{-x^{2}}=e^{-x^{2}}$, $\partial_{x}^{1}e^{-x^{2}}=-2xe^{-x^{2}}$ and $\partial_{x}^{2}e^{-x^{2}}=4x^{2}e^{-x^{2}}-2e^{-x^{2}}$, so (21) holds for $k=1$. Suppose (21) holds for $k$, we prove the inductive step as follows:

[TABLE]

Similarly, observe that (22) holds for $k=0$ and $k=1$. Suppose (22) holds for $k$, then we show the induction step as follows:

[TABLE]

Finally, using the previous two statements we can prove (23) by the following calculation:

[TABLE]

∎

Lemma 47**.**

Suppose we know a $K>1$ such that $\left\lVert H\right\rVert\leq K$. If $\varepsilon\in(0,1/3)$, then we can implement a unitary $\tilde{U}$ such that $\left\lVert(\langle 0|\otimes I)\tilde{U}(|0\rangle\otimes I)-e^{-H^{2}}/2\right\rVert\leq\varepsilon$ using

[TABLE]

Proof.

We apply Corollary 43 to the function $e^{-x^{2}}$. For this let $L_{0}:=\lceil 32K\rceil$, $L:=2L_{0}+1$ and let $x_{\ell}:=(\ell-1-L_{0})/32$ for all $\ell\in[L]$. We choose $r:=1/32$, $\delta:=\theta:=1/128$ and $B=2$ so that the conditions of Corollary 43 are satisfied, as shown by Lemma 46. Indeed $r+\delta\leq 1/(8e)$, hence for all $\ell\in[L]$ we have $\sum_{k=0}^{\infty}a_{k}^{(\ell)}(r+\delta)^{k}\leq 2=B$ as we can see by (23). Since $\delta=\Theta(1)$, Corollary 43 provides the desired complexity. ∎

We can use the above lemma to prepare Gibbs states in a similar way to Lemma 45. In case we have access to $\sqrt{H}$, the advantage of this method is that the dependence on $\left\lVert H\right\rVert$ is reduced to $\sqrt{\left\lVert H\right\rVert}$.

Improved HHL algorithm.

Our techniques can also be applied to the Harrow-Hassidim-Lloyd (HHL) algorithm [HHL09] in order to gain improvements in a similar manner to Childs et al. [CKS17]. The problem the HHL algorithm solves is the following. Suppose we have a circuit $U$ preparing a quantum state $|b\rangle$ (say, starting from $|0\rangle$), and have $d$-sparse oracle access to a non-singular Hamiltonian $H$. The task is to prepare a quantum state, that $\varepsilon$-approximates $H^{-1}|b\rangle/\left\lVert H^{-1}|b\rangle\right\rVert$. For simplicity here we only count the number of uses of $U$ and the number of queries to $H$. Childs et al. [CKS17] present two different methods for achieving this, one based on Hamiltonian simulation, and another directly based on quantum walks. Under the conditions $\left\lVert H\right\rVert\leq 1$ and $\left\lVert H^{-1}\right\rVert\leq\kappa$, the former makes ${\mathcal{O}}\left(\kappa\sqrt{\log(\kappa/\varepsilon)}\right)$ uses of $U$ and has query complexity ${\mathcal{O}}(d\kappa^{2}\log^{2.5}(\kappa/\varepsilon))$. The latter makes ${\mathcal{O}}\left(\kappa\log(d\kappa/\varepsilon)\right)$ uses of $U$ and has query complexity ${\mathcal{O}}(d\kappa^{2}\log^{2}(d\kappa/\varepsilon))$.

Now we provide a sketch of how to solve the HHL problem with ${\mathcal{O}}(\kappa)$ uses of $U$ and with query complexity ${\mathcal{O}}(d\kappa^{2}\log^{2}(\kappa/\varepsilon))$ using our techniques. The improvement on both previously mentioned results is not very large, but our proof is significantly shorter thanks to our general Theorem 41 and Corollary 43.

To solve the HHL problem we need to implement the function $H^{-1}$, i.e., apply the function $f(x)=1/x$ to $H$. Due to the constraints on $H$, the eigenvalues of $H$ lie in the union of the intervals $[-1,-1/\kappa]$ and $[1/\kappa,1]$. We first assume that the eigenvalues actually lie in $[1/\kappa,1]$. In this case we can easily implement the function $1/x$ by Theorem 41 using the Taylor series around $1$:

[TABLE]

As $H^{-1}=(I+(H-I))^{-1}$, we are interested in the eigenvalues of $H-I$. The eigenvalues of $H-I$ lie in the interval $[-1+1/\kappa,0]$, so we choose $r:=1-1/\kappa$ and $\delta:=1/(2\kappa)$. By substituting $z:=-1+1/(2\kappa)$ in (24), we can see that $B:=2\kappa$ satisfies the conditions of Theorem 41. Let $\varepsilon^{\prime}\in(0,1/2)$, then Theorem 41 provides an ${\mathcal{O}}\left(d\kappa\log\left(\kappa/\varepsilon^{\prime}\right)\log\left(1/\varepsilon^{\prime}\right)\right)$-query implementation of an $\varepsilon^{\prime}$-approximation of the operator $H^{-1}/(2\kappa)$, since $\left\lVert H\right\rVert\leq 1$. We can proceed similarly when the eigenvalues of $H-I$ lie in the interval $[-1,-1/\kappa]$, and we can combine the two cases using Corollary 43.

Setting $\varepsilon^{\prime}:=c\varepsilon/\kappa$ for an appropriate constant $c$ and using amplitude amplification, we can prepare an $\varepsilon$-approximation of the state $H^{-1}|b\rangle/\left\lVert H^{-1}|b\rangle\right\rVert$ as required by HHL using ${\mathcal{O}}(\kappa)$ amplitude amplification steps. Therefore we use $U$ at most ${\mathcal{O}}(\kappa)$ times and make ${\mathcal{O}}\left(d\kappa^{2}\log^{2}\left(\kappa/\varepsilon\right)\right)$ queries.

Appendix C Generalized minimum-finding algorithm

In this appendix we describe our generalized quantum minimum-finding algorithm, which we are going to apply to finding an approximation of the ground state energy of a Hamiltonian. This algorithm generalizes the results of Dürr and Høyer [DH96] in a manner similar to the way amplitude amplification [BHMT02] generalizes Grover search: we do not need to assume the ability to query individual elements of the search space, we just need to be able to generate a superposition over the search space. The algorithm also has the benefit over binary search that it removes a logarithmic factor from the complexity.

The backbone of our analysis will be the meta-algorithm below from [DH96]. The meta-algorithm finds the minimal element in the range of the random variable $X$ by sampling, where by “range” we mean the values which occur with non-zero probability. We assume $X$ has finite range.

Note that the above algorithm will always find the minimum, since the obtained samples are strictly decreasing.

Lemma 48**.**

Let $X$ be a finite discrete random variable whose range of values is $x_{1}<x_{2}<\ldots<x_{N}$. Let $S(X)=\{s_{1},s_{2},\dots\}$ denote the random set of values obtained via sampling during a run of Meta-Algorithm 2 with input random variable $X$. If $k\in[N]$, then

[TABLE]

Proof.

The intuition of the proof is to show that, whenever $\mathrm{Pr}(s_{t-1}>x_{k})>0$, we have

[TABLE]

To formulate the statement more precisely131313Throughout this proof, whenever a fraction is $0/0$, we simply interpret it as [math]. Therefore we also interpret conditional probabilities, conditioned on events that happen with probability [math], as [math]. we consider a fixed value $t\in[N]$. For notational convenience let $x:=x_{k}$, $x_{N+1}:=\infty$, we prove (25) by:

[TABLE]

This is enough to conclude the proof, since there is always a smallest $t\in[N]$ such that $s_{t}\leq x$, as the algorithm always finds the minimum in at most $N$ steps. So we can finish the proof by

[TABLE]

∎

Now we describe our generalized minimum-finding algorithm which is based on Meta-Algorithm 2. We take some unitary $U$, and replace $X$ by the distribution obtained if we measured the second register of $U|0\rangle$. We implement conditional sampling via amplitude amplification and the use of the exponential search algorithm of Boyer et al. [BBHT98]. If a unitary prepares the state $|0\rangle|\phi\rangle+|1\rangle|\psi\rangle$ where $\left\lVert\phi\right\rVert^{2}+\left\lVert\psi\right\rVert^{2}=1$, then this exponential search algorithm built on top of amplitude amplification prepares the state $|1\rangle|\psi\rangle$ probabilistically using an expected number of $\mathcal{O}\left(1/\left\lVert\psi\right\rVert\right)$ applications of $U$ and $U^{-1}$ (we will skip the details here, which are straightforward modifications of [BBHT98]).

Lemma 49**.**

There exists $C\in\mathbb{R_{+}}$, such that if we run Algorithm 3 indefinitely (setting $M=\infty$), then for every $U$ and $x_{k}$ the expected number of uses of $U$ and $\kern 0.56905ptU^{-1}$ before obtaining a sample $x\leq x_{k}$ is at most $\frac{C}{\sqrt{\mathrm{Pr}(X\leq x_{k})}}$.

Proof.

Let $X_{<x_{\ell}}$ denote the random variable for which $\mathrm{Pr}(X_{<x_{\ell}}=x)=\mathrm{Pr}(X=x\mid X<x_{\ell})$. The expected number of uses of $U^{\pm 1}$ in Algorithm 3 before obtaining a value $x\leq x_{k}$ is

[TABLE]

where the last equality follows from Equation (30) below. The constant $C$ from the lemma is the constant hidden by the ${\mathcal{O}}$. The remainder of this proof consists of proving (30) using elementary calculus. Let us introduce the notation $p_{0}:=\mathrm{Pr}(X\leq x_{k})$ and for all $j\in[N-k]$ let $p_{j}:=\mathrm{Pr}(X=x_{k+j})$. Then the expression inside the ${\mathcal{O}}$ on the second-to-last line above becomes

[TABLE]

The basic idea is that we treat the expression on the right-hand side of (26) as an integral approximation sum for the integral $\int_{p_{0}}^{1}z^{-3/2}dz$, and show that it is actually always less than the value of this integral. We proceed by showing that subdivision always increases the sum.

Let us fix some $\ell\in[N-k]$ and define

[TABLE]

and observe that

[TABLE]

We show that (27) is $\geq 0$ after simplifying the expression by substituting $a:=\sum_{i=0}^{\ell-1}p_{i}$ and $b:=p_{\ell}/2$:

[TABLE]

Let us fix some parameter $\delta>0$. Recursively applying this halving procedure for different indices, we can find some $J\in\mathbb{N}$ and $\tilde{p}\in\mathbb{R}_{+}^{J+1}$ such that $\sum_{j=0}^{J}\tilde{p}_{j}=1$, $\tilde{p}_{0}=p_{0}$ and $\tilde{p}_{j}\leq\delta$ for all $j\in[J]$, moreover

[TABLE]

Observe that for all $j\in[J]$

[TABLE]

Therefore

[TABLE]

Since this inequality holds for every $\delta>0$, we can conclude using (26) that

[TABLE]

∎

It is not too hard to work out the constant by following the proof of [BBHT98] providing something like $C\approx 25$. The following theorem works with any $C$ satisfying Lemma 49.

Theorem 50** (Generalized Minimum-Finding).**

If we run Algorithm 3 with input satisfying $M\geq 4C/\sqrt{\mathrm{Pr}(X\leq x)}$ for $C$ as in Lemma 49 and a unitary $U$ that acts on $q$ qubits, then at termination we obtain an $x_{i}$ from the range of $X$ that satisfies $x_{i}\leq x$ with probability at least $\frac{3}{4}$. Moreover the success probability can be boosted to at least $1-\delta$ with ${\mathcal{O}}(\log(1/\delta))$ repetitions. This uses at most $M$ applications of $U$ and $U^{-1}$ and ${\mathcal{O}}(qM)$ other gates.

Proof.

Let $x_{k}$ be the largest value in the range of $X$ such that $x_{k}\leq x$. Then Lemma 49 says that the expected number of applications of $U$ and $U^{-1}$ before finding a value $x_{i}\leq x_{k}$ is at most $C/\sqrt{\mathrm{Pr}(X\leq x_{k})}=C/\sqrt{\mathrm{Pr}(X\leq x)}$, therefore by the Markov inequality we know that the probability that we need to use $U$ and $U^{-1}$ at least $4C/\sqrt{\mathrm{Pr}(X\leq x)}$ times is at most $1/4$. The boosting of the success probability can be done using standard techniques, e.g., by repeating the whole procedure ${\mathcal{O}}(\log(1/\delta))$ times and taking the minimum of the outputs.

The number of applications of $U$ and $U^{-1}$ follows directly form the algorithms description. Then, for the number of other gates, each amplitude amplification step needs to implement a binary comparison and a reflection through the $|0\rangle$ state, both of which can be constructed using ${\mathcal{O}}(q)$ elementary gates, giving a total of ${\mathcal{O}}(qM)$ gates. ∎

Note that this result is a generalization of Dürr and Høyer [DH96]: if we can create a uniform superposition over $N$ values $x_{1}<x_{2}<\ldots<x_{N}$, then $\mathrm{Pr}(X\leq x_{1})=1/N$ and therefore Theorem 50 guarantees that we can find the minimum with high probability with ${\mathcal{O}}(\sqrt{N})$ steps.

Now we describe an application of this generalized search algorithm that we need in the paper.

This final lemma in this appendix describes how to estimate the smallest eigenvalue of a Hamiltonian. A similar result was shown by Poulin and Wocjan [PW09a], but we improve on the analysis to fit our framework better. We assume sparse oracle access to the Hamiltonian $H$ as described in Section 2, and will count queries to these oracles. We use some of the techniques introduced in Appendix B.

Lemma 51**.**

If $H\!=\!\sum_{j=1}^{n}E_{j}|\phi_{j}\rangle\!\langle\phi_{j}|$, with eigenvalues $E_{1}\leq E_{2}\leq\ldots\leq E_{n}$, is such that $\left\lVert H\right\rVert\leq K$, $\varepsilon\leq K/2$, and $H$ is given in $d$-sparse oracle form, then we can obtain an estimate $E$ such that $\left|E_{1}-E\right|\leq\varepsilon$, with probability at least $2/3$, using

[TABLE]

Proof.

The general idea is as follows: we prepare a maximally entangled state on two registers, and apply phase estimation [NC00, Section 5.2][CEMM98] to the first register with respect to the unitary $e^{\pi iH/K}$. We then use Theorem 50 to find the minimal phase. In order to guarantee correctness we need to account for all the approximation errors coming from approximate implementations. This causes some technical difficulty, since the approximation errors can introduce phase estimates that are much less than the true minimum. We need to make sure that the minimum-finding algorithm finds these faulty estimates only with a tiny probability.

We first initialize two $\log(n)$-qubit registers in a maximally entangled state $\smash{\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}|j\rangle|j\rangle}$. This can be done for example using $\log(n)$ Hadamard and CNOT gates, when $n$ is a power of two. Due to the invariance of maximally entangled states under transformations of the form $W\otimes W^{*}$ for unitary $W$, we have that

[TABLE]

Let $T:=2^{\left\lceil\log\left(\frac{K}{\varepsilon}\right)+2\right\rceil}$ and first assume that we have access to a perfect unitary $V$ which implements $V=\sum_{t=0}^{T-1}|t\rangle\!\langle t|\otimes e^{\pi tiH/K}$. Let $e_{j}:=E_{j}T/2K$. If we apply phase estimation to the quantum state $|\phi_{j}\rangle$ using $V$, then we get some superposition of phase estimates $|e\rangle$ in the first register. This superposition has the property, that a measurement in the computational basis reveals an $e$ such that $|e-e_{j}|\leq 3$ with high probability, so that the estimate $E:=e2K/T$ would satisfy $|E-E_{j}|=|e-e_{j}|2K/T\leq 3\varepsilon/4<\varepsilon$. If we repeat phase estimation ${\mathcal{O}}(\log(n))$ times (on the same state $|\phi_{j}\rangle$), and take the median of the estimates (in superposition), then we obtain a more concentrated superposition of estimates $e$ such that a measurement would reveal an $|e-e_{j}|\leq 3$ with probability at least $1-b/n$, for some $b=\Theta(1)$.

Since in our maximally entangled state $|\phi_{j}\rangle$ is entangled with $|\phi^{*}_{j}\rangle$ on the second register, applying phase estimation to the first register in superposition does not cause interference. Let us denote by $U$ the unitary corresponding to the above preparation-estimation-boost procedure. Define $\Pi$ to be the projector which projects to the subspace of estimation values $|e\rangle$ such that there is a $j\in[n]$ with $|e-e_{j}|\leq\varepsilon$. By the non-interference argument we can see that, after applying $U$, the probability that we get an estimation $e$ such that $|e-e_{j}|>3$ is at most $b/n$ for all $j\in[n]$, moreover $\left\lVert(I-\Pi)U|0\rangle\right\rVert^{2}\leq b/n$. Also let $\Pi_{1}$ denote the projector which projects to phase estimates that yield $e$ such that $|e-e_{1}|\leq 3$. It is easy to see that $\left\lVert\Pi_{1}U|0\rangle\right\rVert^{2}\geq 1/n-b/n^{2}$.

Now let us replace $V$ by $\tilde{V}$ implemented via Lemma 37, such that $\left\lVert\smash{V-\tilde{V}}\right\rVert\leq c^{\prime}/(n\log(n))$ for some $c^{\prime}=\Theta(1)$. Let $\tilde{U}$ denote the circuit that we obtain from $U$ by replacing $V$ with $\tilde{V}$. Since in the repeated phase-estimation procedure we use $V$ in total ${\mathcal{O}}(\log(n))$ times, by using the triangle inequality we see that $\left\lVert\smash{U-\tilde{U}}\right\rVert\leq c/(2n)$, where $c=\Theta(1)$. We use the well-known fact that if two unitaries are $\delta$-close in operator norm, and they are applied to the same quantum state, then the measurement statistics of the resulting states are $2\delta$-close. Therefore we can upper bound the difference in probability of getting outcome $(I-\Pi)$:

[TABLE]

hence $\left\lVert\smash{(I-\Pi)\tilde{U}|0\rangle}\right\rVert^{2}\leq(b+c)/n$, and we can prove similarly that $\left\lVert\smash{\Pi_{1}\tilde{U}|0\rangle}\right\rVert^{2}\geq 1/n-(b+c)/n$.

Now let $|\psi\rangle:=(I-\Pi)\tilde{U}|0\rangle/\left\lVert\smash{(I-\Pi)\tilde{U}|0\rangle}\right\rVert$ be the state that we would get after post-selecting on the $(I-\Pi)$-outcome of the projective measurement $\Pi$. For small enough $b,c$ we have that $\left\lVert\smash{|\psi\rangle-\tilde{U}|0\rangle}\right\rVert={\mathcal{O}}(\sqrt{(b+c)/n})$ by the triangle inequality. Thus there exists an idealized unitary $U^{\prime}$ such that $|\psi\rangle=U^{\prime}|0\rangle$, and $\left\lVert\smash{\tilde{U}-U^{\prime}}\right\rVert={\mathcal{O}}(\sqrt{(b+c)/n})$. Observe that $\left\lVert\Pi_{1}U^{\prime}|0\rangle\right\rVert^{2}=\left\lVert|\psi\rangle\right\rVert^{2}\geq\left\lVert\smash{\Pi_{1}\tilde{U}|0\rangle}\right\rVert^{2}\geq 1/n-(b+c)/n$.

Now suppose $(b+c)\leq 1/2$ and we run the generalized minimum-finding algorithm of Theorem 50 using $U^{\prime}$ with $M=6C\sqrt{n}$. Since

[TABLE]

we will obtain an estimate $e$ such that $e\leq e_{1}+3$, with probability at least $3/4$. But since $\Pi|\psi\rangle=|\psi\rangle$, we find that any estimate that we might obtain satisfies $e\geq e_{1}-3$. So an estimate $e\leq e_{1}+3$ always satisfies $|e-e_{1}|\leq 3$.

The problem is that we only have access to $\tilde{U}$ as a quantum circuit. Let $C_{MF}(\tilde{U})$ denote the circuit that we get from Theorem 50 when using it with $\tilde{U}$ and define similarly $C_{MF}(U^{\prime})$ for $U^{\prime}$. Since we use $\tilde{U}$ a total of ${\mathcal{O}}(\sqrt{n})$ times in $C_{MF}(\tilde{U})$ and

[TABLE]

Therefore the measurement statistics of the two circuits differ by at most ${\mathcal{O}}(\sqrt{b+c})$. Choosing $b,c$ small enough constants ensures that $C_{MF}(\tilde{U})$ outputs a proper estimate $e$ such that $|e-e_{1}|\leq 3$ with probability at least $2/3$. As we have shown at the beginning of the proof, such an $e$ yields an $\varepsilon$-approximation of $E_{1}$ via $E:=e2K/T$.

The query complexity has an ${\mathcal{O}}(Td\log(Tn))={\mathcal{O}}(Kd/\varepsilon\log(Kn/\varepsilon))$ factor coming from the implementation of $\tilde{V}$ by Lemma 37. This gets multiplied with ${\mathcal{O}}(\log(n))$ by the boosting of phase estimation, and by ${\mathcal{O}}(\sqrt{n})$ due to the minimum-finding algorithm. The gate complexity is dominated by the cost ${\mathcal{O}}(Kd/\varepsilon\log^{7/2}(Kn/\varepsilon))$ of implementing $\tilde{V}$, multiplied with the ${\mathcal{O}}(\sqrt{n}\log(n))$ factor as for the query complexity. ∎

Note that the minimum-finding algorithm of Theorem 50 can also be used for state preparation. If we choose $2\varepsilon$ less than the energy-gap of the Hamiltonian, then upon finding the approximation of the ground state energy we also prepare an approximate ground state. The precision of this state preparation can be improved with logarithmic cost, as can be seen from the proof of Lemma 51.

Appendix D Sparse matrix summation

As seen in Section 2, the Arora-Kale algorithm requires an approximation of $\exp(-\eta H^{(t)})$ where $H^{(t)}$ is a sum of matrices. To keep this section general we simplify the notation. Let $H$ be the sum of $k$ different $d$-sparse matrices $M$:

[TABLE]

In this section we study the complexity of one oracle call to $H$, given access to respective oracles for the matrices $M_{1},\ldots,M_{k}$. Here we assume that the oracles for the $M_{i}$ are given in sparse matrix form, as defined in Section 2. In particular, the goal is to construct a procedure that acts as a similar sparse matrix oracle for $H$. We will only focus on the oracle that computes the non-zero indices of $H$, since the oracle that gives element access is easy to compute by summing the separate oracles.

In the remainder of this section we only consider one row of $H$. We denote this row by $R_{H}$ and the corresponding rows of the matrices $M_{i}$ by $R_{i}$. Notice that such a row is given as an ordered list of integers, where the integers are the non-zero indices in $R_{i}$. Then $R_{H}$ will again be an ordered list of integers, containing all integers in the $R_{i}$ lists once (i.e., $R_{H}$ does not contain duplicates).

D.1 A lower bound

We show a lower bound on the query complexity of the oracle $O^{I}_{H}$ described above by observing that determining the number of elements in a row of $H$ solves the majority function. Notice that, given access to $O^{I}_{H}$, we can decide whether there are at least a certain number of non-zero elements in a row of $H$.

Lemma 52**.**

Given $k+1$ ordered lists of integers $R_{0},\dots,R_{k}$, each of length at most $d$. Let $R_{H}$ be the merged list that is ordered and contains every element in the lists $R_{i}$ only once (i.e., we remove duplicates). Deciding whether $|R_{H}|\leq d+\frac{dk}{2}$ or $|R_{H}|\geq d+\frac{dk}{2}+1$ takes $\Omega\left(dk\right)$ quantum queries to the input lists in general.

Proof.

We prove this by a reduction from MAJ on $dk$ elements. Let $Z\in\{0,1\}^{d\times k}$ be a Boolean string. It is known that it takes at least $\Omega(dk)$ quantum queries to $Z$ to decide whether $|Z|\leq\frac{dk}{2}$ or $|Z|\geq\frac{dk}{2}+1$. Now let $R_{0},R_{1},\dots,R_{k}$ be lists of length $d$ defined as follows:

• •

$R_{0}[j]=j(k+1)$ for $j=1,\ldots,d$.

• •

$R_{i}[j]=j(k+1)+jZ_{ij}$ for $j=1,\ldots,r$ and $i=1,\ldots,k$.

By construction, if $Z_{ij}=1$, then the value of the entry $R_{i}[j]$ is unique in the lists $R_{0},\ldots,R_{k}$, and if $Z_{ij}=0$ then $R_{i}[j]=R_{0}[j]$. So in $R_{H}$ there will be one element for each element in $R_{0}$ and one element for each bit in $Z_{ij}$ that is one. The length of $R_{H}$ is therefore $d+|Z|$. Hence, distinguishing between $|R_{H}|\leq d+\frac{dk}{2}$ and $|R_{H}|\geq d+\frac{dk}{2}+1$ would solve the MAJ problem and therefore requires at least $\Omega(dk)$ queries to the lists in general. ∎

Corollary 53**.**

Implementing a query to a sparse matrix oracle $O^{I}_{H}$ for

[TABLE]

where each $M_{j}$ is $d$-sparse, requires $\Omega\left(dk\right)$ queries to the $O^{I}_{M_{j}}$ in general.

D.2 An upper bound

We first show that an oracle for the non-zero indices of $H$ can be constructed efficiently classically. The important observation is that in the classical case we can write down the oracle once and store it in memory. Hence, we can create an oracle for $H$ as follows. We start from the oracle of $M_{1}$ and then we “add” the oracle of $M_{2}$, then that of $M_{2}$, etc. By “adding” the oracle $M_{i}$ to $H=\sum_{j=1}^{i-1}M_{j}$, we mean that, per row, we insert the non-zero indices in the list of $M_{i}$ into that of $H$ (if it is not already there). When an efficient data structure (for example a binary heap) is used, then such insertions can be done in polylog time. This shows that in the classical case such an oracle can be made in time $\widetilde{\mathcal{O}}\left(ndk\right)$. Note that in the application we are interested in, Meta-algorithm 1, in each iteration $t$ only one new matrix $M^{(t)}$ ‘arrives’, hence from the oracle for $H^{(t-1)}$ an oracle for $H^{(t)}$ can be constructed in time $\widetilde{\mathcal{O}}\left(nd\right)$.

The quantum case is similar, but we need to add all the matrices together each time a query to $O^{I}_{H}$ is made, since writing down each row of $H$ in every iteration would take $\Omega(n)$ operations.

To implement one such query to $O^{I}_{H}$, in particular to the $t$th entry of $R_{H}$, start with an empty heap and add all elements of $R_{1}$ to it. Continue with the elements of $R_{2}$, but this time, for each element first check if it is already present, if not, add it, if it is, just continue. Overall this will take ${\mathcal{O}}(dk)$ insertions and searches in the data structure and hence $\widetilde{\mathcal{O}}\left(dk\right)$ operations. We end up with a full description of $R_{H}$. We can then find the index of the $t$th non-zero element from this and uncompute the whole description. Similarly, we need to be able to compute the inverse function since we need an in-place calculation. Given an index $i$ of a non-zero element in $R_{H}$, we can compute all the indices for $R_{H}$ as above, and find where $i$ is in the heap to find the corresponding value of $t$.

Appendix E Equivalence of $R$, $r$, and $\varepsilon^{-1}$

In this section we will prove the equivalence of the three parameters $R$, $r$ and $\varepsilon^{-1}$ in the Arora-Kale meta-algorithm. That is, we will show any two of the three parameters can be made constant by increasing the third. Therefore, $\frac{Rr}{\varepsilon}$ as a whole is the interesting parameter. This appendix is structured as a set of reductions, in each case we will denote the parameters of the new SDP with a tilde.

Lemma 54**.**

Let $0<\varepsilon\leq 1$. For every SDP with $R,r\geq 1$, there is an SDP with parameters $\tilde{R}=1$ and $\tilde{r}=r$ such that a solution to that SDP with precision $\tilde{\varepsilon}=\frac{\varepsilon}{R}$ provides a solution with precision $\varepsilon$ to the original SDP.

Proof.

Let $\tilde{A}_{j}=A_{j}$, $\tilde{C}=C$ and $\tilde{b}=\frac{b}{R}$. Now clearly $\tilde{R}=1$, but $\widetilde{\mbox{\rm OPT}}=\mbox{\rm OPT}/R$. Hence determining $\widetilde{\mbox{\rm OPT}}$ up to additive error $\tilde{\varepsilon}=\frac{\varepsilon}{R}$ will determine OPT up to additive error $\varepsilon$. Notice that the feasible region of the dual did not change, so $\tilde{r}=r$. ∎

Lemma 55**.**

Let $0<\varepsilon\leq 1$. For every SDP with $R,r\geq 1$, there is an SDP with parameters $\tilde{R}=\frac{R}{\varepsilon}$ and $\tilde{r}=r$ such that a solution to that SDP with precision $\tilde{\varepsilon}=1$ provides a solution with precision $\varepsilon$ to the original SDP.

Proof.

Let $\tilde{A}_{j}=A_{j}$, $\tilde{C}=C$ and $\tilde{b}=\frac{b}{\varepsilon}$. Now $\tilde{R}=\frac{R}{\varepsilon}$ and $\widetilde{\mbox{\rm OPT}}=\mbox{\rm OPT}/\varepsilon$. Hence determining $\widetilde{\mbox{\rm OPT}}$ up to additive error $\tilde{\varepsilon}=1$ will determine OPT up to additive error $\varepsilon$. Notice that again the feasible region of the dual did not change, so $\tilde{r}=r$. ∎

Lemma 56**.**

Let $0<\varepsilon\leq 1$. For every SDP with $R,r\geq 1$, there is an SDP with parameters $\tilde{R}=R$ and $\tilde{r}=1$ such that a solution to that SDP with precision $\tilde{\varepsilon}=\frac{\varepsilon}{r}$ provides a solution with precision $\varepsilon$ to the original SDP.

Proof.

Let $\tilde{A}_{j}=A_{j}$, $\tilde{b}=b$ and $\tilde{C}=\frac{1}{r}C$. Now $\tilde{r}=1$ and $\widetilde{\mbox{\rm OPT}}=\mbox{\rm OPT}/r$. Hence determining $\widetilde{\mbox{\rm OPT}}$ up to additive error $\tilde{\varepsilon}=\frac{\varepsilon}{r}$ will determine OPT up to additive error $\varepsilon$. Since $r\geq 1$ and $\left\lVert C\right\rVert\leq 1$, we find $\left\lVert\smash{\tilde{C}}\right\rVert\leq 1$ as required. Notice that the feasible region of the primal did not change, so $\tilde{R}=R$. ∎

At this point we would like to state the last reduction by setting $\tilde{C}=\frac{1}{\varepsilon}C$, but this would not guarantee that $\left\lVert\smash{\tilde{C}}\right\rVert\leq 1$. Instead we give an algorithm that performs multiple calls to an SDP-solver, each of which has constant $\varepsilon$ but higher $r$.

Lemma 57**.**

Assume an SDP-solver that costs

[TABLE]

to solve one SDP with precision $\varepsilon$, and assume that $C$ is non-decreasing in $r$. Every SDP with parameters $R,r\geq 1$ can be solved up to precision $0<\varepsilon\leq 1$ with cost

[TABLE]

by solving $\log\left(\frac{1}{\varepsilon}\right)$ SDPs, where the $k$-th SDP is solved up to precision $\tilde{\varepsilon}=1$ and has parameters

[TABLE]

and input matrices with elements described by bitstrings of length $\text{poly}(\log n,\log m,\log\left(\frac{1}{\varepsilon}\right))$. Furthermore, if $C(n,m,s,R,1,r)=\text{poly}(n,m,s,R,r)$, then the above cost is $\text{poly}(n,m,s,R,r/eps)$.

Proof.

The high-level idea is that we want to learn a small interval in which the optimum lies, whilst using a “big” precision of $1$. We do so as follows: given an interval $[L,U]$ with the promise that $\mbox{\rm OPT}\in[L,U]$, we formulate another SDP for which a $1$-approximation of the optimum learns us a new, smaller, interval $[L^{\prime},U^{\prime}]$ such that $\mbox{\rm OPT}\in[L^{\prime},U^{\prime}]$. We will moreover have $U^{\prime}-L^{\prime}\leq\frac{1}{2}(U-L)$. In the remainder of the proof we first show how to do this reformulation, we then use this technique to prove the lemma.

Given an SDP $p$, given in the form of Equation (1), of which we know an interval $[L,U]$ such that $\mbox{\rm OPT}\in[L,U]$ (with $0<U-L\leq 1$), we can write down an equivalent SDP $p^{\prime}$ such that an optimal solution of $p$ corresponds one-to-one to an optimal solution of $p^{\prime}$, and the optimum of $p^{\prime}$ lies in $[0,4]$:

[TABLE]

here the variable $\tilde{X}$ is of size $(n+1)\times(n+1)$, and should be thought of as $\tilde{X}=\begin{pmatrix}X&\cdot\\ \cdot&z\end{pmatrix}$, where $X$ is the variable of the original SDP: an $n\times n$ positive semidefinite matrix. Observe that by assumption, for every feasible $X$ of the original SDP, $L\leq\mbox{\rm Tr}\left(CX\right)\leq U$. Therefore, the first constraint implies $0\leq z\leq 4$ and hence the new optimum lies between [math] and $4$, and the new trace bound is $\tilde{R}=R+4$. We now determine $\tilde{r}$. The dual of the above program is given by:

[TABLE]

Claim 58**.**

An optimal solution $\tilde{y}$ to $d^{\prime}$ is of the form $\tilde{y}=\frac{4}{U-L}(1,y)$ where $y$ is an optimal solution to $d$, the dual of $p$.

The proof of the claim is deferred to the end of this section. The claim implies that $\tilde{r}=\frac{4}{U-L}(1+r)$. We also have

[TABLE]

and hence, a $1$-approximation to $\mbox{\rm OPT}^{\prime}$ gives a $\frac{U-L}{4}$-approximation to OPT.

We now use the above technique to prove the lemma. Assume an SDP of the form (1) is given. By assumption, $\left\lVert C\right\rVert\leq 1$ and therefore $\mbox{\rm OPT}\in[-R,R]$. Calling the SDP-solver on this problem with $\varepsilon=1$ will give us an estimate of OPT up to additive error $1$. Call this estimate $\mbox{\rm OPT}_{0}$, then $\mbox{\rm OPT}\in[\mbox{\rm OPT}_{0}-1,\mbox{\rm OPT}_{0}+1]=:[L_{0},U_{0}]$. We now define a new SDP $p^{\prime}$ as above, with $U=U_{0},L=L_{0}$ (notice $U_{0}-L_{0}\leq 2$). By the above, solving $p^{\prime}$ with $\tilde{\varepsilon}=1$ determines a new interval $[L_{1},U_{1}]$ of length at most $2\frac{U_{0}-L_{0}}{4}$ such that $\mbox{\rm OPT}\in[L_{1},U_{1}]$. We use the interval $[L_{1},U_{1}]$ to build a new SDP $p^{\prime}$ and solve that with $\tilde{\varepsilon}=1$ to get an interval $[L_{2},U_{2}]$. Repeating this procedure $k=\log\left(\frac{1}{\varepsilon}\right)+1$ times determines an interval $[L_{k},U_{k}]$ of length at most $\frac{1}{2^{k}}(U_{0}-L_{0})\leq\varepsilon$. Hence, we have determined the optimum of $p$ up to an additive error of $\varepsilon$. The total time needed for this procedure is at most

[TABLE]

Proof of Claim 58.

First observe that the linear matrix inequality in $d^{\prime}$ implies the inequality $y_{0}\frac{U-L}{4}\geq 1$ and hence $y_{0}\geq\frac{4}{U-L}$. Suppose we fix a value $y_{0}$ (with $y_{0}\geq\frac{4}{U-L}$) in $d^{\prime}$, then the resulting SDP is of the form

[TABLE]

and hence, an optimal solution $\tilde{y}$ to $d^{\prime\prime}$ is of the form $\tilde{y}=y_{0}y$ where $y$ is an optimal solution to $d$. It follows that an optimal solution to $d^{\prime}$ is of the form $(y_{0},y_{0}y)$ where $y$ is an optimal solution to $d$. Observe that the optimal value of $d^{\prime}$ as a function of $y_{0}$ is of the form $y_{0}\cdot(-L+\mbox{\rm OPT})$. Since by assumption $\mbox{\rm OPT}\geq L$, the objective is increasing (linearly) with $y_{0}$ and hence $y_{0}=\frac{4}{U-L}$ is optimal. ∎