A Quantum Interior-Point Predictor-Corrector Algorithm for Linear Programming

TL;DR

This paper presents a quantum algorithm for dense linear programming that accelerates classical interior-point methods by leveraging quantum linear system algorithms, providing a quantum speed-up in the number of variables.

Contribution

It introduces a quantum interior-point predictor-corrector algorithm that outputs classical solutions and achieves a speed-up over classical methods for dense LP problems.

Findings

Quantum speed-up in the number of variables n.

Algorithm outputs classical solution and optimal value.

Complexity improves over classical interior-point methods.

Abstract

We introduce a new quantum optimization algorithm for dense Linear Programming problems, which can be seen as the quantization of the Interior Point Predictor-Corrector algorithm \cite{Predictor-Corrector} using a Quantum Linear System Algorithm \cite{DenseHHL}. The (worst case) work complexity of our method is, up to polylogarithmic factors, $O(L\sqrt{n}(n+m)\overline{||M||_F}\bar{\kappa}^2\epsilon^{-2})$ for $n$ the number of variables in the cost function, $m$ the number of constraints, $\epsilon^{-1}$ the target precision, $L$ the bit length of the input data, $\overline{||M||_F}$ an upper bound to the Frobenius norm of the linear systems of equations that appear, $||M||_F$, and $\bar{\kappa}$ an upper bound to the condition number $\kappa$ of those systems of equations. This represents a quantum speed-up in the number $n$ of variables in the cost function with respect to the comparable classical Interior Point algorithms when the initial matrix of the problem $A$ is dense: if we substitute the quantum part of the algorithm by classical algorithms such as Conjugate Gradient Descent, that would mean the whole algorithm has complexity $O(L\sqrt{n}(n+m)^2\bar{\kappa} \log(\epsilon^{-1}))$, or with exact methods, at least $O(L\sqrt{n}(n+m)^{2.373})$. Also, in contrast with any Quantum Linear System Algorithm, the algorithm described in this article outputs a classical description of the solution vector, and the value of the optimal solution.

Tables

Table 1. Table 1: Comparison of complexity of different algorithms 2 2 2 For completitude we want to mention that recently nannicini2019fast has proposed quantum subroutines to speedup the simplex method. that can be used for solving dense Linear Programming problems; the first three are purely classical whereas the following three are hybrid quantum-classical. It includes only leading-order terms. QLSA stands for a dense Quantum Linear System Algorithm block-encoding_1 . Note that the algorithms QuantumSDP3 and kerenidis2018quantum can be applied to more general problems concerning Semidefinite Programming. In the table μ ≤ ‖ M ‖ F ¯ = O ​ ( n ) 𝜇 ¯ subscript norm 𝑀 𝐹 𝑂 𝑛 \mu\leq\overline{||M||_{F}}=O(\sqrt{n}) , which reflects the fact that both kerenidis2018quantum and our algorithm share common features even if they were developed independently. In reference QuantumSDP4 , theorem 24, it is proven that for many combinatorial problems the term ( R ​ r / ϵ ) = O ​ ( n + m ) 𝑅 𝑟 italic-ϵ 𝑂 𝑛 𝑚 (Rr/\epsilon)=O(n+m) , affecting the complexity of QuantumSDP3 . With respect to L 𝐿 L , for many cases the complexity of the Predictor-Corrector method does not depend on L 𝐿 L Average_performance ; but as an upper bound O ​ ( L ) 𝑂 𝐿 O(L) complexity will be common to all Interior Point methods Potra-Wright . In contrast, the Multiplicative weight method will not depend on it. Finally, the column ‘Parallelizable?’ gives the number of quantum or classical processors that can be used in parallel to solve the problem; the time complexity is divided by the corresponding amount.

Algorithms for Linear Programming	Work complexity	Parallelizable?
Pred-Corr. Predictor-Corrector + Conjugate Gradient conjugate_gradient	$O\left(L\sqrt{n}{\left(n+m\right)}^2\overline{\kappa }\log \left({ϵ}^{-1}\right)\right)$	$O\left({\left(n+m\right)}^2\right)$
Pred-Corr. Predictor-Corrector + Cholesky decomposition numerical_recipies	$O\left(L\sqrt{n}{\left(n+m\right)}^3\right)$	$O\left({\left(n+m\right)}^2\right)$
Pred-Corr. Predictor-Corrector + Optimal exact dense_exact	$O\left(L\sqrt{n}{\left(n+m\right)}^{2.373}\right)$	$O\left({\left(n+m\right)}^{2.373}\right)$
Multiplicative weights QuantumSDP3	$O\left(\left(\sqrt{n}\left(\frac{Rr}{ϵ}\right)+\sqrt{m}\right){\left(\frac{Rr}{ϵ}\right)}^4\right)$	O(1)
Quantum Interior Point with Newton system kerenidis2018quantum	$O\left(L\sqrt{n}\left(n+m\right)\mu {\overline{\kappa }}^3{ϵ}^{-2}\log \left({ϵ}^{\prime -1}\right)\right)$	$O\left(\left(n+m\right){ϵ}^{-2}\right)$
Pred-Corr. Predictor-Corrector + Block-encoding block-encoding_1 (This algorithm)	$O\left(L\sqrt{n}\left(n+m\right)\overline{{\Vert M\Vert }_F}\overline{\kappa }{ϵ}^{-2}\right)$	$O\left(\left(n+m\right){ϵ}^{-2}\right)$

Equations

$$L:=\sum _i^n\sum _j^m\lceil {\log }_2(|a_{ij}|+1)+1\rceil ,$$

$$\text{minimizes }{\mathbf{c}}^T\mathbf{x}$$

$$\text{subject to }A\mathbf{x}\ge \mathbf{b},\mathbf{x}\ge 0.$$

$$\text{maximizes }{\mathbf{b}}^T\mathbf{y}$$

$$\text{subject to }{A}^T\mathbf{y}\le \mathbf{c}.$$

$${\mathbf{b}}^T\mathbf{y}-{\mathbf{c}}^T\mathbf{x}=0.$$

$$\min \theta$$

$$\displaystyle\begin{array}[]{r@{\,}lr@{\,}lr@{\,}lr@{\,}lcl}&&+A&\bm{x}&-\bm{b}&\tau&+\bar{\bm{b}}&\theta&=&\bm{0}\\ -A^{T}&\bm{y}&&&+\bm{c}&\tau&-\bar{\bm{c}}&\theta&\geq&\bm{0}\\ +\bm{b}^{T}&\bm{y}&-\bm{c}^{T}&\bm{x}&&&+\bar{z}&\theta&\geq&0\\ -\bar{\bm{b}}^{T}&\bm{y}&+\bar{\bm{c}}^{T}&\bm{x}&-\bar{z}&\tau&&&=&-(\bm{x}^{0})^{T}\bm{s}^{0}-1\end{array}$$

$$\begin{array}{rl}& \bar{\mathbf{b}}:=\mathbf{b}-A{\mathbf{x}}^0,\bar{\mathbf{c}}:=\mathbf{c}-A^T{\mathbf{y}}^0-{\mathbf{s}}^0,\\ & \bar{z}:={\mathbf{c}}^T{\mathbf{x}}^0+1-{\mathbf{b}}^T{\mathbf{y}}^0.\end{array}$$

$$-A^T\mathbf{y}+\mathbf{c}\tau -\bar{\mathbf{c}}\theta -\mathbf{s}=0;\mathbf{s}\ge 0$$

$${\mathbf{b}}^T\mathbf{y}-{\mathbf{c}}^T\mathbf{x}+\bar{z}\theta -k=0;k\ge 0$$

$$({\mathbf{s}}^0{)}^T\mathbf{x}+({\mathbf{x}}^0{)}^T\mathbf{s}+\tau +k-(({\mathbf{x}}^0{)}^T{\mathbf{s}}^0+1)\theta =({\mathbf{x}}^0{)}^T{\mathbf{s}}^0+1.$$

$$\mathbf{y}={\mathbf{y}}^0,\mathbf{x}={\mathbf{x}}^0>0,\mathbf{s}={\mathbf{s}}^0>0,\tau =k=\theta =1.$$

$${\mathbf{y}}^0=0_{m\times 1},{\mathbf{x}}^0=1_{n\times 1}={\mathbf{s}}^0,$$

$$\begin{array}{rl}\mathcal{C}=\{& (\mathbf{y},\mathbf{x},\tau ,\theta ,\mathbf{s},k)\in {\mathcal{F}}_h^0:\\ & \left(\begin{array}{c}X\mathbf{s}\\ \tau k\end{array}\right)=\frac{{\mathbf{x}}^T\mathbf{s}+\tau k}{n+1}{1}_{(n+1)\times 1}\},\end{array}$$

$$\begin{array}{rl}\mathcal{N}(\beta )=& \{(\mathbf{y},\mathbf{x},\tau ,\theta ,\mathbf{s},k)\in {\mathcal{F}}_h^0:\left|\left|\left(\begin{array}{c}X\mathbf{s}\\ \tau k\end{array}\right)-\mu 1_{(n+1)\times 1}\right|\right|\\ & \le \beta \mu \text{ where }\mu =\frac{{\mathbf{x}}^T\mathbf{s}+\tau k}{n+1}\}.\end{array}$$

