$\tilde{\rho}=\frac{1}{1+n\delta}\left[\left(\sum_{k=0}^{K}\frac{1}{\eta^{(k)}}Q^{(k)}\circ\frac{\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)}{\operatorname{tr}\left(\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)\right)}\right)+\delta
I\right],$ (23)

where $\delta$ is defined according to (18) (see, e.g., Proposition 10). We
point out that this marks a key difference between the output of our algorithm
and other quantum SDO solvers based on Gibbs sampling [9, 10, 11, 60, 61],
which need only return a single state preparation pair. This however does not
increase the cost of the method; the iteration bound in Corollary 5 ensures
that there are only at most $\widetilde{{\cal O}}_{\frac{1}{\zeta}}(1)$ (i.e.,
a poly-logarithmic number) of these tuples to be stored over the course of the
algorithm. Using the QRAM input model, one can use the stored tuples to
construct a block-encoding of the final solution up to error $\theta$ using
$\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon},\frac{1}{\theta}}\left(\sqrt{n}\right)$
queries to the QRAM and $\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}\left(n\right)$ classical operations. This
construction, and the associated time complexity are analyzed later in
Proposition 12. We further demonstrate that provided classical access to an
$s$-sparse matrix $A\in\mathbb{R}^{n\times n}$ (with subnormalization factor
1) and access to QRAM, one can estimate $\operatorname{tr}(A\tilde{\rho})$ to
additive error $\theta$ using $\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}\left(\frac{\sqrt{n}}{\theta}\right)$
queries to the QRAM and $\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}\left(ns\right)$ classical operations. If
$A$ has a subnormalization factor $\alpha_{A}>1$, then $\theta$ must be scaled
down by $\alpha_{A}$, increasing the cost.

Additionally, we require Algorithm 1 to return the estimates
$\tilde{p}^{(k)}\in\mathbb{R}^{n}$ (a classical estimate of the diagonal
elements of the solution to the refining problem) and
$\tilde{c}^{(k)}\in\mathbb{R}$ (a classical estimate of the objective value
attained by the solution of the refining problem) that are used to test
$\xi$-closeness for the accepted state. In this fashion, we can (classically)
prepare the refining problem data for the next iteration without increasing
the cost of the algorithm with respect to $n$; the objective value can be
updated using ${\cal O}(1)$ arithmetic operations using $\tilde{c}^{(k)}$,
while updating the residuals along the diagonal of $\rho$ requires ${\cal
O}(n)$ arithmetic operations provided classical access to $\tilde{p}^{(k)}$.

If the current solution is indistinguishable up to precision $\zeta$ from the
maximally mixed state $n^{-1}I$, and provides an objective value of at least
$\gamma-\zeta$, the algorithm terminates and reports the current solution.
Otherwise, we construct the refining problem associated with our current
solution and proceed to the next iteration.

Input: Error tolerances $\epsilon\in(0,1)$ and
$\zeta=\left(\frac{\epsilon}{n\|C\|_{F}}\right)^{4}$, upper bound on objective
value $\gamma\in[-1,1]$

Output: Tuples $(\eta^{(k)},y^{(k)},q^{(k)},d^{(k)})$ that define a
$2\zeta$-precise solution to (5) using Equation (23)

Initialize: $p\leftarrow\mathbf{0}^{n}$, $Q\leftarrow ee^{\top}$,
$\varepsilon_{i}=\frac{1}{n}$ for $i\in[n]$, $\hat{\gamma}\leftarrow 0$,
$\eta^{(0)}\leftarrow 1$, $k\leftarrow 0$

while  _$\max\left\\{\gamma-\hat{\gamma},\|\varepsilon\|_{1}\right\\} >\zeta$_
do

  1. 1.

$(y^{(k)},D^{(k)},\tilde{p}^{(k)},\tilde{c}^{(k)})\leftarrow$ Solve (14) to
precision $\xi$ using Algorithm 1

  2. 2.

Store diagonal elements of $Q$ and $D^{(k)}$

$q_{i}^{(k)}\leftarrow Q_{ii},\leavevmode\nobreak\ \quad d_{i}^{(k)}\leftarrow
D_{ii}^{(k)}\quad\text{for}\leavevmode\nobreak\ i\in[n]$

  3. 3.

Store description of solution to the refining problem
$(\eta^{(k)},y^{(k)},q^{(k)},d^{(k)})$

  4. 4.

Update the objective value of solution:

$\hat{\gamma}\leftarrow\hat{\gamma}+\frac{1}{\eta^{(k)}}\tilde{c}^{(k)}$

  5. 5.

Update estimate of diagonal entries:

$p_{i}\leftarrow
p_{i}+\frac{Q_{ii}}{\eta^{(k)}}\tilde{p}_{i}^{(k)}\quad\text{for}\leavevmode\nobreak\
i\in[n]$

  6. 6.

Compute element-wise deviations from the maximally mixed state:

$\varepsilon_{i}\leftarrow
p_{i}-\frac{1}{n}\quad\text{for}\leavevmode\nobreak\ i\in[n]$

  7. 7.

Classically update:

$Q_{ii}\leftarrow\textup{sign}(-\varepsilon_{i})\leavevmode\nobreak\
\text{for}\leavevmode\nobreak\
i\in[n],\quad\eta^{(k+1)}\leftarrow\frac{1}{\max\left\\{\gamma-\hat{\gamma},\|\varepsilon\|_{1}\right\\}}$

  8. 8.

$k\leftarrow k+1$

end while

Algorithm 3 Iterative Refinement for SDO Approximations of QUBOs using a
quantum computer

The next result gives the overall running time required to solve (3) to
additive error ${\cal O}(\epsilon)$ using the QRAM input model.

######  Theorem 6.

Let $C\in{\cal S}^{n}$, $\epsilon\in(0,1)$, and set
$\zeta=\left(\frac{\epsilon}{n\|C\|_{F}}\right)^{4}$. Assume we have classical
access to $C$. Then, in the QRAM input model, Algorithm 3 solves (3) up to
additive error ${\cal O}(\epsilon)$ using

${\cal
O}\left(n^{1.5}\cdot\textup{polylog}\left(n,\|C\|_{F},\frac{1}{\epsilon}\right)\right)$

accesses to the QRAM and ${\cal O}(ns)$ classical arithmetic operations.

The output of the algorithm is a collection of tuples
$\\{(\eta^{(k)},y^{(k)},q^{(k)},d^{(k)})\\}_{k=0}^{K}$ such that

$\tilde{\rho}=\frac{1}{1+n\delta}\left[\left(\sum_{k=0}^{K}\frac{1}{\eta^{(k)}}Q^{(k)}\circ\frac{\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)}{\operatorname{tr}\left(\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)\right)}\right)+\delta
I\right],$

is a $2\zeta$-precise solution to (5), where $\delta$ is defined according to
(18). The entries of $\tilde{\rho}$ can be modified to construct a matrix
$\rho^{*}$ at trace distance ${\cal
O}\left(\frac{\epsilon}{n\|C\|_{F}}\right)$ of $\tilde{\rho}$ in time ${\cal
O}(n^{2})$, such that $n\rho^{*}$ is a feasible point of the SDO problem (3).

###### Proof.

Given that $C$ is an $s$-sparse matrix, we can classically load $C$ in ${\cal
O}(ns)$ time. Similarly, for normalization purposes we classically compute
$\|C\|_{F}$, which requires ${\cal O}(ns)$ arithmetic operations. In each
iteration we use Algorithm 1 to solve (14), and use classical estimates of the
diagonal elements of the refining solution, and a classical estimate of the
objective value attained by the refining solution to update the solution and
data for the refining problem we need to solve in the next iteration.

Letting ${\cal T}_{HU}^{\text{quantum}}$ be the cost of using Algorithm 1 as
an approximate SDO subroutine, by Proposition 8, Algorithm 1 solves (14) to
additive error $\xi$ using at most

${\cal T}_{HU}^{\text{quantum}}=\widetilde{{\cal
O}}_{\frac{n}{\xi}}\left(n^{1.5}\xi^{-5}\right)$

accesses to the QRAM. Classically updating the objective value requires ${\cal
O}(1)$ arithmetic operations while updating the vector $p$ which stores a
classical description of the diagonal elements of our solution as

$p_{i}\leftarrow p_{i}+\frac{Q_{ii}}{\eta^{(k)}}\tilde{p}_{i}^{(k)}$

