∎

11institutetext: Mohammadhossein Mohammadisiahroudi, Corresponding author
22institutetext: Lehigh University
Bethlehem, PA 18015, USA
<EMAIL_ADDRESS>33institutetext: Ramin Fakhimi 44institutetext: Lehigh
University
Bethlehem, PA 18015, USA
<EMAIL_ADDRESS>55institutetext: Tamás Terlaky 66institutetext: Lehigh
University
Bethlehem, PA 18015, USA
<EMAIL_ADDRESS>
# Efficient Use of Quantum Linear System Algorithms in Interior Point Methods
for Linear Optimization

Mohammadhossein Mohammadisiahroudi  Ramin Fakhimi  and Tamás Terlaky

(Received: date / Accepted: date)

###### Abstract

Quantum computing has attracted significant interest in the optimization
community because it potentially can solve classes of optimization problems
faster than conventional supercomputers. Several researchers proposed quantum
computing methods, especially Quantum Interior Point Methods (QIPMs), to solve
convex optimization problems, such as Linear Optimization, Semidefinite
Optimization, and Second-order Cone Optimization problems. Most of them have
applied a Quantum Linear System Algorithm at each iteration to compute a
Newton step. However, using quantum linear solvers in QIPMs comes with many
challenges, such as having ill-conditioned systems and the considerable error
of quantum solvers. This paper investigates how one can efficiently use
quantum linear solvers in QIPMs. Accordingly, an Inexact Infeasible Quantum
Interior Point Method is developed to solve linear optimization problems. We
also discuss how we can get an exact solution by Iterative Refinement without
excessive time of quantum solvers. Finally, computational results with a
QISKIT implementation of our QIPM using quantum simulators are analyzed.

###### Keywords:

Quantum Interior Point Method Linear Optimization Quantum Linear System
Algorithm Iterative Refinement

###### MSC:

90C51 90C05 81P68

††journal: JOTA

##  1 Introduction

Linear Optimization (LO) is defined as optimizing a linear function over a set
of linear equality and inequality constraints. Several algorithms were
developed to solve LO problems dantzig2016linear (16, 29). Karmarkar1984_New
(25) developed the foundations of polynomial time and practically efficient
Interior Point Methods (IPMs) for solving LO problems. Since Karmarkar’s
publication, a large class of theoretically and practically efficient IPMs
were developed, see, e.g., Terlaky2013_Interior (40, 36, 43). Contrary to the
Simplex method, a feasible IPM reaches an optimal solution by traversing
through the interior of the feasible region Roos2005_Interior (36).

Contemporary IPMs reach an optimal solution by starting from an interior point
and following the central path Roos2005_Interior (36). The most efficient IPMs
are primal-dual methods, meaning they strive to satisfy the optimality
conditions while maintaining strict primal and dual feasibility. It should be
noted that basic IPMs need an initial feasible interior point. Some current
commercial solvers apply Feasible IPMs (F-IPMs) based on the self-dual
embedding formulation of LO problems, e.g., MOSEK, while Infeasible Interior
Point Methods (I-IPMs) can start with an infeasible but positive solution.
Theoretical analysis shows the best iteration complexity of F-IPMs is ${\cal
O}(\sqrt{n}L)$, where $n$ is the number of variables, and $L$ is the binary
length of the input data. On the other hand, the best iteration complexity of
I-IPMs is ${\cal O}(nL)$. In practice, the performance of both feasible and
infeasible IPM are similar Wright1997_Primal (43).

A linear equation system is solved at each iteration of IPMs to calculate a
Newton direction. There are three choices for the linear equation system: (i)
Full Newton System, (ii) Augmented System, and (iii) Normal Equation System
(NES). In classical computers, a general approach is applying Cholesky
Factorization to solve the NES, because it has a symmetric positive definite
coefficient matrix. A partial update technique improved the complexity of
solving the system in each iteration. This approach leads to the best total
complexity of ${\cal O}(n^{3}L)$ arithmetic operations for solving LO
problems. However, this is not always efficient in practice Roos2005_Interior
(36). Since Cholesky Factorization requires $\mathcal{O}(n^{3})$ arithmetic
operations for large dense matrices, several researchers studied inexact
solution methodologies for solving Newton systems. Bellavia1998_Inexact (9)
proved the convergence of an Inexact I-IPM (II-IPM) for general convex
optimization problems. Mizuno and his colleagues studied the convergence of
II-IPMs MizunoJarre1999_Global (33, 20). Korzak2000_Convergence (31) and
Baryamureeba2006_convergence (8) proved the convergence of II-IPM proposed by
Kojima1993_primal (30). Korzak2000_Convergence (31) also showed that the total
time complexity of his algorithm is polynomial. Gondzio2009_Convergence (24),
and Monteiro2003_Convergence (35) investigated the use of Preconditioned
Conjugate Gradient (PCG) methods in II-IPMs. Zhou2004_Polynomiality (46) also
proved the convergence of an II-IPM for SDO problems. The best iteration
complexity of II-IPMs is ${\cal O}(n^{2}L)$ which is $\Theta(n^{1.5})$ weaker
than the best iteration complexity of exact F-IPMs.

Quantum computers have recently emerged as a powerful alternative to classic
computers Childs2017_Quantum (15). Starting from Deutsch’s Problem, a series
of problems and algorithms have demonstrated theoretically exponential speedup
compared to their classical counterparts Deutsch1985 (18, 19, 39). One of the
promising quantum algorithms is the HHL method Harrow2009_Quantum (23) to
solve linear system problems. The HHL method showed quantum advantage with
respect to the dimension compared to classical linear systems solutions
algorithm. However, this method has unfavorable dependence on the condition
number, the sparsity of the matrix, and the inverse precision of the solution.
Several researchers attempted to improve the performance of Quantum Linear
System Algorithms (QLSAs), Wossnig2018_Quantum (42, 15). This paper explores
an efficient use of the QLSAs for solving the Newton system at each iteration
of IPMs.

To investigate quantum speedup for continuous optimization, Brandão and Svore
brandao2017quantumFOCS (11) proposed a non-IPM quantum algorithm based on the
Multiplicative Weight Update Method (MWUM) of arora2012multiplicative (6) to
solve Semidefinite Optimization (SDO) problems. After this papers, many
improved versions of Quantum MWUMs (QMWUMs) are proposed for SDO
brandao2017quantum (12, 4), and LO van2019quantum (5). The advantage of QMWUMs
is that their complexity has linear dependence on the dimension for SDO and
sublinear dependence for LO. The major issue with QMWUMs is that they are
highly dependent on inverse precision and an upper bound for the norm of the
optimal solution, which are exponentially large for LO.

In another direction, a few studies were proposing Quantum IPMs (QIPMs).
Kerenidis and his colleagues presented a series of papers on QIPMs for solving
LO, SDO kerenidisParkas2020_quantum (26), and Second-order Cone Optimization
(SOCO) problems Kerenidis2021_Quantum (28). Note that for both SDO and SOCO,
they missed some important parts, such as the need for symmetrization/scaling.
So, their QIPMs are invalid for their conic optimization problems.
kerenidisParkas2020_quantum (26) claimed that their algorithm has ${\cal
O}(\frac{n^{2}}{\epsilon^{2}}\kappa^{3}\log(\frac{1}{\zeta}))$ complexity for
LO problems where $\zeta$ is the final optimality gap, $\epsilon$ is the
precision of the QLSA. Observe that, to reach final precision $\zeta$, the
Newton system needs to be solved with precision $\epsilon<\zeta$. Further,
$\kappa$ is an upper bound for the condition number of the Newton systems at
each iteration. They used block encoding to construct the Newton system and
Quantum Tomography Algorithm (QTA) to extract the classical solution.
Casares2020_quantum (13) also provided a hybrid predictor-corrector QIPM
scheme to solve LO problems using the well-known predictor-corrector method
proposed by Ye1994_iteration (45). In each step, it uses the QLSA proposed by
Chakraborty2018_power (14) to solve the Newton systems. They claimed that the
complexity of their method is ${\cal
O}(L\sqrt{n}(n+m)\overline{\|M\|}_{F}\kappa^{2}\epsilon^{-2})$, where $m$ is
the number of constraints, and $\overline{\|M\|}_{F}$ is an upper bound to the
Frobenius norm of the coefficient matrix in the Newton system at each
iteration. The authors claimed a quantum speedup with respect to the dimension
$n$ compared to classical algorithms. However, they overlooked some elements
in their time complexity. Thus, the actual complexity has higher dependence on
$n$ and the required precision.

In both mentioned QIPMs, an exact F-IPM framework was used, regardless of the
inherent inexactness of QLSAs. So, the proposed time complexities are not
attainable, since Inexact IPMs have higher iteration complexity. Augustino et
al. augustino2021quantum (7) addressed this issue and proposed two convergent
QIPMs for SDO. Their complexity shows polynomial speed-up with respect to
dimension, but similar to former QIPMs, their complexity suffers from linear
dependence on condition number of the Newton system and precision. A major
reason is that QLSAs’ time complexities depend on the condition number of the
linear system, and it is linear in inverse precision and dimension to extract
an accurate classical solution. In IPMs, the condition number of the Newton
system typically goes to infinity as the algorithm approaches an optimal
solution Roos2005_Interior (36). It is also worth mentioning that $\epsilon$,
the required precision for solution of the Newton system, needs to be
significantly smaller than the final precision $\zeta$. Consequently, these
complexity bounds are not polynomial in the classical sense. This paper
explores an efficient use of the QLSAs for solving the Newton system at each
iteration of QIPMs. We propose an II-QIPM to find an exact solution. We also
employ an iterative refinement scheme to avoid exponential complexity in
finding an exact optimal solution.

This paper is structured as follows. Section 2 discusses the performance of
existing QLSAs. Section 3 introduces the LO problem and its characteristics.
In Section 4, we present an II-QIPM to solve LO problems, in which a QLSA is
used for solving the NES. We also analyze the complexity of the proposed II-
QIPM. In Section 5, we employ an Iterative Refinement scheme to find an exact
optimal solution of a LO problem without excessive time of QLSAs. Section 6
presents the first implementation of the proposed II-QIPM with Iterative
Refinement (IR-II-QIPM) and evaluates the algorithms through computational
experiments. Finally, conclusions are presented in Section 7.

##  2 Quantum Linear Algebra

This section reviews the use of Quantum algorithms to solve Linear System
Problems (LSP).

######  Definition 1 (LSP).

The Linear System Problem: Find a vector $z\in\mathbb{R}^{p}$ such that it
satisfies equation $Mz=\sigma$ with coefficient matrix
$M\in\mathbb{R}^{p\times p}$ and right-hand side (RHS) vector
$\sigma\in\mathbb{R}^{p}$.

In the rest of the paper, $\|M\|=\|M\|_{2}$ is the 2-norm of matrix $M$, and
$\|M\|_{F}$ is the Frobenius norm of $M$. We sometimes use $\tilde{{\cal O}}$
which suppresses the polylogarithmic factors in the ”Big-O” notation. The
quantities of the polylogarithmic factors are indicated as subscripts of
$\tilde{{\cal O}}$. We use $\mathbb{R}^{n}$ for the set of $n$-dimensional
vectors of real numbers and $\mathbb{C}^{n}$ for the set of $n$-dimensional
vectors of complex numbers. It should be mentioned that the complexity of a
classical algorithm when solving LSP means the number of arithmetic operations
and the complexity of a quantum algorithm is the number of quantum gates.

The LSP can have either one, many, or no solutions. A basic approach for
solving an LSP is Gaussian elimination, or LU factorization, with ${\cal
O}(p^{3})$ arithmetic operations. If $M$ is a square symmetric positive
semidefinite (PSD) matrix, we can also apply Cholesky factorization with
${\cal O}(p^{3})$ arithmetic operations. The best complexity for an iterative
algorithm with respect to $p$ is ${\cal O}(pd\sqrt{\kappa}\log(1/\epsilon))$
arithmetic operations for the Conjugate Gradient method solving systems with
symetric PSD matrices, where $d$ is the maximum number of non-zero elements in
any row or column of $M$, $\kappa$ is the condition number of $M$, and
$\epsilon$ is the error allowed. If matrix $M$ is just symmetric, one can use
Lanczos algorithm with higher complexity. For an LSP with general square
matrix $M$, the best iterative method is the GMRES algorithm, which has ${\cal
O}(n^{3})$ worst-case complexity saad2003iterative (37). To sum up, LSPs in
general form are solvable in polynomial time in classical computing setting,
but the worst-case complexity is ${\cal O}(n^{3})$ for either matrix
decomposition methods or iterative methods.

Before discussing QLSAs, we should mention that the $\ket{z}$ notation
represent the quantum state corresponding to the unit classical vector $z$. We
denote the basis state $\ket{i}$, which is a column vector with dimension $p$,
one in coordinate $i$ and zero in other coordinates Childs2017_Quantum (15).
QLSAs have different approaches, while all of them are solving Quantum Linear
System Problems (QLSPs) defined as follows.

######  Definition 2 (QLSP).

Let $M\in\mathbb{C}^{p\times p}$ be a Hermitian matrix with $\|M\|=1$,
$\sigma\in\mathbb{C}^{p}$, and $z\mathrel{\mathop{\mathchar
58\relax}}=M^{-1}\sigma$. We define quantum states

$\ket{\sigma}=\frac{\sum_{i=1}^{p}\sigma_{i}\ket{i}}{\|\sum_{i=1}^{p}\sigma_{i}\ket{i}\|}\quad\text{
and
}\quad\ket{z}=\frac{\sum_{i=1}^{p}z_{i}\ket{i}}{\|\sum_{i=1}^{p}z_{i}\ket{i}\|}.$

For target precision $0<\epsilon_{QLSP}$, the goal is to find
$\ket{\tilde{z}}$ such that $\|\ket{\tilde{z}}-\ket{z}\|\leq\epsilon_{QLSP}$,
succeeding with probability $\Omega(1)$.

Based on Definition 1 and Definition 2, QLSP is a different form of the LSP.
It should be noted that at each iteration of QIPMs, instead of an LSP, we need
to use a QLSA to solve a QLSP. Thus, we need to translate LSP to QLSP, solve
the QLSP by QLSA and extract the solution by a QTA. Here, we analyze the
details and costs of the process of translating LSPs to QLSPs, encoding in the
quantum setting, solving them with a QLSA, and extracting classical solutions
with a QTA as follows.

  1. (i)

Model of Computation: the first important step is determining how to encode
the input data in the quantum setting. There are two major input models. One
is the sparse-access model which is used in the HHL algorithm
Harrow2009_Quantum (23) and then in other QLSAs Childs2017_Quantum (15, 41).
This is a quantum version of classical sparse matrix computation, and we
assume access to unitaries that calculate the index of the $l$th non-zero
element of the $k$th row of a matrix $M$ when given $(k,l)$ as input. A
different input model, now known as the quantum operator input model, is
proposed in low2019hamiltonian (32), which is based on the idea of block-
encoded matrices. In this input model, one has access to unitaries that store
the coefficient matrix:

$U=\begin{pmatrix}M/\alpha&\cdot\\\ \cdot&\cdot\end{pmatrix},$

where $\alpha\geq\|M\|$ is a normalization factor chosen to ensure that $U$
has norm at most 1. Chakraborty2018_power (14) showed that this quantum
operator input model is more efficient than the sparse-access model and
oracles to encode input data using block-encoding has favorable complexity
compared to the sparse-access model. On the other hand, most of the block-
encoding approaches use Quantum Random Access Memory (QRAM). However, the
sparse-access model can be implemented in the standard gate-based quantum
circuit model. Despite efficient encoding procedures, the quantum operator
input model can not be implemented with current quantum computers since there
is no physical implementation of QRAM. In our analysis, we assume that the
data is stored in QRAM, and we use the quantum operator model by
Chakraborty2018_power (14) for the QLSA, which enjoys the best complexity to
date. Using the QRAM structure, kerenidis2016quantum (27) showed that one can
implement $\epsilon$-approximate block-encoding of $M$ with ${\cal
O}(\text{polylog}(\frac{p}{\epsilon}))$ complexity. Further, given $M$ in the
sparse-access input model, there is an $\epsilon$-approximate block-encoding
of $M$ that can be implemented in complexity ${\cal
O}(\text{polylog}(\frac{p\alpha}{\epsilon}))$ Chakraborty2018_power (14).
Thus, our results also apply to the sparse-access input model if we have the
data in that form. Using these results and assuming access to QRAM,
Chakraborty2018_power (14) proposed a QLSA, in which they construct state
$\ket{\sigma}$, build and implement a block-encoding of matrix $M$ with ${\cal
O}(\text{polylog}(\frac{p}{\epsilon}))$ complexity. Since we are interested in
using QLSAs in a hybrid approach, we need to consider the cost of storing data
in a classical form to QRAM which is ${\cal O}(p^{2})$ for a fully dense
matrix Chakraborty2018_power (14). This cost should be paid once, and in the
next section, it is shown that this cost will be dominated by the classical
operations in each iteration of QIPMs.

  2. (ii)

Translating LSP to QLSP: Based on the definition of QLSP, the coefficient
matrix of the system must be Hermitian. If $M$ is not Hermitian, one can
construct $\bar{M}\bar{z}=\bar{\sigma}$, where

$\displaystyle\bar{M}=\begin{bmatrix}0&M\\\
M^{\dagger}&0\end{bmatrix},\quad\bar{\sigma}=\begin{pmatrix}\sigma\\\
0\end{pmatrix},$

and find the vector $\bar{z}=\begin{pmatrix}0\\\ z\end{pmatrix}$, where
$M^{\dagger}$ denotes the conjugate transpose of $M$. The size of the problem
increases from $p$ to $2p$. QLSP assumes $\|M\|=1$. In the structure of block
encoding, we address this normalization, but for sparse encoding, we need to
normalize the system $\bar{M}\bar{z}=\bar{\sigma}$ where

$\displaystyle\bar{M}=\frac{M}{\|M\|},\quad\bar{\sigma}=\frac{\sigma}{\|M\|},\quad\text{and}\quad\bar{z}=z.$

Although we are using block-encoding, we scale the matrix in advance and let
$\alpha=1$ in the block encoding. This scaling will affect the precision and
complexity, but in this way, we can have better complexity analysis since, in
QIPMs, we need to bound the residual. We also have a similar scaling in the
definitions of states $\ket{z}$ and $\ket{\sigma}$. The corresponding QLSP
problem is ${\overline{M\mkern
4.0mu}\mkern-4.0mu}{}\ket{z}=\ket{\bar{\sigma}}$, and we can find an inexact
solution $\ket{\tilde{z}}$ with
$\|\ket{\tilde{z}}-\ket{z}\|\leq\epsilon_{QLSP}$. Since
$\|\ket{\tilde{z}}\|=1$, we need to scale back the solution with
$\tilde{z}=\|\bar{\sigma}\|\ket{\tilde{z}}$. This scaling affects the target
precision. In IPMs, we are interested in finding a solution where the residual
is bounded as $\|\sigma-M\tilde{z}\|\leq\epsilon_{LSP}$. Thus, we have

$\displaystyle\|\sigma-M\tilde{z}\|$
$\displaystyle\leq\|M\|\|\bar{\sigma}-{\overline{M\mkern
4.0mu}\mkern-4.0mu}{}\tilde{z}\|$
$\displaystyle\leq\|M\|\|\bar{\sigma}\|\|\ket{\bar{\sigma}}-{\overline{M\mkern
4.0mu}\mkern-4.0mu}{}\ket{\tilde{z}}\|$
$\displaystyle\leq\|M\|\frac{\|\sigma\|}{\|M\|}\|{\overline{M\mkern
4.0mu}\mkern-4.0mu}{}\ket{z}-{\overline{M\mkern
4.0mu}\mkern-4.0mu}{}\ket{\tilde{z}}\|$
$\displaystyle\leq\|\sigma\|\|\ket{z}-\ket{\tilde{z}}\|$
$\displaystyle\leq\|\sigma\|\epsilon_{QLSP}.$

Thus, we must set the target error of the QLSP as
$\epsilon_{QLSP}=\epsilon_{LSP}/\|\sigma\|.$

  3. (iii)

QLSA: After preprocessing and encoding, we can apply the QLSA to solve the
QLSP. The HHL algorithm Harrow2009_Quantum (23) solves QLSPs with
$\tilde{{\cal O}}_{p}(\frac{d^{2}\kappa^{2}}{\epsilon_{QLSP}})$ complexity.
Several researchers attempted to improve the performance of the HHL algorithm.
As the first attempt, Amplitude Amplification decreases the dependence on
$\kappa^{2}$ to $\kappa$ ambainis2012variable (1). Wossnig2018_Quantum (42)
proposed a QLSA with $\tilde{{\cal
O}}_{p}(\|M\|_{F}\frac{\kappa}{\epsilon_{QLSP}})$ complexity by using the
Quantum Singular Value Estimation. In another direction, Childs2017_Quantum
(15) developed two QLSAs with exponentially better dependence on error with
$\tilde{{\cal O}}_{p,\kappa,\frac{1}{\epsilon_{QLSP}}}(d\kappa)$ complexity.
They proposed two approaches using Fourier and Chebyshev series
representations. The best QLSA with respect to complexity uses block encoding
and QRAM with $\tilde{{\cal O}}_{p,\frac{1}{\epsilon_{QLSP}}}(\alpha\kappa)$
complexity Chakraborty2018_power (14). The normalization factor $\alpha$ is
equal to $1$ for our approach since we apply the QLSA to QLSP, which is
normalized in advance. The details of these methods are out of the scope of
this paper. For further details, see Dervovic2018_Quantum (17).

  4. (iv)

QTA: QLSAs provide a quantum state proportional to the solution. We cannot
extract the classical solution by a single measurement. We need Quantum
Tomography Algorithms (QTAs) to extract the classical solution. There are
several papers improving QTAs, see, e.g., Kerenidis2021_Quantum (28). We used
the best QTA by van2022quantum (3), with complexity ${\cal
O}(\frac{p\varrho}{\epsilon_{QTA}})$, where $\varrho$ is bound for the norm of
the solution. In our approach, $\varrho$ is equal to $1$ since we applied QTA
to QLSP and scaled back the solution after QTA. Since the error is additive,
we may choose

$\epsilon_{QTA}=\epsilon_{QLSA}=\frac{\epsilon_{LSP}}{2\|\sigma\|}.$

Table 1 presents the complexity of different classical and quantum algorithms
for solving an LSP. Here, the complexity of a QTA is considered in the
complexity of QLSAs. As shown in Table 1, although the complexity of solving
LSP using QLSA+QTA will have a similar dependence on $p$ compared to Conjugate
Gradient (CG), the complexity of CG method is valid for solving LSPs with PSD
matrices and the complexity of QLSA+QTA in Table 1 are valid for solving LSPs
in a general form. QLSAs have better dependence on dimension $p$ compared to
factorization and elimination techniques. Generally, QLSA has worse dependence
on $\kappa$, $\frac{1}{\epsilon_{LSP}}$, $\|M\|$, and $\|\sigma\|$. In our
QIPM, we solve a modified NES which has symmetric PSD matrix and both CG and
QLSA+QTA are applicable. Although CG has better complexity than QLSA+QTA based
on Table 1, we use QLSA+QTA to solve the NES since it enables us to also build
the Newton system in quantum computer which has polylogaritimic dependence on
$n$. This approach enable us to avoid matrix-matrix product, which is
inevitable in classical IPM using CG. Thus, the cost per iteration of IPM can
leverage the quadratic dependence on dimension but the linear dependence on
condition number and inverse precision can lead to exponential complexity. The
following sections will discuss how we can deal with errors and condition
numbers when we use QLSA in IPMs.

Table 1: Complexity of solving an LSP. Algorithm | Complexity | Comment
---|---|---
Factorization methods (e.g. LU) | ${\cal O}\big{(}p^{3}\big{)}$ |
Conjugate Gradient | ${\cal O}\big{(}pd\sqrt{\kappa}\log(\frac{1}{\epsilon})\big{)}$ | For LSPs with PSD matrices.
HHL Harrow2009_Quantum (23) \+ QTA van2022quantum (3) | ${\cal O}\Big{(}pd+p\frac{d^{2}\kappa^{2}\|\sigma\|^{2}}{\|M\|^{2}\epsilon^{2}}\text{polylog}(\frac{p\kappa}{\epsilon})\Big{)}$ | Using sparse encoding.
QLSA Childs2017_Quantum (15) \+ QTA van2022quantum (3) | ${\cal O}\Big{(}pd+p\frac{d\kappa\|\sigma\|}{\|M\|\epsilon}\text{polylog}(\frac{p\kappa}{\epsilon})\Big{)}$ | Using sparse encoding.
QLSA Chakraborty2018_power (14) \+ QTA van2022quantum (3) | ${\cal O}\Big{(}p^{2}+p\frac{\kappa\|\sigma\|}{\epsilon}\text{polylog}(\frac{p}{\epsilon})\Big{)}$ | Using block encoding.

##  3 The Linear Optimization Problem

Here we consider the standard form of Linear Optimization (LO) problem as
follows:

$\text{(P)}\quad\begin{aligned} \min\ c^{T}&x\\\ \text{s.t. }Ax&=b,\\\ x&\geq
0,\end{aligned}\qquad\qquad\text{(D)}\quad\begin{aligned} \max\ b^{T}y\ \ &\\\
\text{s.t. }A^{T}y+&s=c,\\\ &s\geq 0,\end{aligned}$ (1)

where $A\mathrel{\mathop{\mathchar 58\relax}}m\times n$ matrix with
rank$(A)=m$, vectors $y,b\in\mathbb{R}^{m}$, and $x,s,c\in\mathbb{R}^{n}$.
Problem (P) is called the primal problem and (D) is called the dual problem.
Due to the Strong Duality Theorem Roos2005_Interior (36), all optimal
solutions, if exist, belong to the set $\mathcal{PD}^{*}$, which is defined as

$\mathcal{PD}^{*}=\left\\{(x,y,s)\in\mathbb{R}^{n+m+n}\mathrel{\mathop{\mathchar
58\relax}}\ Ax=b,\ A^{T}y+s=c,\ x^{T}s=0,\ (x,s)\geq 0\right\\}.$

Now, we can define the optimal partition of the LO problem as

$\displaystyle B$
$\displaystyle=\\{j\in\\{1,\dots,n\\}\mathrel{\mathop{\mathchar
58\relax}}x^{*}_{j}>0\text{ for some
}(x^{*},y^{*},s^{*})\in\mathcal{PD}^{*}\\},$ $\displaystyle N$
$\displaystyle=\\{j\in\\{1,\dots,n\\}\mathrel{\mathop{\mathchar
58\relax}}s^{*}_{j}>0\text{ for some
}(x^{*},y^{*},s^{*})\in\mathcal{PD}^{*}\\}.$

It is shown that $B\cup N=\\{1,\dots,n\\}$, and $B\cap N=\varnothing$
Roos2005_Interior (36).

######  Assumption 1.

From now on, without loss of generality Roos2005_Interior (36), we assume that
the Interior Point Condition (IPC) holds, i.e., there exists a solution
$(x,y,s)$ such that

$Ax=b,\ x>0,\ A^{T}y+s=c\text{, and }s>0.$

The IPC warranties Roos2005_Interior (36) that the optimal set
$\mathcal{PD}^{*}$ is bounded, so there exists $\omega\geq 1$ such that

$\omega\geq\max\\{\|(x^{*},s^{*})\|_{\infty}\mathrel{\mathop{\mathchar
58\relax}}(x^{*},s^{*})\in\mathcal{PD}^{*}\\}.$

The central path is the curve defined by

$\mathcal{CP}=\left\\{(x,y,s)\in\mathbb{R}^{n+m+n}\mathrel{\mathop{\mathchar
58\relax}}\ Ax=b,\ A^{T}y+s=c,\ x_{i}s_{i}=\mu\text{ for
}i\in\\{1,\dots,n\\},\ x,s,\mu>0\right\\}.$

By the IPC, the central path is well-defined, and an interior feasible
solution $(x,y,s)$ exists for all $\mu>0$. Now, for any $0<\gamma_{1}<1$ and
$1\leq\gamma_{2}$, we define an infeasible neighborhood of the central path
for $\mu=\frac{x^{T}s}{n}$ as the following definition in Wright1997_Primal
(43).

$\displaystyle\mathcal{N}(\gamma_{1},\gamma_{2})=\\{(x,y,s)\in\mathbb{R}^{n+m+n}\mathrel{\mathop{\mathchar
58\relax}}\ $ $\displaystyle(x,s)>0,\ x_{i}s_{i}\geq\gamma_{1}\mu\ \text{ for
}\ i\in\\{1,\dots,n\\},\|(R_{P},R_{D})\|\leq\gamma_{2}\mu\\},$

where $R_{P}=b-Ax$, and $R_{D}=c-A^{T}y-s$. Assuming that the input data is
integral, we denote the binary length of the input data by

$L=mn+m+n+\sum_{i,j}\lceil\log(|a_{ij}|+1)\rceil+\sum_{i}\lceil\log(|c_{i}|+1)\rceil+\sum_{j}\lceil\log(|b_{j}|+1)\rceil.$

The following lemma is a classical result first proved by
Khachiyan1979_polynomial (29).

######  Lemma 1.

Let $(x^{*},y^{*},s^{*})\in\mathcal{PD}^{*}$ be a basic solution. If
$x_{i}^{*}>0$, then we have $x_{i}^{*}\geq 2^{-L}$. If $s_{i}^{*}>0$, then we
have $s_{i}^{*}\geq 2^{-L}$.

