Quantum Linear System Solver Based on Time-optimal Adiabatic Quantum Computing and Quantum Approximate Optimization Algorithm

Linear system solvers are used ubiquitously in scientific computing. Quantum algorithms for solving large systems of linear equations, also called the quantum linear system problem (QLSP), have received much attention recently [8, 10, 11, 12, 16, 17, 30, 33, 34]. The goal of QLSP is to efficiently compute \( \mathinner {|{x}\rangle }=A^{-1}\mathinner {|{b}\rangle }/{\parallel }{A^{-1}\mathinner {|{b}\rangle }}_2 \) on a quantum computer, where \( A\in {\mathbb {C}}^{N\times N} \) , and \( \mathinner {|{b}\rangle }\in {\mathbb {C}}^N \) is a normalized vector (for simplicity we assume \( N=2^n \) , and \( \Vert A\Vert _2 = 1 \) ). The ground-breaking Harrow, Hassidim, and Lloyd (HHL) algorithm obtains \( \mathinner {|{x}\rangle } \) with cost \( {\mathcal {O}}(\text{poly}(n) \kappa ^2 /\epsilon) \) , where \( \kappa = \Vert A\Vert \Vert A^{-1}\Vert \) is the condition number of A, and \( \epsilon \) is the target accuracy. On the other hand, the best classical iterative algorithm is achieved by the conjugate gradient method, where the cost is at least \( {\mathcal {O}}(N \sqrt {\kappa }\log (1/\epsilon)) \) , with the additional assumptions that A should be Hermitian positive definite and a matrix-vector product can be done with \( {\mathcal {O}}(N) \) cost [29]. The complexity of direct methods based on the Gaussian elimination procedure removes the dependence on \( \kappa \) , but the dependence on N is typically super-linear even for sparse matrices [21]. Therefore, the HHL algorithm can potentially be exponentially faster than classical algorithms with respect to N. The undesirable dependence with respect to \( \epsilon \) is due to the usage of the quantum phase estimation (QPE) algorithm. Recent progresses based on linear combination of unitaries (LCU) [12] and quantum signal processing (QSP) [16, 22] have further improved the scaling to \( {\mathcal {O}}(\kappa ^2 \text{poly}(\log (\kappa /\epsilon))) \) under different query models, without using QPE. However, the \( {\mathcal {O}}(\kappa ^2) \) scaling can be rather intrinsic to the methods, at least before complex techniques such as variable time amplitude amplification (VTAA) algorithm [2] are applied.

The VTAA algorithm is a generalization of the conventional amplitude amplification algorithm, and allows to quadratically amplify the success probability of quantum algorithms in which different branches stop at different time. In [2], VTAA was first used to successfully improve the complexity of HHL algorithm to \( \widetilde{{\mathcal {O}}}(\kappa /\epsilon ^3) \) . In [12], the authors further combine VTAA algorithm and a low-precision phase estimate to improve the complexity of LCU to \( \widetilde{{\mathcal {O}}}(\kappa \, \text{poly}(\log (\kappa /\epsilon))) \) , which is near-optimal with respect to both \( \kappa \) and \( \epsilon \) . It is worth noting that the VTAA algorithm is a complicated procedure and can be difficult to implement. Thus, it remains of great interest to obtain alternative algorithms to solve QLSP with near-optimal complexity scaling without resorting to VTAA.

Some of the alternative routes for solving QLSP are provided by the adiabatic quantum computing (AQC) [1, 19] and a closely related method called the randomization method (RM) [6, 30]. The key idea of both AQC and RM is to solve QLSP as an eigenvalue problem with respect to a transformed matrix. Assume that a Hamiltonian simulation can be efficiently performed on a quantum computer, it is shown that the runtime of RM scales as \( {\mathcal {O}}(\kappa \log (\kappa)/\epsilon) \) [30], which achieves near-optimal complexity with respect to \( \kappa \) without using VTAA algorithm as a subroutine. The key idea of the RM is to approximately follow the adiabatic path based on the quantum Zeno effect (QZE) using a Monte Carlo method. Although RM is inspired by AQC, the runtime complexity of the (vanilla) AQC is at least \( {\mathcal {O}}(\kappa ^2/\epsilon) \) [1, 7, 30]. Therefore, the RM is found to be at least quadratically faster than AQC with respect to \( \kappa \) .

In this article, we find that with a simple modification of the scheduling function to traverse the adiabatic path, the gap between AQC and RM can be fully closed, along with the following two aspects. (1) We propose a family of rescheduled AQC algorithms called AQC(p). Assuming \( \kappa \) (or its upper bound) is known, we demonstrate that for any matrix A (possibly non-Hermitian or dense), when \( 1\lt p\lt 2 \) , the runtime complexity of AQC(p) can be only \( {\mathcal {O}}(\kappa /\epsilon) \) . Thus, AQC(p) removes a logarithmic factor with respect to \( \kappa \) compared to RM. (2) We propose another rescheduled algorithm called AQC(exp), of which the runtime is \( {\mathcal {O}}(\kappa ~\text{poly}(\log (\kappa /\epsilon))) \) . The main benefit of AQC(exp) is the improved dependence with respect to the accuracy \( \epsilon \) , and this is the near-optimal complexity (up to logarithmic factors) with respect to both \( \kappa \) and \( \epsilon \) . The scheduling function of AQC(exp) is also universal because we do not even need the knowledge of an upper bound of \( \kappa \) . Existing works along this line [15, 25] only suggest that runtime complexity is \( {\mathcal {O}}(\kappa ^3~\text{poly}(\log (\kappa /\epsilon))) \) , which improves the dependence with respect to \( \epsilon \) at the expense of a much weaker dependence on \( \kappa \) . Our main technical contribution is to again improve the dependence on \( \kappa . \) Since the cost of any generic QLSP solver can not be less than \( {\mathcal {O}}(\kappa) \) [17], our result achieves the near-optimal complexity up to logarithmic factors. We remark that in the AQC based algorithm, only the total runtime T depends on \( \kappa \) .

Beyond the runtime complexity, we also discuss gate-efficient approaches to implement our AQC(p) and AQC(exp) methods. In particular, assume that we are given access to the same query models as those in [12]: the sparse input model of a d-sparse matrix A and the prepare oracle of the state \( \mathinner {|{b}\rangle } \) . We demonstrate that, when the adiabatic dynamics is simulated using the truncated Dyson series method [23], the query complexity of the AQC(p) method scales \( {\mathcal {O}}(d\kappa /\epsilon \log (d\kappa /\epsilon)) \) , and that of the AQC(exp) method scales \( {\mathcal {O}}(d\kappa ~\text{poly}\log (d\kappa /\epsilon)) \) . Both algorithms scale almost linearly in terms of \( \kappa \) , and the AQC(exp) method can achieve near-optimal scaling in both \( \kappa \) and \( \epsilon \) . Furthermore, the asymptotic scaling of the AQC(exp) method is the same as that of LCU with VTAA method [Theorem 5][12]. However, the AQC(exp) method avoids the usage of complex VTAA routine, which significantly simplifies its practical implementation.