requires ${\cal O}(n)$ classical arithmetic operations. Likewise,
$\varepsilon$ and $Q$ can each be updated using ${\cal O}(n)$ classical
arithmetic operations, as we only need to store the diagonal elements of $Q$.
This also implies that we can update $Q\circ\frac{C}{\|C\|_{F}}$ using
$\widetilde{{\cal O}}_{n}(n)$ operations, for only the diagonal elements need
to be updated. When compared to loading and normalizing the coefficient matrix
$C$, or our use of Algorithm 1 as a subroutine for solving (14), these
intermediate computation steps are negligible and do not factor into the
overall running time using ${\cal O}$ notation.

By Corollary 5, Algorithm 3 terminates in at most $\widetilde{{\cal
O}}_{\frac{1}{\zeta}}\left(1\right)$ iterations. Therefore, the worst case
complexity of Algorithm 3 can be bounded by

${\cal
O}\left(n^{1.5}\xi^{-5}\cdot\textup{polylog}\left(n,\|C\|_{F},\frac{1}{\epsilon}\right)\right)$

accesses to the QRAM, and ${\cal O}\left(ns\right)$ classical arithmetic
operations to load and normalize $C$. Further, it suffices to use fixed
precision $\xi$ for the every call to Algorithm 1 to reach a final solution
that solves (5) to additive error $\zeta$, as the final solution can always be
arbitrarily precise using a $\widetilde{{\cal O}}_{\frac{1}{\zeta}}(1)$ calls
to Algorithm 1. Since $\xi$ is a fixed constant in the context of Algorithm 3,
the overall running time of Algorithm 3 simplifies to

${\cal
O}\left(n^{1.5}\cdot\textup{polylog}\left(n,\|C\|_{F},\frac{1}{\epsilon}\right)\right).$

accesses to the QRAM, and ${\cal O}\left(ns\right)$ classical arithmetic
operations.

Just as in the proof of Theorem 5, applying Proposition 11 with our choice of
$\zeta=\left(\frac{\epsilon}{n\|C\|_{F}}\right)^{4}$ implies that the above
running time is sufficient to obtain a solution that can be used to solve (3)
up to additive error ${\cal O}(\epsilon)$, and the proof is complete. ∎

We analyze the cost of Algorithm 3 without access to QRAM in Appendix A. Using
the sparse-access input model, one can show that the resulting scheme exhibits
an oracle complexity of

${\cal
O}\left(n^{1.5}s^{0.5+o(1)}\cdot\textup{polylog}\left(n,\|C\|_{F},\frac{1}{\epsilon}\right)\right),$

and requires ${\cal
O}\left(n^{2.5}s^{0.5+o(1)}\cdot\textup{polylog}\left(n,\|C\|_{F},\frac{1}{\epsilon}\right)\right)$
additional gates. To summarize, in the absence of QRAM, the number of oracle
accesses is a factor $\sqrt{s}$ larger due to the Hamiltonian simulation, and
the gate complexity increases by a factor $n$ due to the cost of constructing
$O_{D}$ without QRAM.

We conclude this section by establishing the costs of preparing a block-
encoding of the final solution, and estimating trace inner products of the
form $\operatorname{tr}(A\tilde{\rho})$ for a given matrix $A$.

######  Proposition 12.

Suppose that Algorithm 3 is run with
$\zeta=\left(\frac{\epsilon}{n\|C\|_{F}}\right)^{4}$ for some
$\epsilon\in(0,1)$, and terminates after $K$ iterations, classically
outputting the tuples $\\{(\eta^{(k)},y^{(k)},q^{(k)},d^{(k)})\\}_{k=0}^{K}$.
Then, letting $C\|C\|_{F}^{-1}$ be stored in QRAM, and denoting the refining
problem at iteration $k$ by $\rho^{(k)}$, one can use
$\\{(\eta^{(k)},y^{(k)},q^{(k)},d^{(k)})\\}_{k=0}^{K}$ to implement an
$(n,{\cal O}(\log(n),\theta)$-block-encoding of

$\tilde{\rho}=\frac{1}{1+n\delta}\left[\left(\sum_{k=0}^{K}\frac{1}{\eta^{(k)}}Q^{(k)}\circ\frac{\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)}{\operatorname{tr}\left(\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)\right)}\right)+\delta
I\right],$

with at most $\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon},\frac{1}{\theta}}\left(\sqrt{n}\right)$
queries to the QRAM and $\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon},\frac{1}{\theta}}\left(n\right)$ classical
operations.

###### Proof.

First, note that

$\frac{\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)}{\operatorname{tr}\left(\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)\right)}=n^{-1}I$

whenever $y=(0,0)^{\top}$. Thus, by choosing $y^{(K+1)}=(0,0)^{\top}$,
$\eta^{(K+1)}=\frac{1}{n\delta}$, and $Q^{(K+1)}=ee^{\top}$ we can simplify
the expression of the final solution to

$\frac{1}{1+n\delta}\left[\sum_{k=0}^{K+1}\frac{1}{\eta^{(k)}}Q^{(k)}\circ\frac{\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)}{\operatorname{tr}\left(\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)\right)}\right].$

To ensure that the stated complexity holds, for each $k\in[K+1]$, we block-
encode

$A^{(k)}=\frac{Q^{(k)}\circ C}{\|C\|_{F}}+D^{(k)}.$

First, note that with classical access to $C$ and $q^{(k)}$, one can store
$\frac{Q^{(k)}\circ C}{\|C\|_{F}}$ in the QRAM by properly updating
$C\|C\|_{F}^{-1}$ in the QRAM. This step requires ${\cal O}(n)$ classical
operations, as the only non-trivial computation that is performed is limited
to the diagonal elements of the involved matrices. Then, with
$\frac{Q^{(k)}\circ C}{\|C\|_{F}}$ stored in QRAM, noting that
$\left\|\frac{Q^{(k)}\circ C}{\|C\|_{F}}\right\|_{F}\leq 1$ holds for every
$k\in[K+1]$, we apply Lemma 3 to construct a $(1,\log(n)+2,\theta_{1})$-block-
encoding of $\frac{Q^{(k)}\circ C}{\|C\|_{F}}$ in time ${\cal
O}\left(\textup{polylog}\left(\frac{n}{\theta_{1}}\right)\right)$. Similarly,
as we saw in the proof of Proposition 7, classical access to $d^{(k)}$ and
access to QRAM implies one can implement a $(1,\log(n)+3,\theta_{1})$-block-
encoding of $D^{(k)}$ can be constructed in time $\widetilde{{\cal
O}}_{\frac{n}{\theta_{1}}}(1)$.

Again following the proof of Proposition 7, applying Corollary 3 with
$y^{(k)}$ satisfying $\|y^{(k)}\|_{1}=\widetilde{{\cal O}}_{n}(\xi^{-1})$
implies that we can construct a unitary which prepares a copy of the Gibbs
state $\rho^{(k)}$ encoding the solution to the refining problem at iteration
$k$ with at most

$\widetilde{{\cal
O}}_{\frac{n}{\xi}}\left(\sqrt{n}\alpha\xi^{-1}\right)=\widetilde{{\cal
O}}_{n}\left(\sqrt{n}\right),$

accesses to the QRAM, as $\alpha=1$ and $\xi$ is a fixed constant. Therefore,
by Lemma 7, preparing a $(1,\log(n)+a,\theta_{1})$ block-encoding of a
purification of $\rho^{(k)}$ thus requires $\widetilde{{\cal
O}}_{\frac{n}{\theta_{1}}}(\sqrt{n})$ queries to the QRAM.

Next, provided classical access to the vector $q^{(k)}$ that store the
diagonal elements of $Q^{(k)}$, access to QRAM implies that we can efficiently
implement an oracle $O_{Q^{(k)}}$ that returns the entries of $Q^{(k)}$ in a
binary description:

$O_{Q^{(k)}}:\left|{i}\right\rangle\left|{j}\right\rangle\left|{0}\right\rangle^{\otimes
p}\mapsto\left|{i}\right\rangle\left|{j}\right\rangle\left|{q^{(k)}_{ij}}\right\rangle,\quad\forall
i,j\in[2^{\log n}]-1,$

where $q^{(k)}_{ij}$ is a $p$-bit binary description of the $ij$-matrix
element of $Q^{(k)}$ for $k=0,\dots,K+1$. By construction each matrix
$Q^{(k)}$ may be fully dense, and hence an application of Lemma 4 with
$s_{r}=s_{c}=n$ asserts that in the presence of QRAM, one can construct a
$(n,\log(n)+3,\theta_{2})$-block-encoding of $Q^{(k)}$ in time
$\widetilde{{\cal O}}_{\frac{n}{\theta_{2}}}(1)$.

