# A quantum central path algorithm for linear optimization

Brandon Augustino Corresponding Author<EMAIL_ADDRESS>Sloan School of
Management, Massachusetts Institute of Technology Jiaqi Leng  Department of
Mathematics, University of Maryland, College Park Joint Center for Quantum
Information and Computer Science, University of Maryland Giacomo Nannicini
Department of Industrial and Systems Engineering, University of Southern
California Tamás Terlaky  Department of Industrial and Systems Engineering,
Lehigh University Quantum Computing and Optimization Lab, Lehigh University
Xiaodi Wu  Department of Computer Science, University of Maryland, College
Park Joint Center for Quantum Information and Computer Science, University of
Maryland

###### Abstract

We propose a novel quantum algorithm for solving linear optimization problems
by quantum-mechanical simulation of the central path. While interior point
methods follow the central path with an iterative algorithm that works with
successive linearizations of the perturbed KKT conditions, we perform a single
simulation working directly with the nonlinear complementarity equations.
Combining our approach with iterative refinement techniques, we obtain an
exact solution to a linear optimization problem involving $m$ constraints and
$n$ variables using at most ${\cal O}\left((m+n)\textup{nnz}(A)\kappa({\cal
M})L\cdot\textup{polylog}\left(m,n,\kappa({\cal M})\right)\right)$ elementary
gates and ${\cal O}\left(\textup{nnz}(A)L\right)$ classical arithmetic
operations, where $\textup{nnz}(A)$ is the total number of non-zero elements
found in the constraint matrix, $L$ denotes binary input length of the problem
data, and $\kappa({\cal M})$ is a condition number that depends only on the
problem data.

###### Contents

  1. 1 Introduction
    1. 1.1 Related work
    2. 1.2 Contributions
    3. 1.3 Technical summary and roadmap
  2. 2 Preliminaries
    1. 2.1 Interior Point Methods for Linear Optimization
      1. 2.1.1 Primal and dual LO problems
      2. 2.1.2 The self-dual embedding model
      3. 2.1.3 The central path and its neighborhood
    2. 2.2 Quantum adiabatic theorem for unbounded Hamiltonians
    3. 2.3 Quantum algorithms for Schrödinger equations
  3. 3 A quantum algorithm that traces the central path
    1. 3.1 Quantum representation of the central path
    2. 3.2 Quantum simulation of the central path
    3. 3.3 Quantum central path algorithm
    4. 3.4 Constructing an evaluation oracle for $f_{\mu}$
  4. 4 Recovering polylogarithmic dependence on precision
    1. 4.1 Exactly solving linear optimization problems
    2. 4.2 Iterative refinement for linear optimization
    3. 4.3 Comparison to existing LO solvers
  5. 5 Conclusion and outlook
  6. A Proof of Proposition 1
  7. B Additional Proofs

##  1 Introduction

Given a matrix $A\in\mathbb{R}^{m\times n}$ and vectors $b\in\mathbb{R}^{m}$
and $c\in\mathbb{R}^{n}$, we are interested in solving the following Linear
Optimization (LO) problem:

$\displaystyle\min_{x\in\mathbb{R}^{n}}$ $\displaystyle c^{\top}x$ (P) subject
to $\displaystyle Ax\geq b,$ $\displaystyle\leavevmode\nobreak\
\leavevmode\nobreak\ x\geq 0,$

where $x\in\mathbb{R}^{n}$ is the primal variable. The primal problem (P) has
an associated dual problem, given by

$\displaystyle\max_{y\in\mathbb{R}^{m}}$ $\displaystyle b^{\top}y$ (D) subject
to $\displaystyle A^{\top}y\leq c,$ $\displaystyle\quad\leavevmode\nobreak\
y\geq 0,$

where $A^{\top}\in\mathbb{R}^{n\times m}$ and $y\in\mathbb{R}^{m}$ is the dual
variable. Whenever the canonical primal and dual LO problems (P)-(D) are
feasible, strong duality holds: there exists a primal-dual optimal solution
$(x_{*},y_{*})$ to (P)-(D) with vanishing duality gap, i.e.,
$c^{\top}x_{*}-b^{\top}y_{*}=0$.

Linear Optimization is routinely used to model and solve fundamental problems
in economics, finance, and engineering. LO was brought to the forefront of
computer science and applied mathematics following the advent of digital
computer along with Dantzig’s Simplex method [Dan48, DOW+55]. The first
polynomial time algorithm for LO was Kachiyan’s ellipsoid method [Kha80],
though it failed to elicit an efficient practical implementation. Shortly
thereafter, Karmarkar introduced his projective method for LO problems in the
seminal work [Kar84], which improved the complexity of the ellipsoid method.
Although interior point algorithms had been studied since at least the 1950’s
[Dik67, Dik74, FM64, FM90, Fri54, Fri55, Fri56], Karmarkar’s algorithm was the
first to run in polynomial time, and the so-called Interior Point Method (IPM)
revolution was underway.

Enhancements to Karmarkar’s algorithm followed closely; Renegar [Ren88]
improved the iteration complexity, and Nesterov and Nemirovskii [NN88, NN94]
introduced the paradigm of self-concordant barrier functions. Today, the IPM
literature constitutes a vast and celebrated line of research, making it a
fruitless endeavor to provide a comprehensive review here. Rather, we point
the reader to the excellent texts [Ren01, RTV05, Ter13, Wri97, Ye11] and the
references therein. Many recent developments can be found in [vdB20, vdBLSS20,
CLS21, JKL+20, JSWZ20, LS14, LS15, Vla23], and the current state of the art
running time results are obtained by the randomized IPM of Cohen, Lee and Song
[CLS21] and the deterministic IPM from van den Brand [vdB20]. These IPMs can
solve LO problems to precision $\varepsilon\in(0,1)$ in $\widetilde{{\cal
O}}_{\frac{n}{\varepsilon}}((m+n)^{\omega})$ time111The authors in [vdB20,
CLS21] solve the standard form LO problem which requires the assumption
$\operatorname{rank}(A)=m$ (and in particular, $m\leq n$), so their complexity
result is typically written as $\widetilde{{\cal
O}}_{\frac{n}{\varepsilon}}(n^{\omega})$. We make the dependence on $m$
explicit because this assumption is unnecessary in the cannonical (or,
symmetric) form we consider: introducing $m$ slack variables transforms the
cannonical problem into standard form, and the augmented coefficient matrix is
of full row-rank. We also remark that the complexity of the algorithm in
[CLS21] is more accurately expressed as $\widetilde{{\cal
O}}_{\frac{n}{\varepsilon}}\left(n^{\omega}+n^{2.5-\frac{\alpha}{2}}+n^{2+\frac{1}{6}}\right)$,
where $\alpha$ is the dual-exponent of $\omega$. This simplifies to
$\widetilde{{\cal O}}_{\frac{n}{\varepsilon}}\left(n^{\omega}\right)$ for the
current value $\omega\approx 2.38$. Jiang et al. [JSWZ20] improved the
$n^{2+\frac{1}{6}}$ term to $n^{2+\frac{1}{18}}$., where $\omega\in[2,2.38)$
is the exponent of the running time for matrix multiplication (i.e., we can
multiply two $n\times n$ matrices in time ${\cal O}(n^{\omega}))$, and the
notation $\widetilde{{\cal O}}_{\alpha,\beta}(f(x))$ suppresses
polylogarithmic factors in $f(x)$, $\alpha$ and $\beta$ appearing in the
overall running time.

One can view IPMs as a homotopy approach for Newton’s method. IPMs are
initialized to a strictly feasible point
$x^{(0)}\in\\{x\in\mathbb{R}^{n}:Ax>b,\leavevmode\nobreak\ x>0\\}$, and
approximately track an analytic curve known as the central path towards
optimality. Specifically, letting $\mu>0$ and defining

$\phi:\\{x\in\mathbb{R}^{n}:Ax>b,\leavevmode\nobreak\
x>0\\}\mapsto\mathbb{R},$

a generic path-following scheme solves a sequence of barrier problems of the
form

$x(\mu)\coloneqq\operatorname*{arg\,min}_{\\{x\in\mathbb{R}^{n}:Ax>b,\leavevmode\nobreak\
x>0\\}}\quad c^{\top}x+\mu\phi(x),$ (1)

using Newton’s method. The barrier function $\phi$ is chosen to be self-
concordant: the value $\phi(x)$ diverges to $\infty$ upon approaching the
boundary of $\\{x\in\mathbb{R}^{n}:Ax>b,\leavevmode\nobreak\ x>0\\}$, and at a
high level, the norms of higher order derivatives of $\phi$ can be bounded in
terms of its Hessian.222To put it another way, the second-order Taylor series
approximation of (1) is highly accurate. The value of $\mu$ is decreased in
each iteration, and the optimal solution is reached upon tracing the central
path $\\{x(\mu):\mu>0\\}$ as $\mu\to 0$. The use of the self-concordant
barrier function $\phi$ in (1) ensures that we are always in the region of
rapid local convergence enjoyed by Newton’s method.

It is also standard in the literature to define the central path as the set of
minimizers associated with

$x(t)\coloneqq\operatorname*{arg\,min}_{\\{x\in\mathbb{R}^{n}:Ax>b,\leavevmode\nobreak\
x>0\\}}\quad tc^{\top}x+\phi(x).$ (2)

The value of $t$ is increased in each iteration, and so tracing the central
path $\\{x(t):t>0\\}$, we approach an optimal solution to (P) as $t\to\infty$.
Note however, there is no meaningful difference between (1) and (2). Indeed,
taking $\mu=1/t$ and applying Newton’s method to (1) yields the same sequence
of minimizers as one would obtain from applying Newton’s method to (2). Yet,
the perspective offered by (1) can be attractive for both theoretical and
practical reasons; tracking $t\to\infty$ in (2) would require computation
involving huge numbers.

In every iteration of the classical IPM, one must obtain the solution to a
linearized set of perturbed Karush-Kuhn-Tucker (KKT) optimality conditions
known as the Newton linear system. The computation of this so-called Newton
step is usually performed by solving dense linear systems of equations on a
very large scale. The IPMs that achieve the best complexity results [vdB20,
CLS21] perform this step exactly in the first iteration, and subsequently
utilize sophisticated data structures to maintain an approximation of the
Hessian inverse, amortizing the cost of solving the linear systems over the
run of the algorithm. To the best of our knowledge, the improved theoretical
guarantees of the IPMs in [vdB20, CLS21] have not led to a practical
implementation at the time of writing. Additionally, while the previous decade
saw improvements in the iteration complexity for the specific applications of
maximum flow [M1̨3, M1̨6], minimum cost flow [CMSV17], matrix scaling
[CMTV17], and $\ell_{p}$-regression [BCLL18], whether the asymptotic iteration
count of ${\cal O}(\sqrt{n}\log(1/\varepsilon))$ (or, ${\cal
O}(\sqrt{m}\log(1/\varepsilon))$ if $m\ll n$) can be reduced in the general
setting remains an important open questions in optimization. Quantum IPMs
(QIPMs) currently found in the literature seek to accelerate the solution of
the Newton linear system via quantum linear systems algorithms (QLSAs) [CW12,
CGJ19, HHL09], hence they do not reduce the iteration count, and it is unclear
whether or not this approach can provide an overall speedup. We review these
works in detail next.

###  1.1 Related work

QIPMs were first introduced by Kerenidis and Prakash [KP20], who proposed a
quantum algorithm for solving LO and Semidefinite Optimization (SDO) problems.
The main idea of their approach is to solve the Newton linear system at each
iterate using a QLSA, and obtain a classical estimate of the resulting quantum
state via quantum state tomography. This approach limits one to inexactly
solving the Newton system at each iterate, and therefore, additional
safeguards need to be taken in order to guarantee convergence. To reconcile
the quantum noise introduced into the Newton steps, Augustino et al. [ANTZ23b]
proposed two convergent QIPMs for SDO and LO. The first closely quantized the
classical Inexact-Infeasible IPM of [TK02]. The second framework is a novel
Inexact-Feasible QIPM (IF-QIPM) that uses a nullspace representation of the
Newton system to ensure the sequence of iterates maintain primal-dual
feasibility, in spite of using an inexact linear system subroutine. These
ideas were specialized to LO by Mohammadisiahroudi et al. [MFT22, MFWT23], who
use the iterative refinement algorithm for LO from Gleixner et al. [GSW12,
GSW16, GS20] to exponentially improve the dependence on $\varepsilon^{-1}$,
the inverse precision to which we seek to solve (P)-(D).

While quantizing IPMs in this manner has led to polynomial speedups in the
problem dimension $n$, the current approach has two major drawbacks: (i)
existing QIPMs assume a strong input model (i.e., the quantum RAM model) that
is not always justified in a practical setting, and speedups are heavily
dependent on this model; (ii) the iterative uses of QLSAs rely on an expensive
tomography procedure, prohibiting a conclusive overall speedup. While it is
possible that there is a gap between the algorithm’s theoretical worst-case
running time and its practical performance, the latter cannot yet be studied:
the data structures and computational primitives utilized by existing QIPMs
are beyond the capabilities of near-term quantum devices. A detailed resource
analysis by Dalzell et al. [DCS+22] for the specific application of portfolio
optimization found that a QIPM based on the combined used of QLSA and
tomography would require a $T$-gate count of $10^{28}$ when $n=100$.
Informally, this implies that it would take a quantum computer millions of
years to solve a problem that could be solved in seconds on a personal laptop.
Thus, existing QIPMs do not provide a convincing end-to-end speedup over
classical IPMs, and their estimated resource cost is a negative result in
relation to their viability even on large-scale devices.

These issues aside, the use of quantum subroutines in QIPMs only serves as an
attempt to accelerate a single step of the classical IPM. An inherent drawback
of this approach is that QIPMs produce a sequence of iterates that exactly
mimics the classical IPM trajectory, and therefore lack the potential to
exploit quantum effects or improve on the worst-case iteration complexity
bound. This raises the question as to whether there exists a “more naturally
quantum” algorithm that solves optimization problems by tracking the central
path.

In the recent work [LHLW23], Leng, Hickman, Li and Wu proposed a truly quantum
analogue to classical (accelerated) gradient descent named Quantum Hamiltonian
Descent (QHD). Their work can be viewed as a quantization of the Bregman-
Lagrangian framework for analyzing the continuous-time dynamics of gradient
descent introduced by Wibisono, Wilson, and Jordan [WWJ16]. The QHD framework
admits an implementation on current quantum devices; this is a (practical)
advantage over existing QIPMs, because the quantum primitives underlying QLSAs
and state tomography are out of reach for current quantum computers. The
underlying foundation of QHD is to first recast the task of optimization as a
dynamical system, and then quantize the resulting continuous-time dynamics.
This novel approach offers a new paradigm for the design of quantum algorithms
for optimization, which previously, amounted to using quantum subroutines to
accelerate a single step of a classical optimization algorithm. Interestingly,
the study of dynamical systems arising from IPMs was an active area of
research following the introduction of Karmarkar’s algorithm in 1984. Many
works [Ans88, BL89a, BL89b, Fay91, Fay95, Meg89] studied the trajectories
generated by the so-called Newton barrier flow. Megiddo [Meg89] and Bayer and
Lagarias [BL89a, BL89b] provided (classical) Lagrangian and Hamiltonian
dynamical systems for path-following methods that use a logarithmic barrier
function. Faybusovich later generalized these ideas to convex optimization in
[Fay95]. Yet, the design and complexity analysis of an algorithm based on
these ideas was not considered in these works: although these papers precisely
characterize the central path, these efforts did not translate into new
algorithms.

###  1.2 Contributions

We use a previously unexplored connection between the central path and the
Schrödinger equation to develop a new quantum algorithm for LO, that we call
the Quantum Central Path Method (QCPM). This is achieved by proposing a
1-parameter family of Hamiltonian operators over the positive orthant that
encodes the behavior of the central path. Specifically, we show that one can
approximately follow the central path by simulating a Shrödinger equation
associated to a certain Hamiltonian. The ground state of this Hamiltonian
encodes a squeezed Gaussian distribution centered at an $\varepsilon$-optimal
solution to the primal dual LO pair (P)-(D). The QCPM inherits many desirable
properties from IPM theory (such as provable convergence to the optimal
solution), but we emphasize that the QCPM, unlike (Q)IPMs, is not an iterative
algorithm: the desired solution can be obtained in _one shot_ by measuring the
final state. This could potentially provide a significant advantage over the
state of the art (Q)IPMs, which require ${\cal
O}\left(\min\left\\{\sqrt{m},\sqrt{n}\right\\}\log\left(\frac{1}{\varepsilon}\right)\right)$
iterations to reach an $\varepsilon$-optimal solution. Of course simulating
the Schrödinger equation does not come for free, and this determines the
complexity of our algorithm. For this simulation task we use a slightly
improved version of a result in [CLL+22a], tailored to our specific setup:
this provides a fast quantum algorithm that has near-optimal dependence (with
respect to the simulation task) on the dimension of the problem and the
precision parameter.