The quantum approximate optimization algorithm (QAOA) [14], as a quantum variational algorithm, has received much attention recently thanks to the feasibility of being implemented on near-term quantum devices. Due to the natural connection between AQC and QAOA, our result immediately suggests that the time-complexity for solving QLSP with QAOA is also at most \( {\mathcal {O}}(\kappa ~\text{poly}(\log (\kappa /\epsilon))) \) , which is also confirmed by numerical results. We also remark that both QAOA and AQC schemes prepare an approximate solution to the QLSP in a pure state, while RM prepares a mixed state. Moreover, all methods above can be efficiently implemented on gate-based computers and are much simpler than those using the VTAA algorithm as a subroutine.

Assume \( A\in {\mathbb {C}}^{N\times N} \) is an invertible matrix with condition number \( \kappa \) and \( \Vert A\Vert _2 = 1 \) . Let \( \mathinner {|{b}\rangle }\in {\mathbb {C}}^{N} \) be a normalized vector. Given a target error \( \epsilon \) , the goal of QLSP is to prepare a normalized state \( \mathinner {|{x_{\text{a}}}\rangle } \) , which is an \( \epsilon \) -approximation of the normalized solution of the linear system \( \mathinner {|{x}\rangle }=A^{-1}\mathinner {|{b}\rangle }/{\parallel }{A^{-1}\mathinner {|{b}\rangle }}_2 \) , in the sense that \( \Vert \mathinner {|{x_{\text{a}}}\rangle }\mathinner {\langle {x_{\text{a}}}|}-\mathinner {|{x}\rangle }\mathinner {\langle {x}|}\Vert _2 \le \epsilon \) .

For simplicity, we first assume A is Hermitian and positive definite and will discuss the generalization to non-Hermitian case later.

The first step to design an AQC-based algorithm for solving QLSP is to transform the QLSP to an equivalent eigenvalue problem. Here, we follow the procedure introduced in [30]. Let \( Q_{b}=I_N-\mathinner {|{b}\rangle }\mathinner {\langle {b}|} \) . We introducethen \( H_0 \) is a Hermitian matrix and the null space of \( H_0 \) is \( \text{Null}(H_0)=\text{span}\lbrace \mathinner {|{{\widetilde{b}}}\rangle },\mathinner {|{\bar{b}}\rangle }\rbrace \) . Here, \( \mathinner {|{{\widetilde{b}}}\rangle }=\mathinner {|{0,b}\rangle }:=(b,0)^{\top },\mathinner {|{\bar{b}}\rangle }=\mathinner {|{1,b}\rangle }:=(0,b)^{\top } \) . The dimension of \( H_0 \) is \( 2N \) and one ancilla qubit is needed to enlarge the matrix block. We also defineHere, \( \sigma _{\pm }=\frac{1}{2}(\sigma _x\pm {\mathrm{i}}\sigma _y) \) . Note that if \( \mathinner {|{x}\rangle } \) satisfies \( A\mathinner {|{x}\rangle }\propto \mathinner {|{b}\rangle } \) , we have \( Q_bA\mathinner {|{x}\rangle }=Q_b\mathinner {|{b}\rangle }=0 \) . Then \( \text{Null}(H_1)=\text{span}\lbrace \mathinner {|{{\widetilde{x}}}\rangle },\mathinner {|{\bar{b}}\rangle }\rbrace \) with \( \mathinner {|{{\widetilde{x}}}\rangle }=\mathinner {|{0,x}\rangle } \) . Since \( Q_b \) is a projection operator, the gap between 0 and the rest of the eigenvalues of \( H_0 \) is 1. The gap between 0 and the rest of the eigenvalues of \( H_1 \) is bounded from below by \( 1/\kappa \) (see Appendix A).

QLSP can be solved if we can prepare the zero-energy state \( \mathinner {|{{\widetilde{x}}}\rangle } \) of \( H_1 \) , which can be achieved by the AQC approach. Let \( H(f(s)) = (1-f(s))H_0 + f(s)H_1, 0\le s\le 1 \) . The function \( f:[0,1]\rightarrow [0,1] \) is called a scheduling function, and is a strictly increasing mapping with \( f(0) = 0, f(1) = 1 \) . The simplest choice is \( f(s)=s \) , which gives the “vanilla AQC”. We sometimes omit the s-dependence as \( H(f) \) to emphasize the dependence on f. Note that for any s, \( \mathinner {|{\bar{b}}\rangle } \) is always in \( \text{Null}(H(f(s))) \) , and there exists a state \( \mathinner {|{{\widetilde{x}}(s)}\rangle }=\mathinner {|{0,x(s)}\rangle } \) , such that \( \text{Null}(H(f(s)))=\lbrace \mathinner {|{{\widetilde{x}}(s)}\rangle },\mathinner {|{\bar{b}}\rangle }\rbrace \) . In particular, \( \mathinner {|{{\widetilde{x}}(0)}\rangle }=\mathinner {|{{\widetilde{b}}}\rangle } \) and \( \mathinner {|{{\widetilde{x}}(1)}\rangle }=\mathinner {|{{\widetilde{x}}}\rangle } \) ; and therefore, \( \mathinner {|{{\widetilde{x}}(s)}\rangle } \) is the desired adiabatic path. Let \( P_0(s) \) be the projection to the subspace \( \text{Null}(H(f(s))) \) , which is a rank-2 projection operator \( P_0(s)=\mathinner {|{{\widetilde{x}}(s)}\rangle }\mathinner {\langle {{\widetilde{x}}(s)}|}+\mathinner {|{\bar{b}}\rangle }\mathinner {\langle {\bar{b}}|} \) . Furthermore, the eigenvalue 0 is separated from the rest of the eigenvalues of \( H(f(s)) \) by a gapWe refer to the Appendix A for the derivation.

Consider the adiabatic evolutionwhere \( 0 \le s \le 1 \) , and the parameter T is called the runtime of AQC. The quantum adiabatic theorem [Theorem 3][19] states that for any \( 0\le s\le 1 \) ,whereThe derivatives of H are taken with respect to s, i.e., \( H^{(k)}(s) := \frac{d^k}{ds^k} H(f(s)), k = 1,2 \) . Throughout the article, we shall use a generic symbol C to denote constants independent of \( s,\Delta ,T \) .