######  Theorem 3.1 (Chapter 3 in Wright1997_Primal (43)).

An exact optimal solution can be obtained by a strongly polynomial rounding
procedure when $(x,y,s)\in\mathcal{N}(\gamma_{1},\gamma_{2})$ and $\mu\leq
2^{-2L}.$

If the IPC holds, then the optimal set $\mathcal{PD}^{*}$ is bounded, and we
can find the upper bound for all the coordinates of all optimal solutions as
described in the following lemma.

######  Lemma 2 (Chapter 5 in ye2011interior (44)).

Assuming the IPC, then for any $(x^{*},y^{*},s^{*})\in\mathcal{PD}^{*}$,
$\max_{i}\\{x_{i}^{*}\\}\leq 2^{L},\text{ and }\max_{i}\\{s_{i}^{*}\\}\leq
2^{L}.$

For the theoretical purpose, we can use $\omega=2^{L}$ from Lemma 2, but in
practice, for concrete LO problems, we may find a smaller bound. We define the
set of $\zeta$-optimal solutions as

${\cal P}{\cal
D}(\zeta)=\\{(x,y,s)\in\mathbb{R}^{n+m+n}\mathrel{\mathop{\mathchar
58\relax}}\ (x,s)\geq 0,\
\frac{x^{T}s}{n}\leq\zeta,\|(R_{P},R_{D})\|\leq\zeta\\}$

##  4 An Inexact Infeasible Quantum IPM

To speed up IPMs, we use QLSAs to solve the Newton system at each iteration of
IPMs. As discussed in Section 2, QLSAs inherently produce inexact solutions.
Thus, one approach to use QLSA efficiently is to develop an Inexact Infeasible
QIPM (II-QIPM). In this paper, we utilize the KMM method proposed by
Kojima1993_primal (30) with the inexact Newton steps calculated by a QLSA.
Given $(x^{k},y^{k},s^{k})\in\mathcal{N}(\gamma_{1},\gamma_{2})$, let
$\mu^{k}=\frac{(x^{k})^{T}s^{k}}{n}$ and $0<\beta_{1}<1$ be the centering
parameter, then the Newton system is defined as

$\displaystyle A\Delta x^{k}$ $\displaystyle=b-Ax^{k},$ (2) $\displaystyle
A^{T}\Delta y^{k}+\Delta s^{k}$ $\displaystyle=c-A^{T}y^{k}-s^{k},$
$\displaystyle X^{k}\Delta s^{k}+S^{k}\Delta x^{k}$
$\displaystyle=\beta_{1}\mu^{k}e-X^{k}s^{k},$

where $e$ is all one vector with appropriate dimension,
$X^{k}=\text{diag}(x^{k})$, and $S^{k}=\text{diag}(s^{k})$. Instead of solving
the full Newton system, we may solve the Augmented system or the Normal
Equation System (NES). From the Newton system (2), the NES is formulated as

$M^{k}\Delta y^{k}=\sigma^{k},$ (NES)

where

$\displaystyle D^{k}$ $\displaystyle=(X^{k})^{1/2}(S^{k})^{-1/2},$
$\displaystyle M^{k}$ $\displaystyle=A(D^{k})^{2}A^{T},$
$\displaystyle\sigma^{k}$
$\displaystyle=A(D^{k})^{2}c-A(D^{k})^{2}A^{T}y^{k}-\beta_{1}\mu^{k}A(S^{k})^{-1}e+b-Ax^{k}$
$\displaystyle=b-\beta_{1}\mu^{k}A(S^{k})^{-1}e+A(D^{k})^{2}(c-A^{T}y^{k}-s^{k}).$

As we can see, the NES has a smaller size, $m$, than the full Newton system.
Further, the coefficient matrix of the NES is symmetric and positive definite,
thus Hermitian. Consequently, QLSAs can solve the NES efficiently. By its
nature, a QLSA generates an inexact solution $\widetilde{\Delta y}{}^{k}$ with
error bound $\|\Delta y^{k}-\widetilde{\Delta y}{}^{k}\|\leq\epsilon^{k}$.
This error leads to residual $r^{k}$ as

$M^{k}\widetilde{\Delta y}{}^{k}=\sigma^{k}+r^{k},$

where $r^{k}=M^{k}(\widetilde{\Delta y}{}^{k}-\Delta y^{k})$. After finding
$\widetilde{\Delta y}{}^{k}$ inexactly by solving the NES using QLSA, we
compute the inexact $\widetilde{\Delta x}{}^{k}$ and $\widetilde{\Delta
s}{}^{k}$ classically as

$\displaystyle\widetilde{\Delta s}{}^{k}$
$\displaystyle=c-A^{T}y^{k}-s^{k}-A^{T}\widetilde{\Delta y}{}^{k},$ (3)
$\displaystyle\widetilde{\Delta x}{}^{k}$
$\displaystyle=\beta_{1}\mu^{k}(S^{k})^{-1}e-x^{k}-(D^{k})^{2}\widetilde{\Delta
s}{}^{k}.$

As $\widetilde{\Delta s}{}^{k}$ and $\widetilde{\Delta x}{}^{k}$ are directly
calculated by equations (3), one can verify that $(\widetilde{\Delta
x}{}^{k},\widetilde{\Delta s}{}^{k},\widetilde{\Delta y}{}^{k})$ satisfies

$\displaystyle A\widetilde{\Delta x}{}^{k}$ $\displaystyle=b-Ax^{k}+r^{k},$
(4) $\displaystyle A^{T}\widetilde{\Delta y}{}^{k}+\widetilde{\Delta s}{}^{k}$
$\displaystyle=c-A^{T}y^{k}-s^{k},$ $\displaystyle X^{k}\widetilde{\Delta
s}{}^{k}+S^{k}\widetilde{\Delta x}{}^{k}$
$\displaystyle=\beta_{1}\mu^{k}e-X^{k}s^{k}.$

To have an II-IPM using (NES) with iteration complexity ${\cal O}(n^{2}L)$,
the residual norm must decrease at least ${\cal
O}(\lambda_{\text{min}}(A)\sqrt{n}\log n)$ time faster than $(x^{k})^{T}s^{k}$
where $\lambda_{\text{min}}(A)$ is the smallest singular value of $A$
Zhou2004_Polynomiality (46). We can have wider residual bound but with higher
iteration complexity of II-IPM bellavia2004convergence (10). In the literature
of preconditioning the NES, some papers modified the equations (3) and the
(NES) to transfer the residual from the first equation of (4) to its last
equation. By these changes, we can get much better bounds
Gondzio2009_Convergence (24, 35). Since tight residual bound leads to the high
complexity of QLSA+QTA, in this paper, we use a modification of the NES, which
leads to ${\cal O}(n^{2}L)$ iteration complexity of II-QIPM, where the
residual is decreasing with the rate of ${\cal O}(\sqrt{\mu^{k}})$.

Since $A$ has full row rank, one can choose an arbitrary basis $\hat{B}$, and
calculate $A_{\hat{B}}^{-1}$, $\hat{A}=A_{\hat{B}}^{-1}A$, and
$\hat{b}=A_{\hat{B}}^{-1}b$. This calculation needs ${\cal O}(m^{2}n)$
arithmetic operations and happens just one time before the iterations of IPM.
The cost of this preprocessing is dominated by the cost of II-QIPM, but it can
be reduced by using the structure of $A$. For example, if the problem is in
the canonical form, there is no need for this preprocessing. In the rest of
this paper, all methodology is applied to the preprocessed problem with input
data $(\hat{A},\hat{b},c)$. Now, we can modify the (NES) to

$\hat{M}^{k}z^{k}=\hat{\sigma}^{k}$ (MNES)

where

$\displaystyle\hat{M}^{k}$
$\displaystyle=(D_{\hat{B}}^{k})^{-1}A_{\hat{B}}^{-1}M^{k}((D_{\hat{B}}^{k})^{-1}A_{\hat{B}}^{-1})^{T}=(D_{\hat{B}}^{k})^{-1}\hat{A}(D^{k})^{2}((D_{\hat{B}}^{k})^{-1}\hat{A})^{T},$
$\displaystyle\hat{\sigma}^{k}$
$\displaystyle=(D_{\hat{B}}^{k})^{-1}A_{\hat{B}}^{-1}\sigma^{k}=(D_{\hat{B}}^{k})^{-1}\hat{b}-\beta_{1}\mu^{k}(D_{\hat{B}}^{k})^{-1}\hat{A}(S^{k})^{-1}e+(D_{\hat{B}}^{k})^{-1}\hat{A}(D^{k})^{2}(c-A^{T}y^{k}-s^{k}),$

where $D_{\hat{B}}^{k}=(X_{\hat{B}}^{k})^{1/2}(S_{\hat{B}}^{k})^{-1/2}$. We
use the following procedure to find the Newton direction by solving (MNES)
inexactly with QLSA+QTA.

  1. Step 1. Find $\tilde{z}^{k}$ such that $\hat{M}^{k}\tilde{z}^{k}=\hat{\sigma}^{k}+\hat{r}^{k}$ and $\|\hat{r}^{k}\|\leq\eta\frac{\sqrt{\mu^{k}}}{\sqrt{n}}$.

  2. Step 2. Calculate $\widetilde{\Delta y}{}^{k}=((D_{\hat{B}}^{k})^{-1}A_{\hat{B}}^{-1})^{T}\tilde{z}^{k}$.

  3. Step 3. Calculate $v^{k}=(v^{k}_{\hat{B}},v^{k}_{\hat{N}})=(D_{\hat{B}}^{k}\hat{r}^{k},0)$.

  4. Step 4. Calculate $\widetilde{\Delta s}{}^{k}=c-A^{T}y^{k}-s^{k}-A^{T}\widetilde{\Delta y}{}^{k}$.

  5. Step 5. Calculate $\widetilde{\Delta x}{}^{k}=\beta_{1}\mu^{k}(S^{k})^{-1}e-x^{k}-(D^{k})^{2}\widetilde{\Delta s}{}^{k}-v^{k}$.

The following Lemma shows how the inexact solution of (MNES) leads to residual
only in the last equation of the Newton system.

######  Lemma 3.

For the Newton direction $(\widetilde{\Delta x}{}^{k},\widetilde{\Delta
y}{}^{k},\widetilde{\Delta s}{}^{k})$, we have

$\displaystyle A\widetilde{\Delta x}{}^{k}$ $\displaystyle=b-Ax^{k},$ (5)
$\displaystyle A^{T}\widetilde{\Delta y}{}^{k}+\widetilde{\Delta s}{}^{k}$
$\displaystyle=c-A^{T}y^{k}-s^{k},$ $\displaystyle X^{k}\widetilde{\Delta
s}{}^{k}+S^{k}\widetilde{\Delta x}{}^{k}$
$\displaystyle=\beta_{1}\mu^{k}e-X^{k}s^{k}-S^{k}v^{k}.$

###### Proof.

For the Newton direction $(\widetilde{\Delta x}{}^{k},\widetilde{\Delta
y}{}^{k},\widetilde{\Delta s}{}^{k})$, one can verify that

$\displaystyle\hat{M}^{k}\tilde{z}^{k}$
$\displaystyle=\hat{\sigma}^{k}+\hat{r}^{k}$ $\displaystyle
M^{k}\widetilde{\Delta y}{}^{k}$
$\displaystyle=\sigma^{k}+A_{\hat{B}}D_{\hat{B}}^{k}\hat{r}^{k}.$

For the first equation of (5), we can write

$\displaystyle A\widetilde{\Delta x}{}^{k}=$ $\displaystyle
A(\beta_{1}\mu^{k}(S^{k})^{-1}e-x^{k}-(S^{k})^{-1}X^{k}\widetilde{\Delta
s}{}^{k}-v^{k})$ $\displaystyle=$ $\displaystyle
A(\beta_{1}\mu^{k}(S^{k})^{-1}e-x^{k}-(S^{k})^{-1}X^{k}(c-s^{k}-A^{T}y^{k}-A^{T}\widetilde{\Delta
y}{}^{k})-v^{k})$ $\displaystyle=$
$\displaystyle\beta_{1}\mu^{k}A(S^{k})^{-1}e-Ax^{k}-A(S^{k})^{-1}X^{k}c+A(S^{k})^{-1}X^{k}s^{k}+A(S^{k})^{-1}X^{k}A^{T}y^{k}$
$\displaystyle+A(S^{k})^{-1}X^{k}A^{T}\widetilde{\Delta y}{}^{k}-Av^{k}$
$\displaystyle=$
$\displaystyle\beta_{1}\mu^{k}A(S^{k})^{-1}e-A(S^{k})^{-1}X^{k}c+A(S^{k})^{-1}X^{k}A^{T}y^{k}+\sigma^{k}+A_{\hat{B}}D_{\hat{B}}^{k}\hat{r}^{k}-A_{\hat{B}}D_{\hat{B}}^{k}\hat{r}^{k}$
$\displaystyle=$ $\displaystyle b-Ax^{k}.$

The second and third equations of (5) are obtained by Steps 4 and 5. ∎

To have a convergent IPM, we need $\|S^{k}v^{k}\|_{\infty}\leq\eta\mu^{k}$,
where $0\leq\eta<1$ is an enforcing parameter.

######  Lemma 4.

For the Newton direction $(\widetilde{\Delta x}{}^{k},\widetilde{\Delta
y}{}^{k},\widetilde{\Delta s}{}^{k})$, if the residual
$\|\hat{r}^{k}\|\leq\eta\frac{\sqrt{\mu^{k}}}{\sqrt{n}}$, then
$\|S^{k}v^{k}\|_{\infty}\leq\eta\mu^{k}$.

###### Proof.

We have

$\displaystyle\|S^{k}v^{k}\|_{\infty}=\|S^{k}_{\hat{B}}v^{k}_{\hat{B}}\|_{\infty}=\|S^{k}_{\hat{B}}D_{\hat{B}}^{k}\hat{r}^{k}\|_{\infty}=\|(S^{k}_{\hat{B}})^{1/2}(X_{\hat{B}}^{k})^{1/2}\|_{\infty}\|\hat{r}^{k}\|_{\infty}\leq\sqrt{n\mu^{k}}\|\hat{r}^{k}\|\leq\eta\mu^{k}.$

∎

In the following, we show that by satisfying
$\|\hat{r}^{k}\|\leq\eta\frac{\sqrt{\mu^{k}}}{\sqrt{n}}$, then the iterations
of the II-QIPM remain in the ${\cal N}(\gamma_{1},\gamma_{2})$ neighborhood of
the central path. The following theorem presents the complexity of solving the
(MNES) system by utilizing the QLSA of Chakraborty2018_power (14). We can also
use other QLSAs discussed in Section 2, leading to different complexity
bounds.

######  Lemma 5.

