# Solving the semidefinite relaxation of QUBOs in matrix multiplication time,
and faster with a quantum computer

Brandon Augustino Corresponding Author<EMAIL_ADDRESS>of
Industrial and Systems Engineering, Quantum Computing and Optimization Lab,
Lehigh University  Giacomo Nannicini Department of Industrial and Systems
Engineering, University of Southern California  $\text{Tam\'{a}s
Terlaky}^{\dagger}$ and $\text{Luis F. Zuluaga}^{\dagger}$

###### Abstract

Recent works on quantum algorithms for solving semidefinite optimization (SDO)
problems have leveraged a quantum-mechanical interpretation of positive
semidefinite matrices to develop methods that obtain quantum speedups with
respect to the dimension $n$ and number of constraints $m$. While their
dependence on other parameters suggests no overall speedup over classical
methodologies, some quantum SDO solvers provide speedups in the low-precision
regime. We exploit this fact to our advantage, and present an iterative
refinement scheme for the Hamiltonian Updates algorithm of Brandão et al.
(_Quantum_ 6, 625 (2022)) to exponentially improve the dependence of their
algorithm on the precision $\epsilon$, defined as the absolute gap between
primal and dual solution. As a result, we obtain a classical algorithm to
solve the semidefinite relaxation of Quadratic Unconstrained Binary
Optimization problems (QUBOs) in matrix multiplication time. Provided access
to a quantum read/classical write random access memory (QRAM), a quantum
implementation of our algorithm exhibits ${\cal
O}\left(ns+n^{1.5}\cdot\textup{polylog}\left(n,\|C\|_{F},\frac{1}{\epsilon}\right)\right)$
running time, where $C$ is the cost matrix and $s$ is its sparsity parameter
(maximum number of nonzero elements per row).

##  1 Introduction

We consider optimization problems of the form:

$\displaystyle\max\leavevmode\nobreak$ $\displaystyle x^{\top}Cx$ (1) s.t.
$\displaystyle x\in\\{-1,1\\}^{n},$

where $C\in{\cal S}^{n}$ is the problem data and ${\cal S}^{n}$ is the space
of symmetric matrices in $\mathbb{R}^{n\times n}$. Solving (1) can be viewed
as computing the $\infty\to 1$ norm of the coefficient matrix $C$. This
particular norm is intrinsically related the cut norm of a matrix, which plays
a crucial role in developing efficient approximation algorithms for dense
graph and matrix problems [1, 21], with perhaps the most well-known
application being the task of finding the largest cut in a graph (MaxCut).
These problems also play an important role in quantum information sciences;
the Ising model belongs to this class of problems [52], and quantum algorithms
such as the Quantum Approximate Optimization Algorithm (QAOA) [18] and quantum
annealing [19] can address its solution.

Computing the cut norm corresponds to replacing $x\in\\{-1,1\\}^{n}$ with
$z\in\\{0,1\\}^{n}$ in (1), giving rise to quadratic unconstrained binary
optimization (QUBO) problems. A standard QUBO is of the form

$\displaystyle\max\leavevmode\nobreak$ $\displaystyle z^{\top}Cz$ (2) s.t.
$\displaystyle z\in\\{0,1\\}^{n}.$

Provided that we allow for linear terms (in both formulations), it is well
known that solutions to (1) can be used to compute a solution to (2) which
differs only by a constant factor, and vice-versa, due to the equivalence
$z=\frac{x+e}{2}$ if $z\in\\{0,1\\}^{n}$ and $x\in\\{-1,1\\}^{n}$, where
$e\in\mathbb{R}^{n}$ is the all ones vector of dimension $n$.

Although (1) and (2) cover many applications of interest, they are
intrinsically difficult to solve; computing optimal solutions to either (1) or
(2) is NP-Hard in general. Following the seminal work of Lovász [43] and the
theoretical and practical development of Interior Point Methods (IPMs) for
solving semidefinite optimization (SDO) problems [46, 47, 48, 49, 50, 55, 56],
a prevailing approach has been to obtain approximate solutions to (1) and (2)
by relaxing integrality and lifting the problem from a vector space of
dimension $n$, to the space of $n\times n$ symmetric matrices. The quadratic
form $x^{\top}Cx$ can be equivalently expressed by
$\operatorname{tr}{(Cxx^{\top})}$, where $\operatorname{tr}{(U)}$ denotes the
sum of the diagonal elements (or, trace) of a matrix $U\in\mathbb{R}^{n\times
n}$. To deal with the bilinear term $xx^{\top}$, we introduce a matrix
variable $X\in\mathbb{R}^{n\times n}$, and require that $X$ satisfies the
following:

$\textup{diag}(X)=e,\quad X\succeq 0,\quad\textup{rank}(X)=1,$

where the notation $U\succeq V$ means that the matrix $U-V$ is a symmetric
positive semidefinite matrix. Under these requirements, $X$ is guaranteed to
be of the form $X=xx^{\top}$ for $x\in\\{-1,1\\}^{n}$. The rank constraint,
however, is not convex, and thus dropping it yields the following (convex) SDO
relaxation of (1):

$\displaystyle\max$ $\displaystyle\operatorname{tr}{(CX)}$ (3) s.t.
$\displaystyle\operatorname{diag}{(X)}=e,\quad X\succeq 0.$

Although the optimal solution $X^{*}$ to (3) is no longer guaranteed to
satisfy $X^{*}=x^{*}x^{*^{\top}}$ and may not be integral in general, the
approximation of $x^{*}$ provided by $X^{*}$ is of sufficient quality to
justify its use. In fact, SDO approximations cover some of the most celebrated
results in optimization, such as the 0.878-approximation guarantee of Goemans
and Williamson for MaxCut [28] and the Lovász-$\vartheta$ number [43].

###  1.1 Literature Review

More generally, a (primal) SDO problem involving $n\times n$ matrices and $m$
constraints is of the form

$\displaystyle\max_{X}\quad$ $\displaystyle\operatorname{tr}{(CX)}$ s.t.
$\displaystyle\operatorname{tr}\left(A_{i}X\right)=b_{i}\quad\text{for}\leavevmode\nobreak\
i\in[m],$ $\displaystyle\quad\quad\quad X\succeq 0,$

where $[m]=\\{1,\dots,m\\}$ and $A_{1},\dots,A_{m},C\in{\cal S}^{n}$, and
$b\in\mathbb{R}^{m}$ are the (given) problem data. The dual SDO problem
associated with the primal is given by

$\displaystyle\min_{(u,S)}\quad$ $\displaystyle b^{\top}u$ s.t. $\displaystyle
S=\sum_{i=1}^{m}u_{i}A_{i}-C\succeq 0.$

where $S$ is the dual slack matrix.111While the dual variable is typically
denoted by $y$ rather than $u$, it is also customary in the literature to use
$y$ to denote a certain state preparation pair, and we do so later in this
paper. The classical literature on algorithms for solving SDO problems is rich
and can be categorized into two classes; algorithms that depend poly-
logarithimically on the inverse precision to which we solve the problem and
the size of the minimally inscribed ellipsoid, and algorithms that depend
polynomially on these quantities but exhibit an advantage with respect to $n$
and $m$. For instances with $m\leq\sqrt{n}$, the cutting plane methods (CPMs)
of [35, 42] are the best performing classical algorithms,222We remark that the
running time in [35] does however exhibit improved dependence with respect to
poly-logarithmic factors compared to the running time of [42]. and can solve
SDO problems in time

${\cal
O}\left(m(mns+m^{2}+n^{\omega})\cdot\textup{polylog}\left(m,n,R,\frac{1}{\epsilon}\right)\right),$

where $\omega\in[2,2.38]$ is the matrix multiplication exponent, $R$ is an
upper bound on the trace of a primal optimal solution $X$ (which can be
exponentially large), $\epsilon$ is the precision parameter, $s$ denotes the
maximum number of nonzeros per row of the input matrices and hence, ${\cal
O}(mns)$ is the total number of nonzeros in the constraints of SDO problem.
However, we typically have $m\in[\Omega(n),{\cal O}(n^{2})]$, in which case
the CPMs given in [35, 42] are outperformed by the IPM for SDO from Jiang et
al. [34]. Their IPM exhibits a worst case running time of

${\cal
O}\left(\sqrt{n}(mns+m^{\omega}+n^{\omega})\cdot\textup{polylog}\left(m,n,\frac{1}{\epsilon}\right)\right),$

where the term $m^{\omega}+n^{\omega}$ represents the per-iteration cost of
inverting the Hessian and matrices of the variables.

While quantum SDO solvers could also be categorized in a somewhat similar
fashion, it is perhaps more natural to do so according to how they attempt to
obtain quantum speedups. In this case we also have two classes; at a high
level, all proposed quantum SDO solution methodologies quantize a classical
algorithm by either using quantum linear system algorithms (QLSAs) [12, 14,
31], or a quantum mechanical interpretation of normalized positive
semidefinite matrices. We now review these works in detail.

The former class is comprised of algorithms that quantize IPMs, giving rise to
quantum IPMs (QIPMs). QIPMs attempt to speedup the bottleneck of the classical
IPM by substituting the classical solution of the Newton linear system with
the combined use of QLSA and quantum state tomography (with some classical
computation between iterates). Augustino et al. [6] present a convergent QIPM
for SDO, avoiding the shortcomings prevalent in early works on QIPMs (see,
e.g., [39]), by properly symmetrizing the Newton linear system, and utilizing
an orthogonal subspace representation of the search directions. This
representation guarantees that primal and dual feasibility are satisfied
exactly by all the iterates generated by inexact solutions of the Newton
linear system obtained via quantum subroutines. The worst case complexity of
their algorithm is

$\widetilde{{\cal
O}}_{n,\kappa,\frac{1}{\epsilon}}\left(\sqrt{n}\left(\frac{n^{3}\kappa^{2}}{\epsilon}+n^{4}\right)\right),$

where $\kappa$ is an upper bound on the condition numbers of the intermediate
Newton linear system coefficient matrices that arise over the course of the
algorithm. Here, the notation $\widetilde{{\cal O}}_{a,b}(f(x))$ suppresses
poly-logarithmic factors in $f(x)$, $a$ and $b$ that appear in the overall
running time, i.e., $\widetilde{{\cal O}}_{a,b}(f(x))\equiv{\cal
O}(f(x)\cdot\textup{polylog}(a,b,f(x)))$. While this QIPM achieves a speedup
in $n$ over the IPM from [34] when $m={\cal O}(n^{2})$, its dependence on
$\kappa$ and $\epsilon$ suggest no quantum advantage overall: the complexity
of the classical IPM does not depend on $\kappa$ and its dependence on
$\epsilon^{-1}$ is logarithmic. As the authors in [6] note, dependence on the
condition number bound $\kappa$ is particularly problematic in the context of
IPMs.

The second class of quantum SDO solvers are those that quantize algorithms
based on matrix exponentials and Gibbs states. The most prominent example is
the Matrix Multiplicative Weights Update (MMWU) Method of Arora and Kale [3],
which can solve SDO problems in time

$\widetilde{{\cal
O}}_{n,R,\frac{1}{\epsilon}}\left(nms\left(\frac{Rr}{\epsilon}\right)^{4}+ns\left(\frac{Rr}{\epsilon}\right)^{7}\right).$

where $r$ is a known $\ell_{1}$-norm upper bound333It is also assumed that
$R,r\geq 1$. on a dual optimal solution $u$. Unlike IPMs, the MMWU framework
does not involve the solution of linear systems; rather, these algorithms
alternate between candidate solutions to the primal and dual SDO problems.
IPMs and MMWUs also employ different definitions of optimality; for IPMs,
$\epsilon$-optimality implies that the primal and dual feasible solutions
exhibit a normalized duality gap bounded by $\epsilon$, i.e.:

$\frac{\operatorname{tr}\left(XS\right)}{n}\leq\epsilon,$

whereas an $\epsilon$-optimal solution obtained using an MMWU approximates the
optimal objective value to additive error $\epsilon$ (via binary search).
Finally, we point out a distinction between these algorithms with respect to
output. While primal-dual IPMs return the primal-dual optimal solution
$(X,u,S)$, MMWUs report $u$, but may avoid explicitly reporting $X$ and $S$ to
maintain the speedups they offer with respect to $n$. Reporting $X$ or $S$
under the MMWU framework necessitates the computation of matrix exponentials,
which may impose a considerable overhead because it generally resorts to
matrix multiplication.

The MMWU framework has been specialized to solve SDO problems of the form in
(3) (see, e.g., [4]), and the current state of the art is attributed to Lee
and Padmanabhan [41], who give an algorithm that can solve (3) to additive
error $\|C\|_{\ell_{1}}\epsilon$ with overall complexity

$\widetilde{{\cal O}}_{n,\frac{1}{\epsilon}}\left(ns\epsilon^{-3.5}\right),$

where $\|C\|_{\ell_{1}}=\sum_{i,j}|C_{ij}|$. It is important to note however,
that to achieve the stated complexity their methodology does not explicitly
report444Alternatively, they report a “gradient” $G\in{\cal S}^{n}$ such that
$X=W\exp(G)W$ for a diagonal matrix $W$. the solution $X$ and the authors
assume $\sum_{i,j}|C_{ij}|=n$. To achieve the same error scaling as the
algorithms we present in this work, the algorithm in [41] would have overall
cost $\widetilde{{\cal
O}}_{n,\frac{1}{\epsilon}}\left(\|C\|_{\ell_{1}}^{3.5}ns\epsilon^{-3.5}\right)$,
see Section 5.3.

Brandão and Svore [11] and van Apeldoorn et al. [61] were the first to
quantize the MMWU framework, utilizing a clever interpretation of the primal
variables: Gibbs states, which can be efficiently prepared on a quantum
computer, naturally correspond to trace-normalized positive definite matrices.
The running time of these MMWU-based algorithms was subsequently improved [30,
60], and the current state of the art running time of the quantum MMWU (QMMWU)
algorithm for SDO problems is:

$\widetilde{{\cal
O}}_{n,s,R,\frac{1}{\epsilon}}\left(\left(\sqrt{m}+\sqrt{n}\frac{Rr}{\epsilon}\right)s\left(\frac{Rr}{\epsilon}\right)^{4}\right).$

Similar to the complexity of QIPMs, QMMWU algorithms are faster with respect
to $m$ and $n$ when compared to their classical counterparts, but these
algorithms still exhibit a non-polynomial running time, due to their
polynomial dependence on the scale invariant parameter $\frac{Rr}{\epsilon}$,
whereas the natural input size depends on the logarithm of this quantity.

Seeking to improve the performance of quantum SDO solvers, Brandão et al. [10]
present an algorithm, which they call Hamiltonian Updates (HU), for solving
the SDO approximation (3) of (1). The HU method is a primal-only algorithm
closely related to the QMMWU framework, in that it leverages a Gibbs state
representation of the primal variable and progression towards the optimal
solution is made via matrix-exponentiated gradient updates. Specifically, the
authors in [10] are interested in solving an SDO feasibility problem that
arises upon renormalizing and relaxing (3):

find $\displaystyle X$ (4) s.t.
$\displaystyle\operatorname{tr}{\left(\frac{C}{\|C\|}X\right)}\geq\gamma-\epsilon$
$\displaystyle\sum_{i\in[n]}\left|\langle
i|X|i\rangle-\frac{1}{n}\right|\leq\epsilon$
$\displaystyle\operatorname{tr}{(X)}=1,\quad X\succeq 0.$

Here, $\gamma$ is an upper bound on the absolute value of the optimal
objective value of (3) when the cost matrix $C$ is normalized, obtained via
binary search over $[-1,1]$, and $\left|{i}\right\rangle$ for
$i\in\\{1,\dots,n\\}$ are the computational basis states. Since any
$\log(n)$-qubit Gibbs state is an element of the set
$\\{X\in\mathbb{R}^{n\times n}:\operatorname{tr}(X)=1,X\succeq 0\\}$ by
definition, solutions to (4) can be naturally be expressed as a Gibbs state

$\rho=\frac{\exp(-H)}{\operatorname{tr}(\exp(-H))},$

where $H$ is the Hamiltonian associated with $\rho$. The key observation in
[10] is that upon using the Gibbs state change of variables in (4), one can
model the $n$ constraints on the diagonal elements as single constraint which
requires that the distribution on the diagonal elements of a feasible solution
$\rho$ to (4) be at most $\epsilon$ in total variation distance to the uniform
distribution. In other words, the task of solving (4) reduces to finding a
$\log(n)$-qubit mixed quantum state that upon measurement in the computational
basis is approximately indistinguishable from the maximally-mixed state, and
whose trace inner product with the normalized cost matrix $C\|C\|^{-1}$ is at
least $\gamma-\epsilon$.

Using a quantum computer, the HU method of [10] solves (3) to additive error
${\cal O}\left(n\|C\|\epsilon\right)$ in time

$\widetilde{{\cal
O}}_{n,\frac{1}{\epsilon}}\left(n^{1.5}\sqrt{s}^{1+o(1)}\epsilon^{-28+o(1)}\exp\left(1.6\sqrt{\log(\epsilon^{-1})}\right)\right).$

The authors in [10] also provide an analysis of essentially the same algorithm
when using a classical computer, and show that the classical algorithm has a
complexity of

$\widetilde{{\cal
O}}_{n}\left(\min\\{n^{2}s,n^{\omega}\\}\epsilon^{-12}\right).$

The quantum algorithm yields a speedup in $n$ over classical algorithms, for a
specific class of SDO problems. However, as we have already seen with QIPMs
and QMMWU algorithms, its dependence on other parameters (in this case the
inverse precision) is prohibitive unless a very low precision solution is
acceptable. This raises the question as to whether the poor scaling in the
inverse precision can be mitigated without incurring additional cost in $n$
and $s$. We answer this question in the affirmative using iterative refinement
techniques.

Iterative Refinement (IR) is a methodology for computing high-precision
solutions to linear system of equations [29], as well as linear [25, 26, 27]
and mixed integer optimization problems [2, 17]. We summarize the methodology
at a high level as follows, and present a detailed discussion for the case of
convex feasibility problems later in the paper. Given an initial solution
$x^{(0)}\in\mathbb{R}^{d}$, at each iteration $k$ IR produces a refined
solution $x^{(k+1)}\leftarrow x^{(k)}+u^{(k)}$, where $u^{(k)}$ acts as a
correction of the error $r^{(k)}$ associated with $x^{(k)}$, and is determined
by solving a refining problem induced by the current solution. These
operations can all be carried out using the same level of accuracy, called the
fixed-precision approach. Alternatively, one may increase the accuracy with
which the residuals $r^{(k)}$ are computed as compared to $u^{(k)}$, and this
approach is called a mixed precision approach [29, 62]. In this paper, we
utilize the fixed precision approach.

###  1.2 Contributions

In this paper we develop an IR scheme for SDO approximations of QUBO problems
that uses the HU algorithm of [10] as a subroutine. We show that proceeding in
this way allows one to exponentially improve the dependence on the inverse
precision for both the quantum and classical algorithms.

With the proposed IR scheme, the classical algorithm solves the SDO problem
(3) up to absolute error ${\cal O}(\epsilon)$ with worst-case complexity

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