Intuitively, the quantum adiabatic theorem in Equation (3) says that, if the initial state is an eigenstate corresponding to the eigenvalue 0, then for large enough T the state \( \mathinner {|{\psi _T(s)}\rangle } \) will almost stay in the eigenspace of \( H(s) \) corresponding to the eigenvalue 0, where there is a double degeneracy and only one of the eigenstate \( \mathinner {|{{\widetilde{x}}(s)}\rangle } \) is on the desired adiabatic path. However, such degeneracy will not break the effectiveness of AQC for the following reason. Note that \( \mathinner {\langle {\bar{b}|\psi _T(0)}\rangle }=0 \) , and \( H(f(s))\mathinner {|{\bar{b}}\rangle }=0 \) for all \( 0\le s\le 1 \) , so the Schrödinger dynamics Equation (2) implies \( \mathinner {\langle {\bar{b}|\psi _T(s)}\rangle }=0 \) , which prevents any transition of \( \mathinner {|{\psi _T(s)}\rangle } \) to \( \mathinner {|{\bar{b}}\rangle } \) . Therefore, the adiabatic path will stay along \( \mathinner {|{{\widetilde{x}}(s)}\rangle } \) . Using \( \mathinner {\langle {\bar{b}|\psi _T(s)}\rangle }=0 \) , we have \( P_0(s)\mathinner {|{\psi _T(s)}\rangle }=\mathinner {|{{\widetilde{x}}(s)}\rangle }\mathinner {\langle {{\widetilde{x}}(s)|\psi _T(s)}\rangle } \) . Therefore, the estimate in Equation (3) becomesThis also implies that (see appendix)Therefore \( \eta (1) \) can be an upper bound of the distance of the density matrix.

If we simply assume \( {\parallel }{H^{(1)}}_2,{\parallel }{H^{(2)}}_2 \) are bounded by constants, and use the worst case bound that \( \Delta \ge \kappa ^{-1} \) , we arrive at the conclusion that in order to have \( \eta (1)\le \epsilon \) , the runtime of vanilla AQC is \( T\gtrsim \kappa ^3/\epsilon \) .

Our goal is to reduce the runtime by choosing a proper scheduling function. The key observation is that the accuracy of AQC depends not only on the gap \( \Delta (f(s)) \) but also on the derivatives of \( H(f(s)) \) , as revealed in the estimate in Equation (4). Therefore, it is possible to improve the accuracy if a proper time schedule allows the Hamiltonian \( H(f(s)) \) to slow down when the gap is close to 0. We consider the following schedule [1, 19]:Here, \( \Delta _* \) is defined in Equation (1), and \( c_p = \int _0^1 \Delta _*^{-p}(u) du \) is a normalization constant chosen so that \( f(1) = 1 \) . When \( 1 \lt p \le 2 \) , Equation (5) can be explicitly solved asNote that as \( s\rightarrow 1 \) , \( \Delta _{*}(f(s))\rightarrow \kappa ^{-1} \) , and therefore the dynamics of \( f(s) \) slows down as \( f\rightarrow 1 \) and the gap decreases towards \( \kappa ^{-1} \) . We refer to the adiabatic dynamics Equation \( (2) \) with the schedule in Equation \( (5) \) as the AQC(p) scheme. Our main result is given in Theorem 1 (See Appendix D for the proof).

The runtime complexity of the AQC(p) method with respect to \( \kappa \) is only \( {\mathcal {O}}(\kappa) \) . Compared to Reference [30], AQC(p) further removes the \( \log (\kappa) \) dependence when \( 1\lt p\lt 2 \) , and hence reaches the optimal complexity with respect to \( \kappa \) . Interestingly, though not explicitly mentioned in [30], the success of RM for solving QLSP relies on a proper choice of the scheduling function, which approximately corresponds to AQC(p=1). It is this scheduling function, rather than the QZE or its Monte Carlo approximation per se that achieves the desired \( {\mathcal {O}}(\kappa \log \kappa) \) scaling with respect to \( \kappa \) . Furthermore, the scheduling function in RM is similar to the choice of the schedule in the AQC(p=1) scheme. The speedup of AQC(p) versus the vanilla AQC is closely related to the quadratic speedup of the optimal time complexity of AQC for Grover’s search [1, 19, 27, 28], in which the optimal time scheduling reduces the runtime from \( T\sim {\mathcal {O}}(N) \) (i.e., no speedup compared to classical algorithms) to \( T\sim {\mathcal {O}}(\sqrt {N}) \) (i.e., Grover speedup). In fact, the choice of the scheduling function in Reference [28] corresponds to AQC(p=2) and that in Reference [19] corresponds to AQC(1<p<2).

Although AQC(p) achieves the optimal runtime complexity with respect to \( \kappa \) , the dependence on \( \epsilon \) is still \( {\mathcal {O}}(\epsilon ^{-1}) \) , which limits the method from achieving high accuracy. It turns out that when T is sufficiently large, the dependence on \( \epsilon \) could be improved to \( {\mathcal {O}}(\text{poly} \log (1/\epsilon)) \) , by choosing an alternative scheduling function.

The basic observation is as follows. In the AQC(p) method, the adiabatic error bound we consider, i.e., Equation (4), is the so-called instantaneous adiabatic error bound, which holds true for all \( s\in [0,1] \) . However, when using AQC for solving QLSP, it suffices only to focus on the error bound at the final time \( s=1 \) . It turns out that this allows us to obtain a tighter error bound. In fact, such an error bound can be exponentially small with respect to the runtime [1, 15, 25, 32]. Roughly speaking, with an additional assumption for the Hamiltonian \( H(f(s)) \) that the derivatives of any order vanish at \( s=0,1 \) , the adiabatic error can be bounded by \( c_1\exp (-c_2T^{\alpha }) \) for some positive constants \( c_1, c_2, \alpha \) . Furthermore, it is proved in [15] that if the target eigenvalue is simple, then \( c_1 = {\mathcal {O}}(\Delta _*^{-1}) \) and \( c_2 = {\mathcal {O}}(\Delta _*^3) \) . Note that \( \Delta _* \ge \kappa ^{-1} \) for QLSP, thus, according to this bound, to obtain an \( \epsilon \) -approximation, it suffices to choose \( T = {\mathcal {O}}(\kappa ^3~\text{poly}(\log (\kappa /\epsilon))) \) . This is an exponential speedup with respect to \( \epsilon \) , but the dependence on the condition number becomes cubic again.