The idea for the QCPM is rooted in the early works on the Newton barrier flow
discussed above, with some significant differences: we use the _self-dual
embedding_ model [YTM94] to guarantee an easy-to-prepare ground state for an
initial Hamiltonian, and this is crucial to the eventual application of the
adiabatic theorem to show convergence to the optimal solution (i.e., ground
state of the final Hamiltonian). We also work _directly_ with the nonlinear
complementarity equations in the perturbed KKT conditions that define the
central path, rather than using a local linearization as is done in IPMs. To
the best of our knowledge, the QCPM is the first algorithm that takes this
perspective.

By combining the QCPM with the iterative refinement algorithm introduced in
[GSW12, GSW16, GS20], we obtain an algorithm that can identify an exact
solution of a linear optimization problem involving $m$ constraints and $n$
variables using

$\widetilde{{\cal O}}_{m,n,\kappa({\cal
M}),\frac{1}{\delta}}\left((m+n)\textup{nnz}(A)\kappa({\cal M})L\right)$

elementary gates and ${\cal O}\left(\textup{nnz}(A)L\right)$ classical
operations, where $\textup{nnz}(A)$ is the total number of non-zero entries
found in $A$, $L$ denotes binary input length of the problem data $(A,b,c)$,
$\kappa({\cal M})$ is a condition number that depends only on the problem
data333This contrasts with the $\kappa$ found in the complexity analysis of
QIPMs: there, $\kappa$ is the condition number of the Newton linear system,
which changes at every iteration and in fact $\kappa\to\infty$ as we approach
the set of optimal solutions., and $1-\delta$ is the probability of success.
Note that by exact, we mean a solution obtained to precision $\varepsilon\leq
2^{-L}$. We also stress that our algorithm does not utilize QRAM, and the
stated complexity result only requires access to the sparse binary
representation of $(A,b,c)$.

To contextualize our results, we provide a comparison of our algorithm to the
current state of the art algorithms for solving Linear Optimization problems
in both the classical and quantum models of computation in Table 1. One can
observe that the QCPM outperforms the state of the art classical and quantum
IPMs with respect to $m$ and $n$ when $\textup{nnz}(A)<(m+n)^{\omega-1}$, and
achieves an end-to-end speedup over these algorithms whenever
$\textup{nnz}(A)\kappa({\cal M})<(m+n)^{\omega-1}\cdot\textup{polylog}(m,n)$.
Strikingly, these speedups are not reliant on access to a classical-
write/quantum-read RAM (QRAM): our algorithm only requires access to the
natural binary description for the LO problem data $(A,b,c)$. Like the QIPM in
[MFWT23], our algorithm exhibits a dependence on a condition number that
depends only on the problem input data, but our condition number dependence is
linear, while it is quadratic in [MFWT23]. We also emphasize that the QIPM
from Mohammadisiahroudi et al. [MFWT23] assumes access to an exponentially
large (in the problem dimension $n$) QRAM, for tomography. Our algorithm is
outperformed in $m$ and $n$ by quantum algorithms for zero-sum games [vAG19b,
BGJ+23]: these algorithms achieve sublinear running times in these quantities,
but depend polynomially on the inverse precision (our algorithm has
polylogarithmic dependence) and assume access to QRAM. The algorithms found in
[vAG19b, BGJ+23] are also sensitive to scale: they rely on the problem data
being normalized and their running times scale polynomially in the norm of the
solution. We point out that algorithms for zero-sum games employ different
definitions of optimality than (Q)IPMs and our QCPM, and thus, we refer the
reader to Section 4.3 for a more detailed comparison.

Table 1: Complexity to solve the primal-dual pair (P)-(D) to precision $\varepsilon$ Algorithm | Complexity | QRAM | Notes
---|---|---|---
IPM [CLS21, vdB20] | $\widetilde{\mathcal{O}}_{n,\frac{1}{\varepsilon}}((m+n)^{\omega})$ | - |
QMMWU [BGJ+23] | $\widetilde{{\cal O}}\left(\sqrt{m+n}r^{2.5}\varepsilon^{-2.5}+\varepsilon^{-3}\right)$ | ✓ | $r\geq\|x_{*}\|_{1}$
IR-QIPM [MFWT23] | $\widetilde{\mathcal{O}}_{n,\kappa(Q),\frac{1}{\varepsilon}}\left((m+n)^{2.5}\kappa(Q)^{2}\|Q\|\|x_{*}\|^{5}\right)$ | ✓ |
IR-QCPM (this work) | $\widetilde{\mathcal{O}}_{m,n,\kappa({\cal M}),\frac{1}{\varepsilon}}\left(\left(m+n\right)\textrm{nnz}(A)\kappa(\mathcal{M})\right)$ | ✗ |

The QCPM proposed in this paper provides (for the first time, as far as we are
aware) convincing evidence that an end-to-end quantum speedup for LO is
possible even in the circuit model without QRAM. Previous works on quantum
algorithms for LO or SDO rely on QRAM, and are slower than the best classical
algorithms in at least some problem parameters, such as precision or solution
size [vAG19a, vAG19b, vAGGdW20, BS17, BGJ+23]. An exception to this is
[ANTZ23a], which gives a quantum algorithm for the semidefinite relaxation of
MaxCut that is faster than the best known classical algorithm, but the speedup
in [ANTZ23a] vanishes without QRAM. We view the main contribution of our paper
as conceptual, rather than technical: the technical ingredients for our
algorithm and its analysis were largely present in the literature, but we
believe that the QCPM idea could provide a new, potentially interesting avenue
for constrained convex optimization on quantum computers, and the usefulness
of this idea is already evidenced by the end-to-end quantum speedup that we
obtain for some LO problems. It may be useful to perform a detailed resource
analysis of the QCPM framework, in the same spirit as the analysis of Dalzell
et al. [DCS+22] for QIPMs; we leave this for future work.

Next, we provide a more detailed summary of the technical overview of the
paper’s main results.

###  1.3 Technical summary and roadmap

Hamiltonian simulation is a fundamental task in quantum computation as it
separates the computational power between quantum and classical [Osb12]. For a
given time-dependent Hamiltonian operator ${\cal H}(t)$ and an initial state
$\left|{\psi_{0}}\right\rangle$, the Hamiltonian simulation task is to
simulate the dynamics governed by the Schrödinger equation:

$\mathbf{i}\frac{\text{d}}{\text{d}t}\left|{\psi(t)}\right\rangle={\cal
H}(t)\left|{\psi(t)}\right\rangle,\text{ subject to
}\left|{\psi(0)}\right\rangle=\left|{\psi_{0}}\right\rangle.$

In this paper, we propose a Quantum Central Path Method (QCPM) relying on the
efficient simulation of Schrödinger dynamics on quantum computers (see Theorem
3). Specifically, we design a quantum Hamiltonian (known as the central-path
Hamiltonian) that _encodes_ the behavior of the central path in its ground
state. By preparing an initial state centered at a point on the central path
and simulating the associated Schrödinger equation for a sufficiently long
evolution time $T$, the quantum adiabatic theorem (Theorem 2) guarantees that
the quantum register remains in the instantaneous ground state of the central-
path Hamiltonian. At the end of the Hamiltonian simulation, we obtain an
(approximate) optimal solution to the primal-dual LO pair (P)-(D) by simply
measuring the final state. Our QCPM shares some features with the Quantum
Hamiltonian Descent (QHD) algorithm [LHLW23]: (i) both can be regarded as
_quantizations_ of classical dynamical processes (path-following in QCPM,
accelerated gradient descent in QHD), (ii) both are formulated as quantum
dynamics governed by Schrödinger equations, and (iii) in both cases the ground
state of the Hamiltonian operator encodes the solution to an optimization
problem. The central-path Hamiltonian can be interpreted as a _global_
linearization of the central path: this is essentially different from any
existing QIPMs, where each Newton step represents a _local_ linearization of
the central path. Consequently, our algorithm is more efficient because there
is no need for iterative uses of QLSA and state tomography.

The adiabatic theorem is predicated on the assumption that we start in a
ground state of the initial Hamiltonian ${\cal H}(0)$. In general, this ground
state is hard to prepare because it would require the exact location of the
starting point on the central path for any LO problem, a problem known as hard
as the linear optimization itself. We overcome this challenge by casting the
primal-dual pair (P)-(D) as a slightly larger problem (of dimension $m+n+2$)
known as the self-dual embedding model (first proposed by Ye, Todd and Mizuno
[YTM94]), a standard approach employed by classical IPMs (both in theory and
in practice). This allows us to make a step towards providing a practical
algorithm because the self-dual embedding model always admits the all-ones
vector as an interior-feasible solution. By using the self-dual embedding
model as a foundation for the design of our central-path Hamiltonian, we
ensure that the ground state of the initial Hamiltonian is trivial to prepare.
A review of the self-dual embedding model and some important lemmas are
provided in Section 2.1.2.

Under the self-dual embedding framework, the central path is (uniquely)
characterized by element-wise complementarity and positivity constraints on
the variable $z$ and its associated slack $s(z)$:

${\cal
C}\coloneqq\left\\{z(\mu):\mu\in(0,1]\right\\}=\left\\{z\in\mathbb{R}^{m+n+2}:z_{i}\cdot
s(z)_{i}=\mu,\leavevmode\nobreak\ z_{i},s(z)_{i}>0,\leavevmode\nobreak\
\forall i\in 1,\dots,m+n+2\right\\}.$

Classical IPMs track the central path by iteratively applying Newton’s method,
which requires locally linearizing the nonlinear complementarity equation

$z_{i}\cdot s(z)_{i}=\mu\quad\forall i\in 1,\dots,m+n+2,$

and successively decreasing $\mu$ by a constant factor after each Newton step.
We propose a different approach: a time-dependent Schrödinger operator on the
positive orthant that directly incorporates the nonlinear complementarity
equation as a potential field,

${\cal H}(\mu(t))\coloneqq-\frac{h(t)}{2}\nabla^{2}+\left\|z\odot
s(z)-\mu(t)e\right\|^{2},$

where for a fixed $t>0$, the global minimum of the potential function
$\left\|z\odot s(z)-\mu(t)e\right\|^{2}$ corresponds to a point on the central
path. Here, $h(t)$ and $\mu(t)$ are monotonically decreasing functions in $t$,
$e$ is the all-ones vector of dimension $m+n+2$ and $\odot$ denotes the
Hadamard (or, element-wise) product. The central-path Hamiltonian ${\cal
H}(\mu(t))$ is inspired by the canonical time-independent Hamiltonian in
quantum mechanics:

${\cal H}\coloneqq-\nabla^{2}+V,$

where $\nabla^{2}$ denotes the Laplacian operator in $\mathbb{R}^{n}$,
accounting for kinetic energy, and $V\colon\mathbb{R}^{n}\mapsto\mathbb{R}$ is
a potential function, accounting for potential energy. In our time-dependent
Hamiltonian ${\cal H}(\mu(t))$, the role of $h(t)$ is to control the variance
of the probability distribution encoded in the ground state of the central-
path Hamiltonian. In the usual setting, a (small) fixed value $\hbar$ known as
Planck’s constant is used, but this choice is insufficient to guarantee
convergence in our setting: by decreasing $h(t)$ over the course of the time
evolution, we ensure that we concentrate on a probability mass centered at the
optimal solution. The detailed construction of the central-path Hamiltonian is
provided in Section 3.1. Properly designing $h(t)$ is a crucial technical
insight of our approach, and $\mu(t)$ is constructed using ideas introduced in
[AL22].

Equipped with the Hamiltonian ${\cal H}(\mu(t))$ defined above, the remaining
challenge is two-fold: we need to (i) verify that simulating the Schrödinger
equation associated to our central-path Hamiltonian does in fact solve our
optimization problem; and (ii) rigorously analyze the cost of doing so. We
establish the correctness of our approach through a combination of results
already found in the literature and new results of our own. In Section 3.1, we
apply the harmonic approximation theory of Schrödinger operators in order to
certify that, for fixed $\mu\in(0,1]$, the ground state of the central-path
Hamiltonian gives rise to a Gaussian distribution centered at the
corresponding point on the central path $z(\mu)$. This suggests another
interesting connection to classical IPM theory: the harmonic approximation of
the low-energy spectrum is obtained via a second-order Taylor series
approximation of ${\cal H}(\mu)$ at $z(\mu)$, which is precisely how one
applies Newton’s method to barrier problems of the form (1) in the classical
IPM framework. We also provide a lower bound on the minimum spectral gap: this
is an important component of the complexity of simulating our Schrödinger
equation, and we derive the lower bound by again interpreting the system
governed by ${\cal H}(\mu(t))$ as a quantum harmonic oscillator, whose
spectral gap is well understood. In Section 3.2, we demonstrate that properly
defining $h(t)$ and $\mu(t)$ allows one to solve (P)-(D) via a single
simulation of our Schrödinger equation: the final output is a classical
description of an $\varepsilon$-optimal solution, and no intermediate
measurements are required. This is in contrast with existing quantum
algorithms for LO [vAG19b, ANTZ23b, BGJ+23, KP20, MFWT23], which require
either ${\cal O}\left(\sqrt{n}\log\frac{1}{\varepsilon}\right)$ or ${\cal
O}\left(\frac{\log(n)}{\varepsilon^{2}}\right)$ calls to a Hamiltonian
simulation subroutine (to either solve quantum linear systems or prepare Gibbs
states).

We analyze the complexity of the QCPM in Section 3.3. To achieve our
complexity result, we combine the requirements for our convergence guarantees
with the cost of the simulation algorithm from [CLL+22a]. We prove that the
QCPM returns a classical description of an $\varepsilon$-precise solution to
(P)-(D) using at most

$\widetilde{{\cal O}}_{m,n,\kappa({\cal
M}),\frac{1}{\varepsilon},\frac{1}{\delta}}\left((m+n)\textup{nnz}(A)\kappa({\cal
M})\varepsilon^{-1}\right)$

elementary gates. This implies a polynomial speedup in $m$ and $n$ over the
classical state-of-the-art. We then exponentially improve the dependence on
precision using iterative refining techniques from [GSW12, GSW16, GS20, MFT22,
MFWT23] in Section 4. The basic idea is to treat the QCPM as a low-precision
oracle for LO and solve a sequence of related LO problems to constant
precision, rather than solve one problem to extended precision. This improves
the dependence on precision from $\varepsilon^{-1}$ to
$\log\left(\varepsilon^{-1}\right)$, at the cost of constant overhead and
classical matrix-vector products. The overall complexity of the resulting
scheme suggests that the regime for end-to-end quantum speedup is largely
determined by the condition number dependence. We stress that the speedups
achieved here are independent of QRAM. It may be possible to improve the
condition number dependence by analyzing the algorithm in terms of the local
Riemannian geometry induced by the Hessian of our potential function, similar
to the classical IPM analysis. We also believe our techniques should readily
generalize to other classes of constrained convex optimization problems.

The rest of this paper is organized in the following manner. In Section 2 we
define some notation and review the theory of classical interior point methods
applied to solving LO problems, before providing results on the adiabatic
theorem for unbounded Hamiltonians and algorithms for simulating Schrödinger
equations. Section 3 concerns the design and analysis of the Quantum Central
Path Method, and Section 4 provides some technical improvements and compares
our overall complexity result to the current state-of-the-art. Section 5
concludes the paper.

##  2 Preliminaries

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. 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. For a vector $x\in\mathbb{R}^{n}$, the matrix
$\operatorname{diag}(x)\in\mathbb{R}^{n\times n}$ takes the values of $x$
along its diagonal and zero elsewhere. For vectors $u,v\in\mathbb{R}^{n}$,
$u\odot v$ denotes their Hadamard (or, entry-wise) product:

$u\odot v=\begin{bmatrix}u_{1}\cdot v_{1}\\\ \vdots\\\ u_{n}\cdot
v_{n}\end{bmatrix}.$

We write $0_{n}\in\mathbb{R}^{n}$ when referring to an all-zeros vector of
length $n$, and $0_{m\times n}\in\mathbb{R}^{m\times n}$ denotes the $m\times
n$ all-zeros matrix.

The smallest and largest singular values of a matrix $A$ are denoted
$\sigma_{\min}(A),\sigma_{\max}(A)$, and the smallest and largest eigenvalues
are denoted $\lambda_{\min}(A),\lambda_{\max}(A)$. The operator norm of $A$ is
defined as $\|A\|\coloneqq\sigma_{\max}(A)$. We denote the condition number of
$A$ by $\kappa(A)=\|A\|\cdot\|A^{+}\|$, where $A^{+}$ is the pseudo-inverse of
the matrix $A$.

We let ${\cal S}_{+}^{n}$ and ${\cal S}_{++}^{n}$ represent the spaces of
symmetric positive semidefinite, and symmetric positive definite matrices in
$\mathbb{R}^{n\times n}$, respectively. For $U,V\in{\cal S}^{n}$, we write
$U\succeq V$ ($U\succ V$) to indicate that the matrix $U-V$ is symmetric
positive semidefinite (symmetric positive definite), i.e., $U-V\in{\cal
S}^{n}_{+}$ ($U-V\in{\cal S}^{n}_{++}$).

### Order estimates

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

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

We write $f(x)=\Omega(g(x))\iff g(x)={\cal O}(f(x))$. If there exists positive
constants $\alpha_{1}$ and $\alpha_{2}$ such that