This is a significant speedup compared to general-purpose SDO solvers, such as
IPMs. This algorithm can be quantized following a similar strategy to [10].
When provided access to quantum random access memory (QRAM), the quantum
algorithm takes

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

accesses to the QRAM and additional quantum gates (this is the standard way of
describing complexity in the QRAM model of computation), plus ${\cal O}(ns)$
classical arithmetic operations — note that simply reading the cost matrix $C$
takes ${\cal O}(ns)$ time.

Summarizing, the combination of HU with IR described in this paper provides
exponential speedups over the methodology proposed in [10] with respect to the
precision parameter $\epsilon$. To the best of our knowledge, our classical
and quantum algorithms are the fastest known algorithms in their respective
model of computation for this class of problems, and our quantum algorithm
provides a genuine asymptotic speedup over known classical solution
methodologies, provided that we have access to QRAM. In the sparse-access
input model (without QRAM), the algorithm takes $\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}\left(n^{1.5}s^{0.5+o(1)}\right)$ accesses
to an oracle describing the coefficient matrix $C$ and $\widetilde{{\cal
O}}_{n,\|C\|_{F},\frac{1}{\epsilon}}\left(n^{2.5}s^{0.5+o(1)}\right)$
additional gates, therefore yielding no quantum speedup (the quantum gate
complexity is asymptotically larger than the classical complexity).

The remainder of this paper is organized in the following manner. Section 2
introduces notation, as well as the relevant input models and quantum
subroutines. In Section 3 we introduce the Hamiltonian Updates (HU) algorithm
from [10], and our Iterative Refinement scheme for SDO approximations of QUBOs
is presented in Section 4. The running time analysis is performed in Section
5, and Section 6 concludes the manuscript.

##  2 Preliminaries

We write $[n]$ to represent the set of elements $\\{1,\dots,n\\}$. We denote
the $i$-th element of a vector $x\in\mathbb{R}^{n}$ by $x_{i}$ for $i\in[n]$,
and the $ij$-th element of a matrix $A\in\mathbb{R}^{m\times n}$ by $A_{ij}$
for $i\in[m]$ and $j\in[n]$. To refer to the $i$-th row of a matrix $A$, we
write $A_{i,\cdot}$ and write $A_{\cdot,j}$ when referring to its $j$-th
column. We distinguish the quantity $a$ to the $k$-th power and the value of
$a$ at iterate $k$ using round brackets, writing $a^{k}$ and $a^{(k)}$ to
denote these quantities, respectively.

The smallest and largest singular values of a matrix $A$ are denoted
$\sigma_{\min}(A),\sigma_{\max}(A)$, and if $A\in{\cal S}^{n}$, then the
smallest and largest eigenvalues are denoted
$\lambda_{\min}(A),\lambda_{\max}(A)$. We let ${\cal S}_{+}^{n}$ and ${\cal
S}_{++}^{n}$ represent the cones of symmetric positive semidefinite, and
symmetric positive definite matrices, respectively. For $A,B\in{\cal S}^{n}$,
we write $A\succeq B$ ($A\succ B$) to indicate that the matrix $A-B$ is
symmetric positive semidefinite (symmetric positive definite), i.e.,
$A-B\in{\cal S}^{n}_{+}$ ($A-B\in{\cal S}^{n}_{++}$). The matrix exponential
$\exp(A)$, which is defined by the power series

$\exp(A)=I+A+\frac{1}{2!}A^{2}+\frac{1}{3!}A^{3}+\cdots,$

maps symmetric matrices to the space of symmetric positive definite matrices.
Given the spectral decomposition $A=V\Lambda V^{\top}$, then
$\exp(A)=V\exp(\Lambda)V^{\top}$, where
$\exp(\Lambda)=\operatorname{diag}(\exp(\Lambda_{11}),\exp(\Lambda_{22})\dots,\exp(\Lambda_{nn}))$.

We let $A\circ B$ denote the Hadamard (or element-wise) product of two
matrices, and $A\otimes B$ denotes their tensor product. Later in this work,
we make use of the following facts regarding Hadamard products.

######  Lemma 1 (Lemma 5.1.4 in [33]).

Let $E$, $F$ and $G$ be $m\times n$ matrices. Then, the $i$-th diagonal entry
of the matrix $(E\circ F)G^{\top}$ coincides with the $i$-th diagonal entry of
the matrix $(E\circ G)F^{\top}$. That is,

$[(E\circ F)G^{\top})]_{ii}=[(E\circ G)F^{\top})]_{ii}\quad\forall i\in[m].$

######  Lemma 2 (Theorem 5.3.4 in [33]).

Let $A$ and $B$ be $n\times n$ Hermitian matrices. If $A\in{\cal S}^{n}_{+}$,
then any eigenvalue $\lambda(A\circ B)$ of $A\circ B$ satisfies

$\lambda_{\min}(A)\cdot\lambda_{\min}(B)\leq\lambda(A\circ
B)\leq\lambda_{\max}(A)\cdot\lambda_{\max}(B).$

######  Corollary 1.

Suppose $A$ and $B$ are $n\times n$ Hermitian matrices. If $A\in{\cal
S}^{n}_{+}$ and $B\in{\cal S}^{n}$ with $B$ having at least one negative
eigenvalue, then

$\lambda_{\max}(A)\cdot\lambda_{\min}(B)\leq\lambda_{\min}(A\circ B).$

###### Proof.

First, note that by Lemma 2, we have $\lambda_{\min}(A\circ
B)\geq\lambda_{\min}(A)\cdot\lambda_{\min}(B)$. Moreover, since $B$ has at
least one negative eigenvalue, we have $\lambda_{\min}(B)<0$, which combined
with the fact that $\lambda_{\min}(A)\geq 0$ yields
$\lambda_{\max}(A)\cdot\lambda_{\min}(B)\leq\lambda_{\min}(A)\cdot\lambda_{\min}(B)$.
Hence,

$\lambda_{\max}(A)\cdot\lambda_{\min}(B)\leq\lambda_{\min}(A\circ B),$

and the proof is complete. ∎

We write $e$ to refer to the vector of all ones in $\mathbb{R}^{n}$, and use
the notation $e_{i}$ to refer to the $i$-th unit vector in the standard
orthonormal basis $\\{e_{1},\dots,e_{n}\\}$ for $\mathbb{R}^{n}$. Analogously,
the computational basis states are denoted by $\left|{i}\right\rangle$ for
$i\in[n]$. Hence, for $x\in\mathbb{R}^{n}$, we denote its amplitude encoding
by $\left|{x}\right\rangle$, defined as

$\left|{x}\right\rangle=\frac{1}{\|x\|}\sum_{i\in[n]}x_{i}\left|{i}\right\rangle.$

Observe that $\left|{x}\right\rangle$ is a $\log(n)$-qubit state; for
simplicity, we assume that the dimensions of all spaces are powers of 2. All
logarithms are base 2.

Where appropriate, our analysis makes use of the Schatten $p$-norm, defined
for a bounded linear operator $A$ as

$\|A\|_{p}=\left[\operatorname{tr}\left(|A|^{p}\right)\right]^{\frac{1}{p}},$

where $|A|=(A^{\dagger}A)^{\frac{1}{2}}$ with $A^{\dagger}$ denoting the
conjugate transpose of $A$. Notice that the trace and operator norms
$\|\cdot\|_{\operatorname{tr}}$ and $\|\cdot\|$ are the Schatten-1 and
Schatten-$\infty$ norms, respectively, and the Frobenius norm $\|\cdot\|_{F}$
corresponds to the Schatten-2 norm.

For a scalar $x\in\mathbb{R}$ define the sign function $\textup{sign}(x)$ as

$\textup{sign}(x)=\begin{cases}-1&\text{if}\leavevmode\nobreak\ x<0\\\
0&\text{if}\leavevmode\nobreak\ x=0\\\ 1&\text{if}\leavevmode\nobreak\
x>0.\end{cases}$

When $x\in\mathbb{R}^{n}$,
$\textup{sign}(x)=(\textup{sign}(x_{1}),\dots,\textup{sign}(x_{n}))^{\top}$.

For any positive integer $q$, and binary strings $j,k\in\\{0,1\\}^{q}$, we
denote by $j\oplus k$ the bitwise modulo 2 addition of $q$-digit strings,
defined as

$j\oplus k=h$

where $h\in\\{0,1\\}^{q}$ is the bitstring whose elements $h_{p}$ are defined
for $p\in[q]$ as

$h_{p}=\begin{cases}0&\text{if}\leavevmode\nobreak\ j_{p}=k_{p},\\\
1&\text{otherwise}.\end{cases}$

### “Big-O” notation

We define ${\cal O}(\cdot)$ as

$f(x)={\cal
O}(g(x))\iff\exists\ell\in\mathbb{R},c\in\mathbb{R}_{+},\leavevmode\nobreak\
\text{such that}\leavevmode\nobreak\ f(x)\leq cg(x)\quad\forall x>\ell.$

We write $f(x)=\Omega(g(x))\iff g(x)={\cal O}(f(x))$. We also define
$\widetilde{{\cal O}}(f(x))={\cal O}(f(x)\cdot\textup{polylog}(f(x)))$ and
when the function depends poly-logarithmically on other variables we write

$\widetilde{{\cal O}}_{a,b}\leavevmode\nobreak\ (f(x))={\cal
O}(f(x)\cdot\textup{polylog}(a,b,f(x))).$

###  2.1 Input models and subroutines

For our quantum algorithm, we provide analyses for two distinct models of
input. One model considers a quantum-read/classical-write RAM (QRAM), and the
other is the sparse-access model, which we use to bound the running time
without access to QRAM.

####  2.1.1 Sparse-access model

In the sparse-access model, the input matrix $C$ is assumed to be $s$-row
sparse for some known bound $s\in[n]$. In other words, $C$ has at most $s$
nonzero entries per row. The sparse-access model is closely related to the
classical notion, in that we assume access to an oracle $O_{\text{sparse}}$,
which upon being queried with input $(i,j)$ returns the index of the $j$-th
nonzero entry of the $i$-th row of $C$ by calculating the index function:

$\operatorname{index}:[n]\times[s]\to[n].$

That is, for $i\in[n]$ and $j\in[s]$, $O_{\text{sparse}}$ computes the
position in place:

$O_{\text{sparse}}\left|{i,j}\right\rangle=\left|{i,\operatorname{index}(i,j)}\right\rangle.$

We also assume access to an oracle that returns a bitstring representation of
the individual entries of the normalized cost matrix $C\|C\|_{F}^{-1}$ for
every $i,j\in[n]$:

$O_{C}\left|{i,j,z}\right\rangle=\left|{i,j,z\oplus(C_{ij}\|C\|_{F}^{-1})}\right\rangle.$

####  2.1.2 Quantum random access memory

We consider a quantum-read/classical-write RAM (QRAM), which enables us to
store classical data that our quantum algorithms can make oracle calls to.
This type of storage is the direct quantum analog of classical RAM: it enables
a quantum algorithm to access classical data in superposition. Accessing a
QRAM of size $n$ takes ${\cal O}(n)$ gates [5, 24], but these gates can be
arranged in parallel so that the circuit depth remains ${\cal
O}(\textup{polylog}(n))$. Therefore we make the assumption (standard in the
literature on quantum algorithms) that the cost of accessing a QRAM of size
$n$ is ${\cal O}(\textup{polylog}(n))$.

The next result from Chakraborty et al. [12], is adapted from an earlier
result of Kerenidis and Prakash [38] and summarizes the aspects of the data
structure we utilize.

######  Theorem 1 (Theorem 1 in [12]).

(Implementing quantum operators using an efficient data structure) Let
$A\in\mathbb{R}^{m\times n}$ be a matrix. If $w$ is the number of non-zero
entries of $A$, then there exists a data structure of size ${\cal
O}\left(w\log^{2}(mn)\right)$ that, given the entries $(i,j,A_{ij})$ in an
arbitrary order, stores them such that time taken to store each entry of $A$
is ${\cal O}(\log(mn))$. Once this data structure has been initiated with all
non-zero entries of $A$, there exists a quantum algorithm that can perform the
following maps with $\xi$-precision in time ${\cal
O}\left(\textup{polylog}\left(\frac{mn}{\xi}\right)\right)$:

$\displaystyle\widetilde{U}:$
$\displaystyle\left|{i}\right\rangle\left|{0}\right\rangle\mapsto\left|{i}\right\rangle\frac{1}{\|A_{i,\cdot}\|}\sum_{j=1}^{n}A_{ij}\left|{j}\right\rangle=\left|{i,A_{i,\cdot}}\right\rangle,$
$\displaystyle\widetilde{V}:$
$\displaystyle\left|{0}\right\rangle\left|{j}\right\rangle\mapsto\frac{1}{\|A\|_{F}}\sum_{i=1}^{m}\|A_{i,\cdot}\|\left|{i}\right\rangle\left|{j}\right\rangle=\left|{\widetilde{A},j}\right\rangle,$

where $\left|{A_{i,\cdot}}\right\rangle$ is the normalized quantum state
corresponding to the $i$-th row of $A$ and
$\left|{\widetilde{A}}\right\rangle$ is a normalized quantum state such that
$\langle i|\widetilde{A}\rangle=\|A_{i,\cdot}\|$, i.e., the norm of the $i$-th
row of $A$.

####  2.1.3 Working with block-encoded matrices

We now give a formal definition of a block-encoding from [12].

######  Definition 1 (Block-encoding).

Let $A\in\mathbb{C}^{2^{w}\times 2^{w}}$ be a $w$-qubit operator. Then, a
$(w+a)$-qubit unitary $U$ is an $(\alpha,a,\xi)$-block-encoding of $A$ if
$U=\begin{pmatrix}\widetilde{A}&\cdot\\\ \cdot&\cdot\end{pmatrix}$, with the
property that

$\|\alpha\widetilde{A}-A\|\leq\xi.$

It was shown by Kerenidis and Prakash [38] and Chakraborty et al. [12] how to
efficiently implement block-encodings of matrices that are stored in a QRAM
data structure, which is formalized in the next result.

######  Lemma 3 (Lemma 3.3.7 in [22]).

Let $A\in\mathbb{C}^{2^{w}\times 2^{w}}$ and $\xi>0$.

  1. (i)

Fix $q\in[0,2]$ and define $\mu_{q}(A)=\sqrt{n_{q}(A)n_{(2-q)}(A^{\top})}$
where $n_{q}(A)=\max_{i}\|A_{i,\cdot}\|_{q}^{q}$ is the $q$-th power of the
maximum $q$-norm of the rows of $A$. Defining $A^{\\{q\\}}$ to be the matrix
with elements $A^{\\{q\\}}_{ij}=\sqrt{A_{ij}^{q}}$, if $A^{\\{q\\}}$ and
$(A^{\\{2-q\\}})^{{\dagger}}$ are both stored in QRAM data structures, then
there exist unitaries $U_{R}$ and $U_{L}$ that can be implemented in time
${\cal O}(\textup{poly}(w\log\frac{1}{\xi}))$ and such that
$U^{{\dagger}}_{R}U_{L}$ is a $(\mu_{q}(A),w+2,\xi)$-block-encoding of $A$.

  2. (ii)

If $A$ is stored in a QRAM data structure, then there exist unitaries $U_{R}$
and $U_{L}$ that can be implemented in time ${\cal
O}(\textup{poly}(w\log\frac{1}{\xi}))$ and such that $U^{{\dagger}}_{R}U_{L}$
is an $(\|A\|_{F},w+2,\xi)$-block-encoding of $A$.

Linear combinations of block-encodings can also be constructed at cost that is
merely logarithmic in the dimension.

######  Definition 2 (Definition 3.3.8 in [22]).

(State preparation pair) Let $y\in\mathbb{C}^{m}$ and $\|y\|_{1}\leq\beta$.
The pair of unitaries $(P_{L},P_{R})$ is called a $(\beta,p,\xi)$-state-
preparation-pair if $P_{L}\left|{0}\right\rangle^{\otimes
p}=\sum_{j=0}^{{2^{p}}-1}c_{j}\left|{j}\right\rangle$ and
$P_{R}\left|{0}\right\rangle^{\otimes
p}=\sum_{j=1}^{{2^{p}}-1}d_{j}\left|{j}\right\rangle$ such that
$\sum_{j=0}^{m-1}|\beta(c_{j}^{*}d_{j})-y_{j}|\leq\xi$ and for all $j\in
m,\dots,2^{p}-1$ we have $c_{j}^{*}d_{j}=0$.

######  Proposition 1 (Lemma 52 in [23]).

(Linear combination of block-encoded matrices, with weights given by a state
preparation pair) Let $A=\sum_{j=0}^{m-1}y_{j}A_{j}$ be a $w$-qubit operator,
where $A_{j}$ are matrices. Suppose $P_{L},P_{R}$ is a
$(\beta,p,\xi_{1})$-state-preparation pair for $y$,
$W=\sum_{j=0}^{m-1}\left|{j}\right\rangle\left\langle{j}\right|\otimes
U_{j}+((I-\sum_{j=0}^{m-1}\left|{j}\right\rangle\left\langle{j}\right|)\otimes
I_{a}\otimes I_{s})$ is an $(w+a+p)$-qubit unitary with the property that
$U_{j}$ is an $(\alpha,a,\xi_{2})$-block-encoding of $A_{j}$. Then we can
implement a $(\alpha\beta,a+p,\alpha\xi_{1}+\alpha\beta\xi_{2})$-block-
encoding of $A$ with a single use of $W,P_{R}$ and $P_{L}^{{\dagger}}$.

It turns out that the sparse-access model reduces to the quantum operator
model upon choosing $\alpha=s$ (if row and column sparsity are the same). The
next result from [23] describes how to implement block-encodings using the
sparse-access input model, and the associated costs.

######  Lemma 4 (Lemma 48 in [23]).

Let $A\in\mathbb{C}^{2^{w}\times 2^{w}}$ be a matrix that is $s_{r}$-row-
sparse and $s_{c}$-column-sparse, and each element of $A$ has absolute value
at most 1. Suppose that we have access to the following sparse-access oracles
acting on two $(w+1)$ qubit registers:

$\displaystyle O_{r}:\left|{i}\right\rangle\left|{k}\right\rangle$
$\displaystyle\mapsto\left|{i}\right\rangle\left|{r_{ik}}\right\rangle\quad\forall
i\in[2^{w}]-1,k\in[s_{r}],\leavevmode\nobreak\ \text{and}$ $\displaystyle
O_{c}:\left|{\ell}\right\rangle\left|{j}\right\rangle$
$\displaystyle\mapsto\left|{c_{\ell
j}}\right\rangle\left|{j}\right\rangle\quad\forall\ell\in[s_{c}],j\in[2^{w}]-1,\leavevmode\nobreak\
\text{where}$

$r_{ij}$ is the index for the $j$-th non-zero entry of the $i$-th row of $A$,
or if there are less than $i$ non-zero entries, then it is $j+2^{w}$, and
similarly $c_{ij}$ is the index for the $i$-th non-zero entry of the $j$-th
column of $A$, or if there are less than $j$ non-zero entries, then it is
$i+2^{w}$. Additionally, assume that we have access to an oracle $O_{A}$ that
returns the entries of $A$ in a binary description:

$O_{A}:\left|{i}\right\rangle\left|{j}\right\rangle\left|{0}\right\rangle^{\otimes
p}\mapsto\left|{i}\right\rangle\left|{j}\right\rangle\left|{a_{ij}}\right\rangle,\quad\forall
i,j\in[2^{w}]-1,$

where $a_{ij}$ is a $p$-bit binary description of the $ij$-matrix element of
$A$. Then, we can implement a $(\sqrt{s_{r}s_{c}},w+3,\xi)$-block-encoding of
$A$ with a single use of $O_{r}$, $O_{c}$ and two uses of $O_{A}$, and
additionally using ${\cal
O}\left(w+\log^{2.5}\left(\frac{s_{r}s_{c}}{\xi}\right)\right)$ one and two
qubit gates while using ${\cal
O}\left(p+\log^{2.5}\left(\frac{s_{r}s_{c}}{\xi}\right)\right)$ ancilla
qubits.

The block-encoding framework will be useful in speeding up the overall running
time found in [10], as it allows us to perform matrix computations and
Hamiltonian simulation efficiently.

######  Theorem 2 (Corollary 3.4.7 in [22]).

(Optimal block-Hamiltonian simulation) Suppose that $U$ is an
$(\alpha,a,\xi/|2t|)$-block-encoding of the Hamiltonian $H$. Then, we can
implement a $\xi$-precise Hamiltonian simulation unitary $V$ which is an
$(1,a+2,\xi)$-block-encoding of $e^{itH}$, with ${\cal O}\left(|\alpha
t|+\frac{\log(1/\xi)}{\log\log(1/\xi)}\right)$ uses of controlled-$U$ or its
inverse and with ${\cal O}\left(a|\alpha
t|+a\frac{\log(1/\xi)}{\log\log(1/\xi)}\right)$ two-qubit gates.

Additionally, one can easily take the product of block-encodings.

######  Proposition 2 (Lemma 4 in [12]).

(Product of block-encoded matrices) If $U_{A}$ is an
$(\alpha_{1},a_{1},\xi_{A})$-block-encoding of an $s$-qubit operator $A$, and
$U_{B}$ is an $(\alpha_{2},a_{2},\xi_{B})$-block-encoding of an $s$-qubit
operator $B$, then $(I_{a_{2}}\otimes U_{A})(I_{a_{1}}\otimes U_{B})$ is an
$(\alpha_{1}\alpha_{2},a_{1}+a_{2},\alpha_{1}\xi_{B}+\alpha_{2}\xi_{A})$-block-
encoding of $AB$.

Relevant to our work in the quantum operator input model is the idea of block-
encoding the Hadamard, or element-wise product of two matrices. We will
demonstrate how one can carry out the Hadamard product of block-encodings of
matrices $A$ and $B$ as a reduction of the Kronecker product of block-
encodings, which is straightforward to construct given block encodings of $A$
and $B$.

######  Proposition 3.

(Kronecker product of block-encoded matrices) Suppose that $U_{A}$ is an
$(\alpha_{1},a_{1},\xi_{A})$-block-encoding of $A\in\mathbb{R}^{n\times n}$,
and $U_{B}$ is an $(\alpha_{2},a_{2},\xi_{B})$-block-encoding of
$B\in\mathbb{R}^{n\times n}$. Then, taking the tensor product of $U_{A}$ and
$U_{B}$, we obtain a
$(\alpha_{1}\alpha_{2},a_{1}+a_{2},\xi_{A}+\xi_{B})$-block-encoding of
$A\otimes B$.

We do not give a formal proof here as the result directly follows from the
definition of a block-encoding; to obtain the tensor product of two block-
encoded matrices, it suffices to take the tensor product of their block-
encodings while keeping the ancilla qubits separate.

######  Proposition 4 (Hadamard product of block-encoded matrices).

Suppose that $U_{A}$ is an $(\alpha_{1},a_{1},\xi_{A})$-block-encoding of
$A\in\mathbb{R}^{n\times n}$, and $U_{B}$ is a
$(\alpha_{2},a_{2},\xi_{B})$-block-encoding $B\in\mathbb{R}^{n\times n}$.
Then, using $U_{A}$ and $U_{B}$, we can implement an
$(\alpha_{1}\alpha_{2},a_{1}+a_{2}+8\log(n)+12,5(\xi_{A}+\xi_{B}))$-block-
encoding of $A\circ B$ using one application of $U_{A}$ and $U_{B}$, and
$\widetilde{{\cal O}}_{n}(1)$ additional gates.

###### Proof.

First, note that

$A\circ B=(A\otimes B)[\iota_{A},\iota_{B}],$

where $\iota_{A}=\iota_{B}=\\{1,n+2,2n+3,\dots,n^{2}\\}$ are index sets of
cardinality $n$ (see, e.g., Lemma 5.1.1 in [33]). Our goal is to use the index
sets $\iota_{A}$ and $\iota_{B}$ along with a block encoding of $A\otimes B$
to construct a unitary which block-encodes ${\cal M}\in\mathbb{R}^{n^{2}\times
n^{2}}$, a matrix which contains the elements of $A\circ B$ in its upper left-
most $n\times n$ block, while all other entries are $0$:

${\cal M}_{ij}=\begin{cases}A_{ij}\cdot B_{ij}&\text{for}\leavevmode\nobreak\
i,j=1,\dots,n,\\\ 0&\text{otherwise},\end{cases}$

i.e.,

${\cal M}=\begin{pmatrix}A\circ B&\mathbf{0}^{n\times(n^{2}-n)}\\\
\mathbf{0}^{(n^{2}-n)\times
n}&\mathbf{0}^{(n^{2}-n)\times(n^{2}-n)}\end{pmatrix}.$

We will first show how one can use $\iota_{A}$ and $\iota_{B}$ to construct
sparse matrices that map $A\otimes B$ to ${\cal M}$, and then subsequently
analyze the cost of constructing the corresponding unitary block-encoding.

Consider the matrix $Z\in\mathbb{R}^{n^{2}\times n^{2}}$, whose elements are
defined as

$Z_{ij}=\begin{cases}1&\text{if}\leavevmode\nobreak\ i=j=(k-1)n+k,\quad
k=1,\dots,n,\\\ 0&\text{otherwise}.\end{cases}$

Multiplying $A\otimes B$ on the left by $Z$ sets the rows of $A\otimes B$
which do not contain elements of $A\circ B$ to zero, and subsequently
multiplying $Z(A\otimes B)$ on the right by $Z$ will set the columns of
$Z(A\otimes B)$ which do not appear in $A\circ B$ to zero. As a result, a
block-encoding of $Z(A\otimes B)Z$ corresponds to block-encoding $A\otimes B$,
and setting all terms not appearing in $A\circ B$ to zero:

$[Z(A\otimes B)Z]_{ij}=\begin{cases}[A\otimes
B]_{ij}&\text{if}\leavevmode\nobreak\ i=(k-1)n+k\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ j=(\ell-1)n+\ell\quad k,\ell=1,\dots,n,\\\
0&\text{otherwise}.\end{cases}$

Next, let $G\in\mathbb{R}^{n^{2}\times n^{2}}$ be a matrix whose elements are
defined as follows:

$G_{ij}=\begin{cases}1&\text{if}\leavevmode\nobreak\
i\in[n^{2}]\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\
i=j=(k-1)n+k,\quad k=1,\dots,n,\\\ 1&\text{if}\leavevmode\nobreak\
i\in[n^{2}]\setminus\\{1,n+2,2n+3,\dots,n^{2}\\}\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ j=(i-1)n+i,\\\ 0&\text{otherwise}.\end{cases}$

We will now establish that $GZ(A\otimes B)Z)G^{\top}$ is precisely the matrix
we seek to block-encode, by demonstrating that $G(Z(A\otimes
B)Z)G^{\top}={\cal M}$. First, observe that $G$ is a (partial) permutation
matrix: multiplying $Z(A\otimes B)Z$ on the left by $G$ performs the necessary
row-exchanges, as the elements of $G(Z(A\otimes B)Z)$ are given by

$[G\left(Z(A\otimes B)Z\right)]_{ik}=\begin{cases}A_{ij}\cdot
B_{ij}&\text{for}\leavevmode\nobreak\ k=(j-1)n+j,\quad i,j=1,\dots,n,\\\
0&\text{otherwise}.\end{cases}$

On the other hand, multiplying $Z(A\otimes B)Z)$ on the right by $G^{\top}$
performs this transformation with respect to the columns such that

$[\left(Z(A\otimes B)Z\right)G]_{kj}=\begin{cases}A_{ij}\cdot
B_{ij}&\text{for}\leavevmode\nobreak\ k=(i-1)n+i,\quad i,j=1,\dots,n,\\\
0&\text{otherwise}.\end{cases}$

Hence, multiplying $G\left(Z(A\otimes B)Z\right)$ on the right by $G^{\top}$
conducts the column exchanges to move $A\circ B$ to the top left
$n$-dimensional block of $Z(A\otimes B)Z$, i.e.,

$[G\left(Z(A\otimes B)Z\right)G]_{ij}=\begin{cases}A_{ij}\cdot
B_{ij}&\text{for}\leavevmode\nobreak\ i,j=1,\dots,n,\\\
0&\text{otherwise}.\end{cases}$

Therefore, $G(Z(A\otimes B)Z)G^{\top}={\cal M}$ as desired.

We now analyze the cost associated with block-encoding ${\cal M}$. Under the
stated hypothesis, we have access to an $(\alpha_{1},a_{1},\xi_{A})$-block-
encoding $U_{A}$ of $A$, and an $(\alpha_{2},a_{2},\xi_{B})$-block-encoding
$U_{B}$ of $B$, and thus applying Proposition 3 we can construct an
$(\alpha_{1}\alpha_{2},a_{1}+a_{2},\xi_{A}+\xi_{B})$-block-encoding
$U_{A\otimes B}$ of $A\otimes B$ using one application of $U_{A}$ and of
$U_{B}$, and no additional gates.

Using the description of $Z$, we can construct the sparse-access oracles
$O_{r}$ and $O_{c}$ as defined in Lemma 4 (which act on two $(2\log n+1)$
qubit registers). Additionally, from the definition of $Z$, we can construct
an oracle $O_{Z}$, which returns the entries of $Z$ in a binary description:

$O_{Z}:\left|{i}\right\rangle\left|{j}\right\rangle\left|{0}\right\rangle^{\otimes
p}\mapsto\left|{i}\right\rangle\left|{j}\right\rangle\left|{z_{ij}}\right\rangle,\quad\forall
i,j\in[2^{2\log n}]-1,$

where $z_{ij}$ is a $p$-bit binary description of the $ij$-matrix element of
$Z$. Note that the circuit for the position and value of the nonzero elements
of $Z$ using $\widetilde{{\cal O}}_{n}(1)$ gates because they admit an
efficient description: their value is 1 and we have a compact description of
their position. By construction the matrix $Z$ is 1-row sparse and 1-column
sparse, and hence an application of Lemma 4 with $s_{r}=s_{c}=1$ asserts that
one can construct a $(1,2\log(n)+3,\xi_{Z})$-block-encoding $U_{Z}$ of $Z$.
Given block-encodings $U_{Z}$ and $U_{A\otimes B}$, we can apply Proposition 2
with

$\xi_{Z}=\frac{\xi_{A}+\xi_{B}}{\alpha_{1}\alpha_{2}},\leavevmode\nobreak\
\quad\xi_{A\otimes B}=\xi_{A}+\xi_{B},$

yielding an
$(\alpha_{1}\alpha_{2},a_{1}+a_{2}+2\log(n)+3,2(\xi_{A}+\xi_{B}))$-block-
encoding of $Z(A\otimes B)$. Applying Proposition 2 once more with

$\xi_{Z}=\frac{\xi_{A}+\xi_{B}}{\alpha_{1}\alpha_{2}},\leavevmode\nobreak\
\quad\xi_{Z(A\otimes B)}=2(\xi_{A}+\xi_{B}),$

we obtain an
$(\alpha_{1}\alpha_{2},a_{1}+a_{2}+4\log(n)+6,3(\xi_{A}+\xi_{B}))$-block-
encoding of $Z(A\otimes B)Z$.

Just as was the case with $Z$, we can use the description of $G$ to construct
the sparse-access oracles $O_{r}$ and $O_{c}$ as defined in Lemma 4 (which
again, act on two $(2\log n+1)$ qubit registers), as well as an oracle $O_{G}$
using $\widetilde{{\cal O}}_{n}(1)$ gates, that returns the entries of $G$ in
a binary description:

$O_{G}:\left|{i}\right\rangle\left|{j}\right\rangle\left|{0}\right\rangle^{\otimes
p}\mapsto\left|{i}\right\rangle\left|{j}\right\rangle\left|{g_{ij}}\right\rangle,\quad\forall
i,j\in[2^{2\log n}]-1,$

where $g_{ij}$ is a $p$-bit binary description of $G_{ij}$ (the $ij$-matrix
element of $G$). Noting that $G$ is 1-row sparse and 1-column sparse (and
hence, so its transpose); applying Lemma 4 twice more allows us to construct a
$(1,2\log(n)+3,\xi_{G})$-block-encoding $U_{G}$ of $G$, as well as a
$(1,2\log(n)+3,\xi_{G}^{\top})$-block-encoding $U_{G^{\top}}$ of the transpose
$G^{\top}$. We can then use $U_{G}$ and our
$(\alpha_{1}\alpha_{2},a_{1}+a_{2}+4\log(n)+6,3(\xi_{A}+\xi_{B}))$-block-
encoding $U_{Z(A\otimes B)Z}$ of $Z(A\otimes B)Z$ to construct an
$(\alpha_{1}\alpha_{2},a_{1}+a_{2}+6\log(n)+9,4(\xi_{A}+\xi_{B}))$-block-
encoding of $G(Z(A\otimes B)Z)$ by applying Proposition 2 with

$\xi_{G}=\frac{\xi_{A}+\xi_{B}}{\alpha_{1}\alpha_{2}},\leavevmode\nobreak\
\quad\xi_{Z(A\otimes B)Z}=3(\xi_{A}+\xi_{B}).$

Applying Proposition 2 a final time, with

$\xi_{G^{\top}}=\frac{\xi_{A}+\xi_{B}}{\alpha_{1}\alpha_{2}},\leavevmode\nobreak\
\quad\xi_{G(Z(A\otimes B)Z)}=4(\xi_{A}+\xi_{B}),$

produces an
$(\alpha_{1}\alpha_{2},a_{1}+a_{2}+8\log(n)+12,5(\xi_{A}+\xi_{B}))$-block-
encoding $U_{{\cal M}}$ of ${\cal M}=G(Z(A\otimes B)Z)G^{\top}$.

The stated complexity result follows upon noting that the steps required to
construct the unitary

$U_{{\cal M}}=U_{G}U_{Z}U_{A\otimes B}U_{Z}U_{G^{\top}}$

requires one application of $U_{A\otimes B}$ and one application of each of
the other matrices. In turn, this amounts to 1 application of $U_{A}$ and
$U_{B}$ each, plus the $\widetilde{{\cal O}}_{n}(1)$ gate cost of the
remaining matrices $U_{G}$, $U_{Z}$ and $U_{G^{\top}}$, and the proof is
complete. ∎

We remark that a similar result to Proposition 4 was independently derived and
discussed in the recent paper [13].

####  2.1.4 Gibbs Samplers and Trace Estimators

For clarity, we begin with a formal definition of a subnormalized density
operators and their purifications.

######  Definition 3 (Definition 6.3.1 in [22]).

(Subnormalized density operators & Purification) A subnormalized density
operator $\rho$ is a positive semidefinite matrix of trace at most 1. A
purification $\varrho$ of a subnormalized density operator $\rho$ is a
3-register pure state such that tracing out the third register and projecting
on the subspace where the second register is $\left|{0}\right\rangle$ yields
$\rho$.

The frameworks introduced later in this paper require that we implement a
Gibbs sampler and a trace estimator, which we define next.

######  Definition 4 (Definition 4.11 in [58]).

(Gibbs Sampler) A $\theta$-precise Gibbs-sampler for the input matrix $H$, is
a unitary that takes as input a data structure storing a Hamiltonian $H$ and
creates as output a purification of a $\theta$-approximation (in trace
distance) of the Gibbs state

$\rho=\frac{\exp(-H)}{\operatorname{tr}(\exp(-H))}.$

We will use these approximate Gibbs states in order to check the diagonal
entries of our solutions, as well as compute the trace inner products of
matrices (or, expectation values), i.e., quantities of the form
$\operatorname{tr}(A\rho)$.

######  Definition 5 (Definition 4.12 in [58]).

(Trace Estimator) A $\theta$-precise trace estimator is a unitary that as
input takes a state $\rho$ and a matrix $A$. It outputs a sample from a random
variable $x\in\mathbb{R}$ such that $x$ is an estimator for
$\operatorname{tr}(A\rho)$ that is at most $\theta/4$ biased.

These implementations require polynomial approximations of the exponential
function, which can be obtained using quantum singular value transformation
techniques introduced in [22, 23].

######  Lemma 5 (Lemma 4.14 in [58]).

Let $\xi\in(0,1/6]$ and $\beta\geq 1$. There exists a polynomial $P(x)$ such
that

  * •

For all $x\in[-1,0]$, we have $\left|P(x)-\exp(2\beta x)/4\right|\leq\xi.$

  * •

For all $x\in[-1,1]$, we have $\left|P(x)\right|\leq 1/2.$

  * •

$\operatorname{deg}(P)=\widetilde{{\cal O}}_{\frac{1}{\xi}}(\beta)$.

######  Lemma 6 (Lemma 4.15 in [58]).

Let $\theta\in(0,1/3]$, $\beta>1$, and let $d$ be the degree of the polynomial
from Lemma 5 when we let $\xi=\frac{\theta}{128n}$. Let $U$ be a
$(\beta,a,\frac{\theta^{2}\beta}{1024^{2}d^{2}n^{2}})$-block-encoding of a
Hermitian operator $H\in\mathbb{R}^{n\times n}$, i.e,, a
$(\beta,a,\widetilde{{\cal O}}(\theta/\beta n^{2}))$-block-encoding. Then, we
can create a purification of a state $\tilde{\rho}$ such that

$\left\|\tilde{\rho}-\frac{\exp(H)}{\operatorname{tr}{(\exp(H))}}\right\|_{\operatorname{tr}}\leq\theta$

using $\widetilde{{\cal O}}_{\frac{1}{\theta}}(\sqrt{n}\beta)$ applications of
$U$ and $\widetilde{{\cal O}}_{\frac{1}{\theta}}(\sqrt{n}\beta a)$ elementary
operations.

Provided access to a unitary that prepares a purification of a density
operator, we can also construct a block-encoding of it. This is formalized in
the following lemma from [22], which was based on ideas found in [45,
Corollary 9].

######  Lemma 7 (Lemma 6.4.4 in [22]).

(Block-encoding of a (subnormalized) density operator) Let $G$ be a $(w+a)$
unitary which on the input state
$\left|{0}\right\rangle^{w}\left|{0}\right\rangle^{a}$ prepares a purification
$\left|{\varrho}\right\rangle$ of the subnormalized $w$-qubit density operator
$\rho$. Then we can implement a $(1,w+a,0)$-block-encoding of $\rho$ with a
single use of $G$ and its inverse and with $w+1$ two-qubit gates.