However, it is possible to reduce the runtime if the change of the Hamiltonian is slow when the gap is small, as we have already seen in the AQC(p) method. For QLSP, the gap monotonically decreases, and the smallest gap occurs uniquely at the final time, where the Hamiltonian \( H(s) \) can be set to vary slowly by requiring its derivatives to vanish at the boundary.

We consider the following schedulewhere \( c_e = \int _0^1\exp \left(-1/(s^{\prime }(1-s^{\prime }))\right){\mathrm{d}} s^{\prime } \) is a normalization constant such that \( f(1) = 1 \) . This schedule can assure that \( H^{(k)}(0) = H^{(k)}(1) = 0 \) for all \( k\ge 1 \) . We refer to the adiabatic dynamics Equation (2) with the schedule in Equation (8) as the AQC(exp) scheme. Our main result is given in Theorem 2 (see Appendix E for the proof).

Compared with RM and AQC(p), although the \( \log (\kappa) \) dependence reoccurs, AQC(exp) achieves an exponential speedup over RM and AQC(p) with respect to \( \epsilon \) (and hence giving its name), and thus is more suitable for preparing the solution of QLSP with high fidelity. Furthermore, the time scheduling of AQC(exp) is universal and AQC(exp) does not require knowledge on the bound of \( \kappa \) .

We remark that the performance of the AQC(exp) method is sensitive to the perturbations in the scheduling function, which can affect the final error in the AQC(exp) method. This is similar to the finite precision effect in the adiabatic Grover search reported in [18]. Therefore, the scheduling function should be computed to sufficient accuracy on classical computers using numerical quadrature, and implemented accurately as a control protocol on quantum computers.

We briefly discuss how to implement AQC(p) and AQC(exp) on a gate-based quantum computer. Since \( \mathinner {|{\psi _T(s)}\rangle } = \mathcal {T}\exp (-{\mathrm{i}} T\int _0^sH(f(s^{\prime }))ds^{\prime })\mathinner {|{\psi _T(0)}\rangle } \) , where \( \mathcal {T} \) is the time-ordering operator, it is sufficient to implement an efficient time-dependent Hamiltonian simulation of \( H(f(s)) \) .

One straightforward approach is to use the Trotter splitting method. The lowest order approximation takes the formwhere \( h = s/M, s_m = mh \) . It is proved in [31] that the error of such an approximation is \( {\mathcal {O}}(\text{poly} \) \( (\log (N))T^2/M) \) , which indicates that to achieve an \( \epsilon \) -approximation, it suffices to choose \( M = {\mathcal {O}}(\text{poly}(\log (N))T^2/\epsilon) \) . On a quantum computer, the operations \( e^{-{\mathrm{i}} \tau H_0} \) and \( e^{-{\mathrm{i}} \tau H_1} \) require a time-independent Hamiltonian simulation process, which can be implemented via techniques such as LCU and QSP [4, 22]. For a d-sparse matrix A, according to [5], the query complexity is \( f=\widetilde{{\mathcal {O}}}(g) \) if \( f={{\mathcal {O}}}(g \text{poly}\log (g)) \) for a single step. Here, the \( \widetilde{{\mathcal {O}}} \) means that we neglect the \( \log \log \) factors. Note that the total sum of the simulation time of single steps is exactly T regardless of the choice of M, and the total query complexity is \( \widetilde{{\mathcal {O}}}(dT\log (dT/\epsilon)) \) . Using Theorems 1 and 2, the query complexity of AQC(p) and AQC(exp) is \( \widetilde{{\mathcal {O}}}(d\kappa /\epsilon \log (d\kappa /\epsilon)) \) and \( \widetilde{{\mathcal {O}}}(d\kappa ~\text{poly}(\log (d\kappa /\epsilon))) \) , respectively. Nevertheless, M scales as \( {\mathcal {O}}(T^2) \) with respect to the runtime T, which implies that the number of time slices should be at least \( {\mathcal {O}}(\kappa ^2) \) . Therefore the gate complexity scales superlinearly with respect to \( \kappa \) . The scaling of the Trotter expansion can be improved using high order Trotter–Suzuki formula as well as the recently developed commutator-based error analysis [13], but we will not pursue this direction here.

There is an efficient way to directly perform the time evolution of \( H(f(s)) \) without using the splitting strategy, following the algorithm proposed by Low and Wiebe in [23], where the time-dependent Hamiltonian simulation is performed based on a truncated Dyson expansion. A detailed discussion on how to implement this algorithm in a gate-efficient way is presented in [Appendix C][20], and here we summarize the basic idea as follows. Assume that we are given two input query models: \( \mathcal {P}_A \) that gives the locations and values of the nonzero entries of the matrix A, and \( \mathcal {P}_B \) that produces the quantum state \( \mathinner {|{b}\rangle } \) . Here, the input query models are the same as those in [12]. Then, one can construct the block-encoding representations of the matrix A [16] and the matrix \( Q_b \) with \( {\mathcal {O}}(1) \) additional primitive gates. Next, the block-encodings of A and \( Q_b \) can be applied to build the block-encodings of \( H_0 \) and \( H_1 \) , and then the HAM-T model, which is a block-encoding of the select oracle of the time-dependent Hamiltonian \( H(s) \) evaluated at different time steps and serves as the input model in the truncated Dyson series method [23]. Finally, after the construction of HAM-T, the adiabatic dynamics can be simulated following the procedure for solving time-dependent Schrödinger equations discussed in [23].

The costs of AQC(p) and AQC(exp) are summarized in Table 1, where for both AQC(p) and AQC(exp), almost linear dependence with respect to \( \kappa \) is achieved. The almost linear dependence on \( \kappa \) cannot be expected to be improved to \( {\mathcal {O}}(\kappa ^{1-\delta }) \) with any \( \delta \gt 0 \) [17]. Thus, both AQC(p) and AQC(exp) are almost optimal with respect to \( \kappa \) , and AQC(exp) further achieves an exponential speedup with respect to \( \epsilon \) .