From here, we can utilize Proposition 4 with
$\theta_{1}=\theta_{2}=\frac{\tilde{\theta}}{10}$ to construct an
$(n,a+4\log(n^{2})+12,\tilde{\theta})$-block-encoding of
$Q^{(k)}\circ\rho^{(k)}$ in time $\widetilde{{\cal
O}}_{\frac{n}{\tilde{\theta}}}(1)$. Repeating the above steps for
$k=0,\dots,K+1$, it follows that we can block-encode each of the terms
$Q^{(k)}\circ\rho^{(k)}$ using at most

$\widetilde{{\cal
O}}_{n,\frac{1}{\tilde{\theta}}}\left(K\sqrt{n}\right)=\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon},\frac{1}{\tilde{\theta}}}\left(\sqrt{n}\right)$

queries to the QRAM and $\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon},\frac{1}{\tilde{\theta}}}\left(n\right)$
classical operations, as $K={\cal
O}\left(\textup{polylog}\left(n,\|C\|_{F},\frac{1}{\epsilon}\right)\right)$ by
Corollary 5.

Finally, what remains is to take the linear combination of these terms. To do
so, we choose our weights to be $w_{k}=\frac{1}{2(1+n\delta)\eta^{(k)}}$,
which indeed satisfies $\|w\|_{1}\leq 1$. Then, we can construct a
$(K+2,\log(K+2),0)$-state-preparation pair $P_{L}$, $P_{R}$ for $w$, which can
be constructed by taking a $\log(K+2)$-fold tensor product of the Hadamard
gate, i.e.,

$P_{L}=P_{R}=\frac{1}{\sqrt{2}}\begin{pmatrix}1&1\\\
1&-1\end{pmatrix}^{\otimes\log(K+2)}.$

We are now in a position to apply Proposition 1, and choosing
$\tilde{\theta}=\frac{\theta}{n}$, we can obtain $W$ upon adding a control
qubit to the circuits used to construct the block-encoding of each
$Q^{(k)}\circ\rho^{(k)}$. As a result, we obtain an $(n,{\cal
O}(\log(n),\theta)$-block-encoding of $\tilde{\rho}$ with a single use of
$W,P_{R}$ and $P_{L}^{\dagger}$. Summing the cost of each step in the
construction we arrive at total cost of

$\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon},\frac{1}{\theta}}\left(\sqrt{n}\right)$

queries to the QRAM and $\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon},\frac{1}{\theta}}\left(n\right)$ classical
operations, and proof is complete. ∎

######  Proposition 13.

Suppose that Algorithm 3 is run with
$\zeta=\left(\frac{\epsilon}{n\|C\|_{F}}\right)^{4}$ for some
$\epsilon\in(0,1)$, and terminates after $K$ iterations, classically
outputting the tuples $\\{(\eta^{(k)},y^{(k)},q^{(k)},d^{(k)})\\}_{k=0}^{K}$.
Let $A\in\mathbb{R}^{n\times n}$ be a matrix with $\|A\|_{F}\leq 1$ and assume
classical access to $A$ and $C/\|C\|_{F}$. Then, with access to QRAM, one can
compute a $\theta$-precise estimate of $\operatorname{tr}(A\tilde{\rho})$
using at most $\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}\left(\frac{\sqrt{n}}{\theta}\right)$
queries to the QRAM and $\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}\left(ns\right)$ classical operations.

###### Proof.

See the proof of Theorem 8 in Appendix B. ∎

A QRAM-free version of Proposition 13 is also analyzed in Appendix B, and the
cost is summarized in Corollary 7. Without access to QRAM, the cost increases
with respect to $n$ because computing the Hadamard product of block-encodings
introduces $n$ as a subnormalization factor. This is compounded in the running
time, upon noting that we then have to scale down the error for the amplitude
estimation steps by $n$, and constructing sparse-access oracles for the
intermediate block-encodings of $Q$ and $D$ that arise in the trace estimation
procedure requires $\widetilde{{\cal O}}_{n}(n)$ gates.

###  5.3 Comparison to existing SDO algorithms

Table 1 presents a comparison of the running time results for the algorithms
we have proposed with the running times of the best performing methods from
both the classical and quantum literature when applied to solving (3).

Note that when directly solving (3), $m=n$, and any feasible solution $X$ to
(3) satisfies ${\operatorname{tr}{(X)}=n}$, implying $R=n$ for the algorithms
based on the (Q)MMWU framework. We also point out that the running times in
Table 1 take into account the role of sparsity in context of the algorithms,
which is measured as the maximum number of nonzero entries per row of the
constraint matrices $A_{1},\dots,A_{n}$. When using either an IPM or CPM to
solve (3), the $n$ constraint matrices are $A_{i}=e_{i}e_{i}^{\top}$ (with row
sparsity one) enforcing $X_{ii}=1$ for each diagonal element. On the other
hand, algorithms based on the (Q)MMWU or HU frameworks solve (3) by reducing
the problem to a feasibility problem; $C$ enters into the resulting
formulation as another constraint matrix, and as a result, the relevant
sparsity parameter is the maximum number of non-zeroes per row of $C$, which
we denote by $s$ in Table 1.