We are now in a position to define a trace estimator using the quantum
operator input model.

######  Lemma 8 (Lemma 4.18 in [58]).

Let $\rho$ be an $n$-dimensional quantum state and $U$ an
$(\alpha,a,\theta/2)$-block-encoding of a matrix $A\in\mathbb{R}^{n\times n}$
with $\|A\|\leq 1$. A trace estimator for $\operatorname{tr}{(A\rho)}$ with
bias at most $\theta$ and $\sigma={\cal O}(1)$ can be implemented using
$\widetilde{{\cal O}}(\alpha)$ uses of $U$ and $U^{\dagger}$ and
$\widetilde{{\cal O}}_{\frac{1}{\theta}}(\alpha)$ elementary operations.

####  2.1.5 Computational complexity

When discussing the computational complexity of quantum algorithms we normally
express the cost in terms of the number of calls to some input oracle. Unless
otherwise specified, the gate complexity is at most a poly-logarithmic factor
larger than the stated oracle complexity. The meaning of “input oracle access”
depends on the input model:

  * •

For the sparse-oracle access model, it refers to a query to the oracle
describing $C/\|C\|_{F}$.

  * •

For the QRAM model, it refers to the number of accesses to QRAM. A QRAM of
size ${\cal O}\left(ns\log^{2}(n)\right)$ is sufficient for our algorithms,
and in particular, we only need classical write access to the QRAM, i.e., we
do not write in superposition.

It is straightforward to translate each of these oracle costs into a running
time in the standard gate model without QRAM, by considering the cost of
implementing each oracle.

##  3 Hamiltonian Updates

In this section, we present the algorithm from [10] and relevant results
required to prove its convergence and analyze its cost.

###  3.1 Convex Feasibility Problems

In order to avoid any normalization issues for the problems that arise over
the course of our IR scheme, we deviate slightly from [10] and renormalize the
problem (3) using the Frobenius norm of the cost matrix rather than use its
operator norm:

find $\displaystyle X$ (5) s.t.
$\displaystyle\operatorname{tr}{\left(\frac{C}{\|C\|_{F}}X\right)}\geq\gamma-\epsilon$
$\displaystyle\sum_{i\in[n]}\left|\langle
i|X|i\rangle-\frac{1}{n}\right|\leq\epsilon$
$\displaystyle\operatorname{tr}{(X)}=1,\quad X\succeq 0.$

The relaxed renormalized SDO problem (5) is a specific example of the convex
optimization problem

$\displaystyle\max$ $\displaystyle f(X)$ (6) s.t. $\displaystyle X\in{\cal
P}_{1}\cap{\cal P}_{2}\cap\cdots\cap{\cal P}_{m},$
$\displaystyle\operatorname{tr}(X)=1,\leavevmode\nobreak\ X\succeq 0,$

where ${\cal P}_{1},\dots,{\cal P}_{m}$ are convex sets.

In this context, the trace constraint enforces normalization, but also allows
us to obtain a bound on the optimal objective value. Letting
$\widetilde{C}=C\|C\|_{F}^{-1}$ and invoking the tracial matrix Hölder
inequality [8], it follows that any $X^{*}$ that solves (6) satisfies the
following relation:

$\left|\operatorname{tr}(\widetilde{C}X^{*})\right|\leq\|\widetilde{C}\|\|X^{*}\|_{\operatorname{tr}}=\|\widetilde{C}\|.$

It is well known in the optimization literature that performing binary search
over the range of values

$\gamma\in\left[-\|\widetilde{C}\|,\|\widetilde{C}\|\right]\subseteq[-1,1]$

that the objective can take reduces the task of solving (6) to solving a
sequence of feasibility problems of the form:

find $\displaystyle X\in{\cal S}^{n}_{+}\cap\\{X:\operatorname{tr}(X)=1\\}$
(7) s.t. $\displaystyle\operatorname{tr}(\widetilde{C}X)\geq\gamma$
$\displaystyle X\in{\cal P}_{1}\cap{\cal P}_{2}\cap\cdots\cap{\cal P}_{m}.$

In particular, $\log(\|\widetilde{C}\|\epsilon^{-1})$ queries to (7) are
sufficient to estimate the optimal objective value of of (6) up to additive
error $\epsilon$.

###  3.2 Solving Convex Feasibility Problems via Hamiltonian Updates

Hamiltonian Updates (HU) is a meta-algorithm for solving convex feasibility
problems of the form (7), adapted from the work of Tsdua, Rätsch and Warmuth
[57] as well as [4, 9, 32, 42]. At a high level, HU can be viewed as a mirror
descent algorithm with the von Neumann entropy as the mirror map. In each
iteration, the method uses certain subroutines to test $\epsilon$-closeness to
convex sets ${\cal P}_{1},{\cal P}_{2},\dots,{\cal P}_{m}$, which we formally
define next.

######  Definition 6 (Definition 2.1 in [10]).

Let ${\cal P}\subset\\{X\in{\cal S}^{n}_{+}:\operatorname{tr}(X)=1\\}$ be a
closed, convex subset of quantum states, and $\widetilde{{\cal
P}}\subset\\{X\in\mathbb{C}^{n\times n}:X=X^{\dagger},\|X\|\leq 1\\}$ be a
closed, convex subset of observables of operator norm at most 1. For
$\epsilon>0$, an $\epsilon$-separation oracle with respect to
$\widetilde{{\cal P}}$ is a subroutine that either accepts a state $\rho$ (in
the sense that observables from $\widetilde{{\cal P}}$ cannot distinguish
$\rho$ from the elements of ${\cal P}$), or provides a normal vector (in the
matrix space) $P$ of a hyperplane that separates $\rho$ from the set ${\cal
P}$ using a test from $\widetilde{{\cal P}}$:

$O_{{\cal P},\epsilon}(\rho)=\begin{cases}\textit{accept}\leavevmode\nobreak\
\rho&\textit{if}\leavevmode\nobreak\ \min_{Y\in{\cal
P}}\max_{P\in\widetilde{{\cal P}}}\operatorname{tr}(P(\rho-Y))\leq\epsilon,\\\
\textit{output}\leavevmode\nobreak\ P\in\widetilde{{\cal
P}}\leavevmode\nobreak\ \textit{s.t.}\leavevmode\nobreak\
\operatorname{tr}(P(\rho-Y))\geq\frac{\epsilon}{2}\leavevmode\nobreak\
\textit{for all}\leavevmode\nobreak\ Y\in{\cal
P}&\textit{otherwise}.\end{cases}$

The authors in [10] point out that the above oracle construction is well
defined, as we can always choose some hyperplane $P\in\widetilde{{\cal P}}$
such that

$\operatorname{tr}{(P(\rho-Y))}\geq\frac{\epsilon}{2},$

holds for all $Y\in{\cal P}$ whenever

$\min_{Y\in{\cal P}}\max_{P\in\widetilde{{\cal
P}}}\operatorname{tr}(P(\rho-Y))>\epsilon.$

From Sion’s min-max theorem [54], it follows that

$\max_{P\in\widetilde{{\cal P}}}\min_{Y\in{\cal
P}}\operatorname{tr}(P(\rho-Y))=\min_{Y\in{\cal P}}\max_{P\in\widetilde{{\cal
P}}}\operatorname{tr}(P(\rho-Y))>\epsilon,$

and hence there exists a hyperplane which separates $\rho$ from ${\cal P}$ by
$\epsilon$. By relaxing the requirement to $\frac{\epsilon}{2}$-separation,
the algorithm is able to reconcile with the errors that result from
approximating quantities computed with $\rho$, or estimating its entries.

The Hamiltonian Updates (HU) algorithm of Brandão et al. [10] is provided in
full detail in Algorithm 1. The algorithm takes as input the precision
parameter $\epsilon$, and $m$ $\epsilon$-separation oracles
$O_{1,\epsilon},O_{2,\epsilon},\dots,O_{m,\epsilon}$. In the initialization
steps, the starting point is defined to be the maximally mixed state
$\rho\leftarrow n^{-1}I$. This is critical to ensuring the convergence of
mirror descent-based approaches such as Algorithm 1 and the works in [4, 9,
32, 42, 57]; initialization to the maximally mixed state ensures that the
quantum relative entropy between any feasible state and the initial state is
bounded by $\log(n)$ (see, e.g., Theorem 11.8 pt. 2 [51]), and is reduced at
every iteration. Consequently, Algorithm 1 terminates in a finite number of
iterations.

As noted in [10], how we define $\widetilde{{\cal P}}$ determines the number
of closeness conditions that need to be tested. By using the Gibbs state
change of variables, we do not need to test if our candidate solution is trace
normalized or positive semidefinite; any Gibbs state

$\rho_{H}=\frac{\exp(-H)}{\operatorname{tr}(\exp(-H))}$

is an element of the set $\\{X\in{\cal S}^{n}_{+}:\operatorname{tr}(X)=1\\}$
by definition. Our task therefore reduces to finding a $\log(n)$-qubit mixed
state state $\rho$ which is $\epsilon$-close to the convex sets ${\cal P}_{i}$
that arise from any other constraints included in the feasibility problem. At
each iteration, $\epsilon$-closeness is tested by querying
$\epsilon$-separation oracles which are constructed using observables in
$\widetilde{{\cal P}}_{i}$. If each of our oracles accepts the candidate
state, the algorithm terminates and reports $(\rho,H)$ as an
$\epsilon$-precise solution. Otherwise, upon detecting infeasibility the
matrix exponent is updated to penalize the infeasible directions using the
rule

$H\leftarrow H+\frac{\epsilon}{16}P,$

where $P$ is a normal vector in the matrix space of a hyperplane that
witnesses infeasibility.

Input: Query access to $m$ $\epsilon$-separation oracles
$O_{1,\epsilon}(\cdot),\dots,O_{m,\epsilon}(\cdot)$

Initialize $\rho\leftarrow n^{-1}I$ and $H\leftarrow\mathbf{0}^{n\times n}$

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

for  _$i=1,\dots,m$_ do

if  _$O_{i,\epsilon}(\rho)=P$_ then

$H\leftarrow H+\frac{\epsilon}{16}P$

$\rho\leftarrow\frac{\exp(-H)}{\operatorname{tr}(\exp(-H))}$

break

end if

end for

return $(\rho,H)$ and exit

end for

Algorithm 1 Hamiltonian Updates for Convex Feasibility Problems

The following result establishes the iteration complexity of Algorithm 1.

######  Theorem 3 (Theorem 2.1 in [10]).

Algorithm 1 requires at most $T=\lceil 64\log(n)\epsilon^{-2}\rceil+1$
iterations to certify that (7) is infeasible or output a state $\rho$
satisfying

$\text{for all}\leavevmode\nobreak\ 1\leq i\leq m:\leavevmode\nobreak\
\max_{P_{i}\in\widetilde{{\cal P}}_{i}}\min_{Y_{i}\in{\cal
P}_{i}}\operatorname{tr}(P_{i}(\rho-Y_{i}))\leq\epsilon.$

Note that Theorem 3 applies to _any_ convex feasibility problem (on density
operators, i.e., trace-normalized positive semidefinite matrices) for which we
have separation oracles as outlined in Definition 6. This is crucial for the
development of an iterative refinement scheme.

There is an important distinction with respect to output across the models of
computation we study. A classical implementation of Algorithm 1 outputs an
explicit description of an $\epsilon$-precise solution $\rho^{*}$ to (5) and
its associated Hamiltonian $H^{*}$, whereas a quantum implementation reports a
real valued vector $y\in\mathbb{R}^{2}$ along with a diagonal matrix $D$ (with
$\|D\|\leq 1$) such that $H^{*}=y_{1}\frac{C}{\|C\|_{F}}+y_{2}D$. The vector
$y=(y_{1},y_{2})^{\top}$ is the state preparation pair of $\rho^{*}$, in
particular:

$\rho^{*}=\frac{\exp\left(-\left(y_{1}\frac{C}{\|C\|_{F}}+y_{2}D\right)\right)}{\operatorname{tr}\left[\exp\left(-\left(y_{1}\frac{C}{\|C\|_{F}}+y_{2}D\right)\right)\right]},$

and we refer to this type of output as a state preparation pair description of
$\rho$. This choice of output is used in all quantum SDO solvers based on
Gibbs sampling techniques (see, e.g., [9, 10, 11, 60, 61]), and is motivated
by the fact that it is difficult to develop quantum algorithms that are
substantially faster than classical algorithms if we still have to output each
entry of the solution (an $n\times n$ matrix).

The Gibbs sampling approaches that we apply later exhibit a cost that depends
on a norm bound for $y$. Observe that we initialize $y$ to the all zeros
vector of appropriate dimension, and in every iteration, at most one entry of
$y$ changes by a magnitude of $\frac{\epsilon}{16}$ (specifically, an entry
$y_{i}$, where the oracle $O_{i,\epsilon}$ has detected infeasibility). As a
consequence, the vector $y$ satisfies the inequality

$\left\|y^{(t+1)}-y^{(t)}\right\|\leq\frac{\epsilon}{16}$ (8)

for each iteration $t$. In view of the iteration bound for Algorithm 1
provided in Theorem 3, it is easy to see that for any $y$ obtained from
Algorithm 1 we have

$\|y\|_{1}\leq\lceil
64\log(n)\epsilon^{-2}\rceil\left\|y^{(t+1)}-y^{(t)}\right\|\leq\lceil
64\log(n)\epsilon^{-2}\rceil\frac{\epsilon}{16}\leq 4\log(n)\epsilon^{-1}.$
(9)

To instantiate the algorithm to solve problem (3) we need to choose the sets
${\cal P}_{i}$, and provide separation oracles for them. This is what we do in
the following section.

####  3.2.1 Oracle Construction

The goal of Hamiltonian Updates is to solve, for fixed $\gamma\in[-1,1]$, the
following feasibility problem:

find $\displaystyle\rho\in\\{X\in{\cal
S}^{n}_{+}:\operatorname{tr}(X)=1\\}\cap{\cal C}_{\gamma}\cap{\cal D}_{n}$
(10) where $\displaystyle{\cal
C}_{\gamma}=\left\\{X:\operatorname{tr}\left(\frac{C}{\|C\|_{F}}X\right)\geq\gamma\right\\},$
$\displaystyle{\cal D}_{n}=\left\\{X:\langle
i|X|i\rangle=\frac{1}{n},i\in[n]\right\\}.$

One can observe that the set ${\cal C}_{\gamma}$ constitutes a halfspace,
while ${\cal D}_{n}$ is an affine space of codimension $n$. The sets of
observables for ${\cal C}_{\gamma}$ and ${\cal D}_{n}$ are given by
$\widetilde{{\cal C}}_{\gamma}$ and $\widetilde{{\cal D}}_{n}$ respectively,
with

$\widetilde{{\cal C}}_{\gamma}=\\{-C\|C\|_{F}^{-1}\\},\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ \widetilde{{\cal
D}}_{n}=\\{D\in\mathbb{R}^{n\times n}:\|D\|\leq 1,\leavevmode\nobreak\
D\leavevmode\nobreak\ \text{is diagonal}\\}.$

As noted in [10], it follows

$\max_{P\in\widetilde{{\cal C}}_{\gamma}}\min_{Y\in{\cal
C}_{\gamma}}\operatorname{tr}(P(\rho-Y))\leq\epsilon\iff-\operatorname{tr}\left(C\|C\|_{F}^{-1}(\rho-Y)\right)\leq\epsilon\quad\text{for
some}\leavevmode\nobreak\ Y\in{\cal C}_{\gamma},$

which in turn implies
$\operatorname{tr}{(C\|C\|_{F}^{-1}\rho)}\geq\gamma-\epsilon$.

Given the structure of ${\cal C}_{\gamma}$ and ${\cal D}_{n}$, the authors in
[10] suggest the following two separation oracles:

$\displaystyle O_{{\cal C}_{\gamma}}:\leavevmode\nobreak\ $
$\displaystyle\text{compute an approximation}\leavevmode\nobreak\
\tilde{c}\leavevmode\nobreak\ \text{of}\leavevmode\nobreak\
\operatorname{tr}\left(C\|C\|_{F}^{-1}\rho\right)\leavevmode\nobreak\ \text{up
to additive error}\leavevmode\nobreak\ \frac{\epsilon}{4}.\leavevmode\nobreak\
\text{Check if}\leavevmode\nobreak\
\tilde{c}\geq\gamma-\frac{3\epsilon}{4}\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ $
$\displaystyle\text{output}\leavevmode\nobreak\
P=-C\|C\|_{F}^{-1}\leavevmode\nobreak\ \text{if the inequality is violated}.$
$\displaystyle O_{{\cal D}_{n}}:\leavevmode\nobreak\ $
$\displaystyle\text{compute an approximation}\leavevmode\nobreak\
\tilde{p}\in\mathbb{R}^{n}\leavevmode\nobreak\ \text{of}\leavevmode\nobreak\
p_{i}=\langle i|\rho|i\rangle\leavevmode\nobreak\
\text{satisfying}\leavevmode\nobreak\
\sum_{i=1}^{n}|p_{i}-\tilde{p}_{i}|\leq\frac{\epsilon}{4}.$
$\displaystyle\text{Check if}\leavevmode\nobreak\
\sum_{i=1}^{n}\left|\tilde{p}_{i}-\frac{1}{n}\right|\leq\frac{3\epsilon}{4}\leavevmode\nobreak\
\text{and output}\leavevmode\nobreak\
P=\sum_{i=1}^{n}\left(\mathbb{I}\left\\{\tilde{p}_{i}>\frac{1}{n}\right\\}-\mathbb{I}\left\\{\tilde{p}_{i}<\frac{1}{n}\right\\}\right)\left|{i}\right\rangle\left\langle{i}\right|$
$\displaystyle\text{if the inequality is violated}.$

For any given

$\rho_{H}=\frac{\exp(-H)}{\operatorname{tr}(\exp(-H))},$

the required separation oracles are straightforward to implement on a
classical computer that has access to $\rho_{H}$. Thus, classically we only
need to prepare $\rho_{H}$ once and store it to build the separation oracles.
The next result from [10] establishes that computing an ${\cal
O}(\log(n)\epsilon^{-1})$-degree Taylor series suffices to produce accurate
approximations.

######  Lemma 9 (Lemma 3.2 in [10]).

Fix a Hermitian $n\times n$ matrix $H$, an accuracy $\epsilon$, and let $\ell$
be the smallest even number satisfying $(\ell+1)(\log(\ell+1)-1)\geq
2\|H\|+\log(n)+\log\left(\frac{1}{\epsilon}\right)$. Then, the truncated
matrix exponential $T_{\ell}=\sum_{k=0}^{\ell}\frac{1}{k!}(-H)^{k}$ satisfies

$\left\|\frac{\exp(-H)}{\operatorname{tr}\left(\exp(-H)\right)}-\frac{T_{\ell}}{\operatorname{tr}(T_{\ell})}\right\|_{\operatorname{tr}}\leq\epsilon.$

