# Quantum algorithms for Second-Order Cone Programming and Support Vector
Machines

Iordanis Kerenidis  QCWare, Palo Alto, California  Université de Paris, CNRS,
IRIF, F-75006, Paris, France<EMAIL_ADDRESS>Anupam Prakash  QCWare, Palo
Alto, California  Université de Paris, CNRS, IRIF, F-75006, Paris, France
<EMAIL_ADDRESS>Dániel Szilágyi  Université de Paris, CNRS, IRIF,
F-75006, Paris, France<EMAIL_ADDRESS>
###### Abstract

We present a quantum interior-point method (IPM) for second-order cone
programming (SOCP) that runs in time
$\widetilde{O}\left(n\sqrt{r}\frac{\zeta\kappa}{\delta^{2}}\log\left(1/\epsilon\right)\right)$
where $r$ is the rank and $n$ the dimension of the SOCP, $\delta$ bounds the
distance of intermediate solutions from the cone boundary, $\zeta$ is a
parameter upper bounded by $\sqrt{n}$, and $\kappa$ is an upper bound on the
condition number of matrices arising in the classical IPM for SOCP. The
algorithm takes as its input a suitable quantum description of an arbitrary
SOCP and outputs a classical description of a $\delta$-approximate
$\epsilon$-optimal solution of the given problem.

Furthermore, we perform numerical simulations to determine the values of the
aforementioned parameters when solving the SOCP up to a fixed precision
$\epsilon$. We present experimental evidence that in this case our quantum
algorithm exhibits a polynomial speedup over the best classical algorithms for
solving general SOCPs that run in time $O(n^{\omega+0.5})$ (here, $\omega$ is
the matrix multiplication exponent, with a value of roughly $2.37$ in theory,
and up to $3$ in practice). For the case of random SVM (support vector
machine) instances of size $O(n)$, the quantum algorithm scales as $O(n^{k})$,
where the exponent $k$ is estimated to be $2.59$ using a least-squares power
law. On the same family random instances, the estimated scaling exponent for
an external SOCP solver is $3.31$ while that for a state-of-the-art SVM solver
is $3.11$.

##  1 Introduction

It is well known that many interesting and relevant optimization problems in
the domain of Machine Learning can be expressed in the framework of convex
optimization [6, 10]. The landmark result in this area was the discovery of
interior-point methods (IPM) by [23], and their subsequent generalization to
all “self-scaled” (i.e. symmetric) cones by [29, 30]. Very recently, [13] have
shown that it is possible to solve linear programs (LP) in
$\widetilde{O}(n^{\omega})$, the time it takes to multiply two matrices (as
long as $\omega\geq 2+1/6$, which is currently the case). This result has been
further extended in [27] to a slightly more general class of cones, however,
their techniques did not yield improved complexities for second-order (SOCP)
and semidefinite programming (SDP). An efficient algorithm for SOCP would also
yield efficient algorithms for many interesting problems, such as (standard
and quadratically-constrained) convex quadratic programming, portfolio
optimization, and many others [1].

Starting with the landmark results of [20, 35], and, more recently, [21], it
has been demonstrated that quantum computers offer significant (sometimes even
exponential) asymptotic speedups for a number of important problems. More
recently, there has been substantial work in the area of convex optimization.
Quantum speedups for gradient descent were investigated by [17], whereas [9,
7, 41, 40] presented quantum algorithms for SDP based on the the
multiplicative weights framework of [4]. However, it has been difficult to
establish asymptotic speedups for this family of quantum SDP solvers as their
the running time depends on problem-specific parameters, including a 5th-power
dependence on the width of the SDP. Interestingly, the recent result of [8]
suggests that such a speedup might be obtained when applying an SDP algorithm
of this type to some low-precision instances of quadratic binary optimization.

In an orthogonal approach, [26] proposed a quantum algorithm for LPs and SDPs
by quantizing a variant of the classical interior point method and using the
state of the art quantum linear algebra tools [11, 18] – in particular, the
matrix multiplication and inversion algorithms whose running time is sub-
linear in the input size. However, the complexity of this algorithm depends on
the condition number of $O(n^{2})$-sized matrices that is difficult to bound
theoretically. It therefore remains an open question to find an end-to-end
optimization problems for which quantum SDP solvers achieve an asymptotic
speedup over state of the art classical algorithms. In this paper, we propose
support vector machines (SVM) as a candidate for such an end-to-end quantum
speedup using a quantum interior point method based algorithm.

###  1.1 Our results and techniques

In this section, we provide a high level sketch of our results and the
techniques used for the quantum interior point method for SOCPs, we begin by
discussing the differences between classical and quantum interior point
methods.

A classical interior point method solves an optimization problem over
symmetric cones by starting with a feasible solution and iteratively finding
solutions with a smaller duality gap while maintaining feasibility. A single
iterative step consists of solving a system of linear equations called the
Newton linear system and updating the current iterate using the solutions of
the Newton system. The analysis of the classical IPM shows that in each
iteration, the updated solutions remain feasible and the duality gap is
decreased by a factor of $(1-\alpha/\sqrt{n})$ where $n$ is the dimension of
the optimization problem and $\alpha>0$ is a constant. The algorithm therefore
converges to a feasible solution with duality gap $\epsilon$ in
$O(\sqrt{n}\log(1/\epsilon))$ iterations.

A quantum interior point method [26] uses a quantum linear system solver
instead of classical one in each iteration of the IPM. However, there is an
important difference between classical and quantum linear algebra procedures
for solving the linear system $A\bm{x}=\bm{b}$. Unlike classical linear system
solvers which return an exact description of $\bm{x}$, quantum tomography
procedures can return an $\epsilon$-accurate solution $\widetilde{\bm{x}}$
such that $\lVert\overline{\bm{x}}-\bm{x}\rVert\leq\epsilon\lVert\bm{x}\rVert$
with $O(n/\epsilon^{2})$ runs of the quantum linear system solver.
Additionally, these linear system solvers require $A$ and $\bm{b}$ to be given
as block-encodings [11], so this input model is used by our algorithm as well.
One of the main technical challenges in developing a quantum interior point
method (QIPM) is to establish convergence of the classical IPM which uses
$\epsilon$-approximate solutions of the Newton linear system (in the
$\ell_{2}$ norm) instead of the exact solutions in the classical IPM.

The quantum interior point method for second-order cone programs requires
additional ideas going beyond those used for the quantum interior point
methods for SDP [26]. Second order cone programs are optimization problems
over the product of second-order or Lorentz cones (see section 2.1 for
definitions), interior point methods for SOCP can be described using the
Euclidean Jordan algebra framework [28]. The Euclidean Jordan algebra
framework provides analogs of concepts like eigenvalues, spectral and
Frobenius norms for matrices and positive semidefinite constraints for the
case of SOCPs. Using these conceptual ideas from the Euclidean Jordan algebra
framework [28] and the analysis of the approximate SDP interior point method
[26] we provide an approximate IPM for SOCP that converges in
$O(\sqrt{n}\log(1/\epsilon))$ iterations. Approximate IPMs for SOCP have not
been previously investigated in the classical or the quantum optimization
literature, this analysis is one of the main technical contributions of this
paper.

From an algorithmic perspective, SOCPs are much closer to LPs (Linear
Programs) than to SDPs, since for cones of dimension $n$, the Newton linear
systems arising in LP and SOCP IPMs are of size $O(n)$, whereas in the SDP
case they are of size $O(n^{2})$. Namely, a second-order conic constraint of
dimension $n$ can be expressed as a single PSD constraint on a (sparse)
$n\times n$ matrix [1] – this allows us to embed an SOCP with $n$ variables
and $m$ constraints in an SDP with an $n\times n$ matrix and $m$ constraints.
The cost of solving that SDP would have a worse dependence on the error [40]
or the input size [26]. On the other hand, the quantum representations (block
encodings) of the Newton linear systems for SOCP are also much simpler to
construct than those for SDP. The smaller size of the SOCP linear system also
makes it feasible to empirically estimate the condition number for these
linear systems in a reasonable amount of time allowing us to carry out
extensive numerical experiments to validate the running time of the quantum
algorithm.

The theoretical analysis of the quantum algorithm for SOCP shows that its
worst-case running time is

$\widetilde{O}\left(n\sqrt{r}\frac{\zeta\kappa}{\delta^{2}}\log\left(\frac{1}{\epsilon}\right)\right),$
(1)

where $r$ is the rank and $n$ the dimension of the SOCP, $\delta$ bounds the
distance of intermediate solutions from the cone boundary, $\zeta$ is a
parameter bounded by $\sqrt{n}$, $\kappa$ is an upper bound on the condition
number of matrices arising in the interior-point method for SOCPs, and
$\epsilon$ is the target duality gap. The running time of the algorithm
depends on problem dependent parameters like $\kappa$ and $\delta$ that are
difficult to bound in terms of the problem dimension $n$ – this is also the
case with previous quantum SDP solvers [26, 40] and makes it important to
validate the quantum optimization speedups empirically. Interestingly, since
we require a classical solution of the Newton system, the linear system solver
could also be replaced by a classical iterative solver [34] which would yield
a complexity of $O(n^{2}\sqrt{r}\kappa\log(n/\epsilon))$.

Let us make a remark about the complexity: as it is the case with all
approximation algorithms, (1) depends on the inverse of the target duality gap
$\epsilon$, thus the running time of our algorithm grows to infinity as
$\epsilon$ approaches zero, as in the case of classical IPM. Our running time
also depends on $\kappa$, which in turn is empirically observed to grow
inversely with the duality gap (in particular as $O(1/\epsilon)$) which again
makes the running time go to infinity as $\epsilon$ approaches zero. The
quantum IPM is a low precision method, unlike the classical IPM, and it can
offer speedups for settings where the desired precision $\epsilon$ is moderate
or low. Thus although at first glance it seems that the $\epsilon$-dependence
in (1) is logarithmic, experimental evidence suggests that the factor
$\frac{\zeta\kappa}{\delta^{2}}$ depends polynomially on $1/\epsilon$.

Support Vector Machines (SVM) is an important application in machine learning,
where even a modest value of $\epsilon=0.1$ yields an almost optimal
classifier. Since the SVM training problem can be reduced to SOCP, the quantum
IPM for SOCP can be used to obtain an efficient quantum SVM algorithm. We
perform extensive numerical experiments to evaluate our algorithm on random
SVM instances and compare it against state of the art classical SOCP and SVM
solvers.

The numerical experiments on random SVM instances indicate that the running
time of the quantum algorithm scales as roughly $O(n^{2.591})$, where all the
parameters in the running time are taken into account and the exponent is
estimated using a least squares fit. We also benchmarked the exponent for
classical SVM algorithms on the same instances and for a comparable accuracy,
the scaling exponent was found to be $3.31$ for general SOCP solvers and
$3.11$ for state-of-the-art SVM solvers. We note that this does not prove a
worst-case asymptotic speedup, but the experiments on unstructured SVM
instances provide strong evidence for a significant polynomial speedup of the
quantum SVM algorithm over state-of-the-art classical SVM algorithms. We can
therefore view SVMs as a candidate problem for which quantum optimization
algorithms can achieve a polynomial speedup over state of the art classical
algorithms for an end-to-end application.

###  1.2 Related work

Our main result is the first specialized quantum algorithm for training
support-vector machines (SVM). While several quantum SVM algorithms have been
proposed, they are unlikely to offer general speedups for the most widely used
formulation of SVM – the soft-margin ($\ell_{1}$-)SVM (defined in eq. (10)).
On one hand, papers such as [3] formulate the SVM as a SDP, and solve that
using existing quantum SDP solvers such as [9, 7, 41, 40] – with the
conclusion being that a speedup is observed only for very specific sparse
instances. On the other hand, [33] solves an easier related problem – the
least-squares SVM ($\ell_{2}$-SVM or LS-SVM, see eq. (11) for its
formulation), thus losing the desirable sparsity properties of $\ell_{1}$-SVM
[38]. It turns out that applying our algorithm to this problem also yields the
same complexity as in [33] for the $\ell_{2}$-SVM. Very recently, a quantum
algorithm for yet another variant of SVM (SVM-perf, [22]) has been presented
in [2].

##  2 Preliminaries

###  2.1 Second-order cone programming

For the sake of completeness, in this section we present the most important
results about classical SOCP IPMs, from [28, 1]. We start by defining SOCP as
the optimization problem over the product of second-order (or Lorentz) cones
$\mathcal{L}=\mathcal{L}^{n_{1}}\times\cdots\times\mathcal{L}^{n_{r}}$, where
$\mathcal{L}_{k}\subseteq\mathbb{R}^{k}$ is defined as
$\mathcal{L}^{k}=\left\\{\left.\bm{x}=(x_{0};\widetilde{\bm{x}})\in\mathbb{R}^{k}\;\right\rvert\;\lVert\widetilde{\bm{x}}\rVert\leq
x_{0}\right\\}$. In this paper we consider the problem (2) and its dual (3):

$\begin{array}[]{ll}\min&\bm{c}^{\top}\bm{x}\\\ \text{s.t.}&A\bm{x}=\bm{b}\\\ &\bm{x}\in\mathcal{L},\end{array}$ (2) |  $\begin{array}[]{ll}\max&\bm{b}^{\top}\bm{y}\\\ \text{s.t.}&A^{\top}\bm{y}+\bm{s}=\bm{c}\\\ &\bm{s}\in\mathcal{L},\;\bm{y}\in\mathbb{R}^{m}.\end{array}$ (3)
---|---

We call $n:=\sum_{i=1}^{r}n_{i}$ the _size_ of the SOCP (2), and $r$ is its
_rank_.