References | Method | Runtime | Error Scaling
---|---|---|---
[34] | IPM | $\widetilde{{\cal O}}_{n,\frac{1}{\epsilon}}\left(n^{\omega+0.5}\right)$ | $\epsilon$
[6] | QIPM | $\widetilde{{\cal O}}_{n,\kappa,\frac{1}{\epsilon}}\left(\sqrt{n}(n^{3}\kappa\epsilon^{-1}+n^{4})\right)$ | $\epsilon$
[41] | MMWU | $\widetilde{{\cal O}}_{n,\frac{1}{\epsilon}}\left(ns\epsilon^{-3.5}\right)$ | $\|C\|_{\ell_{1}}\epsilon$
[60] | QMMWU | $\widetilde{{\cal O}}_{n,\frac{1}{\epsilon}}\left(n^{5.5}s\epsilon^{-4}\right)$ | $n\|C\|\epsilon$
[10] (Classical) | HU | $\widetilde{{\cal O}}_{n,\|C\|}\left(\min\\{n^{2}s,n^{\omega}\\}\epsilon^{-12}\right)$ | $n\|C\|\epsilon$
[10] (Quantum) | HU | $\widetilde{{\cal O}}_{n,\|C\|,\frac{1}{\epsilon}}\left(n^{2.5}s^{0.5+o(1)}\epsilon^{-28+o(1)}\exp\left(1.6\sqrt{\log(\epsilon^{-1})}\right)\right)$ | $n\|C\|\epsilon$
[10] (Quantum) | HU-QRAM | $\widetilde{{\cal O}}_{n,\|C\|,\frac{1}{\epsilon}}\left(n^{1.5}s^{0.5+o(1)}\epsilon^{-28+o(1)}\exp\left(1.6\sqrt{\log(\epsilon^{-1})}\right)\right)$ | $n\|C\|\epsilon$
This work (Classical) | IR-HU | $\widetilde{{\cal O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}\left(\min\\{n^{2}s,n^{\omega}\\}\right)$ | $\epsilon$
This work (Quantum) | IR-HU | $\widetilde{{\cal O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}\left(n^{2.5}s^{0.5+o(1)}\right)$ | $\epsilon$
This work (Quantum) | IR-HU-QRAM | $\widetilde{{\cal O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}(n^{1.5})+ns$ | $\epsilon$
Table 1: Total running times for classical and quantum algorithms to solve
(3).

There are additional considerations that need to be taken into account when
making comparisons across methodologies listed in Table 1. Broadly speaking,
both (Q)MMWUs and HU require normalizing the problem by an upper bound on the
trace of a primal solution, and in the case of (3), we have the natural bound
$\operatorname{tr}(X)=n$. Moreover, (Q)MMWUs and HU additionally normalize the
cost matrix so that it exhibits unit norm with respect to some norm. While
these modifications amount to scaling the optimal objective value of (3) by a
fixed quantity, without employing any safeguards such as IR, these
modifications impact the scaling of the error as reflected in the fourth
column of Table 1. On the contrary, (Q)IPMs do not require the SDO problem to
be normalized in any way. Finally there is a distinction with regard to
output; (Q)IPMs explicitly report a classical description of the solution $X$,
whereas only the classical HU algorithm of [10] and our own classical IR-HU
method do so; the primal QMMWU of [60] reports a state-preparation pair $y$,
and the MMWU algorithm found in [41] reports a “gradient” $G\in{\cal S}^{n}$
such that $X=W\exp(G)W$ for a diagonal matrix $W$. As we noted earlier,
(Q)IPMs and (Q)MMWUs also utilize different definitions of optimality.

It can be easily seen that both the classical and quantum implementations of
our proposed methodology outperform all existing algorithms that exhibit poly-
logarithmic dependence on the precision $\epsilon$. Our classical algorithm is
only outperformed with respect to dimension by our own quantum algorithms, and
the algorithm from [41], which has an exponentially worse dependence on the
inverse prevision. Moreover, to achieve the same error scaling as our
algorithms, the algorithm from [41] would require time $\widetilde{{\cal
O}}_{n,\frac{1}{\epsilon}}\left(n^{4.5}s\epsilon^{-3.5}\right)$ (as it is
assumed $\|C\|_{\ell_{1}}=n$ in [41]). Up to poly-logarithmic factors, our
quantum algorithms outperform each of the classical and quantum solvers in
every parameter, suggesting the first evidence of quantum advantage for
solving a special class of SDO problems. Moreover, our implementation with
access to QRAM dominates all other algorithms. We therefore conclude that our
proposed algorithms are respectively, the fastest both in the classical and
quantum regimes.

##  6 Conclusion

In this work we devised an iterative refinement scheme for a particular class
of semidefinite optimization problems. The key to our idea behind our speedup
is to solve a sequence of related SDO problems in fixed low precision, rather
than solve one SDO problem using high accuracy requirements. Moreover, our
solutions satisfy a far stronger approximation guarantee over previous quantum
solution methodologies for this class of problem. We show that, provided
access to QRAM, a quantum implementation of our algorithm can produce accurate
solutions to SDO approximations of QUBO problems in time ${\cal
O}\left(ns+n^{1.5}\cdot\textup{polylog}\left(n,\|C\|_{F},\frac{1}{\epsilon}\right)\right)$
in the worst case. In the absence of QRAM, one can bound the running time of
the quantum algorithm using using the sparse-access input model, in which case
the algorithm exhibits an oracle complexity of ${\cal
O}\left(n^{2.5}s^{0.5+o(1)}\cdot\textup{polylog}\left(n,\|C\|_{F},\frac{1}{\epsilon}\right)\right)$.
A classical implementation of the algorithm exhibits worst case running time
of ${\cal
O}\left(\min\\{n^{2}s,n^{\omega}\\}\cdot\textup{polylog}\left(n,\|C\|_{F},\frac{1}{\epsilon}\right)\right)$,
which is at least a $\sqrt{n}$ factor better than classical IPMs.

When compared to the best performing algorithms in the literature, our
algorithms are the fastest in both the quantum and classical regimes,
respectively. This work indicates that there could be a genuine quantum
advantage (in the QRAM model) for this specific class of SDO problems; to
establish such an advantage, one would have to show that no classical
algorithm can beat the quantum running time. At the moment, we can only make
the weaker claim that our quantum algorithm is faster than any currently known
classical algorithm. We believe one can improve the theoretical performance of
our classical algorithm by not explicitly computing the density operator in
our subroutines. In particular, it may be possible to construct the separation
oracles as we do in the quantum setting using techniques to classically
estimate trace inner products of the form $\operatorname{tr}(A\rho)$ (see,
e.g., Appendix A in [61]), and applying ideas developed in [3, 41] to estimate
the diagonal elements of matrix exponentials via randomized projection [36].
It remains an open question as to whether our techniques can be applied to
general SDO problems using the matrix-multiplicative weights update framework
as a subroutine.

## Acknowledgements

This project has been carried out thanks to funding by the Defense Advanced
Research Projects Agency (DARPA), ONISQ grant W911NF2010022, titled The
Quantum Computing Revolution and Optimization: Challenges and Opportunities.

Giacomo Nannicini is partially supported by the Army Research Office under
grant number W911NF-20-1-0014.

##  Appendix A Running time of Algorithm 3 without QRAM

The following result from [10] gives the sample complexity of implementing the
oracles in the sparse-access model.

######  Lemma 16 (see, proof of Lemma 3.3 in [10]).

We can implement the oracle $O_{{\cal C}_{\gamma}}$ on a quantum computer
given access to ${\cal O}(\epsilon^{-2})$ copies of a state that is an
$\frac{\epsilon}{8}$-approximation of the input state $\rho$ in trace
distance. The oracle $O_{{\cal D}_{n}}$ can be implemented using ${\cal
O}(n\epsilon^{-2})$ $\frac{\epsilon}{8}$-approximate copies of the input, and
the classical post-processing time needed to implement the oracle is ${\cal
O}(n\epsilon^{-2})$.

Next, we bound the overall complexity of Algorithm 1 without access to QRAM.

######  Proposition 14.

Suppose that $C\in{\cal S}^{n}$ has row sparsity $s$ and $\xi\in(0,1)$. Then,
in the sparse-access input model, the complexity of solving (5) up to additive
error $\xi$ using Algorithm 1 on a quantum computer requires

$\widetilde{{\cal
O}}_{n}\left(n^{1.5}\sqrt{s}^{1+o(1)}\xi^{-7+o(1)}\exp\left(1.6\sqrt{\log(\xi^{-1})}\right)\right)$

queries to the input oracle $O_{C}$ and $\widetilde{{\cal
O}}_{n}\left(n^{2.5}\sqrt{s}^{1+o(1)}\xi^{-7+o(1)}\exp\left(1.6\sqrt{\log(\xi^{-1})}\right)\right)$
additional gates.

###### Proof.

Our proof can be viewed as the QRAM-free analogue of the discussion found in
[10, Section 3.4], and we repeat it here for completeness. In order to derive
an appropriate bound on the per-iteration cost, we need to evaluate the cost
of constructing our separation oracles. By Lemma 16, we can conclude that the
time to construct the oracle $O_{{\cal D}_{n}}$ for the diagonal elements
dominates that of constructing the oracle $O_{{\cal C}_{\gamma}}$ to test the
objective value.

We now turn our attention to the cost of simulating our Hamiltonian $H$. From
the results in [53, Appendix] it follows that we can produce a state that is
$\frac{\xi}{8}$ close to $\rho$ using $\widetilde{{\cal O}}(\sqrt{n}\xi^{-3})$
invocations of a controlled $U$ which satisfies

$\left\|U-e^{it_{0}H}\right\|\leq{\cal O}\left(\xi^{3}\right),$

with $t_{0}=\frac{\pi}{4\|H\|}$. Further, the authors in [10] note that each
of the Hamiltonians we seek to simulate are of the form
$H=y_{1}C\|C\|_{F}^{-1}+y_{2}D$ where $y_{1},y_{2}={\cal O}(\log(n)\xi^{-1})$
and $D$ is a diagonal matrix which satisfies $\|D\|\leq 1$. Invoking [15,
Theorem 1], we can simulate $H$ for time $t$ up to error $\xi^{3}$ using

$\widetilde{{\cal
O}}\left(t(a+b)\exp\left(1.6\sqrt{\log\left(\log(n)t\xi^{-3}\right)}\right)\right)$

separate simulations of $y_{1}C\|C\|_{F}$ and $y_{2}D$.

As noted in [10], access to the oracles $O_{\text{sparse}}$ and $O_{C}$ we
described in Section 2.1.1 allows us to simulate $\exp(itC\|C\|_{F}^{-1})$ in
time ${\cal O}\left((t\sqrt{s})^{1+o(1)}\xi^{o(1)}\right)$ if we utilize the
algorithm in [44]. Similarly, we follow [10] in constructing an oracle $O_{D}$
acting on $\mathbb{C}\otimes(\mathbb{C}^{2})^{\otimes a}$, where $a$ is a
sufficiently large constant such that we can represent the diagonal elements
of $D$ as

$O_{D}\left|{i,z}\right\rangle\mapsto\left|{i,z\oplus D_{ii}}\right\rangle$

to the desired level of precision in binary. Accordingly, we can simulate
$e^{iDt}$ for $t=\widetilde{{\cal O}}(\xi^{-1})$ using $\widetilde{{\cal
O}}_{n}(1)$ queries to $O_{D}$ and $\widetilde{{\cal O}}_{n}(1)$ elementary
operations [7], and we can implement $O_{D}$ using $\widetilde{{\cal
O}}_{n}(n)$ gates.

To summarize, the Gibbs sampler from [53] requires $\widetilde{{\cal
O}}(\sqrt{n}\xi^{-3})$ Hamiltonian simulation steps, each of which requires
time

$\widetilde{{\cal
O}}\left(\sqrt{s}^{1+o(1)}\xi^{o(1)}\exp\left(1.6\sqrt{\log(\xi^{-1})}\right)\right).$

Hence, each iteration of Algorithm 1 requires a total of

$\widetilde{{\cal
O}}_{n}\left(n^{1.5}\sqrt{s}^{1+o(1)}\xi^{-5+o(1)}\exp\left(1.6\sqrt{\log(\xi^{-1})}\right)\right)$

sparse-access oracle queries. Combining the above per-iteration cost with the
iteration bound ${\cal O}(\log(n)\xi^{-2})$ provided in Theorem 3, it follows
that Algorithm 1 solves (5) up to additive error $\xi$ with at most

$\widetilde{{\cal
O}}_{n}\left(n^{1.5}\sqrt{s}^{1+o(1)}\xi^{-7+o(1)}\exp\left(1.6\sqrt{\log(\xi^{-1})}\right)\right)$

queries to the input oracle $O_{C}$ and $\widetilde{{\cal
O}}_{n}\left(n^{2.5}\sqrt{s}^{1+o(1)}\xi^{-7+o(1)}\exp\left(1.6\sqrt{\log(\xi^{-1})}\right)\right)$
additional gates. ∎

Theorem 7 formalizes the complexity of of Algorithm 3 in the quantum setting
without access to QRAM. In our analysis, we employ the same Hamiltonian
simulation subroutines and Gibbs sampler used in [10] to construct our
separation oracles.

######  Theorem 7.

Let $C\in{\cal S}^{n}$ with row sparsity $s$ and $\epsilon\in(0,1)$. Then,
setting $\zeta=\left(\frac{\epsilon}{n\|C\|_{F}}\right)^{4}$ and fixing
$\xi=10^{-2}$, a quantum implementation of Algorithm 3 using the sparse-access
input model solves (3) up to additive error ${\cal O}(\epsilon)$ using

${\cal
O}\left(n^{1.5}s^{0.5+o(1)}\cdot\textup{polylog}\left(n,\|C\|_{F},\frac{1}{\epsilon}\right)\right)$

queries to the input oracle $O_{C}$ and ${\cal
O}\left(n^{2.5}s^{0.5+o(1)}\cdot\textup{polylog}\left(n,\|C\|_{F},\frac{1}{\epsilon}\right)\right)$
additional gates.

The output of the algorithm is a collection of tuples
$\\{(\eta^{(k)},y^{(k)},q^{(k)},d^{(k)})\\}_{k=0}^{K}$ such that

$\tilde{\rho}=\frac{1}{1+n\delta}\left[\left(\sum_{k=0}^{K}\frac{1}{\eta^{(k)}}Q^{(k)}\circ\frac{\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)}{\operatorname{tr}\left(\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)\right)}\right)+\delta
I\right],$

is a $2\zeta$-precise solution to (5), where $\delta$ is defined according to
(18). The entries of $\tilde{\rho}$ can be modified to construct a matrix
$\rho^{*}$ at trace distance ${\cal
O}\left(\frac{\epsilon}{n\|C\|_{F}}\right)$ of $\tilde{\rho}$ in time ${\cal
O}(n^{2})$, such that $n\rho^{*}$ is a feasible point of the SDO problem (3).

###### Proof.

Given that $C$ is an $s$-sparse matrix, we can load $C$ in ${\cal O}(ns)$
time. Similarly, for normalization purposes we classically compute
$\|C\|_{F}$, which requires ${\cal O}(ns)$ arithmetic operations. In each
iteration we use Algorithm 1 to solve (14), and use classical estimates of the
diagonal elements of the refining solution, and a classical estimate of the
objective value attained by the refining solution to update the solution and
data for the refining problem we need to solve in the next iteration.

Letting ${\cal T}_{HU}^{\text{sparse}}$ be the cost of using Algorithm 1 as an
approximate SDO subroutine, we saw in Proposition 14, Algorithm 1 solves (14)
to additive error $\xi$ using

${\cal T}_{HU}^{\text{sparse}}=\widetilde{{\cal
O}}_{n}\left(n^{1.5}\sqrt{s}^{1+o(1)}\xi^{-7+o(1)}\exp\left(1.6\sqrt{\log(\xi^{-1})}\right)\right)$

queries to the oracle describing the problem data and $\widetilde{{\cal
O}}_{n}\left(n^{2.5}\sqrt{s}^{1+o(1)}\xi^{-7+o(1)}\exp\left(1.6\sqrt{\log(\xi^{-1})}\right)\right)$
additional gates. In the context of Algorithm 3, $\xi$ is a fixed constant, so
the cost of our oracle call to Algorithm 1 simplifies to

${\cal T}_{HU}^{\text{sparse}}=\widetilde{{\cal
O}}_{n}\left(n^{1.5}\sqrt{s}^{1+o(1)}\right)$

queries to the oracle describing the problem data and $\widetilde{{\cal
O}}_{n}\left(n^{2.5}\sqrt{s}^{1+o(1)}\right)$ additional gates.

Classically updating the objective value requires ${\cal O}(1)$ arithmetic
operations while updating the vector $p$ which stores a classical description
of the diagonal elements of our solution as

$p_{i}\leftarrow p_{i}+\frac{Q_{ii}}{\eta^{(k)}}\tilde{p}_{i}^{(k)}$

requires ${\cal O}(n)$ arithmetic operations. Again, $\varepsilon$ and $Q$ can
each be updated using ${\cal O}(n)$ arithmetic operations, as we only need to
store the diagonal elements of $Q$. This also implies that we can also
calculate $Q\circ\frac{C}{\|C\|_{F}}$ in time ${\cal O}(n)$, for only the
element-wise products along the diagonal are non-trivial. When compared to
loading and normalizing the data or our use of Algorithm 1 as a subroutine for
solving (14), these intermediate computation steps are negligible and do not
factor into the overall running time using ${\cal O}$ notation.

Factoring in the ${\cal
O}\left(\text{polylog}\left(\frac{1}{\zeta}\right)\right)={\cal
O}\left(\text{polylog}\left(n,\|C\|_{F},\frac{1}{\epsilon}\right)\right)$ from
Corollary 5, it follows that a quantum implementation of Algorithm 3 requires
at most

${\cal
O}\left(n^{1.5}s^{0.5+o(1)}\cdot\textup{polylog}\left(n,\|C\|_{F},\frac{1}{\epsilon}\right)\right)$

queries to the input oracle $O_{C}$ and ${\cal
O}\left(n^{2.5}s^{0.5+o(1)}\cdot\textup{polylog}\left(n,\|C\|_{F},\frac{1}{\epsilon}\right)\right)$
additional gates. Just as in the proof of Theorem 5, applying Proposition 11
with our choice of $\zeta=\left(\frac{\epsilon}{n\|C\|_{F}}\right)^{4}$
implies that the above running time is sufficient to obtain a solution that
can be used to solve (3) up to additive error ${\cal O}(\epsilon)$, and the
proof is complete. ∎

##  Appendix B Estimating trace inner products with the final solution

Given that we do not explicitly report a classical description of the final
solution $\tilde{\rho}$ defined in equation (23), it may be of interest to
understand how, for a user specified matrix $A$, one can compute the trace
inner product $\operatorname{tr}(A\tilde{\rho})$. We outline a procedure for
doing so using the state preparation pair description of solution
$\\{(\eta^{(k)},y^{(k)},q^{(k)},d^{(k)})\\}_{k=0}^{K}$ in Algorithm 4, and
subsequently analyze the complexity of doing so.

Input: Access to an $s$-sparse matrix $A\in\mathbb{R}^{n\times n}$ with
$\|A\|_{F}\leq 1$, state preparation pair description of solution
$\\{(\eta^{(k)},y^{(k)},q^{(k)},d^{(k)})\\}_{k=0}^{K}$, precision
$\theta\in(0,1)$, $\zeta=\left(\frac{\epsilon}{n\|C\|_{F}}\right)^{4}$

Output: A $\theta$-precise classical estimate of
$\operatorname{tr}(A\tilde{\rho})$

Initialize: $a\leftarrow 0$, $k\leftarrow 0$,
$y^{(K+1)}\leftarrow(0,0)^{\top}$, $\eta^{(K+1)}\leftarrow\frac{1}{n\delta}$,
$Q^{(K+1)}\leftarrow ee^{\top}$

for  _$k=0,\dots,K+1$_ do

  1. 1.

Implement an $(\alpha,a,\zeta/2(K+2))$-block-encoding of $Q^{(k)}\circ A$

  2. 2.

Use block-encoding of $Q^{(k)}\circ A$ to implement a trace estimator for

$a^{(k)}=\operatorname{tr}\left[\left(Q^{(k)}\circ
A\right)\left(\frac{\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)}{\operatorname{tr}\left(\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)\right)}\right)\right]$

  3. 3.