The QLSA by Chakraborty2018_power (14) and the QTA by van2022quantum (3) can
build the (MNES) system, and produce a solution ${\widetilde{z}}^{k}$ for the
(MNES) system satisfying
$\|\hat{r}^{k}\|\leq\eta\frac{\sqrt{\mu^{k}}}{\sqrt{n}}$ with $\tilde{{\cal
O}}_{n,\frac{1}{\mu^{k}},\|\hat{\sigma}^{k}\|}(mn+m\sqrt{n}\frac{\kappa_{E}^{k}\|\hat{\sigma}^{k}\|}{\sqrt{\mu^{k}}})$
complexity, where $E^{k}=(D^{k}_{\hat{B}})^{-1}\hat{A}D^{k}$, and
$\kappa_{E}^{k}$ is the condition number of $E^{k}$.

###### Proof.

Building the (MNES) system in classical computer needs some matrix
multiplications, which costs ${\cal O}(m^{2}n)$ arithmetic operations. We can
write (MNES) as $E^{k}(E^{k})^{T}\tilde{z}{}^{k}=\hat{\sigma}^{k}$. As we can
see, calculating $E^{k}$ and $\hat{\sigma}^{k}$ needs just ${\cal O}(mn)$
arithmetic operations. Chakraborty2018_power (14) proposed an efficient way to
build and solve a linear system in the form
$E^{k}(E^{k})^{T}\tilde{z}{}^{k}=\hat{\sigma}^{k}$, with $\tilde{{\cal
O}}(\text{polylog}(\frac{n}{\epsilon_{QLSA}^{k}})\kappa_{E}^{k})$ complexity.
Also, we need to find the target precision for QLSA and QTA such that
$\|\hat{r}^{k}\|\leq\eta\frac{\sqrt{\mu^{k}}}{\sqrt{n}}$ is satisfied. Thus,
to have $\|\hat{r}^{k}\|\leq\eta\frac{\sqrt{\mu^{k}}}{\sqrt{n}}$, it is
sufficient to require
$\epsilon_{LSP}^{k}\leq\eta\frac{\sqrt{\mu^{k}}}{\sqrt{n}\|\hat{M}^{k}\|}$.
Based on the discussion of Section 2, we need to have

$\displaystyle\epsilon_{QLSP}^{k}\leq\eta\frac{\sqrt{\mu^{k}}}{\sqrt{n}\|\hat{\sigma}^{k}\|},\quad\epsilon_{QLSA}^{k}\leq\eta\frac{\sqrt{\mu^{k}}}{2\sqrt{n}\|\hat{\sigma}^{k}\|},\text{
and
}\epsilon_{QTA}^{k}\leq\eta\frac{\sqrt{\mu^{k}}}{2\sqrt{n}\|\hat{\sigma}^{k}\|}.$

With this target precision, the QLSA by Chakraborty2018_power (14) has ${\cal
O}(\text{polylog}(\frac{n^{1.5}\|\hat{\sigma}^{k}\|}{\sqrt{\mu^{k}}})\kappa_{E}^{k})$
complexity, and the QTA by van2022quantum (3) has ${\cal
O}(\frac{m\sqrt{n}\|\hat{\sigma}^{k}\|}{\sqrt{\mu^{k}}})$ complexity. Since
calculating $E^{k}$ and $\hat{\sigma}^{k}$ classicaly needs ${\cal O}(mn)$
arithmetic operations and the cost of storing them in QRAM is also ${\cal
O}(mn)$, the total cost of building and solving the (NES) system is

$\tilde{{\cal
O}}_{n,\frac{1}{\mu^{k}},\|\hat{\sigma}^{k}\|}(mn+m\sqrt{n}\frac{\kappa_{E}^{k}\|\hat{\sigma}^{k}\|}{\sqrt{\mu^{k}}}).$

The proof is complete. ∎

We present the II-QIPM as Algorithm 1 for solving LO problems. In this
algorithm, we use QLSA and QTA to solve the NES.

Algorithm 1 II-QIPM

1:Choose $\zeta>0$,
$\gamma_{1}\in(0,1)$,$\gamma_{2}>0$,$0<\eta<\beta_{1}<\beta_{2}<1$,

2:Choose $\omega\geq\max\\{1,\|x^{*},s^{*}\|_{\infty}\\}$.

3:$k\leftarrow 0$, $(x^{0},y^{0},s^{0})\leftarrow(\omega e,0e,\omega e)$, and
$\gamma_{2}\leftarrow\max\left\\{1,\frac{\|(R_{p}^{0},R_{D}^{0})\|}{\mu^{0}}\right\\}$

4:while $(x^{k},y^{k},s^{k})\notin\mathcal{PD}_{\zeta}$ do

5: $\mu^{k}\leftarrow\frac{(x^{k})^{T}s^{k}}{n}$

6: $E^{k}\leftarrow(D_{\hat{B}}^{k})^{-1}\hat{A}D^{k}$ and
$\hat{\sigma}^{k}\leftarrow(D_{\hat{B}}^{k})^{-1}\hat{b}-\beta_{1}\mu^{k}(D_{\hat{B}}^{k})^{-1}\hat{A}(S^{k})^{-1}e+(D_{\hat{B}}^{k})^{-1}\hat{A}(D^{k})^{2}(c-A^{T}y^{k}-s^{k})$

7:
$\epsilon_{QLSA}^{k}\leftarrow\eta\frac{\sqrt{\mu^{k}}}{2\sqrt{n}\|\sigma^{k}\|}$
and
$\epsilon_{QTA}^{k}\leftarrow\eta\frac{\sqrt{\mu^{k}}}{2\sqrt{n}\|\sigma^{k}\|}$

8: $(\Delta x^{k},\Delta y^{k},\Delta s^{k})\leftarrow$ solve
MNES($\beta_{1}$) by QLSA+QTA with precision $\epsilon_{QLSA}^{k}\text{ and
}\epsilon_{QTA}^{k}$

9: $\begin{aligned}
\hat{\alpha}^{k}\leftarrow\max\Big{\\{}\bar{\alpha}\in[0,1]\ |\ &\text{for all
}\alpha\in[0,\bar{\alpha}]\text{ we have }\\\
&\big{(}(x^{k},y^{k},s^{k})+\alpha(\Delta x^{k},\Delta y^{k},\Delta
s^{k})\big{)}\in\mathcal{N}(\gamma_{1},\gamma_{2})\text{ and }\\\
&(x^{k}+\alpha\Delta x^{k})^{T}(s^{k}+\alpha\Delta
s^{k})\leq\big{(}1-\alpha(1-\beta_{2})\big{)}(x^{k})^{T}s^{k}\Big{\\}}\end{aligned}$

10:
$(x^{k+1},y^{k+1},s^{k+1})\leftarrow(x^{k},y^{k},s^{k})+\hat{\alpha}^{k}(\Delta
x^{k},\Delta y^{k},\Delta s^{k})$

11: if $\|x_{k+1},s_{k+1}\|_{\infty}>\omega$ then

12: return Primal or dual is infeasible.

13: $k\leftarrow k+1$

14:return $(x^{k},y^{k},s^{k})$

It can be easily verified that e.g., $\beta_{1}=0.5$, $\beta_{2}=0.9995$,
$\eta=0.4$ and $\gamma_{1}=0.5$ yield a valid choice, i.e. satisfying the
conditions in the first line of Algorithm 1. In the following, we prove the
polynomial complexity of Algorithm 1.

###  4.1 Convergence of the II-QIPM

In this section, a convegence proof and iteration complexity bound are
provided for the proposed II-QIPM. The analysis closest to ours can be found
in Gondzio2009_Convergence (24) and Monteiro2003_Convergence (35) where they
are preconditioning the Newton systems which are then solved by CG. However,
we modified the NES and solved it by QLSA+QTA under different conditions and
assumptions. First, in Lemma 6, we study basic properties of the proposed II-
QIPM as presented in Algorithm 1. Then, Lemma 9 shows that the sequence
$\\{\hat{\alpha}^{k}\\}$ is strictly positive for all $k$. The iteration
complexity of Algorithm 1 is proved in Theorem 4.1. In the proof, we need to
analyze values of $\alpha^{k}$ such that at iteration $k$ of Algorithm 1, the
Newton step satisfies all the conditions in line 9. To ease notation, we use

$\displaystyle x^{k}(\alpha)=x^{k}+\alpha\Delta x^{k},\
y^{k}(\alpha)=y^{k}+\alpha\Delta y^{k},$ $\displaystyle\
s^{k}(\alpha)=s^{k}+\alpha\Delta s^{k},\
\mu^{k}(\alpha)=\frac{x^{k}(\alpha)^{T}s^{k}(\alpha)}{n},$ $\displaystyle
R_{P}^{k}(\alpha)=b-Ax^{k}(\alpha),\text{ and }$ $\displaystyle
R_{D}^{k}(\alpha)=c-A^{T}y^{k}(\alpha)-s^{k}(\alpha).$

The following lemma shows some properties of the proposed II-QIPM.

######  Lemma 6.

At iteration $k$ of Algorithm 1, for any $\alpha\in[0,1]$, with
$\mu^{k}=\frac{(x^{k})^{T}s^{k}}{n}$, we have

$\displaystyle R_{P}^{k}(\alpha)$ $\displaystyle=(1-\alpha)R_{P}^{k},$ (6a)
$\displaystyle R_{D}^{k}(\alpha)$ $\displaystyle=(1-\alpha)R_{D}^{k},$ (6b)
$\displaystyle(x^{k}(\alpha))^{T}s^{k}(\alpha)$
$\displaystyle\geq(1+\alpha(\beta_{1}-\eta-1))n\mu^{k}+\alpha^{2}(\Delta
x^{k})^{T}\Delta s^{k},$ (6c) $\displaystyle
x_{i}^{k}(\alpha)s_{i}^{k}(\alpha)$
$\displaystyle\geq(1-\alpha)x^{k}_{i}s^{k}_{i}+\alpha(\beta_{1}-\eta)\mu^{k}+\alpha^{2}\Delta
x^{k}_{i}\Delta s^{k}_{i}\text{ for }i\in\\{1,2,\dots,n\\}.$ (6d)

###### Proof.

To prove (6a) and (6b), for any $\alpha\in[0,1]$, by (2) we have

$\displaystyle R_{P}^{k}(\alpha)$ $\displaystyle=b-A(x^{k}+\alpha\Delta
x^{k})=b-Ax^{k}-\alpha A\Delta
x^{k}=b-Ax^{k}-\alpha(b-Ax^{k})=(1-\alpha)R_{P}^{k},$ $\displaystyle
R_{D}^{k}(\alpha)$ $\displaystyle=c-A^{T}(y^{k}+\alpha\Delta
y^{k})-s^{k}-\alpha\Delta s^{k}=c-A^{T}y^{k}-s^{k}-\alpha(A^{T}\Delta y+\Delta
s)=(1-\alpha)R_{D}^{k}.$

To prove (6c), using (5), we have

$\displaystyle(x^{k}+\alpha\Delta x^{k})^{T}(s^{k}+\alpha\Delta s^{k})$
$\displaystyle=(x^{k})^{T}s^{k}+\alpha[(x^{k})^{T}\Delta
s^{k}+(s^{k})^{T}\Delta x^{k}]+\alpha^{2}(\Delta x^{k})^{T}\Delta s^{k}$
$\displaystyle\geq(x^{k})^{T}s^{k}+\alpha[n\beta_{1}\mu^{k}-(x^{k})^{T}s^{k}-n\eta\mu^{k}]+\alpha^{2}(\Delta
x^{k})^{T}\Delta s^{k}$
$\displaystyle=[1+\alpha(\beta_{1}-\eta-1)](x^{k})^{T}s^{k}+\alpha^{2}(\Delta
x^{k})^{T}\Delta s^{k}.$

Using (5) again, we can similarly prove (6d) for all $i\in\\{1,2,\dots,n\\}$
as follows:

$\displaystyle(x^{k}_{i}+\alpha\Delta x^{k}_{i})(s^{k}+\alpha\Delta
s^{k}_{i})$ $\displaystyle=x^{k}_{i}s^{k}_{i}+\alpha(x^{k}_{i}\Delta
s^{k}_{i}+s^{k}_{i}\Delta x^{k}_{i})+\alpha^{2}\Delta x^{k}_{i}\Delta
s^{k}_{i}$ $\displaystyle\geq
x^{k}_{i}s^{k}_{i}+\alpha(\beta_{1}\mu^{k}-x^{k}_{i}s^{k}_{i}-\eta\mu^{k})+\alpha^{2}\Delta
x^{k}_{i}\Delta s^{k}_{i}$
$\displaystyle=(1-\alpha)x^{k}_{i}s^{k}_{i}+\alpha(\beta_{1}-\eta)\mu^{k}+\alpha^{2}\Delta
x^{k}_{i}\Delta s^{k}_{i}.$

Thus, the proof is complete. ∎

Let us define the following functions

$\displaystyle G^{k}_{i}(\alpha)$
$\displaystyle=x^{k}_{i}(\alpha)s^{k}_{i}(\alpha)-\gamma_{1}\mu^{k}(\alpha)\text{
for }i\in\\{1,\dots,n\\},$ $\displaystyle g^{k}(\alpha)$
$\displaystyle=x^{k}(\alpha)^{T}s^{k}(\alpha)-(1-\alpha)(x^{k})^{T}s^{k},$
$\displaystyle h^{k}(\alpha)$
$\displaystyle=\big{(}1-\alpha(1-\beta_{2})\big{)}(x^{k})^{T}s^{k}-x^{k}(\alpha)^{T}s^{k}(\alpha).$

We use the defined functions to check if a step with length $\alpha$
sufficiently reduces the complementarity gap and keeps the next iterate in the
neighborhood ${\cal N}(\gamma_{1},\gamma_{2})$. It is obvious that
$h^{k}(\alpha)\geq 0$ means the Armijo condition
$(x^{k}(\alpha))^{T}s^{k}(\alpha)\leq\big{(}1-\alpha(1-\beta_{2})\big{)}(x^{k})^{T}s^{k}$
holds, and the next lemma shows how we check that an iterate is in the
neighborhood of the central path.

######  Lemma 7.

For step length $0<\alpha\leq 1$, if $G^{k}_{i}(\alpha)\geq 0$ and
$g^{k}(\alpha)\geq 0$ then
$\big{(}x^{k}(\alpha),y^{k}(\alpha),s^{k}(\alpha)\big{)}\in{\cal
N}(\gamma_{1},\gamma_{2})$.

###### Proof.

It is easy to verify that that conditions $G^{k}_{i}(\alpha)\geq 0$ and
$g^{k}(\alpha)\geq 0$ lead to

$\displaystyle x^{k}_{i}(\alpha)s^{k}_{i}(\alpha)$
$\displaystyle\geq\gamma_{1}\mu^{k}(\alpha)\text{ for }i\in\\{1,\dots,n\\},$
$\displaystyle(x^{k}+\alpha\Delta x^{k})^{T}(s^{k}+\alpha\Delta s^{k})$
$\displaystyle\leq\big{(}1-\alpha(1-\beta_{2})\big{)}(x^{k})^{T}s^{k},$