The QAOA scheme [14] considers the following parameterized wavefunctionHere, \( \theta \) denotes the set of \( 2P \) adjustable real parameters \( \lbrace \beta _i,\gamma _i\rbrace _{i=1}^{P} \) , and \( \mathinner {|{\psi _i}\rangle } \) is an initial wavefunction. The goal of QAOA is to choose \( \mathinner {|{\psi _i}\rangle } \) and to tune \( \theta \) , so that \( \mathinner {|{\psi _\theta }\rangle } \) approximates a target state. In the context of QLSP, we may choose \( \mathinner {|{\psi _i}\rangle }=\mathinner {|{{\widetilde{b}}}\rangle } \) , and each step of the QAOA ansatz in Equation (11) can be efficiently implemented using the quantum singular value transformation [16]. More specifically, as discussed in Section 5 and in [20], the block-encodings of \( H_0 \) and \( H_1 \) can be efficiently constructed via the input models for the matrix A and the vector \( \mathinner {|{b}\rangle } \) . Then, the quantum singular value transformation can be directly applied to simulate \( e^{-{\mathrm{i}} \beta H_0} \) and \( e^{-{\mathrm{i}} \gamma H_1} \) . According to [Corollary 62][16], the cost of each single simulation scales linearly in time and logarithmically in precision, and hence the total complexity of implementing a QAOA ansatz scales linearly in total runtime of QAOA, defined to be \( T:=\sum _{i=1}^P(|\beta _i|+|\gamma _i|) \) , and logarithmically in precision. Notice that with a sufficiently large P, the optimal Trotter splitting method becomes a special form of Equation (11). Hence, Theorem 2 implies that with an optimal choice of \( \lbrace \beta _i,\gamma _i\rbrace _{i=1}^P \) , the QAOA runtime T is at most \( {\mathcal {O}}(\kappa ~\text{poly}(\log (\kappa /\epsilon))) \) . We remark that the validity of such an upper bound requires a sufficiently large P and an optimal choice of \( \theta \) . On the other hand, our numerical results suggest that the same scaling can be achieved with a much smaller P.

For a given P, the optimal \( \theta \) maximizes the fidelity asHowever, the maximization of the fidelity requires the knowledge of the exact solution \( \mathinner {|{{\widetilde{x}}}\rangle }, \) which is not practical. We may instead solve the following minimization problemSince the null space of \( H_1 \) is of dimension 2, the unconstrained minimizer \( \mathinner {|{\psi _{\theta }}\rangle } \) seems possible to only have a small overlap with \( \mathinner {|{{\widetilde{x}}}\rangle } \) . However, this is not a problem due to the choice of the initial state \( \mathinner {|{\psi _i}\rangle }=\mathinner {|{{\widetilde{b}}}\rangle } \) . Notice that by the variational principle the minimizer \( \mathinner {|{\psi _\theta }\rangle } \) maximizes \( \mathinner {\langle {\psi _\theta |P_0(1)|\psi _\theta }\rangle } \) . Using the fact that \( e^{-{\mathrm{i}} \beta H_0}\mathinner {|{\bar{b}}\rangle }=e^{-{\mathrm{i}} \gamma H_1}\mathinner {|{\bar{b}}\rangle }=\mathinner {|{\bar{b}}\rangle } \) for any \( \beta ,\gamma \) , we obtain \( \mathinner {\langle {\bar{b}|\psi _\theta }\rangle }=\mathinner {\langle {\bar{b}|{\widetilde{b}}}\rangle }=0, \) which means the QAOA ansatz prevents the transition to \( \mathinner {|{\bar{b}}\rangle } \) , similar to AQC. Then, \( \mathinner {\langle {\psi _\theta |P_0(1)|\psi _\theta }\rangle }=\mathinner {\langle {\psi _\theta |{\widetilde{x}}}\rangle }\mathinner {\langle {{\widetilde{x}}|\psi _\theta }\rangle }=F_\theta \) , so the minimizer of Equation \( (12) \) indeed maximizes the fidelity.

For every choice of \( \theta \) , we evaluate the expectation value \( \mathinner {\langle {\psi _\theta |H_1^2|\psi _\theta }\rangle } \) . Then, the next \( \theta \) is adjusted on a classical computer towards minimizing the objective function. The process is repeated till convergence. Efficient classical algorithms for the optimization of parameters in QAOA are currently an active topic of research, including methods using gradient optimization [24, 36], Pontryagin’s maximum principle [3, 35], reinforcement learning [9, 26], to name a few. Algorithm 1 describes the procedure using QAOA to solve QLSP.

Compared to AQC(p) and AQC(exp), the QAOA method has the following two potential advantages. The first advantage is that QAOA provides the possibility of going beyond AQC-based algorithms. Notice that the Trotter splitting method is a special form of the QAOA ansatz in Equation (11). If the angles \( \lbrace \beta _i,\gamma _i\rbrace _{i=1}^P \) have been properly optimized (which is a very strong assumption and will be further discussed later), the total QAOA runtime T will be by definition comparable to or even shorter than the runtime of AQC with the best scheduling function (after discretization). Second, one way of implementing AQC(p) and AQC(exp) using an operator splitting method requires the time interval to be explicitly split into a large number of intervals, while numerical results indicate that the number of intervals P in QAOA can be much smaller. This could reduce the depth of the quantum circuit. Compared to AQC, QAOA has the additional advantage that it only consists of \( 2P \) time-independent Hamiltonian simulation problem, once \( \theta \) is known.

Despite the potential advantages, several severe caveats of using QAOA for QLSP arise when we consider beyond the time complexity. The first is that classical optimization of the angles \( \lbrace \beta _i,\gamma _i\rbrace _{i=1}^{P} \) can be difficult. Commonly used classical optimization algorithms, such as the gradient descent method, are likely to be stuck at local optimizers and thus result in sub-optimal performance. The cost for the classical optimization is also hard to known a priori. The optimization may require many iterations, which can diminish the gain of the runtime reduction. The second is related to the accurate computation of the objective function \( O(\theta ^{(k)}) \) . Note that the minimal spectrum gap of \( H_1 \) is \( {\mathcal {O}}(\kappa ^{-1}) \) . In order to obtain an \( \epsilon \) -approximation, the precision of measuring \( O(\theta)=\mathinner {\langle {\psi _{\theta }|H_1^2|\psi _\theta }\rangle } \) should be at least \( {\mathcal {O}}(\epsilon ^2/\kappa ^2) \) . Hence \( {\mathcal {O}}(\kappa ^4/\epsilon ^4) \) repeated measurements can be needed to achieve the desired accuracy.

Now we discuss the case when A is not Hermitian positive definite. First, we still assume that A is Hermitian (but not necessarily positive definite). In this case, we adopt the family of Hamiltonians introduced in [30], which overcomes the difficulty brought by the indefiniteness of A at the expense of enlarging the Hilbert space to dimension \( 4N \) (so two ancilla qubits are needed to enlarge the matrix block). Here, we definewhere \( Q_{+,b} = I_{2N}-\mathinner {|{+,b}\rangle }\mathinner {\langle {+,b}|} \) , and \( \mathinner {|{\pm }\rangle }=\frac{1}{\sqrt {2}}(\mathinner {|{0}\rangle }\pm \mathinner {|{1}\rangle }) \) . The null space of \( H_0 \) is \( \text{Null}(H_0)=\text{span}\lbrace \mathinner {|{0,-,b}\rangle }, \mathinner {|{1,+,b}\rangle }\rbrace \) . We also defineNote that \( \text{Null}(H_1)=\text{span}\lbrace \mathinner {|{0,+,x}\rangle },\mathinner {|{1,+,b}\rangle }\rbrace \) . Therefore, the solution of the QLSP can be obtained if we can prepare the zero-energy state \( \mathinner {|{0,+,x}\rangle } \) of \( H_1 \) .