A solution $(\bm{x},\bm{y},\bm{s})$ satisfying the constraints of both (2) and
(3) is _feasible_ , and if in addition it satisfies
$\bm{x}\in\operatorname{int}\mathcal{L}$ and
$\bm{s}\in\operatorname{int}\mathcal{L}$, it is _strictly feasible_. If at
least one constraint of (2) or (3), is violated, the solution is _infeasible_.
The duality gap of a feasible solution $(\bm{x},\bm{y},\bm{s})$ is defined as
$\mu:=\mu(\bm{x},\bm{s}):=\frac{1}{r}\bm{x}^{\top}\bm{s}$. As opposed to LP,
and similarly to SDP, strict feasibility is required for _strong duality_ to
hold [1]. From now on, we assume that our SOCP has a strictly feasible primal-
dual solution (this assumption is valid, since the homogeneous self-dual
embedding technique from [42] allows us to embed (2) and (3) in a slightly
larger SOCP where this condition is satisfied).

###  2.2 Euclidean Jordan algebras

The cone $\mathcal{L}^{n}$ has an algebraic structure similar to that of
symmetric matrices under the matrix product. Here, we consider the Jordan
product of $(x_{0},\widetilde{\bm{x}})\in\mathbb{R}^{n}$ and
$(y_{0},\widetilde{\bm{y}})\in\mathbb{R}^{n}$, defined as

$\bm{x}\circ\bm{y}:=\begin{bmatrix}\bm{x}^{T}\bm{y}\\\
x_{0}\widetilde{\bm{y}}+y_{0}\widetilde{\bm{x}}\end{bmatrix},\text{ and its
identity element }\bm{e}:=\begin{bmatrix}1\\\ 0^{n-1}\end{bmatrix}.$

This product is closely related to the (linear) matrix representation
$\operatorname{Arw}(\bm{x}):=\begin{bmatrix}x_{0}&\widetilde{\bm{x}}^{T}\\\
\widetilde{\bm{x}}&x_{0}I_{n-1}\end{bmatrix}$, which in turn satisfies the
following equality:

$\bm{x}\circ\bm{y}=\operatorname{Arw}(\bm{x})\bm{y}=\operatorname{Arw}(\bm{x})\operatorname{Arw}(\bm{y})\bm{e}.$

The key observation is that this product induces a spectral decomposition of
any vector $\bm{x}$, that has similar properties as its matrix relative.
Namely, for any vector $\bm{x}$ we define

$\displaystyle\lambda_{1}(\bm{x})$
$\displaystyle:=x_{0}+\lVert\widetilde{\bm{x}}\rVert,\;\bm{c}_{1}(\bm{x}):=\frac{1}{2}\begin{bmatrix}1\\\
\frac{\widetilde{\bm{x}}}{\lVert\widetilde{\bm{x}}\rVert}\end{bmatrix},$
$\displaystyle\lambda_{2}(\bm{x})$
$\displaystyle:=x_{0}-\lVert\widetilde{\bm{x}}\rVert,\;\bm{c}_{2}(\bm{x}):=\frac{1}{2}\begin{bmatrix}1\\\
\frac{-\widetilde{\bm{x}}}{\lVert\widetilde{\bm{x}}\rVert}\end{bmatrix}.$ (4)

We use the shorthands $\lambda_{1}:=\lambda_{1}(\bm{x})$,
$\lambda_{2}:=\lambda_{2}(\bm{x})$, $\bm{c}_{1}:=\bm{c}_{1}(\bm{x})$ and
$\bm{c}_{2}:=\bm{c}_{2}(\bm{x})$ whenever $\bm{x}$ is clear from the context,
so we observe that $\bm{x}=\lambda_{1}\bm{c}_{1}+\lambda_{2}\bm{c}_{2}$. The
set of eigenvectors $\\{\bm{c}_{1},\bm{c}_{2}\\}$ is called the _Jordan frame_
of $\bm{x}$, and satisfies several properties:

######  Proposition 1 (Properties of Jordan frames).

Let $\bm{x}\in\mathbb{R}^{n}$ and let $\\{\bm{c}_{1},\bm{c}_{2}\\}$ be its
Jordan frame. Then, the following holds:

  1. 1.

$\bm{c}_{1}\circ\bm{c}_{2}=0$ (the eigenvectors are “orthogonal”)

  2. 2.

$\bm{c}_{1}\circ\bm{c}_{1}=\bm{c}_{1}$ and
$\bm{c}_{2}\circ\bm{c}_{2}=\bm{c}_{2}$

  3. 3.

$\bm{c}_{1}$, $\bm{c}_{2}$ are of the form
$\left(\frac{1}{2};\pm\widetilde{\bm{c}}\right)$ with
$\lVert\widetilde{\bm{c}}\rVert=\frac{1}{2}$

On the other hand, just like a given matrix is positive (semi)definite if and
only if all of its eigenvalues are positive (nonnegative), a similar result
holds for $\mathcal{L}^{n}$ and $\operatorname{int}\mathcal{L}^{n}$ (the
Lorentz cone and its interior):

######  Proposition 2.

Let $\bm{x}\in\mathbb{R}^{n}$ have eigenvalues $\lambda_{1},\lambda_{2}$.
Then, the following holds:

  1. 1.

$\bm{x}\in\mathcal{L}^{n}$ if and only if $\lambda_{1}\geq 0$ and
$\lambda_{2}\geq 0$.

  2. 2.

$\bm{x}\in\operatorname{int}\mathcal{L}^{n}$ if and only if $\lambda_{1}>0$
and $\lambda_{2}>0$.

Now, using this decomposition, we can define arbitrary real powers
$\bm{x}^{p}$ for $p\in\mathbb{R}$ as
$\bm{x}^{p}:=\lambda_{1}^{p}\bm{c}_{1}+\lambda_{2}^{p}\bm{c}_{2}$, and in
particular the “inverse” and the “square root”

$\displaystyle\bm{x}^{-1}$
$\displaystyle=\frac{1}{\lambda_{1}}\bm{c}_{1}+\frac{1}{\lambda_{2}}\bm{c}_{2},\text{
if }\lambda_{1}\lambda_{2}\neq 0,$ $\displaystyle\bm{x}^{1/2}$
$\displaystyle=\sqrt{\lambda_{1}}\bm{c}_{1}+\sqrt{\lambda_{2}}\bm{c}_{2},\text{
if }\bm{x}\in\mathcal{L}^{n}.$

Moreover, we can also define some operator norms, namely the Frobenius and the
spectral one:

$\displaystyle\lVert\bm{x}\rVert_{F}$
$\displaystyle=\sqrt{\lambda_{1}^{2}+\lambda_{2}^{2}}=\sqrt{2}\lVert\bm{x}\rVert,$
$\displaystyle\lVert\bm{x}\rVert_{2}$
$\displaystyle=\max\\{|\lambda_{1}|,|\lambda_{2}|\\}=|x_{0}|+\lVert\widetilde{\bm{x}}\rVert.$

Finally, we define an analogue to the operation $Y\mapsto XYX$. It turns out
that for this we need another matrix representation (_quadratic
representation_) $Q_{\bm{x}}$, defined as

$\displaystyle Q_{\bm{x}}$
$\displaystyle:=2\operatorname{Arw}^{2}(\bm{x})-\operatorname{Arw}(\bm{x}^{2})=\begin{bmatrix}\lVert\bm{x}\rVert^{2}&2x_{0}\widetilde{\bm{x}}^{T}\\\
2x_{0}\widetilde{\bm{x}}&\lambda_{1}\lambda_{2}I_{n}+2\widetilde{\bm{x}}\widetilde{\bm{x}}^{T}\end{bmatrix}.$
(5)

Now, the matrix-vector product $Q_{\bm{x}}\bm{y}$ will behave as the quantity
$XYX$. To simplify the notation, we also define the matrix
$T_{\bm{x}}:=Q_{\bm{x}^{1/2}}$.

The definitions that we introduced so far are suitable for dealing with a
single constraint $\bm{x}\in\mathcal{L}^{n}$. For dealing with multiple
constraints
$\bm{x}_{1}\in\mathcal{L}^{n_{1}},\dots,\bm{x}_{r}\in\mathcal{L}^{n_{r}}$, we
need to deal with block-vectors
$\bm{x}=(\bm{x}_{1};\bm{x}_{2};\dots;\bm{x}_{r})$ and
$\bm{y}=(\bm{y}_{1};\bm{y}_{2};\dots;\bm{y}_{r})$. We call the number of
blocks $r$ the _rank_ of the vector (thus, up to now, we were only considering
rank-1 vectors). Now, we extend all our definitions to rank-$r$ vectors.

  1. 1.

$\bm{x}\circ\bm{y}:=(\bm{x}_{1}\circ\bm{y}_{1};\dots;\bm{x}_{r}\circ\bm{y}_{r})$

  2. 2.

The matrix representations $\operatorname{Arw}(\bm{x})$ and $Q_{\bm{x}}$ are
the block-diagonal matrices containing the representations of the blocks:

$\displaystyle\operatorname{Arw}(\bm{x})$
$\displaystyle:=\operatorname{Arw}(\bm{x}_{1})\oplus\cdots\oplus\operatorname{Arw}(\bm{x}_{r})\text{
and }Q_{\bm{x}}:=Q_{\bm{x}_{1}}\oplus\cdots\oplus Q_{\bm{x}_{r}}$

  3. 3.

$\bm{x}$ has $2r$ eigenvalues (with multiplicities) – the union of the
eigenvalues of the blocks $\bm{x}_{i}$. The eigenvectors of $\bm{x}$
corresponding to block $i$ contain the eigenvectors of $\bm{x}_{i}$ as block
$i$, and are zero everywhere else.

  4. 4.

The identity element is $\bm{e}=(\bm{e}_{1};\dots;\bm{e}_{r})$, where
$\bm{e}_{i}$’s are the identity elements for the corresponding blocks.

Thus, all things defined using eigenvalues can also be defined for rank-$r$
vectors:

  1. 1.

The norms are extended as
$\lVert\bm{x}\rVert^{2}_{F}:=\sum_{i=1}^{r}\lVert\bm{x}_{i}\rVert_{F}^{2}$ and
$\lVert\bm{x}\rVert_{2}:=\max_{i}\lVert\bm{x}_{i}\rVert_{2}$, and

  2. 2.

Powers are computed blockwise as
$\bm{x}^{p}:=(\bm{x}_{1}^{p};\dots;\bm{x}_{r}^{p})$ whenever the corresponding
blocks are defined.

Some further matrix-algebra inspired properties of block vectors are stated in
the following two claims:

######  Claim 1 (Algebraic properties).

Let $\bm{x},\bm{y}$ be two arbitrary block-vectors. Then, the following holds:

  1. 1.

The spectral norm is subadditive:
$\lVert\bm{x}+\bm{y}\rVert_{2}\leq\lVert\bm{x}\rVert_{2}+\lVert\bm{y}\rVert_{2}$.

  2. 2.

The spectral norm is less than the Frobenius norm:
$\lVert\bm{x}\rVert_{2}\leq\lVert\bm{x}\rVert_{F}$.

  3. 3.

If $A$ is a matrix with minimum and maximum singular values
$\sigma_{\text{min}}$ and $\sigma_{\text{max}}$ respectively, then the norm
$\lVert A\bm{x}\rVert$ is bounded as
$\sigma_{\text{min}}\lVert\bm{x}\rVert\leq\lVert
A\bm{x}\rVert\leq\sigma_{\text{max}}\lVert\bm{x}\rVert$.

  4. 4.

The minimum eigenvalue of $\bm{x}+\bm{y}$ is bounded as
$\lambda_{\text{min}}(\bm{x}+\bm{y})\geq\lambda_{\text{min}}(\bm{x})-\lVert\bm{y}\rVert_{2}$.

  5. 5.

The following submultiplicativity property holds:
$\lVert\bm{x}\circ\bm{y}\rVert_{F}\leq\lVert\bm{x}\rVert_{2}\cdot\lVert\bm{y}\rVert_{F}$.

In general, the proofs of these statements are analogous to the matrix case,
with a few notable differences: First, the vector spectral norm
$\lVert\cdot\rVert_{2}$ is not actually a norm, since there exist nonzero
vectors outside $\mathcal{L}$ which have zero norm. It is, however, still
bounded by the Frobenius norm (just like in the matrix case), which is in fact
a proper norm. Secondly, the minimum eigenvalue bound also holds for matrix
spectral norms, with the exact same statement. Finally, the last property is
reminiscent of the matrix submultiplicativity property $\lVert A\cdot
B\rVert_{F}\leq\lVert A\rVert_{2}\lVert B\rVert_{F}$.

We finish with several well-known properties of the quadratic representation
$Q_{\bm{x}}$ and $T_{\bm{x}}$.

######  Proposition 3 (Properties of $Q_{\bm{x}}$, from [1]).

Let $\bm{x}\in\operatorname{int}\mathcal{L}$. Then, the following holds:

  1. 1.

$Q_{\bm{x}}\bm{e}=\bm{x}^{2}$, and thus $T_{\bm{x}}\bm{e}=\bm{x}$.

  2. 2.

$Q_{\bm{x}^{-1}}=Q_{\bm{x}}^{-1}$, and more generally
$Q_{\bm{x}^{p}}=Q_{\bm{x}}^{p}$ for all $p\in\mathbb{R}$.

  3. 3.

$\lVert Q_{\bm{x}}\rVert_{2}=\lVert\bm{x}\rVert_{2}^{2}$, and thus $\lVert
T_{\bm{x}}\rVert_{2}=\lVert\bm{x}\rVert_{2}$.

  4. 4.

$Q_{\bm{x}}$ preserves $\mathcal{L}$, i.e.
$Q_{\bm{x}}(\mathcal{L})=\mathcal{L}$ and
$Q_{\bm{x}}(\operatorname{int}\mathcal{L})=\operatorname{int}\mathcal{L}$.

###  2.3 Interior-point methods

Our algorithm follows the general IPM structure, i.e. it uses Newton’s method
to solve a sequence of increasingly strict relaxations of the Karush-Kuhn-
Tucker (KKT) optimality conditions:

$\displaystyle A\bm{x}=\bm{b},\;A^{\top}\bm{y}+\bm{s}=\bm{c}$ (6)
$\displaystyle\bm{x}\circ\bm{s}=\nu\bm{e},\;\bm{x}\in\mathcal{L},\bm{s}\in\mathcal{L},$

where the parameter $\nu>0$ is decreased by a factor $\sigma<1$ in each
iteration. Since $\bm{x}\circ\bm{s}=\nu\bm{e}$ implies that the duality gap is
$\mu=\nu$, by letting $\nu\to 0$ the IPM converges towards the optimal
solution. The curve traced by (feasible) solutions $(\bm{x},\bm{y},\bm{s})$ of
(6) for $\nu>0$ is called the _central path_ , and we note that all points on
it are strictly feasible.

More precisely, in each iteration we need to find
$\Delta\bm{x},\Delta\bm{y},\Delta\bm{s}$ such that
$\bm{x}_{\text{next}}:=\bm{x}+\Delta\bm{x}$,
$\bm{y}_{\text{next}}:=\bm{y}+\Delta\bm{y}$ and
$\bm{s}_{\text{next}}:=\bm{s}+\Delta\bm{s}$ satisfy (6) for $\nu=\sigma\mu$.
After linearizing the product $\bm{x}_{\text{next}}\circ\bm{s}_{\text{next}}$,
we obtain the following linear system – the _Newton system_ :

$\displaystyle\begin{bmatrix}A&0&0\\\ 0&A^{\top}&I\\\
\operatorname{Arw}(\bm{s})&0&\operatorname{Arw}(\bm{x})\end{bmatrix}\begin{bmatrix}\Delta\bm{x}\\\
\Delta\bm{y}\\\ \Delta\bm{s}\end{bmatrix}=\begin{bmatrix}\bm{b}-A\bm{x}\\\
\bm{c}-\bm{s}-A^{\top}\bm{y}\\\
\sigma\mu\bm{e}-\bm{x}\circ\bm{s}\end{bmatrix}.$ (7)

As a final remark, it is not guaranteed that
$(\bm{x}_{\text{next}},\bm{y}_{\text{next}},\bm{s}_{\text{next}})$ is on the
central path (6), or even that it is still strictly feasible. Luckily, it can
be shown that long as $(\bm{x},\bm{y},\bm{s})$ starts out in a neighborhood
$\mathcal{N}$ of the central path,
$(\bm{x}_{\text{next}},\bm{y}_{\text{next}},\bm{s}_{\text{next}})$ will remain
both strictly feasible and in $\mathcal{N}$. It can be shown that this
algorithm halves the duality gap every $O(\sqrt{r})$ iterations, so indeed,
after $O(\sqrt{r}\log(\mu_{0}/\epsilon))$ it will converge to a (feasible)
solution with duality gap at most $\epsilon$ (given that the initial duality
gap was $\mu_{0}$).

###  2.4 Quantum linear algebra

As it was touched upon in the paper, the main speedup in our algorithms comes
from the fact that we use the quantum linear algebra algorithms from [11, 18]
whose complexity is sublinear in the dimension. Of course, we need to change
our computational model for this sentence to make any sense: namely, we encode
$n$-dimensional unit vectors as quantum states of a
$\lceil\log_{2}(n)\rceil$-qubit system. In other words, for a vector
$\bm{z}\in\mathbb{R}^{2^{k}}$ with $\lVert\bm{z}\rVert=1$, we use the notation
$\mathinner{\left\lvert{\bm{z}}\right\rangle}$ to refer to the $k$-qubit state
$\mathinner{\left\lvert{\bm{z}}\right\rangle}:=\sum_{i=0}^{2^{k}-1}z_{i}\mathinner{\left\lvert{i}\right\rangle}$,
where $\mathinner{\left\lvert{i}\right\rangle}$ are the standard basis vectors
of $\mathbb{C}^{2^{k}}$, the state space of a $k$-qubit system [32].

Given a quantum state $\mathinner{\left\lvert{\bm{z}}\right\rangle}$, we have
very limited ways of interacting with it: we can either apply a unitary
transformation $U:\mathbb{C}^{2^{k}}\to\mathbb{C}^{2^{k}}$, or we can
_measure_ it, which means that we discard the state and obtain a single random
integer $0\leq i\leq 2^{k}-1$, with the probability of measuring $i$ being
$z_{i}^{2}$. In particular this means that we can neither observe the
_amplitudes_ $z_{i}$ directly, nor can we create a copy of
$\mathinner{\left\lvert{\bm{z}}\right\rangle}$ for an arbitrary
$\mathinner{\left\lvert{\bm{z}}\right\rangle}$. In addition to this, it is _a
priori_ not clear how (and whether it is even possible) to implement the state
$\mathinner{\left\lvert{\bm{z}}\right\rangle}$ or an arbitrary unitary $U$
using the gates of a quantum computer. Luckily, there exists a quantum-
classical framework using _QRAM_ data structures described in [24] that
provides a positive answer to both of these questions.

The QRAM can be thought of as the quantum analogue to RAM, i.e. an array
$[\bm{b}^{(1)},\dots,\bm{b}^{(m)}]$ of $w$-bit bitstrings, whose elements we
can access in poly-logarithmic time given their address (position in the
array). More precisely, QRAM is just an efficient implementation of the
unitary transformation

$\mathinner{\left\lvert{i}\right\rangle}\mathinner{\left\lvert{0}\right\rangle}^{\otimes
w}\mapsto\mathinner{\left\lvert{i}\right\rangle}\mathinner{\left\lvert{b^{(i)}_{1}\dots
b^{(i)}_{w}}\right\rangle},\text{ for }i\in[m].$

The usefulness of QRAM data structures becomes clear when we consider the
_block encoding_ framework:

######  Definition.

Let $A\in\mathbb{R}^{n\times n}$ be a symmetric matrix. Then, the $\ell$-qubit
unitary matrix $U\in\mathbb{C}^{2^{\ell}\times 2^{\ell}}$ is a $(\zeta,\ell)$
block encoding of $A$ if $U=\begin{bmatrix}A/\zeta&\cdot\\\
\cdot&\cdot\end{bmatrix}$. For an arbitrary matrix $B\in\mathbb{R}^{n\times
m}$, a block encoding of $B$ is any block encoding of its symmetrized version
$\operatorname{sym}(B):=\begin{bmatrix}0&B\\\ B^{\top}&0\end{bmatrix}$.

We want $U$ to be implemented efficiently, i.e. using an $\ell$-qubit quantum
circuit of depth (poly-)logarithmic in $n$. Such a circuit would allow us to
efficiently create states $\mathinner{\left\lvert{A_{i}}\right\rangle}$
corresponding to rows (or columns) of $A$. Moreover, we need to be able to
construct such a data structure efficiently from the classical description of
$A$. It turns out that we are able to fulfill both of these requirements using
a data structure built on top of QRAM.

######  Theorem 1 (Block encodings using QRAM [24, 25]).

There exist QRAM data structures for storing vectors
$\bm{v}_{i}\in\mathbb{R}^{n}$, $i\in[m]$ and matrices $A\in\mathbb{R}^{n\times
n}$ such that with access to these data structures one can do the following:

  1. 1.

Given $i\in[m]$, prepare the state
$\mathinner{\left\lvert{\bm{v}_{i}}\right\rangle}$ in time $\widetilde{O}(1)$.
In other words, the unitary
$\mathinner{\left\lvert{i}\right\rangle}\mathinner{\left\lvert{0}\right\rangle}\mapsto\mathinner{\left\lvert{i}\right\rangle}\mathinner{\left\lvert{\bm{v}_{i}}\right\rangle}$
can be implemented efficiently.

  2. 2.

A $(\zeta(A),2\log n)$ unitary block encoding for $A$ with $\zeta(A)=\lVert
A\rVert_{2}^{-1}\min(\lVert A\rVert_{F},s_{1}(A))$, where
$s_{1}(A)=\max_{i}\sum_{j}|A_{i,j}|$ can be implemented in time
$\widetilde{O}(\log n)$. Moreover, this block encoding can be constructed in a
single pass over the matrix $A$, and it can be updated in $O(\log^{2}n)$ time
per entry.

From now on, we will also refer to storing vectors and matrices in QRAM,
meaning that we use the data structure from Theorem 1. Note that this is the
same quantum oracle model that has been used to solve SDPs in [26] and [40].

Once we have these block encodings, we may use them to perform linear algebra.
In particular, we want to construct the quantum states
$\mathinner{\left\lvert{A\bm{b}}\right\rangle}$ and
$\mathinner{\left\lvert{A^{-1}\bm{b}}\right\rangle}$, corresponding to the
matrix-vector product $A\bm{b}$ and the solution of the linear system
$A\bm{x}=\bm{b}$:

######  Theorem 2 (Quantum linear algebra with block encodings [11, 18]).

Let $A\in\mathbb{R}^{n\times n}$ be a matrix with non-zero eigenvalues in the
interval $[-1,-1/\kappa]\cup[1/\kappa,1]$, and let $\epsilon>0$. Given an
implementation of an $(\zeta,O(\log n))$ block encoding for $A$ in time
$T_{U}$ and a procedure for preparing state
$\mathinner{\left\lvert{b}\right\rangle}$ in time $T_{b}$,

  1. 1.

A state $\epsilon$-close to $\mathinner{\left\lvert{A^{-1}b}\right\rangle}$
can be generated in time
$O((T_{U}\kappa\zeta+T_{b}\kappa)\operatorname{polylog}(\kappa\zeta/\epsilon))$.

  2. 2.

A state $\epsilon$-close to $\mathinner{\left\lvert{Ab}\right\rangle}$ can be
generated in time
$O((T_{U}\kappa\zeta+T_{b}\kappa)\operatorname{polylog}(\kappa\zeta/\epsilon))$.

  3. 3.

For $\mathcal{A}\in\\{A,A^{-1}\\}$, an estimate $\Lambda$ such that
$\Lambda\in(1\pm\epsilon)\lVert\mathcal{A}b\rVert$ can be generated in time
$O((T_{U}+T_{b})\frac{\kappa\zeta}{\epsilon}\operatorname{polylog}(\kappa\zeta/\epsilon))$.

Finally, in order to recover classical information from the outputs of a
linear system solver, we require an efficient procedure for _quantum state
tomography_. The tomography procedure is linear in the dimension of the
quantum state.

######  Theorem 3 (Efficient vector state tomography, [26]).

There exists an algorithm that given a procedure for constructing
$\mathinner{\left\lvert{\bm{x}}\right\rangle}$ (i.e. a unitary mapping
$U:\mathinner{\left\lvert{0}\right\rangle}\mapsto\mathinner{\left\lvert{\bm{x}}\right\rangle}$
and its controlled version in time $T_{U}$) and precision $\delta>0$ produces
an estimate $\overline{\bm{x}}\in\mathbb{R}^{n}$ with
$\lVert\overline{\bm{x}}\rVert=1$ such that
$\lVert\bm{x}-\overline{\bm{x}}\rVert\leq\sqrt{7}\delta$ with probability at
least $(1-1/n^{0.83})$. The algorithm runs in time $O\left(T_{U}\frac{n\log
n}{\delta^{2}}\right)$.

Of course, repeating this algorithm $\widetilde{O}(1)$ times allows us to
increase the success probability to at least $1-1/\operatorname{poly}(n)$.
Putting Theorems 1, 2 and 3 together, assuming that $A$ and $\bm{b}$ are
already in QRAM, we obtain that the complexity of a completely self-contained
algorithm for solving the system $A\bm{x}=\bm{b}$ with error $\delta$ is
$\widetilde{O}\left(n\cdot\frac{\kappa\zeta}{\delta^{2}}\right)$. For well-
conditioned matrices, this presents a significant improvement over
$O(n^{\omega})$ (or, in practice, $O(n^{3})$) needed for solving linear
systems classically, especially when $n$ is large and the desired precision is
not too high. This can be compared with classical iterative linear algebra
algorithms [34] that have $O(\operatorname{nnz}(A))$ complexity per iteration,
as well as the new quantum-inspired solvers [19] that have high-degree
polynomial dependence on the rank, error, and the condition number.

##  3 A quantum interior-point method

Having introduced the classical IPM for SOCP, we can finally introduce our
quantum IPM (Algorithm 1). The main idea is to use quantum linear algebra as
much as possible (including solving the Newton system – the most expensive
part of each iteration), and falling back to classical computation when
needed. Since quantum linear algebra introduces inexactness, we need to deal
with it in the analysis. The inspiration for the “quantum part” of the
analysis is [26], whereas the “classical part” is based on [28]. Nevertheless,
the SOCP analysis is unique in many aspects and a number of hurdles had to be
overcome to the make the analysis go through. In the rest of the section, we
give a brief introduction of the most important quantum building blocks we
use, as well as present a sketch of the analysis.

Require: Matrix $A$ and vectors $\bm{b},\bm{c}$ in QRAM, precision parameter
$\epsilon$

  1. 1.

Find feasible initial point $(\bm{x},\bm{y},\bm{s},\mu)$ and store it in QRAM.

  2. 2.

Repeat the following steps for $O(\sqrt{r}\log(\mu_{0}/\epsilon))$ iterations:

    1. (a)

Compute the vector $\sigma\mu\bm{e}-\bm{x}\circ\bm{s}$ classically and store
it in QRAM.

    2. (b)

Prepare and update the block encodings of the LHS and the RHS of the Newton
system (7)

    3. (c)

Solve the Newton system to obtain
$\mathinner{\left\lvert{(\Delta\bm{x};\Delta\bm{y};\Delta\bm{s})}\right\rangle}$,
and obtain a classical approximate solution
$\overline{\left(\Delta\bm{x};\Delta\bm{y};\Delta\bm{s}\right)}$ using
tomography.

    4. (d)

Update $\bm{x}\leftarrow\bm{x}+\overline{\Delta\bm{x}}$,
$\bm{s}\leftarrow\bm{s}+\overline{\Delta\bm{s}}$ and store in QRAM.

    5. (e)