$$\left(\begin{array}{cccccc}& +A& -\mathbf{b}& +\bar{\mathbf{b}}& & \\ -A^T& & +\mathbf{c}& -\bar{\mathbf{c}}& -1& \\ +{\mathbf{b}}^T& -{\mathbf{c}}^T& & +\bar{z}& & -1\\ -{\bar{\mathbf{b}}}^T& +{\bar{\mathbf{c}}}^T& -\bar{z}& & & \end{array}\right)\left(\begin{array}{c}{\mathbf{d}}_{\mathbf{y}}\\ {\mathbf{d}}_{\mathbf{x}}\\ d_{\tau }\\ d_{\theta }\\ {\mathbf{d}}_{\mathbf{s}}\\ d_k\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right)$$

$$\left(\begin{array}{c}{X}^t{\mathbf{d}}_{\mathbf{s}}+S^t{\mathbf{d}}_{\mathbf{x}}\\ {\tau }^t{d}_k+k^t{d}_{\tau }\end{array}\right)={\gamma }^t{\mu }^t{1}_{(n+1)\times 1}-\left(\begin{array}{c}{X}^t{\mathbf{s}}^t\\ {\tau }^t{k}^t\end{array}\right),$$

$$\blockarray{ccccccc}&mn11n1\\ \block{c(cccccc)}m0A-\bm{b}\bar{\bm{b}}00\\ n-A^{T}0\bm{c}-\bar{\bm{c}}-10\\ 1\bm{b}^{T}-\bm{c}^{T}0\bar{z}0-1\\ 1-\bar{\bm{b}}^{T}\bar{\bm{c}}^{T}-\bar{z}000\\ n0S^{t}00X^{t}0\\ 100k^{t}00\tau^{t}\\ \qquad\blockarray{c}\\ \block{(c)}\bm{d_{y}}\\ \bm{d_{x}}\\ d_{\tau}\\ d_{\theta}\\ \bm{d_{s}}\\ d_{k}\\ \quad=\quad\blockarray{c}\\ \block{(c)}\bm{0}\\ \bm{0}\\ 0\\ 0\\ \gamma^{t}\mu^{t}1_{n\times 1}-X^{t}\bm{s}^{t}\\ \gamma^{t}\mu^{t}-\tau^{t}k^{t}\\ .$$

$${\mathbf{v}}^{t+1}={\mathbf{v}}^t+\delta {\mathbf{d}}_{{\mathbf{v}}^t}$$

$${\mathbf{v}}^{t+1}={\mathbf{v}}^t+{\mathbf{d}}_{{\mathbf{v}}^t}$$

$$\left(\begin{array}{c}{\mathbf{x}}^{\ast }+{\mathbf{s}}^{\ast }\\ {\tau }^{\ast }+k^{\ast }\end{array}\right)>0.$$

$$\begin{array}{r}({\mathbf{x}}^t/{\tau }^t{)}^T({\mathbf{s}}^t/{\tau }^t)\le {ϵ}_1\text{ and }\\ ({\theta }^t/{\tau }^t)||({\bar{\mathbf{b}}}^T,{\bar{\mathbf{c}}}^T)||\le {ϵ}_2,\end{array}$$

$${\tau }^t\le {ϵ}_3.$$

$$\underset{\mathbf{y},{\mathbf{x}}_B,\tau }{\min }||{\mathbf{y}}^t-\mathbf{y}|{|}^2+||{\mathbf{x}}_B^t-{\mathbf{x}}_B|{|}^2+({\tau }^t-\tau {)}^2$$

$$B{\mathbf{x}}_B-\mathbf{b}\tau =0;-B^T\mathbf{y}+{\mathbf{c}}_B\tau =0;{\mathbf{b}}^T\mathbf{y}-{\mathbf{c}}_B^T{\mathbf{x}}_B=0;$$

$$\underset{\mathbf{y},{\mathbf{x}}_B,k}{\min }||{\mathbf{y}}^t-\mathbf{y}|{|}^2+||{\mathbf{x}}_B^t-{\mathbf{x}}_B|{|}^2+(k^t-k{)}^2$$

$$B{\mathbf{x}}_B=0;-B^T\mathbf{y}=0;{\mathbf{b}}^T\mathbf{y}-{\mathbf{c}}_B^T{\mathbf{x}}_B-k=0.$$

$${\mathbf{d}}_{\mathbf{v}}=M^{-1}\mathbf{f},|d⟩=\frac{{\mathbf{d}}_{\mathbf{v}}}{||{\mathbf{d}}_{\mathbf{v}}||}.$$

$$\tilde{O}((||M|{|}_F+T_f)\kappa \text{ poly}\log (n^{\prime }{ϵ}^{-1})).$$

$$\tilde{O}((||M|{|}_F+T_f)\kappa {ϵ}^{-1}\text{ poly}\log (n^{\prime }{ϵ}^{-1})).$$

$$U_{\mathcal{M}}:|i⟩|0⟩\to \frac{1}{||M_{i\cdot }||}\sum _j{M}_{ij}|ij⟩;$$

$$U_{\mathcal{N}}:|0⟩|j⟩\to \frac{1}{||M|{|}_F}\sum _i||M_{i\cdot }|||ij⟩;$$

$$B_{i,l}=\sum _{j_1,...,j_d=l}{M}_{ij}^2.$$

Full text

Linear Programming Problem, Quantum Algorithms, Quantum Linear Approximation, Interior Point Method, Iteration Complexity, Strong Polynomiality.

††preprint: APS/123-QED

I Introduction

Linear Programming problems are among the most fundamental optimization problems Nering-Tucker ; Padberg ; Murty . Applications abound both at personal and professional fronts: improving a project delivery, scheduling of tasks, analyzing supply chain operations, shelf space optimization, designing better strategies and logistics and scheduling problems in general. Linear Programming is also used in Machine Learning where Supervised Learning works on the basis of linear programming. A system is trained to fit a mathematical model of an objective (cost) function from the labeled input data that later can predict values from unknown test data Russell-Norvig ; MRT . More specifically, linear programming is a method to find the best outcome from a linear function, such as maximum profit or lowest cost, in a mathematical model whose requirements are represented by linear constraints of the variables. Semi-Definite Programming (SDP) is an extension of Linear Programming where the objective or cost function is formulated with a non-diagonal matrix and constraints contain more general inequalities Vandernberghe-Boyd ; Todd ; Laurent-Rendl ; deKlerk .

We are in the time of small quantum computers with reduced computational capabilities due to noisy physical qubits Preskill ; Ions ; SCQ1 ; SCQ2 . The challenge of surpassing the power of current and foreseeable classical computers is attracting a lot of attention in the academia Preskill2 ; Aaronson-Arkhipov and in technological companies. This motivates the endeavour of searching for new quantum algorithms beyond the standard ones that spurred the field of quantum computation in the mid 90s (Shor, Grover, etc.) Shor ; Grover ; Nielsen-Chuang ; GMD . Only recently, a quantum algorithm for solving SDP problems has been proposed by Brandão and Svore providing us with the first quantum advantage for these optimization problems Brandao-Svore ; QuantumSDP1 ; QuantumSDP2 ; QuantumSDP3 ; QuantumSDP4 .

I.1 Background on Linear Programming.

The development of methods to solve Linear Programming problems has a long tradition starting with the Simplex Method Murty , which is simple and widely used in practice, but has (in the worst case) exponential time complexity in the number of variables. In 1979 Khachiyan proved that the ellipsoid method ensured (weak) polynomial complexity the number of variables, $O(n^{6}L)$ Khachiyan . However, in practice the ellipsoid algorithm is complicated and not competitive. In 1984 Karamark proposed the first Interior Point algorithm Karmarkar , with complexity $O(n^{3.5}L)$. It was more practical than the ellipsoid method and gave rise to a large variety of available Interior Point methods Potra-Wright . The best advantage of these methods is that, contrary to what happens in the Simplex Method, Interior Point algorithms have a worst case runtime polynomial in the number of variables. Among them, the Predictor-Corrector Method Predictor-CorrectorI ; Predictor-Corrector is arguably one of the best procedures to achieve an extremely well-behaved solution, and requires just $O(\sqrt{n}L)$ iterations. However, the Predictor-Corrector method does not explicitly indicate what method should be used to solve the linear system of equations that appear in each iteration, so the solution of the system will depend on what method is used, as can be seen in table 2. For further background on Interior Point methods we refer the reader to the review Potra-Wright .

I.2 Our algorithm

Here we present a quantum algorithm that relies on the quantization of this method. One important feature of our quantum Interior Point algorithm is that it is a hybrid algorithm: partially classical, partially quantum. This feature has become very common and a similar situation occurs with the Brandão-Svore algorithm in SDP, the Quantum Eigen-Solver for quantum chemistry Aspuru ; Mezzacapo ; Solano ; Innsbruck ; ReviewQCh , and many others, and has the advantage of requiring shorter coherence times. The core of the quantization of the Interior Point algorithm relies on the use of the block encoding techniques block-encoding_1 , that extend the applicability of the Quantum Linear System Algorithm proposed by Harrow, Hassadim and Lloyd (HHL) HHL in the case where $A$ is dense, to solve the linear system of equations that appear in the Predictor-Corrector steps.