The family of Hamiltonians for AQC(p) is still given by \( H(f(s)) = (1-f(s))H_0 + f(s)H_1, 0\le s\le 1 \) . Similar to the case of Hermitian positive definite matrices, there is a double degeneracy of the eigenvalue 0, and we aim at preparing one of the eigenstate via time-optimal adiabatic evolution. More precisely, for any s, \( \mathinner {|{1,+,b}\rangle } \) is always in \( \text{Null}(H(f(s))) \) , and there exists a state \( \mathinner {|{{\widetilde{x}}(s)}\rangle } \) with \( \mathinner {|{{\widetilde{x}}(0)}\rangle } = \mathinner {|{0,-,b}\rangle }, \mathinner {|{{\widetilde{x}}(1)}\rangle } = \mathinner {|{0,+,x}\rangle } \) , such that \( \text{Null}(H(f(s)))=\lbrace \mathinner {|{{\widetilde{x}}(s)}\rangle },\mathinner {|{1,+,b}\rangle }\rbrace \) . Such degeneracy will not influence the adiabatic computation starting with \( \mathinner {|{0,-,b}\rangle } \) for the same reason we discussed for Hermitian positive definite case (also discussed in [30]), and the error of AQC(p) is still bounded by \( \eta (s) \) given in Equation (4).

Furthermore, the eigenvalue 0 is separated from the rest of the eigenvalues of \( H(f(s)) \) by a gap \( \Delta (f(s))\ge \sqrt {(1-f(s))^2+(f(s)/\kappa)^2} \) [30]. For technical simplicity, note that \( \sqrt {(1-f)^2+(f/\kappa)^2} \ge (1-f+f/\kappa)/\sqrt {2} \) for all \( 0 \le f \le 1 \) , we define the lower bound of the gap to be \( \Delta _*(f) = (1-f+f/\kappa)/\sqrt {2} \) , which is exactly proportional to that for the Hermitian positive definite case. Therefore, we can use exactly the same time schedules as the Hermitian positive definite case to perform AQC(p) and AQC(exp) schemes, and properties of AQC(p) and AQC(exp) are stated in the following theorems (see Appendices D and E for the proof).

For a most general square matrix \( A\in {\mathbb {C}}^{N\times N} \) , following [17] we may transform it into the Hermitian case at the expense of further doubling the dimension of the Hilbert space. Consider an extended QLSP \( \mathfrak {A}\mathinner {|{\mathfrak {x}}\rangle } = \mathinner {|{\mathfrak {b}}\rangle } \) in dimension \( 2N \) whereNote that \( \mathfrak {A} \) is a Hermitian matrix of dimension \( 2N \) , with condition number \( \kappa \) and \( \Vert \mathfrak {A}\Vert _2 = 1 \) , and \( \mathinner {|{\mathfrak {x}}\rangle } := \mathinner {|{1,x}\rangle } \) solves the extended QLSP. Therefore, we can directly apply AQC(p) and AQC(exp) for Hermitian matrix \( \mathfrak {A} \) to prepare an \( \epsilon \) -approximation of x. The total dimension of the Hilbert space becomes \( 8N \) for non-Hermitian matrix A (therefore three ancilla qubits are needed).

We first report the performance of AQC(p), AQC(exp), and QAOA for a series of Hermitian positive definite dense matrices with varying condition numbers, together with the performance of RM and vanilla AQC. The details of the setup of the numerical experiments are given in Appendix F.

Figure 1 shows how the total runtime T depends on the condition number \( \kappa \) and the accuracy \( \epsilon \) for the Hermitian positive definite case. The numerical scaling is reported in Table 2. For the \( \kappa \) dependence, despite that RM and AQC(1) share the same asymptotic linear complexity with respect to \( \kappa \) , we observe that the preconstant of RM is larger due to its Monte Carlo strategy and the mixed state nature resulting in the same scaling of errors in fidelity and density (see Appendix C for a detailed explanation). The asymptotic scaling of the vanilla AQC is at least \( {\mathcal {O}}(\kappa ^2) \) . When higher fidelity (0.999) is desired, the cost of vanilla AQC becomes too expensive, and we only report the timing of AQC(p), AQC(exp), and QAOA. For the \( \kappa \) dependence tests, the depth of QAOA ranges from 8 to 60. For the \( \epsilon \) dependence test, the depth of QAOA is fixed to be 20. We find that the runtime for AQC(p), AQC(exp), and QAOA depends approximately linearly on \( \kappa \) , while QAOA has the smallest runtime overall. It is also interesting to observe that although the asymptotic scalings of AQC(1) and AQC(2) are both bounded by \( {\mathcal {O}}(\kappa \log \kappa) \) instead of \( {\mathcal {O}}(\kappa) \) , the numerical performance of AQC(2) is much better than AQC(1); in fact, the scaling is very close to that with the optimal value of p. For the \( \epsilon \) dependence, the scaling of RM and AQC(p) is \( {\mathcal {O}}(1/\epsilon) \) , which agrees with the error bound. Again the preconstant of RM is slightly larger. Our results also confirm that AQC(exp) only depends poly logarithmically on \( \epsilon \) . Note that when \( \epsilon \) is relatively large, AQC(exp) requires a longer runtime than that of AQC(p), and it eventually outperforms AQC(p) when \( \epsilon \) is small enough. The numerical scaling of QAOA with respect to \( \epsilon \) is found to be only \( {\mathcal {O}}(\log ^{1.5}(1/\epsilon)) \) together with the smallest preconstant.

Figure 2 and Table 3 demonstrate the simulation results for non-Hermitian matrices. We find that numerical performances of RM, AQC(p), AQC(exp), and QAOA are similar to that of the Hermitian positive definite case. Again QAOA obtains the optimal performance in terms of the runtime. The numerical scaling of the optimal AQC(p) is found to be \( {\mathcal {O}}(\kappa /\epsilon) \) , while the time complexity of QAOA and AQC(exp) is only \( {\mathcal {O}}(\kappa ~\text{poly}(\log (\kappa /\epsilon))) \) .