Use ${\cal O}\left(\frac{K}{\theta}\right)$ samples from the trace estimator
to produce $\frac{\theta}{K+2}$-precise estimate $\tilde{a}^{(k)}$ of
$a^{(k)}$

  4. 4.

Update solution:

$a\leftarrow a+\frac{1}{\eta^{(k)}}\tilde{a}^{(k)}$

  5. 5.

$k\leftarrow k+1$

end for

Scale down estimate to account for spectrum shift:
$a\leftarrow\frac{1}{1+n\delta}a$

Algorithm 4 Trace estimation procedure for the final solution

######  Theorem 8.

Let $A\in\mathbb{R}^{n\times n}$, and $\frac{C}{\|C\|_{F}}\in{\cal S}^{n}$ be
stored in QRAM, $\theta\in(0,1)$, and
$\\{(\eta^{(k)},y^{(k)},q^{(k)},d^{(k)})\\}_{k=0}^{K}$ be a state preparation
pair description of the solution obtained from running Algorithm 3 to final
precision $\zeta=\left(\frac{\epsilon}{n\|C\|_{F}}\right)^{4}$. Suppose $A$ is
an $s$-sparse matrix with $\|A\|_{F}\leq 1$, and assume classical access to
$A$ and $\frac{C}{\|C\|_{F}}\in{\cal S}^{n}$. Then, Algorithm 4 outputs a
$\theta$-precise estimate of