However, in order to apply the QLSA in the context of Linear Programming, we have to solve several caveats since the straightforward application of it is doomed to failure.

The quantum Interior Point algorithm we propose benefits from several fundamental properties inherited from the classical Predictor-Corrector algorithm, and has a better performance than other classical Interior Point algorithms. In particular Predictor-Corrector :

The Predictor-Corrector method can solve the Linear Programming problem without assuming the existence of feasible or optimal solutions. 2. 2.

If the Linear Programming problem has solution, the loop of this interior point algorithm approaches feasibility and optimality at the same time for both the primal and dual problem, and if the problem is infeasible or unbounded the algorithm detects infeasibility for either the primal or dual problem. 3. 3.

The algorithm can start from any point near the center of the positive orthant.

The notions of feasible, optimal solutions etc. are defined in Sec. II where a self-contained review of the Predictor-Corrector method is presented.

The work complexity of the algorithm proposed here is $O(L\sqrt{n}(n+m)\overline{||M||_{F}}\bar{\kappa}\epsilon^{-2})$, where $n$ is the number of variables of the cost function, $m$ is the number of constraints, $L$ is the bit length of the input data (see Eq. (1)), $\overline{||M||_{F}}$ is an upper bound to the Frobenius norm of the linear systems of equations that appear, $\bar{\kappa}$ is an upper bound to the condition numbers of the linear systems of equations that appear in the Predictor-Corrector steps, and $\epsilon^{-1}$ is the precision with which one wants to solve the linear system of equations. To avoid confusion notice that in the text we will call $\epsilon$ the error and $\epsilon^{-1}$ its associated precision, because a low error means high precision and viceversa. The time complexity of the proposed quantum Interior Point algorithm can be reduced from $O(L\sqrt{n}(n+m)\overline{||M||_{F}}\bar{\kappa}\epsilon^{-2})$ to $O(L\sqrt{n}\overline{||M||_{F}}\bar{\kappa})$ distributing the work of each iteration between $O((n+m)\epsilon^{-2})$ quantum processors.

If we substituted the QLSA by a classical Linear System Algorithm, the price to pay would be, at least, an $O(\sqrt{n+m})$ increase in the work complexity, as $||M||_{F}=O(\sqrt{n+m})$ if the spectral norm of $M$ is bounded DenseHHL . For example, if we used conjugate gradient descent, the overall algorithm complexity would be $O(L\sqrt{n}(n+m)^{2}\bar{\kappa}\log(\epsilon^{-1}))$. Also, if we wanted to use an exact Linear System Algorithm the best we could hope for is the complexity it takes to exactly invert a matrix dense_exact , $O((n+m)^{2.373})$ coppersmith1990matrix , thus implying an overall work complexity for the algorithm $O(\sqrt{n}(n+m)^{2.373}L)$, that could be parallelized in $(n+m)^{2.373}$ processors to lower the time complexity $O(\sqrt{n}L)$ up to polylogarithmic terms dense_exact . A summary of these results is presented in table 2.

It is worth mentioning that our quantization approach to Linear Programming problems is radically different from the method of Brandão and Svore and this comes with several benefits. Namely, the problem of quantising linear programming using multiplicative weight methods Arora-Kale as in Brandão-Svore is that they yield an efficiency depending on parameters $R$ and $r$ of the primal and dual problems. In fact, these parameters might depend on the sizes $n,m$ of the cost function, thereby the real time complexity of the algorithm remains hidden. For instance, for some classes of problems $Rr/\epsilon=O(n)$ according to theorem 24 of QuantumSDP4 . Moreover and generically, unless specified, these $R,r$ parameters cannot be computed beforehand, but after running the algorithm (we will have a similar situation with $\bar{\kappa})$. Thus, the real efficiency of the quantum algorithm is masqueraded by overhead factors behaving badly on $R$ and $r$. Their algorithm has nevertheless a good complexity on $n$, $O(\sqrt{n+m})$, but much worse complexity on the precision, $O(\epsilon^{-5})$ for the most recent improvement of the Brandão-Svore algorithm QuantumSDP3 .

I.3 Structure of the paper

The paper has two main sections. The first reviews the Predictor-Corrector algorithm from Predictor-Corrector . It is itself divided in subsections where we explain how to initialize and terminate the algorithm, and the main loop.

In the second section we explain the changes we carry out to be able to use the QLSA from block-encoding_1 . In particular we start with a subsection discussing the condition number and then we focus on how to prepare the initial quantum states for the QLSA and read out the results using a tomography protocol from kerenidis2017recommendation . Finally we explain the QLSA, comment on the possibility of quantizing the termination of the algorithm, and devote two subsections to the complexity of the overall algorithm and its comparison with other alternatives, and the possibility of failure.

We recommend the reader to first understand the Predictor-Corrector algorithm in its classical form, and then take a look at figures 1 and 2 in order to get an overall impression of the algorithm before trying to understand the technical details.

II The Predictor-corrector algorithm

In this section we review the Predictor-Corrector algorithm of Mizuno, Todd and Ye for solving Linear Programming problems Predictor-Corrector . As stated in the original article, we will see that it performs $O(\sqrt{n}L)$ iterations of the main loop in the worst case scenario, where $n$ is the number of variables and $L$ the size of the encoding the input data in bits:

[TABLE]

for $a_{ij}$ elements of the matrix defining the problem, $A$, that is defined in the following equations (2b) and (3b). Note that $L=O(mn)$. However, in a typical case the number of iterations will not depend on $L$, but will rather be $O(\sqrt{n}\log{n})$ Average_performance .

The linear programming problem we want to solve is called primal problem (Linear Problem, LP): Given $A\in\mathbb{R}^{m\times n}$, $\bm{c}\in\mathbb{R}^{n}$ and $\bm{b}\in\mathbb{R}^{m}$, find $\bm{x}\in\mathbb{R}^{n}$ such that:

[TABLE]

The dual problem (Dual Problem, DP) has the same solution: finding $y\in\mathbb{R}^{m}$ such that

[TABLE]

Then, for linear programming problems, the primal-dual gap is [math]:

[TABLE]

A usual strategy is to use slack variables to turn all inequality constraints into equality constraints, at the cost of additional constraints. Thus, we can substitute $A^{T}\bm{y}\leq\bm{c}$ by $A^{T}\bm{y}+\bm{s}=\bm{c}$, $\bm{s}\geq\bm{0}\in\mathbb{R}^{n}$ being the slack (dual) variable to the constraint (3b).

II.1 Initialization

According to the prescription of Predictor-Corrector , one way to solve the previous problems (2b) and (3b) is to set another problem from which the solution of (2b) and (3b) can be easily obtained. This new problem is homogeneous, in the sense that there is a single non-zero constraint, and self-dual, as its dual problem is itself. Therefore, let $\bm{x}^{0}>\bm{0}\in\mathbb{R}^{n}$, $\bm{s}^{0}>\bm{0}\in\mathbb{R}^{n}$, and $\bm{y}^{0}\in\mathbb{R}^{m}$ be arbitrary initialization variables which will be chosen later on. Then, formulate the Homogeneous Linear Problem (HLP) as

[TABLE]

such that ($\bm{x}\geq\bm{0},\tau\geq 0,\tau\in\mathbb{R}$):

[TABLE]

with

[TABLE]

The last constraint from (10) is used to impose self-duality. It is also important to remark that $\bar{\bm{b}}$, $\bar{\bm{c}}$ and $\bar{z}$ indicate the infeasibility of the initial primal and dual points, and the dual gap, respectively.

Recall also that we use slack variables to convert inequality constraints into equality constraints. Those slack variables indicate the amount by which the original constraint deviates from an equality. As we have two inequality constraints, we introduce slack variables $\bm{s}\in\mathbb{R}^{n}$ for the second constraint in (10) and $k\in\mathbb{R}$ (in Predictor-Corrector denoted $\kappa$) for the third:

[TABLE]

This implies that we can rewrite the last constraint in (10) as

[TABLE]

Once we have defined these variables, theorem 2 of Predictor-Corrector proves that any point fulfilling

[TABLE]

is a feasible point, and therefore a suitable set of initialization parameters for our algorithm. A particularly simple one can choose is

[TABLE]

where $1_{n\times 1}=[1,...,1]^{T}$, and $0_{m\times 1}=[0,...,0]^{T}$.

II.2 Main loop

In this section we explain how to set up an iterative method that allows us to get close to the optimal point, following a path along the interior of the feasible region. The original references are Predictor-CorrectorI ; Predictor-Corrector . Begin defining $X:=\text{diag}(\bm{x})$ and $S:=\text{diag}(\bm{s})$. Define also $\mathcal{F}_{h}$ the set of feasible points of (HLP) $\bm{v}=(\bm{y},\bm{x},\tau,\theta,\bm{s},k)$; and $\mathcal{F}_{h}^{0}\subset\mathcal{F}_{h}$ those such that $(\bm{x},\tau,\bm{s},k)>\bm{0}$.

