An Inexact Feasible Quantum Interior Point Method for Linearly Constrained Quadratic Optimization

Quantum linear system algorithms (QLSAs) have the potential to speed up algorithms that rely on solving linear systems. Interior point methods (IPMs) yield a fundamental family of polynomial-time algorithms for solving optimization problems. IPMs solve a Newton linear system at each iteration to compute the search direction; thus, QLSAs can potentially speed up IPMs. Due to the noise in contemporary quantum computers, quantum-assisted IPMs (QIPMs) only admit an inexact solution to the Newton linear system. Typically, an inexact search direction leads to an infeasible solution, so, to overcome this, we propose an inexact-feasible QIPM (IF-QIPM) for solving linearly constrained quadratic optimization problems. We also apply the algorithm to ℓ1-norm soft margin support vector machine (SVM) problems, and demonstrate that our algorithm enjoys a speedup in the dimension over existing approaches. This complexity bound is better than any existing classical or quantum algorithm that produces a classical solution.


Introduction
Linearly constrained quadratic optimization (LCQO) is defined as optimizing a convex quadratic objective function over a set of linear constraints. Linear optimization is a special case of LCQO that corresponds to the case where the objective function is linear. LCQO has rich theory, algorithms, and applications. Many problems in machine learning can be formulated as LCQO problems, including variants of least square problems and variants of support vector machine training [1,2]. Some important optimization algorithms also have LCQO subproblems, e.g., sequential quadratic programming [1].
The modern age of IPMs was launched by Karmarkar's projective method for linear optimization (LO). Since then, many variants of IPMs have also been applied to nonlinear optimization problems, including LCQO problems [3,4]. Contemporary IPMs progress towards the set of optimal solutions by moving within a neighbourhood of an analytic curve known as the central path. IPMs can be categorized according to whether or not the the sequence of iterates produced by the algorithm satisfies feasibility. Feasible IPMs are initialized with a strictly feasible solution and maintain feasibility in each iteration, whereas infeasible IPMs start from an infeasible interior solution and do not require feasibility to be exactly satisfied at any point of the algorithm. For LCQO problems with n variables, feasible IPMs can produce an -approximate solution using O( √ n log(1/ )) iterations, whereas infeasible IPMs require O(n 2 log(1/ )) IPM iterations to converge to an -approximate solution [5,6].
At each IPM iteration, a linear system needs to be solved to obtain the search direction, called the Newton direction. This so-called Newton linear system is traditionally in the form of the augmented system or the normal equation system. Classically, these linear systems can be solved exactly using Bunch-Parlett factorization if the matrices in the systems are symmetric indefinite [7], or Cholesky factorization if the matrices are symmetric positive definite. Solving the Newton linear systems using direct factorization approaches requires the use of O(n 3 ) arithmetic operations, which suggests that feasible IPMs based on factoring methods cannot exhibit complexity better than O(n 3.5 log(1/ )), whereas, with the partial update, they achieve O(n 3 log(1/ )) arithmetic operation complexity. The linear systems can also be solved inexactly using some inexact methods, e.g., Krylov subspace methods, which may require fewer iterations if the desired accuracy of the solutions to the linear systems is not high. However, inaccurately solving the Newton linear systems (i.e., the inaccuracy of the search directions) may result in the infeasibility of the sequence of solutions generated by IPMs; therefore, they have only been used in infeasible IPMs.
The advent of quantum technology has led to the development of many quantumassisted algorithms for optimization and machine learning applications, such as linear regression [8] and the support vector machine training problem [9]. Following the seminal work on quantum algorithms for solving linear systems of equations [10], researchers have been studying whether QLSAs could yield quantum speedups in classical optimization algorithms. In particular, quantum IPMs (QIPMs) that utilize QLSAs to solve the Newton linear system arising in each iteration have been proposed for LO problems [11,12] and semidefinite optimization problems [13]. To maintain the feasibility of the iterates using quantum subroutines, the authors of [13,14] introduce the so-called orthogonal subspace system (OSS) for SDO and LO problems, and, in particular, demonstrate that a feasible solution to the original Newton system can be recovered from an inexact solution to the OSS. However, linearly constrained quadratic optimization problems, which are fundamental to both optimization and machine learning, have yet to be formally studied in the quantum literature.
In this work, we generalize the OSS for LO problems in [14] to LCQO problems and provide an efficient method for constructing the OSS using a quantum computer. Using the OSS, we can obtain an inexact feasible IPM, solving for the search directions inexactly but maintaining the feasibility of the iterates throughout the process of our IPM. The feasibility of the iterates gives better IPM iteration complexity and the bottleneck becomes solving the linear system, OSS. In particular, we show that a quantum implementation of our algorithm with access to quantum RAM (QRAM) obtains an -approximate solution to a given LCQO problem with worst-case complexity O n,ω, 1 √ n n ω 2 + σ max (Q) κ VAQ + n 2 , whereω = max k ω k , σ max (Q) is the maximum singular value of the Hessian of the objective function and κ VAQ is the condition number of a matrix determined by initial data; see Lemma 3. We also consider the application of 1 -norm soft margin SVM problems, in which case, an -approximate solution is obtained with complexity Here, m is the number of features and n is the number of data points.ω, Q, and κ VAQ are defined similarly from the LCQO formulation of the SVM problem; see Section 4. The dependence on dimension is better than any existing quantum or classical algorithm.
The rest of this paper is organized as follows: in Section 2, we introduce IPMs for LCQO and the OSS system; in Section 3, we discuss how to use quantum algorithms to find the Newton directions and analyze the complexity of our IF-QIPM; in Section 4, we apply our IF-QIPM to the support vector machine problem. Discussions are provided in Section 5, and some technical proofs are moved to the Appendixes A and B.