By reformulating QLSP into an eigenvalue problem, AQC provides an alternative route to solve QLSP other than those based on phase estimation (such as HHL) and those based on approximation of matrix functions (such as LCU and QSP). However, the scaling of the vanilla AQC is at least \( {\mathcal {O}}(\kappa ^2/\epsilon) \) , which is unfavorable with respect to both \( \kappa \) and \( \epsilon \) . Thanks to the explicit information of the energy gap along the adiabatic path, we demonstrate that we may reschedule the AQC and dramatically improve the performance of AQC for solving QLSP. When the target accuracy \( \epsilon \) is relatively large, the runtime complexity of the AQC(p) method ( \( 1\lt p\lt 2 \) ) is reduced to \( {\mathcal {O}}(\kappa /\epsilon) \) ; for highly accurate calculations with a small \( \epsilon \) , the AQC(exp) method is more advantageous, and its runtime complexity is \( {\mathcal {O}}(\kappa ~\text{poly}(\log (\kappa /\epsilon))) \) . To our knowledge, our ACP(exp) method provides the first example that an adiabatic algorithm can simultaneously achieve near linear dependence on the spectral gap, and poly-logarithmic dependence on the precision.

Due to the close connection between AQC and QAOA, the runtime complexity of QAOA for solving QLSP is also bounded by \( {\mathcal {O}}(\kappa ~\text{poly}(\log (\kappa /\epsilon))) \) . Both AQC and QAOA can be implemented on gate-based quantum computers. Our numerical results can be summarized using the following relation:Here, \( A\lt B \) means that the runtime of A is smaller than that of B. The two exceptions are: \( \text{QAOA}\lesssim \text{AQC(exp)} \) means that the runtime of QAOA is smaller only when the optimizer \( \theta \) is found, while \( \text{AQC(exp)}\lesssim \text{AQC}(p) \) holds only when \( \epsilon \) is sufficiently small. While the runtime complexity of AQC(exp) readily provides an upper bound of the runtime complexity of QAOA, numerical results indicate that the optimizer of QAOA often involves a much smaller depth, and hence the dynamics of QAOA does not necessarily follow the adiabatic path. Therefore, it is of interest to find alternative routes to directly prove the scaling of the QAOA runtime without relying on AQC. In the work [20], our AQC based algorithm has been combined with the eigenvector filtering technique. Reference [20] also proposed another AQC inspired quantum linear system solver, which is based on the quantum Zeno effect. Both methods can scale linearly in \( \kappa \) and logarithmically in \( 1/\epsilon \) . We expect our AQC based QLSP solvers may serve as useful subroutines in other quantum algorithms as well.

The Hamiltonian \( H(f) \) can be written in the block matrix form as

Let \( \lambda \) be an eigenvalue of H, thenwhere the second equality holds because the bottom two blocks are commutable. Thus, \( \lambda ^2 \) is an eigenvalue of \( ((1-f)I + fA)Q_b^2((1-f)I + fA) \) , and \( \Delta ^2(f) \) equals the smallest non-zero eigenvalue of \( ((1-f)I + fA)Q_b^2((1-f)I + fA) \) . Applying a proposition of matrices that XY and YX have the same non-zero eigenvalues, \( \Delta ^2(f) \) also equals the smallest non-zero eigenvalue of \( Q_b((1-f)I + fA)^2Q_b \) .

Now we focus on the matrix \( Q_b((1-f)I + fA)^2Q_b \) . Note that \( \left|b\right\gt \) is the unique eigenstate corresponding to the eigenvalue 0, and all eigenstates corresponding to non-zero eigenvalues must be orthogonal to with \( \left|b\right\gt \) . Therefore,and \( {\Delta }(f) \ge {\Delta }_*(f) = 1-f+f/\kappa \) .

We show two relations that connect the error of density matrix distance and the error of fidelity, which are used in our proof for AQC(p) and AQC(exp). Following the notations in the main text, let \( \mathinner {|{{\widetilde{x}}(s)}\rangle } \) denote the desired eigenpath of \( H(f(s)) \) corresponding to the 0 eigenvalue, and \( \text{Null}(H(f(s)))=\lbrace \mathinner {|{{\widetilde{x}}(s)}\rangle },\mathinner {|{\bar{b}}\rangle }\rbrace \) . \( P_0(s) \) denotes the projector onto the eigenspace corresponding to the 0 eigenvalue.

In our numerical test, we observe that to reach a desired fidelity, RM encounters a much larger pre-constant than AQC(p). This is due to the following reason. Although the runtime of both RM and AQC(p) scales as \( {\mathcal {O}}(1/\epsilon) \) where \( \epsilon \) is the 2-norm error of the density matrix, the scalings with respect to the infidelity are different. More specifically, Lemma 5 shows that for AQC, the square of the 2-norm error is exactly equal to the infidelity. Thus, in order to reach infidelity \( 1-F \) using AQC(p), it suffices to choose the runtime to be \( T = {\mathcal {O}}(\kappa /\sqrt {1-F}) \) . Meanwhile, it has been proved in [6] that the runtime complexity of RM is \( {\widetilde{\mathcal {O}}}(\kappa /(1-F)) \) . Therefore, to reach a desired fidelity, the runtime of AQC(p) will be smaller than that of RM, as shown in our numerical examples.

We further verify the statement above by studying the relation between the 2-norm error of the density matrix and the infidelity for AQC(p), AQC(exp), and RM using the positive definite example with \( \kappa = 10 \) . In AQC(p) and AQC(exp), we change the runtime to obtain approximations with different errors and infidelity. In RM we vary the number of exponential operators to obtain different levels of accuracy. All other numerical treatments remain unchanged. As shown in Figure 3, the infidelity is exactly the square of 2-norm error in the case of AQC(p) and AQC(exp), while the infidelity of RM only scales approximately linearly with respect to 2-norm error. This also verifies that although the runtime of both AQC(p) and RM scales linearly with respect to \( \epsilon \) , the runtime of AQC(p) can be much smaller to reach desired fidelity.

The proof of Theorems 1 and 3 rests on some delicate cancellation of the time derivatives \( {\parallel }{H^{(1)}}_2, \) \( {\parallel }{H^{(2)}}_2 \) and the gap \( \Delta (f(s)) \) in the error bound, and can be completed by carefully analyzing the \( \kappa \) -dependence of each term in \( \eta (s) \) given in Equation (4). Note that in both cases \( {H}(f) = (1-f)H_0+fH_1 \) , and we let \( \Delta _*(f) = (1-f+f/\kappa)/\sqrt {2} \) since such a choice of \( \Delta _* \) can serve as a lower bound of the spectrum gap for both the cases of Theorems 1 and 3. We first compute the derivatives of \( H(f(s)) \) by chain rule asandThen the first two terms of \( \eta (s) \) can be rewritten asHere, C stands for a general positive constant independent of \( s,\Delta , \) and T. To compute the remaining two terms of \( \eta (s) \) , we use the following change of variableand the last two terms of \( \eta (s) \) becomeand similarlySummarize all terms above, an upper bound of \( \eta (s) \) isFinally, since for \( 1\lt p\lt 2 \) andwe haveThe leading term of the bound is \( {\mathcal {O}}(\kappa /T) \) when \( 1\lt p\lt 2 \) .