Update $\mu\leftarrow\frac{1}{r}\bm{x}^{\top}\bm{s}$.

  3. 3.

Output $(\bm{x},\bm{y},\bm{s})$.

Algorithm 1 A quantum IPM for SOCP

First, we note that the algorithms from the previous section allow us to
“forget” that Algorithm 1 is quantum, and treat it as a small modification of
the classical IPM, where the system (7) is solved up to an $\ell_{2}$-error
$\delta$. Since Algorithm 1 is iterative, the main part of the analysis is
proving that a single iteration preserves closeness to the central path,
strict feasibility, and improves the duality gap. In the remainder of this
section, we state our main results informally, while the exact statements and
proofs of all claims can be found in the supplementary material.

######  Theorem 4 (Per-iteration correctness, informal).

Let $(\bm{x},\bm{y},\bm{s})$ be a strictly feasible primal-dual solution that
is close to the central path, with duality gap $\mu$, and at distance at least
$\delta$ from the boundary of $\mathcal{L}$. Then, the Newton system (7) has a
unique solution $(\Delta\bm{x},\Delta\bm{y},\Delta\bm{s})$. There exist
positive constants $\xi,\alpha$ such that the following holds: If we let
$\overline{\Delta\bm{x}},\overline{\Delta\bm{s}}$ be approximate solutions of
(7) that satisfy

$\lVert\Delta\bm{x}-\overline{\Delta\bm{x}}\rVert_{F}\leq\xi\delta\text{ and
}\lVert\Delta\bm{s}-\overline{\Delta\bm{s}}\rVert_{F}\leq\xi\delta,$

and let $\bm{x}_{\text{next}}:=\bm{x}+\overline{\Delta\bm{x}}$ and
$\bm{s}_{\text{next}}:=\bm{s}+\overline{\Delta\bm{s}}$ be the updated
solution, then:

  1. 1.

The updated solution is strictly feasible, i.e.
$\bm{x}_{\text{next}}\in\operatorname{int}\mathcal{L}$ and
$\bm{s}_{\text{next}}\in\operatorname{int}\mathcal{L}$.

  2. 2.

The updated solution is close to the central path, and the new duality gap is
less than $(1-\alpha/\sqrt{r})\mu$.

The proof of this theorem consists of 3 main parts:

  1. 1.

Rescaling $\bm{x}$ and $\bm{s}$ so that they commute in the Jordan-algebraic
sense [1]. This part can be reused from the classical analysis [28].

  2. 2.

Bounding the norms of $\overline{\Delta\bm{x}}$ and $\overline{\Delta\bm{s}}$,
and proving that $\bm{x}+\overline{\Delta\bm{x}}$ and
$\bm{s}+\overline{\Delta\bm{s}}$ are still strictly feasible (in the sense of
belonging to $\operatorname{int}\mathcal{L}$). This part of the analysis is
also inspired by the classical analysis, but it has to take into account the
inexactness of the Newton system solution.

  3. 3.

Proving that the new solution
$(\bm{x}+\overline{\Delta\bm{x}},\bm{y}+\overline{\Delta\bm{y}},\bm{s}+\overline{\Delta\bm{s}})$
is in the neighborhood of the central path, and the duality gap/central path
parameter have decreased by a factor of $1-\alpha/\sqrt{r}$, where $\alpha$ is
constant. This part is the most technical, and while it is inspired by [26],
it required using many of the Jordan-algebraic tools from [1, 28].

Theorem 4 formalizes the fact that the Algorithm 1 has the same iteration
invariant as the classical IPM. Since the duality gap is reduced by the same
factor in both algorithms, their iteration complexity is the same, and a
simple calculation shows that they need $O(\sqrt{r})$ iterations to halve the
duality gap. On the other hand, the cost of each iteration varies, since the
complexity of the quantum linear system solver depends on the precision
$\xi\delta$, the condition number $\kappa$ of the Newton matrix, as well as
its $\zeta$-parameter. While exactly bounding these quantities is the subject
of future research, it is worth noting that for $\zeta$, we have the trivial
bound $\zeta\leq\sqrt{n}$, and research on iterative linear algebra methods
[14] suggests that $\kappa=O(1/\mu)=O(1/\epsilon)$. The final complexity is
summarized in the following theorem:

######  Theorem 5.

Let (2) be a SOCP with $A\in\mathbb{R}^{m\times n}$, $m\leq n$, and
$\mathcal{L}=\mathcal{L}^{n_{1}}\times\cdots\times\mathcal{L}^{n_{r}}$. Then,
Algorithm 1 achieves duality gap $\epsilon$ in time

$T=\widetilde{O}\left(\sqrt{r}\log\left(\mu_{0}/\epsilon\right)\cdot\frac{n\kappa\zeta}{\delta^{2}}\log\left(\frac{\kappa\zeta}{\delta}\right)\right),$

where the $\widetilde{O}(\cdot)$ notation hides the factors that are poly-
logarithmic in $n$ and $m$.

Finally, the quality of the resulting (classical) solution is characterized by
the following theorem:

######  Theorem 6.

Let (2) be a SOCP as in Theorem 5. Then, after $T$ iterations, the (linear)
infeasibility of the final iterate $\bm{x},\bm{y},\bm{s}$ is bounded as

$\displaystyle\lVert A\bm{x}-\bm{b}\rVert$ $\displaystyle\leq\delta\lVert
A\rVert,$ $\displaystyle\lVert A^{\top}\bm{y}+\bm{s}-\bm{c}\rVert$
$\displaystyle\leq\delta\left(\lVert A\rVert+1\right).$

##  4 Technical results

In this section, we present our main technical results – the proofs of
Theorems 4, 5 and 6.

###  4.1 Central path

In addition to the central path defined in (6), we define the distance from
the central path as $d(\bm{x},\bm{s},\nu)=\lVert
T_{\bm{x}}\bm{s}-\nu\bm{e}\rVert_{F}$, so the corresponding
$\eta$-neighborhood is given by

$\displaystyle\mathcal{N}_{\eta}(\nu):=\\{$
$\displaystyle(\bm{x},\bm{y},\bm{s})\;|\;(\bm{x},\bm{y},\bm{s})\;\text{strictly
feasible}\text{ and }d(\bm{x},\bm{s},\nu)\leq\eta\nu\\}.$

Using this neighborhood definition, we can specify what exactly do we mean
when we claim that the important properties of the central path are valid in
its neighborhood as well.

######  Lemma 1 (Properties of the central path).

Let $\nu>0$ be arbitrary and let
$\bm{x},\bm{s}\in\operatorname{int}\mathcal{L}$. Then, $\bm{x}$ and $\bm{s}$
satisfy the following properties:

  1. 1.

For all $\nu>0$, the duality gap and distance from the central path are
related as

$\displaystyle|\bm{x}^{\top}\bm{s}-r\nu|\leq\sqrt{\frac{r}{2}}\cdot
d(\bm{x},\bm{s},\nu).$

  2. 2.

The distance from the central path is symmetric in its arguments i.e.
$d(\bm{x},\bm{s},\nu)=d(\bm{s},\bm{x},\nu)$.

  3. 3.

Let $\mu=\frac{1}{r}\bm{x}^{\top}\bm{s}$. If
$d(\bm{x},\bm{s},\mu)\leq\eta\mu$, then
$(1+\eta)\lVert\bm{s}^{-1}\rVert_{2}\geq\lVert\mu^{-1}\bm{x}\rVert_{2}$.

###### Proof.

For part 1, let $\\{\lambda_{i}\\}_{i=1}^{2r}$ be the eigenvalues of
$T_{\bm{x}}\bm{s}$, note that $T_{x}$ is invertible as $x\in\mathcal{L}$. Then
using the properties of $T_{\bm{x}}$, we have

$\displaystyle\bm{x}^{\top}\bm{s}$
$\displaystyle=\bm{x}^{\top}T_{\bm{x}}^{-1}T_{\bm{x}}\bm{s}=(T_{\bm{x}^{-1}}\bm{x})^{\top}T_{\bm{x}}\bm{s}=\bm{e}^{\top}T_{\bm{x}}\bm{s}=\frac{1}{2}\sum_{i=1}^{2r}\lambda_{i}.$

We can therefore bound the duality gap $x^{\top}s$ as follows,

$\displaystyle\bm{x}^{\top}\bm{s}=\frac{1}{2}\sum_{i=1}^{2r}\lambda_{i}$
$\displaystyle\leq r\nu+\frac{1}{2}\sum_{i=1}^{2r}|\lambda_{i}-\nu|\leq
r\nu+\sqrt{\frac{r}{2}}\sqrt{\sum_{i=1}^{2r}(\lambda_{i}-\nu)^{2}}=r\nu+\sqrt{\frac{r}{2}}\cdot
d(\bm{x},\bm{s},\nu).$

The second step used the Cauchy-Schwarz inequality while the third follows
from the definition
$d(\bm{x},\bm{s},\nu)^{2}=\sum_{i=1}^{2r}(\lambda_{i}-\nu)^{2}$. The proof of
the lower bound is similar, but starts instead with the inequality

$\frac{1}{2}\sum_{i=1}^{2r}\lambda_{i}\geq
r\nu-\frac{1}{2}\sum_{i=1}^{2r}|\nu-\lambda_{i}|.$

For part 2, it suffices to prove that $T_{\bm{x}}\bm{s}$ and
$T_{\bm{s}}\bm{x}$ have the same eigenvalues. This follows from part 2 of
Theorem 10 in [1]. Finally for part 3, as $d(\bm{s},\bm{x},\mu)\leq\eta\mu$ we
have,

$\displaystyle\eta\mu$ $\displaystyle\geq\lVert
T_{\bm{s}}\bm{x}-\mu\bm{e}\rVert_{F}$ $\displaystyle=\lVert
T_{\bm{s}}\bm{x}-\mu\left(T_{\bm{s}}T_{\bm{s}^{-1}}\right)\bm{e}\rVert_{F}$
$\displaystyle=\lVert T_{\bm{s}}\left(\bm{x}-\mu
T_{\bm{s}^{-1}}\bm{e}\right)\rVert_{F}$
$\displaystyle\geq\lambda_{\text{min}}(T_{\bm{s}})\lVert\bm{x}-\mu
T_{\bm{s}^{-1}}\bm{e}\rVert_{F}$
$\displaystyle\geq\lambda_{\text{min}}(T_{\bm{s}})\lVert\bm{x}-\mu\bm{s}^{-1}\rVert_{2}$
$\displaystyle=\frac{1}{\lVert\bm{s}^{-1}\rVert_{2}}\cdot\mu\cdot\lVert\mu^{-1}\bm{x}-\bm{s}^{-1}\rVert_{2}$

Therefore,
$\eta\lVert\bm{s}^{-1}\rVert_{2}\geq\lVert\mu^{-1}\bm{x}-\bm{s}^{-1}\rVert_{2}$.
Finally, by the triangle inequality for the spectral norm,

$\eta\lVert\bm{s}^{-1}\rVert_{2}\geq\lVert\mu^{-1}\bm{x}\rVert_{2}-\lVert\bm{s}^{-1}\rVert_{2},$

so we can conclude that
$\lVert\bm{s}^{-1}\rVert_{2}\geq\frac{1}{1+\eta}\lVert\mu^{-1}\bm{x}\rVert_{2}$.
∎

###  4.2 A single quantum IPM iteration

Recall that the essence of our quantum algorithm is repeated solution of the
Newton system (7) using quantum linear algebra. As such, our goal is to prove
the following theorem:

######  Theorem 4 (Per-iteration correctness, formal).

Let $\chi=\eta=0.01$ and $\xi=0.001$ be positive constants and let
$(\bm{x},\bm{y},\bm{s})$ be a feasible solution of (2) and (3) with
$\mu=\frac{1}{r}\bm{x}^{\top}\bm{s}$ and $d(\bm{x},\bm{s},\mu)\leq\eta\mu$.
Then, for $\sigma=1-\chi/\sqrt{n}$, the Newton system (7) has a unique
solution $(\Delta\bm{x},\Delta\bm{y},\Delta\bm{s})$. Let
$\overline{\Delta\bm{x}},\overline{\Delta\bm{s}}$ be approximate solutions of
(7) that satisfy

$\lVert\Delta\bm{x}-\overline{\Delta\bm{x}}\rVert_{F}\leq\frac{\xi}{\lVert
T_{\bm{x}^{-1}}\rVert},\;\lVert\Delta\bm{s}-\overline{\Delta\bm{s}}\rVert_{F}\leq\frac{\xi}{2\lVert
T_{\bm{s}^{-1}}\rVert},$

where $T_{\bm{x}}$ and $T_{\bm{x}}$ are the square roots of the quadratic
representation matrices in equation (5). If we let
$\bm{x}_{\text{next}}:=\bm{x}+\overline{\Delta\bm{x}}$ and
$\bm{s}_{\text{next}}:=\bm{s}+\overline{\Delta\bm{s}}$, the following holds:

  1. 1.

The updated solution is strictly feasible, i.e.
$\bm{x}_{\text{next}}\in\operatorname{int}\mathcal{L}$ and
$\bm{s}_{\text{next}}\in\operatorname{int}\mathcal{L}$.

  2. 2.

The updated solution satisfies
$d(\bm{x}_{\text{next}},\bm{s}_{\text{next}},\overline{\mu})\leq\eta\overline{\mu}$
and
$\frac{1}{r}\bm{x}_{\text{next}}^{\top}\bm{s}_{\text{next}}=\overline{\mu}$
for $\overline{\mu}=\overline{\sigma}\mu$,
$\overline{\sigma}=1-\frac{\alpha}{\sqrt{r}}$ and a constant
$0<\alpha\leq\chi$.