Preliminaries
In this section, we introduce notations before reviewing the theory of IPMs applied to LCQO, and derive the OSS system for the class of problems.

Notation
Vectors are typically represented by lower-case letters. We write 0 n when referring to the n-dimensional all-zeros vector, and the n-dimensional all-ones vector is denoted by e n . When the dimension is obvious from the context, we may write 0 or e, respectively. Matrices are typically represented with upper-case letters. The identity of dimension n is denoted by I n×n , and 0 n×m represents the n × m-dimensional all-zero matrix, again, dropping these subscripts when the dimension is obvious from the context. For a general n × m-dimensional matrix H, we write H i· to refer to its ith row, and, similarly, denote the jth column by H ·j . For the (i, j)th element of H, we write H ij or H i,j .
For real-valued functions f 1 , f 2 , and f 3 , we write if there exists a positive number k 4 such that f 1 ≤ k 4 f 2 . We write if there exists a positive number k 5 such that f 1 ≤ k 5 f 2 × poly log( f 3 ).

IPMs for LCQO
In this work, LCQO is defined as follows.
Definition 1 (LCQO Problem). For vectors b ∈ R m , c ∈ R n , and matrices A ∈ R m×n and Q ∈ R n×n with rank(A) = m ≤ n and Q symmetric positive semidefinite, we define the primal and dual LCQO problems as: where x ∈ R n is the vector of primal variables, and y ∈ R m , s ∈ R n are vectors of the dual variables. Problem (P) is called the primal problem and (D) is called the dual problem.
Since A is of full row-rank, A does not contain any null rows, and we further make the following assumption on matrix A.

Assumption 1.
Matrix A has no all-zero columns.

Remark 1.
Suppose that A has zero columns. Without a loss of generality, assume that the nth column is all-zero. Introducing a new variable x n+1 , we can rewrite the problem as The new LCQO problem is equivalent to the original one, and contains fewer all-zero columns. Iterating this procedure to eliminate each of the all-zero columns, we obtain a new LCQO problem satisfying Assumption 1 with no more than 2n − m variables and n constraints in the worst case.