respectively. Since $g^{k}(\alpha)\geq 0$, i.e.
$x^{k}(\alpha)^{T}s^{k}(\alpha)\geq(1-\alpha)(x^{k})^{T}s^{k}$, we have

$\frac{\big{\|}\big{(}R_{P}^{k}(\alpha),R_{D}^{k}(\alpha)\big{)}\big{\|}}{\|(R_{P}^{0},R_{D}^{0})\|}=\frac{(1-\alpha)\|(R_{P}^{k},R_{D}^{k})\|}{\|(R_{P}^{0},R_{D}^{0})\|}\leq\frac{(1-\alpha)\mu^{k}}{\mu^{0}}\leq\frac{\mu^{k}(\alpha)}{\mu^{0}}.$

Further, as $\gamma_{2}=\frac{\|(R_{P}^{0},R_{D}^{0})\|}{\mu^{0}}$, we can
conclude that $\big{(}x^{k}(\alpha),y^{k}(\alpha),s^{k}(\alpha)\big{)}\in{\cal
N}(\gamma_{1},\gamma_{2})$. ∎

In order to prove polynomial complexity of II-QIPM we need to find a positive
lower bound for the step length $\hat{\alpha}^{k}$. The following lemma is
bounding some remaining elements to get the step length bound.

######  Lemma 8.

There exist $0\leq\nu^{k}={\cal O}(n^{2}\mu^{k})$ such that $\left|\Delta
x^{k}_{i}\Delta s^{k}_{i}-\gamma_{1}\frac{(\Delta x^{k})^{T}\Delta
s^{k}}{n}\right|\leq\nu^{k}$ for $i\in\\{1,2,\dots,n\\}$ and $|(\Delta
x^{k})^{T}\Delta s^{k}|\leq\nu^{k}$.

###### Proof.

To prove this lemma, we need to do the following steps

  1. Step 1. finding bound ${\cal C}_{1}^{k}={\cal O}(n\mu^{k})$ such that $\omega\theta^{k-1}\|x^{k},s^{k}\|_{1}\leq{\cal C}_{1}^{k}$;

  2. Step 2. finding bound ${\cal C}_{2}^{k}={\cal O}(n\sqrt{\mu^{k}})$ such that $\|D^{-1}\Delta x^{k}\|\leq{\cal C}_{2}^{k}$ and $\|D\Delta s^{k}\|\leq{\cal C}_{2}^{k}$;

  3. Step 3. finding bound $0\leq\nu^{k}={\cal O}(n^{2}\mu^{k})$ such that

$|(\Delta x^{k})^{T}\Delta s^{k}|\leq\nu^{k}\text{ and }\left|\Delta
x^{k}_{i}\Delta s^{k}_{i}-\gamma_{1}\frac{(\Delta x^{k})^{T}\Delta
s^{k}}{n}\right|\leq\nu^{k}\text{ for }i\in\\{1,2,\dots,n\\}.$

Step 1. Let us define $\theta^{k}=\prod^{k}_{i=0}(1-\hat{\alpha}^{i})$. One
can verify that

$R_{P}^{k}=\theta^{k-1}R_{p}^{0},\text{ and }R_{D}^{k}=\theta^{k-1}R_{D}^{0}.$
(7)

Based on the definition of ${\cal N}(\gamma_{1},\gamma_{2})$ and the choice of
$\gamma_{2}$, we have

$\frac{\|(R_{P}^{k},R_{D}^{k})\|}{\mu^{k}}\leq\frac{\|(R_{P}^{0},R_{D}^{0})\|}{\mu^{0}},$
(8)

implying $\mu^{k}\geq\theta^{k-1}\mu^{0}$. We also define

$\displaystyle(\mkern 1.5mu\overline{\mkern-1.5mux}^{k},\mkern
1.5mu\overline{\mkern-1.5muy}^{k},\mkern 1.5mu\overline{\mkern-1.5mus}^{k})=$
$\displaystyle\theta^{k-1}(x^{0},y^{0},s^{0})+(1-\theta^{k-1})(x^{*},y^{*},s^{*})-(x^{k},y^{k},s^{k}),$

where, $(x^{*},y^{*},s^{*})\in{\cal P}{\cal D}^{*}$. One can verify that

$\displaystyle A^{T}\mkern 1.5mu\overline{\mkern-1.5muy}^{k}+\mkern
1.5mu\overline{\mkern-1.5mus}^{k}$ $\displaystyle=0,$ $\displaystyle A\mkern
1.5mu\overline{\mkern-1.5mux}^{k}$ $\displaystyle=0.$

Since $\mkern 1.5mu\overline{\mkern-1.5mus}^{k}$ is in the row space of $A$
and $\mkern 1.5mu\overline{\mkern-1.5mux}^{k}$ is in the null space of $A$, we
have $(\mkern 1.5mu\overline{\mkern-1.5mux}^{k})^{T}\mkern
1.5mu\overline{\mkern-1.5mus}^{k}=0$, or equivalently,

$\displaystyle[\theta^{k-1}x^{0}+(1-\theta^{k-1})x^{*}-x^{k}]^{T}[\theta^{k-1}s^{0}+(1-\theta^{k-1})s^{*}-s^{k}]$
$\displaystyle=0.$

Since $(x^{*},s^{*},x^{k},s^{k})\geq 0$, $x^{*}s^{*}=0$, and
$(x^{0},s^{0})=(\omega e,\omega e)$, we can write

$\displaystyle[\theta^{k-1}x^{0}+(1-\theta^{k-1})x^{*}-x^{k}]^{T}[\theta^{k-1}s^{0}+(1-\theta^{k-1})s^{*}-s^{k}]$
$\displaystyle=0,$ (9)

$(\theta^{k-1})^{2}(s^{0})^{T}x^{0}+\theta^{k-1}(1-\theta^{k-1})[(s^{0})^{T}x^{*}+(s^{*})^{T}x^{0}]+(s^{k})^{T}x^{k}$

$\displaystyle\geq\theta^{k-1}[(s^{0})^{T}x^{k}+(x^{0})^{T}s^{k}],$ (10)
$\displaystyle(\theta^{k-1})^{2}n\mu^{0}+2\theta^{k-1}(1-\theta^{k-1})n\mu^{0}+(s^{k})^{T}x^{k}$
$\displaystyle\geq\theta^{k-1}\omega[(e)^{T}x^{k}+(e)^{T}s^{k}],$ (11)
$\displaystyle\theta^{k-1}n\mu^{k}+2(1-\theta^{k-1})\mu^{k}+n\mu^{k}$
$\displaystyle\geq\theta^{k-1}\omega\|x^{k},s^{k}\|_{1}.$ (12)

Inequality (11) is obtained by using $(x^{0},s^{0})=(\omega e,\omega e)$ and
$\|(x^{*},s^{*})\|_{\infty}\leq\omega$, and the last inequality is obtained by
using (7). Thus, let ${\cal C}_{1}^{k}$ be defined as the left-hand side of
(12), we have ${\cal C}_{1}^{k}={\cal O}(n\mu^{k})$.

Step 2. In this step we define

$\displaystyle(\mkern 1.5mu\overline{\mkern-1.5mux}^{k},\mkern
1.5mu\overline{\mkern-1.5muy}^{k},\mkern 1.5mu\overline{\mkern-1.5mus}^{k})=$
$\displaystyle(\Delta x^{k},\Delta y^{k},\Delta
s^{k})+\theta^{k-1}(x^{0},y^{0},s^{0})-\theta^{k-1}(x^{*},y^{*},s^{*}),$

where, $(x^{*},y^{*},s^{*})\in{\cal P}{\cal D}^{*}$. Similar to Step 1, one
can verify that

$\displaystyle A^{T}\mkern 1.5mu\overline{\mkern-1.5muy}^{k}+\mkern
1.5mu\overline{\mkern-1.5mus}^{k}$ $\displaystyle=0,$ $\displaystyle A\mkern
1.5mu\overline{\mkern-1.5mux}^{k}$ $\displaystyle=0,$ $\displaystyle(\mkern
1.5mu\overline{\mkern-1.5mux}^{k})^{T}\mkern
1.5mu\overline{\mkern-1.5mus}^{k}$ $\displaystyle=0.$

Consequently, one can verify

$\displaystyle\|D^{-1}(\Delta x^{k}+\theta^{k-1}(x^{0}-x^{*}))+D(\Delta
s^{k}+\theta^{k-1}(s^{0}-s^{*}))\|$ $\displaystyle\qquad=\|D^{-1}(\Delta
x^{k}+\theta^{k-1}(x^{0}-x^{*}))\|+\|D(\Delta
s^{k}+\theta^{k-1}(s^{0}-s^{*}))\|$
$\displaystyle\qquad\leq\|(XS)^{-1/2}\|(\|XSe-\beta_{1}\mu^{k}e\|+\eta\mu^{k})+\theta^{k-1}\|D^{-1}(x^{0}-x^{*}))\|+\theta^{k-1}\|D(s^{0}-s^{*})\|.$

Now, we have

$\displaystyle\|D^{-1}\Delta
x^{k}\|\leq\|(XS)^{-1/2}\|(\|XSe-\beta_{1}\mu^{k}e\|+\eta\mu^{k})+2\theta^{k-1}\|D^{-1}(x^{0}-x^{*}))\|+2\theta^{k-1}\|D(s^{0}-s^{*})\|.$

According to pages 116-118 of Wright1997_Primal (43), we can derive the
following inequalities:

$\displaystyle\|(XS)^{-1/2}\|$
$\displaystyle\leq\frac{1}{\sqrt{\gamma_{1}}\sqrt{\mu^{k}}},$
$\displaystyle\|XSe-\beta_{1}\mu^{k}e\|$ $\displaystyle\leq n\mu^{k},$
$\displaystyle\theta^{k-1}\|D^{-1}(x^{0}-x^{*}))\|+\theta^{k-1}\|D(s^{0}-s^{*})\|$
$\displaystyle\leq\theta^{k-1}\|x^{k},s^{k}\|_{1}\|(XS)^{-1/2}\|\max(\|x^{0}-x^{*}\|,\|s^{0}-s^{*}\|)$
$\displaystyle\leq\frac{\omega\theta^{k-1}\|x^{k},s^{k}\|_{1}}{\sqrt{\gamma_{1}}\sqrt{\mu^{k}}}\leq\frac{C_{1}}{\sqrt{\gamma_{1}}\sqrt{\mu^{k}}}.$

Thus, $\|D^{-1}\Delta x^{k}\|\leq\frac{n\mu^{k}+\eta\mu^{k}+{\cal
C}_{1}^{k}}{\sqrt{\gamma_{1}}\sqrt{\mu^{k}}}={\cal C}_{2}^{k}$, where ${\cal
C}_{2}^{k}={\cal O}(n\sqrt{\mu^{k}})$. We can similarly show that $\|D\Delta
s^{k}\|\leq{\cal C}_{2}^{k}$.

Step 3. Based on Steps 1 and 2, we have

$\displaystyle(\Delta x^{k})^{T}\Delta s^{k}$ $\displaystyle\leq\|D^{-1}\Delta
x^{k}\|\|D\Delta s^{k}\|\leq({\cal C}_{2}^{k})^{2},$ $\displaystyle|\Delta
x^{k}_{i}\Delta s^{k}_{i}|$ $\displaystyle\leq\|D^{-1}\Delta x^{k}\|\|D\Delta
s^{k}\|\leq({\cal C}_{2}^{k})^{2},$ $\displaystyle\left|\Delta x^{k}_{i}\Delta
s^{k}_{i}-\gamma_{1}\frac{(\Delta x^{k})^{T}\Delta s^{k}}{n}\right|$
$\displaystyle\leq(1+\frac{\gamma_{1}}{n})({\cal C}_{2}^{k})^{2}\leq 2({\cal
C}_{2}^{k})^{2}.$

Thus, $\nu^{k}=2({\cal C}_{2}^{k})^{2}={\cal O}(n^{2}\mu^{k})$. ∎

In the next lemma, we present a strictly positive lower bound for
$\hat{\alpha}^{k}$. In the proof, we use parameters $\delta_{1}$,
$\delta_{2}$, and $\delta_{3}$ defined as follows:

$\displaystyle\delta_{1}=\frac{(1-\gamma_{1})(\beta_{1}-\eta)}{n}>0,\quad\delta_{2}=\beta_{1}-\eta>0,\quad\delta_{3}=\beta_{2}-\beta_{1}+\eta>0.$
(13)

######  Lemma 9.

At line 9 of Algorithm 1, at iteration $k$, we have

$\hat{\alpha}^{k}\geq\tilde{\alpha}^{k}\coloneqq\min\left\\{1,\ \
\min\\{\delta_{1},\delta_{2},\delta_{3}\\}\frac{(x^{k})^{T}s^{k}}{\nu^{k}}\right\\}>0.$

###### Proof.

It is enough to show that the conditions of Lemma 7 hold for all
$\alpha\in[0,\tilde{\alpha}^{k}]$. Based on Lemma 6, for any
$\alpha\in[0,\tilde{\alpha}^{k}]$, we have

$\displaystyle G^{k}_{i}(\alpha)$ $\displaystyle=(x^{k}_{i}+\alpha\Delta
x^{k}_{i})(s^{k}_{i}+\alpha\Delta
s^{k}_{i})-\gamma_{1}\frac{(x^{k}+\alpha\Delta x^{k})^{T}(s^{k}+\alpha\Delta
s^{k})}{n}$ (14a)
$\displaystyle=\textstyle(1-\alpha)x^{k}_{i}s^{k}_{i}+\alpha(\beta_{1}-\eta)\mu^{k}+\alpha^{2}\Delta
x^{k}_{i}\Delta
s^{k}_{i}-\gamma_{1}\frac{(1+\alpha(\beta_{1}-\eta-1))(x^{k})^{T}s^{k}+\alpha^{2}(\Delta
x^{k})^{T}\Delta s^{k}}{n}$ (14b)
$\displaystyle\geq\scalebox{0.95}{$\alpha^{2}(\Delta x^{k}_{i}\Delta
s^{k}_{i}-\frac{\gamma_{1}}{n}(\Delta x^{k})^{T}\Delta
s^{k})+(1-\alpha)(x^{k}_{i}s^{k}_{i}-\frac{\gamma_{1}}{n}(x^{k})^{T}s^{k})+\alpha(\beta_{1}-\eta)(1-\gamma_{1})\mu^{k}$}$
(14c)
$\displaystyle\geq-\alpha^{2}\nu^{k}+\alpha\delta_{1}(x^{k})^{T}s^{k}\geq\alpha\big{(}\delta_{1}(x^{k})^{T}s^{k}-\nu^{k}\tilde{\alpha}^{k}\big{)}\geq
0.$ (14d)