$\alpha_{1}g(x)\leq f(x)\leq\alpha_{2}g(x)\quad\forall x>0,$

then we write $f(x)=\Theta(g(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}}_{\alpha,\beta}\leavevmode\nobreak\ (f(x))={\cal
O}(f(x)\cdot\textup{polylog}(\alpha,\beta,f(x))).$

###  2.1 Interior Point Methods for Linear Optimization

In this section we outline the classical IPM theory applied to solving linear
optimization problems.

####  2.1.1 Primal and dual LO problems

Recall that we are interested in the primal LO problem and its dual

$\min_{x\in{\cal P}}\quad c^{\top}x,\quad\quad\max_{y\in{\cal D}}\quad
b^{\top}y,$

where ${\cal P}$ and ${\cal D}$ are the primal and dual feasible sets, defined
as

${\cal P}\coloneqq\left\\{x\in\mathbb{R}^{n}:Ax\geq b,\leavevmode\nobreak\
x\geq 0\right\\},\quad{\cal
D}\coloneqq\left\\{y\in\mathbb{R}^{m}:A^{\top}y\leq c,\leavevmode\nobreak\
y\geq 0\right\\}.$

Likewise, the sets of interior feasible solutions of (P) and (D) are given by

$\operatorname{int}\left({\cal
P}\right)\coloneqq\left\\{x\in\mathbb{R}^{n}:Ax>b,\leavevmode\nobreak\
x>0\right\\},\quad\operatorname{int}\left({\cal
D}\right)\coloneqq\left\\{y\in\mathbb{R}^{m}:A^{\top}y<c,\leavevmode\nobreak\
y>0\right\\}.$

The well known property of weak duality always holds, and asserts that any
primal feasible $x\in{\cal P}$ provides an upper bound $c^{\top}x$ on the
value $b^{\top}y$, and conversely, any dual feasible $y\in{\cal D}$ provides a
lower bound $b^{\top}y$ on the value $c^{\top}x$. The nonnegative quantity

$c^{\top}x-b^{\top}y$

is referred to as the duality gap associated to the pair $(x,y)\in{\cal
P}\times{\cal D}$.

Whenever $(x_{*},y_{*})\in{\cal P}\times{\cal D}$ exhibit vanishing duality
gap, i.e., $c^{\top}x_{*}-b^{\top}y_{*}=0$, then $x_{*}$ is an optimal
solution to (P), and $y_{*}$ is an optimal solution to (D). Hence, we define
the optimal set of (P)-(D) to be

${\cal P}{\cal D}_{*}\coloneqq\left\\{(x,y)\in{\cal P}\times{\cal
D}:c^{\top}x=b^{\top}y\right\\}.$

Clearly, ${\cal P}{\cal D}_{*}$ is nonempty if and only if the inequality
system

$\displaystyle Ax$ $\displaystyle\geq b,\quad$ $\displaystyle x\geq 0,$ (3)
$\displaystyle-A^{\top}y$ $\displaystyle\geq-c,\quad$ $\displaystyle y\geq 0,$
$\displaystyle b^{\top}y-c^{\top}x$ $\displaystyle\geq 0,$

is solvable. Introducing a homogenizing variable $\beta$, system (3) is
equivalent to the homogeneous system

$\begin{bmatrix}0_{m\times m}&A&-b\\\ -A^{\top}&0_{n\times n}&c\\\
b^{\top}&-c^{\top}&0\end{bmatrix}\begin{bmatrix}y\\\ x\\\
\beta\end{bmatrix}\geq\begin{bmatrix}0_{m}\\\ 0_{n}\\\ 0\end{bmatrix},\quad
x\geq 0,\leavevmode\nobreak\ y\geq 0,\leavevmode\nobreak\ \beta\geq 0,$ (4)

when $\beta=1$. Indeed, one can verify that any $(x,y,\beta)$ which solves (4)
with $\beta>0$ gives rise to a solution
$\left(\frac{x}{\beta},\frac{y}{\beta},1\right)$ to (3). More concisely,
letting

${\overline{M\mkern 4.0mu}\mkern-4.0mu}{}\coloneqq\begin{bmatrix}0_{m\times
m}&A&-b\\\ -A^{\top}&0_{n\times n}&c\\\
b^{\top}&-c^{\top}&0\end{bmatrix},\quad\mkern
1.5mu\overline{\mkern-1.5muz}\coloneqq\begin{bmatrix}y\\\ x\\\
\beta\end{bmatrix},$

solving (P) and (D) to optimality is equivalent [RTV05, see, Theorem I.3] to
obtaining a feasible solution to the system

${\overline{M\mkern 4.0mu}\mkern-4.0mu}{}\mkern
1.5mu\overline{\mkern-1.5muz}\geq 0,\quad\mkern
1.5mu\overline{\mkern-1.5muz}\geq 0,\quad\beta>0.$ (5)

Next, we discuss the condition under which (5) is solvable, which is
intimately related to the practical concern of starting from an interior
feasible solution.

####  2.1.2 The self-dual embedding model

We begin with a formal definition of the interior point condition (IPC) from
[RTV05].

######  Definition 1 (Definition I.4 in [RTV05]).

We say that any system of (linear) equalities and (linear) inequalities
satisfies the interior-point condition (IPC) if there exists a solution that
satisfies all inequality constraints in the system.

The IPC ensures that the primal-dual pair (P) and (D) has an optimal solution
$(x_{*},y_{*})\in{\cal P}{\cal D}_{*}$ whenever there exists a strictly
feasible solution $(x,y)\in\operatorname{int}({\cal
P})\times\operatorname{int}({\cal D})$. Being able to easily determine an
interior point $(x,y)\in\operatorname{int}({\cal
P})\times\operatorname{int}({\cal D})$ is also a practical concern: (feasible)
IPMs are initialized to a strictly feasible starting point, and naïvely
determining a strictly feasible solution to (P) and (D) is as challenging as
solving these problems to optimality. Here we discuss a way to embed the
problems (P) and (D) into a slightly larger one that always has a trivial
strictly feasible solution, and can be readily solved with an IPM, following a
strategy first introduced by Ye, Todd and Mizuno [YTM94].

While system (5) does not satisfy the IPC, this can be reconciled upon
introducing another auxiliary variable $\vartheta\geq 0$, and adding one
additional row and column to ${\overline{M\mkern 4.0mu}\mkern-4.0mu}{}$ as
follows:

$M\coloneqq\begin{bmatrix}{\overline{M\mkern 4.0mu}\mkern-4.0mu}{}&r\\\
-r^{\top}&0\end{bmatrix}=\begin{bmatrix}\begin{bmatrix}0_{m\times m}&A&-b\\\
-A^{\top}&0_{n\times n}&c\\\ b^{\top}&-c^{\top}&0\end{bmatrix}&r\\\
-r^{\top}&0\end{bmatrix},\quad z\coloneqq\begin{bmatrix}\mkern
1.5mu\overline{\mkern-1.5muz}\\\ \vartheta\end{bmatrix}=\begin{bmatrix}y\\\
x\\\ \beta\\\ \vartheta\end{bmatrix},$

where

$r=\bar{e}-{\overline{M\mkern 4.0mu}\mkern-4.0mu}{}\bar{e},$

and $\bar{e}$ is the all-ones vector of length $m+n+1$. Note that $M$ is a
skew-symmetric matrix of dimension $\mkern
1.5mu\overline{\mkern-1.5mun}\coloneqq n+m+2$, i.e., $M^{\top}=-M$. Defining
$q\in\mathbb{R}^{\mkern 1.5mu\overline{\mkern-1.5mun}}$ to be the vector

$q=\begin{bmatrix}0_{\mkern 1.5mu\overline{\mkern-1.5mun}-1}\\\ \mkern
1.5mu\overline{\mkern-1.5mun}\end{bmatrix},$

one can also see that the system

$Mz\geq-q,\quad z\geq 0,$ (6)

admits the all-ones vector of length $\mkern 1.5mu\overline{\mkern-1.5mun}$ as
a strictly feasible solution, which we denote by $e$. If $z$ solves (6) with
$\vartheta=0$, i.e., we have $z=(\mkern 1.5mu\overline{\mkern-1.5muz},0)$,
then $\mkern 1.5mu\overline{\mkern-1.5muz}$ must be a solution to (5). Hence,
solving (P) and (D) can be reduced to finding a solution to (6) with
$\vartheta=0$ and $\beta>0$.

The self-dual embedding of (P) and (D) is defined as

$\min\left\\{q^{\top}z:Mz\geq-q,\leavevmode\nobreak\ z\geq 0\right\\}.$ (SP)

This problem is called “self-dual” because its dual problem

$\max\left\\{-q^{\top}u:M^{\top}u\leq q,\leavevmode\nobreak\ u\geq
0\right\\},$

can be recognized as equivalent to (SP), upon recalling that $M$ is a skew-
symmetric matrix with $M^{\top}=-M$. Like system (6), the self-dual embedding
problem (SP) also admits the strictly feasible solution $e$, and thus
trivially satisfies the IPC. Moving forward, we will write ${\cal S}{\cal P}$
when referring to the set of feasible solutions to (SP). We denote the sets of
interior feasible, and optimal solutions by $\operatorname{int}({\cal S}{\cal
P})$ and ${\cal S}{\cal P}_{*}$, respectively.

The next result from [RTV05] asserts that we obtain optimal solutions to (P)
and (D) by solving (SP).

######  Theorem 1 (Theorem I.6 in [RTV05]).

The system (6) has a solution with $\vartheta=0$ and $\beta>0$ if and only if
the problem (SP) has an optimal solution with $\beta=z_{\mkern
1.5mu\overline{\mkern-1.5mun}-1}>0$.

####  2.1.3 The central path and its neighborhood

For any vector $z\in\mathbb{R}^{\mkern 1.5mu\overline{\mkern-1.5mun}}$, we
define its slack vector444As we have defined it here, $s$ really is a surplus
variable.

$s(z)\coloneqq Mz+q.$

It follows that

$z\leavevmode\nobreak\ \textit{is a feasible solution to \eqref{e:self-
dual}}\leavevmode\nobreak\ \iff\leavevmode\nobreak\ z\geq
0\leavevmode\nobreak\ \textit{and}\leavevmode\nobreak\ s(z)\geq 0.$

Recalling that $e$ denotes the all-one vector of length $\mkern
1.5mu\overline{\mkern-1.5mun}$, note that $s(e)=e$ and $e\odot s(e)=e$. Hence,
$\mu=1$ for the point $\left(z,s(z)\right)=(e,e)\in\operatorname{int}({\cal
S}{\cal P})$.

For every positive $\mu$ there exists a unique non-negative vector $z$ such
that

$\displaystyle z\odot s(z)=\mu e,\quad z\geq 0,\leavevmode\nobreak\ s(z)\geq
0,$ (7)

see, [RTV05, Lemma I.13]. We denote the unique non-negative solution to the
quadratic system (7) as $z(\mu)$, which we refer to as the $\mu$-center. Using
this notation, we may write $z(1)=e$. The set of $\mu$-centers ${\cal
C}=\left\\{z(\mu):\mu>0\right\\}$ is the central path of the problem (SP), and
an optimal solution of (SP) is the limiting point of the central path as
$\mu\to 0^{+}$. The central path constitutes an analytic curve
$\xi(\mu):(0,\infty)\mapsto\mathbb{R}^{\mkern 1.5mu\overline{\mkern-1.5mun}}$,
and its graph $(\mu,\xi(\mu))$ satisfies

$\displaystyle z\odot s(z)=\mu e,\quad z>0,\leavevmode\nobreak\ s(z)>0.$ (8)

Since IPC is satisfied, the central path is guaranteed to exist and is
uniquely determined from the starting point $\left(z(1),s(z(1))\right)=(e,e)$.
Put another way, $(e,e)$ is the point on the central path corresponding to
$\mu=1.$

The standard IPM can be interpreted as an algorithm that approximately follows
the central path by locally linearizing the quadratic system (7). The solution
of each locally linearized system is the so-called Newton step. To see this,
suppose that $z$ is a positive solution to (SP) such that its slack vector
$s(z)$ is also positive, i.e., $z$ and $s(z)$ satisfy (8). To find the
displacement $\Delta z$ such that $z^{+}\coloneqq z+\Delta z$ is the
$\mu$-center, we want to solve the equation:

$\displaystyle(z+\Delta z)\odot s(z+\Delta z)=\mu e.$

Defining $s=s(z)$ and $\Delta s=M\Delta z$, we obtain a nonlinear equation
with quadratic term $\Delta z\odot\Delta s$, namely,

$\displaystyle Zs+Z\Delta s+S\Delta z+\operatorname{diag}\left(\Delta
z\right)\Delta s=\mu e,$ (9)

where $Z\coloneqq\operatorname{diag}(z)$ and
$S\coloneqq\operatorname{diag}(s)$. If we further assume $z^{+}$ is in a small
neighborhood of $z$, i.e., $\Delta z$ is small (and hence, so is $\Delta
s=M\Delta z$), the quadratic term in (9) can be omitted, and the nonlinear
equation reduces to a linear system in the unknowns $\Delta z$ and $\Delta s$.

The foundation for IPM theory is that performing the local linearization we
just described at each iterate, and solving the resulting equation system for
$\Delta z$ and $\Delta s$ allows us to make sufficient progress towards the
optimal solution (and remain in $\operatorname{int}({\cal S}{\cal P})$),
provided that the current iterate $(z,s(z))$ is in some sense close to the
central path. To quantify a notion of closeness, one defines a proximity
measure that gives rise to a neighborhood of the central path. In this paper,
we consider the distance metric

$d_{2}(x,s;\mu)\coloneqq\left\|z\odot s(z)-\mu e\right\|_{2},$

which for fixed $\gamma\in(0,1)$ gives rise to a narrow neighborhood of the
central path:

${\cal N}_{2}(\gamma)\coloneqq\left\\{(z,s(z))\in\operatorname{int}({\cal
S}{\cal P}):d_{2}(x,s;\mu)\leq\gamma\mu\right\\}.$

One can observe that for every $\gamma\in(0,1)$, the following set of
inclusions hold:

${\cal C}\subset{\cal N}_{2}(\gamma)\subset\operatorname{int}({\cal S}{\cal
P}).$

Before proceeding further, we provide two results that will be useful in our
analysis later in the paper. For ease of notation, we define $F(z)\coloneqq
z\odot s(z)\in\mathbb{R}^{\mkern 1.5mu\overline{\mkern-1.5mun}}$. We begin by
establishing that the Jacobian of $F(z)$ is non-singular over the positive
orthant.

######  Lemma 1.

Let $Z\coloneqq\operatorname{diag}(z)$, $S\coloneqq\operatorname{diag}(s(z))$
and ${\cal J}(z)\coloneqq\frac{\partial F(z)}{\partial z}$. Then,

$\displaystyle{\cal J}(z)=ZM+S$ (10)

is non-singular whenever $z,s(z)>0$.

###### Proof.

We use the following identities:

$F(z)=z\odot s(z)=z\odot(Mz+q)=\operatorname{diag}(z)Mz+q\odot z.$

To show that ${\cal J}(z)=ZM+S$ is non-singular whenever $z,s>0$, note that in
this case $D=Z^{-1}S$ is a positive diagonal matrix, and hence

$u^{\top}(M+D)u=u^{\top}Mu+u^{\top}Du=0\iff u^{\top}Du=0\iff u=0,$

as $u^{\top}Mu=0$ for all $u\in\mathbb{R}^{\mkern
1.5mu\overline{\mkern-1.5mun}}$ since $M$ is skew-symmetric. ∎

The result in Lemma 1 allows us to analyze the singular values of the Jacobian
${\cal J}(z)$ of $F(z)$, and thus the spectrum of the matrix ${\cal
J}(z)^{\top}{\cal J}(z)$, which will play a crucial role in the quantum
Hamiltonian we define later in Section 3.

######  Lemma 2.

Let $\mu>0$ and $z(\mu)$ denote the $\mu$-center. Then, the smallest singular
value of

${\cal J}\left(z(\mu)\right)\coloneqq\left.\frac{\partial F(z)}{\partial
z}\right|_{z=z(\mu)}$

satisfies

$\sigma_{\min}\left({\cal J}\left(z(\mu)\right)\right)\geq\eta\mu>0,$

where $\eta>0$.

###### Proof.

Let $s(\mu)\coloneqq s(z(\mu))$. Note that if $z(\mu)$ is the $\mu$-center,
then we must have

$\operatorname{diag}\left(z(\mu)\right)\operatorname{diag}\left(s\left(z(\mu)\right)\right)=Z(\mu)S(\mu)=\mu
I,$

with $z(\mu),s\left(\mu\right)>0$. Since $Z(\mu)$ and $S(\mu)$ are positive
diagonal matrices, for any positive integer $k$ we may write

$\prod_{j=1}^{k}Z(\mu)=Z(\mu)^{k}\coloneqq\operatorname{diag}\left(\left(z(\mu)\right)_{1}^{k},\dots,\left(z(\mu)\right)_{n}^{k}\right),\leavevmode\nobreak\
\prod_{j=1}^{k}S(\mu)=S(\mu)^{k}\coloneqq\operatorname{diag}\left(\left(s(\mu)\right)_{1}^{k},\dots,\left(s(\mu)\right)_{n}^{k}\right).$