states [10], which are used to test closeness to the sets ${\cal C}_{\gamma}$
and ${\cal D}_{n}$ via quantum measurements. While in Lemma 9 we bound the
number of required Taylor series steps for computing $\rho$ via a matrix
exponential, in the quantum case we bound the number of copies of $\rho$
required to estimate its diagonal entries and expectation values
$\operatorname{tr}(A\rho)$.

######  Lemma 10.

Fix $\epsilon\in(0,1).$ Let $\rho$ be a $\log(n)$-qubit quantum state and $U$
a $(1,\log(n)+2,\epsilon/(2n))$-block-encoding of $C\|C\|_{F}^{-1}$. Then, we
can implement the oracle $O_{{\cal C}_{\gamma}}$ on a quantum computer given
access to ${\cal O}(\epsilon^{-1})$ copies of a state that is an
$\frac{\epsilon}{8}$-approximation of the input state $\rho$ in trace distance
and ${\cal O}(\epsilon^{-1})$ applications of $U$ and $U^{\dagger}$. 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})$.

###### Proof.

First, note that we can obtain an estimate $\tilde{p}$ of the diagonal
elements of $\rho$ whose total variation distance from $p$ is no more than
$\frac{\epsilon}{8}$ using $\widetilde{{\cal
O}}_{n}\left(n\epsilon^{-2}\right)$ copies of $\rho$ to measure $\rho$ in the
computational basis. Further, provided accesses to $\rho$ and a
$(1,\log(n)+2,\epsilon/(2n))$-block-encoding $U$ of $C\|C\|_{F}^{-1}$, by
Lemma 8, a trace estimator for
$\operatorname{tr}{\left(C\|C\|_{F}^{-1}\rho\right)}$ with bias at most
$\frac{\epsilon}{n}$ can be implemented using $\widetilde{{\cal O}}(1)$ uses
of $U$ and $U^{\dagger}$ and $\widetilde{{\cal O}}_{\frac{n}{\epsilon}}(1)$
elementary operations. From here, applying amplitude estimation using ${\cal
O}(\epsilon^{-1})$ quantum samples (i.e., state preparation unitaries) from
the trace estimator to suffice to compute an approximation
$\operatorname{tr}{\left(C\|C\|_{F}^{-1}\rho\right)}$ up to additive
$\frac{\epsilon}{8}$ to implement $O_{{\cal C}_{\gamma}}$. The rest of the
proof exactly follows the proof of [10, Lemma 3.3]. ∎

We remark that multidimensional phase estimation techniques from [59] could
improve the dependence on $\epsilon^{-1}$ for estimating the diagonal elements
of $\rho$ to linear, which is a factor $\epsilon^{-1}$ better than a näive
application of computational basis measurements. However, in the context of
the iterative refinement scheme we present later, the improvement would only
reduce the amount of constant overhead in the overall running time, and
multidimensional phase estimation has a larger gate complexity (which can be
reduced with QRAM). There are also numerous ways to prepare Gibbs states using
a quantum computer [16, 20, 37, 53, 60, 61, 63]. Following [10], we utilize
the Gibbs sampler from [53] when working with the sparse-access input model,
and for the QRAM input model we consider Gibbs sampling techniques introduced
in [60].

###  3.3 Complexity

Having understood the cost of constructing the oracles in both the classical
and quantum settings, we are now in a position to analyze the complexity
associated with using Algorithm 1 to obtain solutions to (5) and
approximations to (3). Relevant to this discussion is the following result,
which imposes precision requirements on solving (3) to an additive error of
the order ${\cal O}\left(n\|C\|_{F}\epsilon\right)$ using Algorithm 1.

######  Proposition 5 (Proposition 3.1 in [10]).

Let $\rho$ be an $\epsilon^{4}$-accurate solution to the relaxed SDO problem
(5) with input matrix $C$. Let
$\gamma_{\epsilon^{4}}=\operatorname{tr}{(C\rho)}$ be the value attained by
$\rho$. Then, there is a quantum state $\rho^{*}$ at trace distance ${\cal
O}(\epsilon)$ of $\rho$ such that $n\rho^{*}$ is a feasible point of SDO
problem (3). In particular

$\left|\gamma_{\epsilon^{4}}n\|C\|_{F}-\operatorname{tr}\left(n\rho^{*}C\right)\right|={\cal
O}\left(n\|C\|_{F}\epsilon\right).$

Moreover, it is possible to construct $\rho^{*}$ in time ${\cal O}(n^{2})$
given the entries of $\rho$.

We do not provide a proof of this result here, as later we will provide an
improved approximation guarantee and a proof of the improved statement.

####  3.3.1 Classical running time

Using Lemma 9 in combination with Theorem 3, we can bound the running time
required to solve (5) to additive error $\epsilon$ using a classical
implementation of Algorithm 1.

######  Proposition 6.

Suppose that $C$ has row sparsity $s$. Then, the classical cost of solving (5)
up to additive error $\epsilon$ using Algorithm 1 is ${\cal
O}\left(\min\\{n^{2}s,n^{\omega}\\}\log^{2}(n)\epsilon^{-3}\right)$.

###### Proof.

The result follows directly from the proof of Corollary 3.1 in [10], but we
repeat the argument here for completeness.

First, observe that over the course of the iterations $t=0,\dots,T$, the
operator norms $\|H^{(t)}\|$ do not become prohibitively large. This follows
from initializing $H^{(0)}=\mathbf{0}^{n\times n}$, and that by (8), the
inequality

$\left\|H^{(t+1)}-H^{(t)}\right\|\leq\frac{\epsilon}{16}\left\|P^{(t)}\right\|\leq\frac{\epsilon}{16}$

holds for all $t$. By Theorem 3, Algorithm 1 requires at most $T=\lceil
64\log(n)\epsilon^{-2}\rceil$ iterations, which implies $\|H^{(t)}\|\leq
4\log(n)\epsilon^{-1}$ for all $t$.

By Lemma 9, it suffices to compute ${\cal O}(\log(n)\epsilon^{-1})$ steps of
the Taylor series corresponding to $\exp(-H^{(t)})$ in order to obtain a
matrix $\tilde{\rho}^{(t)}$ that is at most a trace distance of
$\frac{\epsilon}{4}$ from $\rho^{(t)}$. Moreover, given that $H^{(t)}$ is
defined as a linear combination of $C\|C\|_{F}^{-1}$ with a diagonal matrix,
matrix multiplication involving $H^{(t)}$ can be carried out in ${\cal
O}(\min\\{n^{2}s,n^{\omega}\\})$ arithmetic operations. Given classical access
to $\tilde{\rho}^{(t)}$, the diagonal constraints comprising ${\cal D}_{n}$
can be checked in time ${\cal O}(n)$, whereas computing
$\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\tilde{\rho}^{(t)}\right)$ requires
${\cal O}(ns)$ arithmetic operations. Thus, the dominant operation at each
iteration is computing the matrix exponential and the classical per-iteration
cost of Algorithm 1 is given by

${\cal O}\left(\min\\{n^{2}s,n^{\omega}\\}\log(n)\epsilon^{-1}\right).$

Taking into account the iteration bound ${\cal O}(\log(n)\epsilon^{-2})$
provided in Theorem 3, we arrive at an overall running time of

${\cal O}\left(\min\\{n^{2}s,n^{\omega}\\}\log^{2}(n)\epsilon^{-3}\right).$

The proof is complete. ∎

The next corollary from [10] follows from Proposition 5 in the context of the
previous result, and provides the overall running time of Algorithm 1 to solve
(3) to additive error ${\cal O}\left(n\|C\|_{F}\epsilon\right)$ in the
classical setting.

######  Corollary 2.

Suppose that $C$ has row-sparsity $s$. Then, the classical cost of solving (3)
up to an additive error ${\cal O}\left(n\|C\|_{F}\epsilon\right)$ using
Algorithm 1 is ${\cal
O}\left(\min\\{n^{2}s,n^{\omega}\\}\log^{2}(n)\epsilon^{-12}\right)$.

###### Proof.

By Proposition 6, Algorithm 1 requires time

${\cal
O}\left(\min\\{n^{2}s,n^{\omega}\\}\log^{2}(n)\tilde{\epsilon}^{-3}\right),$

to solve (5) up to additive error $\tilde{\epsilon}$. In order to satisfy the
approximation guarantee for (3) given in Proposition 5, it suffices to solve
(5) to error $\tilde{\epsilon}=\epsilon^{4}$. Plugging in this value for the
precision parameter, the total cost required to solve (3) up to an additive
error ${\cal O}\left(n\|C\|_{F}\epsilon\right)$ using Algorithm 1 is

${\cal
O}\left(\min\\{n^{2}s,n^{\omega}\\}\log^{2}(n)\tilde{\epsilon}^{-3}\right)={\cal
O}\left(\min\\{n^{2}s,n^{\omega}\\}\log^{2}(n)(\epsilon^{4})^{-3}\right)={\cal
O}\left(\min\\{n^{2}s,n^{\omega}\\}\log^{2}(n)\epsilon^{-12}\right).$

∎

####  3.3.2 Quantum running time

Combining the sampling requirements provided in Lemma 10 with the cost of
preparing a single Gibbs state and the iteration bound from Theorem 3 gives
the complexity of Algorithm 1 when run on a quantum computer. However, Gibbs
samplers based on the block-encoding framework depend only poly-
logarithmically on the inverse precision, therefore they are exponentially
faster (in the parameter $\epsilon^{-1}$) compared to the Gibbs sampling
algorithm from [53] utilized in [10]. It thus makes sense to analyze the
running time in the more efficient model. This will require an efficient data
structure for storing $y$ so that we can efficiently prepare linear
combinations of block-encodings.

######  Lemma 11 (Lemma 15 in [60]).

There is a data structure that can store an $m$-dimensional $\chi$-sparse
vector $y$ with $\theta$-precision using a QRAM of size $\widetilde{{\cal
O}}_{\frac{m}{\theta}}(\chi)$. Furthermore:

  * •

Given a classical ${\cal O}(1)$-sparse vector, adding it to the stored vector
has classical cost $\widetilde{{\cal O}}_{\frac{m}{\theta}}(1)$.

  * •

Given that $\beta\geq\|y\|_{1}$, we can implement a (symmetric)
$(\beta,\widetilde{{\cal O}}_{\frac{m}{\theta}}(1),\theta)$-state preparation
pair for $y$ with $\widetilde{{\cal O}}_{\frac{m}{\theta}}(1)$ queries to the
QRAM.

######  Corollary 3 (Corollary 16 in [60]).

Suppose $A_{1},\dots,A_{m}$ are Hermitian matrices with operator norm at most
1, and that $y\in\mathbb{R}^{m}$ satisfies $\|y\|_{1}\leq\beta$. Having access
to the above data structure for $y$, we can prepare one copy of the Gibbs
state

$\rho=\frac{\exp\left(-\sum_{i=1}^{m}y_{i}A_{i}\right)}{\operatorname{tr}\left(\exp\left(-\sum_{i=1}^{m}y_{i}A_{i}\right)\right)}$

using $\widetilde{{\cal O}}_{\theta}(\sqrt{n}\alpha\beta)$ accesses to the
data structure for $y$ and block-encodings of $A_{1},\dots,A_{m}$.

We can now use Corollary 3 in combination with results from Sections 2.1.3 and
2.1.4 to establish the running time of Algorithm 1 in the QRAM input model.

######  Proposition 7.

Let $\frac{C}{\|C\|_{F}}\in{\cal S}^{n}$ be stored in QRAM. Then, the
complexity of solving (5) up to additive error $\epsilon$ with Algorithm 1
using the QRAM input model is

$\widetilde{{\cal O}}_{\frac{n}{\epsilon}}\left(n^{1.5}\epsilon^{-5}\right).$

Here, the complexity corresponds to the number of accesses to the QRAM.

###### Proof.

Given that $\frac{C}{\|C\|_{F}}$ is stored in QRAM, Lemma 3(ii) asserts that
when constructing a block-encoding of $\frac{C}{\|C\|_{F}}$, one can set the
subnormalization factor to be
$\alpha_{C}=\left\|\frac{C}{\|C\|_{F}}\right\|_{F}=1$. Hence, one can
construct a $(1,\log(n)+2,\epsilon/(2n))$-block-encoding of $C\|C\|_{F}^{-1}$
in time $\widetilde{{\cal O}}_{\frac{n}{\epsilon}}(1)$.

Next, recall that in iteration $t\in[T]$ of Algorithm 1, our Hamiltonian is
defined as

$H^{(t)}=y_{1}^{(t)}\frac{C}{\|C\|_{F}}+y_{2}^{(t)}D^{(t)},$

where $D^{(t)}$ is a diagonal matrix with the diagonal entries taking value
$-1$, $0$ or $1$. The diagonal elements of $D$ change in each iteration, and
therefore, a new $D$ must be block-encoded in each iteration. For this, we use
the QRAM model described in Section 2.1.2, which allows for insertions to be
made in time $\widetilde{{\cal O}}_{n}(1)$ to keep the cost of this step
negligible. Provided a classical description of $D$, we can store $D$ in the
QRAM in time ${\cal O}(n\log(n))$. Applying Lemma 4, a
$(1,\log(n)+3,\epsilon)$-block-encoding of $D^{(t)}$ can be constructed in
time $\widetilde{{\cal O}}_{\frac{n}{\epsilon}}(1)$.

In an earlier discussion we saw that any $y$ obtained from a call to Algorithm
1 will satisfy $\|y\|_{1}=\widetilde{{\cal O}}_{n}(\epsilon^{-1})$ if we call
Algorithm 1 using precision $\epsilon$ (see, e.g., equation (9)). Hence, an
application of Corollary 3 with $\beta=\widetilde{{\cal
O}}_{n}(\epsilon^{-1})$ implies that we can prepare one copy of our Gibbs
state using

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

accesses to the data structure for $y$ and the block-encodings of
$C\|C\|_{F}^{-1}$ and $D$, where $\alpha$ is defined as the maximum over the
subnormalization factors used to block-encode $\frac{C}{\|C\|_{F}}$ and $D$.
Since $\alpha=\max\\{\alpha_{C},\alpha_{D}\\}=1$, it follows

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

Now, one can see from Lemma 10 that the cost of constructing $O_{{\cal
D}_{n}}$ dominates that of constructing $O_{{\cal C}_{\gamma}}$. Noting that
$O_{{\cal D}_{n}}$ can be implemented using ${\cal O}(n\epsilon^{-2})$ copies
of a state that is an $\frac{\epsilon}{8}$-approximation of the input state
$\rho$ in trace distance and its inverse, the per-iteration cost of Algorithm
1 in the QRAM input model is given by

$\widetilde{{\cal O}}_{\frac{n}{\epsilon}}\left(n^{1.5}\epsilon^{-3}\right).$

Factoring in the iteration bound of $\widetilde{{\cal O}}_{n}(\epsilon^{-2})$
from Theorem 3, it follows that when provided access to QRAM, Algorithm 1
solves (5) up to additive error $\epsilon$ using

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

accesses to the QRAM. The proof is complete. ∎

######  Corollary 4.

Let $\frac{C}{\|C\|_{F}}\in{\cal S}^{n}$ be stored in QRAM. Then, the
complexity of solving (3) up to additive error ${\cal O}(n\|C\|_{F}\epsilon)$
with Algorithm 1 using the QRAM input model is

$\widetilde{{\cal O}}_{\frac{n}{\epsilon}}\left(n^{1.5}\epsilon^{-20}\right).$

Here, the complexity corresponds to the number of accesses to the QRAM.

###### Proof.

By Proposition 7, Algorithm 1 requires

$\widetilde{{\cal
O}}_{\frac{n}{\tilde{\epsilon}}}\left(n^{1.5}\tilde{\epsilon}^{-5}\right),$

accesses to the QRAM to solve (5) up to additive error $\tilde{\epsilon}$. In
order to satisfy the approximation guarantee for (3) given in Proposition 5,
it suffices to solve (5) to error $\tilde{\epsilon}=\epsilon^{4}$. Plugging in
this value for the precision parameter, the total cost required to solve (3)
up to an additive error ${\cal O}\left(n\|C\|_{F}\epsilon\right)$ using
Algorithm 1 is

$\widetilde{{\cal
O}}_{\frac{n}{\tilde{\epsilon}}}\left(n^{1.5}\tilde{\epsilon}^{-5}\right)=\widetilde{{\cal
O}}_{\frac{n}{\epsilon}}\left(n^{1.5}(\epsilon^{4})^{-5}\right)=\widetilde{{\cal
O}}_{\frac{n}{\epsilon}}\left(n^{1.5}\epsilon^{-20}\right).$

The proof is complete. ∎

Corollary 4 establishes that utilizing Gibbs samplers and trace estimators
based on the block-encoding framework for our oracle construction in Algorithm
1 leads to an

${\cal
O}\left(\sqrt{s}^{1+o(1)}\epsilon^{-8+o(1)}\exp\left(1.6\sqrt{\log(\epsilon^{-4})}\right)\right)$

speedup over the running time result provided in [10, Corollary 3.2] when
applied to solving (3). Yet, the costly accuracy requirements for the rounding
procedure (see, e.g., Proposition 5) lead to a prohibitive scaling in the
inverse precision for the overall running time. Given the advantageous
dependence on the dimension, as compared to classical algorithms, we study how
to improve the dependence on the precision parameter. This is discussed next.

##  4 Iterative Refinement for SDO approximations of QUBOs

In this section, we introduce an iterative refinement method for obtaining
accurate solutions to the renormalized relaxed SDO problem (5), that at a high
level can be viewed as solving a series of problems related to the feasibility
problem (10) associated with (5). We then discuss how to test
$\epsilon$-closeness to the convex sets which comprise the feasible regions of
the intermediate refining problems before presenting our algorithm in full
detail. We conclude the section by proving our algorithm’s correctness and
iteration complexity, and use these results to provide an improved
approximation guarantee.

###  4.1 The refining problem

To develop an iterative refinement scheme for (5), we need to design a problem
whose solution can be used to improve the quality of solutions to (5). Suppose
we run Algorithm 1 and obtain an $\epsilon$-precise solution $\hat{\rho}$ to
(5). Letting
$\hat{\gamma}=\operatorname{tr}\left(C\|C\|_{F}^{-1}\hat{\rho}\right)$,
$\hat{\rho}$ must satisfy

$\displaystyle\operatorname{tr}\left(C\|C\|_{F}^{-1}\hat{\rho}\right)=\hat{\gamma}$
$\displaystyle\geq\gamma-\epsilon,$ $\displaystyle\sum_{i=1}^{n}\left|\langle
i|\hat{\rho}|i\rangle-\frac{1}{n}\right|$ $\displaystyle\leq\epsilon.$

In refining our solution to (5), we should aim to reduce the trace distance to
the maximally mixed state $n^{-1}I$, while also improving the precision to
which the optimal objective value is approximated. Thus, an improved solution
$\rho^{\prime}$ should obey

$\displaystyle\operatorname{tr}\left(C\|C\|_{F}^{-1}\rho^{\prime}\right)$
$\displaystyle\geq\gamma-\epsilon^{\prime},$
$\displaystyle\sum_{i=1}^{n}\left|\langle
i|\rho^{\prime}|i\rangle-\frac{1}{n}\right|$
$\displaystyle\leq\epsilon^{\prime},$