Assumption 2.
There exists a solution (x, y, s) ∈ R n × R m × R n such that Ax = b, x > 0, A T y + s − Qx = c, and s > 0.
The set of primal-dual feasible solutions is defined as and, similarly, the set of interior feasible primal-dual solutions is given by By strong duality, the set of optimal solutions can be characterized as P D * := {(x, y, s) ∈ P D : xs = 0}, where xs denotes the Hadamard, i.e., component-wise product of x and s. Let > 0; then, the set of -approximate solutions to Problem (1) can be defined as P D := (x, y, s) ∈ P D : x T s ≤ n . ( Let X and S be diagonal matrices of x and s, respectively. Under Assumption 2, for all µ > 0, the perturbed system of optimality conditions has a unique solution (x(µ), y(µ), s(µ)), and this set of solutions gives rise to the primaldual central path CP := (x, y, s) ∈ P D 0 |x i s i = µ for i ∈ {1, . . . , n}; for µ > 0 .
IPMs apply Newton's method to solve system (3). At each iteration of infeasible IPMs, a candidate solution to the primal-dual LCQO pair in (1) is updated by solving the following linear system to find the Newton direction: where are residuals, and σ ∈ (0, 1) is the barrier reduction parameter. If r p = 0 and r d = 0, then the solution (x, y, s) exactly satisfies primal-dual feasibility. We can also define residuals in different ways as we will show later. Once the Newton direction is found, one can move along the direction but has to stay in a neighbourhood of the central path, which is defined as where θ ∈ (0, 1).
Until relatively recently, inexact solution approaches to solve the Newton linear system (4) had only been utilized in inexact infeasible IPMs (II-IPMs). For LCQO problems, ref. [6] proposes an II-IPM using an iterative method to solve the Newton systems and obtains a worst-case iteration complexity O(n 2 log( 1 )). On the other hand, feasible IPMs for LCQO problems enjoy O( √ n log( 1 )) iteration complexity [15][16][17]. In [5], the author provides a general inexact feasible IPM for LCQO problems but does not discuss how the sequence of iterates could be guaranteed to maintain primal-dual feasibility exactly when using inexact linear system solvers. This is a vital consideration, as the feasible neighborhood of the central path as outlined in (5) is a subset of the primal-dual feasible set; if primal and dual feasibility are not satisfied exactly at any point in the algorithm, the iterates leave this neighborhood and the method fails. Our work fills this gap by using a method inspired by the QIPMs of [13,14].

Orthogonal Subspaces System
Assume that (x, y, s) ∈ P D 0 . To maintain the feasibility of the primal and dual variables, the first two linear equations in system (4) need to be solved with r p = 0 and r d = 0 exactly, which can be guaranteed if ∆x lies in the null space of A, denoted as Null(A), and ∆s = Q∆x − A T ∆y. Accordingly, we can rewrite system (4) by representing ∆x as a linear combination of basis elements of Null(A). To achieve this, we partition A as Then, we construct the following matrix: Matrix V has a full column rank and satisfies AV = 0, i.e., the columns of V span the null space of A. Let ∆x = Vλ, where λ ∈ R n−m is the unknown coefficient vector used to determine ∆x. Subsequently, we can rewrite system (4) by substituting ∆x and ∆s in the third equation as A similar system was proposed and called "Orthogonal Subspaces System" (OSS) in [13,14], and we use the same name in this work. The matrix in the OSS system (6) is of size n × n, and it is nonsingular. Even if the OSS system is solved inexactly, primal and dual feasibility are preserved by computing ∆x = Vλ and ∆s = QVλ − A T ∆y. Thus, we can conclude that any inexactness will only impact the third equation of (4), i.e., r p = 0 and r d = 0. This property of the OSS system is very convenient when analyzing the proposed inexact IPM, and allows us to obtain the best known iteration complexity for IPMs.

Inexact Feasible IPM with QLSAs
In this section, we propose our IF-QIPM for LCQO problems. We begin with the IF-IPM structure introduced by [5] and describe how to quantize it into an IF-QIPM. Then, we analyze the construction of the OSS system and conclude by analyzing the overall complexity of our IF-QIPM.

IF-IPM for LCQO
In [5], the author studies a general conceptual form IF-IPM for QCLO problems by assuming the feasibility of the primal and dual iterates, which induces the following system: where r c = σµe − XSe, with σ ∈ (0, 1) being the reduction factor of the central path parameter µ, i.e., µ new = σµ. When system (7) is solved with r c = σµe − XSe inexactly yielding an error r, if r 2 ≤ δ r c 2 for some δ ∈ (0, 1), the inexact IPM converges to an -approximate solution to Problem (1) in at most O( √ n log(1/ )) iterations. As we mentioned earlier, it is not specified in [5] how to preserve primal and dual feasibility when system (7) is solved inexactly. Thus, it is presently not clear whether one could recover the convergence conditions described in [5] using inexact approaches, which are reliant on the assumption of primal-dual feasibility (see, e.g., system (7)). Now, we present a general procedure of how to solve system (7) inexactly, while the inexactness error occurs only in the third equation of system (7). Let (λ, ∆y) be an inexact solution for system (6) and r be the error at this solution, i.e., The corresponding Newton step Recall that once (λ, ∆y) is determined, then (∆x, ∆s) is also (uniquely) determined. An interesting property is that, if (λ, ∆y) and (∆x, ∆y, ∆s) can be deduced from each other, then the OSS system and system (7) yield the same error term r. Hence, the convergence conditions built upon system (7) can be directly examined using the residual r c and error r of the OSS system. Let OSS be the target accuracy of the OSS system (6), i.e., where (λ * , ∆y * ) is the accurate solution. According to [5], in order to guarantee that the IF-IPM converges, we must have where δ ∈ (0, 1) is a constant parameter. Therefore, to ensure the convergence of the IF-IPM, it suffices to set The IF-IPM is presented in full detail in Algorithm 1. In each iteration, we build and solve system (6) classically. We solve system (6) to the accuracy just introduced above and then compute the feasible Newton step from the inexact solution and take a full Newton step.
In the quantum-assisted IF-IPM, or IF-QIPM, we propose accelerating Step 7 using quantum subroutines. In the next sections, we investigate how to use quantum algorithms to build and solve the OSS system and obtain the Newton direction.

IF-QIPM for LCQO
The pseudocode of our IF-QIPM is presented in Algorithm 2. At each iteration of the IF-QIPM, we construct and solve system (6) and compute the Newton direction using quantum algorithms. To obtain an OSS -approximate solution of system (6), we first block encode system (8); see Appendix A. Then, we use quantum algorithms to solve for an QLSA -approximate solution of system (8). This solution is normalized but we can rescale it to obtain an OSS -approximate solution of system (6). Details are discussed later in this section.
(λ k , ∆y k ) ← solve system (6) with accuracy k OSS quantumly 8: ∆x k = Vλ k and ∆s k = −A T ∆y k 9: Here, θ 0 < 1 and its value will be discussed later. First, we introduce some notations to simplify the OSS system. In the kth iteration of Algorithm 2, let Then, the OSS system can be rewritten as As discussed in [14], to solve the OSS system (6) using quantum algorithms, we can first rewrite it as the normalized Hermitian OSS system To use the QLSAs mentioned earlier, we need to turn the linear system (8) into a quantum linear system using the block encoding introduced in [18]. To this end, we first decompose the coefficient matrix in linear system (8) as where To compute matrix V, we need to find a basis matrix A B of matrix A and we need to compute the inverse matrix A −1 B . Both steps are nontrivial and can be expensive. However, we can reformulate the LCQO problem as follows: In this case, we have an obvious basis A B = I 0 0 I and matrix V can be constructed efficiently Since matrix A has no all-zero rows, matrix V has no all-zero rows either. This property of the reformulation is useful in the analysis of the proposed IF-QIPM but we do not want to build the complexity analysis on the reformulated problem. Thus, without a loss of generality we may make the following assumption.

Assumption 3. Matrix A is of the form A = I A N .
To simplify the analysis, we further assume that the input data are integers. Based on the two assumptions above, we have the following lemma.
where V i· and (A N ) i· are the ith row of V and A N , respectively. Now, we are ready to give θ 0 in our definition of the central path neighborhood; see (5). We set We also define ω k as the maximum of the values of primal variables and dual slack variables in the kth iteration.

Definition 2.
Let (x k , y k , s k ) be a candidate solution for Problem (1); then, As is standard in the literature on quantum algorithms, in this work, we assume access to quantum random access memory (QRAM). Then, Step 7 of Algorithm 2 consists of three parts: (1). use block encoding to build system (8); (2). use QLSAs to solve system (8); (3). use quantum tomography algorithms (QTAs) to extract the classical solution. We use the block-encoding methods introduced in [18] to block-encode linear system (8).

Proposition 1.
In the kth iteration of Algorithm 2, using the block-encoding methods introduced in [18] and the decomposition described in Equations (9) and (10), a -block-encoding of the matrix in system (8) can be implemented efficiently and the complexity will be dominated by the complexity of the QLSA step. Here, QLSA is the accuracy required for the QLSA step and κ M k is the condition number of matrix M k .
Proof. See Appendix A for proof.
Provided access to QRAM, the complexity associated with block encoding the OSS system coefficient matrix and preparing a quantum state encoding the right hand side amounts to polylogarithmic overhead. The cost of these steps is therefore negligible when compared with the complexity contributed by QLSAs and QTAs, so we ignore it here. To bound the total complexity contributed by QLSAs and QTAs, we first need to analyze the accuracy of QLSA characterized by QLSA , the accuracy of QTA characterized by QTA , and their relationship.
In each iteration, we use a QLSA to solve the block-encoded version of system (8) and obtain an QLSA -approximate solution. Then, we use a QTA to extract an QTAapproximate solution from the quantum machine. In the context of QLSAs and QTAs, ifz is an -approximate solution of z, thenz satisfies Observe that this definition of accuracy differs from the concept of -approximate solutions defined in (2). Similar to [12,13], the QLSA we use is proposed by [19] and the QTA we use is proposed by [20]. Following the argument in Section 2 in [12], we can establish the relationship among QLSA , QTA , and k OSS as where k OSS is defined as the 2 norm of the residual when solving system (8) inexactly in the kth iteration. This coefficient is also used to rescale the solution. According to [12], we rescale the normalized solution obtained from QLSA and QTA by to obtain the k OSS -approximate solution for system (6). Here, we did not add superscript to QLSA and QTA , and the reason shall be revealed later. Let 0 k z k be an inexact solution for system (8) in the kth iteration. Then, the norm of residual of system (8), which is k OSS , and the norm of residual of system (6), which is M kzk − r k c 2 , satisfies Recall that the error arising from the OSS system (6) is the same as the error in the full Newton system (7); then, we can directly use the convergence condition in [5], i.e., We can require and it follows that ensures the convergence of the IF-QIPM. The complexities for each step are also available now. Using the QLSA from [19] and QTA from [20], we have the complexity for QLSA and QTA: Since we have QTA = δ 2 and δ ∈ (0, 1) is a constant parameter, we omit QTA in the Big-O notation. Note that the complexity of the block-encoding procedure is dominated by that of QLSA and QTA and thus we ignore the complexity contributed by block encoding. In Step 8, the complexity contributed by computing Newton step from OSS solution is O(n 2 ). The total complexity for the kth iteration of IF-QIPM will be 3.2.1. Bound for ω k / M k F In this section, all of the quantities that we consider are from the kth iteration. For simplicity, we omit the superscript k in this section unless we need it. Using the property of trace, we have Recalling the central path neighborhood that we defined in (5), we define a matrix E such that It is obvious that E is a diagonal matrix and satisfies With this, we can have tr XQVV T S = tr SXQVV T = tr (θµE + µI)QVV T = tr θµEQVV T + tr µQVV T .
For the second term, we know that Q and V T QV are both positive semidefinite. Thus, we can have tr QVV T = tr V T QV ≥ 0 because of the cyclic invariant property of trace. According to the Cauchy-Schwarz inequality, we have Thus, we have where the last inequality holds due to condition (11). Thus, we can bound M F by Since XQVV T QX 0, we have Since X and S are both positive diagonal matrices, we have As we mentioned in the very beginning of this section, at each iteration, ω is indeed ω k , but the superscript is ignored here. Now, we aim to find a bound for µ so we can further bound M 2 F . Since ω is the upper bound for the magnitude of the primal and dual slack variables, we have Recall the definition of matrix E; see (14). Thus, we have Thus, where the last inequality follows from the bound for θ; see (11). Thus, we have

Bound for κ M k
Similar to the previous section, we ignore the superscript k unless we need it. We will start with a general result and then work on the matrix M k . The following lemma is a well-known result regarding condition numbers of matrices and can be proven using Courant-Fischer-Weyl min-max principle [21].

Lemma 2.
For any full row rank matrix P ∈ R m×n and symmetric positive definite matrix D ∈ R n×n , their condition number satisfies Next, we analyze the matrix in the OSS system (8). Specifically, we focus on M T M since we are interested in the spectral property of the OSS system (8). Using the matrix E defined in (14), we have the following decomposition: The second equality holds because as AV = 0 and Q is symmetric. Then, plugging (14) into the first diagonal block of the decomposition we obtained earlier, we have The first two matrices are nonsingular, so we can apply the Lemma 2, and thus we only need to study the middle matrix. Denote the middle matrix by Ψ. Observe that Ψ is almost the same as its counterpart in [14]. Subsequently, we have the following result regarding the spectral property of M k .
where κ VAQ is the condition number of the matrix Putting all of these together, we have the complexity for our IF-QIPM for LCQO problems.
Theorem 1. The IF-QIPM for LCQO problems stops with the final duality gap less than in at most O √ n log(1/ ) IPM iterations and, in each IPM iteration, the Newton direction can be obtained with complexity O n,ω, 1 n ω 2 + σ max (Q) κ VAQ + n 2 , whereω = max k ω k .
Proof. The complexity bound for the IPM iterations comes from the result in [5]. According to (13), the complexity for obtaining the Newton direction is Combining this with the result in Sec. 3.2.1, the bound in Lemma 3, and µ k ≥ , we have O n,ω, 1 nω k κ M k M k F + n 2 = O n,ω, 1 n ω 2 + σ max (Q) κ VAQ + n 2 .

Application in Support Vector Machine Problems
In this section, we discuss how to use our IF-QIPM to solve SVM problems. We show that our algorithm can solve 1 -norm soft margin SVM problems faster than any existing classical or quantum algorithms with respect to dimension.
The ordinary SVM problem works on a linearly separable dataset, in which the data points have binary labels. The ordinary SVM aims to find a hyperplane correctly separating the data points with a maximum margin. However, in practice, the data points are not necessarily linearly separable. To allow for mislabelling, the concept of a soft margin SVM was introduced in [22]. Let {(φ i , ζ i ) ∈ R m × {−1, +1}|i = 1, . . . , n} be the set of data points, Φ be a matrix with the ith column being φ i , and Z be a diagonal matrix with the ith diagonal element being ζ i . The SVM problem with an l 1 -norm soft margin can be formulated as below. min (ξ,w,t)∈R n ×R m ×R Here, (w, t) determines a hyperplane and C is a penalty parameter. In [9], the authors rewrote the SVM problem as a second-order conic optimization (SOCO) problem and used the quantum algorithm that they proposed to solve the resulting SOCO problem. They claim the complexity of their algorithm has O(n 2 ) dependence on the dimension, which is better than any classical algorithm. However, the algorithm in [9] is invalid. Their algorithm is an inexact infeasible-QIPM (II-QIPM), while they used the IPM complexity for the feasible-QIPM, which ignores at least O(n 1.5 ) dependence on n. They also missed the symmetrization of the Newton step, which is necessary for SOCO problems and makes their Newton step invalid.
Aside from [9], some pure quantum algorithms for SVM problems are also proposed. In [23], the authors propose a pure quantum algorithm for SVM problems. They claim the complexity is O(κ 3 eff −3 log(mn)), where κ eff is the condition number of a matrix involving the kernel matrix and is the accuracy. In the worst case, κ eff = O(m). Their complexity is worse than ours regarding the dependence of dimension and accuracy. In addition, their algorithm does not provide classical solutions. Namely, the solution is in the quantum machine and we cannot read or use it in a classical computer. However, our algorithm produces a classical solution.
To convert the problem into standard-form LCQO, we introduce (w + , w − ) ∈ R m + × R m + , (t + , t − ) ∈ R + × R + , and a slack variable ρ ∈ R n + . Then, we can obtain the following formulation: This is a standard-form LCQO problem with non-negative variables (w + , Thus, we can use the proposed IF-QIPM for LCQO problems to solve the 1 -norm soft margin SVM problems and obtain an -approximate solution with complexity O m,n,ω, 1 (m + n) 1.5 ω 2 This dependence on dimension is better than any existing quantum or classical algorithm.

Discussion
In this work, we present an IF-QIPM for LCQO problems by combining the IF-IPM framework proposed in [5] and the OSS system introduced in [14]. Our algorithm has n 1.5 dependence on n, which is better than any existing algorithms for LCQO problems. The dependence on the accuracy is polynomial, which is worse than classic IPMs. Iterative refinement techniques might help to improve the dependence on the accuracy but they are beyond the discussion of this work.

Conflicts of Interest:
The funder had no role in the design of the study; in the writing of the manuscript; or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Block Encoding of the OSS System
In this section, we ignore the superscript k for simplicity. As described in Equation (9), we first block encode each of the matrices involved in (10). We assume that V, A, S, and X are given and are stored in a quantum accessible data structure (we ignore the complexity to store the classical information into the quantum machine). For the first matrix -block-encoding of M 1 can be implemented efficiently according to Lemma 50 from [18]. The second matrix Then, we can block encode the two matrices first, and then apply a linear combination to obtain M 3 . In fact, a ( Q F , O(poly log(n)), 3 ) -block-encoding of the left matrix can be implemented efficiently according to Lemma 50 from [18] and a (1, O(poly log(n)), 3 ) -block-encoding of the right matrix can be implemented efficiently according to Lemma 48 in [18]. With the state-preparation cost of the linear combination coefficient vector (1, 1) neglected, a ( Q F + 1, O(poly log(n)), ( Q F + 1) 3 ) -block-encoding of M 3 can be implemented efficiently according to Lemma 52 from [18]. The fourth matrix is one-row-sparse and two-columns-sparse. After being scaled by 1 ω , each element of M 4 /ω has an absolute value of at most 1. According to Lemma 48 in [18], a -block-encoding can be implemented efficiently according to Lemma 53 from [18]. For the linear combination M 2 /ω + M 3 M 4 /ω, the cost for the state-preparation of the coefficient vector (1, 1) is negligible and thus a -block-encoding can be implemented efficiently according to Lemma 52 from [18]. For the matrix multiplication of O(poly log(n)), -block-encoding can be implemented efficiently according to Lemma 53 from [18].
Finally, considering that the complexity of the state-preparation of the vector can be neglected, a O(poly log(n)), -block-encoding of the coefficient matrix of system (8) can be implemented efficiently according to Lemma 52 from [18]. We can choose where K depends on the initial data Now, considering that the complexity for all of the block-encoding algorithms that we have used so far has poly-logarithmic dependence on the dimension and accuracy, and that, for i = 1, 2, 3, 4 O poly log( 1 i ) = O(poly log(κ M )), the complexity for block encoding will be dominated by the complexity for QLSA because QLSA has linear dependence on κ M , we can ignore the complexity of block encoding.

Appendix B. Spectral Analysis for Matrix Ψ
In this section, we provide the spectral analysis for the matrix Just like in the previous section, for simplicity, we ignore the superscript k. We can perform the following decomposition: Let us use the following notation: It can be proven that Ψ 1 is positive definite. The majority of the proof of this conclusion comes from the paper [14]. For the reader's convenience, we provide the complete proof here.
Matrix Ψ 1 is a block diagonal matrix, with all four blocks being diagonal matrices. Thus, we can easily compute the eigenvalues using the characteristic polynomial Clearly, det(Ψ 1 − qI) = 0 gives n quadratic equations and each quadratic equation gives two eigenvalues. The two eigenvalues from the ith quadratic equation are Recalling the definition of E in (14), we can write One can verify that the square root always exists because With θ ∈ 0, min 1 3 √ n , 1 4 QVV T F +1 , we have This means that matrix Ψ 1 is positive definite and its eigenvalues coincide with its singular values because Ψ 1 is also real and symmetric. Analogously, we have Thus, the condition number of Ψ satisfies κ(Ψ) ≤ σ max (Ψ 1 ) + σ max (Ψ 2 ) σ min (Ψ 1 ) + σ min (Ψ 2 ) where the last inequality comes from the definition of ω. Since ω 2 ≥ x i s i ≥ (1 − θ)µ, we have κ(Ψ) = O ω 2 (ω 2 + µσ max (Q)) µ 2 .
Using Lemma 2, we can also bound the condition number of matrix M by