Since $z(\mu)$ and $s\left(z(\mu)\right)$ are strictly positive, Lemma 1
ensures that ${\cal J}\left(z(\mu)\right)$ is non-singular, and as a
consequence, ${\cal J}\left(z(\mu)\right)^{\top}{\cal J}\left(z(\mu)\right)$
is positive definite. Combining these facts, and using that $M$ is skew-
symmetric, we obtain

$\displaystyle{\cal J}\left(z(\mu)\right)^{\top}{\cal J}\left(z(\mu)\right)$
$\displaystyle=\left(Z(\mu)M+S(\mu)\right)^{\top}\left(Z(\mu)M+S(\mu))\right)$
$\displaystyle=M^{\top}Z(\mu)^{2}M+M^{\top}Z(\mu)S(\mu)+S(\mu)Z(\mu)M+S(\mu)^{2}$
$\displaystyle=M^{\top}Z(\mu)^{2}M-MZ(\mu)S(\mu)+Z(\mu)S(\mu)M+S(\mu)^{2}$
$\displaystyle=M^{\top}Z(\mu)^{2}M-\mu MI+\mu IM+S(\mu)^{2}$
$\displaystyle=M^{\top}Z(\mu)^{2}M+S(\mu)^{2}\succ 0.$

Since $S^{2}(\mu)$ is a positive diagonal matrix, we may write

${\cal J}\left(z(\mu)\right)^{\top}{\cal
J}\left(z(\mu)\right)=M^{\top}Z(\mu)^{2}M+S(\mu)^{2}\succeq S(\mu)^{2}\succ
0.$

To see this, note that for any $u\in\mathbb{R}^{\mkern
1.5mu\overline{\mkern-1.5mun}}$ we have

$u^{\top}\left(M^{\top}Z(\mu)^{2}M\right)u=(uM)^{\top}Z(\mu)^{2}(Mu)\geq 0.$

Hence, an application of Weyl’s inequality yields

$\lambda_{\min}\left(M^{\top}Z(\mu)^{2}M+S(\mu)^{2}\right)\geq\lambda_{\min}\left(M^{\top}Z(\mu)^{2}M\right)+\lambda_{\min}\left(S(\mu)^{2}\right)\geq\lambda_{\min}\left(S(\mu)^{2}\right)=\min_{i\in[\mkern
1.5mu\overline{\mkern-1.5mun}]}\left\\{\left(S(\mu)^{2}\right)_{ii}\right\\}.$

The diagonal elements of $S(\mu)^{2}$ are bounded away from zero by a positive
constant times $\mu^{2}$ (see, e.g., [RTV05, Section 3.3.3]), and therefore we
clearly have

$\displaystyle\sigma_{\min}\left({\cal
J}\left(z(\mu)\right)\right)=\sqrt{\lambda_{\min}\left(M^{\top}Z(\mu)^{2}M+S(\mu)^{2}\right)}$
$\displaystyle\geq\sqrt{\lambda_{\min}\left(S(\mu)^{2}\right)}$
$\displaystyle=\min_{i\in[\mkern
1.5mu\overline{\mkern-1.5mun}]}\left\\{\left(S(\mu)\right)_{ii}\right\\}=\eta\mu>0,$
(11)

for some constant $\eta>0$. The proof is complete. ∎

###  2.2 Quantum adiabatic theorem for unbounded Hamiltonians

Given a quantum Hamiltonian ${\cal H}(t)$, $t\in[0,1]$, we consider the
dynamics described by the Schrödinger equation:

$\displaystyle\mathbf{i}\eta\frac{\text{d}}{\text{d}t}\left|{\psi(t)}\right\rangle={\cal
H}(t)\left|{\psi(t)}\right\rangle,$ (12)

where $\mathbf{i}$ denotes the imaginary unit and $\eta$ is a positive real
number. We suppose that $U(t)$ is the propagator of the dynamics, such that
the solution to (12) at time $t$ is given by

$\left|{\psi(t)}\right\rangle=U(t)\left|{\psi(0)}\right\rangle.$

If $\left|{\psi(0)}\right\rangle$ is a nondegenerate ground state of ${\cal
H}(0)$, then the quantum adiabatic theorem [Mes58] asserts that in the limit
$T\to\infty$, the state $\left|{\psi(T)}\right\rangle$ obtained from (12),
will be close to the ground state of ${\cal H}(T)$.

We assume ${\cal H}(t)$ has a non-degenerate ground state for all $t\in[0,1]$,
and we define $P(t)$ as a rank-$1$ projector onto the ground-energy subspace
of ${\cal H}(t)$.

######  Theorem 2 (Exponential estimate).

Define ${\cal H}^{(k)}(t)\coloneqq\frac{\text{d}^{k}}{\text{d}t^{k}}{\cal
H}(t)$. Let ${\cal H}(t)$ be a quantum Hamiltonian for $t\in[0,1]$ such that
the following conditions hold:

  1. (a)

${\cal H}(t)$ admits an analytic continuation to some strip on $\mathbb{C}$
containing $[0,1]$.

  2. (b)

For $t\in[0,1]$, the spectral gap of ${\cal H}(t)$ is greater than a constant
$\Delta_{0}>0$.

  3. (c)

For any $k=1,2,\dots$, we have that ${\cal H}^{(k)}(0)={\cal H}^{(k)}(1)=0$.

Then, there exists a constant $C$ such that

$\|\psi(1)-\phi_{*}\|\leq Ce^{-1/\eta},$

where $\psi(1)$ is the solution to (12) at $t=1$ and $\phi_{*}$ is a ground
state of ${\cal H}(1)$, i.e., $P(1)\phi_{*}=\phi_{*}$.

###### Proof.

Combining the main theorem in [HJ02] and the vanishing derivative condition at
$t=0,1$. ∎

###  2.3 Quantum algorithms for Schrödinger equations

The primary subroutine of our algorithm is simulating the Schrödinger
equation, otherwise known as quantum simulation. Our work here will consider
the Schrödinger equation over the time interval $[t_{0},t_{1}]$ for a given
time-dependent potential $V(t,x)$,

$\displaystyle\mathbf{i}\frac{\partial}{\partial
t}\left|{\Psi(x,t)}\right\rangle=\left[-\frac{1}{2}\nabla^{2}+V(x,t)\right]\left|{\Psi(x,t)}\right\rangle,$
(13)

where we specify $\Omega=[-R,R]^{d}$ for a sufficiently large $R$ and
$V(x,t)\colon\Omega\times[t_{0},t_{1}]\to\mathbb{R}$ is a time-dependent
potential function. Accordingly,
$\Psi(t,x)\colon\Omega\times[t_{0},t_{1}]\to\mathbb{C}$ is the wave function
subject to certain initial data $\Psi(x,t_{0})=\Psi_{0}(x)$ and the periodic
boundary condition.

We utilize the quantum simulation algorithm of Childs, Leng, Li, Liu, and
Zhang [CLL+22a], which exhibits near-optimal dependence in the dimension $d$,
and precision $\epsilon$ to which the simulation is carried out. However, the
gate complexity of that quantum simulation algorithm also involves a parameter
$g^{\prime}$ that depends on the higher-order derivatives of the wave
function. An accurate upper bound of $g^{\prime}$ relies on a refined _a
priori_ estimate of the wave function. By leveraging the fact that the initial
condition is always analytic in our setting, we eschew the regularity
parameter $g^{\prime}$ in [CLL+22a, Theorem 8], resulting in an enhanced
complexity result.

######  Theorem 3 (Improved version of Theorem 8 in [CLL+22a]).

Suppose the potential field $V(x,t)$ is bounded, smooth in $x$ and $t$, and
periodic in $x$. Define the $\|\cdot\|_{\infty,1}$-norm of $V(x,t)$ as

$\displaystyle\|V\|_{\infty,1}\coloneqq\int^{t_{1}}_{t_{0}}\|V(\cdot,t)\|_{\infty}\leavevmode\nobreak\
\textup{d}t.$

We assume that we have access to the zeroth-order oracle of $V$, which is a
unitary map $O_{V}$ on $\Omega\otimes[t_{0},t_{1}]\to\mathbb{R}$ such that for
any $\left|{x}\right\rangle\in\Omega$ and
$\left|{s}\right\rangle\in[t_{0},t_{1}]$,

$O_{V}\left(\left|{x}\right\rangle\otimes\left|{s}\right\rangle\otimes\left|{0}\right\rangle\right)=\left|{x}\right\rangle\otimes\left|{s}\right\rangle\otimes\left|{V(x,s)}\right\rangle.$

Moreover, we assume the initial data $\Psi_{0}(x)$ is analytic on $\Omega$ and
$V$ is $L$-Lipschitz in $t$. Then, the Schrödinger equation (13) can be
simulated for time $t\in[t_{0},t_{1}]$ up to accuracy $\epsilon$ with the
following cost:

  1. 1.

Queries to $O_{V}$: ${\cal
O}\left(\|V\|_{\infty,1}\frac{\log(\|V\|_{\infty,1}/\epsilon)}{\log\log(\|V\|_{\infty,1}/\epsilon)}\right)$,

  2. 2.

1- and 2-qubit gates:

${\cal
O}\left(\|V\|_{\infty,1}\left(\operatorname{poly}(z)+\log^{2.5}\left(L\|V\|_{\infty,1}/\epsilon\right)+d\log\log(1/\epsilon)\right)\frac{\log(\|V\|_{\infty,1}/\epsilon)}{\log\log(\|V\|_{\infty,1}/\epsilon)}\right).$

We point out that the simulation accuracy referred to in Theorem 3 corresponds
to the $\ell_{2}$-distance between the actual final state and the state
returned by the quantum simulation algorithm. This theorem is an improved
version of the original Theorem 8 in [CLL+22a]. We provide the proof of this
theorem in Appendix B. We also refer the readers to [LZW23, Section 2.4] for a
detailed discussion.

##  3 A quantum algorithm that traces the central path

In this section, we propose a Hamiltonian formalism for the central path of
linear optimization problems. We show how one can solve LO problems by
simulating the Schrödinger equation associated with our Hamiltonian, and
provide a rigorous complexity analysis of the resulting scheme.

###  3.1 Quantum representation of the central path

In this subsection, we construct a family of quantum Hamiltonian operators
${\cal H}(\mu)$, where $\mu\in(0,1]$ is a positive parameter. We will show
that, for each $0<\mu\leq 1$, the ground state of the operator ${\cal H}(\mu)$
gives rise to a Gaussian distribution centered at the $\mu$-center $z(\mu)$.
Therefore, by measuring the ground state of ${\cal H}(\mu)$, we can
approximate the $\mu$-center.

Recall that we define $F(z)\coloneqq z\odot s(z)\in\mathbb{R}^{\mkern
1.5mu\overline{\mkern-1.5mun}}$, where $s(z)=Mz+q$. For any $0<\mu\leq 1$, we
define the function

$\displaystyle f_{\mu}(z)\coloneqq\frac{1}{2}\left(F(z)-\mu
e\right)^{\top}\left(F(z)-\mu e\right)=\frac{1}{2}\left\|F(z)-\mu
e\right\|^{2}.$ (14)

######  Definition 2 (Central-path Hamiltonian).

For $0<\mu\leq 1$ and $h>0$, we define a 1-parameter family of elliptic
operators over the positive orthant $\mathbb{R}^{\mkern
1.5mu\overline{\mkern-1.5mun}}_{++}$,

$\displaystyle{\cal H}(\mu)=-\frac{h^{2}}{2}\nabla^{2}+f_{\mu}(z),$ (15)

where $\nabla^{2}=\sum^{\mkern
1.5mu\overline{\mkern-1.5mun}}_{j=1}\frac{\partial^{2}}{\partial z_{j}^{2}}$
is the Laplacian operator.

As defined in (14), the function $f_{\mu}(z)$ is non-negative. The
$\mu$-center $z(\mu)$ is the unique zero of $f_{\mu}(z)$ in the positive
orthant $z>0$. When the parameter $h$ is sufficiently small, the harmonic
approximation theory of Schrödinger operators [HS12, Section 11] allows us to
approximate the low-energy spectrum of ${\cal H}(\mu)$ by that of the
following operator,

$\displaystyle\widetilde{{\cal
H}}(\mu)=-\frac{h^{2}}{2}\nabla^{2}+\frac{1}{2}\left[z-z(\mu)\right]^{\top}H(\mu)\left[z-z(\mu)\right],$
(16)

where $H(\mu)$ is the Hessian of $f_{\mu}(z)$ at $z=z(\mu)$. A quick
calculation yields that

$H(\mu)={\cal J}(z(\mu))^{\top}{\cal J}(z(\mu))\succ 0.$ (17)

In the following lemma, we summarize some important properties of the Hessian
matrix $H(\mu)$.

######  Lemma 3.

Let $0<\mu\leq 1$ and define ${\cal
M}\coloneqq\begin{bmatrix}M&I\end{bmatrix}^{\top}\in\mathbb{R}^{2\mkern
1.5mu\overline{\mkern-1.5mun}\times\mkern 1.5mu\overline{\mkern-1.5mun}}$.
Then, the condition number of $H(\mu)$ is bounded by as

$\kappa\left(H(\mu)\right)={\cal O}\left(\kappa({\cal
M})^{2}\frac{1}{\mu^{2}}\right),$

where $\kappa({\cal M})$ is an absolute constant that only depends on the LO
problem data $(A,b,c)$. Moreover, there exists constants $\eta_{1},\eta_{2}>0$
such that

$\eta_{1}\mu^{2}\leq\lambda_{0}(H(\mu))\leq\eta_{2}\mu^{2}.$

###### Proof.

In order to bound the condition number of $H(\mu)$, recall that

$H(\mu)={\cal J}(z(\mu))^{\top}{\cal
J}(z(\mu))=M^{\top}Z(\mu)^{2}M+S(\mu)^{2}.$

We can therefore express $H(\mu)$ as

$H(\mu)=\begin{bmatrix}M\\\
I\end{bmatrix}^{\top}\begin{bmatrix}Z(\mu)^{2}&0\\\
0&S(\mu)^{2}\end{bmatrix}\begin{bmatrix}M\\\ I\end{bmatrix}.$

Taking ${\cal M}\coloneqq\begin{bmatrix}M\\\ I\end{bmatrix}$ and
$D\coloneqq\begin{bmatrix}Z(\mu)^{2}&0\\\ 0&S(\mu)^{2}\end{bmatrix}$, it can
be readily seen that $H(\mu)={\cal M}^{\top}D{\cal M}$, from which it follows

$\kappa\left(H(\mu)\right)={\cal O}\left(\kappa\left({\cal M}^{\top}{\cal
M}\right)\kappa(D)\right)={\cal O}\left(\kappa({\cal
M})^{2}\kappa(D)\right)={\cal O}\left(\kappa({\cal
M})^{2}\frac{1}{\mu^{2}}\right).$

To understand why ${\cal O}\left(\kappa(D)\right)={\cal
O}\left(\frac{1}{\mu^{2}}\right)$, recall that (see, the proof of Lemma 2) the
smallest elements of $Z(\mu)^{2}$ and $S(\mu)^{2}$ are bounded away from zero
by a positive constant times $\mu^{2}$ [RTV05, Chapter 3.3.3], and the
smallest eigenvalue of $D$ is given by

$\min\left\\{\min_{i\in[\mkern
1.5mu\overline{\mkern-1.5mun}]}\left\\{\left(Z(\mu)^{2}\right)_{ii}\right\\},\min_{i\in[\mkern
1.5mu\overline{\mkern-1.5mun}]}\left\\{\left(S(\mu)^{2}\right)_{ii}\right\\}\right\\}.$

In order to prove the second part of the statement, recall from the proof of
Lemma 2 that the smallest eigenvalue of is of the form

$\lambda_{\min}\left(H(\mu)\right)=\lambda_{\min}\left({\cal
J}\left(z(\mu)\right)^{\top}{\cal
J}\left(z(\mu)\right)\right)\geq\eta^{2}\mu^{2},$

where $\eta>0$ is defined according to (11). Therefore, there must exist
$\eta_{1},\eta_{2}>0$ such that

$\eta_{1}\mu^{2}\leq\lambda_{\min}\left(H(\mu)\right)\leq\eta_{2}\mu^{2},$

for example, we can choose

$\eta_{1}=\eta^{2},\quad\eta_{2}=2\eta^{2}.$

More succinctly,

$\lambda_{\min}\left(H(\mu)\right)=\Theta\left(\mu^{2}\right),$

as desired. ∎

The Hamiltonian operator $\widetilde{{\cal H}}(\mu)$ describes a quantum
harmonic oscillator. The eigenvalues and eigenstates of a quantum harmonic
oscillator are well understood [GS18]. Let
$0<\lambda_{0}(\mu)\leq\lambda_{1}(\mu)\leq\cdots\leq\lambda_{\mkern
1.5mu\overline{\mkern-1.5mun}-1}(\mu)$ be the eigenvalues of the Hessian
matrix $H(\mu)$. The spectral gap (i.e., the difference between the first two
eigenvalues) of the operator $\widetilde{{\cal H}}(\mu)$ is
$\Delta=h\lambda^{1/2}_{0}(\mu)$. Moreover, let $\Phi_{0}(\mu)$ be the ground
state of the operator $\widetilde{{\cal H}}(\mu)$. It turns out that
$|\Phi_{0}(\mu)|^{2}$ is the probability density function of the multivariate
normal distribution