Since the Newton system (7) is the same as in the classical case, we can reuse
Theorem 1 from [28] for the uniqueness part of Theorem 4. Therefore, we just
need to prove the two parts about strict feasibility and improving the duality
gap. Our analysis is inspired by the general case analysis from [5], the
derived SDP analysis from [26], and uses some technical results from the SOCP
analysis in [28]. The proof of Theorem 4 consists of three main steps:

  1. 1.

Rescaling $\bm{x}$ and $\bm{s}$ so that they share the same Jordan frame.

  2. 2.

Bounding the norms of $\Delta\bm{x}$ and $\Delta\bm{s}$, and proving that
$\bm{x}+\Delta\bm{x}$ and $\bm{s}+\Delta\bm{s}$ are still strictly feasible
(in the sense of belonging to $\operatorname{int}\mathcal{L}$).

  3. 3.

Proving that the new solution
$(\bm{x}+\Delta\bm{x},\bm{y}+\Delta\bm{y},\bm{s}+\Delta\bm{s})$ is in the
$\eta$-neighborhood of the central path, and the duality gap/central path
parameter have decreased by a factor of $1-\alpha/\sqrt{n}$, where $\alpha$ is
constant.

###  4.3 Rescaling $\bm{x}$x and $\bm{s}$s

As in the case of SDPs, the first step of the proof uses the symmetries of the
Lorentz cone to perform a commutative scaling, that is to reduce the analysis
to the case when $\bm{x}$ and $\bm{s}$ share the same Jordan frame. Although
$\circ$ is commutative by definition, two vectors sharing a Jordan frame are
akin to two matrices sharing a system of eigenvectors, and thus commuting
(some authors [1] say that the vectors _operator commute_ in this case). The
easiest way to achieve this is to scale by $T_{\bm{x}}=Q_{\bm{x}^{1/2}}$ and
$\mu^{-1}$, i.e. to change our variables as

$\bm{x}\mapsto\bm{x}^{\prime}:=T_{\bm{x}}^{-1}\bm{x}=\bm{e}\text{ and
}\bm{s}\mapsto\bm{s}^{\prime}:=\mu^{-1}T_{\bm{x}}\bm{s}.$

Note that for convenience, we have also rescaled the duality gap to 1. Recall
also that in the matrix case, the equivalent of this scaling was $X\mapsto
X^{-1/2}XX^{-1/2}=I$ and $S\mapsto\mu^{-1}X^{1/2}SX^{1/2}$. We use the
notation $\bm{z}^{\prime}$ to denote the appropriately-scaled vector $\bm{z}$,
so that we have

$\Delta\bm{x}^{\prime}:=T_{\bm{x}}^{-1}\Delta\bm{x},\quad\Delta\bm{s}^{\prime}:=\mu^{-1}T_{\bm{x}}\Delta\bm{s}$

For approximate quantities (e.g. the ones obtained using tomography, or any
other approximate linear system solver), we use the notation
$\overline{\phantom{i}\cdot\phantom{i}}$, so that the increments become
$\overline{\Delta\bm{x}}$ and $\overline{\Delta\bm{s}}$, and their scaled
counterparts are
$\overline{\Delta\bm{x}}^{\prime}:=T_{\bm{x}}^{-1}\overline{\Delta\bm{x}}$ and
$\overline{\Delta\bm{s}}^{\prime}:=\mu^{-1}T_{\bm{x}}\overline{\Delta\bm{s}}$.
Finally, we denote the scaled version of the next iterate as
$\bm{x}_{\text{next}}^{\prime}:=\bm{e}+\overline{\Delta\bm{x}}^{\prime}$ and
$\bm{s}_{\text{next}}^{\prime}:=\bm{s}^{\prime}+\overline{\Delta\bm{s}}^{\prime}$.
Now, we see that the statement of Theorem 4 implies the following bounds on
$\lVert\Delta\bm{x}^{\prime}-\overline{\Delta\bm{x}}^{\prime}\rVert_{F}$ and
$\lVert\Delta\bm{s}^{\prime}-\overline{\Delta\bm{s}}^{\prime}\rVert_{F}$:

$\displaystyle\lVert\Delta\bm{x}^{\prime}-\overline{\Delta\bm{x}}^{\prime}\rVert_{F}$
$\displaystyle=\lVert
T_{\bm{x}^{-1}}\Delta\bm{x}-T_{\bm{x}^{-1}}\overline{\Delta\bm{x}}\rVert_{F}$
$\displaystyle\leq\lVert
T_{\bm{x}^{-1}}\rVert\cdot\lVert\Delta\bm{x}-\overline{\Delta\bm{x}}\rVert_{F}\leq\xi,\text{
and }$
$\displaystyle\lVert\Delta\bm{s}^{\prime}-\overline{\Delta\bm{s}}^{\prime}\rVert_{F}$
$\displaystyle=\mu^{-1}\lVert
T_{\bm{x}}\Delta\bm{s}-T_{\bm{x}}\overline{\Delta\bm{s}}\rVert_{F}$
$\displaystyle\leq\mu^{-1}\lVert
T_{\bm{x}}\rVert\lVert\Delta\bm{s}-\overline{\Delta\bm{s}}\rVert_{F}$
$\displaystyle=\mu^{-1}\lVert\bm{x}\rVert_{2}\lVert\Delta\bm{s}-\overline{\Delta\bm{s}}\rVert_{F}$
$\displaystyle\leq(1+\eta)\lVert\bm{s}^{-1}\rVert_{2}\lVert\Delta\bm{s}-\overline{\Delta\bm{s}}\rVert_{F}\text{
by Lemma \ref{lemma:Properties of the central path}}$ $\displaystyle\leq
2\lVert
T_{\bm{s}^{-1}}\rVert\lVert\Delta\bm{s}-\overline{\Delta\bm{s}}\rVert_{F}\leq\xi.$

Throughout the analysis, we will make use of several constants: $\eta>0$ is
the distance from the central path, i.e. we ensure that our iterates stay in
the $\eta$-neighborhood $\mathcal{N}_{\eta}$ of the central path. The constant
$\sigma=1-\chi/\sqrt{r}$ is the factor by which we aim to decrease our duality
gap, for some constant $\chi>0$. Finally constant $\xi>0$ is the approximation
error for the scaled increments
$\overline{\Delta\bm{x}}^{\prime},\overline{\Delta\bm{s}}^{\prime}$. Having
this notation in mind, we can state several facts about the relation between
the duality gap and the central path distance for the original and scaled
vectors.

######  Claim 2.

The following holds for the scaled vectors $\bm{x}^{\prime}$ and
$\bm{s}^{\prime}$:

  1. 1.

The scaled duality gap is $\frac{1}{r}\bm{x}^{\prime\top}\bm{s}^{\prime}=1$.

  2. 2.

$d(\bm{x},\bm{s},\mu)\leq\eta\mu$ is equivalent to
$\lVert\bm{s}^{\prime}-\bm{e}\rVert\leq\eta$.

  3. 3.

$d(\bm{x},\bm{s},\mu\sigma)=\mu\cdot
d(\bm{x}^{\prime},\bm{s}^{\prime},\sigma)$, for all $\sigma>0$.

At this point, we claim that it suffices to prove the two parts of Theorem 4
in the scaled case. Namely, assuming that
$\bm{x}_{\text{next}}^{\prime}\in\operatorname{int}\mathcal{L}$ and
$\bm{s}_{\text{next}}^{\prime}\in\operatorname{int}\mathcal{L}$, by
construction we get

$\bm{x}_{\text{next}}=T_{\bm{x}}\bm{x}_{\text{next}}^{\prime}\text{ and
}\bm{s}_{\text{next}}=T_{\bm{x}^{-1}}\bm{s}_{\text{next}}^{\prime}$

and thus
$\bm{x}_{\text{next}},\bm{s}_{\text{next}}\in\operatorname{int}\mathcal{L}$.

On the other hand, if $\mu
d(\bm{x}_{\text{next}}^{\prime},\bm{s}_{\text{next}}^{\prime},\overline{\sigma})\leq\eta\overline{\mu}$,
then
$d(\bm{x}_{\text{next}},\bm{s}_{\text{next}},\overline{\mu})\leq\eta\overline{\mu}$
follows by Claim 2. Similarly, from
$\frac{1}{r}\bm{x}_{\text{next}}^{\prime\top}\bm{s}_{\text{next}}^{\prime}=\overline{\sigma}$,
we also get
$\frac{1}{r}\bm{x}_{\text{next}}^{\top}\bm{s}_{\text{next}}=\overline{\mu}$.
We conclude this part with two technical results from [28], that use the
auxiliary matrix $R_{xs}$ defined as
$R_{xs}:=T_{\bm{x}}\operatorname{Arw}(\bm{x})^{-1}\operatorname{Arw}(\bm{s})T_{\bm{x}}$.
These results are useful for the later parts of the proof of Theorem 4.

######  Claim 3 ([28], Lemma 3).

Let $\eta$ be the distance from the central path, and let $\nu>0$ be
arbitrary. Then, $R_{xs}$ is bounded as

$\lVert R_{xs}-\nu I\rVert\leq 3\eta\nu.$

######  Claim 4 ([28], Lemma 5, proof).

Let $\mu$ be the duality gap. Then, the scaled increment
$\Delta\bm{s}^{\prime}$ is

$\Delta\bm{s}^{\prime}=\sigma\bm{e}-\bm{s}^{\prime}-\mu^{-1}R_{xs}\Delta\bm{x}^{\prime}.$

###  4.4 Maintaining strict feasibility

The main tool for showing that strict feasibility is conserved is the
following bound on the increments $\Delta\bm{x}^{\prime}$ and
$\Delta\bm{s}^{\prime}$:

######  Lemma 2 ([28], Lemma 6).

Let $\eta$ be the distance from the central path and let $\mu$ be the duality
gap. Then, we have the following bounds for the scaled direction:

$\begin{array}[]{rl}\lVert\Delta\bm{x}^{\prime}\rVert_{F}&\leq\frac{\Theta}{\sqrt{2}}\\\
\lVert\Delta\bm{s}^{\prime}\rVert_{F}&\leq\Theta\sqrt{2}\end{array},\quad\text{where}\quad\Theta=\frac{2\sqrt{\eta^{2}/2+(1-\sigma)^{2}r}}{1-3\eta}$

Moreover, if we substitute $\sigma$ with its actual value $1-\chi/\sqrt{r}$,
we get $\Theta=\frac{\sqrt{2\eta^{2}+4\chi^{2}}}{1-3\eta}$, which we can make
arbitrarily small by tuning the constants. Now, we can immediately use this
result to prove
$\bm{x}_{\text{next}}^{\prime},\bm{s}_{\text{next}}^{\prime}\in\operatorname{int}\mathcal{L}$.

######  Lemma 3.

Let $\eta=\chi=0.01$ and $\xi=0.001$. Then, $\bm{x}_{\text{next}}^{\prime}$
and $\bm{s}_{\text{next}}^{\prime}$ are strictly feasible, i.e.
$\bm{x}_{\text{next}}^{\prime},\bm{s}_{\text{next}}^{\prime}\in\operatorname{int}\mathcal{L}$.

###### Proof.

By Lemma 2, $\lambda_{\text{min}}(\overline{\bm{x}})\geq
1-\lVert\overline{\Delta\bm{x}}^{\prime}\rVert_{F}\geq
1-\frac{\Theta}{\sqrt{2}}-\xi$. On the other hand, since
$d(\bm{x},\bm{s},\mu)\leq\eta\mu$, we have
$d(\bm{x}^{\prime},\bm{s}^{\prime},1)\leq\eta$, and thus

$\displaystyle\eta^{2}$
$\displaystyle\geq\lVert\bm{s}^{\prime}-e\rVert_{F}^{2}=\sum_{i=1}^{2r}(\lambda_{i}(\bm{s}^{\prime})-1)^{2}$

The above equation implies that
$\lambda_{i}(\bm{s}^{\prime})\in\left[1-\eta,1+\eta\right],\forall i\in[2r]$.
Now, since $\lVert\bm{z}\rVert_{2}\leq\lVert\bm{z}\rVert_{F}$,

$\displaystyle\lambda_{\text{min}}(\overline{\bm{s}})$
$\displaystyle\geq\lambda_{\text{min}}(\bm{s}^{\prime}+\Delta\bm{s}^{\prime})-\lVert\overline{\Delta\bm{s}}^{\prime}-\Delta\bm{s}^{\prime}\rVert_{F}$
$\displaystyle\geq\lambda_{\text{min}}(\bm{s}^{\prime})-\lVert\Delta\bm{s}^{\prime}\rVert_{F}-\lVert\overline{\Delta\bm{s}}^{\prime}-\Delta\bm{s}^{\prime}\rVert_{F}$
$\displaystyle\geq 1-\eta-\Theta\sqrt{2}-\xi,$

where we used Lemma 2 for the last inequality. Substituting $\eta=\chi=0.01$
and $\xi=0.001$, we get that $\lambda_{\text{min}}(\overline{\bm{x}})\geq 0.8$
and $\lambda_{\text{min}}(\overline{\bm{s}})\geq 0.8$. ∎

###  4.5 Maintaining closeness to central path

Finally, we move on to the most technical part of the proof of Theorem 4,
where we prove that
$\bm{x}_{\text{next}}^{\prime},\bm{s}_{\text{next}}^{\prime}$ is still close
to the central path, and the duality gap has decreased by a constant factor.
We split this into two lemmas.

######  Lemma 4.

Let $\eta=\chi=0.01$, $\xi=0.001$, and let $\alpha$ be any value satisfying
$0<\alpha\leq\chi$. Then, for $\overline{\sigma}=1-\alpha/\sqrt{r}$, the
distance to the central path is maintained, that is,
$d(\bm{x}_{\text{next}}^{\prime},\bm{s}_{\text{next}}^{\prime},\overline{\sigma})<\eta\overline{\sigma}$.

###### Proof.

By Claim 2, the distance of the next iterate from the central path is