Now we consider the limiting case when \( p = 1,2 \) . Note that the bound for \( \eta (s) \) can still be written asStraightforward computation shows thatandHence, when \( p=1,2 \) ,This completes the proof of Theorems 1 and 3.

We provide a rigorous proof of the error bound for the AQC(exp) scheme. We mainly follow the methodology of [25] and a part of technical treatments of [15]. Our main contribution is carefully revealing an explicit constant dependence in the adiabatic theorem, which is the key to obtain the \( \widetilde{{\mathcal {O}}}(\kappa) \) scaling. In the AQC(exp) scheme, the Hamiltonian \( H(s) = (1-f(s))H_0 + f(s)H_1 \) with \( \Vert H_0\Vert , \Vert H_1\Vert \le 1 \) andThe normalization constant \( c_e = \int _0^1 \exp (-\frac{1}{t(1-t)})dt \approx 0.0070 \) . Let \( U_T(s) \) denote the corresponding unitary evolution operator, and \( P_0(s) \) denote the projector onto the eigenspace corresponding to 0. We use \( \Delta _*(f) = (1-f+f/\kappa)/\sqrt {2} \) since this can serve as a lower bound of the spectrum gap for both the cases of Theorems 2 and 4.

We first restate the theorems universally with more technical details as following.

The remaining part is devoted to some preliminary results regarding \( H, E, \) and the technical estimates for the growth of \( E_j \) . It is worth mentioning in advance that in the proof we will encounter many derivatives taken on a contour integral. In fact all such derivatives taken on a contour integral will not involve derivatives on the contour. Specifically, since \( (H(s)-z)^{-1} \) is analytic for any \( 0 \lt |z| \lt \Delta (s) \) , for any \( s_0 \in (0,1) \) , there exists a small enough neighborhood \( B_{\delta }(s_0) \) such that \( \forall s \in B_{\delta }(s_0) \) , \( \oint _{\Gamma (s)} G(s,(H(s)-z)^{-1})dz = \oint _{\Gamma (s_0)} G(s,(H(s)-z)^{-1})dz \) for any smooth mapping G. This means locally the contour integral does not depend on the smooth change of the contour, and thus the derivatives will not involve derivatives on the contour. In the spirit of this trick, we write the resolvent \( R(z,s,s_0) = (H(s)-z)^{-1} \) for \( 0 \le s \le 1, 0 \le s_0 \le 1, z \in {\mathbb {C}} \) and \( |z| = \Delta (s_0)/2 \) and let \( R^{(k)} \) denote the partial derivative with respect to s, i.e., \( \frac{\partial }{\partial s}R(z,s,s_0) \) , which means by writing \( R^{(k)} \) we only consider the explicit time derivatives brought by H.

The following lemma gives the bound for the derivatives of H.

The following result demonstrates that \( E_{j} \) ’s for all \( j \ge 1 \) vanish on the boundary.

Before we process, we recall three technical lemmas introduced in [15, 25]. Throughout let \( c_f = 4\pi ^2/3 \) denote an absolute constant.

Next we bound the derivatives of the resolvent. This bound provides the most important improvement of the general adiabatic bound.

The next lemma is the main technical result, which gives the bound of derivatives of \( E_{j} \) defined in Equation (20).

Proof. (a) By Lemma 14,

(b) We prove by induction with respect to j. For \( j = 1 \) , Equation (22) tellsBy Lemmas 13 and 14,Now assume \( \lt j \) the estimate holds, for j, by Lemmas 12, 13, and the induction hypothesis,Again by Lemmas 13, 14, and the induction hypothesis,

For simulation purpose, the AQC schemes are carried out using and the induction hypothesis with a time step size 0.2. We use the gradient descent method to optimize QAOA and record the running time corresponding to the lowest error in each case. In QAOA, we also use the true fidelity to measure the error. RM is a Monte Carlo method, and each RM calculation involves performing 200 independent runs to obtain the density matrix \( \rho ^{(i)} \) for i-th repetition, then we use the averaged density \( \bar{\rho } = 1/n_{\text{rep}}\sum \rho ^{(i)} \) to compute the error. We report the averaged runtime of each single RM calculation. We perform calculations for a series of 64-dimensional Hermitian positive definite dense matrices \( A_1 \) , and 32-dimensional non-Hermitian dense matrices \( A_2 \) with varying condition number \( \kappa \) .

For concreteness, for the Hermitian positive definite example, we choose \( A = U\Lambda U^{\dagger } \) . Here, U is an orthogonal matrix obtained by Gram–Schmidt orthogonalization (implemented via a QR factorization) of the discretized periodic Laplacian operator given by \( \Lambda \) is chosen to be a diagonal matrix with diagonals uniformly distributed in \( [1/\kappa ,1] \) . More precisely, \( \Lambda = \text{diag}(\lambda _1,\lambda _2,\ldots ,\lambda _{N}) \) with \( \lambda _k = 1/\kappa + (k-1)h, h = (1-1/\kappa)/(N-1) \) . Such construction ensures A to be a Hermitian positive definite matrix which satisfies \( \Vert A\Vert _2 = 1 \) and the condition number of A is \( \kappa \) . We choose \( \mathinner {|{b}\rangle } = \sum _{k=1}^{N}u_k / \Vert \sum _{k=1}^{N}u_k\Vert _2 \) where \( \lbrace u_k\rbrace \) is the set of the column vectors of U. Here, \( N=64 \) .

For the non-Hermitian positive definite example, we choose \( A = U\Lambda V^{\dagger } \) . Here, U is the same as those in the Hermitian positive definite case, except that the dimension is reduced to \( N=32 \) . \( \Lambda = \text{diag}(\lambda _1,\lambda _2,\ldots ,\lambda _N) \) with \( \lambda _k = (-1)^k(1/\kappa +(k-1)h), h = (1-1/\kappa)/(N-1) \) . V is an orthogonal matrix obtained by Gram–Schmidt orthogonalization of the matrixSuch construction ensures A to be non-Hermitian, satisfying \( \Vert A\Vert _2 = 1 \) and the condition number of A is \( \kappa \) . We choose the same \( \mathinner {|{b}\rangle } \) as that in the Hermitian positive definite example.

We thank Rolando Somma, Yu Tong and Nathan Wiebe for helpful discussions.