Finally, define the following central path in (HLP)

[TABLE]

and its neighbourhood

[TABLE]

Then, theorem 5 of Predictor-Corrector ensures that the central path lies in the feasibility region of (HLP).

In consequence, the algorithm proceeds as follows: start from an interior feasible point $\bm{v}^{0}=(\bm{y}^{0},\bm{x}^{0},\tau^{0},\theta^{0},\bm{s}^{0},k^{0})\in\mathcal{F}_{h}^{0}$. Then, recursively, form the following system of equations for variables $\bm{d}_{\bm{v}}=(\bm{d_{y}},\bm{d_{x}},d_{\tau},d_{\theta},\bm{d_{s}},d_{k})$ and $t=0,1,...\in\mathbb{N}$:

[TABLE]

where $\gamma^{t}$ takes values 0 and 1 for even and odd steps respectively, starting in $t=0$. The linear system of equations can be written in matrix form as $M^{t}\bm{d}_{\bm{v}^{t}}=\bm{f}^{t}$, i.e.

[TABLE]

So, the main loop of the algorithm consists in performing the following steps iteratively:

Predictor step: Solve (20) with $\gamma^{t}=0$ for $\bm{d}_{\bm{v}^{t}}$ where $\bm{v}^{t}=(\bm{y}^{t},\bm{x}^{t},\tau^{t},\theta^{t},\bm{s}^{t},k^{t})\in\mathcal{N}(1/4)$. Then find the biggest step length $\delta$ such that

[TABLE]

is in $\mathcal{N}(1/2)$, and update the values accordingly. Then $t\leftarrow t+1$.

Corrector step: Solve (20) with $\gamma^{t}=1$ and set

[TABLE]

that will be back in $\mathcal{N}(1/4)$. Update $t\leftarrow t+1$.

II.3 Termination

Define a strictly self-complementary solution of (HLP) $\bm{v}^{*}=(\bm{y}^{*},\bm{x}^{*},\tau^{*},\theta^{*}=0,\bm{s}^{*},k^{*})$ as an optimal solution to (HLP) that fulfills

[TABLE]

Theorem 3 in Predictor-Corrector tells us that if we have a strictly self-complementary solution to (HLP), then a solution to (LP) and (LD) exits whenever $\tau^{*}>0$, in which case $\bm{x}^{*}/\tau^{*}$ and $(\bm{y}^{*}/\tau^{*},\bm{s}^{*}/\tau^{*})$ are the solutions respectively. On the other hand, if $\tau^{*}=0$ at least one of two things will happen: $\bm{c}^{T}\bm{x}^{*}<0$, meaning that (LD) is not feasible, or $-\bm{b}^{T}\bm{y}^{*}<0$ in which case (LP) is not feasible.

The loop from the previous section will run over $t$ until one of the following two criteria are fulfilled: For $\epsilon_{1},\epsilon_{2},\epsilon_{3}$ small numbers, either

[TABLE]

We can see that the two equations in (24a) are related to the dual gap being 0, and $\theta^{*}=0$ (needed conditions for the solution to be optimal); supposing $\tau^{*}>0$. The equation (24b) is the procedure to detect $\tau^{*}=0$. $\epsilon_{1}$ and $\epsilon_{2}$ should therefore be chosen taking into account the precision we are seeking in the optimality of the solution, and the error our calculations will have. In particular, $\epsilon_{1}$ and $\epsilon_{2}$ can be taken to be the target error of the algorithm, $\epsilon$.

To get to this point we will have to iterate up to $O(L\bar{t}\sqrt{n})$ times, with $\bar{t}=\max[\log((\bm{x}^{0})^{T}(\bm{s}^{0})/(\epsilon_{1}\epsilon_{3}^{2})),\log(||(\bar{\bm{b}}^{T},\bar{\bm{c}}^{T})||/\epsilon_{2}\epsilon_{3})]$.

If the termination is due to condition (24b), then we know that there is no solution fulfilling $||(\bm{x},^{T}\bm{s}^{T})||\leq 1/(2\epsilon_{3})-1$. Therefore one should choose $\epsilon_{3}$ small enough so that the region we are exploring is reasonable. We will then consider, following Predictor-Corrector , that either (LP) or (LD) are infeasible or unbounded.

However, if termination is due to (24a), denote by $\zeta^{t}$ the index set $\{j\in 0,...,n:\bm{x}_{j}^{t}\geq\bm{s}^{t}_{j}\}$. Let also $B$ the columns of $M^{t}$ such that their index is in $\zeta^{t}$, and the rest by $C$.

Case 1: If $\tau^{t}\geq k^{t}$ solve for $\bm{y},\bm{x}_{B},\tau$

[TABLE]

Case 2: If $\tau^{t}<k^{t}$ and we solve for $\bm{y},\bm{x}_{B},$ and $k$ from

[TABLE]

The result of either of these two calculations will be the output of our algorithm, and the estimate of the solution of the (HLP) problem. In particular, $\bm{x}$ will be the calculated $\bm{x}_{B}$ in the least square projection together with $\bm{x}_{C}$, and $\bm{y}$ will be the calculated $\bm{y}$ again in the least square projection. Calculating the solution to (LP) and (LD) is then straightforward: $\bm{x}^{*}/\tau^{*}$ and $(\bm{y}^{*}/\tau^{*},\bm{s}^{*}/\tau^{*})$ respectively.

III The quantum algorithm

The aim of this section is to explain how the Quantum Linear System Algorithm (QLSA) can help us efficiently run this algorithm, in the same spirit of, for example, FEM solving the problem of the Finite Element Method. This is due to the fact that solving (20) is the most computationally expensive part of each step for large matrices.

To solve this system we propose leveraging the use of a technique called block encoding or qubitization, introduced in qubitization . Using such technique, which we shall expand further, we have the following result (theorem 34 from block-encoding_1 ).

Theorem 1**.**

block-encoding_1 : Let $M$ be an $n^{\prime}\times n^{\prime}$ Hermitian matrix (if the matrix is not Hermitian it can be included as a submatrix of a Hermitian one) with condition number $\kappa$, Frobenius norm $||M||_{F}=\sqrt{\sum_{ij}|M_{ij}|^{2}}$ and spectral norm $||M||\leq 1$. Let $\bm{f}$ be an $n^{\prime}$-dimensional unit vector, and assume that there is an oracle $\mathcal{P}_{f}$ which produces the state $\ket{f}=\bm{f}/||\bm{f}||$ in time $T_{f}$. Let also $M$ be encoded in the quantum accessible data structure indicated in III.2. Let

[TABLE]

Then,

We can prepare state $\ket{d}$ in time complexity

[TABLE] 2. 2.

We can get an $\epsilon$-multiplicative estimate of $||M^{-1}\ket{f}||$ in time

[TABLE]

Proof omited. ∎

In our case the variable $n^{\prime}$ is the size of the matrix of (20), that is $n^{\prime}=2(m+2n+3)$, the $2$ coming from symmetrisation as in the HHL algorithm. For spectral norm bounded matrices, $||M||\leq C$ constant, $||M||_{F}=O(\sqrt{n^{\prime}})$ DenseHHL . Thus, the time complexity of running the algorithm would be $O(\sqrt{n^{\prime}})$. We are also assuming $m=O(n)$. Notice that since $n$ appears in the number of iterations but $m$ does not, it is convenient to set $m\geq n$ by exchanging the primal and dual problems if needed.

An alternative could be to use a combination of the Hamiltonian simulation for dense matrices of Dense_hamiltonian_simulation with the Quantum Linear System Algorithm from QLSAchilds , that would result in slightly different complexity factors.

Let us know study how to integrate this algorithm within the Interior-Point algorithm.

III.1 The condition number $\kappa$.

We have seen that the QLSA is linear in $\kappa$. Therefore it is important to check that $\kappa$ is as low as possible.

However, preconditioning a dense matrix is much more complicated than a sparse matrix. In fact, we are not aware of any method that allows us to do it without incurring in expensive computations in the worst case. For example, the method proposed in SPAI-QLSA is only useful for sparse matrices.

Thus, as preconditioning does not seem possible, we might attempt setting an upper bound to $\kappa$ for all steps of the iteration, taking into account that only a small part of the matrix $M^{t}$ depends on $t$. The entries that depend on $t$ are the $n+1$ last rows, $2(n+1)$ entries, see (20). However, if we try doing that we will see that even if it is possible to upper bound the maximum singular value $\sigma_{\max}(M^{t})$ knowing the entries of the last rows, we cannot see a way to lower bound $\sigma_{\min}(M^{t})$, so we cannot bound the condition number.

In conclusion, we have not been able to bound the condition number from the start, so we have to rely on a rather unknown upper bound $\bar{\kappa}$. We remark that this shortcoming of the algorithm is common to both our algorithm, and the Predictor-Corrector if we substituted QLSA by other iterative methods.

III.2 Quantum state preparation and quantum-accessible data structure

In order to prepare quantum states there are many options that include State_prep_1 ; State_prep_grover ; State_prep_2 . However we are interested here in some method that can allow us to prove some quantum advantage for the case.