$\displaystyle|\Phi_{0}(\mu)|^{2}\sim{\cal
N}\left(z(\mu),\frac{h}{2}(H(\mu))^{-1/2}\right).$ (18)

######  Proposition 1.

Fix $\delta>0$. If we choose

$\displaystyle h\coloneqq\frac{\gamma^{2}\mu}{2\sqrt{\eta_{2}\kappa}(\mkern
1.5mu\overline{\mkern-1.5mun}/2+3\log(1/\delta)/4)},$ (19)

then, for any $\mu\in(0,1]$, we have that

$\displaystyle\Pr_{x\sim|\Phi_{0}(\mu)|^{2}}\left[x\not\in{\cal
N}_{2}(\gamma)\right]\leq\delta.$

The proof of Proposition 1 is available in Appendix A.

###  3.2 Quantum simulation of the central path

The key idea of our quantum algorithm is to simulate a quantum evolution for
some $t\in[0,1]$ in which the quantum state $\left|{\Psi(t)}\right\rangle$ is
(approximately) the ground state of ${\cal H}(\mu(t))$. Here, the function
$\mu(t)$ is a monotonically decreasing function such that

$\mu(0)=1,\quad\mu(1)=\mu_{f}\ll 1.$

In this way, the quantum state $\left|{\Psi(t)}\right\rangle$ follows the
central path as $\mu(t)$ decreases. At $t=1$, if we measure the final state
$\left|{\Psi(1)}\right\rangle$, Proposition 1 guarantees that we will obtain
an approximate solution that is in the neighborhood of $z(\mu_{f})$ with high
probability.

First, we introduce a function $g\colon[0,1]\mapsto[0,1]$ such that (i) $g$ is
analytic in $(0,1)$, and (ii) for any $k=1,2,\dots$, we have
$g^{(k)}(0)=g^{(k)}(1)=0$, where
$g^{(k)}(t)\coloneqq\frac{\text{d}^{k}}{\text{d}t^{k}}g(t)$. The construction
of this function $g(t)$ follows [AL22, Equation 8]. Namely, we choose

$\displaystyle
g(t)=c^{-1}_{e}\int^{t}_{0}\exp\left(-\frac{1}{\tau(1-\tau)}\right)\leavevmode\nobreak\
\text{d}\tau,$ (20)

where
$c_{e}=\int^{1}_{0}\exp\left(-\frac{1}{\tau(1-\tau)}\right)\leavevmode\nobreak\
\text{d}\tau$ is a normalization constant such that $g(1)=1$.

Then, given a fixed $\mu_{f}\in(0,1)$, we define the function $\mu(t)$,

$\displaystyle\mu(t)=1-(1-\mu_{f})g(t).$ (21)

Clearly, $\mu(t)$ is a monotonically decreasing function such that $\mu(0)=1$,
$\mu(1)=\mu_{f}$, and we have $\mu^{(k)}(0)=\mu^{(k)}(1)=0$ for any positive
integer $k$.

######  Lemma 4.

Let $\mu(t)$ be the same as in (21). Then, the Hamiltonian ${\cal
H}(t)\coloneqq{\cal H}(\mu(t))$ is analytic in $t\in(0,1)$. Moreover, for any
$k=1,2,\dots$, we have that ${\cal H}^{(k)}(0)={\cal H}^{(k)}(1)=0$, where
${\cal H}^{(k)}(t)\coloneqq\frac{\text{d}^{k}}{\text{d}t^{k}}{\cal H}(t)$.

###### Proof.

The analyticity immediately follows from our definition of ${\cal H}$ (see
(15)) because $f_{\mu}$ is a quadratic function in $\mu$. By the chain rule,
we have

$\frac{\text{d}{\cal H}(t)}{\text{d}t}=-2\dot{\mu}\left(F(z)-\mu
e\right)^{\top}e.$

Since $\dot{\mu}(0)=\dot{\mu}(1)=0$, we have ${\cal H}^{\prime}(0)={\cal
H}^{\prime}(1)=0$. Similarly, using an induction argument, we can prove ${\cal
H}^{(k)}(0)={\cal H}^{(k)}(1)=0$ for any $k\geq 1$. ∎

######  Proposition 2.

Let ${\cal H}(t)\coloneqq{\cal H}(\mu(t))$, where $\mu(t)$ is defined in (21).
For sufficiently small $h$, we consider the following Schrödinger equation,

$\displaystyle\mathbf{i}\eta\frac{\partial}{\partial
t}\Psi^{\eta}(t)=\frac{1}{h\mu(t)}{\cal H}(t)\Psi^{\eta}(t),$ (22)

where the initial state $\Psi(0)$ is the ground state of ${\cal H}(0)$. Then,
there exists a constant $C$ such that

$\displaystyle\|\Psi^{\eta}(1)-\Phi_{0}(\mu_{f})\|\leq Ce^{-1/\eta},$ (23)

where $\Phi_{0}(\mu_{f})$ is a ground state of the final Hamiltonian ${\cal
H}(1)$.

###### Proof.

This is a direct consequence of Theorem 2. In Lemma 4, we proved that ${\cal
H}(t)$ is analytic in $t\in(0,1)$ and and ${\cal H}^{(k)}(0)={\cal
H}^{(k)}(1)=0$ for any positive integer $k$. It remains to show that the
spectral gap of ${\cal H}(t)$ has a lower bound $\Delta_{0}$. For sufficiently
small $h$, the spectral gap of ${\cal H}(t)$ is approximately the same as its
harmonic approximation $\widetilde{{\cal H}}(\mu(t))$ (see (16)). The operator
$\widetilde{{\cal H}}(\mu(t))$ describes a quantum harmonic oscillator and its
spectral gap is precisely $\Delta(t)=h\lambda^{1/2}_{0}(\mu(t))$. It follows
from Lemma 3 that

$\Delta(t)\geq h\eta^{1/2}_{1}\mu(t).$

Therefore, the spectral gap of the Hamiltonian operator in (22) has a lower
bound $\Delta_{0}=\eta^{1/2}_{1}$. ∎

###  3.3 Quantum central path algorithm

We present the quantum central path method for linear optimization in full
detail in Algorithm 1. The algorithm takes as input: (i) the LO problem data
$(A,b,c)$; (ii) the optimality tolerance $\varepsilon\in(0,1)$ to which we
seek to solve (P)-(D); (iii) the neighborhood opening parameter
$\gamma\in(0,1)$ that specifies ${\cal N}_{2}(\gamma)$, the
$d_{2}$-neighborhood of the central path; and (iv) the failure rate
$\delta\in(0,1)$ to which we allow the quantum algorithm returns a point that
is not in ${\cal N}_{2}(\gamma)$. In our initialization steps, we form the
self-dual embedding formulation associated to the problem data $(A,b,c)$, see,
e.g., Section 2.1.2. Proceeding in this way, we may provide our quantum
central path algorithm with the initial quantum state centered at the strictly
feasible point starting point $e\in\mathbb{R}^{\mkern
1.5mu\overline{\mkern-1.5mun}}$ with associated duality gap $\mu=1$.

In Algorithm 1, we simulate the Schrödinger equation (22) to error
$\frac{\delta}{8}$, which can be accomplished through the use of

$\widetilde{{\cal O}}_{m,n,\frac{1}{\delta}}\left((m+n)\kappa({\cal
M})\frac{1}{\varepsilon}\right)$

queries to an evaluation oracle for $f_{\mu}(z)$ and

$\widetilde{{\cal O}}_{m,n,\kappa({\cal
M}),\frac{1}{\delta}}\left((m+n)^{2}\kappa({\cal
M})\frac{1}{\varepsilon}\right)$

elementary gates. We address the construction an evaluation oracle for
$f_{\mu}(z)$ and its associated cost in Section 3.4. Letting
$\Phi_{0}(\mu_{f})$ denote the ground state associated to the Hamiltonian
${\cal H}(\mu_{f})$ and choosing $h$ according to Equation (19) simulating the
Shrödinger equation, suffices to ensure that $|\Phi_{0}(\mu_{f})|^{2}$ is the
probability density function of the multivariate normal distribution ${\cal
N}\left(z\left(\mu_{f}\right),\frac{h}{2}(H(\mu_{f}))^{-1/2}\right)$, with
$\mu_{f}\leq\varepsilon$. This process can be visualized in Figure 1.

Input: Matrix $A\in\mathbb{R}^{m\times n}$, vectors $b\in\mathbb{R}^{m}$ and
$c\in\mathbb{R}^{n}$, optimality tolerance $\varepsilon\in(0,1)$, neighborhood
opening $\gamma\in(0,1)$, failure rate $\delta\in(0,1)$

Output: $\varepsilon$-optimal solution $z\in\mathbb{R}^{\mkern
1.5mu\overline{\mkern-1.5mun}}$ to the LO problem (SP)

Initialize: Form self-dual embedding $(M,q)$ for $(A,b,c)$,
$\mu_{f}\leftarrow\varepsilon$

  1. 1.

Initialize the quantum register to $\left|{\Phi(0)}\right\rangle$, the ground
state of ${\cal H}(0)$.

  2. 2.

$\left|{\Psi^{\eta}_{\mathrm{sim}}(1)}\right\rangle\leftarrow$ simulate the
Schrödinger equation (22) to error $\delta/8$ with

$h(t)=\frac{\gamma^{2}\mu(t)}{2\sqrt{\eta_{2}\kappa}(\mkern
1.5mu\overline{\mkern-1.5mun}/2+3\log(2/\delta)/4)},\quad\eta=\frac{1}{\log(8C/\delta)},$

where $\mu(t)$ is defined in (21) and the constant $C$ is independent of $m,n$
and $\delta$ (see Proposition 2).

  3. 3.

$z\leftarrow$ sample from the quantum state
$\left|{\Psi^{\eta}(1)}\right\rangle$.

Algorithm 1 Quantum central path algorithm for linear optimization

We may now establish the correctness and complexity of Algorithm 1.

######  Theorem 4.

Let $\varepsilon\in(0,1)$. Choose $\gamma\in(0,1)$ and $\delta\in(0,1)$.
Define $\kappa({\cal M})$ to be the condition number of the matrix ${\cal
M}=\begin{bmatrix}M&I\end{bmatrix}^{\top}$, which is an absolute constant that
depends only on the LO problem data $(A,b,c)$. Assume access to a quantum
oracle $O_{f_{\mu}}$ evaluates the potential function $f_{\mu}(z,t)$ in a
superposition:

$O_{f_{\mu}}\left(\left|{z}\right\rangle\otimes\left|{t}\right\rangle\otimes\left|{0}\right\rangle\right)=\left|{z}\right\rangle\otimes\left|{t}\right\rangle\otimes\left|{f_{\mu}(z,t)}\right\rangle.$

Then, with probability at least $1-\delta$, Algorithm 1 returns a classical
vector $z\in\mathbb{R}^{m+n}_{++}$ such that $(z,s(z))\leavevmode\nobreak\
\in\leavevmode\nobreak\ {\cal N}_{2}(\gamma)$ with
$d_{2}(z,s(z);\varepsilon)\leq\gamma\varepsilon$.

Moreover, Algorithm 1 can be implemented with

${\cal O}\left((m+n)\kappa({\cal
M})\frac{1}{\varepsilon}\cdot\textup{polylog}\left(m,n,\kappa({\cal
M}),\frac{1}{\delta}\right)\right)$

queries to $O_{f_{\mu}}$ and

${\cal O}\left((m+n)^{2}\kappa({\cal
M})\frac{1}{\varepsilon}\cdot\textup{polylog}\left(m,n,\kappa({\cal
M}),\frac{1}{\delta}\right)\right)$

elementary gates.

###### Proof.

We obtain a quantum state $\left|{\Psi^{\eta}_{\mathrm{sim}}(1)}\right\rangle$
by simulating the Schrödinger equation (12) to error $\delta/8$. By
Proposition 2, if we choose $\eta=1/\log(8C/\delta)$ (where $C$ is a constant
independent of $d$ and $\delta$), we have
$\left\|\Psi^{\eta}(1)-\Phi_{0}(\mu_{f})\right\|\leq\delta/8$. By the triangle
inequality, it turns out that

$\displaystyle\left\|\Psi^{\eta}_{\mathrm{sim}}(1)-\Phi_{0}(\mu_{f})\right\|\leq\frac{\delta}{4}.$
(24)

We define $B_{\gamma,\varepsilon}\coloneqq\\{z\in\mathbb{R}^{\mkern
1.5mu\overline{\mkern-1.5mun}}:f(z)\leq\frac{1}{2}\gamma^{2}\varepsilon^{2}\\}$,
which is a compact set containing the point $z(\varepsilon)$. Let
$\mathbbm{1}_{B_{\gamma,\varepsilon}}(z)$ be the indicator function of the set
$B_{{\gamma,\varepsilon}}$, then we have

$\displaystyle\Pr_{z\sim|\Psi^{\eta}_{\mathrm{sim}}(1)|^{2}}\left[f_{\mu}(z)>\frac{1}{2}\gamma^{2}\varepsilon^{2}\right]$
$\displaystyle=1-\left\langle{\Psi^{\eta}_{\mathrm{sim}}(1)}\right|\mathbbm{1}_{B_{\gamma,\varepsilon}}\left|{\Psi^{\eta}_{\mathrm{sim}}(1)}\right\rangle,$
$\displaystyle\Pr_{z\sim|\Phi_{0}(\mu_{f})|^{2}}\left[f_{\mu}(z)>\frac{1}{2}\gamma^{2}\varepsilon^{2}\right]$
$\displaystyle=1-\left\langle{\Phi_{0}(\mu_{f})}\right|\mathbbm{1}_{B_{\gamma,\varepsilon}}\left|{\Phi_{0}(\mu_{f})}\right\rangle.$

Then, it follows from Lemma 7 that

$\displaystyle\Pr_{z\sim|\Psi^{\eta}_{\mathrm{sim}}(1)|^{2}}\left[f_{\mu}(z)>\frac{1}{2}\gamma^{2}\varepsilon^{2}\right]\leq\Pr_{z\sim|\Phi_{0}(\mu_{f})|^{2}}\left[f_{\mu}(z)>\frac{1}{2}\gamma^{2}\varepsilon^{2}\right]+\frac{\delta}{2}\leq\delta,$
(25)

where the last step follows from Proposition 1. This proves the first part of
the theorem.

Next, we discuss the complexity of Algorithm 1. The Schrödinger equation (22)
can be written as follows,

$\mathbf{i}\eta\frac{\partial}{\partial
t}\Psi^{\eta}(t)=\left[-\frac{\theta}{2}\nabla^{2}+\frac{1}{\mu(t)h(t)}f_{\mu}(z)\right]\Psi^{\eta}(t),$

where the coefficient

$\theta=\frac{h(t)}{\mu(t)}=\frac{\gamma^{2}}{\sqrt{\eta_{2}\kappa}\left(\frac{\mkern
1.5mu\overline{\mkern-1.5mun}}{2}+3\log\left(\frac{2/\delta}{4}\right)\right)}.$
(26)

We introduce the change of variable $t\mapsto\frac{\eta}{\theta}\tau$ for
$\tau\in[0,\theta/\eta]$, and we define a new wave function

$\widetilde{\Psi}^{\eta}(\tau)\coloneqq\Psi^{\eta}\left(\frac{\eta\tau}{\theta}\right).$

It turns out that simulating the Schrödinger equation (22) is equivalent to
simulate the following time-dilated Schrödinger equation

$\displaystyle\mathbf{i}\frac{\partial}{\partial\tau}\widetilde{\Psi}^{\eta}(\tau)=\left[-\frac{1}{2}\nabla^{2}+\frac{1}{h^{2}(\tau)}f_{\mu}(z)\right]\widetilde{\Psi}^{\eta}(\tau)$
(27)

for $0\leq\tau\leq\frac{\theta}{\eta}$, where $h(\tau)\coloneqq
h\left(\frac{\theta t}{\eta}\right)$.

The time-dilated Schrödinger equation (27) can be simulated using Theorem 3.
Note that in the quantum simulation, we need to truncate the time-dependent
potential function $f_{\mu}(z)/h^{2}(\tau)$ over a compact domain
$\Omega=[0,R]^{\mkern 1.5mu\overline{\mkern-1.5mun}}$, where $R$ is a large
number such that $\Omega$ contains the neighborhood of the central path (e.g.,
${\cal N}_{2}(\gamma)$). Usually, $R$ can be fixed as an absolute constant. We
denote the truncated time-dependent potential function as

$V(z,\tau)=\frac{1}{h^{2}(\tau)}V_{\mu}(z),$

where $V_{\mu}(z)=f_{\mu}(z)$ for $z\in\Omega$. The
$\|\cdot\|_{\infty,1}$-norm of $V(z,\tau)$ is defined as