Equality (14b) follows from equation (6a) of Lemma 6, and inequality (14d) is
due to the definition of $\tilde{\alpha}^{k}$ and the neighborhood
$\mathcal{N}(\gamma_{1},\gamma_{2})$. Similarly, we show that
$g^{k}_{i}(\alpha)\geq 0$ as

$\displaystyle g^{k}_{i}(\alpha)$
$\displaystyle=x^{k}(\alpha)^{T}s^{k}(\alpha)-(1-\alpha)(x^{k})^{T}s^{k}$
$\displaystyle\geq\big{(}1+\alpha(\beta_{1}-\eta-1)\big{)}(x^{k})^{T}s^{k}+\alpha^{2}(\Delta
x^{k})^{T}\Delta s^{k}-(1-\alpha)(x^{k})^{T}s^{k}$
$\displaystyle\geq\alpha(\beta_{1}-\eta)(x^{k})^{T}s^{k}-\alpha^{2}\nu^{k}$
$\displaystyle\geq\alpha\delta_{2}(x^{k})^{T}s^{k}-\alpha^{2}\nu^{k}$
$\displaystyle\geq\alpha\big{(}\delta_{2}(x^{k})^{T}s^{k}-\nu^{k}\tilde{\alpha}^{k}\big{)}\geq
0.$

Again, by (6c) of Lemma 6, we have for all $\alpha\in[0,\tilde{\alpha}^{k}]$

$\displaystyle h^{k}(\alpha)$
$\displaystyle=\big{(}1-\alpha(1-\beta_{2})\big{)}(x^{k})^{T}s^{k}-(x^{k}+\alpha\Delta
x^{k})^{T}(s^{k}+\alpha\Delta s^{k})$
$\displaystyle=\big{(}1-\alpha(1-\beta_{2})\big{)}(x^{k})^{T}s^{k}-\big{(}1+\alpha(\beta_{1}-\eta-1)\big{)}(x^{k})^{T}s^{k}-\alpha^{2}(\Delta
x^{k})^{T}\Delta s^{k}$
$\displaystyle\geq\alpha(\beta_{2}-\beta_{1}+\eta)(x^{k})^{T}s^{k}-\alpha^{2}\nu^{k}$
$\displaystyle\geq\alpha\delta_{3}(x^{k})^{T}s^{k}-\alpha^{2}\nu^{k}$
$\displaystyle\geq\alpha\big{(}\delta_{3}(x^{k})^{T}s^{k}-\nu^{k}\tilde{\alpha}^{k}\big{)}\geq
0.$

We showed that for all $\alpha\in[0,\tilde{\alpha}^{k}]$, all the conditions
of Lemma 7 hold. Thus, we can conclude that
$\tilde{\alpha}^{k}\leq\hat{\alpha}^{k}$ and the proof is complete. ∎

By Lemma 9, we have a strictly positive lower bound for step length
$\tilde{\alpha}^{k}$ to remain in the neighborhood of the central path while
we decrease the optimality gap. In what follows, using the results of the
previous lemmas, we establish the iteration complexity of Algorithm 1.

######  Theorem 4.1.

If Algorithm 1 does not terminate in line 11, then it reaches a
$\zeta$-optimal solution in at most ${\cal O}(n^{2}\log\frac{1}{\zeta})$
iterations.

###### Proof.

Based on Lemma 9, we have

$\hat{\alpha}^{k}\geq\alpha^{k}\geq\tilde{\alpha}^{k}\coloneqq\min\left\\{1,\min\left\\{\delta_{1},\delta_{2},\delta_{3}\right\\}\frac{n\mu^{k}}{\nu^{k}}\right\\}\in(0,1].$

Hence, by the definition of the neighborhood, we have

$(x^{k})^{T}s^{k}\leq\big{(}1-\tilde{\alpha}^{k}(1-\beta_{2})\big{)}^{k}(x^{0})^{T}s^{0}.$

This implies that $\lim_{k\to\infty}(x^{k})^{T}s^{k}=0$. Thus, the algorithm
terminates in finite number of steps. Since the iterates are in the ${\cal
N}(\gamma_{1},\gamma_{2})$ neighborhood of the central path, a $\zeta$-optimal
solution is obtained when $\mu^{k}\leq\frac{\zeta}{\gamma_{1}}$. The algorithm
stops when

$\mu^{k}\leq\big{(}1-\tilde{\alpha}^{k}(1-\beta_{2})\big{)}^{k}\mu^{0}\leq\frac{\zeta\mu^{0}}{\|R_{p}^{0},R_{D}^{0}\|}.$

By the definition of $\tilde{\alpha}^{k}$ and Lemma 8, we have
$\frac{1}{\tilde{\alpha}^{k}}={\cal O}(n^{2})$. We can conclude that $k={\cal
O}(n^{2}\log\frac{\omega}{\zeta})$. The proof is complete. ∎

######  Remark 1.

The sequences of $\left\\{\mu^{k}\right\\}$, primal infeasibility
$\left\\{\|Ax^{k}-b\|\right\\}$, and dual infeasibility
$\left\\{\|A^{T}y^{k}+s^{k}-c\|\right\\}$, generated by II-QIPM, converge
linearly to zero.

######  Remark 2.

Based on Theorem 3.1, we can calculate an exact solution by a rounding
procedure if $\zeta\leq 2^{-4L}.$ Consequently, the iteration complexity of
II-QIPM for finding an exact solution is ${\cal O}(n^{2}L)$.

In the next section, we provide the total time complexity of II-QIPM.

###  4.2 Total time complexity of the II-QIPM

As discussed in Lemma 5, we need to solve the NES at each iteration of the II-
QIPM by subsequent application of QLSA and QTA, which requires $\tilde{{\cal
O}}_{n,\frac{1}{\mu^{k}},\|\hat{\sigma}^{k}\|}(mn+m\sqrt{n}\frac{\kappa_{E}^{k}\|\hat{\sigma}^{k}\|}{\sqrt{\mu^{k}}})$
computational cost. We can calculate the total time complexity of the II-QIPM
as the product of the complexity of the QLSA and the number of iterations of
the II-QIPM. The computational cost of the QLSA depends on $\kappa_{E}^{k}$,
$\|\hat{\sigma}^{k}\|$, $\|E^{k}\|_{F}$, and $\mu^{k}$ which change through
the algorithm. In the following theorem, we bound them properly and obtain the
detailed total time complexity of the proposed II-QIPM algorithm.

######  Theorem 4.2.

The total time complexity of the proposed II-QIPM with QLSA by
Chakraborty2018_power (14) and QTA by van2022quantum (3) is

$\tilde{{\cal
O}}_{n,\frac{1}{\zeta},\omega,\|\hat{A}\|,\|\hat{b}\|}\bigg{(}n^{2}\Big{[}mn+\frac{m\sqrt{n}\kappa_{\hat{A}}(\|\hat{A}\|+\|\hat{b}\|)}{\zeta^{3}}\Big{]}\bigg{)}.$

###### Proof.

To establish the total time complexity of II-QIPM, we need to analyze how the
matrices $M^{k}$ and $E^{k}$ evolve through the iterations. As in
Roos2005_Interior (36), considering the optimal partition $B$ and $N$, we have

$\frac{x_{i}^{k}}{s_{i}^{k}}={\cal O}(\frac{1}{\mu^{k}})\to\infty\text{ for
}i\in B\qquad\text{and}\qquad\frac{x_{i}^{k}}{s_{i}^{k}}={\cal O}(\mu^{k})\to
0\text{ for }i\in N.$ (15)

Appropriate bounds are provided in the following results.

  1. (i)

Based on Theorem 4.1, we have at the termination $\frac{1}{\mu^{k}}={\cal
O}(\frac{1}{\zeta})$ and $\mu^{k}\leq{\cal O}(\mu^{0})={\cal O}(\omega^{2})$.

  2. (ii)

Since $\|E^{k}\|\leq\|(D^{k})^{-1}\|\|\hat{A}\|\|D^{k}\|$, and
$\|D^{k}\|={\cal O}(\frac{1}{\sqrt{\mu^{k}}})={\cal
O}(\frac{1}{\sqrt{\zeta}})$ by (15). Similarly, $\|(D^{k})^{-1}\|_{F}={\cal
O}(\frac{1}{\sqrt{\zeta}})$, and we have $\|E^{k}\|={\cal
O}\big{(}\frac{\|\hat{A}\|}{\zeta}\big{)}$. Let $\kappa_{\hat{A}}$ be the
condition number of $\hat{A}$. Thus, we have

$\kappa^{k}_{E}={\cal O}\big{(}\kappa_{\hat{A}}\zeta^{-2}\big{)}.$

  3. (iii)

In the time complexity of II-QIPM, we also have
$\frac{\|\sigma^{k}\|}{\sqrt{\mu^{k}}}$ coming from the precision of QLSA and
QTA. We can easily verify that

$\displaystyle\frac{\|\hat{\sigma}^{k}\|}{\sqrt{\mu^{k}}}$
$\displaystyle\leq\frac{\|(D^{k})^{-1}\|}{\sqrt{\mu^{k}}}(\|\hat{b}\|+\|\hat{A}X^{k}(S^{k})^{-1}\|\|c-A^{T}y^{k}-s^{k}\|+\beta_{1}\|\mu^{k}\hat{A}(S^{k})^{-1}e\|)$
$\displaystyle\leq\frac{1}{\mu^{k}}(\|\hat{b}\|+\|\hat{A}X^{k}(S^{k})^{-1}\|\|R_{D}^{k}\|+\beta_{1}\|\hat{A}(S^{k})^{-1}e\|\mu^{k}).$

Based on the definition of the neighborhood of the central path,
$\|R_{D}^{k}\|\leq\gamma_{2}\mu^{k}$. Thus, we can get bounds ${\cal
O}(\frac{\|\sigma^{k}\|}{\sqrt{\mu^{k}}})={\cal
O}(\frac{\|\hat{A}\|+\|\hat{b}\|}{\zeta})$.

  4. (iv)

Based on Lemma 5, the complexity of QLSA by Chakraborty2018_power (14) and QTA
by van2022quantum (3) for building and solving the MNES is

$\tilde{{\cal
O}}_{n,\frac{1}{\zeta},\|\hat{A}\|,\|\hat{b}\|}\big{(}mn+m\sqrt{n}\frac{\kappa_{\hat{A}}(\|\hat{A}\|+\|\hat{b}\|)}{\zeta^{3}}\big{)}.$

Thus, the detailed time complexity of the proposed II-QIPM with QLSA by
Chakraborty2018_power (14) and QTA by van2022quantum (3) is

$\tilde{{\cal
O}}_{n,\frac{1}{\zeta},\omega,\|\hat{A}\|,\|\hat{b}\|}\bigg{(}n^{2}\Big{[}mn+\frac{m\sqrt{n}\kappa_{\hat{A}}(\|\hat{A}\|+\|\hat{b}\|)}{\zeta^{3}}\Big{]}\bigg{)}.$
(16)

The time complexity is achieved by multiplying the number of iterations of II-
QIPM and the total cost of each iteration, including building and solving it
by QLSA+QTA. Thus, the proof is complete. ∎

Let $\phi\mathrel{\mathop{\mathchar 58\relax}}=\|\hat{A}\|+\|\hat{b}\|$, the
complexity of II-QIPM can be simplified as

$\tilde{{\cal
O}}_{n,\frac{1}{\zeta},\omega,\phi}\Big{(}n^{4}\frac{\phi\kappa_{\hat{A}}}{\zeta^{3}}\Big{)}.$
(17)

In the complexity of II-QIPM, the $\frac{1}{\zeta}$ factors come from bounding
the condition number and from QTA. Based on Theorem 3.1,
$\frac{1}{\zeta}=2^{{\cal O}(L)}$ leading to exponential complexity. We
discuss how we can solve this and improve the complexity of the algorithm by
stopping II-QIPM early, e.g., $\zeta=10^{-2}$, and using the Iterative
Refinement scheme discussed in the Section 5 to improve the precision.

######  Remark 3.

Some QLSAs, such as Harrow2009_Quantum (23) and Childs2017_Quantum (15), take
advantage of the sparsity of the NES. Consider (MNES), the sparsity of
$\hat{M}^{k}=\hat{A}X^{k}(S^{k})^{-1}\hat{A}^{T}$ is independent of $X^{k}$
and $S^{k}$. So, let $d$ be the maximum number of nonzero elements in any row
or column of the matrix $\hat{A}\hat{A}^{T}$. Since matrix $\hat{A}$ has two
blocks $\begin{bmatrix}I&A_{\hat{B}}^{-1}A_{\hat{N}}\end{bmatrix}$, we have
$d\leq\min\\{m,n-m+1\\}$. As matrix $\hat{A}$ determines the sparsity of
matrix $M^{k}$, so we can take advantage of the sparsity structure of matrix
$\hat{A}$. In case $\hat{A}$ is mostly sparse, e.g., $n-m\ll m$, but has a few
dense columns, then this structure can be exploited. As described in
Andersen2004_use (2), the sparse part can be separated to solve a sparse
linear system by QLSA, and then the use of the Sherman-Morrison-Woodbury
sherman1950adjustment (38) formula allows calculating the solution of the
original linear system efficiently.

##  5 Iterative Refinement Method

Based on Remark 2, we need $\zeta\leq 2^{-4L}$ to have an exact solution for
the LO problem with integer data. The proposed II-QIPM has exponential
complexity in finding an exact solution. An Iterative Refinement (IR) scheme
can be employed to achieve polynomial complexity. By this scheme, a series of
LO problems are solved by the II-QIPM with low precision $\hat{\zeta}$, e.g.,
$\hat{\zeta}=10^{-2}$, and an IR method improves the precision to reach an
exact solution. In the classical IPM literature, the IR method by
Gleixner2016_Iterative (22) and the Rational Reconstruction method by
Gleixner2020_Linear (21) employed low-precision methods to generate a high-
precision solution. Here, we adopt the IR method to generate a high precision
solution that allows the identification of an exact optimal solution. While
using only low-precision II-QIPM, Theorem 5.1 is the foundation of the IR
method.

######  Theorem 5.1 (Gleixner2016_Iterative (22)).

Let the primal problem $(P)$ be given as (1). For $\tilde{x}\in\mathbb{R}^{n}$
and $\tilde{y}\in\mathbb{R}^{m}$ and scaling factor $\nabla>1$ consider the
refining problem $(\bar{P})$

$\min\left\\{\nabla\bar{c}^{T}x|Ax=\nabla\bar{b}\text{ and
}x\geq-\nabla\tilde{x}\right\\},$

where $\bar{c}=c-A^{T}\tilde{y}$ and $\bar{b}=b-A\tilde{x}$. Then $\bar{x}$
and $\bar{y}$ are $\hat{\zeta}$-optimal solution for problem $(\bar{P})$ if
and only if $\tilde{x}+\frac{1}{\nabla}\bar{x}$ and
$\tilde{y}+\frac{1}{\nabla}\bar{y}$ are the
$\frac{\hat{\zeta}}{\nabla}$-optimal solution for problem $(P)$.