$\operatorname{tr}(A\tilde{\rho})=\frac{1}{1+n\delta}\operatorname{tr}\left(A\left[\left(\sum_{k=0}^{K}\frac{1}{\eta^{(k)}}Q^{(k)}\circ\frac{\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)}{\operatorname{tr}\left(\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)\right)}\right)+\delta
I\right]\right)$

using at most

$\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}\left(\frac{\sqrt{n}}{\theta}\right)$

queries to the QRAM and $\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}\left(ns\right)$ classical operations.

###### Proof.

We begin by establishing the correctness of Algorithm 4. First, note that
following the proof of Proposition 12, we can simplify the expression of the
final solution to

$\tilde{\rho}=\frac{1}{1+n\delta}\left[\sum_{k=0}^{K+1}\frac{1}{\eta^{(k)}}Q^{(k)}\circ\frac{\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)}{\operatorname{tr}\left(\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)\right)}\right].$

by setting $y^{(K+1)}=(0,0)^{\top}$, $\eta^{(K+1)}=\frac{1}{n\delta}$, and
$Q^{(K+1)}=ee^{\top}$. Then, by linearity of the trace and Lemma 1, one has:

$\displaystyle\operatorname{tr}(A\tilde{\rho})$
$\displaystyle=\operatorname{tr}\left(A\frac{1}{1+n\delta}\left[\sum_{k=0}^{K+1}\frac{1}{\eta^{(k)}}Q^{(k)}\circ\frac{\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)}{\operatorname{tr}\left(\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)\right)}\right]\right)$
$\displaystyle=\frac{1}{1+n\delta}\sum_{k=0}^{K+1}\frac{1}{\eta^{(k)}}\operatorname{tr}\left(A\left[Q^{(k)}\circ\frac{\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)}{\operatorname{tr}\left(\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)\right)}\right]\right)$
$\displaystyle=\frac{1}{1+n\delta}\sum_{k=0}^{K+1}\frac{1}{\eta^{(k)}}\operatorname{tr}\left(\left(Q^{(k)}\circ
A\right)\frac{\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)}{\operatorname{tr}\left(\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)\right)}\right).$

In other words, the output of Algorithm 4 is indeed an estimate of
$\operatorname{tr}(A\tilde{\rho})$.

Next, we analyze the complexity of the procedure. If $A$ is classically known,
one can store $Q^{(k)}\circ A$ in the QRAM using ${\cal O}(ns)$ classical
operations, as $A$ is $s$-sparse. With $Q\circ A$ stored in a QRAM data
structure, one can apply Lemma 3 to implement an
$(1,\log(n)+2,\zeta/2(K+2))$-block-encoding of $Q\circ A$ in time
$\widetilde{{\cal O}}_{\frac{nK}{\zeta}}(1)$ (as $\|Q\circ
A\|_{F}\leq\|A\|_{F}\leq 1$ for any $Q$ defined according to (13)). As we saw
in the proof of Proposition 12, with $\frac{C}{\|C\|_{F}}$ stored in QRAM, one
can implement the state

$\rho^{(k)}=\frac{\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)}{\operatorname{tr}\left(\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)\right)}$

using at most

$\widetilde{{\cal O}}_{n}\left(\sqrt{n}\right),$

accesses to the QRAM and ${\cal O}(n)$ classical operations.

Having prepared the state $\rho^{(k)}$ and a
$(1,\log(n)+2,\zeta/2(K+2))$-block-encoding $U_{k}$ of $Q^{(k)}\circ A$, Lemma
8 asserts that one can implement a trace estimator for

$\operatorname{tr}\left[\left(Q^{(k)}\circ A\right)\rho^{(k)}\right]$

with bias at most $\frac{\zeta}{K+2}$ using $\widetilde{{\cal O}}(1)$
applications of $U_{k}$ and $U_{k}^{\dagger}$. Applying amplitude estimation
using ${\cal O}\left(\frac{K}{\theta}\right)=\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}\left(\frac{1}{\theta}\right)$ samples
from the estimator, we obtain a $\frac{\theta}{K+2}$-precise classical
estimate $\tilde{a}^{(k)}$ of $a^{(k)}$, as $K={\cal
O}\left(\textup{polylog}\left(n,\|C\|_{F},\frac{1}{\epsilon}\right)\right)$.

From here, we classically update $a$ using ${\cal O}(1)$ arithmetic
operations. Therefore, each iteration of Algorithm 4 requires at most

$\widetilde{{\cal O}}_{n,\frac{K}{\zeta}}\left(\frac{\sqrt{n}}{\theta}\right)$

accesses to the QRAM and ${\cal O}(ns)$ classical operations. Summing over
$K+2$ iterations implies a total of

$\widetilde{{\cal
O}}_{n,\frac{K}{\zeta}}\left(K\left(\frac{\sqrt{n}}{\theta}\right)\right)=\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}\left(\frac{\sqrt{n}}{\theta}\right)$

accesses to the QRAM and

${\cal O}\left(Kns\right)=\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}\left(ns\right)$

classical operations. The proof is complete. ∎

Note that if $\|A\|_{F}>1$, because of the subnormalization to block-encode
$A$ we need to increase precision of the estimation procedure: the cost
increases by a factor proportional to $\|A\|_{F}$.

######  Corollary 7.

Let $A\in\mathbb{R}^{n\times n}$, $\theta\in(0,1)$, and
$\\{(\eta^{(k)},y^{(k)},q^{(k)},d^{(k)})\\}_{k=0}^{K}$ be a state preparation
pair description of the solution obtained from running Algorithm 3 to final
precision $\zeta=\left(\frac{\epsilon}{n\|C\|_{F}}\right)^{4}$. Suppose $A$ is
an $s$-sparse matrix with $\|A\|_{F}\leq 1$, and assume sparse oracle access
to $A$ and $\frac{C}{\|C\|_{F}}\in{\cal S}^{n}$. Then, Algorithm 4 outputs a
$\theta$-precise estimate of

$\operatorname{tr}(A\tilde{\rho})=\frac{1}{1+n\delta}\operatorname{tr}\left(A\left[\left(\sum_{k=0}^{K}\frac{1}{\eta^{(k)}}Q^{(k)}\circ\frac{\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)}{\operatorname{tr}\left(\exp\left(-\left[y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}\operatorname{diag}\left(d^{(k)}\right)\right]\right)\right)}\right)+\delta
I\right]\right)$

using at most

$\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}\left(\frac{n^{2.5}s^{2}}{\theta}\right)$

queries to $O_{A}$, $O_{C}$, and $\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}\left(\frac{n^{3.5}s^{2}}{\theta}\right)$
additional gates.

###### Proof.