$d(\bm{x}_{\text{next}}^{\prime},\bm{s}_{\text{next}}^{\prime},\overline{\sigma})=\lVert
T_{\bm{x}_{\text{next}}^{\prime}}\bm{s}_{\text{next}}^{\prime}-\overline{\sigma}\bm{e}\rVert_{F},$

and we can bound it from above as

$\displaystyle
d(\bm{x}_{\text{next}}^{\prime},\bm{s}_{\text{next}}^{\prime},\overline{\sigma})$
$\displaystyle=\lVert
T_{\bm{x}_{\text{next}}^{\prime}}\bm{s}_{\text{next}}^{\prime}-\overline{\sigma}\bm{e}\rVert_{F}$
$\displaystyle=\lVert
T_{\bm{x}_{\text{next}}^{\prime}}\bm{s}_{\text{next}}^{\prime}-\overline{\sigma}T_{\bm{x}_{\text{next}}^{\prime}}T_{\bm{x}_{\text{next}}^{\prime-1}}\bm{e}\rVert_{F}$
$\displaystyle\leq\lVert
T_{\bm{x}_{\text{next}}^{\prime}}\rVert\cdot\lVert\bm{s}_{\text{next}}^{\prime}-\overline{\sigma}\cdot(\bm{x}_{\text{next}}^{\prime})^{-1}\rVert.$

So, it is enough to bound
$\lVert\bm{z}\rVert_{F}:=\lVert\bm{s}_{\text{next}}^{\prime}-\overline{\sigma}\cdot\bm{x}_{\text{next}}^{\prime-1}\rVert_{F}$
from above, since

$\displaystyle\lVert T_{\bm{x}_{\text{next}}^{\prime}}\rVert$
$\displaystyle=\lVert\bm{x}_{\text{next}}^{\prime}\rVert_{2}\leq
1+\lVert\overline{\Delta\bm{x}}^{\prime}\rVert_{2}\leq
1+\lVert\Delta\bm{x}^{\prime}\rVert_{2}+\xi\leq
1+\frac{\Theta}{\sqrt{2}}+\xi.$

We split $\bm{z}$ as

$\displaystyle\bm{z}$
$\displaystyle=\underbrace{\left(\bm{s}^{\prime}+\overline{\Delta\bm{s}}^{\prime}-\overline{\sigma}e+\overline{\Delta\bm{x}}^{\prime}\right)}_{\bm{z}_{1}}+\underbrace{(\overline{\sigma}-1)\overline{\Delta\bm{x}}^{\prime}}_{\bm{z}_{2}}+\underbrace{\overline{\sigma}\left(e-\overline{\Delta\bm{x}}^{\prime}-(e+\overline{\Delta\bm{x}}^{\prime})^{-1}\right)}_{\bm{z}_{3}},$

and we bound $\lVert\bm{z}_{1}\rVert_{F},\lVert\bm{z}_{2}\rVert_{F}$, and
$\lVert\bm{z}_{3}\rVert_{F}$ separately.

  1. 1.

By the triangle inequality,
$\lVert\bm{z}_{1}\rVert_{F}\leq\lVert\bm{s}^{\prime}+\Delta\bm{s}^{\prime}-\overline{\sigma}e+\Delta\bm{x}^{\prime}\rVert_{F}+2\xi$.
Furthermore, after substituting $\Delta\bm{s}^{\prime}$ from Claim 4, we get

$\displaystyle\bm{s}^{\prime}+\Delta\bm{s}^{\prime}-\overline{\sigma}e+\Delta\bm{x}^{\prime}$
$\displaystyle=\sigma
e-\mu^{-1}R_{xs}\Delta\bm{x}^{\prime}-\overline{\sigma}e+\Delta\bm{x}^{\prime}$
$\displaystyle=\frac{\alpha-\chi}{\sqrt{r}}e+\mu^{-1}(\mu
I-R_{xs})\Delta\bm{x}^{\prime}.$

Using the bound for $\lVert\mu I-R_{xs}\rVert$ from Claim 3 as well as the
bound for $\lVert\Delta\bm{x}^{\prime}\rVert_{F}$ from Lemma 2, we obtain

$\displaystyle\lVert\bm{z}_{1}\rVert_{F}\leq
2\xi+\frac{\chi}{\sqrt{r}}+\frac{3}{\sqrt{2}}\eta\Theta.$

  2. 2.

$\lVert\bm{z}_{2}\rVert_{F}\leq\frac{\chi}{\sqrt{r}}\left(\frac{\Theta}{\sqrt{2}}+\xi\right)$,
where we used the bound from Lemma 2 again.

  3. 3.

Here, we first need to bound
$\lVert(e+\Delta\bm{x}^{\prime})^{-1}-(e+\overline{\Delta\bm{x}}^{\prime})^{-1}\rVert_{F}$.
For this, we use the submultiplicativity of $\lVert\cdot\rVert_{F}$:

$\displaystyle\lVert$
$\displaystyle(e+\Delta\bm{x}^{\prime})^{-1}-(e+\overline{\Delta\bm{x}}^{\prime})^{-1}\rVert_{F}=\lVert(e+\Delta\bm{x}^{\prime})^{-1}\circ\left(e-(e+\Delta\bm{x}^{\prime})\circ(e+\overline{\Delta\bm{x}}^{\prime})^{-1}\right)\rVert_{F}$
$\displaystyle\leq\lVert(e+\Delta\bm{x}^{\prime})^{-1}\rVert_{2}\cdot\lVert
e-(e+\overline{\Delta\bm{x}}^{\prime}+\Delta\bm{x}^{\prime}-\overline{\Delta\bm{x}}^{\prime})\circ(e+\overline{\Delta\bm{x}}^{\prime})^{-1}\rVert_{F}$
$\displaystyle=\lVert(e+\Delta\bm{x}^{\prime})^{-1}\rVert_{2}\cdot\lVert(\Delta\bm{x}^{\prime}-\overline{\Delta\bm{x}}^{\prime})\circ(e+\overline{\Delta\bm{x}}^{\prime})^{-1}\rVert_{F}$
$\displaystyle\leq\lVert(e+\Delta\bm{x}^{\prime})^{-1}\rVert_{2}\cdot\lVert\Delta\bm{x}^{\prime}-\overline{\Delta\bm{x}}^{\prime}\rVert_{F}\cdot\lVert(e+\overline{\Delta\bm{x}}^{\prime})^{-1}\rVert_{2}$
$\displaystyle\leq\xi\cdot\lVert(e+\Delta\bm{x}^{\prime})^{-1}\rVert_{2}\cdot\lVert(e+\overline{\Delta\bm{x}}^{\prime})^{-1}\rVert_{2}.$

Now, we have the bound
$\lVert(e+\Delta\bm{x}^{\prime})^{-1}\rVert_{2}\leq\frac{1}{1-\lVert\Delta\bm{x}^{\prime}\rVert_{F}}$
and similarly
$\lVert(e+\overline{\Delta\bm{x}}^{\prime})^{-1}\rVert_{2}\leq\frac{1}{1-\lVert\Delta\bm{x}^{\prime}\rVert_{F}-\xi}$,
so we get

$\lVert(e+\Delta\bm{x}^{\prime})^{-1}-(e+\overline{\Delta\bm{x}}^{\prime})^{-1}\rVert_{F}\leq\frac{\xi}{(1-\lVert\Delta\bm{x}^{\prime}\rVert_{F}-\xi)^{2}}.$

Using this, we can bound $\lVert\bm{z}_{3}\rVert_{F}$:

$\displaystyle\lVert\bm{z}_{3}\rVert_{F}\leq\overline{\sigma}$
$\displaystyle\left(\lVert
e-\Delta\bm{x}^{\prime}-(e+\Delta\bm{x}^{\prime})^{-1}\rVert_{F}+\xi+\frac{\xi}{(1-\lVert\Delta\bm{x}^{\prime}\rVert_{F}-\xi)^{2}}\right).$

If we let $\lambda_{i}$ be the eigenvalues of $\Delta\bm{x}^{\prime}$, then by
Lemma 2, we have

$\displaystyle\lVert$ $\displaystyle
e-\Delta\bm{x}^{\prime}-(e+\Delta\bm{x}^{\prime})^{-1}\rVert_{F}$
$\displaystyle=\sqrt{\sum_{i=1}^{2r}\left((1-\lambda_{i})-\frac{1}{1+\lambda_{i}}\right)^{2}}$
$\displaystyle=\sqrt{\sum_{i=1}^{2r}\frac{\lambda_{i}^{4}}{(1+\lambda_{i})^{2}}}\leq\frac{\Theta}{\sqrt{2}-\Theta}\sqrt{\sum_{i=1}^{2r}\lambda_{i}^{2}}$
$\displaystyle\leq\frac{\Theta^{2}}{2-\sqrt{2}\Theta}.$

Combining all bound from above, we obtain

$\displaystyle
d(\bm{x}_{\text{next}}^{\prime},\bm{s}_{\text{next}}^{\prime},\overline{\sigma})\leq\left(1+\frac{\Theta}{\sqrt{2}}+\xi\right)\cdot$
$\displaystyle\left(2\xi+\frac{\chi}{\sqrt{r}}+\frac{3}{\sqrt{2}}\eta\Theta\right.$
$\displaystyle+\frac{\chi}{\sqrt{r}}\left(\frac{\Theta}{\sqrt{2}}+\xi\right)$
$\displaystyle+\left.\overline{\sigma}\left(\frac{\Theta^{2}}{2-\sqrt{2}\Theta}+\xi+\frac{\xi}{(1-\Theta/\sqrt{2}-\xi)^{2}}\right)\right).$

Finally, if we plug in $\chi=0.01$, $\eta=0.01$, $\xi=0.001$, we get
$d(\overline{\bm{x}},\overline{\bm{s}},\overline{\sigma})\leq
0.005\overline{\sigma}\leq\eta\overline{\sigma}$. ∎

Now, we prove that the duality gap decreases.

######  Lemma 5.

For the same constants, the updated solution satisfies
$\frac{1}{r}\bm{x}_{\text{next}}^{\prime\top}\bm{s}_{\text{next}}^{\prime}=\left(1-\frac{\alpha}{\sqrt{r}}\right)$
for $\alpha=0.005$.

###### Proof.

Since $\bm{x}_{\text{next}}^{\prime}$ and $\bm{s}_{\text{next}}^{\prime}$ are
scaled quantities, the duality gap between their unscaled counterparts is
$\frac{\mu}{r}\bm{x}_{\text{next}}^{\prime\top}\bm{s}_{\text{next}}^{\prime}$.
Applying Lemma 1 (and Claim 2) with $\nu=\sigma\mu$ and
$\bm{x}_{\text{next}}^{\prime},\bm{s}_{\text{next}}^{\prime}$, we obtain the
upper bound

$\displaystyle\mu\bm{x}_{\text{next}}^{\prime\top}\bm{s}_{\text{next}}^{\prime}\leq
r\sigma\mu+\sqrt{\frac{r}{2}}\mu
d(\bm{x}_{\text{next}}^{\prime},\bm{s}_{\text{next}}^{\prime},\sigma),$

which in turn implies

$\frac{1}{r}\bm{x}_{\text{next}}^{\prime\top}\bm{s}_{\text{next}}^{\prime}\leq\left(1-\frac{0.01}{\sqrt{r}}\right)\left(1+\frac{d(\bm{x}_{\text{next}}^{\prime},\bm{s}_{\text{next}}^{\prime},\sigma)}{\sigma\sqrt{2r}}\right).$

By instantiating Lemma 1 for $\alpha=\chi$, from its proof, we obtain
$d(\bm{x}_{\text{next}}^{\prime},\bm{s}_{\text{next}}^{\prime},\sigma)\leq
0.005\sigma$, and thus

$\frac{1}{r}\bm{x}_{\text{next}}^{\prime\top}\bm{s}_{\text{next}}^{\prime}\leq
1-\frac{0.005}{\sqrt{r}}$

Therefore, the final $\alpha$ for this Lemma is $0.005$. ∎

###  4.6 Final complexity and feasibility

In every iteration we need to solve the Newton system to a precision dependent
on the norms of $T_{\bm{x}^{-1}}$ and $T_{\bm{s}^{-1}}$. Thus, to bound the
running time of the algorithm (since the complexity of the vector state
tomography procedure depends on the desired precision), we need to bound
$\lVert T_{\bm{x}^{-1}}\rVert$ and $\lVert T_{\bm{s}^{-1}}\rVert$. Indeed, by
the properties of the quadratic representation, we get

$\displaystyle\lVert T_{\bm{x}^{-1}}\rVert$
$\displaystyle=\lVert\bm{x}^{-1}\rVert=\lambda_{\text{min}}(\bm{x})^{-1}\text{
and }$ $\displaystyle\lVert T_{\bm{s}^{-1}}\rVert$
$\displaystyle=\lVert\bm{s}^{-1}\rVert=\lambda_{\text{min}}(\bm{s})^{-1}.$

If the tomography precision for iteration $i$ is chosen to be at least (i.e.
smaller than)

$\delta_{i}:=\frac{\xi}{4}\min\left\\{\lambda_{\text{min}}(\bm{x}_{i}),\lambda_{\text{min}}(\bm{s}_{i})\right\\},$
(8)

then the premises of Theorem 4 are satisfied. The tomography precision for the
entire algorithm can therefore be chosen to be $\delta:=\min_{i}\delta_{i}$.
Note that these minimum eigenvalues are related to how close the current
iterate is to the boundary of $\mathcal{L}$ – as long as
$\bm{x}_{i},\bm{s}_{i}$ are not “too close” to the boundary of $\mathcal{L}$,
their minimal eigenvalues should not be “too small”.

There are two more parameters that impact the complexity of the quantum linear
system solver: the condition number of the Newton matrix $\kappa_{i}$ and the
matrix parameter $\zeta_{i}$ of the QRAM encoding in iteration $i$. For both
of these quantities we define their global versions as
$\kappa=\max_{i}\kappa_{i}$ and $\zeta=\max_{i}\zeta_{i}$. Therefore, we
arrive to the following statement about the complexity of Algorithm 1.