with $\epsilon^{\prime}<\epsilon$. The basic idea behind constructing the
refining problem is to use our current solution $\hat{\rho}$ to first shift
the renormalized relaxed SDO problem (5) to the origin, and then scale the
shifted problem back to the domain of the original problem. In particular, we
solve a series of problems related to the feasibility problem (10).

Let $\varepsilon\in\mathbb{R}^{n}$ be a vector whose elements are the
residuals along the diagonal $\varepsilon_{i}=\hat{\rho}_{ii}-\frac{1}{n}$ for
$i\in[n]$, and $\eta\geq 1$ to be a scalar defined as

$\eta=\frac{1}{\max\left\\{\gamma-\operatorname{tr}\left(C\|C\|_{F}^{-1}\hat{\rho}\right),\sum_{i=1}^{n}|\varepsilon_{i}|\right\\}}=\frac{1}{\max\left\\{\gamma-\operatorname{tr}\left(C\|C\|_{F}^{-1}\hat{\rho}\right),\left\|\hat{\rho}-n^{-1}I\right\|_{\operatorname{tr}}\right\\}}.$

Using these quantities, the refining problem is given by:

find $\displaystyle\rho^{r}\in\\{X\in{\cal
S}^{n}_{+}:\operatorname{tr}(X)=1\\}\cap{\cal
C}_{\eta(\gamma-\hat{\gamma})}\cap{\cal D}_{\eta\varepsilon}$ (11) where
$\displaystyle{\cal
C}_{\eta(\gamma-\hat{\gamma})}=\left\\{X:\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\left(Q\circ
X\right)\right)\geq\eta(\gamma-\hat{\gamma})\right\\},$ $\displaystyle{\cal
D}_{\eta\varepsilon}=\left\\{X:\langle
i|X|i\rangle=\eta|\varepsilon_{i}|,\leavevmode\nobreak\ \forall
i\in[n]\right\\},$

where $Q\in{\cal S}^{n}$ is a matrix whose diagonal elements are chosen such
that for any $X\in{\cal D}_{\eta\varepsilon}$, we have

$(Q\circ X)_{ii}=\text{sign}(-\varepsilon_{i})\eta|\varepsilon_{i}|$

for $i\in[n]$. Further details and requirements on the structure of $Q$ are
specified later in this section. We refer to solutions $\rho^{r}$ to (11) as
refining solutions, which we use to update our current solution $\hat{\rho}$
to (5).

The set ${\cal D}_{\eta\varepsilon}$ is comprised of the diagonal constraints

$\langle i|X|i\rangle=\eta|\varepsilon_{i}|,\quad\forall i\in[n],$

and similar to ${\cal D}_{n}$, is an affine space with codimension $n$. Our
use of the absolute value function of the residuals and scaling by $\eta$
ensures the viability of applying Gibbs sampling techniques to solve the
refining problem (11); the diagonal terms of any density matrix must be
nonnegative and sum to 1. Whenever

$\sum_{i=1}^{n}|\varepsilon_{i}|>\gamma-\operatorname{tr}\left(C\|C\|_{F}^{-1}\hat{\rho}\right),$

then $\eta\|\varepsilon\|_{1}=1$, and the parameter $\eta$ therefore scales
the shifted problem back to the space of the $\log(n)$-qubit mixed states,
ensuring that any solution $\rho^{r}$ to (11) is indeed a (trace normalized)
Gibbs state.

On the other hand, should it be the case that

$\sum_{i=1}^{n}|\varepsilon_{i}|\leq\gamma-\operatorname{tr}\left(C\|C\|_{F}^{-1}\hat{\rho}\right),$

then for any $X\in{\cal D}_{\eta\varepsilon}$ we have
$\operatorname{tr}(X)\leq 1$, rather than $\operatorname{tr}(X)=1$. Our primal
SDO oracle in Algorithm 1 solves feasibility problems in which the trace upper
bound is tight, i.e., $\operatorname{tr}(X)=1$. The authors in [60] note that
this can be dealt with adding one extra variable $w$ such that

$\bar{\rho}^{r}:=\begin{bmatrix}\rho^{r}&0\\\ 0&w\end{bmatrix}.$

Then, $\operatorname{tr}{(\bar{\rho}^{r})}=1$ and $\bar{\rho}^{r}\succeq 0$
imply that $\operatorname{tr}(\rho^{r})\leq 1$, and as a result we obtain an
SDO problem that is equivalent to (11). Since we know exactly the amount of
subnormalization, we can also get rid of the extra variable in subsequent
calculations and rescale the trace back to 1 when necessary (e.g., when
combining solutions from multiple iterative refinement iterations for trace
estimations). Crucially, using the input models described in Section 2.1,
these modifications do not introduce more than constant overhead in the
overall complexity, as the problem data in this case is simply given by

$\mkern 1.5mu\overline{\mkern-1.5muC\mkern-1.5mu}\mkern
1.5mu=\begin{bmatrix}\frac{C}{\|C\|_{F}}&0\\\ 0&0\end{bmatrix},\quad\mkern
1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern 1.5mu=\begin{bmatrix}Q&0\\\
0&0\end{bmatrix},$

with $(\mkern 1.5mu\overline{\mkern-1.5muC\mkern-1.5mu}\mkern 1.5mu,\mkern
1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern 1.5mu)\in{\cal
S}^{n+1}\times{\cal S}^{n+1}$.

The Hadamard product $Q\circ\rho^{r}$ that appears in the definition of ${\cal
C}_{\eta(\gamma-\hat{\gamma})}$ is required for similar reasons; properly
setting $Q$ allows us to drive the trace distance to the maximally mixed state
to zero using the solutions to the refining problem. Later, in Section 4.3 we
demonstrate that this can be achieved by updating the current solution
$\hat{\rho}$ using the rule

$\hat{\rho}=\hat{\rho}+\frac{1}{\eta}Q\circ\rho^{r},$ (12)

with a suitable choice for $Q$ being

$Q=(ee^{\top}-I)+\operatorname{diag}{(\text{sign}(-\varepsilon))}=\begin{pmatrix}\text{sign}(-\varepsilon_{1})&1&\ldots&1\\\
1&\text{sign}(-\varepsilon_{2})&\ddots&\vdots\\\ \vdots&\ddots&\ddots&1\\\
1&\ldots&1&\text{sign}(-\varepsilon_{n})\end{pmatrix}.$ (13)

Choosing $Q$ in this manner also implies that the Hadamard product $Q\circ A$
can be carried out classically using ${\cal O}(n)$ arithmetic operations for
any $A\in\mathbb{R}^{n\times n}$, as the element-wise products
$Q_{ij}A_{ij}=A_{ij}$ for $i\neq j$. Similarly, updating $Q$ at each iterate
only requires updating its diagonal elements, an ${\cal O}(n)$ operation.

It is important to note that the update we propose in (12) does not preserve
positive semidefiniteness in general. However, later in our analysis, we
demonstrate that the eigenvalues of the final solution $\hat{\rho}$ are only
slightly negative in the worst case, i.e.,
$\lambda_{\min}(\hat{\rho})\geq-\delta$ for a small constant $\delta$; one can
restore positive semidefiniteness by adding $\delta$ to the diagonal elements
of the final solution, and we renormalize by $(1+n\delta)$ to obtain unit
trace. It turns out that the trace distance from resulting matrix to
$\hat{\rho}$ is bounded by the final precision tolerance parameter of our
refining scheme. We show that these modifications required to restore positive
semidefiniteness have only a mild (in fact, constant) impact on feasibility.
To do this, we will bound the eigenvalues of $Q$. We first state a special
instance of Weyl’s inequality.

######  Lemma 12.

Suppose that $A\in\mathbb{R}^{n\times n}$ and $B\in\mathbb{R}^{n\times n}$ are
Hermitian matrices. Then

$\lambda_{\min}(A+B)\geq\lambda_{\min}(A)+\lambda_{\min}(B).$

Using the preceding lemma, the following result bounds the minimum eigenvalue
of $Q$.

######  Lemma 13.

Suppose that $Q\in{\cal S}^{n}$ is defined according to Equation (13). Then,
$\lambda_{\min}(Q)\geq-2$.

###### Proof.

Let $A=(ee^{\top}-I)$ and
$B=\operatorname{diag}{(\text{sign}(-\varepsilon))}$, such that $Q=A+B$. Now,
it can be easily seen from the definition of $A$ that $A+I$ is an all-ones
matrix of dimension $n$. Upon performing row-reduction (via, e.g., Guassian
elimination) on $A$, it is trivial to observe that the resulting row-echelon
form will have $n-1$ zero rows, and as a consequence, $A$ has the eigenvalue
$-1$, repeated (at least) $n-1$ times. Further, since
$\operatorname{tr}{(A)}=0$, the other eigenvalue is $n-1$. Therefore, we have
$\lambda_{\min}(A)\geq-1$. On the other hand, $B$ is a diagonal matrix whose
diagonal elements can take value $-1$, $0$, or $1$, from which
$\lambda_{\min}(B)\geq-1$ readily follows.

Applying Lemma 12, we obtain

$\lambda_{\min}(Q)=\lambda_{\min}(A+B)\geq\lambda_{\min}(A)+\lambda_{\min}(B)\geq-2.$

The proof is complete. ∎

###  4.2 Oracle construction for the refining problem

In order to construct separation oracles for testing closeness to ${\cal
C}_{\eta(\gamma-\hat{\gamma})}$, we rely on the following result.

######  Lemma 14.

Let $E$, $F$ and $G\in{\cal S}^{n}$. We have

$\operatorname{tr}{\left(G(E\circ F)\right)}=\operatorname{tr}{\left((E\circ
G)F\right)}.$

###### Proof.

Applying Lemma 1 with $m=n$, we have

$\displaystyle\left[(E\circ F)G\right]_{ii}=\left[(E\circ
G)F\right]_{ii}\quad\forall i\in[n].$

Note that we have dropped the transpose terms, as $E$, $F$ and $G$ are
symmetric matrices, and hence, so are $E\circ F$ and $E\circ G$. It follows

$\operatorname{tr}{\left(G(E\circ F)\right)}=\operatorname{tr}{\left((E\circ
F)G\right)}=\sum_{i\in[n]}\left[(E\circ
F)G\right]_{ii}=\sum_{i\in[n]}\left[(E\circ
G)F\right]_{ii}=\operatorname{tr}{\left((E\circ G)F\right)}.$

∎

In addition to $Q\in{\cal S}^{n}$, we also require
$\max_{i,j\in[n]}\\{|Q_{ij}|\\}\leq 1$ to avoid any normalization issues with
respect to $Q\circ\frac{C}{\|C\|_{F}}$. Note that defining of $Q$ according to
equation (13) satisfies both of these properties trivially, as each of the
diagonal elements are $1$, $0$, or $-1$, while the off-diagonal elements are
all set to 1. This idea is formalized next.

######  Lemma 15.

Let $A,Q\in{\cal S}^{n}$ be matrices satisfying
$\max_{i,j\in[n]}\\{|Q_{ij}|\\}\leq 1$ and $\|A\|_{F}\leq 1$. Then,

$\left\|Q\circ A\right\|\leq\left\|Q\circ A\right\|_{F}\leq 1.$

###### Proof.

Under the stated conditions for $Q$, it follows

$\displaystyle\left\|Q\circ
A\right\|_{F}^{2}=\sum_{i\in[n]}\sum_{j\in[n]}\left(\left[Q\circ
A\right]_{ij}\right)^{2}=\sum_{i\in[n]}\sum_{j\in[n]}\left(Q_{ij}\cdot
A_{ij}\right)^{2}$
$\displaystyle=\sum_{i\in[n]}\sum_{j\in[n]}\left(Q_{ij}\right)^{2}\left(A_{ij}\right)^{2}$
$\displaystyle\leq\sum_{i\in[n]}\sum_{j\in[n]}\left(A_{ij}\right)^{2}=\left\|A\right\|_{F}^{2},$

and applying the square root throughout the above we obtain $\left\|Q\circ
A\right\|_{F}\leq\left\|A\right\|_{F}$. From here, the result follows upon
noting $\left\|A\right\|_{F}\leq 1$ and $\|A\|\leq\|A\|_{F}$ is true for any
$A\in\mathbb{R}^{n\times n}$. ∎

Although the sets ${\cal C}_{\gamma}$ and ${\cal D}_{n}$ differ from their
refining counterparts ${\cal C}_{\eta(\gamma-\hat{\gamma})}$ and ${\cal
D}_{\eta\varepsilon}$, their dissimilarity merely affects the right hand side
of the inequality defining the sets, and are thus no more difficult to
construct. Just as in the case of (10), the task of obtaining separation
oracles for the refining problem (11) in the quantum regime reduces to
preparing many copies of Gibbs states. Likewise, these oracles can also be
implemented on a classical computer, given access to $\rho^{r}$.

The similarities between (10) and (11) become transparent when we demonstrate
that they are specific instances of the same problem. In particular, it is
easy to see that solving (10) corresponds to solving

find $\displaystyle\rho\in\\{X\in{\cal
S}^{n}_{+}:\operatorname{tr}(X)=1\\}\cap{\cal
C}_{\eta(\gamma-\hat{\gamma})}\cap{\cal D}_{\eta\varepsilon}$ (14) where
$\displaystyle{\cal
C}_{\eta(\hat{\gamma}-\gamma)}=\left\\{X:\operatorname{tr}\left(\frac{C}{\|C\|_{F}}Q\circ
X\right)\geq\eta(\gamma-\hat{\gamma})\right\\},$ $\displaystyle{\cal
D}_{\eta\varepsilon}=\left\\{X:\langle
i|X|i\rangle=\eta|\varepsilon_{i}|,\leavevmode\nobreak\ \forall
i\in[n]\right\\},$

with $\varepsilon_{i}=\frac{1}{n}$, $\eta=1$, $Q=ee^{\top}$, and
$\hat{\gamma}=0$. In view of this relationship, we can unify the oracle
construction for (10) and (11) as follows:

$\displaystyle O_{{\cal C}_{\eta(\gamma-\hat{\gamma})}}:\leavevmode\nobreak\ $
$\displaystyle\text{Compute an approximation}\leavevmode\nobreak\
\tilde{c}\leavevmode\nobreak\ \text{of}\leavevmode\nobreak\
\operatorname{tr}\left(Q\circ C\|C\|_{F}^{-1}\rho\right)\leavevmode\nobreak\
\text{up to additive error}\leavevmode\nobreak\
\frac{\epsilon}{4}.\leavevmode\nobreak\ $ $\displaystyle\text{Check
if}\leavevmode\nobreak\
\tilde{c}\geq\eta(\gamma-\hat{\gamma})+\frac{3\epsilon}{4}\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ \text{output}\leavevmode\nobreak\ P=-Q\circ
C\|C\|_{F}^{-1}\leavevmode\nobreak\ \text{if the inequality is violated}.$
$\displaystyle O_{{\cal D}_{\eta\varepsilon}}:\leavevmode\nobreak\ $
$\displaystyle\text{Compute an approximation}\leavevmode\nobreak\
\tilde{p}\in\mathbb{R}^{n}\leavevmode\nobreak\ \text{of}\leavevmode\nobreak\
p_{i}=\langle i|\rho|i\rangle\leavevmode\nobreak\
\text{satisfying}\leavevmode\nobreak\
\sum_{i\in[n]}|p_{i}-\tilde{p}_{i}|\leq\frac{\epsilon}{4}.$
$\displaystyle\text{Check if}\leavevmode\nobreak\
\sum_{i\in[n]}\left|\tilde{p}_{i}-\eta|\varepsilon_{i}|\right|\leq\frac{3\epsilon}{4}\leavevmode\nobreak\
\text{and output}\leavevmode\nobreak\
P=\sum_{i\in[n]}\left(\mathbb{I}\\{\tilde{p}_{i}>\eta|\varepsilon_{i}|\\}-\mathbb{I}\\{\tilde{p}_{i}<\eta|\varepsilon_{i}|\\}\right)\left|{i}\right\rangle\left\langle{i}\right|$
$\displaystyle\text{if the inequality is violated}.$

Again, the sets of observables for ${\cal C}_{\eta(\gamma-\hat{\gamma})}$ and
${\cal D}_{\eta\varepsilon}$ are given by

$\widetilde{{\cal C}}_{\eta(\gamma-\hat{\gamma})}=\\{-Q\circ
C\|C\|_{F}^{-1}\\},\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\
\widetilde{{\cal D}}_{\eta\varepsilon}=\\{D\in\mathbb{R}^{n\times n}:\|D\|\leq
1,D\leavevmode\nobreak\ \text{is diagonal}\\}.$

Although these observations are straightforward, they justify our use of
Algorithm 1 as a semidefinite optimization oracle that solves a convex
feasibility problem at hand in every iteration for different values of $Q$. In
particular, these facts, along with Lemmas 14 and 15 ensure that the
complexity results in Propositions 6 and 7 hold when applying Algorithm 1 to
solve (14).

######  Proposition 8.

Let $Q\circ\frac{C}{\|C\|_{F}}\in{\cal S}^{n}$ be stored in QRAM. Algorithm 1
solves (14) up to additive error $\epsilon$ using

$\widetilde{{\cal O}}_{\frac{n}{\epsilon}}\left(n^{1.5}\epsilon^{-5}\right)$

accesses to the QRAM.

###### Proof.

Given that $Q\circ\frac{C}{\|C\|_{F}}$ is stored in QRAM, Lemma 3(ii) asserts
that when constructing a block-encoding of $Q\circ\frac{C}{\|C\|_{F}}$, one
can set the subnormalization factor to be
$\alpha_{C}=\left\|Q\circ\frac{C}{\|C\|_{F}}\right\|_{F}$. In particular, one
can always choose $\alpha_{C}=1$, as it can be seen from the proof of Lemma 15
that the inequality

$\left\|Q\circ\frac{C}{\|C\|_{F}}\right\|_{F}\leq\left\|\frac{C}{\|C\|_{F}}\right\|_{F}=1$

always holds for any $Q$ defined according to equation (13). Collecting these
facts, one can construct a $(1,{\cal O}(\log(n)),\epsilon/(2n))$-block-
encoding of $Q\circ C\|C\|_{F}^{-1}$ in time $\widetilde{{\cal
O}}_{\frac{n}{\epsilon}}(1)$. Note that the quantity
$Q\circ\frac{C}{\|C\|_{F}}$ remains unchanged for the duration of Algorithm 1.
From here, the rest of the proof follows exactly that of Proposition 7 upon
replacing $\frac{C}{\|C\|_{F}}$, $O_{{\cal C}_{\gamma}}$ and $O_{{\cal
D}_{n}}$ with $Q\circ\frac{C}{\|C\|_{F}}$, $O_{{\cal
C}_{\eta(\gamma-\hat{\gamma})}}$ and $O_{{\cal D}_{\eta\varepsilon}}$,
respectively, in what remains. ∎

###  4.3 Iterative Refinement using Hamiltonian Updates

We are now in a position to provide our iterative refinement method for SDO
approximations of QUBOs presented in full detail in Algorithm 2.

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: A matrix $\hat{\rho}\in{\cal S}^{n}$ satisfying
$\max\left\\{\gamma-\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\hat{\rho}\right),\|\hat{\rho}-n^{-1}I\|_{\operatorname{tr}}\right\\}\leq\zeta$