$\|V\|_{\infty,1}\coloneqq\int^{\theta/\eta}_{0}\frac{1}{h^{2}(\tau)}\|V_{\mu}(\cdot)\|_{\infty}\leavevmode\nobreak\
\text{d}\tau.$

At each $\tau$, the $\ell_{\infty}$-norm of $V_{\mu}(z)$ could be large
because the function increases when $z$ is far away from the central path.
However, as we are simulating the quantum dynamics in the adiabatic regime and
the quantum state is approximately the ground state of the quantum
Hamiltonian, we can give an improved estimate on $\|V_{\mu}\|_{\infty}$ based
on our knowledge of the ground state. Thanks to our choice of $h(\cdot)$, the
ground state is concentrated in a narrow neighborhood ${\cal N}_{2}(\gamma)$
on which the function $f_{\mu}(z)$ is controlled by $\gamma^{2}\mu^{2}/2$.
Thus, we may assume there is a small constant $K>0$ such that for any
$\mu\in[\varepsilon,1]$,

$\|V_{\mu}\|_{\infty}\leq K\gamma^{2}\mu^{2}.$

It follows that

$\displaystyle\|V\|_{\infty,1}\leq\int^{\theta/\eta}_{0}\frac{K\gamma^{2}\mu^{2}(\tau)}{h^{2}(\tau)}\leavevmode\nobreak\
\text{d}\tau\leq\frac{K\gamma^{2}}{\theta\eta}={\cal
O}\left(\sqrt{\kappa}\log\left(\frac{1}{\delta}\right)\cdot\left(\mkern
1.5mu\overline{\mkern-1.5mun}+\log\left(\frac{1}{\delta}\right)\right)\right),$
(28)

where $\theta$ is defined according to (26). Similarly, we estimate the
Lipschitz constant of $V(x,s)$,

$\displaystyle\left\|\dot{V}\right\|_{\infty}\leq\frac{2\dot{\mu}}{\theta^{2}\mu^{2}}\|V_{\mu}\|_{\infty}\leq\frac{2K\gamma^{2}}{\theta^{2}}={\cal
O}\left(\kappa\left(\mkern
1.5mu\overline{\mkern-1.5mun}+\log\left(\frac{1}{\delta}\right)\right)^{2}\right).$
(29)

Therefore, by Theorem 3, we can simulate the Schrödinger equation (27) using

${\cal O}\left(\sqrt{\kappa}\left(\mkern
1.5mu\overline{\mkern-1.5mun}+\log\left(\frac{1}{\delta}\right)\right)\log\left(\frac{1}{\delta}\right)\log\left(\frac{\mkern
1.5mu\overline{\mkern-1.5mun}}{\delta}\right)\right)$

queries to $O_{V_{\mu}}$, and

$\widetilde{{\cal O}}_{\kappa}\left(\sqrt{\kappa}\left(\mkern
1.5mu\overline{\mkern-1.5mun}+\log^{2.5}\left(\frac{\mkern
1.5mu\overline{\mkern-1.5mun}}{\delta}\right)\right)\left(\mkern
1.5mu\overline{\mkern-1.5mun}+\log\left(\frac{1}{\delta}\right)\right)\log\left(\frac{1}{\delta}\right)\log\left(\frac{\mkern
1.5mu\overline{\mkern-1.5mun}}{\delta}\right)\right)=\widetilde{{\cal
O}}_{m,n,\kappa,\frac{1}{\delta}}\left(\sqrt{\kappa}(m+n)^{2}\right)$

additional gates. Note that we ignore $\log\log(\mkern
1.5mu\overline{\mkern-1.5mun})$ and $\log\log(1/\delta)$ factors in the
big-${\cal O}$ notation.

The stated result then follows from Lemma 3, from which we have

$\kappa={\cal O}\left(\kappa({\cal M})^{2}\frac{1}{\mu^{2}}\right)={\cal
O}\left(\kappa({\cal M})^{2}\left(\frac{1}{\varepsilon}\right)^{2}\right).$

∎

neighborhoodfeasible regioncentral path $\mu\to
0$$z{(1)}$$z(\mu(T))$$z_{*}$$\mu=0$optimal solution Figure 1: Visualization of
the quantum central path algorithm. The dotted lines define the boundary of a
neighborhood of the central path. The dashed circles indicate the progression
of the wave packet from time $t=0$ to $t=T$. The wave packet begins to
concentrate on a small ball centered at $z(\mu(T))$, near the optimal solution
to the linear optimization problem $z_{*}$.

###  3.4 Constructing an evaluation oracle for $f_{\mu}$

In this subsection, we discuss the cost of constructing a quantum oracle
$O_{f_{\mu}}$ that queries the function value of $f_{\mu}(z)$ for
$z\in\Omega$. To construct $O_{f_{\mu}}$, we only need the _classical_ binary
description of the matrix $A\in\mathbb{R}^{m\times n}$ and the vectors
$b\in\mathbb{R}^{m}$, $c\in\mathbb{R}^{n}$.

######  Lemma 5.

Suppose that the total number of non-zeros in $A$ is $\textup{nnz}(A)\geq
m+n$. Then, the quantum oracle $O_{f_{\mu}}$ can be constructed with $\ell$
bits of precision using ${\cal
O}\left(\textup{nnz}(A)\cdot\operatorname{poly}(\ell)\right)$ elementary
quantum gates and at most ${\cal
O}(\textup{nnz}(A)\cdot\operatorname{poly}(\ell))$ ancilla qubits.

###### Proof.

For any $\left|{z}\right\rangle\in\Omega$ and
$\left|{t}\right\rangle\in[0,1]$, the quantum oracle is defined as

$O_{f_{\mu}}\left(\left|{z}\right\rangle\otimes\left|{t}\right\rangle\otimes\left|{0}\right\rangle\right)=\left|{z}\right\rangle\otimes\left|{t}\right\rangle\otimes\left|{f_{\mu}(z,t)}\right\rangle,$

where the function $f_{\mu}(z,t)$ is given by

$f_{\mu}(z,t)=\frac{1}{2}\sum^{\mkern
1.5mu\overline{\mkern-1.5mun}}_{j=1}\left(z_{j}s_{j}-\mu(t)\right)^{2}.$

Recall that vector $s$ is given by $s(z)=Mz+q$. We write $M=(m_{j,k})^{\mkern
1.5mu\overline{\mkern-1.5mun}}_{j,k=1}$, so

$s_{j}=\sum^{\mkern
1.5mu\overline{\mkern-1.5mun}}_{k=1}m_{j,k}z_{k}+q_{j}\quad\forall j\in[\mkern
1.5mu\overline{\mkern-1.5mun}].$

Note that the input $\left|{z}\right\rangle$ and $\left|{t}\right\rangle$ are
in binary representation. To maintain $\ell$ bit of precision,
$\left|{z}\right\rangle$ is represented using $\ell\mkern
1.5mu\overline{\mkern-1.5mun}$ qubits (in which each entry uses $\ell$
qubits), and $\left|{t}\right\rangle$ is represented by $\ell$ qubits.

For $j=1,\dots,\mkern 1.5mu\overline{\mkern-1.5mun}$, we denote
$\tilde{N}_{j}$ as the number of non-zeros in the $j$-th row of $M$. Due to
our definition of $M$ (see Section 2.1.2), we have that $\sum^{\mkern
1.5mu\overline{\mkern-1.5mun}}_{j=1}\tilde{N}_{j}={\cal O}(\textup{nnz}(A))$.
To compute each $s_{j}$ (for $j=1,\dots,\mkern
1.5mu\overline{\mkern-1.5mun})$, we need $\tilde{N}_{j}$ uses of quantum
multipliers and quantum adders, respectively.

Since we have a closed-form formula for the function
$\mu(t)\colon[0,1]\to\mathbb{R}$, we assume the value of $\mu(t)$ can be
computed by a quantum circuit with ${\cal O}(\operatorname{poly}(\ell))$
elementary gates and at most ${\cal O}(\operatorname{poly}(\ell))$ ancilla
qubits, where $\ell$ is the digit length of the floating point number $t$.

Counting the arithmetic operations in the definition of $f_{\mu}$, we find
that the function value of $f_{\mu}(z,t)$ can be computed using

$1+\mkern 1.5mu\overline{\mkern-1.5mun}+\sum^{\mkern
1.5mu\overline{\mkern-1.5mun}}_{j=1}\tilde{N}_{j}={\cal O}(\textup{nnz}(A))$

quantum multipliers and $2\mkern 1.5mu\overline{\mkern-1.5mun}$ quantum
adders. In this process, we also need to query $\mu(t)$ for $\mkern
1.5mu\overline{\mkern-1.5mun}$ times. A quantum adder with input size $\ell$
can be implemented using ${\cal O}(\ell)$ elementary gates and ${\cal O}(1)$
ancilla qubits [VBE96]. A quantum multiplier with input size $\ell$ can be
implemented using ${\cal O}(\ell^{\log(3)})$ elementary gates and ${\cal
O}(\ell)$ ancilla qubits [Gid19]. Therefore, we can construct the oracle
$O_{f_{\mu}}$ using ${\cal
O}\left(\textup{nnz}(A)\cdot\operatorname{poly}(\ell)\right)$ elementary gates
and ${\cal O}\left(\textup{nnz}(A)\cdot\operatorname{poly}(\ell)\right)$
ancilla qubits. ∎

Now, we are ready to give an end-to-end complexity result of our quantum
algorithm.

######  Corollary 1.

Let $\varepsilon\in(0,1)$. Choose $\gamma\in(0,1)$ and $\delta\in(0,1)$.
Define $\kappa({\cal M})$ to be the condition number of the matrix ${\cal
M}=\begin{bmatrix}M&I\end{bmatrix}^{\top}$, which depends only on the LO
problem data $(A,b,c)$. Then, with probability at least $1-\delta$, Algorithm
1 returns a classical vector $z\in\mathbb{R}^{\mkern
1.5mu\overline{\mkern-1.5mun}}_{++}$ such that $(z,s(z))\in{\cal
N}_{2}(\gamma)$ with $d_{2}(z,s(z);\varepsilon)\leq\gamma\varepsilon$.
Moreover, if the LO problem data $(A,b,c)$ is specified in a classical data
structure (e.g., a sparse matrix list, etc.), Algorithm 1 can be implemented
with

${\cal O}\left((m+n)\textup{nnz}(A)\kappa({\cal
M})\frac{1}{\varepsilon}\cdot\textup{polylog}\left(m,n,\kappa({\cal
M}),\frac{1}{\delta}\right)\right)$

elementary gates (up to constant digit precision), where $\textup{nnz}(A)$ is
the number of non-zeros in the matrix $A$.

###### Proof.

Without loss of generality, we assume the number of non-zeros
$\textup{nnz}(A)\geq m+n$. By Lemma 5, the quantum evaluation oracle
$O_{f_{\mu}}$ can be constructed up to fixed digit precision using ${\cal
O}(\textup{nnz}(A))$ elementary gates. Therefore, the overall gate complexity
follows from Theorem 4. ∎

##  4 Recovering polylogarithmic dependence on precision

For both practical and theoretical reasons, algorithms for solving
mathematical optimization problems leverage Iterative Refinement (IR) [Wil63]
through which extended precision solutions to a given optimization problem can
be computed using low-precision oracle calls. In this section, we discuss how
IR can be leveraged in combination with the approach described in Section 3 to
exponentially improve the dependence on precision. This leads to a quantum
algorithm for LO problems that exhibits favorable asymptotic scaling compared
to existing state-of-the-art methodologies.

###  4.1 Exactly solving linear optimization problems

Consider the standard setting in which the input data of the LO problem is
integer. The binary encoding of the LO problem (P)-(D) then requires at most

$L=\sum_{i\in[m]}\sum_{j\in[n]}\left(1+\log_{2}(|A_{ij}|+1)\right)+\sum_{i\in[m]}\left(1+\log_{2}(|b_{i}|+1)\right)+\sum_{j\in[n]}\left(1+\log_{2}(|c_{j}|+1)\right)$
(30)

bits. It is well known (see, e.g., [RTV05, Wri97]) that determining an exact
solution to (P)-(D) requires that the IPM solve the primal-dual pair to
precision $\varepsilon=2^{-L}$. In other words, we terminate precisely when we
obtain an iterate $(z,s(z))$ with normalized duality gap

$\mu\leq 2^{-L}.$

Noting the linear dependence on $\varepsilon^{-1}$ in the complexity result
provided in Theorem 4, this would imply that directly identifying an exact
solution would require

$\widetilde{{\cal O}}_{m,n,\kappa({\cal
M})}\left(\left(m+n\right)\textup{nnz}(A)\kappa({\cal M})2^{L}\right)$

gates in the worst case, where $\textup{nnz}(A)$ is the total number of non-
zero entries found in $A$ and $\kappa({\cal M})$ is a condition number that
depends on the LO problem data $\left(A,b,c\right)$. Though this obviously
indicates an exponential worst-case running time, fortunately, one can recover
an algorithm with polylogarithmic dependence on precision through the use of
IR. This is discussed next.

###  4.2 Iterative refinement for linear optimization

We consider a specialization of the techniques first introduced in [GSW12,
GSW16, GS20] to the definition of IPM optimality.

The IR algorithm for LO is presented in full detail in Algorithm 2. The
algorithm takes as input: (i) the classically specified LO problem data
$\left(A,b,c\right)$, (ii) the fixed precision constant $\xi\in(0,1)$ used for
each of our LO oracle calls and (iii) the final precision $\zeta\in(0,1)$ to
which we seek to solve the primal and dual LO problems (P)-(D). Note that the
only requirements on our precision parameters are that $\xi$ is a positive
constant less than one (e.g., it suffices to fix $\xi=10^{-2}$), and $\zeta$
must be less than or equal to $2^{-L}$, where $L$ is the length of the binary
encoding of the LO problem data $\left(A,b,c\right)$ as defined in equation
(30).

Input: Matrix $A\in\mathbb{R}^{m\times n}$, vectors $b\in\mathbb{R}^{m}$ and
$c\in\mathbb{R}^{n}$, error tolerances $0<\zeta\ll\xi<1$

Output: $\zeta$-precise primal solution to (P)-(D)

Initialize: $x^{(0)}\leftarrow 0\in\mathbb{R}^{n}$, $y^{(0)}\leftarrow
0\in\mathbb{R}^{m}$, $\eta^{(0)}\leftarrow 1,k\leftarrow 0$

while  _$x^{\top}\left(c-A^{\top}y\right) >\zeta$_ do

  1. 1.

Compute residuals

$\hat{b}\leftarrow b-Ax^{(k)},\quad\hat{c}\leftarrow c-A^{\top}y^{(k)}$

  2. 2.

Scale shifted problem data

$(\bar{b},\bar{c})=(\eta^{(k)}\hat{b},\eta^{(k)}\hat{c})$

  3. 3.

$\bar{z}^{(k)}\leftarrow$ solve $(A,\bar{b},\bar{c})$ to $\xi$-optimality
using Algorithm 1

  4. 4.

Compute $(\bar{x}^{(k)},\bar{y}^{(k)})$ from $\bar{z}^{(k)}$

  5. 5.

Update solution:

$x^{(k+1)}\leftarrow x^{(k)}+\frac{1}{\eta^{(k)}}\bar{x}^{(k)},\quad
y^{(k+1)}\leftarrow y^{(k)}+\frac{1}{\eta^{(k)}}\bar{y}^{(k)}$

  6. 6.

Update complementarity gap:

$\hat{\mu}^{(k+1)}=\frac{x^{(k)\top}\left(c-A^{\top}y^{(k)}\right)}{n}$

  7. 7.

Update scaling factor

$\eta^{(k+1)}\leftarrow\frac{1}{\hat{\mu}^{(k+1)}}$

  8. 8.

$k\leftarrow k+1$

end while

Algorithm 2 Iterative Refinement for Linear Optimization

The following result originally appeared in [GSW12, GSW16, GS20] and
establishes that the sequence of iterates produced by Algorithm (2) are
increasingly accurate solutions to the primal and dual LO problem. We state a
version of the original statement which is specialized to the definition of
IPM optimality, and is due to Mohammadisiahroudi et al. [MFWT23].

######  Theorem 5 (Corollary 1 in [GS20], Lemma 6.1 in [MFWT23]).

Suppose we are given a primal and dual feasible LO problem with constraint
matrix $A$. Let $x^{(k)}$, $y^{(k)}$ and $\eta^{(k)}$ be the sequences of
primal-dual solutions and scaling factors produced by Algorithm 2. Then for
all $k$,

  * (a)

$\eta^{(k)}\geq\frac{1}{\xi^{k}}$

  * (b)

The primal-dual pair $\left(x^{(k)},y^{(k)}\right)$ is a $\xi^{k+1}$-optimal
solution to (P)-(D). That is:

$\displaystyle Ax^{(k)}$ $\displaystyle\geq b,\quad x^{(k)}>0,$ $\displaystyle
A^{\top}y^{(k)}$ $\displaystyle\leq c,\quad y^{(k)}>0,$

with

$\frac{x^{(k)\top}\left(c-A^{\top}y^{(k)}\right)}{n}\leq\xi^{k+1}.$

Thus, Algorithm 2 terminates after at most $\log(1/\zeta)$ calls to the
limited-precision LO oracle.