######  Theorem 5 (restatement).

Let (2) be a SOCP with $A\in\mathbb{R}^{m\times n}$, $m\leq n$, and
$\mathcal{L}=\mathcal{L}^{n_{1}}\times\cdots\times\mathcal{L}^{n_{r}}$. Then,
our algorithm achieves duality gap $\epsilon$ in time

$T=\widetilde{O}\left(\sqrt{r}\log\left(\mu_{0}/\epsilon\right)\cdot\frac{n\kappa\zeta}{\delta^{2}}\log\left(\frac{\kappa\zeta}{\delta}\right)\right).$

This complexity can be easily interpreted as product of the number of
iterations and the cost of $n$-dimensional vector tomography with error
$\delta$. So, improving the complexity of the tomography algorithm would
improve the running time of our algorithm as well.

Note that up to now, we cared mostly about strict (conic) feasibility of
$\bm{x}$ and $\bm{s}$. Now, we address the fact that the linear constraints
$A\bm{x}=\bm{b}$ and $A^{\top}\bm{y}+\bm{s}=\bm{c}$ are not exactly satisfied
during the execution of the algorithm. Luckily, it turns out that this error
is not accumulated, but is instead determined just by the final tomography
precision:

######  Theorem 6 (restatement).

Let (2) be a SOCP as in Theorem 5. Then, after $T$ iterations, the (linear)
infeasibility of the final iterate $\bm{x},\bm{y},\bm{s}$ is bounded as

$\displaystyle\lVert A\bm{x}_{T}-\bm{b}\rVert$ $\displaystyle\leq\delta\lVert
A\rVert,$ $\displaystyle\lVert A^{\top}\bm{y}_{T}+\bm{s}_{T}-\bm{c}\rVert$
$\displaystyle\leq\delta\left(\lVert A\rVert+1\right).$

###### Proof.

Let $(\bm{x}_{T},\bm{y}_{T},\bm{s}_{T})$ be the $T$-th iterate. Then, the
following holds for $A\bm{x}_{T}-\bm{b}$:

$A\bm{x}_{T}-\bm{b}=A\bm{x}_{0}+A\sum_{t=1}^{\top}\overline{\Delta\bm{x}}_{t}-\bm{b}=A\sum_{t=1}^{\top}\overline{\Delta\bm{x}}_{t}.$
(9)

On the other hand, the Newton system at iteration $T$ has the constraint
$A\Delta\bm{x}_{T}=\bm{b}-A\bm{x}_{T-1}$, which we can further recursively
transform as,

$\displaystyle A\Delta\bm{x}_{T}$
$\displaystyle=\bm{b}-A\bm{x}_{T-1}=\bm{b}-A\left(\bm{x}_{T-2}+\overline{\Delta\bm{x}}_{T-1}\right)$
$\displaystyle=\bm{b}-A\bm{x}_{0}-\sum_{t=1}^{T-1}\overline{\Delta\bm{x}}_{t}=-\sum_{t=1}^{T-1}\overline{\Delta\bm{x}}_{t}.$

Substituting this into equation (9), we get

$A\bm{x}_{T}-\bm{b}=A\left(\overline{\Delta\bm{x}}_{T}-\Delta\bm{x}_{T}\right).$

Similarly, using the constraint
$A^{\top}\Delta\bm{y}_{T}+\Delta\bm{s}_{T}=\bm{c}-\bm{s}_{T-1}-A^{\top}\bm{y}_{T-1}$
we obtain that

$A^{\top}\bm{y}_{T}+\bm{s}_{T}-\bm{c}=A^{\top}\left(\overline{\Delta\bm{y}}_{T}-\Delta\bm{y}_{T}\right)+\left(\overline{\Delta\bm{s}}_{T}-\Delta\bm{s}_{T}\right).$

Finally, we can bound the norms of these two quantities,

$\displaystyle\lVert A\bm{x}_{T}-\bm{b}\rVert$ $\displaystyle\leq\delta\lVert
A\rVert,$ $\displaystyle\lVert A^{\top}\bm{y}_{T}+\bm{s}_{T}-\bm{c}\rVert$
$\displaystyle\leq\delta\left(\lVert A\rVert+1\right).$

∎

##  5 Quantum Support-Vector Machines

In this section we present our quantum support vector machine (SVM) algorithm
as an application of our SOCP solver. Given a set of vectors
$\mathcal{X}=\\{\bm{x}^{(i)}\in\mathbb{R}^{n}\;|\;i\in[m]\\}$ (_training
examples_) and their _labels_ $y^{(i)}\in\\{-1,1\\}$, the objective of the SVM
training process is to find the “best” hyperplane that separates training
examples with label $1$ from those with label $-1$. In this paper we focus on
the (traditional) _soft-margin_ ($\ell_{1}$-)SVM, which can be expressed as
the following optimization problem:

$\begin{array}[]{ll}\min\limits_{\bm{w},b,\bm{\xi}}&\lVert\bm{w}\rVert^{2}+C\lVert\bm{\xi}\rVert_{1}\\\
\text{s.t.}&y^{(i)}(\bm{w}^{\top}\bm{x}^{(i)}+b)\geq 1-\xi_{i},\;\forall
i\in[m]\\\ &\bm{\xi}\geq 0.\end{array}$ (10)

Here, the variables $\bm{w}\in\mathbb{R}^{n}$ and $b\in\mathbb{R}$ correspond
to the hyperplane, $\bm{\xi}\in\mathbb{R}^{m}$ corresponds to the “linear
inseparability” of each point, and the constant $C>0$ is a hyperparameter that
quantifies the tradeoff between maximizing the margin and minimizing the
constraint violations.

As a slightly less traditional alternative, one might also consider the
$\ell_{2}$-SVM (or least-squares SVM, LS-SVM) [37], where the
$\lVert\bm{\xi}\rVert_{1}$ regularization term is replaced by
$\lVert\bm{\xi}\rVert^{2}$. This formulation arises from considering the
least-squares regression problem with the constraints
$y^{(i)}(\bm{w}^{\top}\bm{x}^{(i)}+b)=1$, which we solve by minimizing the
squared $2$-norm of the residuals:

$\begin{array}[]{ll}\min\limits_{\bm{w},b,\bm{\xi}}&\lVert\bm{w}\rVert^{2}+C\lVert\bm{\xi}\rVert^{2}\\\
\text{s.t.}&y^{(i)}(\bm{w}^{\top}\bm{x}^{(i)}+b)=1-\xi_{i},\;\forall
i\in[m]\end{array}$ (11)

Since this is a least-squares problem, the optimal $\bm{w},b$ and $\bm{\xi}$
can be obtained by solving a linear system. In [33], a quantum algorithm for
LS-SVM is presented, which uses a single quantum linear system solver.
Unfortunately, replacing the $\ell_{1}$-norm with $\ell_{2}$ in the objective
of (11) leads to the loss of a key property of ($\ell_{1}$-)SVM – weight
sparsity [38].

###  5.1 Reducing SVM to SOCP

Finally, we are going to reduce the SVM problem (10) to SOCP. In order to do
that, we define an auxiliary vector $\bm{t}=\left(t+1;t;\bm{w}\right)$, where
$t\in\mathbb{R}$ – this allows us to “compute” $\lVert\bm{w}\rVert^{2}$ using
the constraint $\bm{t}\in\mathcal{L}^{n+2}$ since

$\bm{t}\in\mathcal{L}^{n+2}\Leftrightarrow(t+1)^{2}\geq
t^{2}+\lVert\bm{w}\rVert^{2}\Leftrightarrow 2t+1\geq\lVert\bm{w}\rVert^{2}.$

Thus, minimizing $\lVert\bm{w}\rVert^{2}$ is equivalent to minimizing $t$.
Note we can restrict our bias $b$ to be nonnegative without any loss in
generality, since the case $b<0$ can be equivalently described by a bias
$-b>0$ and weights $-\bm{w}$. Using these transformations, we can restate (10)
as the following SOCP:

$\begin{array}[]{ll}\min\limits_{\bm{t},b,\bm{\xi}}&\begin{bmatrix}0&1&0^{n}&0&C^{m}\end{bmatrix}\begin{bmatrix}\bm{t}&b&\bm{\xi}\end{bmatrix}^{\top}\\\
\text{s.t.}&\begin{bmatrix}0&0&&1&\\\
\vdots&\vdots&X^{\top}&\vdots&\operatorname{diag}(\bm{y})\\\ 0&0&&1&\\\
1&-1&0^{n}&0&0^{m}\end{bmatrix}\begin{bmatrix}\bm{t}\\\ b\\\
\bm{\xi}\end{bmatrix}=\begin{bmatrix}\bm{y}\\\ 1\end{bmatrix}\\\
&\bm{t}\in\mathcal{L}^{n+2},\;b\in\mathcal{L}^{1},\;\xi_{i}\in\mathcal{L}^{1}\quad\forall
i\in[m]\end{array}$ (12)

Here, we use the notation $X\in\mathbb{R}^{n\times m}$ for the matrix whose
columns are the training examples $\bm{x}^{(i)}$, and
$\bm{y}\in\mathbb{R}^{m}$ for the vector of labels. This problem has $O(n+m)$
variables, and $O(m)$ conic constraints (i.e. its rank is $r=O(m)$).
Therefore, in the interesting case of $m=\Theta(n)$, it can be solved in
$\widetilde{O}(\sqrt{n})$ iterations. More precisely, if we consider both the
primal and the dual, in total they have $3m+2n+7$ scalar variables and $2m+4$
conic constraints.

In practice (as evidenced by the LIBSVM and LIBLINEAR libraries [12, 16]), a
small modification is made to the formulations (10) and (11): instead of
treating the bias separately, all data points are extended with a constant
unit coordinate. In this case, the SOCP formulation remains almost identical,
with the only difference being that the constraints
$\bm{t}\in\mathcal{L}^{n+2}$ and $b\in\mathcal{L}^{1}$ are replaced by a
single conic constraint $(\bm{t};b)\in\mathcal{L}^{n+3}$. This change allows
us to come up with a simple feasible initial solution in our numerical
experiments, without going through the homogeneous self-dual formalism of
[42].

Note also that we can solve the LS-SVM problem (11), by reducing it to a SOCP
in a similar manner. In fact, this would have resulted in just $O(1)$ conic
constraints, so an IPM would converge to a solution in $\widetilde{O}(1)$
iterations, which is comparable with the result from [33].

###  5.2 Experimental results

We next present some experimental results to assess the running time
parameters and the performance of our algorithm for random instances of SVM.
If an algorithm demonstrates a speedup on unstructured instances like these,
it is reasonable to extrapolate that the speedup is generic, as it could not
have used any special properties of the instance to derive an advantage. For a
given dimension $n$ and number of training points $m$, we denote our
distribution of random SVMs with $\mathcal{SVM}(n,m,p)$, where $p$ denotes the
probability that a datapoint is misclassified by the optimal separating
hyperplane. Additionally, for every training set sampled from
$\mathcal{SVM}(n,m,p)$, a corresponding test set of size $\lfloor m/3\rfloor$
was also sampled from the same distribution. These test sets are used to
evaluate the generalization error of SVMs trained in various ways.

Our experiments consist of generating roughly 16000 instances of
$\mathcal{SVM}(n,2n,p)$, where $n$ is chosen to be uniform between $2^{2}$ and
$2^{9}$ and $p$ is chosen uniformly from the discrete set
$\\{0,0.1,\dots,0.9,1\\}$. The instances are then solved using a simulation of
Algorithm 1 (with the target duality gap of $\epsilon=0.1$) as well as using
the ECOS SOCP solver [15] (with the default target duality gap). We simulate
the execution of Algorithm 1 by implementing the classical IPM and adding
noise to the solution of the Newton system (7). The noise added to each
coordinate is uniform, from an interval selected so that the noisy increment
$(\Delta\bm{x},\Delta\bm{y},\Delta\bm{s})$ simulates the outputs of the
tomography algorithm with precision determined by Theorem 4. The SVM parameter
$C$ is set to be equal to $1$ in all experiments. Additionally, a separate,
smaller experiment with roughly 1000 instances following the same distribution
is performed for comparing Algorithm 1 with LIBSVM [12] using a linear kernel.

The experiments are performed on a Dell Precision 7820T workstation with two
Intel Xeon Silver 4110 CPUs and 64GB of RAM, and experiment logs are available
at [39].

Figure 1: Observed complexity of Algorithm 1, its power law fit, and 95%
confidence interval.
Figure 2: Empirical CDF of the difference in accuracy between SVMs trained in
different ways.

By finding the least-squares fit of the power law $y=ax^{b}$ through the
observed values of the quantity $\frac{n^{1.5}\kappa\zeta}{\delta^{2}}$, we
obtain the exponent $b=2.591$, and its 95% confidence interval $[2.564,2.619]$
(this interval is computed in the standard way using Student’s
$t$-distribution, as described in [31]). These observations, the power law,
and its confidence interval are show on Figure 1. Thus, we can say that for
random $\mathcal{SVM}(n,2n,p)$-instances, and fixed $\epsilon=0.1$, the
complexity of Algorithm 1 scales as $O(n^{2.591})$. This represents a
polynomial improvement over general dense SOCP solvers with complexity
$O(n^{\omega+0.5})$. In practice, the polynomial speedup is conserved when
compared to ECOS [15], that has a measured running time scaling of
$O(n^{3.314})$, with a 95% confidence interval for the exponent of [3.297,
3.330] (this is consistent with the internal usage of a [36]-like matrix
multiplication algorithm, with complexity $O(n^{2.807})$). Neglecting constant
factors, this gives us a speedup of $10^{4}$ for $n=10^{6}$. The results from
the LIBSVM solver indicate that the training time with a linear kernel has a
complexity of $O(n^{3.112})$, with a 95% confidence interval for the exponent
of [2.799, 3.425]. These results suggest Algorithm 1 retains its advantage
even when compared to state-of-the-art specialized classical algorithms.