Provided classical access to $A$, we use Lemma 4 with $s_{r}=s_{c}$ to
construct an $(s,\log(n)+3,\theta/n)$-block-encoding of $A$ with two uses of
$O_{A}$ (an oracle describing the elements of $A$ in binary), and additionally
using $\widetilde{{\cal O}}_{n}\left(1\right)$ one and two qubit gates.

Likewise, with access to the oracle $O_{C}$ describing the elements of
$C\|C\|_{F}^{-1}$, one can construct an $(s,\log(n)+3,\theta/n)$-block-
encoding of $C\|C\|_{F}^{-1}$ with two uses of $O_{C}$, and additionally using
$\widetilde{{\cal O}}_{n}\left(1\right)$ one and two qubit gates. Note that
without access to QRAM, we must compute the Hadamard products by taking the
Hadamard products of block-encodings, which causes the subnormalization factor
for the Hadamard product $Q^{(k)}\circ C\|C\|_{F}^{-1}$ to be $ns$, as
$Q^{(k)}$ may be fully dense and $C$ is $s$-sparse. It follows that preparing
one copy of each Gibbs state requires

$\widetilde{{\cal O}}_{n}\left(\sqrt{n}(ns)\right)=\widetilde{{\cal
O}}_{n}\left(n^{1.5}s\right)$

accesses to block-encodings of $Q^{(k)}\circ C\|C\|_{F}^{-1}$ and $D$, which
each require an additional $\widetilde{{\cal O}}_{n}(n)$ gates (to construct
sparse-access oracles for $Q^{(k)}$ and $D$).

Similarly, the subnormalization factor for a block-encoding $U_{k}$ of
$Q^{(k)}\circ A$ will be $ns$. Having prepared the state $\rho^{(k)}$ and a
block-encoding $Q^{(k)}\circ A$, Lemma 8 asserts that one can implement a
trace estimator for

$\operatorname{tr}\left[\left(Q^{(k)}\circ A\right)\rho^{(k)}\right]$

with bias at most $\frac{\zeta}{K+2}$ using $\widetilde{{\cal O}}(ns)$
applications of $U_{k}$ and $U_{k}^{\dagger}$. Applying amplitude estimation
using ${\cal O}\left(\frac{K}{\theta}\right)=\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}\left(\frac{1}{\theta}\right)$ samples
from the estimator to obtain a $\frac{\theta}{K+2}$-precise classical estimate
$\tilde{a}^{(k)}$ of $a^{(k)}$, as $K={\cal
O}\left(\textup{polylog}\left(n,\|C\|_{F},\frac{1}{\epsilon}\right)\right)$.

Just as in the QRAM setting, classically updating $a$ requires ${\cal O}(1)$
arithmetic operations. Therefore, without access to QRAM, each iteration of
Algorithm 4 requires at most

$\widetilde{{\cal
O}}_{n,\frac{K}{\zeta}}\left(\frac{n^{2.5}s^{2}}{\theta}\right)=\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}\left(\frac{n^{2.5}s^{2}}{\theta}\right)$

applications of block-encodings for $Q^{(k)}\circ C\|C\|_{F}^{-1}$, $D^{(k)}$
and $Q^{(k)}\circ A$ and $\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}\left(\frac{n^{3.5}s^{2}}{\theta}\right)$
additional gates. This corresponds to $\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}\left(\frac{n^{2.5}s^{2}}{\theta}\right)$
queries to $O_{A}$ and $O_{C}$ in each iteration, and $\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}\left(\frac{n^{3.5}s^{2}}{\theta}\right)$
additional gates. Summing over the $K+2=\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}(1)$ iterations yields the stated
complexity. ∎