So, in order to do this, we introduce the method of block-encoding_1 , which additionally we will need in order to apply the QLSA.

Theorem 2**.**

kerenidis2017recommendation ; block-encoding_1 : Let $M\in\mathbb{R}^{n^{\prime}\times n^{\prime}}$ be a matrix. If $w$ is the number of nonzero entries, there is a quantum accessible data structure of size $O(w\log^{2}(n^{\prime 2}))$, which takes time $O(\log(n^{\prime 2}))$ to store or update a single entry. Once the data structure is set up, there are quantum algorithms that can perform the following maps to precision $\epsilon^{-1}$ in time $O(\text{poly}\log(n^{\prime}/\epsilon))$:

[TABLE]

[TABLE]

where $||M_{i\cdot}||$ is the $l_{2}-$norm of row $i$ of $M$. This means in particular that given a vector $f$ in this data structure, we can prepare an $\epsilon$ approximation of it, $1/||v||_{2}\sum_{i}v_{i}\ket{i}$, in time $O(\text{poly}\log(n^{\prime}/\epsilon))$.

Proof: To construct the classical data structure, create $n^{\prime}$ trees, one for each row of $M$. Then, in leaf $j$ of tree $B_{i}$ one saves the tuple ($M_{ij}^{2}$, sgn($M_{ij}$)). Also, intermediate nodes are created (that join nearby branches) so that node $l$ of tree $B_{i}$ at depth $d$ contains the value

[TABLE]

Notice that $j_{1},...,j_{d}$ is a string of values [math] and $1$, as is $l$. The root node contains the value $||M_{i\cdot}||^{2}$.

An additional tree is created taking the root nodes of all the other trees, as the leaves of the former. One can see that the depth of the structure is polylogarithmic on $n^{\prime}$, and so a single entry of $M$ can be found or updated in time polylogarithmic on $n^{\prime}$.

Now, to apply $U_{\mathcal{M}}$, we perform the following kind of controlled rotations

[TABLE]

except for the last rotation, where the sign of the leaf is included in the coefficients. It is simple to see that $U_{\mathcal{N}}$ is the same algorithm applied with the last tree, the one that contains $||M_{i\cdot}||$ for each $i$. Finally, for a vector, we have just one tree, and the procedure is the same. ∎

One may worry about two things: the first is that setting up the database might take too long, since our matrices are dense. However, notice that in $M^{t}$ only $O(n+m)$ entries depend on $t$, so the rest can be prepared at the beginning of the algorithm with an overall time complexity of $O((n+m)^{2})$, up to polylogarithmic factors. This is the same complexity as the overall algorithm when matrices are spectrally bounded.

On the other hand twice per iteration one must update the entries in the last $n+1$ rows of $M^{t}$, and prepare the data structures for the preparation of the quantum states, which will take time $O(n+m)$, but that is fine since the work complexity on $n+m$ is the same as needed to read out the result, and so will not add any complexity to the result, and it has to be done just once for each linear system of equations.

Finally, preparing the states $\ket{f}$ themselves comes at a polylogarithmic cost in both $n+m$ and $\epsilon$, so we do not need to care about it. This quantum state preparation protocol seems particularly useful when we want to solve the same linear system of equations multiple times to read out the entire solution.

III.3 Readout of the solution of QLSA.

In the same way that we need some procedure to prepare the quantum state that feeds in the QLSA, we need some way to read out the information in $\ket{d}$, defined as in equation (27).

This is a tomography problem where the solution is a pure state with real amplitudes. Thus, one may be inclined to use Amplitude Estimation to read each of the entries. However, this has the problem that entries of large normalized vectors will be small, making it costly to determine each amplitude with the additive precision provided by Amplitude Estimation.

So, instead of that procedure, we propose using the tomography algorithm explained in section 4 of kerenidis2018quantum . The complexity of the tomography is $O(n^{\prime}\epsilon^{-2})$. Theorem 4.3 in kerenidis2018quantum proves that if $\ket{\tilde{d}}$ is the output of the tomography algorithm 1, then $||\ket{\tilde{d}}-\ket{d}||_{2}\leq\sqrt{7}\epsilon$ with probability greater or equal to $1-1/n^{\prime 0.83}$. Notice that theorem 1 allows us to recover the norm of the solution. The global sign of the solution might be recovered multiplying a single row of $A$ with the solution, and comparing the result against the corresponding entry of $\bm{f}^{t}$.

III.4 Block encoding and Quantum Linear System Algorithm (QLSA)

In this subsection let us briefly review some of the results that are needed for the Quantum Linear system algorithm. In the first place, let us review what we mean by block encoding (definition 3 in block-encoding_1 ).

Definition 1**.**

Block encoding block-encoding_1 Given a $s$-qubit operator $A$, the unitary operator $U$ is an $(a,\alpha,\epsilon)$-block encoding of $A$ when

[TABLE]

That means that one can write

[TABLE]

Block encodings are useful because they allow for a variety of linear algebra techniques including multiplication of block encodings, exponentiation to positive and negative powers, and Hamiltonian simulation; efficiently. For more information we refer the reader to block-encoding_1 ; block-encoding_2 .

Another important algorithm we are using is Variable Time Amplitude Amplification Ambainis and Variable Time Amplitude Estimation block-encoding_1 . These algorithms were developed to reduce the quadratic complexity in the condition number $\kappa$ in the original HHL algorithm HHL . In that algorithm the quadratic complexity appears because of the $\kappa$ complexity in phase estimation and the $\kappa$ cost in amplitude amplification (or estimation).

The insight of Ambainis was to realise that for the eigenvalues $\lambda$ for which the phase estimation was costly ($\lambda$ small) the cost of amplitude amplification is low, and viceversa. Thus, he proposed variable time algorithms, that allow to stop some branches of the algorithm before others. Since the more expensive branches require less Amplitude Amplification, they stop earlier. The result is a reduction of the cost to linear in $\kappa$. For more information on this complex algorithm we refer the reader to those references. These two elements are the key ideas needed to obtain the result of theorem 1.

III.5 On quantizing the termination.

If there exist a feasible and optimal solution, we have seen that the loop should terminate with either procedures (25) or (26b). However it is unclear how to carry out this minimization. What we know is that it can be efficiently calculated, or substituted by any more modern and efficient method if found.

The cost of carrying out this termination by classical procedures should be not too big. In fact, according to ye_termination the overall cost is around that of one iteration of the main loop.

However, we can also propose a quantum method to finish this. It would consist on using a small Grover subroutine Grover to find all solutions of (25b) or (26b) in a small neighbourhood of the latest calculated point. After that, without reading out the state, one could apply the algorithm described in GroverMin to calculate the one with the smallest distance to the calculated point, as in (25a) or (26a).

III.6 Complexity

In this section we will indicate the complexity of our algorithm against other algorithms that can be used to solve Linear Programming problems. In particular, in table 2 we compare against the same Predictor-Corrector algorithm but using one iterative Classical Linear System Algorithm (conjugate gradient descent conjugate_gradient ), two exact classical methods (like Gauss or Cholesky decomposition, or the optimal exact algorithm dense_exact ), and against the recent quantum algorithms proposed by Brandão and Svore Brandao-Svore , and Keredinis and Prakash kerenidis2018quantum for solving Semi Definite Programming problems, a more general class of problems than those studied here (Linear Programming).

Firstly, we must take into account that, as we are using the Predictor-Corrector algorithm Predictor-Corrector , that means by construction $O(\sqrt{n}L)$ iterations of the main loop. For dense problems (as those we are considering), we should also take into account the complexity of solving two Linear Systems of Equations. The QLSA we are using is described in block-encoding_1 , with complexity $O(\overline{||M||_{F}}\bar{\kappa}\text{poly}\log((n+m)\epsilon^{-1}))$. In contrast, the fastest comparable Classical Linear System Algorithm is the conjugate gradient method conjugate_gradient , which has time complexity $O((n+m)^{2}\bar{\kappa}\log(\epsilon^{-1}))$ for general (not symmetric, positive semidefinite) dense matrices.

But we also have to take into account other procedures. Those are: the preparation of quantum states has work complexity $O(n+m)$ if we take into account the preparation of the classical data structure, and the tomography requires $O((n+m)\epsilon^{-2})$ complexity, multiplied by the complexity of QLSA.

In general we have for our algorithm a runtime of $O(L\sqrt{n}(n+m)\overline{||M||_{F}}\bar{\kappa}\epsilon^{-2})$, where each iteration comes at a runtime cost of $O((n+m)\overline{||M||_{F}}\bar{\kappa}\epsilon^{-2})$, up to polylogarithmic terms. All of this is quantum work, since the only classical operations we perform are multiplicating the vectors needed to find $\delta$ in the Predictor step, and recalculating and updating the data-base for $M^{t}$ and $f^{t}$ in each round. In the next section, III.7, we will see that some times a small number of steps of a gradient descent will be necessary after the Corrector steps, with cost $O(n+m)\text{poly}\log(n+m)$, which we already had. Finally, an additional inexpensive classical postprocessing step after the Predictor step will be needed to ensure the same convergence guarantees as the original Predictor Corrector algorithm.