Additionally, we use the gathered data to verify that the accuracy of our
quantum (or approximate) SVM is close to the optimum: Figure 2 shows that both
at train- and at test-time the accuracies of all three classifiers are most of
the time within a few percent of each other, with Algorithm 1 often
outperforming the exact SOCP SVM classifier.

In conclusion, the performed numerical experiments indicate that Algorithm 1
provides a polynomial speedup for solving SOCPs with low- and medium precision
requirements. In particular, for SVM, we achieve a polynomial speedup with no
detriment to the quality of the trained classifier.

## Acknowledgmenets

This work was partly supported by IdEx Université de Paris ANR-18-IDEX-0001,
as well the French National Research Agency (ANR) projects QuBIC and QuDATA.

## References

  * Alizadeh and Goldfarb [2003] F. Alizadeh and D. Goldfarb.  Second-order cone programming.  _Math. Program._ , 95(1, Ser. B):3–51, 2003.  ISSN 0025-5610.  doi: 10.1007/s10107-002-0339-5.  ISMP 2000, Part 3 (Atlanta, GA).
  * Allcock and Hsieh [2020] Jonathan Allcock and Chang-Yu Hsieh.  A quantum extension of SVM-perf for training nonlinear SVMs in almost linear time.  _Quantum_ , 4:342, October 2020.  ISSN 2521-327X.  doi: 10.22331/q-2020-10-15-342.
  * Arodz and Saeedi [2019] Tomasz Arodz and Seyran Saeedi.  Quantum sparse support vector machines, 2019.  URL https://arxiv.org/abs/1902.01879.
  * Arora et al. [2012] Sanjeev Arora, Elad Hazan, and Satyen Kale.  The multiplicative weights update method: a meta-algorithm and applications.  _Theory Comput._ , 8:121–164, 2012.  doi: 10.4086/toc.2012.v008a006.
  * Ben-Tal and Nemirovski [2001] Aharon Ben-Tal and Arkadi Nemirovski.  _Lectures on modern convex optimization_.  MPS/SIAM Series on Optimization. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA; Mathematical Programming Society (MPS), Philadelphia, PA, 2001.  ISBN 0-89871-491-5.  doi: 10.1137/1.9780898718829.  Analysis, algorithms, and engineering applications.
  * Boyd and Vandenberghe [2004] Stephen Boyd and Lieven Vandenberghe.  _Convex optimization_.  Cambridge University Press, Cambridge, 2004.  ISBN 0-521-83378-7.  doi: 10.1017/CBO9780511804441.
  * Brandão et al. [2019] Fernando G. S. L. Brandão, Amir Kalev, Tongyang Li, Cedric Yen-Yu Lin, Krysta M. Svore, and Xiaodi Wu.  Quantum SDP solvers: large speed-ups, optimality, and applications to quantum learning.  In _46th International Colloquium on Automata, Languages, and Programming_ , volume 132 of _LIPIcs. Leibniz Int. Proc. Inform._ , pages Art. No. 27, 14. Schloss Dagstuhl. Leibniz-Zent. Inform., Wadern, 2019.  doi: 10.4230/LIPICS.ICALP.2019.27.
  * Brandão et al. [2020] Fernando G. S L. Brandão, Richard Kueng, and Daniel Stilck França.  Faster quantum and classical sdp approximations for quadratic binary optimization, 2020.  URL https://arxiv.org/abs/1909.04613.
  * Brandão and Svore [2017] Fernando G.S.L. Brandão and Krysta M. Svore.  Quantum speed-ups for solving semidefinite programs.  In _58th Annual IEEE Symposium on Foundations of Computer Science—FOCS 2017_ , pages 415–426. IEEE Computer Soc., Los Alamitos, CA, 2017.  doi: 10.1109/focs.2017.45.
  * Bubeck [2015] Sébastien Bubeck.  Convex optimization: Algorithms and complexity.  _Foundations and Trends® in Machine Learning_ , 8(3-4):231–357, 2015.  doi: 10.1561/2200000050.
  * Chakraborty et al. [2019] Shantanav Chakraborty, András Gilyén, and Stacey Jeffery.  The power of block-encoded matrix powers: improved regression techniques via faster Hamiltonian simulation.  In _46th International Colloquium on Automata, Languages, and Programming_ , volume 132 of _LIPIcs. Leibniz Int. Proc. Inform._ , pages Art. No. 33, 14. Schloss Dagstuhl. Leibniz-Zent. Inform., Wadern, 2019.  doi: 10.4230/LIPICS.ICALP.2019.33.
  * Chang and Lin [2011] Chih-Chung Chang and Chih-Jen Lin.  LIBSVM: A library for support vector machines.  _ACM Transactions on Intelligent Systems and Technology_ , 2(3):1–27, April 2011.  doi: 10.1145/1961189.1961199.
  * Cohen et al. [2019] Michael B. Cohen, Yin Tat Lee, and Zhao Song.  Solving linear programs in the current matrix multiplication time.  In _STOC’19—Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing_ , pages 938–942. ACM, New York, 2019.  doi: 10.1145/3313276.3316303.
  * Dollar [2005] Hilary Dollar.  _Iterative linear algebra for constrained optimization_.  PhD thesis, University of Oxford, 2005.  URL https://www.numerical.rl.ac.uk/people/hsd/thesismain.pdf.
  * Domahidi et al. [2013] Alexander Domahidi, Eric Chu, and Stephen Boyd.  ECOS: An SOCP solver for embedded systems.  In _2013 European Control Conference (ECC)_. IEEE, July 2013\.  doi: 10.23919/ecc.2013.6669541.
  * Fan et al. [2008] Rong-En Fan, Kai-Wei Chang, Cho-Jui Hsieh, Xiang-Rui Wang, and Chih-Jen Lin.  LIBLINEAR: A library for large linear classification.  _Journal of machine learning research_ , 9(Aug):1871–1874, 2008.  URL https://www.jmlr.org/papers/volume9/fan08a/fan08a.
  * Gilyén et al. [2019a] András Gilyén, Srinivasan Arunachalam, and Nathan Wiebe.  Optimizing quantum optimization algorithms via faster quantum gradient computation.  In _Proceedings of the Thirtieth Annual ACM-SIAM Symposium on Discrete Algorithms_ , pages 1425–1444. SIAM, Philadelphia, PA, 2019a.  doi: 10.1137/1.9781611975482.87.
  * Gilyén et al. [2019b] András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe.  Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics.  In _STOC’19—Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing_ , pages 193–204. ACM, New York, 2019b.  doi: 10.1145/3313276.3316366.
  * Gilyén et al. [2018] András Gilyén, Seth Lloyd, and Ewin Tang.  Quantum-inspired low-rank stochastic regression with logarithmic dependence on the dimension, 2018.  URL https://arxiv.org/abs/1811.04909.
  * Grover [1996] Lov K. Grover.  A fast quantum mechanical algorithm for database search.  In _Proceedings of the Twenty-eighth Annual ACM Symposium on the Theory of Computing (Philadelphia, PA, 1996)_ , pages 212–219. ACM, New York, 1996.  doi: 10.1145/237814.237866.
  * Harrow et al. [2009] Aram W. Harrow, Avinatan Hassidim, and Seth Lloyd.  Quantum algorithm for linear systems of equations.  _Physical Review Letters_ , 103(15), October 2009.  doi: 10.1103/physrevlett.103.150502.
  * Joachims [2006] Thorsten Joachims.  Training linear SVMs in linear time.  In _Proceedings of the 12th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’06, page 217–226, New York, NY, USA, 2006. Association for Computing Machinery.  ISBN 1595933395.  doi: 10.1145/1150402.1150429.
  * Karmarkar [1984] N. Karmarkar.  A new polynomial-time algorithm for linear programming.  In _Proceedings of the sixteenth annual ACM symposium on Theory of computing - STOC '84_ , pages 302–311. ACM Press, 1984.  doi: 10.1145/800057.808695.
  * Kerenidis and Prakash [2017] Iordanis Kerenidis and Anupam Prakash.  Quantum recommendation systems.  In Christos H. Papadimitriou, editor, _8th Innovations in Theoretical Computer Science Conference (ITCS 2017)_ , volume 67 of _Leibniz International Proceedings in Informatics (LIPIcs)_ , pages 49:1–49:21, Dagstuhl, Germany, 2017. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik.  ISBN 978-3-95977-029-3.  doi: 10.4230/LIPIcs.ITCS.2017.49.
  * Kerenidis and Prakash [2020a] Iordanis Kerenidis and Anupam Prakash.  Quantum gradient descent for linear systems and least squares.  _Phys. Rev. A_ , 101:022316, Feb 2020a.  doi: 10.1103/PhysRevA.101.022316.
  * Kerenidis and Prakash [2020b] Iordanis Kerenidis and Anupam Prakash.  A quantum interior point method for LPs and SDPs.  _ACM Transactions on Quantum Computing_ , 1(1), October 2020b.  ISSN 2643-6809.  doi: 10.1145/3406306.
  * Lee et al. [2019] Yin Tat Lee, Zhao Song, and Qiuyi Zhang.  Solving empirical risk minimization in the current matrix multiplication time.  In Alina Beygelzimer and Daniel Hsu, editors, _Proceedings of the Thirty-Second Conference on Learning Theory_ , volume 99 of _Proceedings of Machine Learning Research_ , pages 2140–2157, Phoenix, USA, 25–28 Jun 2019. PMLR.  URL http://proceedings.mlr.press/v99/lee19a.html.
  * Monteiro and Tsuchiya [2000] Renato D. C. Monteiro and Takashi Tsuchiya.  Polynomial convergence of primal-dual algorithms for the second-order cone program based on the MZ-family of directions.  _Math. Program._ , 88(1, Ser. A):61–83, 2000\.  ISSN 0025-5610.  doi: 10.1007/PL00011378.
  * Nesterov and Todd [1997] Yu. E. Nesterov and M. J. Todd.  Self-scaled barriers and interior-point methods for convex programming.  _Math. Oper. Res._ , 22(1):1–42, 1997.  ISSN 0364-765X.  doi: 10.1287/moor.22.1.1.
  * Nesterov and Todd [1998] Yu. E. Nesterov and M. J. Todd.  Primal-dual interior-point methods for self-scaled cones.  _SIAM J. Optim._ , 8(2):324–364, 1998.  ISSN 1052-6234.  doi: 10.1137/S1052623495290209.
  * Neter et al. [1996] John Neter, Michael H Kutner, Christopher J Nachtsheim, and William Wasserman.  _Applied linear statistical models_ , volume 4.  Irwin Chicago, 1996.
  * Nielsen and Chuang [2009] Michael A. Nielsen and Isaac L. Chuang.  _Quantum Computation and Quantum Information_.  Cambridge University Press, 2009.  doi: 10.1017/cbo9780511976667.
  * Rebentrost et al. [2014] Patrick Rebentrost, Masoud Mohseni, and Seth Lloyd.  Quantum support vector machine for big data classification.  _Physical Review Letters_ , 113(13), September 2014.  doi: 10.1103/physrevlett.113.130503.
  * Saad [2003] Yousef Saad.  _Iterative methods for sparse linear systems_.  Society for Industrial and Applied Mathematics, Philadelphia, PA, second edition, 2003.  ISBN 0-89871-534-2.  doi: 10.1137/1.9780898718003.
  * Shor [1994] Peter W. Shor.  Algorithms for quantum computation: discrete logarithms and factoring.  In _35th Annual Symposium on Foundations of Computer Science (Santa Fe, NM, 1994)_ , pages 124–134. IEEE Comput. Soc. Press, Los Alamitos, CA, 1994.  doi: 10.1109/SFCS.1994.365700.
  * Strassen [1969] Volker Strassen.  Gaussian elimination is not optimal.  _Numerische Mathematik_ , 13(4):354–356, August 1969.  doi: 10.1007/bf02165411.
  * Suykens and Vandewalle [1999] J.A.K. Suykens and J. Vandewalle.  Least squares support vector machine classifiers.  _Neural Processing Letters_ , 9(3):293–300, 1999\.  doi: 10.1023/a:1018628609742.
  * Suykens et al. [2002] J.A.K. Suykens, J. De Brabanter, L. Lukas, and J. Vandewalle.  Weighted least squares support vector machines: robustness and sparse approximation.  _Neurocomputing_ , 48(1-4):85–105, October 2002\.  doi: 10.1016/s0925-2312(01)00644-0.
  * Szilagyi et al. [2021] Daniel Szilagyi, Iordanis Kerenidis, and Anupam Prakash.  Quantum SVM via SOCP experiment logs.  Mar 2021.  doi: 10.6084/m9.figshare.11778189.v1.
  * van Apeldoorn and Gilyén [2019] Joran van Apeldoorn and András Gilyén.  Improvements in quantum SDP-solving with applications.  In _46th International Colloquium on Automata, Languages, and Programming_ , volume 132 of _LIPIcs. Leibniz Int. Proc. Inform._ , pages Art. No. 99, 15. Schloss Dagstuhl. Leibniz-Zent. Inform., Wadern, 2019.  doi: 10.4230/LIPICS.ICALP.2019.99.
  * van Apeldoorn et al. [2017] Joran van Apeldoorn, András Gilyén, Sander Gribling, and Ronald de Wolf.  Quantum SDP-solvers: Better upper and lower bounds.  In _2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS)_. IEEE, October 2017.  doi: 10.1109/focs.2017.44.
  * Ye et al. [1994] Yinyu Ye, Michael J. Todd, and Shinji Mizuno.  An $O(\sqrt{n}L)$-iteration homogeneous and self-dual linear programming algorithm.  _Math. Oper. Res._ , 19(1):53–67, 1994.  ISSN 0364-765X.  doi: 10.1287/moor.19.1.53.