As presented in Algorithm 2, we use the II-QIPM of Algorithm 1 to solve the
refining model and update the solution with an intelligent scaling procedure.
Theorem 5.2 shows that a polynomial number of iterations are sufficient to
reach an exact optimal solution.

Algorithm 2 IR-II-QIPM

1:$\big{(}A\in\mathbb{Z}^{m\times
n},b\in\mathbb{Z}^{m},c\in\mathbb{Z}^{n}\big{)}$

2:Choose scaling multiplier $\rho\in\mathbb{N}$ such that $\rho>1$

3:$\hat{\zeta}\leftarrow 10^{-2},\zeta\leftarrow 2^{-4L}$

4:$k\leftarrow 0$

5:$(x^{*0},y^{*0},s^{*0})\leftarrow$ Algorithm 1 $(A,b,c)$ with $\hat{\zeta}$
precision

6:while $r^{k}>\zeta$ do

7: $\pi^{k}\leftarrow\max\left\\{r^{k},\frac{1}{\rho\nabla^{k}}\right\\}$

8: $\nabla^{k}\leftarrow
2^{\left\lceil\log(\nicefrac{{1}}{{\pi_{k}}})\right\rceil}$

9: $\bar{b}^{k}\leftarrow b-Ax^{*k}$ and $\bar{c}^{k}\leftarrow c-A^{T}y^{*k}$

10: $(\hat{x}^{k},\hat{y}^{k},\hat{s}^{k})\leftarrow$ Algorithm 1
$(A,\nabla^{k}\bar{b}^{k},\nabla^{k}\bar{c}^{k})$ with $\hat{\zeta}$ precision

11: $x^{*k+1}\leftarrow x^{*k}+\frac{1}{\nabla^{k}}\hat{x}^{k}$ and
$y^{*k+1}\leftarrow y^{*k}+\frac{1}{\nabla^{k}}\hat{y}^{k}$

12: $r^{k+1}\leftarrow\max\left\\{\max_{j}|\bar{b}_{\hskip
1.0ptj}^{k+1}|,\max_{i}(-\bar{c}_{\hskip
1.0pti}^{k+1}),\sum_{i}|\bar{c}_{\hskip 1.0pti}^{k+1}x^{*k+1}_{i}|\right\\}$

13: $k\leftarrow k+1$

######  Theorem 5.2.

Let $(A,b,c)$ be integer. The number of iterations of Algorithm 2 to get
$2^{-4L}$-precise solution is at most ${\cal O}(L)$.

###### Proof.

Based on Corollary 3.6 in Gleixner2016_Iterative (22), we have

$\left\lceil\frac{\log(\zeta)}{\log(\hat{\zeta})}\right\rceil=\left\lceil\frac{-4L\log(2)}{-2}\right\rceil={\cal
O}(L).$

∎

One can observe that, except $\omega$, all parameters in the complexity of II-
QIPM as specified in (17), are constant in all iterations of the IR method.
Lemma 10 provides an upper bound for $\omega^{k}$ at iteration $k$ of the IR
method.

######  Lemma 10.

At the $k$th iteration of the IR method, let
$(\hat{x}^{*k},\hat{y}^{*k},\hat{s}^{*k})$ be an exact optimal solution of the
refining problem $(A,\nabla^{k}\bar{b}_{k-1},\nabla^{k}\bar{c}_{k-1})$, and
$\omega^{k}\geq\max\\{1,\|\hat{x}^{*k},\hat{s}^{*k}\|_{\infty}\\}$. Then,
$\omega^{k}={\cal O}\big{(}(2\rho)^{L}\big{)}$.

###### Proof.

From Theorem 5.2 and last line of Algorithm 2, we have

$\displaystyle x^{*}$ $\displaystyle=x^{*}_{0}+\sum_{k=1}^{{\cal
O}(L)}\frac{1}{\nabla^{k}}\hat{x}^{k},\text{ and
}s^{*}=s^{*}_{0}+\sum_{k=1}^{{\cal O}(L)}\frac{1}{\nabla^{k}}\hat{s}^{k}.$

Based on Lemma 1, we know that $\|x^{*},s^{*}\|_{\infty}\leq 2^{L}$, then we
have $\|\hat{x}^{*k},\hat{s}^{*k}\|_{\infty}\leq\nabla^{k}2^{L}.$

Based on the procedure of updating the scaling factor in Algorithm 2, we can
drive $\nabla^{k}={\cal O}(\rho^{k})$. We can conclude that $\omega^{k}={\cal
O}\big{(}(2\rho)^{L}\big{)}$. ∎

In Theorem 5.3, we have the total time complexity of the IR method using the
proposed II-QIPM to find an exact optimal solution for LO problems.

######  Theorem 5.3.

Let $\hat{\zeta}=10^{-2}$, then the total time complexity of finding an exact
optimal solution using the IR-II-QIPM Algorithm 1 for solving the LO problem
(1) is polynomial with

$\tilde{{\cal
O}}_{n,\|\hat{A}\|,\|\hat{b}\|}\bigg{(}n^{2}L\Big{[}mn+m\sqrt{n}\kappa_{\hat{A}}(\|\hat{A}\|+\|\hat{b}\|)\Big{]}\bigg{)},$

the arithmetic operations, where $\hat{A}$ and $\hat{b}$ are preprocessed $A$
and $b$.

###### Proof.

The proof follows from combining the result of Theorem 5.2, the total time
complexity of the proposed II-QIPM in (16), and Lemma 10. ∎

######  Corollary 1.

The simplified complexity of the proposed IR-II-QIPM using the (MNES) is

$\tilde{{\cal O}}_{n,\phi}(n^{4}L\phi\kappa_{\hat{A}}).$

For finding $\zeta$-optimal solution the complexity of IR-II-QIPM is

$\tilde{{\cal O}}_{n,\phi,\omega,\frac{1}{\zeta}}(n^{4}\phi\kappa_{\hat{A}}).$

Thus, the iterative refinement procedure speeds QIPM up exponentially with
respect to precision. It also addresses the condition number since the growing
condition number of the Newton system is replaced by the condition number of
constant matrix $\hat{A}$.

The following section presents the results of our numerical results.

##  6 Numerical Experiments

This section provides numerical results for the proposed II-QIPM using the
QISKIT AQUA quantum simulator. Due to the limited number of qubits available
in quantum computers and simulators, we use the NES, which has a smaller
dimension. The numerical results are run on a workstation with Dual Intel
Xeon® CPU E5-2630 @ 2.20 GHz (20 cores) and 64 GB of RAM. For the
computational experiments, we have developed a Python qipm package available
for public use111https://github.com/qcol-lu/qipm.

IBM has implemented a QLSA, which is similar to the HHL method, without block-
encoding and QRAM. With the current technology, the number of available qubits
in gate-based quantum computers is limited to about one hundred. One of the
main issues with quantum computers is that they are not scalable compared to
classical computers. Currently, larger NISQ222NISQ: Noisy Intermediate Scale
Quantum devices suffer more from the lack of precision. On the other hand,
quantum simulator algorithms are computationally expensive. The maximum number
of qubits in a quantum simulator is roughly similar to that in an actual
quantum computer. The main advantage of using a quantum simulator is that we
do not need to handle the noise of NISQ devices. Despite this, we still need
to deal with a high error level due to insufficient qubits and the high cost
of QLSA+QTA in order to find high-quality solutions.

We used two post-processing procedures to improve the performance of the QLSA.

  1. (i)

We first scale the linear system solution such that it satisfies the equality
$\|Mz\|=\|\sigma\|$.

  2. (ii)

We also check the sign of the linear system solution by comparing it with its
negate.

The condition number increases the solution time. Notably, the dimension of
the linear system increases the solution time as a step function for building
the quantum circuit. The dimension of the system must be a power of two. The
simulator expands the system size to the smallest possible power of two. Table
2 presents the number of qubits in quantum circuits for achieving the same
precision. The QISKIT simulator error oscillates between zero and the norm of
an actual solution. In some cases, the IBM QISKIT simulator fails, and it
reports a zero vector as a solution. While the simulator has a parameter for
tuning precision, it violates the predefined precision. We could not find any
meaningful relationship between the precision of the simulator and the
condition number of the coefficient matrix, and the dimension of the system.
In the following sections, we discuss the implementation of the proposed II-
QIPM and IR-II-QIPM using the QISKIT simulator of QLSA and evaluate their
performance.

$\kappa$ | $2^{1}$ | $2^{2}$ | $2^{3}$ | $2^{4}$ | $2^{5}$ | $2^{6}$ | $2^{7}$ | $2^{8}$ | $2^{9}$ | $2^{10}$ | $2^{11}$
---|---|---|---|---|---|---|---|---|---|---|---
#qubits | 6 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 12 | 13 | 15
Table 2: Size of the circuit for linear systems with different condition
numbers.

###  6.1 Evaluation of the II-QIPM

As discussed in Section 2, QLSAs have better dependence on the size of the
linear system than classical algorithms. However, unreliable qubits cause a
significant error in the solution of a linear system. Here, we use the IBM
QISKIT simulator to solve linear systems arising in the proposed II-QIPM. For
a fair comparison, we only consider those experiments that reach the desired
precision. The running time of the quantum computers and quantum simulators is
not comparable. Thus, instead of running time, we use the number of iterations
as a performance measure.

We use the random instance generator of generator (34), where the norm of
primal and dual solutions is set to two. The norms of the coefficient matrix
and the RHS vector are set to one and two, respectively. The condition number
of the coefficient matrix is set to two. The desired precision of the II-QIPM
is equal to $0.1$. It is worth noting that we set $\omega$ to 10 and all the
instances are feasible.

Figure 1a illustrates II-QIPM performance on instances solved to the desired
precision. As illustrated in Figure 1a, the number of iterations fluctuates
with increasing the number of variables. The noisy behavior of the QLSA and
the relatively small-scaled instances justifies this observation. Figure 1b
shows the number of iterations increases by increasing the desired precision.
The number of iterations also increases by the norm of the RHS vector, while
$\omega$ has less impact on the number of iterations (see Figures 1c and 1d).
As discussed in Section 4.2, the number of iterations is affected by the error
of QLSAs. Scaling the RHS vector norm increases the error of QLSA and, as a
result, increases the number of iterations. Furthermore, the larger $\omega$
is, the further away the starting point is from the central optimal solution,
i.e., the limit point of the central path. It implies that increasing $\omega$
should increase the number of iterations. However, the fast convergence of
IPMs and the noisy behavior of the QLSA make it difficult to observe a direct
relationship between $\omega$ and the number of iterations.

a Effect of number of variables.

b Effect of the desired precision.

c Effect of norm of RHS vector.

d Effect of $\omega$.

Figure 1: Effect of different characteristics of a LO problem on the number of
iterations.

###  6.2 Evaluation of the IR-II-QIPM

Even with smart parameter tuning and perfect implementation, the HHL simulator
has limited precision. Here, we analyze the performance of the IR method when
combined with the II-QIPM. One important parameter that highly affects the
norm of the RHS vector and, consequently, the HHL simulator’s error is
$\omega$. As the IR progresses, the norm of a solution to the LO subproblem
increases. Thus, we need to set $\omega$ based on the desired precision. Here,
we set $\omega$ to $1000$ to avoid infeasibility reported by the II-QIPM.

Figure 2a shows the logarithmic error of IR-II-QIPM for different desired
precisions. We set the LO precision to 0.01 to evaluate the effect of the IR
on the proposed II-QIPM. Figure 2b illustrates that by using IR, we can reach
higher precision. The final precision in both the pure II-QIPM and the IR-II-
QIPM is set to $10^{-4}$. However, we can still not reach the desired
precision in half of the instances because of the QLSA’s error.

a Effect of the LO precision on the IR-II-QIPM.

b Effect of iterative refinement on the II-QIPM.

Figure 2: Obtained precision versus the desired precision on the IR-II-QIPM.

By reducing the condition number, the iterative refinement method speeds up
the solution time of Newton systems. The condition number of the Newton system
is bounded by ${\cal O}(\kappa_{\hat{A}})$. Figure 3 shows how the condition
number of the solved linear systems in IR-QIPM is bounded, while the condition
number of the solved linear systems in QIPM without IR goes to infinity.

Figure 3: The condition number of linear systems in QIPM and IR-QIPM to get
$10^{-6}$-precision solution for a primal-degenerate LO with 10 variables and
5 constraints

##  7 Conclusion

This paper analyzes in detail the benefits and challenges of novel proposed
Quantum Interior Point Methods. Specifically, we analyze the use of QLSAs
within IPMs and present a convergent II-QIPM. Previous papers overlooked that
when one uses QLSA with QTA, the solution of the Newton system is inexact, and
the Newton system’s condition number goes to infinity as IPMs approach the
optimal set. Here, we also adopt the IR method to find an exact solution in
polynomial time. After addressing issues in earlier QIPMs, we proved the
correctness and convergence of the proposed II-QIPM and analyzed its
performance, both theoretically and empirically.

Algorithm | Simplified Complexity | Comment
---|---|---
II-IPM with Cholesky | ${\cal O}\big{(}n^{5}L\big{)}$ |
IR-II-IPM with CG | ${\cal O}\big{(}n^{5}L\kappa_{\hat{A}}\big{)}$ | Anticipated complexity and needing complete analysis.
QIPM of kerenidisParkas2020_quantum (26) | $\tilde{{\cal O}}_{n}\big{(}n^{2}L\kappa^{3}2^{4L}\big{)}$ | Unattainable due to using exact IPM complexity.
QIPM of Casares2020_quantum (13) | $\tilde{{\cal O}}_{n}\big{(}n^{2}L\kappa^{2}2^{2L}\big{)}$ | Unattainable due to using exact IPM complexity.
Proposed IR-II-QIPM | $\tilde{{\cal O}}_{n,\phi}\big{(}n^{4}L\phi\kappa_{\hat{A}}\big{)}$ |
Table 3: Time complexity of finding an exact solution

Table 3 compares the complexity result of the proposed IR-II-QIPM with an
analogous classical II-IPM and two recent QIPMs. The proposed IR-II-QIPM has
polynomial complexity, while the other QIPMs cannot find an exact optimal
solution in polynomial time. The exponential complexity of those QIPMs is
caused by QLSA’s error and the increasing condition number of the Newton
system. At first glance, one can get the impression that the other two QIPMs
have better complexity with respect to dimension. Still, these time
complexities cannot be attained since they only contain the iteration
complexity of exact IPMs, while these QIPMs solve the Newton system inexactly.
They also need appropriate bounds for condition number $\kappa$ and precision
$\epsilon$ based on their setting of QIPMs. To correct the complexity of the
QIPMs, at least ${\cal O}(n^{1.5})$ must be added for inexact Newton steps and
an appropriate upper bound for QLSAs’ error.