It is also remarkable to mention that, thanks to the de-quantization algorithm of Ewin Tang tang2018quantum , it is possible to solve linear systems of equations (and therefore use our Interior-Point algorithm) in work complexity $O(||M||_{F}^{6}k^{6}\kappa^{16}\epsilon^{-6})$. Therefore, this classical algorithm is only useful if the rank of the matrix $k$ is low compared to $O(n+m)$ chia2018quantum ; arrazola2019quantum . However notice that we have not made any assumption about the rank of the matrix $A$, and for Linear Programming problems we do not expect this to be the case in general.

III.7 Probability of failure

Finally, we want to analyze the probability of failure of the algorithm. The reason for this is because the classical Predictor-Corrector algorithm assumes exact arithmetic, and we have to take care of the error $\epsilon$. What is more, the time complexity of the algorithm is quadratic on the associated precision $\epsilon^{-1}$, so it is computationally expensive to reach high precision. The failure of the algorithm may happen because we get out of $\mathcal{N}(\beta)$ in one of the steps. Let us now analyze if this is in fact possible.

In the Predictor steps we can state that the failure is not possible. This is because we are moving from $\mathcal{N}(1/4)$ to $\mathcal{N}(1/2)$ where we classically calculate $\delta$ such that this step is performed correctly. Therefore there is no chance of failure here, since in the worst case we can always find $\delta$ small enough such that in (21) $\bm{v}^{t+1}\in\mathcal{N}(1/2)$ if $\bm{v}^{t}\in\mathcal{N}(1/4)$.

In the Corrector steps the problem is different since now we are moving from $\mathcal{N}(1/2)$ to $\mathcal{N}(1/4)$ and there is no parameter we can tune. Therefore, one could think that the conditions of a complexity $O(\epsilon^{\prime-2})$ on the error parameter of the quantum subroutine (that forces us to have a loose error) and the possibility of a corrector step getting out of $\mathcal{N}(1/4)$ (that forces a high precision) could be incompatible.

In this section we will see that from a naive point of view they could seem incompatible, but in reality we can avoid this problem. To do this we will introduce a very simple procedure that allows us to ensure that even if we have a relatively high error, we still end up inside $\mathcal{N}(1/4)$.

The first idea one may have to check whether we end up in $\mathcal{N}(1/4)$, is to look at lemma 3 in Predictor-CorrectorI , in which Predictor-Corrector is based. If we analyze the details of the proof we can see that in fact the exact solution of the Corrector Step is not only in $\mathcal{N}(1/4)$ but also in $\mathcal{N}(1/4\sqrt{2})$. This means that we need to lower the error sufficiently so that if the exact arithmetic result is in $\mathcal{N}(1/4\sqrt{2})$, then the approximate solution is in $\mathcal{N}(1/4)$. In the appendix A we see that in fact this would require very low error condition, something that the quantum subroutine cannot provide in reasonable time complexity. There we derive the full calculation for completeness, even if they are a bit tedious. In particular the precision needed would be $\epsilon^{\prime-1}=O(n^{\text{poly}\log n})$. The reason for the $n$ dependence on the precision is related to the fact that the error $\epsilon^{\prime}$ affects each coordinate of $\bm{x}^{t}$ or $\bm{s}^{t}$, and there are $O(n)$ of them.

To solve this problem let us introduce a different strategy. Choose a low precision on the error for the corrector step $\epsilon^{\prime-1}=O(n^{1/a})$, with $a$ sufficiently large. This is ‘cheap’ in terms of the complexity of the quantum subroutine for $a$ relatively large. We will now shift the result slightly ($\epsilon^{\prime\prime}$ close) but in such a way that the new point fulfils $||X\bm{s}-\mathbf{1}\mu||\leq\beta\mu$, thus introducing an error of size $\epsilon^{\prime\prime}$ in each coordinate of the result. To do that we can perform one step of gradient descent (or a small number of them). Therefore we shift each coordinate of the output of the Corrector step at most $\epsilon^{\prime\prime}$. How large is that $\epsilon^{\prime\prime}$? Intuitively one can think that if the error of the Corrector step is at most $\epsilon^{\prime}$ for each coordinate, then $\epsilon^{\prime\prime}=O(\epsilon^{\prime})$, which we already had anyway. As we calculate in appendix B this is in fact true and performing gradient descent with such $\epsilon^{\prime\prime}$ is sufficient to end up inside $\mathcal{N}(1/4)$. We want to emphasise that the legitimacy of this procedure rests on the fact that the precision on the Corrector step is only important because it could cause the algorithm to fail by ending outside $\mathcal{N}(1/4)$, apart from affecting the final error of the algorithm, $\epsilon$. Finally, note that this procedure has time complexity $O(n+m)$ as calculated in appendix B, so it adds no additional complexity to the algorithm. In appendix C we also see that imposing (97), which has no complexity consequences, we can prove that our modification does not affect the convergence of our algorithm, and so the number of iterations of our algorithm is the same as the original Predictor-Corrector.

The conclusion of the previous is that, in order to comply with both the requirements of having a low complexity in $\epsilon^{\prime}$ in the quantum subroutine and still ending up in $\mathcal{N}(1/4)$ in the end of the corrector step, we can $O(\epsilon^{\prime})$-shift the output of the quantum subroutine and make it comply with the second requirement while $\epsilon^{\prime}$ is sufficiently large. This only imposes is that the final precision we want from our algorithm ($\epsilon^{-1}$) cannot be greater than the one induced by the error $\epsilon^{\prime}$ that the Corrector step and ‘shifting’ procedure has in the worst case. Or, in other words, $\epsilon\geq O(\epsilon^{\prime})=O(n^{-1/a})$. This should nevertheless pose no problem since we assume that in general we do not ask for a greater precision in each variable for a larger problem, what means that $\epsilon=O(1)$ in the variable $n$. In other words, even if the target error $\epsilon$ for each variable is small, it will have a priori no relation with the number of variables $n$ or constraints $m$. It is reasonable to assume this since in general one is interested in a given precision for each variable, independent of their number.

Another way our algorithm could fail is because in each step we saw that the tomography has a failure probability of $1/(3+2n+m)^{0.83}$. Let us analyze if that is a problem. First notice that after each step we can easily check that the solution is correct up to $\epsilon$ precision. For example, using the quantum accessible data structure we can prepare the trial solution, and perform a swap test to check that the solution is $\epsilon$-close to the actual solution. The swap test may be combined with Amplitude Estimation and the median lemma from nagaj2009fast to yield exponentially small failure probability. This means that if we perform $c$ attempts at each step, the probability of failure in each iteration decreases to $1/(3+2n+m)^{0.83c}$.

So, we have to choose a small constant $c$ such that after $O(\sqrt{n}L)$ iterations the total probability of failure $(\sqrt{n}L)/(3+2n+m)^{0.83c}\ll 1$. Since in many cases the number of iterations will be $O(\sqrt{n}\log n)$ taking $c=1$ should do, but in general $L=O(mn)$, so we can take $c=4>2.5/0.83$. In conclusion, the probability of failure of the tomography algorithm should pose no problem.

Finally, in order to decrease the final error of the algorithm even further, in practice it could be a good idea to perform a small (constant number) of iterations classically (doing this will not affect the theoretical complexity). Theoretically though, this does not give any guarantee of success, since the number of additional iterations would be $O(\sqrt{n}L(\bar{t}_{\epsilon_{a}}-\bar{t}_{\epsilon_{b}}))$, when we want to lower the error from $\epsilon_{b}$ to $\epsilon_{a}$. Even if the difference is small, there is a dependence on $\sqrt{n}$ that would increase the complexity of the algorithm making it comparable to the complexity of classical algorithms.

We already explained in section II.3 that the number of iterations should be $\bar{t}=\max[\log((\bm{x}^{0})^{T}(\bm{s}^{0})/(\epsilon_{1}\epsilon_{3}^{2})),\log(||(\bar{\bm{b}}^{T},\bar{\bm{c}}^{T})||/\epsilon_{2}\epsilon_{3})]$. In particular, if any one supposes every entry of $\bm{x}^{0}$, $\bm{s}^{0}$, $\bm{b}$ and $\bm{c}$ to be of order $O(1)$, then $\bar{t}_{a}=\log(O(n+1)/\epsilon^{2}_{a})$ and similarly for $\bar{t}_{b}$ and $\epsilon_{b}$. This would mean that to lower the error from $\epsilon_{b}$ to $\epsilon_{a}$, one would need

[TABLE]

what implies that

[TABLE]

Since $\lim_{n\rightarrow\infty}e^{-1/2\sqrt{n}}=1$, we would need the final error $\epsilon_{a}$ to be very similar to the one we achieve with the quantum subroutine $\epsilon_{b}$. Or, with the notation we are using: $\epsilon_{a}=O(\epsilon_{b})=O(n^{-1/a})$.

To summarize all the components in our quantum Predictor-Corrector algorithm and the interrelations among them, we show a diagram in Fig. 1 in the form of a flow chart of actions from the initialization to the termination of the quantum algorithm providing the solution to the given (LP) and (LD) problems in (2b) and (3b).