## References

  * [1] Noga Alon, W. Fernandez De La Vega, Ravi Kannan, and Marek Karpinski.  Random sampling and approximation of MAX-CSP Problems.  Journal of Computer and System Sciences, 67(2):212–243, 2003.
  * [2] David L. Applegate, William Cook, Sanjeeb Dash, and Daniel G. Espinoza.  Exact solutions to linear programming problems.  Operations Research Letters, 35(6):693–699, 2007.
  * [3] Sanjeev Arora, Elad Hazan, and Satyen Kale.  The multiplicative weights method: a meta-algorithm and its applications.  Theory of Computing, 8(6) 121-164, 2012.
  * [4] Sanjeev Arora and Satyen Kale.  A combinatorial, primal-dual approach to semidefinite programs.  Journal of the ACM (JACM), 63(2):1–35, 2016.
  * [5] Srinivasan Arunachalam, Vlad Gheorghiu, Tomas Jochym-O’Connor, Michele Mosca, and Priyaa Varshinee Srinivasan.  On the robustness of bucket brigade quantum ram.  New Journal of Physics, 17(12):123010, 2015.
  * [6] Brandon Augustino, Giacomo Nannicini, Tamás Terlaky, and Luis F. Zuluaga.  Quantum interior point methods for semidefinite optimization.  arXiv preprint arXiv:2112.06025, 2021.
  * [7] Dominic W. Berry, Graeme Ahokas, Richard Cleve, and Barry C. Sanders.  Efficient quantum algorithms for simulating sparse Hamiltonians.  Communications in Mathematical Physics, 270(2):359–371, 2007.
  * [8] Rajendra Bhatia.  Matrix Analysis, volume 169.  Springer Science & Business Media, 2013.
  * [9] Fernando G.S.L. Brandão, Amir Kalev, Tongyang Li, Cedric Yen-Yu Lin, Krysta M. Svore, and Xiaodi Wu.  Quantum SDP Solvers: Large Speed-Ups, Optimality, and Applications to Quantum Learning.  In Christel Baier, Ioannis Chatzigiannakis, Paola Flocchini, and Stefano Leonardi, editors, 46th International Colloquium on Automata, Languages, and Programming (ICALP 2019), volume 132, pages 27:1–27:14, Dagstuhl, Germany, 2019. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik.
  * [10] Fernando G.S.L. Brandão, Richard Kueng, and Daniel Stilck França.  Faster quantum and classical SDP approximations for quadratic binary optimization.  Quantum, 6:625, 2022.
  * [11] Fernando G.S.L. Brandão and Krysta M. Svore.  Quantum speed-ups for solving semidefinite programs.  In Rafail Ostrovsky and Chris Umans, editors, 2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS), pages 415–426. IEEE, 2017.
  * [12] Shantanav Chakraborty, András Gilyén, and Stacey Jeffery.  The power of block-encoded matrix powers: improved regression techniques via faster Hamiltonian simulation.  In Christel Baier, Ioannis Chatzigiannakis, Paola Flocchini, and Stefano Leonardi, editors, 46th International Colloquium on Automata, Languages, and Programming (ICALP 2019), volume 132, pages 33:1–33:14, Dagstuhl, Germany, 2019. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik.
  * [13] Chen-Fu Chiang, Anirban Chowdhury, and Pawel Wocjan.  Space-efficient quantization method for reversible markov chains.  arXiv preprint arXiv:2206.06886, 2022.
  * [14] Andrew M. Childs, Robin Kothari, and Rolando D. Somma.  Quantum algorithm for systems of linear equations with exponentially improved dependence on precision.  SIAM Journal on Computing, 46(6):1920–1950, 2017.
  * [15] Andrew M. Childs and Nathan Wiebe.  Hamiltonian Simulation using linear combinations of unitary operations.  Quantum Information and Computation, 12(11–12):901–924, Nov 2012\.
  * [16] Anirban Narayan Chowdhury and Rolando D. Somma.  Quantum algorithms for Gibbs sampling and Hitting-Time estimation.  Quantum Information & Computing, 17(1–2):41–64, Feb 2017.
  * [17] William Cook, Thorsten Koch, Daniel E. Steffy, and Kati Wolter.  An exact rational mixed-integer programming solver.  In Oktay Günlük and Gerhard J. Woeginger, editors, International Conference on Integer Programming and Combinatorial Optimization, pages 104–116. Springer, 2011.
  * [18] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  A quantum approximate optimization algorithm.  arXiv preprint arXiv:1411.4028, 2014.
  * [19] Aleta Berk Finnila, Maria A. Gomez, C. Sebenik, Catherine Stenson, and Jimmie D. Doll.  Quantum annealing: A new method for minimizing multidimensional functions.  Chemical Physics Letters, 219(5-6):343–348, 1994.
  * [20] Daniel Stilck França.  Perfect sampling for quantum Gibbs states.  Quantum Information and Computation, 18:361–388, 2018.
  * [21] Alan Frieze and Ravi Kannan.  Quick approximation to matrices and applications.  Combinatorica, 19(2):175–220, 1999.
  * [22] András Gilyén.  Quantum singular value transformation & its algorithmic applications.  PhD thesis, University of Amsterdam, 2019.
  * [23] András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe.  Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics.  In Moses Charikar and Edith Cohen, editors, Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, pages 193–204, 2019\.
  * [24] Vittorio Giovannetti, Seth Lloyd, and Lorenzo Maccone.  Quantum random access memory.  Physical Review Letters, 100(16):160501, 2008.
  * [25] Ambros M. Gleixner and Daniel E. Steffy.  Linear programming using limited-precision oracles.  Mathematical Programming, 183(1):525–554, 2020.
  * [26] Ambros M. Gleixner, Daniel E. Steffy, and Kati Wolter.  Improving the accuracy of linear programming solvers with iterative refinement.  In Joris van der Hoeven and Mark van Hoeij, editors, Proceedings of the 37th International Symposium on Symbolic and Algebraic Computation, pages 187–194, 2012.
  * [27] Ambros M. Gleixner, Daniel E. Steffy, and Kati Wolter.  Iterative refinement for linear programming.  INFORMS Journal on Computing, 28(3):449–464, 2016.
  * [28] Michel X. Goemans and David P. Williamson.  Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming.  Journal of the ACM (JACM), 42(6):1115–1145, 1995.
  * [29] Gene H. Golub and Charles F. Van Loan.  Matrix Computations.  Johns Hopkins University Press, 2013.
  * [30] Sander Gribling.  Applications of optimization to factorization ranks and quantum information theory.  PhD thesis, Tilburg University, 2019.
  * [31] Aram W. Harrow, Avinatan Hassidim, and Seth Lloyd.  Quantum algorithm for linear systems of equations.  Physical Review Letters, 103(15):150502, 2009.
  * [32] Elad Hazan.  Introduction to online convex optimization.  Foundations and Trends® in Optimization, 2(3-4):157–325, 2016.
  * [33] Roger Horn and Charles R. Johnson.  Topics in Matrix Analysis.  Cambridge University Press Cambridge, UK, 1994.
  * [34] Haotian Jiang, Tarun Kathuria, Yin Tat Lee, Swati Padmanabhan, and Zhao Song.  A faster interior point method for semidefinite programming.  In Sandy Irani, Lisa O’Conner, and Patrick Kellenberger, editors, 2020 IEEE 61st Annual Symposium on Foundations of Computer Science (FOCS), pages 910–918. IEEE, 2020.
  * [35] Haotian Jiang, Yin Tat Lee, Zhao Song, and Sam Chiu-wai Wong.  An improved cutting plane method for convex optimization, convex-concave games, and its applications.  In Konstantin Makarychev, Yury Makarychev, Madhur Tulsiani, Gautam Kamath, and Julia Chuzhoy, editors, Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing, pages 944–953, 2020.
  * [36] William B. Johnson, Joram Lindenstrauss, and Gideon Schechtman.  Extensions of Lipschitz mappings into a Hilbert space.  Contemporary Mathematics, 26:189–206, 1984.
  * [37] Michael J. Kastoryano and Fernando G.S.L. Brandao.  Quantum Gibbs samplers: The commuting case.  Communications in Mathematical Physics, 344(3):915–957, 2016.
  * [38] Iordanis Kerenidis and Anupam Prakash.  Quantum gradient descent for linear systems and least squares.  Physical Review A, 101(2):022316, 2020.
  * [39] Iordanis Kerenidis and Anupam Prakash.  A quantum interior point method for LPs and SDPs.  ACM Transactions on Quantum Computing, 1(1):1–32, 2020.
  * [40] Christopher King.  Inequalities for trace norms of 2$\times$ 2 block matrices.  Communications in Mathematical Physics, 242(3):531–545, 2003.
  * [41] Yin Tat Lee and Swati Padmanabhan.  An $\widetilde{\mathcal{O}}(m/\varepsilon^{3.5})$-cost algorithm for semidefinite programs with diagonal constraints.  In Jacob Abernethy and Shivani Agarwal, editors, Conference on Learning Theory, pages 3069–3119. PMLR, 2020.
  * [42] Yin Tat Lee, Aaron Sidford, and Sam Chiu-wai Wong.  A faster cutting plane method and its implications for combinatorial and convex optimization.  In Rafail Ostrovsky and Venkatesan Guruswami, editors, 2015 IEEE 56th Annual Symposium on Foundations of Computer Science (FOCS), pages 1049–1065. IEEE, 2015.
  * [43] László Lovász.  On the Shannon capacity of a graph.  IEEE Transactions on Information Theory, 25(1):1–7, 1979.
  * [44] Guang Hao Low.  Hamiltonian simulation with nearly optimal dependence on spectral norm.  In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, pages 491–502, 2019.
  * [45] Guang Hao Low and Isaac L. Chuang.  Hamiltonian simulation by qubitization.  Quantum, 3:163, 2019.
  * [46] Renato D.C. Monteiro.  Polynomial convergence of primal-dual algorithms for semidefinite programming based on the Monteiro and Zhang family of directions.  SIAM Journal on Optimization, 8(3):797–812, 1998.
  * [47] Yurii E. Nesterov and Arkadi Nemirovskii.  A general approach to polynomial-time algorithms design for convex programming.  Report, Central Economical and Mathematical Institute, USSR Academy of Sciences, Moscow, 1988.
  * [48] Yurii E. Nesterov and Arkadi Nemirovskii.  Interior-Point Polynomial Algorithms in Convex Programming, volume 13.  SIAM, 1995.
  * [49] Yurii E. Nesterov and Michael J. Todd.  Self-scaled barriers and interior-point methods for convex programming.  Mathematics of Operations Research, 22(1):1–42, 1997.
  * [50] Yurii E. Nesterov and Michael J. Todd.  Primal-dual interior-point methods for self-scaled cones.  SIAM Journal on Optimization, 8(2):324–364, 1998.
  * [51] Michael A. Nielsen and Isaac Chuang.  Quantum Computation and Quantum Information.  American Association of Physics Teachers, 2002.
  * [52] Foad Mahdavi Pajouh, Balabhaskar Balasundaram, and Oleg A. Prokopyev.  On characterization of maximal independent sets via quadratic optimization.  Journal of Heuristics, 19(4):629–644, 2013.
  * [53] David Poulin and Pawel Wocjan.  Sampling from the thermal quantum Gibbs state and evaluating partition functions with a quantum computer.  Physical Review Letters, 103(22):220502, 2009.
  * [54] Maurice Sion.  On general minimax theorems.  Pacific Journal of Mathematics, 8(1):171–176, 1958.
  * [55] Jos F. Sturm.  Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones.  Optimization Methods and Software, 11(1-4):625–653, 1999.
  * [56] Kim-Chuan Toh, Michael J. Todd, and Reha H. Tütüncü.  SDPT3—a MATLAB software package for semidefinite programming, version 1.3.  Optimization Methods and Software, 11(1-4):545–581, 1999.
  * [57] Koji Tsuda, Gunnar Rätsch, and Manfred K. Warmuth.  Matrix exponentiated gradient updates for on-line learning and Bregman projection.  Journal of Machine Learning Research, 6(Jun):995–1018, 2005.
  * [58] Joran van Apeldoorn.  A quantum view on convex optimization.  PhD thesis, University of Amsterdam, February 2020.
  * [59] Joran van Apeldoorn.  Quantum probability oracles & multidimensional amplitude estimation.  In Min-Hsiu Hsieh, editor, 16th Conference on the Theory of Quantum Computation, Communication and Cryptography (TQC 2021). Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2021.
  * [60] Joran van Apeldoorn and András Gilyén.  Improvements in Quantum SDP-Solving with Applications.  In Christel Baier, Ioannis Chatzigiannakis, Paola Flocchini, and Stefano Leonardi, editors, 46th International Colloquium on Automata, Languages, and Programming (ICALP 2019), volume 132 of Leibniz International Proceedings in Informatics (LIPIcs), pages 99:1–99:15, Dagstuhl, Germany, 2019. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik.
  * [61] Joran van Apeldoorn, András Gilyén, Sander Gribling, and Ronald de Wolf.  Quantum SDP-solvers: Better upper and lower bounds.  Quantum, 4:230, 2020.
  * [62] James Hardy Wilkinson.  Rounding Errors in Algebraic Processes.  Courier Corporation, 1994.
  * [63] Man-Hong Yung and Alán Aspuru-Guzik.  A quantum–quantum Metropolis algorithm.  Proceedings of the National Academy of Sciences, 109(3):754–759, 2012.