The complexity of the proposed IR-II-QIPM has better dependence on $n$ than
both its classical counterparts and the realistic complexity of its quantum
counterparts. Still, the complexity of the proposed method depends on
constants $\kappa_{\hat{A}}$ and $\|\hat{A}\|$. One may apply scaling and
preconditioning techniques for LO problems with large $\kappa_{\hat{A}}$ to
decrease $\kappa_{\hat{A}}$. Monteiro2003_Convergence (35) used the speculated
optimal partition instead of predefined basis $(\hat{B},\hat{N})$ as a
preconditioned NES. The major problem of this approach is that the cost of
calculating the precondition $A_{\hat{B}}^{-1}$ in each iteration will destroy
the quantum speed up in QIPMs. A viable research direction is to explore how
to mitigate the effects of condition number and norm of Newton systems with a
quantum-friendly approach. In addition, the proposed iterative refinement
approach can also be used to mitigate the impact of the condition number in
classical IPMs with CG in solving the NES. The complexity analysis of Inexact
IPM augmented with an iterative refinement method that uses classical
iterative solvers for solving the Newton system is part of ongoing work.
However, the anticipated complexity of such algorithms, as reported in Table
3, will have unfavorable dependence on dimension than the proposed IR-II-QIPM.
Because at each iteration, there are some matrix-matrix products to build the
Newton system. By using QLSA+QTA, we can avoid the costs of classical matrix-
matrix products and achieve polynomial quantum speedups.

We analyzed how the state-of-the-art QLSA+QTA can solve a classical linear
system problem with quadratic dependence on dimension and linear dependence on
the condition number and precision. Although there are classical polynomial-
time algorithms to solve LSPs in general form, QLSA-QTA is more scalable with
regard to dimension than classical iterative and direct approaches. The
performance of QIPMs will be improved if faster QTAs and QLSAs are proposed,
and it is worth investigating polynomial-time quantum algorithms to solve
LSPs.

We can also investigate Inexact-Feasible IPMs that are more adaptable with
QLSAs since Feasible IPMs have better complexity than infeasible IPMs. Another
direction can be developing pure QIPMs where all calculations happen in the
quantum setting. However, current NISQ devices have some limitations that
prevent having a pure QIPM. Such a method will circumvent the use of QTA
within and enjoy the fast QLSAs.

Our computational experiments show that the proposed II-QIPM embedded in the
IR scheme with the QLSA simulator of QISKIT AQUA can solve problems with
hundreds of variables to a user-defined precision. However, a classical
computer can simulate only a limited number of qubits, limiting the number of
constraints. Although there are Feasible IPMs with better complexity than the
proposed IR-II-QIPM, this paper is a significant step toward using quantum
solvers in classical methods correctly and efficiently. In addition, the IR-
II-QIPM leverage the wide neighborhood of the central path and starts with an
infeasible interior solution compared to F-IPMs. We also demonstrated for the
first time that LO problems could practically be solved using quantum solvers.

###### Acknowledgements.

This work is supported by Defense Advanced Research Projects Agency as part of
the project W911NF2010022: The Quantum Computing Revolution and Optimization:
Challenges and Opportunities.

## References

  * (1) Andris Ambainis  “Variable time amplitude amplification and quantum algorithms for linear algebra problems”  In _STACS’12 (29th Symposium on Theoretical Aspects of Computer Science)_ 14, 2012, pp. 636–647  LIPIcs  URL: https://bit.ly/3IldDGR
  * (2) Erling D. Andersen, Cornelis Roos, Tamás Terlaky, Theodore Trafalis and Joost P. Warners  “The use of low-rank updates in interior-point methods”  In _Numerical Linear Algebra and Optimization_ Science Press, Beijing, China, 2004, pp. 3–14  URL: https://bit.ly/3WM0Ly2
  * (3) Joran Apeldoorn, Arjan Cornelissen, András Gilyén and Giacomo Nannicini  “Quantum tomography using state-preparation unitaries”  In _arXiv preprint_ , 2022  DOI: 10.48550/arXiv.2207.08800
  * (4) Joran Apeldoorn and András Gilyén  “Improvements in quantum SDP-solving with applications”  In _arXiv preprint_ , 2018  DOI: 10.48550/arXiv.1804.05058
  * (5) Joran Apeldoorn and András Gilyén  “Quantum algorithms for zero-sum games”  In _arXiv preprint_ , 2019  DOI: 10.48550/arXiv.1904.03180
  * (6) Sanjeev Arora, Elad Hazan and Satyen Kale  “The multiplicative weights update method: a meta-algorithm and applications”  In _Theory of Computing_ 8.1 Theory of Computing Exchange, 2012, pp. 121–164  DOI: 10.4086/toc.2012.v008a006
  * (7) Brandon Augustino, Giacomo Nannicini, Tamás Terlaky and Luis F Zuluaga  “Quantum interior point methods for semidefinite optimization”  In _arXiv preprint_ , 2021  DOI: 10.48550/arXiv.2112.06025
  * (8) Venansius Baryamureeba and Trond Steihaug  “On the convergence of an inexact primal-dual interior point method for linear programming”  In _The Large-Scale Scientific Computing_ Berlin, Heidelberg: Springer Berlin Heidelberg, 2006, pp. 629–637  DOI: 10.1007/11666806˙72
  * (9) Stefania Bellavia  “Inexact interior-point method”  In _Journal of Optimization Theory and Applications_ 96.1 Springer, 1998, pp. 109–121  DOI: 10.1023/A:1022663100715
  * (10) Stefania Bellavia and Sandra Pieraccini  “Convergence analysis of an inexact infeasible interior point method for semidefinite programming”  In _Computational Optimization and Applications_ 29.3 Springer, 2004, pp. 289–313  DOI: 10.1023/B:COAP.0000044183.43473.34
  * (11) Fernando GSL Brandao and Krysta M Svore  “Quantum speed-ups for solving semidefinite programs”  In _2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS)_ , 2017, pp. 415–426  IEEE  DOI: 10.1109/FOCS.2017.45
  * (12) Fernando G.S.L. Brandão, Amir Kalev, Tongyang Li, Cedric Yen-Yu Lin, Krysta M Svore and Xiaodi Wu  “Quantum SDP solvers: Large speed-ups, optimality, and applications to quantum learning”  In _arXiv preprint_ , 2017  DOI: 10.48550/arXiv.1710.02581
  * (13) Pablo AM Casares and Miguel Angel Martin-Delgado  “A quantum interior-point predictor–corrector algorithm for linear programming”  In _Journal of Physics A: Mathematical and Theoretical_ 53.44 IOP Publishing, 2020, pp. 445305  DOI: 10.1088/1751-8121/abb439/meta
  * (14) Shantanav Chakraborty, András Gilyén and Stacey Jeffery  “The power of block-encoded matrix powers: Improved regression techniques via faster Hamiltonian simulation”  In _arXiv preprint_ , 2018  DOI: 10.48550/arXiv.1804.01973
  * (15) Andrew M. Childs, Robin Kothari and Rolando D. Somma  “Quantum algorithm for systems of linear equations with exponentially improved dependence on precision”  In _SIAM Journal on Computing_ 46.6 Society for Industrial & Applied Mathematics (SIAM), 2017, pp. 1920–1950  DOI: 10.1137/16m1087072
  * (16) George B. Dantzig  “Linear Programming and Extensions”  Princeton, NJ, USA: Princeton University Press, 1963  DOI: 10.1515/9781400884179
  * (17) Danial Dervovic, Mark Herbster, Peter Mountney, Simone Severini, Naïri Usher and Leonard Wossnig  “Quantum linear systems algorithms: a primer”  In _arXiv preprint_ , 2018  DOI: 10.48550/arXiv.1802.08227
  * (18) David Deutsch  “Quantum theory, the Church-Turing principle and the universal quantum computer”  In _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 400.1818 The Royal Society, 1985, pp. 97–117  DOI: 10.1098/rspa.1985.0070
  * (19) David Deutsch and Richard Jozsa  “Rapid solution of problems by quantum computation”  In _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 439.1907 The Royal Society, 1992, pp. 553–558  DOI: 10.1098/rspa.1992.0167
  * (20) Roland W. Freund, Florian Jarre and Shinji Mizuno  “Convergence of a class of inexact interior-point algorithms for linear programs”  In _Mathematics of Operations Research_ 24.1, 1999, pp. 50–71  DOI: 10.1287/moor.24.1.50
  * (21) Ambros Gleixner and Daniel E. Steffy  “Linear programming using limited-precision oracles”  In _Mathematical Programming_ 183.1, 2020, pp. 525–554  DOI: 10.1007/s10107-019-01444-6
  * (22) Ambros M. Gleixner, Daniel E. Steffy and Kati Wolter  “Iterative refinement for linear programming”  In _INFORMS Journal on Computing_ 28.3, 2016, pp. 449–464  DOI: 10.1287/ijoc.2016.0692
  * (23) Aram W. Harrow, Avinatan Hassidim and Seth Lloyd  “Quantum algorithm for linear systems of equations”  In _Physical Review Letters_ 103.15 American Physical Society (APS), 2009  DOI: 10.1103/physrevlett.103.150502
  * (24) Ghussoun Al-Jeiroudi and Jacek Gondzio  “Convergence analysis of the inexact infeasible interior-point method for linear optimization”  In _Journal of Optimization Theory and Applications_ 141.2 Springer, 2009, pp. 231–247  DOI: 10.1007/s10957-008-9500-5
  * (25) Narendra Karmarkar  “A new polynomial-time algorithm for linear programming”  In _Proceedings of the Sixteenth Annual ACM Symposium on Theory of Computing_ , STOC ’84  New York, NY, USA: Association for Computing Machinery, 1984, pp. 302–311  DOI: 10.1145/800057.808695
  * (26) Iordanis Kerenidis and Anupam Prakash  “A quantum interior point method for LPs and SDPs”  In _ACM Transactions on Quantum Computing_ 1.1 ACM New York, NY, USA, 2020, pp. 1–32  DOI: 10.1145/3406306
  * (27) Iordanis Kerenidis and Anupam Prakash  “Quantum recommendation systems”  In _arXiv preprint_ , 2016  DOI: 10.48550/arXiv.1603.08675
  * (28) Iordanis Kerenidis, Anupam Prakash and Dániel Szilágyi  “Quantum algorithms for second-order cone programming and support vector machines”  In _Quantum_ 5 Verein zur Förderung des Open Access Publizierens in den Quantenwissenschaften, 2021, pp. 427  DOI: 10.22331/q-2021-04-08-427
  * (29) Leonid Genrikhovich Khachiyan  “A polynomial algorithm in linear programming”  In _Doklady Akademii Nauk_ 244.5, 1979, pp. 1093–1096  Reports of the Russian Academy of Sciences  URL: http://mi.mathnet.ru/eng/dan42319
  * (30) Masakazu Kojima, Nimrod Megiddo and Shinji Mizuno  “A primal—dual infeasible-interior-point algorithm for linear programming”  In _Mathematical Programming_ 61.1 Springer, 1993, pp. 263–280  DOI: 10.1007/BF01582151
  * (31) János Korzak  “Convergence analysis of inexact infeasible-interior-point algorithms for solving linear programming problems”  In _SIAM Journal on Optimization_ 11.1, 2000, pp. 133–148  DOI: 10.1137/S1052623497329993
  * (32) Guang Hao Low and Isaac L Chuang  “Hamiltonian simulation by qubitization”  In _Quantum_ 3 Verein zur Förderung des Open Access Publizierens in den Quantenwissenschaften, 2019, pp. 163  DOI: 10.22331/q-2019-07-12-163
  * (33) Shinji Mizuno and Florian Jarre  “Global and polynomial-time convergence of an infeasible-interior-point algorithm using inexact computation”  In _Mathematical Programming_ 84.1, 1999, pp. 105–122  DOI: 10.1007/s10107980020a
  * (34) Mohammadhossein Mohammadisiahroudi, Ramin Fakhimi, Brandon Augustino and Tamás Terlaky  “Generating Linear, Semidefinite, and Second-order Cone Optimization Problems for Numerical Experiments”  In _arXiv preprint arXiv:2302.00711_ , 2023
  * (35) Renato Monteiro and Jerome W. O’Neal  “Convergence analysis of a long-step primal-dual infeasible interior-point LP algorithm based on iterative linear solvers”  In _Georgia Institute of Technology_ , 2003  URL: https://bit.ly/3ss0kg1
  * (36) Cornelis Roos, Tamás Terlaky and Jean-Philippe Vial  “Interior Point Methods for Linear Optimization”  Springer Science & Business Media, 2005  DOI: 10.1007/b100325
  * (37) Yousef Saad  “Iterative Methods for Sparse Linear Systems”  SIAM, 2003  DOI: 10.1137/1.9780898718003
  * (38) Jack Sherman and Winifred J Morrison  “Adjustment of an inverse matrix corresponding to a change in one element of a given matrix”  In _The Annals of Mathematical Statistics_ 21.1 JSTOR, 1950, pp. 124–127  URL: https://bit.ly/3jJFwy6
  * (39) Daniel R. Simon  “On the power of quantum computation”  In _SIAM Journal on Computing_ 26.5 Society for Industrial & Applied Mathematics (SIAM), 1997, pp. 1474–1483  DOI: 10.1137/s0097539796298637
  * (40) Tamás Terlaky  “Interior Point Methods of Mathematical Programming”  Applied Optimization Kluwer Academic Publisher, 1996  DOI: 10.1007/978-1-4613-3449-1
  * (41) Almudena Carrera Vazquez, Ralf Hiptmair and Stefan Woerner  “Enhancing the quantum linear systems algorithm using Richardson extrapolation”  In _ACM Transactions on Quantum Computing_ 3.1 ACM New York, NY, 2022, pp. 1–37  DOI: 10.1145/3490631
  * (42) Leonard Wossnig, Zhikuan Zhao and Anupam Prakash  “Quantum linear system algorithm for dense matrices”  In _Physical Review Letters_ 120.5 American Physical Society (APS), 2018  DOI: 10.1103/physrevlett.120.050502
  * (43) Stephen J. Wright  “Primal-Dual Interior-Point Methods”  SIAM, 1997  DOI: 10.1137/1.9781611971453
  * (44) Yinyu Ye  “Interior Point Algorithms: Theory and Analysis”  John Wiley & Sons, 2011  URL: https://bit.ly/3i9ztlU
  * (45) Yinyu Ye, Michael J. Todd and Shinji Mizuno  “An $O(\sqrt{n}L)$-iteration homogeneous and self-dual linear programming algorithm”  In _Mathematics of Operations Research_ 19.1 INFORMS, 1994, pp. 53–67  DOI: 10.1287/moor.19.1.53
  * (46) Guanglu Zhou and Kim-Chuan Toh  “Polynomiality of an inexact infeasible interior point algorithm for semidefinite programming”  In _Mathematical Programming_ 99.2, 2004, pp. 261–282  DOI: 10.1007/s10107-003-0431-5