We are now in a position to state our main result.

######  Theorem 6.

Let $\xi=10^{-2}$. Define $\kappa({\cal M})$ to be the condition number of the
matrix ${\cal M}=\begin{bmatrix}M&I\end{bmatrix}^{\top}$, which depends only
on the LO problem data $(A,b,c)$. Let $\textup{nnz}(A)$ denote the number of
non-zeros in the matrix $A$.555Without loss of generality, we assume
$\textup{nnz}(A)\geq m+n$. Then, obtaining an exact solution to (P)-(D) using
Algorithm 2 requires at most

$\widetilde{{\cal O}}_{m,n,\kappa({\cal
M})}\left((m+n)\textup{nnz}(A)\kappa({\cal M})L\right)$

elementary gates and ${\cal O}\left(\textup{nnz}(A)L\right)$ classical
operations.

###### Proof.

First, observe that we can classically load the LO problem data
$\left(A,b,c\right)$ using ${\cal O}(\textup{nnz}(A))$ classical operations.
In each iteration, we make a $\xi$-precise call to Algorithm 1, which by
Corollary 1 requires at most

$\widetilde{{\cal O}}_{m,n,\kappa({\cal
M})}\left(\left(m+n\right)\textup{nnz}(A)\kappa({\cal M})\frac{1}{\xi}\right)$

elementary gates. Since it suffices to perform our LO solves to constant
precision in Algorithm 2, i.e., $\xi=10^{-2}$ is a fixed constant, the term
$\frac{1}{\xi}$ in the complexity simply amounts to constant overhead,
implying that each oracle call in fact requires

$\widetilde{{\cal O}}_{m,n,\kappa({\cal
M})}\left((m+n)\textup{nnz}(A)\kappa({\cal M})\right)$

elementary gates.

From here, preparing the next refining problem requires ${\cal
O}(\textup{nnz}(A))$ classical operations, since these steps amount to
computing matrix-vector products and vector operations. Hence, each iteration
of Algorithm 2 requires at most

$\widetilde{{\cal O}}_{m,n,\kappa({\cal
M})}\left((m+n)\textup{nnz}(A)\kappa({\cal M})\right)$

elementary gates, and ${\cal O}(\textup{nnz}(A))$ classical operations.
Combining these facts with the iteration bound provided in Theorem 5, a
$\zeta$-optimal solution to (P)-(D) can be computed using at most

$\widetilde{{\cal O}}_{m,n,\kappa({\cal
M})}\left(\left(m+n\right)\textup{nnz}(A)\kappa({\cal
M})\log\left(\frac{1}{\zeta}\right)\right)$

elementary gates, and ${\cal
O}\left(\textup{nnz}(A)\log\left(\frac{1}{\zeta}\right)\right)$ classical
operations. Plugging in $\zeta=2^{-L}$, we obtain the stated result. ∎

###  4.3 Comparison to existing LO solvers

In Table 2 we provide a comparison of our algorithms to the current state of
the art algorithms for solving Linear Optimization problems in both the
classical and quantum models of computation. One can observe that the QCPM
outperforms the state of the art classical and quantum IPMs with respect to
$m$ and $n$ when $\textup{nnz}(A)<(m+n)^{\omega-1}$, and achieves an end-to-
end speedup over these algorithms whenever

$\textup{nnz}(A)\kappa({\cal M})<(m+n)^{\omega-1}\cdot\textup{polylog}(m,n).$

Strikingly, these speedups are not reliant on QRAM: our algorithm only
requires access to the natural binary description for the LO problem data
$(A,b,c)$. Like the QIPM in [MFWT23], our algorithm exhibits a dependence on a
condition number that depends only on the problem input data. However, our
condition number dependence is linear, while it is quadratic in [MFWT23]. We
also emphasize that the QIPM from Mohammadisiahroudi et al. [MFWT23] assumes
access to a classical-write/quantum-read RAM (QRAM), and its running time
becomes uncompetitive in the gate model (without QRAM).

There are also classical and quantum algorithms based on zero-sum games and
Gibbs sampling [BGJ+23, vAG19b] that achieve sublinear running times in $m$
and $n$, but exhibit polynomial dependence on precision and the size of the
solution to (P)-(D). These algorithms assume $A\in[-1,1]^{m\times n}$, and
employ different definitions of optimality and feasibility than (Q)IPMs and
our QCPM. The output of (Q)IPMs and the QCPM is a classical primal-dual pair
$(x,y)$ such that

$\displaystyle\left(A_{i,\cdot}\right)^{\top}x$ $\displaystyle\leq
b_{i}\quad\forall i\in[m],\quad x>0,$ $\displaystyle A^{\top}y$
$\displaystyle\leq c,\quad\quad\quad\quad\quad\quad y>0,$

with666To be a bit more technically precise: (Q)IPMs and the QCPM output a
solution satisfying $c^{\top}x-b^{\top}y\leq n\varepsilon$, but these
algorithms exhibit polylogarithmic dependence on precision (provided the
quantum algorithms use iterative refinement), so one can recover duality gap
$\varepsilon$ by taking on additional polylogarithmic overhead in $n$ (i.e.,
simply run the (Q)IPM/QCPM to precision
$\varepsilon^{\prime}=\frac{\varepsilon}{n}.$)

$c^{\top}x-b^{\top}y\leq\varepsilon.$

Conversely, quantum algorithms for zero-sum games output a dual solution
$y\in\mathbb{R}^{m}$ such that the primal-dual pair $(x,y)$ satisfies

$\displaystyle\left(A_{i,\cdot}\right)^{\top}x$ $\displaystyle\leq
b_{i}+\varepsilon_{\text{abs}}\quad\forall i\in[m],\quad x>0,$ $\displaystyle
A^{\top}y$ $\displaystyle\leq
c+\varepsilon_{\text{abs}},\quad\quad\quad\quad\quad\quad y>0,$

and the objective value attained by this solution is
$\text{OPT}\in[\varsigma-{\cal O}(\epsilon_{\text{abs}}),\varsigma+{\cal
O}(\epsilon_{\text{abs}})]$, where $\varsigma$ is a bound on the optimal
objective value determined using binary search. Observe that this is another
difference between the QCPM/(Q)IPMs and quantum algorithms for zero-sum games:
in contrast with the approximate infeasibility of solutions obtained from the
zero-sum games approach, the output of the QCPM/(Q)IPMs always satisfies
primal and dual feasibility exactly. We also point out that, like QIPMs, the
state of the art running times for quantum algorithms for zero-sum games rely
on access to QRAM.

Table 2: Complexity to solve the primal-dual pair (P)-(D) to precision $2^{-L}$ Algorithm | Complexity | QRAM | Notes
---|---|---|---
IPM [CLS21, vdB20] | $\widetilde{{\cal O}}_{n}((m+n)^{\omega}L)$ | - | -
QMMWU [BGJ+23] | $\widetilde{{\cal O}}\left(\sqrt{m+n}r^{2.5}2^{2.5L}+2^{3L}\right)$ | ✓ | $r\geq\|x_{*}\|_{1}$
IR-QIPM [MFWT23] | $\widetilde{{\cal O}}_{n,\kappa(Q)}\left((m+n)^{2.5}L\kappa(Q)^{2}\|Q\|\|x_{*}\|^{5}\right)$ | ✓ |
IR-QCPM (this work) | $\widetilde{{\cal O}}_{m,n,\kappa({\cal M})}\left(\left(m+n\right)\textup{nnz}(A)\kappa({\cal M})L\right)$ | ✗ |

##  5 Conclusion and outlook

In this work, we proposed a new quantum algorithm that solve linear
optimization problems by quantum evolution of the central path. Combining our
approach with iterative refinement techniques, one can obtain an exact
solution to a linear optimization problem involving $m$ constraints and $n$
variables using at most $\widetilde{{\cal O}}_{m,n,\kappa({\cal
M}),\frac{1}{\delta}}\left((m+n)\textup{nnz}(A)\kappa({\cal M})L\right)$
elementary gates and ${\cal O}\left(\textup{nnz}(A)L\right)$ classical
arithmetic operations, where $\textup{nnz}(A)$ is the total number of nonzero
entries found in $A$, $L$ denotes binary input length of the problem data, and
$\kappa({\cal M})$ is the condition number that depends only on the input data
$(A,b,c)$. When the constraint matrix $A$ is sufficiently sparse in the sense
that $\textup{nnz}(A)<(m+n)^{\omega-1}$, our results imply a polynomial
speedup in $m$ and $n$ over all general-purpose classical and quantum
algorithms for solving linear optimization problems that have polylogarithmic
dependence on the error. For problems satisfying $\textup{nnz}(A)\kappa({\cal
M})<(m+n)^{\omega-1}\cdot\textup{polylog}(m,n)$ the QCPM attains an end-to-end
speedup over these algorithms. We stress that these speedups are not reliant
on QRAM, our algorithm only relies on data structures for the sparse binary
representation of $(A,b,c)$.

Though we leave this for future work, our framework should readily generalize
to more complex classes of convex optimization problems, such as semidefinite
optimization and second-order conic optimization. Our work highlights a
previously unexplored connection between Interior Point Methods and the
Quantum Adiabatic Algorithm [FGGS00], and we believe further exploration into
this relationship to be a worthy endeavor.

## Acknowledgements

This project has been carried out thanks to funding by the Defense Advanced
Research Projects Agency (DARPA), ONISQ grant W911NF2010022, titled The
Quantum Computing Revolution and Optimization: Challenges and Opportunities.
G. Nannicini is supported by ONR award # N000142312585. J.L. and X.W. are
partially supported by the U.S. National Science Foundation grant CCF-1816695
and CCF-1942837 (CAREER), and a Sloan research fellowship.