IV Conclusions

Quantization of Linear Programming problems thus far have been achieved by using multiplicative weight methods as in the pioneering work of Brandão and Svore for Semidefinite Programming (SDP) problems Brandao-Svore , which are more general than Linear Programming problems. In this work, we have enlarged the range of applicability of quantum algorithms for Linear Programming problems by using Interior Point methods instead. Specifically, our quantum algorithm relies on a type of Interior Point algorithm known as the Predictor-Corrector method that is very well behaved with respect to the feasibility, optimality conditions of the output solution and the iteration complexity.

The core of our quantum Interior Point algorithm is the application of block encoding techniques as a Quantum Linear System Algorithm block-encoding_1 to an auxiliary system of equations that comprises an homogeneous self-dual primal-dual problem associated to the original Linear Programming problem. This is the basis of the Predictor-Corrector method, from which many of its good properties derive. In particular, the iteration complexity of the classical part scales as the square root of the size $n$ of the cost function. Then, the advantage of the quantum part of the Predictor-Corrector algorithm amounts to a faster solution of the linear system of equations, with complexity $O((n+m)\sqrt{n+m})$ including the readout process, as can be seen in table 2.

Hence, this quantum Predictor Corrector algorithm is an hybrid algorithm, partially classical, partially quantum. Applying the QLSA is not an easy task if we want to achieve a clear advantage. These algorithms come with several shortcomings, some of which have been recently overcome SPAI-QLSA for sparse linear systems. Also, even though the solution to the system of linear equations can be obtained in a quantum state, then it is not easy to extract all the information provided by the solution. One has to be satisfied by obtaining partial information from the encoded solution such as an expectation value of interest or a single entry of the vector solution. Nevertheless this does not stop us from obtaining a polynomial quantum advantage in the number of variables of the problem $n$, if the matrix is dense, well-conditioned, with $m=O(n)$, and with a constant-bounded spectral norm.

Acknowledgements

We acknowledge financial support from the Spanish MINECO grants FIS2015-67411P, and the CAM research consortium QUITEMAD-CM, Grant No.S2018/TCS-4342. The research of M.A.M.-D. has been partially supported by the U.S. Army Research Office through Grant No. W911NF-14-1-0103. P. A. M. C. thanks the support of a FPU MECD Grant.

Appendix A Calculation of the error $\epsilon^{\prime}$ in the Corrector step

In this appendix we want to calculate the size of the error $\epsilon^{\prime}$ we need in order to make the Corrector step successfully output a point within $\mathcal{N}(1/4)$, and its comparison with the complexity of the quantum subroutine.

To prove this suppose we define $\bm{x}$ as the concatenation of $\bm{x}^{t}$ and $\tau^{t}$, and $\bm{s}$ as concatenating $\bm{s}^{t}$ and $k^{t}$. Call $\bm{x}_{0}$ and $\bm{s}_{0}$ the exact arithmetic solution, so that being $\epsilon^{\prime}$ the error of the quantum subroutine

[TABLE]

since each entry will have an error of $\epsilon^{\prime}$ at most, and then each entry in $\bm{x}_{1}$ and $\bm{s}_{1}$ are, say, at most $1$. Using (18), we can see that

[TABLE]

and we want to calculate how small $\epsilon^{\prime}$ needs to be in order to comply with

[TABLE]

Expanding, to leading order $O(\epsilon^{\prime})$

[TABLE]

The first term is clearly our hypothesis (40), so let us calculate one of other two terms (calculating one is the same as calculating the other, they are symmetrical). Let $\bm{x}_{i}=(\bar{x}_{i,1},...,\bar{x}_{i,n+1})^{T}$ for $i\in\{0,1\}$ and similarly for $\bm{s}_{i}$, and let us expand the expression.

[TABLE]

Now we can see that the second term partially cancels out with the third

[TABLE]

Therefore, we can write

[TABLE]

Enforcing (41) can be done if we choose $\epsilon^{\prime}$ such that

[TABLE]

To leading order $O(\epsilon^{\prime})$ that means

[TABLE]

or equivalently

[TABLE]

implying

[TABLE]

Let us start analyzing the denominator. The worst case is that the denominator is very large so it forces the error to be small. The fraction $\frac{1}{4}\frac{\bm{x}_{1}^{T}\bm{s}_{0}+\bm{x}_{0}^{T}\bm{s}_{1}}{n+1}$ could turn negative, thus increasing the denominator. But we can see that its influence will be low due to its denominator $n+1$. Due to (39), we can expect $\bm{x}_{1}^{T}\bm{s}_{0}\leq(n+1)\max_{i}{|\bar{s}_{0,i}|}=O(n)$ and $\bm{x}_{0}^{T}\bm{s}_{1}\leq(n+1)\max_{i}{|\bar{x}_{0,i}|}=O(n)$. Therefore $\frac{1}{4}\frac{\bm{x}_{1}^{T}\bm{s}_{0}+\bm{x}_{0}^{T}\bm{s}_{1}}{n+1}=O(1)$.

Additionally, if all entries in $\bm{x}_{0}$ and $\bm{s}_{0}$ are $O(1)$ it is easy to check that $\left|\left|X_{0}\bm{s}_{1}\right|\right|^{2}=O(n)$ and $\left|\left|X_{1}\bm{s}_{0}\right|\right|^{2}=O(n)$. Therefore the denominator will be $O(n^{1/2})$.

The numerator of (68) is a bit more tricky. In this case the worst case is when it is small. What seems the biggest problem is that Predictor-Corrector tells us that in the exact solution $\bm{x}_{0}^{T}\bm{s}_{0}=0$, and furthermore it is divided by $n+1$. Therefore we can see that in the exact solution $\epsilon$ should be [math] (which is what we should expect). What matters is the convergence speed. At the beginning recall that $\bm{x}_{0}=\bm{s}_{0}=[1,...,1]^{T}$. Therefore, at the very beginning $\frac{\bm{x}_{0}^{T}\bm{s}_{0}}{n+1}=1$. According to Predictor-Corrector , every two iterations $\bm{x}^{T}_{0}\bm{s}_{0}$ should decrease at a rate $(1-1/(\sqrt[4]{8}\sqrt{n+1}))$. The question therefore is what happens after $O(L\sqrt{n}\bar{t})$ iterations (the iterations of the algorithm).

This finally gives us the threshold for $\epsilon^{\prime}$ we would like in order to stay inside $\mathcal{N}(1/4)$:

[TABLE]

Since $L=\log n$ in the usual case and $\bar{t}=\max[\log((\bm{x}^{0})^{T}(\bm{s}^{0})/(\epsilon_{1}\epsilon_{3}^{2})),\log(||(\bar{\bm{b}}^{T},\bar{\bm{c}}^{T})||/\epsilon_{2}\epsilon_{3})]$, the error threshold is $\epsilon^{\prime}=O(n^{-\text{poly}\log n-1/2})$.

On the other hand, the complexity on the precision of the algorithm is $O(\epsilon^{-2})$ and the difference in complexity in $n$ that we stated in tabla 2 between our algorithm and the ‘best classical’ algorithm is $O(\sqrt{n})$. This means that in order to maintain the quantum advantage one would want $\epsilon\leq O(n^{-1/4})$. Therefore we have seen that setting an $\epsilon^{\prime}$ small enough might be too expensive computationally in general for the quantum subroutine.

Appendix B Gradient descent for shifting the output of the corrector step.

We have seen that setting an $\epsilon^{\prime}$ small enough might be too expensive computationally. Therefore let us try something: choose a small precision on the error for the corrector step $\epsilon^{\prime}=O(n^{-1/a})$, and then, once we get the result, perform one step gradient descent of the size $\epsilon^{\prime\prime}$ towards $\mathcal{N}(1/4)$.

For that call, using the definition of $\mathcal{N}(\beta)$ (18),

[TABLE]

Using the equivalent to (54), we can rewrite it like

[TABLE]

Notice that the points that are in $\mathcal{N}(\beta)$ are those for which $g(x,s)\leq 0$. We calculate the gradient, calling $B=\frac{\beta^{2}+(n+1)}{(n+1)^{2}}=O(n^{-1})$:

[TABLE]

Clearly,

[TABLE]

The idea is now, supposing that

[TABLE]

if we can find $\epsilon^{\prime\prime}=O(n^{-1/a})$ such that

[TABLE]

For that, as it will be very useful, first calculate to leading order $O(\epsilon^{\prime\prime})$

[TABLE]

We can write, to leading order $O(\epsilon)$, how much is $||Xs||^{2}$ in the new point:

[TABLE]

The next thing we want to calculate is, to leading order $O(\epsilon)$, the new value of $(\bm{x}^{T}\bm{s})^{2}$

[TABLE]

We now have all the parts we need to calculate the result we were seeking. Let us try to check that we can make

[TABLE]

Let’s check it

[TABLE]