Initialize: $\hat{\rho}\leftarrow\mathbf{0}^{n\times 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-\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\hat{\rho}\right),\|\varepsilon\|_{1}\right\\}
>\zeta$_ do

  1. 1.

Solve (14) to precision $\frac{\xi}{4}$ for $\rho^{(k)}$ using Algorithm 1
with oracles $O_{{\cal C}_{\eta(\gamma-\hat{\gamma})}}$ and $O_{{\cal
D}_{\eta\varepsilon}}$

  2. 2.

Update Solution:

$\hat{\rho}\leftarrow\hat{\rho}+\frac{1}{\eta^{(k)}}Q\circ\rho^{(k)},\quad\hat{\gamma}\leftarrow\hat{\gamma}+\frac{1}{\eta^{(k)}}\operatorname{tr}\left(\frac{C}{\|C\|_{F}}Q\circ\rho^{(k)}\right)$

  3. 3.

Compute element-wise deviations from the maximally mixed state:

$\varepsilon_{i}\leftarrow\hat{\rho}_{ii}-\frac{1}{n}\leavevmode\nobreak\
\text{for}\leavevmode\nobreak\ i\in[n]$

  4. 4.

Update refining problem parameters:

$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-\operatorname{tr}\left(\frac{C}{\|C\|}\hat{\rho}\right),\|\varepsilon\|_{1}\right\\}}$

  5. 5.

$k\leftarrow k+1$

end while

Algorithm 2 Iterative Refinement for SDO Approximations of QUBOs

The algorithm takes three parameters as input; (i) $\xi$, the fixed precision
used to test closeness to the sets ${\cal C}_{\eta(\gamma-\hat{\gamma})}$ and
${\cal D}_{\eta\varepsilon}$ in every iteration, (ii) $\zeta$, the precision
to which the final solution satisfies the functional constraints of (5), and
(iii) $\epsilon$, the additive error to which we seek to solve (3). In our
initialization steps we set the values of $Q$, $\varepsilon$ and $\eta$ such
that the first iteration corresponds to solving the feasibility problem (10).

In each iteration $k$, Algorithm 2 calls Algorithm 1 with separation oracles
$O_{{\cal C}_{\eta(\gamma-\hat{\gamma})}}$ and $O_{{\cal
D}_{\eta\varepsilon}}$ using fixed precision $\xi$ such that every call to
Algorithm 1 produces a $\xi$-precise classical solution $\rho^{(k)}$ to (14).
If $\hat{\rho}$ is indistinguishable up to precision $\zeta$ from the
maximally mixed state $n^{-1}I$ upon measurement in the computational basis,
and satisfies
$\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\hat{\rho}\right)\geq\gamma-\zeta$,
the algorithm terminates and reports $\hat{\rho}$. Otherwise, we construct the
refining problem associated with our current solution, and proceed to the next
iteration.

To define the parameters for the next refining problem, we first calculate the
deviation of the diagonal elements from $\frac{1}{n}$, and the violation with
respect to satisfying our objective value. Then, we define our scaling factor
to be the maximum over the $\ell_{1}$-norm of the diagonal deviations, and the
objective violation. We stress that $\xi$ is a (chosen) constant, and does not
change throughout the algorithm.

We now state a series of results in order to bound the iteration complexity of
Algorithm 2, and use our findings to improve the approximation guarantee given
in Proposition 5. We begin by proving that the iterates generated by Algorithm
2 satisfy the constraints in (5) with increasing accuracy.

######  Theorem 4.

Let $\hat{\rho}$ be the current overall solution, and let $\rho^{(k)}$ be a
solution to (14) obtained from running Algorithm 1 using fixed precision
$\xi\in(0,1)$ in iteration $k$ of Algorithm 2. Then, the following hold:

  * (a)

For $k\geq 0$, $\eta^{(k)}\geq\frac{1}{\xi^{k}}$.

  * (b)

For $k\geq 0$, $\rho=\hat{\rho}+\frac{1}{\eta^{(k)}}Q\circ\rho^{(k)}$
satisfies

$\max\left\\{\gamma-\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\rho\right),\|\rho-n^{-1}I\|_{\operatorname{tr}}\right\\}\leq\xi^{k+1}.$

###### Proof.

We begin by establishing that for $k\geq 0$, the current solution satisfies

$\max\left\\{\gamma-\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\rho\right),\|\rho-n^{-1}I\|_{\operatorname{tr}}\right\\}\leq\frac{\xi}{\eta^{(k)}}.$
(15)

First, observe that when $k=0$, we have $\varepsilon_{i}=\frac{1}{n}$ for
$i\in[n]$, $\hat{\rho}=\mathbf{0}^{n\times n}$, $\eta^{(0)}=1$ and
$Q=ee^{\top}$. Under these conditions, one can observe that if $\rho^{(0)}$
solves (14) to precision $\xi$, then

$\displaystyle\sum_{i=1}^{n}\left|\langle i|\rho|i\rangle-\frac{1}{n}\right|$
$\displaystyle=\sum_{i=1}^{n}\left|\langle
i|\hat{\rho}+\frac{1}{\eta^{(0)}}\rho^{(0)}|i\rangle-\frac{1}{n}\right|=\sum_{i=1}^{n}\left|\frac{1}{\eta^{(0)}}\rho^{(0)}-\varepsilon_{i}\right|=\frac{1}{\eta^{(0)}}\sum_{i=1}^{n}\left|\rho^{(0)}-\eta^{(0)}\varepsilon_{i}\right|\leq\frac{\xi}{\eta^{(0)}}.$

In other words, $\rho=\rho^{(0)}$ satisfies

$\|\rho-n^{-1}I\|_{\operatorname{tr}}\leq\frac{\xi}{\eta^{(0)}}.$

Next, by the definition of $O_{{\cal C}_{\eta(\gamma-\hat{\gamma})}}$ we have

$\displaystyle\operatorname{tr}\left(C\|C\|_{F}^{-1}\rho\right)$
$\displaystyle=\operatorname{tr}\left(C\|C\|_{F}^{-1}\left[\hat{\rho}+\frac{1}{\eta^{(0)}}Q\circ\rho^{(0)}\right]\right)=\operatorname{tr}\left(C\|C\|_{F}^{-1}\hat{\rho}\right)+\frac{1}{\eta^{(0)}}\operatorname{tr}\left(C\|C\|_{F}^{-1}\left(Q\circ\rho^{(0)}\right)\right)$
$\displaystyle=\underbrace{\operatorname{tr}\left(C\|C\|_{F}^{-1}\mathbf{0}^{n\times
n}\right)}_{=0}+\frac{1}{\eta^{(0)}}\operatorname{tr}\left(C\|C\|_{F}^{-1}\left(Q\circ\rho^{(0)}\right)\right)=\frac{1}{\eta^{(0)}}\operatorname{tr}\left(C\|C\|_{F}^{-1}\left(Q\circ\rho^{(0)}\right)\right)$
$\displaystyle=\frac{1}{\eta^{(0)}}\operatorname{tr}\left(C\|C\|_{F}^{-1}\rho^{(0)}\right)\geq\frac{1}{\eta^{(0)}}\left([\eta^{(0)}(\gamma-\hat{\gamma}^{(0)})]-\xi\right)=\gamma-\frac{\xi}{\eta^{(0)}},$

where we used the fact that $\hat{\gamma}^{(0)}=0$.

Next, letting $\hat{\rho}$ be our current solution and
$\hat{\gamma}=\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\hat{\rho}\right)$ be
the objective value attained at this solution. When $k\geq 1$, we have
$\varepsilon_{i}=\hat{\rho}_{ii}-\frac{1}{n}$ for $i\in[n]$ and
$Q=(ee^{\top}-I)+\textup{diag}(\textup{sign}(-\varepsilon))$. For this choice
of parameters, the general feasibility problem (14) reduces to the refining
problem (11) and the solution $\rho^{(k)}$ obtained via Algorithm 1 is
therefore a $\xi$-precise solution to (11). Accordingly, for $k\geq 1$,
setting $\rho=\hat{\rho}+\frac{1}{\eta^{(k)}}Q\circ\rho^{(k)}$ improves the
precision to which the maximally mixed state is approximated:

$\displaystyle\sum_{i=1}^{n}\left|\langle i|\rho|i\rangle-\frac{1}{n}\right|$
$\displaystyle=\sum_{i=1}^{n}\left|\langle
i|\hat{\rho}+\frac{1}{\eta^{(k)}}Q\circ\rho^{(k)}|i\rangle-\frac{1}{n}\right|$
$\displaystyle=\sum_{i=1}^{n}\left|\left(\hat{\rho}_{ii}+\frac{1}{\eta^{(k)}}\left(\textup{sign}(-\varepsilon_{i})\cdot\rho^{(k)}_{ii}\right)\right)-\frac{1}{n}\right|$
$\displaystyle=\sum_{i=1}^{n}\left|\left(\hat{\rho}_{ii}-\frac{1}{n}\right)+\frac{1}{\eta^{(k)}}\textup{sign}(-\varepsilon_{i})\rho^{(k)}_{ii}\right|$
$\displaystyle=\sum_{i=1}^{n}\left|\varepsilon_{i}+\frac{1}{\eta^{(k)}}\textup{sign}(-\varepsilon_{i})\rho^{(k)}_{ii}\right|=\frac{1}{\eta^{(k)}}\sum_{i=1}^{n}\left|\eta^{(k)}\varepsilon_{i}+\textup{sign}(-\varepsilon_{i})\rho^{(k)}_{ii}\right|\leq\frac{\xi}{\eta^{(k)}}.$

Consequently, we can conclude that at iteration $k\geq 1$, the trace distance
from our solution to the maximally mixed state is

$\|\rho-n^{-1}I\|_{\operatorname{tr}}\leq\frac{\xi}{\eta^{(k)}}.$ (16)

Next, letting
$\tilde{\gamma}^{(k)}=\operatorname{tr}{\left(\frac{C}{\|C\|_{F}}Q\circ\rho^{(k)}\right)}$,
one can observe

$\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\rho\right)=\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\left(\hat{\rho}+\frac{1}{\eta^{(k)}}Q\circ\rho^{(k)}\right)\right)=\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\hat{\rho}\right)+\frac{1}{\eta^{(k)}}\operatorname{tr}\left(\frac{C}{\|C\|_{F}}Q\circ\rho^{(k)}\right)=\hat{\gamma}+\frac{\tilde{\gamma}^{(k)}}{\eta^{(k)}}.$

For any $\rho^{(k)}$ which is $\xi$-close to the set ${\cal
C}_{\eta(\gamma-\hat{\gamma})}$ we must have

$\tilde{\gamma}^{(k)}\geq\eta^{(k)}(\gamma-\hat{\gamma})-\xi.$

It follows:

$\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\rho\right)=\hat{\gamma}+\frac{\tilde{\gamma}^{(k)}}{\eta^{(k)}}\geq\hat{\gamma}+\frac{1}{\eta^{(k)}}\left[\eta^{(k)}(\gamma-\hat{\gamma})-\xi\right]=\gamma-\frac{\xi}{\eta^{(k)}}.$
(17)

It therefore follows from (16) and (17) that
$\rho=\hat{\rho}+\frac{1}{\eta^{(k)}}Q\circ\rho^{(k)}$ satisfies inequality
(15) for all $k\geq 0$, and we can now use this fact to establish the lower
bound on $\eta^{(k)}$, which we prove by induction.

For $k=0$, we have $\eta^{(0)}=1$, for which $\eta^{(k)}\geq\frac{1}{\xi^{k}}$
trivially holds. By the induction hypothesis, it assumed that
$\eta^{(\ell)}\geq\frac{1}{\xi^{\ell}}$ is true for $\ell=1,\dots,k$. From
here, applying (15) yields

$\displaystyle\eta^{(k+1)}=\frac{1}{\max\left\\{\gamma-\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\rho\right),\|\rho-n^{-1}I\|_{\operatorname{tr}}\right\\}}$
$\displaystyle\geq\frac{1}{\frac{\xi}{\eta^{(k)}}}\geq\frac{1}{\xi^{k+1}},$

which completes the proof of (a).

Having demonstrated that (a) holds, to prove (b), we can simply combine
inequality (15) with the lower bound $\eta^{(k)}\geq\frac{1}{\xi^{k}}$, which
together imply

$\max\left\\{\gamma-\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\rho\right),\|\rho-n^{-1}I\|_{\operatorname{tr}}\right\\}\leq\frac{\xi}{\eta^{(k)}}\leq\xi^{k+1}.$

That is, (b) holds, and the proof is complete. ∎

The next result establishes polynomial convergence of Algorithm 2.

######  Corollary 5.

Let $0<\zeta\ll\xi<1$, and $\eta^{(0)}=1$. Then, Algorithm 2 terminates in at
most

$K={\cal O}\left(\log\left(\frac{1}{\zeta}\right)\right)$

iterations.

###### Proof.

The result follows from Theorem 4(b). ∎

Observe that in Theorem 4, we do not consider whether the updated solution
remains positive semidefinite. It turns out that, by the nature of our update
scheme, the minimum eigenvalue of the final solution obtained via Algorithm 2
will never fall significantly below zero. Moreover, we employ a rounding
procedure to modify the entries of the solution we obtain from Algorithm 2 so
that we arrive at an exactly feasible solution to (3). Before proceeding
further, we formalize this notion, and derive a lower bound on the smallest
eigenvalue of the final solution. We utilize Lemma 13 to bound the minimum
eigenvalue of the terms $\frac{1}{\eta^{(k)}}Q^{(k)}\circ\rho^{(k)}$ that are
used to update the overall solution in each iteration, from which a lower
bound minimum eigenvalue bound for the final solution obtained by Algorithm 2
readily follows upon applying Lemma 12.

######  Proposition 9.

Let $\rho^{(k)}$ be a solution to (14) obtained from running Algorithm 1 using
fixed precision $\xi\in(0,1)$ in iteration $k$ of Algorithm 2. Then:

  * (a)

For $k\geq 1$,
$\frac{1}{\eta^{(k)}}Q^{(k)}\circ\rho^{(k)}\succeq-2\cdot\xi^{k}I$.

  * (b)

Suppose Algorithm 2 is run with final precision $\zeta$, and terminates after
$K$ iterations. Then, the solution $\rho$ output by Algorithm 2 satisfies

$\lambda_{\min}(\rho)\geq\lambda_{\min}\left(\rho^{(0)}\right)-2\cdot\sum_{k=1}^{K}\xi^{k}.$

###### Proof.

We begin with a proof of (a). In what follows, we assume without loss of
generality that $Q^{(k)}$ has at least one negative eigenvalue (otherwise,
$Q^{(k)}\circ\rho^{(k)}\succeq 0$ trivially holds). Hence, a combined
application of Corollary 1 and Lemma 13 yields

$\lambda_{\min}\left(Q^{(k)}\circ\rho^{(k)}\right)\geq\lambda_{\max}\left(\rho^{(k)}\right)\cdot\lambda_{\min}\left(Q^{(k)}\right)\geq-2.$

Therefore, recalling that by Theorem 4(a) we have
$\eta^{(k)}\geq\frac{1}{\xi^{k}}$, it follows

$\displaystyle\lambda_{\min}\left(\frac{1}{\eta^{(k)}}Q\circ\rho^{(k)}\right)=\frac{1}{\eta^{(k)}}\lambda_{\min}\left(Q\circ\rho^{(k)}\right)\geq\frac{-2}{\eta^{(k)}}\geq-2\cdot\xi^{k},$

which completes the proof of (a).

Noting that the final solution can be expressed as

$\rho=\sum_{k=0}^{K}\frac{1}{\eta^{(k)}}Q^{(k)}\circ\rho^{(k)}=\rho^{(0)}+\sum_{k=1}^{K}\frac{1}{\eta^{(k)}}Q^{(k)}\circ\rho^{(k)},$

the result in (b) follows from (a) and repeated application of Lemma 12. ∎

The next corollary computes bounds for the geometric series that appears in
Proposition 9(b) for different values of the fixed precision parameter $\xi$,
by computing the value of the series to $\infty$.

######  Corollary 6.

Let $\xi=\frac{\tilde{\xi}}{4}$. Then, for any positive integer $K$, the
following hold.

  * (a)

If $\tilde{\xi}=10^{-2}$, then $-2\cdot\sum_{k=1}^{K}\xi^{k}\geq-0.005$.

  * (b)

If $\tilde{\xi}=10^{-4}$, then $-2\cdot\sum_{k=1}^{K}\xi^{k}\geq-0.00005$.

The goal of Corollary 6 is simply to show that with fixed precision we obtain
matrices with eigenvalues that are, in the worst case, slightly negative. A
shift of the spectrum suffices to restore positive semidefiniteness, and it
does not change the constraint violation or the objective function value by a
large amount, as we show next.

######  Proposition 10.

Suppose Algorithm 2 is run with final precision $\zeta$, and terminates after
$K$ iterations. Let $\rho$ be the solution output by Algorithm 2, and let
$\xi$ be the fixed precision used in every iteration. Then, letting

$\delta=2\cdot\sum_{k=1}^{K}\xi^{k},$ (18)

it follows that

$\tilde{\rho}=\frac{1}{1+n\delta}\left(\rho+\delta I\right)$

is a positive semidefinite matrix at trace distance $\zeta$ from $\rho$.
Moreover, $\tilde{\rho}$ satisfies:

$\max\left\\{\gamma-\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\tilde{\rho}\right),\left\|\tilde{\rho}-n^{-1}I\right\|_{\operatorname{tr}}\right\\}\leq
2\zeta.$

In other words, $\tilde{\rho}$ is a $2\zeta$-precise solution to (5).

###### Proof.

First, observe that $\tilde{\rho}\succeq 0$ by the definition of
$\tilde{\rho}$; applying Proposition 9(b), we have
$\lambda_{\min}(\rho)\geq-\delta$, which implies that $\rho+\delta I\succeq
0$. From the definition of $\tilde{\rho}$, we also have:

$\displaystyle\left\|\rho-\tilde{\rho}\right\|_{\operatorname{tr}}=\left\|\rho-\left[\frac{1}{1+n\delta}\left(\rho+\delta
I\right)\right]\right\|_{\operatorname{tr}}$
$\displaystyle=\left\|\frac{1+n\delta-1}{1+n\delta}\rho-\frac{\delta}{1+n\delta}I\right\|_{\operatorname{tr}}$
$\displaystyle=\frac{\delta}{1+n\delta}\left\|n\rho-I\right\|_{\operatorname{tr}}$
$\displaystyle=\frac{n\delta}{1+n\delta}\left\|\rho-n^{-1}I\right\|_{\operatorname{tr}}$
$\displaystyle\leq\frac{n\delta}{1+n\delta}\zeta$ $\displaystyle<\zeta,$

where the second to last inequality follows from the fact that $\rho$ is
obtained from running Algorithm 2 with $\zeta$ as the final precision
parameter.

Next, we leverage our bound on the trace distance from $\tilde{\rho}$ to
$\rho$ to establish that $\tilde{\rho}$ is indeed an accurate solution to the
renormalized SDO problem (5). First, note that