## References

  * [AL22] Dong An and Lin Lin.  Quantum linear system solver based on time-optimal adiabatic quantum computing and quantum approximate optimization algorithm.  ACM Transactions on Quantum Computing, 3(2):1–28, 2022.
  * [Ans88] Kurt M. Anstreicher.  Linear programming and the Newton barrier flow.  Mathematical Programming, 41(1-3):367–373, 1988.
  * [ANTZ23a] Brandon Augustino, Giacomo Nannicini, Tamás Terlaky, and Luis Zuluaga.  Solving the semidefinite relaxation of QUBOs in matrix multiplication time, and faster with a quantum computer, 2023.  arXiv:2301.04237.
  * [ANTZ23b] Brandon Augustino, Giacomo Nannicini, Tamás Terlaky, and Luis F. Zuluaga.  Quantum interior point methods for semidefinite optimization.  Quantum, 7:1110, 2023.
  * [BCLL18] Sébastien Bubeck, Michael B. Cohen, Yin Tat Lee, and Yuanzhi Li.  An homotopy method for $\ell_{p}$-regression provably beyond self-concordance and in input-sparsity time.  In Proceedings of the 50th Annual ACM SIGACT Symposium on Theory of Computing, pages 1130–1137, 2018.
  * [BGJ+23] Adam Bouland, Yosheb M. Getachew, Yujia Jin, Aaron Sidford, and Kevin Tian.  Quantum speedups for zero-sum games via improved dynamic Gibbs sampling.  In International Conference on Machine Learning, pages 2932–2952. PMLR, 2023.
  * [BL89a] Dave A. Bayer and Jeffrey C. Lagarias.  The nonlinear geometry of linear programming. I. Affine and projective scaling trajectories.  Transactions of the American Mathematical Society, 314(2):499–526, 1989.
  * [BL89b] Dave A. Bayer and Jeffrey C. Lagarias.  The nonlinear geometry of linear programming. II. Legendre transform coordinates and central trajectories.  Transactions of the American Mathematical Society, 314(2):527–581, 1989.
  * [Bou99] Jean Bourgain.  On growth of Sobolev norms in linear Schrödinger equations with smooth time dependent potential.  Journal d’Analyse Mathématique, 77(1):315–348, 1999.
  * [BS17] Fernando G.S.L. Brandão and Krysta M. Svore.  Quantum speed-ups for solving semidefinite programs.  In Rafail Ostrovsky and Chris Umans, editors, 2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS), pages 415–426. IEEE, 2017.
  * [CGJ19] Shantanav Chakraborty, András Gilyén, and Stacey Jeffery.  The power of block-encoded matrix powers: improved regression techniques via faster Hamiltonian simulation.  In Christel Baier, Ioannis Chatzigiannakis, Paola Flocchini, and Stefano Leonardi, editors, 46th International Colloquium on Automata, Languages, and Programming (ICALP 2019), volume 132, pages 33:1–33:14, Dagstuhl, Germany, 2019. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik.
  * [CLL+22a] Andrew M. Childs, Jiaqi Leng, Tongyang Li, Jin-Peng Liu, and Chenyi Zhang.  Quantum simulation of real-space dynamics.  Quantum, 6:860, 2022.
  * [CLL+22b] Andrew M Childs, Tongyang Li, Jin-Peng Liu, Chunhao Wang, and Ruizhe Zhang.  Quantum algorithms for sampling log-concave distributions and estimating normalizing constants.  Advances in Neural Information Processing Systems, 35:23205–23217, 2022.
  * [CLS21] Michael B. Cohen, Yin Tat Lee, and Zhao Song.  Solving linear programs in the current matrix multiplication time.  Journal of the ACM (JACM), 68(1):1–39, 2021.
  * [CMSV17] Michael B. Cohen, Aleksander Mądry, Piotr Sankowski, and Adrian Vladu.  Negative-weight shortest paths and unit capacity minimum cost flow in $\widetilde{{\cal O}}(m^{10/7}\log{W})$ time.  In Proceedings of the Twenty-Eighth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 752–771. SIAM, 2017.
  * [CMTV17] Michael B. Cohen, Aleksander Mądry, Dimitris Tsipras, and Adrian Vladu.  Matrix scaling and balancing via box constrained Newton’s method and interior point methods.  In 2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS), pages 902–913. IEEE, 2017.
  * [CW12] Andrew M. Childs and Nathan Wiebe.  Hamiltonian Simulation using linear combinations of unitary operations.  Quantum Information and Computation, 12(11–12):901–924, Nov 2012\.
  * [Dan48] George B. Dantzig.  Programming in a linear structure.  In Bulletin of the American Mathematical Society, volume 54, pages 1074–1074. AMER MATHEMATICAL SOC 201 CHARLES ST, PROVIDENCE, RI 02940-2213, 1948.
  * [DCS+22] Alexander M. Dalzell, B. David Clader, Grant Salton, Mario Berta, Cedric Yen-Yu Lin, David A. Bader, Nikitas Stamatopoulos, Martin J.A. Schuetz, Fernando G.S.L. Brandão, Helmut G. Katzgraber, et al.  End-to-end resource analysis for quantum interior point methods and portfolio optimization, 2022.  arXiv:2211.12489.
  * [Dik67] I.I. Dikin.  Iterative solution of problems of linear and quadratic programming.  In Doklady Akademii Nauk, volume 174, pages 747–748. Russian Academy of Sciences, 1967.
  * [Dik74] I.I. Dikin.  On the convergence of an iteration process.  Upravlyaemya Sistemi, (12):54–60, 1974.
  * [DOW+55] George B. Dantzig, Alex Orden, Philip Wolfe, et al.  The generalized simplex method for minimizing a linear form under linear inequality restraints.  Pacific Journal of Mathematics, 5(2):183–195, 1955.
  * [Fay91] Leonid Faybusovich.  Hamiltonian structure of dynamical systems which solve linear programming problems.  Physica D: Nonlinear Phenomena, 53(2-4):217–232, 1991.
  * [Fay95] Leonid Faybusovich.  A Hamiltonian formalism for optimization problems.  Journal of Mathematical Systems, Estimation, and Control, 5(3):367–370, 1995.
  * [FGGS00] Edward Farhi, Jeffrey Goldstone, Sam Gutmann, and Michael Sipser.  Quantum computation by adiabatic evolution, 2000.  arXiv:quant-ph/0001106.
  * [FM64] Anthony V. Fiacco and Garth P. McCormick.  The sequential unconstrained minimization technique for nonlinear programing, a primal-dual method.  Management Science, 10(2):360–366, 1964.
  * [FM90] Anthony V. Fiacco and Garth P. McCormick.  Nonlinear programming: sequential unconstrained minimization techniques.  SIAM, 1990.  First published in 1968 by Research Analysis Corporation.
  * [Fri54] Ragnar Anton Kittil Frisch.  Principles of linear programming: With particular reference to the double gradient form of the logarithmic potential method.  Universitetets Socialokonomiske Institutt, 1954.
  * [Fri55] Ragnar Anton Kittil Frisch.  The logarithmic potential method of convex programming.  Memorandum, University Institute of Economics, Oslo, 5(6), 1955\.
  * [Fri56] Ragnar Anton Kittil Frisch.  La résolution des problèmes de programme linéaire par la méthode du potentiel logarithmique.  Cahiers du Seminaire D’Econometrie, pages 7–23, 1956.
  * [Gid19] Craig Gidney.  Asymptotically efficient quantum Karatsuba multiplication, 2019.  arXiv:1904.07356.
  * [GS18] David J. Griffiths and Darrell F. Schroeter.  Introduction to quantum mechanics.  Cambridge University Press, 2018.
  * [GS20] Ambros M. Gleixner and Daniel E. Steffy.  Linear programming using limited-precision oracles.  Mathematical Programming, 183(1):525–554, 2020.
  * [GSW12] Ambros M. Gleixner, Daniel E. Steffy, and Kati Wolter.  Improving the accuracy of linear programming solvers with iterative refinement.  In Joris van der Hoeven and Mark van Hoeij, editors, Proceedings of the 37th International Symposium on Symbolic and Algebraic Computation, pages 187–194, 2012.
  * [GSW16] Ambros M. Gleixner, Daniel E. Steffy, and Kati Wolter.  Iterative refinement for linear programming.  INFORMS Journal on Computing, 28(3):449–464, 2016.
  * [HHL09] Aram W. Harrow, Avinatan Hassidim, and Seth Lloyd.  Quantum algorithm for linear systems of equations.  Physical Review Letters, 103(15):150502, 2009.
  * [HJ02] George A. Hagedorn and Alain Joye.  Elementary exponential error estimates for the adiabatic approximation.  Journal of Mathematical Analysis and Applications, 267(1):235–246, 2002.
  * [HKZ12] Daniel Hsu, Sham Kakade, and Tong Zhang.  A tail inequality for quadratic forms of subgaussian random vectors.  Electronic Communications in Probability, 17:1 – 6, 2012.
  * [HS12] Peter D. Hislop and Israel Michael Sigal.  Introduction to spectral theory: With applications to Schrödinger operators, volume 113.  Springer Science & Business Media, 2012.
  * [JKL+20] Haotian Jiang, Tarun Kathuria, Yin Tat Lee, Swati Padmanabhan, and Zhao Song.  A faster interior point method for semidefinite programming.  In Sandy Irani, Lisa O’Conner, and Patrick Kellenberger, editors, 2020 IEEE 61st Annual Symposium on Foundations of Computer Science (FOCS), pages 910–918. IEEE, 2020.
  * [JSWZ20] Shunhua Jiang, Zhao Song, Omri Weinstein, and Hengjie Zhang.  Faster dynamic matrix inverse for faster lps, 2020.  arXiv:2004.07470.
  * [Kar84] Narendra Karmarkar.  A new polynomial-time algorithm for linear programming.  Combinatorica, (4):373–395, 1984.
  * [Kha80] Leonid G. Khachiyan.  Polynomial algorithms in linear programming.  USSR Computational Mathematics and Mathematical Physics, 20(1):53–72, 1980.
  * [KP20] Iordanis Kerenidis and Anupam Prakash.  A quantum interior point method for LPs and SDPs.  ACM Transactions on Quantum Computing, 1(1):1–32, 2020.
  * [LHLW23] Jiaqi Leng, Ethan Hickman, Joseph Li, and Xiaodi Wu.  Quantum Hamiltonian Descent, 2023.  arXiv:2303.01471.
  * [LS14] Yin Tat Lee and Aaron Sidford.  Path finding methods for linear programming: Solving linear programs in $\mathcal{O}(\sqrt{\text{rank}})$ iterations and faster algorithms for maximum flow.  In 2014 IEEE 55th Annual Symposium on Foundations of Computer Science, pages 424–433. IEEE, 2014.
  * [LS15] Yin Tat Lee and Aaron Sidford.  Efficient inverse maintenance and faster algorithms for linear programming.  In 2015 IEEE 56th Annual Symposium on Foundations of Computer Science, pages 230–249. IEEE, 2015.
  * [LZW23] Jiaqi Leng, Yufan Zheng, and Xiaodi Wu.  A quantum-classical performance separation in nonconvex optimization, 2023\.  arXiv:2311.00811.
  * [Meg89] Nimrod Megiddo.  Pathways to the optimal set in linear programming.  Springer, 1989.
  * [Mes58] Albert Messiah.  Quantum Mechanics.  John Wiley & Sons, New York, 1958.
  * [MFT22] Mohammadhossein Mohammadisiahroudi, Ramin Fakhimi, and Tamás Terlaky.  Efficient use of quantum linear system algorithms in interior point methods for linear optimization, 2022.  arXiv:2205.01220.
  * [MFWT23] Mohammadhossein Mohammadisiahroudi, Ramin Fakhimi, Zeguan Wu, and Tamás Terlaky.  An inexact feasible interior point method for linear optimization with high adaptability to quantum computers, 2023.  arXiv:2307.14445.
  * [M1̨3] Aleksander Mądry.  Navigating central path with electrical flows: From flows to matchings, and back.  In 2013 IEEE 54th Annual Symposium on Foundations of Computer Science, pages 253–262. IEEE, 2013.
  * [M1̨6] Aleksander Mądry.  Computing maximum flow with augmenting electrical flows.  In 2016 IEEE 57th Annual Symposium on Foundations of Computer Science (FOCS), pages 593–602. IEEE, 2016.
  * [NN88] Yurii E. Nesterov and Arkadi Nemirovskii.  A general approach to polynomial-time algorithms design for convex programming.  Report, Central Economical and Mathematical Institute, USSR Academy of Sciences, Moscow, 1988.
  * [NN94] Yurii E. Nesterov and Arkadi Nemirovskii.  Interior-Point Polynomial Algorithms in Convex Programming, volume 13.  SIAM, 1994.
  * [Osb12] Tobias J Osborne.  Hamiltonian complexity.  Reports on progress in physics, 75(2):022001, 2012.
  * [Ren88] James Renegar.  A polynomial-time algorithm, based on Newton’s method, for linear programming.  Mathematical Programming, 40(1-3):59–93, 1988.
  * [Ren01] James Renegar.  A mathematical view of interior-point methods in convex optimization.  SIAM, 2001.
  * [RTV05] Cornelis Roos, Tamás Terlaky, and Jean-Phillipe Vial.  Interior Point Methods for Linear Optimization.  Springer Science & Business Media, 2005.
  * [Ter13] Tamás Terlaky.  Interior Point Methods of Mathematical Programming, volume 5.  Springer Science & Business Media, 2013.
  * [TK02] Kim-Chuan Toh and Masakazu Kojima.  Solving some large scale semidefinite programs via the conjugate residual method.  SIAM Journal on Optimization, 12(3):669–691, 2002.
  * [vAG19a] Joran van Apeldoorn and András Gilyén.  Improvements in Quantum SDP-Solving with Applications.  In Christel Baier, Ioannis Chatzigiannakis, Paola Flocchini, and Stefano Leonardi, editors, 46th International Colloquium on Automata, Languages, and Programming (ICALP 2019), volume 132 of Leibniz International Proceedings in Informatics (LIPIcs), pages 99:1–99:15, Dagstuhl, Germany, 2019. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik.
  * [vAG19b] Joran van Apeldoorn and András Gilyén.  Quantum algorithms for zero-sum games, 2019.  arXiv:1904.03180.
  * [vAGGdW20] Joran van Apeldoorn, András Gilyén, Sander Gribling, and Ronald de Wolf.  Quantum SDP-solvers: Better upper and lower bounds.  Quantum, 4:230, 2020.
  * [VBE96] Vlatko Vedral, Adriano Barenco, and Artur Ekert.  Quantum networks for elementary arithmetic operations.  Physical Review A, 54(1):147, 1996.
  * [vdB20] Jan van den Brand.  A deterministic linear program solver in current matrix multiplication time.  In Proceedings of the Fourteenth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 259–278. SIAM, 2020.
  * [vdBLSS20] Jan van den Brand, Yin Tat Lee, Aaron Sidford, and Zhao Song.  Solving tall dense linear programs in nearly linear time.  In Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing, pages 775–788, 2020.
  * [Vla23] Adrian Vladu.  Interior point methods with a gradient oracle, 2023.  arXiv:2304.04550.
  * [Wil63] James Hardy Wilkinson.  Rounding Errors in Algebraic Processes.  London: Her Majesty’s Stationary Office; Eaglewoocl Cliffs, N.J.: Prentice-Hall, 1963.
  * [Wri97] Stephen J. Wright.  Primal-Dual Interior-Point Methods, volume 54.  SIAM, 1997.
  * [WWJ16] Andre Wibisono, Ashia C. Wilson, and Michael I. Jordan.  A variational perspective on accelerated methods in optimization.  Proceedings of the National Academy of Sciences, 113(47):E7351–E7358, 2016.
  * [Ye11] Yinyu Ye.  Interior point algorithms: theory and analysis.  John Wiley & Sons, 2011.
  * [YTM94] Yinyu Ye, Michael J. Todd, and Shinji Mizuno.  An $\mathcal{O}(\sqrt{n}L)$-iteration homogeneous and self-dual linear programming algorithm.  Mathematics of Operations Research, 19(1):53–67, 1994.

Appendices

##  Appendix A Proof of Proposition 1

######  Proposition 3 (Proposition 1, [HKZ12]).

Let $V\in\mathbb{R}^{m\times n}$ be a matrix, and let $\Sigma=V^{T}V$. Let
$x=(x_{1},\dots,x_{n})$ be an isotropic multivariate Gaussian random vector
with mean zero. For all $t>0$,

$\displaystyle\Pr\left[\|Vx\|^{2}>\operatorname{Tr}(\Sigma)+2\sqrt{\operatorname{Tr}(\Sigma^{2})t}+2\|\Sigma\|t\right]\leq
e^{-t}.$

Now, we are ready to prove Proposition 1.

###### Proof.

Recall that we define $f(x;\mu)=\frac{1}{2}\|F(x)-\mu
e\|^{2}=\frac{1}{2}d^{2}_{2}(x;\mu)$. So $x\in{\cal N}_{2}(\gamma)$ is
equivalent to

$f(x)\leq\frac{\gamma^{2}\mu^{2}}{2}.$

Let $W\coloneqq\frac{h}{2}(H(\mu))^{-1/2}$. We define $y\coloneqq
W^{-1/2}\left(x-z(\mu)\right)$, then for $x\sim|\Phi_{0}(\mu)|^{2}$, we have
$y\sim{\cal N}(0,\mathbb{I})$. In a small neighborhood of $z(\mu)$,

$f_{\mu}(x)\approx\frac{1}{2}(x-z(\mu))^{T}H(\mu)(x-z(\mu))=\frac{1}{2}y^{T}(H(\mu)W)y=\|Vy\|^{2},$

where $V\coloneqq\frac{1}{2}\sqrt{h}(H(\mu))^{1/4}$. Let $\Sigma\coloneqq
V^{2}=\frac{h}{4}(H(\mu))^{1/2}$, we have that

$\displaystyle\operatorname{Tr}[\Sigma]=\frac{h}{4}\left(\sum^{\mkern
1.5mu\overline{\mkern-1.5mun}}_{k=1}\lambda^{1/2}_{k}\right)\leq\frac{\sqrt{\eta_{2}\kappa}\mkern
1.5mu\overline{\mkern-1.5mun}h\mu}{4},\quad\operatorname{Tr}[\Sigma^{2}]=\frac{h^{2}}{16}\left(\sum^{\mkern
1.5mu\overline{\mkern-1.5mun}}_{k=1}\lambda_{k}\right)\leq\frac{\eta_{2}\kappa\mkern
1.5mu\overline{\mkern-1.5mun}h^{2}\mu^{2}}{16},$

in which we use the estimate
$\lambda_{k}\leq\lambda_{\max}\leq\kappa\lambda_{\min}\leq\kappa\eta_{2}\mu^{2}$.

Therefore, by Proposition 3,

$\Pr_{x\sim|\Phi_{0}(\mu)|^{2}}\left[f(x;\mu)>\sqrt{\eta_{2}\kappa}h\mu\left(\mkern
1.5mu\overline{\mkern-1.5mun}/4+\sqrt{\mkern
1.5mu\overline{\mkern-1.5mun}\log(1/\delta)}/2+\log(1/\delta)/2\right)\right]\leq\delta.$

Due to our choice of $h$, it follows that

$\Pr_{x\sim|\Phi_{0}(\mu)|^{2}}\left[f_{\mu}(x)>\frac{1}{2}\gamma^{2}\mu^{2}\right]\leq\delta.$

∎

##  Appendix B Additional Proofs

######  Lemma 6.

Let $\Psi(x,t)$ denote the exact solution of (13) and $\widetilde{\Psi}(x,t)$
denote the approximated solution by the Fourier spectral method (truncated up
to frequency $n$).777More details on the Fourier spectral method can be found
in Section 2.2 (in particular Lemma 1) in [CLL+22b]. We assume that the
initial data $\Psi_{0}(x)$ is periodic and analytic in $x\in\Omega$. Then, for
any integer $n\geq 1$, the error from the Fourier spectral method satisfies

$\displaystyle\max_{x,t}|\Psi(x,t)-\widetilde{\Psi}(x,t)|\leq 2r^{n/2+1},$
(31)

where $0<r<\min(\frac{1}{2},\frac{1}{At})$, $A$ is an absolute constant that
only depends on $\Psi_{0}(x)$.

###### Proof.

We give the proof in one dimension, as the same argument is readily
generalized to arbitrary finite dimensions. We assume the initial data
$\Psi_{0}(x)$ is periodic over $[0,2\pi]$. The analyticity implies that, for
any $x\in[0,2\pi]$, there is a constant $C$ such that

$\displaystyle\left|\Psi^{(k)}_{0}(x)\right|\leq C^{k+1}(k!).$ (32)

Therefore, the function $\Psi_{0}(x)$ admits an analytic continuation in the
strip

$\Gamma_{0}=\left\\{z\in\mathbb{C}^{:}|z-x|<\frac{1}{C},x\in[0,2\pi]\right\\}.$

Due to [Bou99, Proposition 1], there is an absolute constant $A$ such that for
any $t\in[0,T]$,

$\displaystyle\left\|\Psi^{(k)}(\cdot,t)\right\|\leq
At\left\|\Psi^{(k)}_{0}\right\|\leq AC^{k+1}t(k!),$ (33)

which implies that the wave function $\Psi(x,t)$ is periodic and analytic for
any finite $t$. The strip on which $\Psi(x,t)$ admits an analytic continuation
is

$\Gamma_{t}=\left\\{z\in\mathbb{C}^{:}|z-x|<\frac{1}{ACt},x\in[0,2\pi]\right\\}.$

Suppose that the function $\Psi(x,t)$ allows an exact, infinite trigonometric
polynomial representation (see [CLL+22b, Lemma 16]),

$\displaystyle\Psi(x,t)=\sum^{\infty}_{k=0}c_{k}(t)e^{ikx}.$ (34)

Let $\tilde{\Psi}(t,x)$ be the truncated Fourier series up to $k=n/2$, then
the error from the Fourier spectral method satisfies

$\displaystyle|\Psi(x,t)-\tilde{\Psi}(x,t)|\leq\sum^{\infty}_{k=n/2+1}|c_{k}(t)|.$
(35)

Meanwhile, the function $\Psi(x,t)$ has an analytic continuation defined by
$\Psi(x,t)=\sum^{\infty}_{k=0}c_{k}(t)z^{k}$ for $z\in\Gamma_{t}$. By Cauchy’s
integral formula, for any simply connected curve $\gamma$ on the strip
$\Gamma_{t}$, we have that

$\displaystyle c_{k}(t)=\frac{1}{k!}\int_{\gamma}\Psi^{(k)}(z,t)\mathrm{d}z.$
(36)

Together with (33), it turns out that $|c_{k}(t)|\leq r^{k}$ for some
$0<r<\frac{1}{At}$. Therefore, the error from the Fourier spectral method
satisfies

$\displaystyle|\Psi(x,t)-\tilde{\Psi}(x,t)|\leq\sum^{\infty}_{k=n/2+1}|c_{k}(t)|\leq\frac{r^{n/2+1}}{1-r}.$
(37)

Moreover, if we force $r<1/2$, the error is bounded by $2r^{n/2+1}$. ∎

##### Proof of Theorem 3.

To ensure that the simulation error is bounded by $\eta$, we may choose the
truncation number

$n=\left\lceil
2\left(\frac{\log(4/\eta)}{\log(1/r)}-1\right)\right\rceil\leq{\cal
O}\left(\log(1/\eta)\right).$

Next, by plugging this truncation number in Equation (113) in [CLL+22b], we
prove Theorem 3.

######  Lemma 7.

Suppose that $\psi_{1}$, $\psi_{2}$ are two unit vectors such that
$\|\psi_{1}-\psi_{2}\|\leq\delta$. Then, we have

$\left|\left\langle{\psi_{1}}\right|O\left|{\psi_{1}}\right\rangle-\left\langle{\psi_{2}}\right|O\left|{\psi_{2}}\right\rangle\right|\leq
2\|O\|\delta.$

###### Proof.

By the triangle inequality,

$\displaystyle\left|\left\langle{\psi_{1}}\right|O\left|{\psi_{1}}\right\rangle-\left\langle{\psi_{2}}\right|O\left|{\psi_{2}}\right\rangle\right|\leq\left|(\left\langle{\psi_{1}}\right|-\left\langle{\psi_{2}}\right|)O\left|{\psi_{1}}\right\rangle\right|+\left|\left\langle{\psi_{2}}\right|O(\left|{\psi_{1}}\right\rangle-\left|{\psi_{2}}\right\rangle)\right|\leq
2\|O\|\delta.$

∎