This allows us to calculate the approximate $\epsilon^{\prime\prime}$ in time $O(n)$ such that the point is in $\mathcal{N}(1/4)$. Is this $\epsilon^{\prime\prime}$ large? To answer that question expand in $\bm{x}=\bm{x}_{0}+\epsilon^{\prime}\bm{x}_{1}$ and $\bm{s}=\bm{s}_{0}+\epsilon^{\prime}\bm{s}_{1}$

[TABLE]

We know by hypothesis that the second line of (86) is less than 0 (the exact point is in $\mathcal{N}(1/4)$). Therefore we want the term of $\epsilon^{\prime\prime}$ to cancel that of $\epsilon^{\prime}$. Notice that this means that $\epsilon^{\prime\prime}=O(\epsilon^{\prime})$, since both terms in the square brackets are of size $O(n)$. Therefore, we are only introducing an error of size $O(\epsilon^{\prime})$, which is what we already had, but we can now ensure that the output is in $\mathcal{N}(1/4)$. (85) explains how to calculate $\epsilon^{\prime\prime}$ such that, to order $O(\epsilon^{\prime\prime})$ we are inside $\mathcal{N}(1/4)$; and (86) certifies that $\epsilon^{\prime\prime}=O(\epsilon^{\prime})$.

Since in the expansion we have only consider terms up to $\epsilon^{\prime\prime}$, there is the possibility that after this shift we are still outside $\mathcal{N}(1/4)$. However, in this case, we will only be $\epsilon^{\prime\prime 2}=O(\epsilon^{\prime 2})$ away from $\mathcal{N}(1/4)$, so iterating this process we can make this error $O(\epsilon^{\prime\prime 2^{i}})$, which decreases quickly on the iteration $i$. Alternatively one could take into account further terms with higher exponents of $\epsilon^{\prime\prime}$ in the expansions (81), (82) and (83).

Further, one could also wonder if the gradient step could converge towards different local minima of $g(x,s)$ than those in the central path. However, it is easy to check using (75) and (76) that all local minima (or maxima or saddle points) fulfil either $\bm{x}_{i}=0$ and $\bm{s}_{i}=0$, or

[TABLE]

for all $i\in\{1,...,n+1\}$. Substituting those points in the definition of $g(x,s)$ (74),

[TABLE]

This means that all those points are within $\mathcal{N}(\beta)$, so the gradient step will be taken towards $\mathcal{N}(\beta)$.

With this we conclude that with an $O(\epsilon^{\prime})$-shift to the output of the already $O(\epsilon^{\prime})$-precise Corrector step we can ensure that the point is in $\mathcal{N}(1/4)$.

Appendix C Convergence of the algorithm

The legitimacy of the previous procedure rests on the fact that we only have to check that we do not get out of the neighbourhood of the central path. However, a question remains: does the algorithm still converges, with this modification, in $O(n^{1/2}L)$ steps? In this appendix we analyse this question and find a positive answer. The references we will be using are Predictor-CorrectorI , mainly, and Predictor-Corrector .

Let us first review how the convergence is studied in the original article Predictor-Corrector . The convergence happens when the duality gap $\mu$ closes. In those references, one proves that in the corrector step $\mu^{k}=\mu^{k+1}$. The convergence happens in the predictor corrector, where the duality gap $\mu^{k+1}=(1-\alpha)\mu^{k}$. Lets see these two things.

For convenience define $\bm{x}^{t}$ concatenating $\bm{x}^{t}$ and $\tau^{t}$; $\bm{s}^{t}$ as the concatenation of $\bm{s}^{t}$ and $\kappa^{t}$; and similarly with $\bm{d_{x}}$, $d_{\tau}$, and $\bm{d_{s}}$ and $d_{\kappa}$. The starting equation is

[TABLE]

Since the exact solution fulfils equation (19b), that would mean that

[TABLE]

In practice, since we get an approximation solution, we can substitute the term multiplied by $\delta$ in (89) with $\gamma^{t}\mu^{t}(n+1)-(\bm{x}^{t})^{T}\bm{s}^{t}+(\bm{x}^{t})^{T}\bm{\epsilon_{x}}+(\bm{s}^{t})^{T}\bm{\epsilon_{s}}$, where $\bm{\epsilon_{x}}$ and $\bm{\epsilon_{s}}$ are the error vectors in the estimation of $\bm{d_{x}}$ and $\bm{d_{s}}$ respectively.

On the other hand, using equation (19a) for the exact solution, and the second part of theorem 5 of Predictor-Corrector , we know that $(\bm{d_{x}})^{T}\bm{d_{s}}=0$. In our case though, the equation will not be exact due to errors, and thus, the term multiplied by $\delta^{2}$ in (89) can be substituted with $\bm{d_{x}}^{T}\bm{\epsilon_{s}}+\bm{d_{s}}^{T}\bm{\epsilon_{x}}$. Therefore, in the approximate case, we can rewrite (89) for the corrector step ($\delta=1=\mu^{t}$) as

[TABLE]

From this it follows that in the worst case

[TABLE]

where $\epsilon$ is the norm of each of the entries of $\bm{\epsilon_{x}}$ and $\bm{\epsilon_{s}}$, if the entries of $\bm{x}^{t}$, $\bm{s}^{t}$, $\bm{d_{x}}$ and $\bm{d_{s}}$ are $O(1)$.

The predictor step is a bit more involved, since $\delta$ is not fixed. The calculation, though, is similar, and we get

[TABLE]

where the term $(1-\delta)\mu^{t}$ is the exact value we would recover, and $O(\epsilon)$ the part due to the error. Then question that remains is, how large can we make $\delta$ in the Predictor step? To answer this question, remember that condition for the Predictor step to end up in $\mathcal{N}(1/2)$ is

[TABLE]

We have already seen that the right hand side of (94) is $O(\epsilon)$-approximate to the exact $\frac{1}{2}\mu^{t+1}$. The left-hand side on the other hand, is

[TABLE]

where $D_{x}$ is the diagonal matrix corresponding to $\bm{d_{x}}$. In order to make this work, we will enforce a correction to make sure that (19b) is fulfilled exactly, not only approximately. This means that we have to make $X^{t}\bm{d_{x}}+S^{t}\bm{d_{s}}=-X^{t}\bm{s}^{t}$. Let us expand it with the errors (recall that $\mu^{t}=0$ for the predictor step).

[TABLE]

The terms not containing either of the error vectors cancel out. To enforce that the previous equation is fulfilled exactly, we go one by one of the entries of $\bm{x}^{t}$ and $\bm{s}^{t}$ (complexity $O(n+1)$), and choose the largest value of both. If the modulus of the entry of $\bm{x}^{t}$ is larger, then substract $\epsilon^{\prime}$ to ${\bm{d_{x}}}_{i}$ such that

[TABLE]

Clearly $\epsilon^{\prime}=O(\epsilon)$, so the introduced error will be of the size of the original one. We do equivalently for the case when the modulus of the entry of $\bm{s}^{t}$ is larger than that of $\bm{x}^{t}$. If both entries are very small, then no correction is necessary. Then, we rewrite (95)

[TABLE]

For the next steps we follow the procedure indicated in lemma 4 and theorem 1 of Predictor-CorrectorI . According to them, one chooses $\delta$ such that

[TABLE]

Recall that in the last line of (98), by hypothesis the first term is smaller than $(1-\delta)\mu^{t}/4$, because the previous step was a corrector step. With the $\delta$ from (99), (94) is fulfilled, as proved in Lemma 4 from Predictor-CorrectorI . Then we must calculate the size of $||D_{x}\bm{d_{s}}||$. Using the same methods as theorem 1 of Predictor-CorrectorI , we get

[TABLE]

for $\bm{\epsilon}$ a vector that also has entries of size $O(\epsilon)$. Thus $||\bm{\epsilon}||^{2}=(n+1)O(\epsilon)$. Using (99) we deduce that $\delta\geq 8^{-1/4}(n+1)^{-1/2}\sqrt{\mu^{t}/(\mu^{t}+O(\epsilon))}$. The fact that $\delta\geq O(n^{-1/2})$ is key to obtain a number of steps that grows as $O(n^{1/2})$ as explained in theorem 1 of Predictor-CorrectorI and theorem 6 of Predictor-Corrector . Consequently, we can see that the step size in the Predictor step gets almost the same convergence guarantees as in the exact arithmetic case. Notice that the only caveat is that we will not be able to reduce the duality gap more than $\mu=O(\epsilon)$, but that is fine. This concludes our algorithm.

Appendix D Overall structure of the algorithm

D.1 Initialization

The initialization procedure consists in preparing the matrix $M$, and the state $f$.

D.2 Termination

In the termination we propose one possible way of using Grover to run the termination explained in Predictor-Corrector . Any other classical termination is also possible.

D.3 Main loop

The main loop consists in two steps called predictor and corrector. The structure of them is very similar:

Update the data structures for $f^{t}$ and $M^{t}$. 2. 2.

Prepare $\ket{f}$ and solve $M\ket{d}=\ket{f}$ with QLSA. 3. 3.

Read $\ket{d}\rightarrow d$ and calculate the new vector $\bm{v}^{t+1}=(\bm{y}^{t+1},\bm{x}^{t+1},\tau^{t+1},\theta^{t+1},\bm{s}^{t+1},k^{t+1})$