$\displaystyle\left\|\tilde{\rho}-n^{-1}I\right\|_{\operatorname{tr}}$
$\displaystyle=\left\|\tilde{\rho}-n^{-1}I+(\rho-\rho)\right\|_{\operatorname{tr}}$
$\displaystyle=\left\|\left(\tilde{\rho}-\rho\right)+\left(\rho-n^{-1}I\right)\right\|_{\operatorname{tr}}\leq\left\|\tilde{\rho}-\rho\right\|_{\operatorname{tr}}+\left\|\rho-n^{-1}I\right\|_{\operatorname{tr}}\leq
2\zeta.$

Further, applying a matrix Hölder inequality, one can observe:

$\displaystyle\left|\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\rho\right)-\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\tilde{\rho}\right)\right|\leq\left\|\frac{C}{\|C\|_{F}}\right\|\|\rho-\tilde{\rho}\|_{\operatorname{tr}}\leq\|\rho-\tilde{\rho}\|_{\operatorname{tr}}<\zeta,$

from which we can conclude

$\displaystyle\gamma-\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\tilde{\rho}\right)$
$\displaystyle=\gamma-\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\tilde{\rho}\right)+\left[\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\rho\right)-\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\rho\right)\right]$
$\displaystyle=\left[\gamma-\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\rho\right)\right]+\left[\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\rho\right)-\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\tilde{\rho}\right)\right]\leq
2\zeta,$

as $\gamma-\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\rho\right)\leq\zeta$. ∎

It is important at this point for us to remark that fixing $\xi\in(0,1)$ does
not limit us with respect to how accurately we can solve (3). We can always
make the final precision parameter arbitrarily small using only
$\widetilde{{\cal O}}_{\frac{1}{\zeta}}(1)$ iterations, as the overall running
time depends only poly-logarithmically on $\zeta^{-1}$. Accordingly, we take
advantage of this fact and revisit the approximation guarantee provided in
Proposition 5.

######  Proposition 11.

Let $\rho$ be a $\zeta$-accurate solution to the renormalized and relaxed SDO
problem (5) with input matrix $C$ and
$\zeta=\left(\frac{\epsilon}{n\|C\|_{F}}\right)^{4}$. Let
$\gamma_{\zeta}=\operatorname{tr}{(C\rho)}$ be the value attained by $\rho$.
Then, there is a quantum state $\rho^{*}$ at trace distance ${\cal
O}\left(\frac{\epsilon}{n\|C\|_{F}}\right)$ of $\rho$ such that $n\rho^{*}$ is
a feasible point of SDO problem (3). In particular

$\left|\gamma_{\zeta}n\|C\|_{F}-\operatorname{tr}\left(n\rho^{*}C\right)\right|={\cal
O}\left(\epsilon\right).$

Moreover, it is possible to construct $\rho^{*}$ in time ${\cal O}(n^{2})$
given the entries of $\rho$.

###### Proof.

The proof almost exactly follows the proof of Proposition 3.1 in [10],
regardless, we present the adjusted proof for completeness. Our aim is to show
that a $\zeta$-precise solution $\rho$ to (5) obtained using Algorithm 2 can
be used to construct $\rho^{*}$ such that $n\rho^{*}$ is an exactly feasible
solution to (3).

We begin shifting $\rho$ in order to ensure that our solution is positive
semidefinite. In particular, choosing $\delta$ according to (18), we set

$\tilde{\rho}=\frac{1}{1+n\delta}\left(\rho+\delta I\right).$

It then follows from Proposition 10 that $\tilde{\rho}$ satisfies

$\tilde{\rho}\succeq
0,\quad\left\|\rho-\tilde{\rho}\right\|_{\operatorname{tr}}\leq\zeta,\quad\max\left\\{\gamma-\operatorname{tr}\left(\frac{C}{\|C\|_{F}}\tilde{\rho}\right),\left\|\tilde{\rho}-n^{-1}I\right\|_{\operatorname{tr}}\right\\}\leq
2\zeta.$ (19)

Next, we examine the diagonal elements of $\tilde{\rho}$ and check whether
modifications need to be made to ensure that our solution is an exactly
feasible point to the renormalized SDO problem (5). Namely, if $|\langle
i|\tilde{\rho}|i\rangle-\frac{1}{n}|>\frac{\sqrt{2\zeta}}{n}$ for $i\in[n]$,
we replace $\tilde{\rho}_{ii}$ with $\frac{1}{n}$ and set all elements in the
$i$-th row and the $i$-th column to 0, and denote the resulting matrix by
$\rho^{\prime}$. From here we introduce another matrix $W$ which we obtain by
replacing each diagonal entry of $\rho^{\prime}$ with $\frac{1}{n}$. In
general we may not have $W\succeq 0$, so the authors in [10] suggest using the
convex combination:

$\rho^{*}=\frac{1}{1+\sqrt{2\zeta}}\left(W+\frac{\sqrt{2\zeta}}{n}I\right).$

Then, $\rho^{*}\succeq 0$ and by construction $\langle
i|\rho^{*}|i\rangle=\frac{1}{n}$ for all $i\in[n]$. Hence, $\rho^{*}$ is a
feasible solution to the renormalized SDO problem (5).

What remains is to show that the above reformulations yield the desired
approximation. Denote by ${\cal B}=\\{i:|n\langle
i|\tilde{\rho}|i\rangle-1|>\sqrt{2\zeta}\\}\subset[n]$ the set of diagonal
entries that deviate substantially from $\frac{1}{n}$. Without loss of
generality, it suffices to assume that such elements are found in the first
$|{\cal B}|$ rows of $\tilde{\rho}$, in which case

$\displaystyle\|\rho^{\prime}-\tilde{\rho}\|_{\operatorname{tr}}=\left\|\begin{pmatrix}n^{-1}I_{{\cal
B}}&0\\\
0&\tilde{\rho}_{22}\end{pmatrix}-\begin{pmatrix}\tilde{\rho}_{11}&\tilde{\rho}_{12}\\\
\tilde{\rho}_{21}&\tilde{\rho}_{22}\end{pmatrix}\right\|_{\operatorname{tr}}$
$\displaystyle=\left\|\begin{pmatrix}n^{-1}I_{{\cal
B}}-\tilde{\rho}_{11}&-\tilde{\rho}_{12}\\\
-\tilde{\rho}_{21}&0\end{pmatrix}\right\|_{\operatorname{tr}}$
$\displaystyle\leq\|\tilde{\rho}_{11}\|_{\operatorname{tr}}+2\|\tilde{\rho}_{12}\|_{\operatorname{tr}}+\|n^{-1}I_{{\cal
B}}\|_{\operatorname{tr}}.$ (20)

Since $\tilde{\rho}$ is a $2\zeta$-precise solution to (5), $\tilde{\rho}$
obeys

$\sum_{i=1}^{n}\left|\langle i|\tilde{\rho}|i\rangle-\frac{1}{n}\right|\leq
2\zeta.$

Therefore, we must have

$|{\cal B}|\frac{\sqrt{2\zeta}}{n}\leq 2\zeta,$

which equates to $|{\cal B}|\leq n\sqrt{2\zeta}$. Now, by the definition of
${\cal B}$, it follows

$\|\tilde{\rho}_{22}\|_{\operatorname{tr}}\geq(n-|{\cal
B}|)\frac{1-\sqrt{2\zeta}}{n}\geq(n-n\sqrt{2\zeta})\frac{1-\sqrt{2\zeta}}{n}=(1-\sqrt{2\zeta})^{2}.$

Following [10], we invoke a result from [40], which states

$\left\|\begin{bmatrix}\|\tilde{\rho}_{11}\|_{\operatorname{tr}}&\|\tilde{\rho}_{12}\|_{\operatorname{tr}}\\\
\|\tilde{\rho}^{\top}_{12}\|_{\operatorname{tr}}&\|\tilde{\rho}_{22}\|_{\operatorname{tr}}\end{bmatrix}\right\|\leq\left\|\begin{bmatrix}\tilde{\rho}_{11}&\tilde{\rho}_{12}\\\
\tilde{\rho}^{\top}_{12}&\tilde{\rho}_{22}\end{bmatrix}\right\|_{\operatorname{tr}}=\|\tilde{\rho}\|_{\operatorname{tr}}=\operatorname{tr}{(\tilde{\rho})}=1.$

Using the fact that $\|\cdot\|_{\operatorname{tr}}\geq\|\cdot\|_{2}$, where
$\|\cdot\|_{2}$ is the Frobenius, or Schatten-2 norm, the above implies

$\|\tilde{\rho}_{11}\|^{2}_{\operatorname{tr}}+2\|\tilde{\rho}_{12}\|^{2}_{\operatorname{tr}}+\|\tilde{\rho}_{22}\|^{2}_{\operatorname{tr}}\leq
1.$

As $\|\tilde{\rho}_{22}\|_{\operatorname{tr}}\geq(1-\sqrt{2\zeta})^{2}$, it
can be seen trivially that
$\|\tilde{\rho}_{22}\|^{2}_{\operatorname{tr}}\geq(1-\sqrt{2\zeta})^{4}$, and
thus

$\|\tilde{\rho}_{11}\|^{2}_{\operatorname{tr}}+2\|\tilde{\rho}_{12}\|^{2}_{\operatorname{tr}}\leq
1-(1-\sqrt{2\zeta})^{4}={\cal O}(\sqrt{\zeta}).$

Consequently
$\|\tilde{\rho}_{11}\|_{\operatorname{tr}}+2\|\tilde{\rho}_{12}\|_{\operatorname{tr}}={\cal
O}\left(\zeta^{\frac{1}{4}}\right)$, and plugging this into equation (20)
asserts

$\|\rho^{\prime}-\tilde{\rho}\|_{\operatorname{tr}}={\cal
O}\left(\zeta^{\frac{1}{4}}\right).$ (21)

Let $R$ be a diagonal matrix whose elements are
$R_{ii}\in\left[-\frac{\sqrt{2\zeta}}{n},\frac{\sqrt{2\zeta}}{n}\right]$ for
$i\in[n]$, such that

$W=\rho^{\prime}+R,$

and note that $R+\sqrt{2\zeta}n^{-1}I\succeq 0$. Upon normalizing the trace,
one can observe

$\rho^{*}=\frac{1}{1+\sqrt{2\zeta}}\left(\rho^{\prime}+R+\sqrt{2\zeta}n^{-1}I\right)\succeq
0,$

with $\rho^{*}_{ii}=\frac{1}{n}n$ for all $i\in[n]$. Thus, $n\rho^{*}$ is a
feasible solution to the SDO problem (3). Further, by a triangle inequality we
have

$\|\rho^{\prime}-\rho^{*}\|_{\operatorname{tr}}=\frac{1}{1+\sqrt{2\zeta}}\left\|\sqrt{2\zeta}\rho^{\prime}+R+\sqrt{2\zeta}n^{-1}I\right\|_{\operatorname{tr}}={\cal
O}(\sqrt{\zeta}).$ (22)

Combining equations (21) and (22) and noting
$\zeta=\left(\frac{\epsilon}{n\|C\|_{F}}\right)^{4}$, applying another
triangle inequality yields

$\|\tilde{\rho}-\rho^{*}\|_{\operatorname{tr}}={\cal
O}\left(\zeta^{\frac{1}{4}}\right)={\cal
O}\left(\left[\left(\frac{\epsilon}{n\|C\|_{F}}\right)^{4}\right]^{\frac{1}{4}}\right)={\cal
O}\left(\frac{\epsilon}{n\|C\|_{F}}\right).$

Then, the result follows from a matrix Hölder inequality:

$\displaystyle\left|\operatorname{tr}\left(nC\rho\right)-\operatorname{tr}\left(nC\rho^{*}\right)\right|\leq
n\|C\|\|\rho-\rho^{*}\|_{\operatorname{tr}}$ $\displaystyle\leq
n\|C\|_{F}\left(\|\rho-\tilde{\rho}\|_{\operatorname{tr}}+\|\tilde{\rho}-\rho^{*}\|_{\operatorname{tr}}\right)$
$\displaystyle={\cal
O}\left(n\|C\|_{F}\left(\zeta+\zeta^{\frac{1}{4}}\right)\right)$
$\displaystyle={\cal
O}\left(n\|C\|_{F}\left[\left(\frac{\epsilon}{n\|C\|_{F}}\right)^{4}+\frac{\epsilon}{n\|C\|_{F}}\right]\right)$
$\displaystyle={\cal
O}\left(\frac{\epsilon^{4}}{\left(n\|C\|_{F}\right)^{3}}+\epsilon\right)$
$\displaystyle={\cal O}\left(\epsilon\right).$

∎

##  5 Complexity

We now analyze the worst case overall running time of our Iterative Refinement
Method given in Algorithm 2 in both the classical and quantum settings.

###  5.1 Classical running time

As we saw in Section 3, the complexity of using Algorithm 1 to solve the SDO
problem (3) scales poorly in the inverse precision, with the classical
algorithm exhibiting an ${\cal O}(\epsilon^{-12})$ dependence. In both the
classical and quantum cases, our iterative refinement scheme reconciles the
poor scaling in $\epsilon$ because it possesses the following two properties.
First, we can obtain an arbitrarily precise solution to (5) in at most
$\widetilde{{\cal O}}_{\frac{1}{\zeta}}(1)$ iterations. Second, it suffices to
treat $\xi$ as fixed for the oracle calls that occur in each iteration, as the
precision of the final solution is a byproduct of how we use these solution of
the refining problems to produce a solution to (5).

The next result formalizes the above argument, and establishes the complexity
of Algorithm 2 for the classical case.

######  Theorem 5.

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

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

The output of the algorithm is a classical description of a matrix
$\hat{\rho}\in{\cal S}^{n}$ such that

$\tilde{\rho}=\frac{1}{1+n\delta}\left(\hat{\rho}+\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, and from here we must compute $\|C\|_{F}$, which requires ${\cal O}(ns)$
arithmetic operations. In every iteration of Algorithm 2, we make a call to
our subroutine in Algorithm 1, before updating the solution and preparing the
next refining problem. Updating the solution involves matrix addition between
two $n\times n$ matrices and requires ${\cal O}(n^{2})$ arithmetic operations,
whereas updating $Q$ and $\varepsilon$ for the next refining problem can be
accomplished using ${\cal O}(n)$ arithmetic operations, as only the diagonal
entries of $Q$ need to be stored and maintained.

In view of Proposition 6, the dominant operation at each iteration is the use
of Algorithm 1 to solve the SDO problem at hand. By Proposition 6, Algorithm 1
can be used to solve (14) to additive error $\xi$ in time

${\cal T}_{HU}^{\text{classical}}={\cal
O}\left(\min\\{n^{2}s,n^{\omega}\\}\log^{2}(n)\xi^{-3}\right).$

If every call to Algorithm 1 is made using precision $\xi$, then by Corollary
5, Algorithm 2 converges in at most ${\cal O}\left(\log(\zeta^{-1})\right)$
iterations, and we can thus express the overall running time of Algorithm 2 as

${\cal
O}\left(\left(\min\\{n^{2}s,n^{\omega}\\}\log^{2}(n)\xi^{-3}\right)\log(\zeta^{-1})\right).$

In the context of Algorithm 2, it suffices to carry out each of the calls to
the SDO subroutine (i.e., calls to Algorithm 1) using fixed precision $\xi$ to
obtain a $2\zeta$-precise solution to (5) (see, e.g., Proposition 10). The
above complexity thus reduces to

${\cal O}\left(\min\\{n^{2}s,n^{\omega}\\}\log^{2}(n)\log(\zeta^{-1})\right).$

For our choice of $\zeta=\left(\frac{\epsilon}{n\|C\|_{F}}\right)^{4}$, one
can observe

${\cal
O}\left(\min\\{n^{2}s,n^{\omega}\\}\log^{2}(n)\log(\zeta^{-1})\right)={\cal
O}\left(\min\\{n^{2}s,n^{\omega}\\}\cdot\textup{polylog}\left(n,\|C\|_{F},\frac{1}{\epsilon}\right)\right).$

Proposition 11 certifies that the above running time suffices to obtain a
$\rho$ from which we can construct $\rho^{*}$ in time ${\cal O}(n^{2})$, such
that $n\rho^{*}$ is a feasible point of the SDO problem (3) satisfying

$\left|\gamma_{\zeta}n\|C\|_{F}-\operatorname{tr}\left(n\rho^{*}C\right)\right|={\cal
O}(\epsilon),$

and the proof is complete. ∎

###  5.2 Quantum running time

Just as in the classical case, we show that a quantum implementation of
Algorithm 2 mitigates the poor scaling in the running time with respect to the
inverse precision.

Our quantum implementation of Algorithm 2 is provided in Algorithm 3. The
relevant error parameters are the same as those appearing in Algorithm 2: (i)
$\xi$, the fixed precision used to test closeness to the sets ${\cal
C}_{\eta(\gamma-\hat{\gamma})}$ and ${\cal D}_{\eta\varepsilon}$ in every
iteration, (ii) $\zeta$, the precision to which the final solution satisfies
the functional constraints in (5), and (iii) $\epsilon$, the additive error to
which we seek to solve (3). In our initialization steps we set the values of
$Q$, $\varepsilon$ and $\eta$ such that the first iteration corresponds to
solving the feasibility problem (10). We also create a vector
$p=\mathbf{0}^{n\times 1}$ that will be used to maintain a classical
description of the diagonal elements of our solution over the course of the
algorithm.

At every iteration $k$, a call is made to Algorithm 1 with separation oracles
$O_{{\cal C}_{\eta(\gamma-\hat{\gamma})}}$ and $O_{{\cal
D}_{\eta\varepsilon}}$ to solve (14) using fixed precision $\xi$. If the
oracles accept the candidate state, then Algorithm 1 returns a real-valued
vector $y^{(k)}\in\mathbb{R}^{2}$ along with a diagonal matrix $D^{(k)}$ such
that the Hamiltonian associated with the Gibbs state that solves the refining
problem is

$H^{(k)}=y^{(k)}_{1}Q^{(k)}\circ\frac{C}{\|C\|_{F}}+y^{(k)}_{2}D^{(k)},$

with $\|y^{(k)}\|_{1}\leq 4\log(n)\xi^{-1}$ and $\|D^{(k)}\|\leq 1$ for every
$k\geq 0$. This allows us to efficiently describe the solution to each
refining problem, and once the algorithm has terminated, it facilitates an
efficient way to describe the final solution as well.555Requiring an explicit
classical description of the solution would in fact lead to a worse running
time overall when compared to the classical implementation we studied in
Section 5.1. First, observe that the matrices $Q^{(k)}$ and $D^{(k)}$ can be
completely described by their diagonal elements; letting
$q^{(k)}\in\mathbb{R}^{n}$ and $d^{(k)}\in\mathbb{R}^{n}$ be the vectors that
store the diagonal elements of $Q^{(k)}$ and $D^{(k)}$, respectively, we have

$\displaystyle Q^{(k)}$
$\displaystyle=(ee^{\top}-I)+\text{diag}\left(q^{(k)}\right),$ $\displaystyle
D^{(k)}$ $\displaystyle=\text{diag}\left(d^{(k)}\right).$

Therefore, we store the solution to the refining problem at iteration $k$ as
the tuple

$(\eta^{(k)},y^{(k)},q^{(k)},d^{(k)}),$

and the final solution to (5) is defined as