HHL: solving linear systems, and what the condition number does to the speedup

Definition

The HHL algorithm, due to Harrow, Hassidim, and Lloyd, estimates features of the solution of a linear system. Given an N by N matrix A and a vector b, the equation is A x = b. The 2009 paper treats the case in which the desired output is not the list of entries of x. It is an approximation of an expectation x† M x for some operator M, or a quantum state that encodes x.

The algorithm assumes A is Hermitian. If it is not, the paper embeds it in a larger Hermitian matrix whose off-diagonal blocks are A and its adjoint, and solves that system for a right-hand side whose second block is zero. The solution of the larger system then carries x in its second block. The rest of the argument is written for Hermitian A.

Two access assumptions come before any runtime. The matrix is s-sparse and efficiently row computable: at most s nonzero entries in a row, and, given a row index, those entries can be computed in time proportional to s. The vector b is presented as a quantum state, a superposition whose amplitudes are the entries of a unit vector. The paper gives two ways that state can arise. The entries and the partial sums of their squared magnitudes are efficiently computable, in which case a state-preparation procedure of Grover and Rudolph applies, or some other part of a larger quantum algorithm produces the state. The number of gates in that preparation is called T_B in the appendix. It is not free.

The output, after a successful run, is a quantum state proportional to A inverse applied to the input state, not a classical vector. The introduction is explicit that reading every component would take at least N runs. What one measurement can return is a feature of that state: an expectation, a normalization, a weight on a subset of basis states.

Quantum phase estimation is the subroutine that reads the eigenvalues. Hamiltonian simulation is the subroutine that supplies the unitary e to the power i A t. The cost used below is the one Harrow, Hassidim, and Lloyd write down for the simulation paper they cite. The simulation page is where the later product-formula and qubitization bounds are stated.

Intuition

Simulate the matrix as a Hamiltonian, read its eigenvalues, and write their reciprocals onto an ancilla. Phase estimation, given a superposition of evolution times, decomposes the input state in the eigenbasis and stores an estimate of each eigenvalue. A controlled rotation then maps an eigenvalue estimate λ to an amplitude proportional to 1/λ on an ancilla, with a normalizing constant on the order of 1/κ. The eigenvalue register is uncomputed. Measuring the ancilla and keeping the outcome that marks success leaves a state proportional to the inverted eigenvalues applied to the input.

The rotation is not unitary on its own, so the postselection can fail. If every singular value is at most 1 and the constant in the rotation is on the order of 1/κ, the success probability is at least on the order of 1/κ². Amplitude amplification turns that into on the order of κ repetitions rather than on the order of κ². The appendix makes the schedule explicit: a geometrically increasing number of repetitions, up to a power of two that is at least κ, succeeds with constant probability using at most 4κ repetitions.

The same appendix does not invert an ill-conditioned direction blindly. It chooses a cutoff κ and a looser cutoff twice as large. Eigenvalues at least 1/κ are inverted. Eigenvalues at most 1 over 2κ are flagged as ill-conditioned and left uninverted. Eigenvalues in between are interpolated by a pair of filter functions. The flag is part of the output, so a user can estimate how much of the input sat in the bad subspace. The guarantee ||A|| ≤ 1 and ||A inverse|| ≤ κ is what makes the cutoff honest. If the smallest singular value is zero, there is no inverse on the whole space, and the flag is the paper's way of saying so.

Complexity

This is a gate-complexity bound under the sparse-access and state-preparation assumptions. It is not a query-complexity bound in the sense of Grover, and it is not a promise that a classical list of N solution entries can be printed in polylogarithmic time.

The condition number κ is the ratio of the largest eigenvalue magnitude to the smallest. The introduction normalizes so that the singular values lie between 1/κ and 1, which is the same as κ to the minus 2 times the identity sitting below A† A, which sits below the identity. In that normalization the introduction says the runtime scales as κ² log(N) / ε, where ε is the additive error in the output state. The largest advantage over classical algorithms, the paper says, is when both κ and 1/ε are polylogarithmic in N.

The section headed "Run-time and error analysis" puts the sparsity back in. Simulating e to the power i A t for an s-sparse matrix, to error ε, costs Õ of t s², with the tilde hiding factors that grow more slowly than any positive power of 1/ε. Phase estimation is the dominant error. An evolution time t₀ on the order of κ/ε makes the error in the inverted eigenvalue at most ε when every eigenvalue is at least 1/κ. Together with the κ repetitions from amplification, that section states the runtime Õ of log(N) s² κ² / ε.

Appendix A.1 is the version that keeps the state-preparation cost visible. The run time is Õ of κ T_B plus κ² s² log(N) / ε. The tilde suppresses a squared iterated logarithm of N, a sub-polynomial factor from the simulation error, and a polylogarithmic cost for preparing the phase-estimation register. Theorem 1 is the error bound those choices of t₀ are aimed at. The paper's distance between two pure states is the square root of 2 times 1 minus the real part of their inner product. With no postselection, the implemented inversion is at most O(κ / t₀) from the ideal in operator distance. With postselection onto the flag register, or onto the "well" flag when the input already lies in the well-conditioned subspace, the output state is at most O(κ / t₀) from the ideal normalized state. Taking t₀ on the order of κ/ε makes that distance at most ε.

The simulation cost inside those expressions is not re-derived on this page. Section II states that, for s-sparse efficiently row-computable A, the algorithm of Berry, Ahokas, Cleve, and Sanders simulates e to the power i A t in time Õ of log(N) s² t. That quotation is the one to keep next to the 2009 runtime. Hamiltonian simulation is the later comparison of product formulas, the truncated Taylor series, and qubitization. The thing to cite for the inversion runtime is still HHL.

Classical conjugate gradient, as HHL compare it, uses on the order of the square root of κ times log(1/ε) matrix-vector multiplications when A is positive definite, each multiplication costing on the order of N s, for a total on the order of N s times the square root of κ times log(1/ε). Without positive definiteness the multiplication count becomes on the order of κ log(1/ε). The dependence on 1/ε is logarithmic. HHL's is not. The paper treats that gap as real, and then argues it cannot be removed for the estimation problem.

The estimation problem is defined just before Theorem 6. The output is a number within ε of the expectation of a fixed observable on the normalized solution state, with probability at least 2/3. Producing the state to trace distance ε/2, in time Õ of log(N) κ² s² / ε, and then taking on the order of 1/ε² samples of the observable, solves the estimation problem in time Õ of log(N) κ² s² / ε³. Theorem 6 says a quantum algorithm for that estimation problem running in time polynomial in κ, log(N), and log(1/ε) would imply BQP equals PP. The polylogarithmic precision that later papers achieve is for the state, not for this classical number.

Theorem 4 is the optimality claim for the condition number and for classical algorithms. A quantum matrix-inversion algorithm running in time κ to the power 1 minus δ, times polylog(N), for some positive δ, would imply BQP equals PSPACE. No relativizing quantum algorithm can run in that time either. A classical algorithm running in time polynomial in κ and log(N) would imply BPP equals BQP. The reduction behind those claims, in the body of the letter, simulates an n-qubit, T-gate circuit by inverting an O(1)-sparse matrix of dimension on the order of 2 to the n times κ. The condition number is on the order of T if the matrix need not be positive definite, and on the order of T squared if it must. Measuring a clock register in a window that the paper shows has probability at least 1/10 leaves the data register in the state of the simulated computation. The letter calls the resulting problem BQP-complete. The speedup over classical algorithms is, on those instances, the speedup of the encoded quantum computation.

Childs, Kothari, and Somma study a different accounting, and a better dependence on ε for the state. Their Problem 1, the quantum linear systems problem, asks for a state whose vector distance from the normalized solution is at most ε, with success probability at least 1/2 and a flag. The matrix is Hermitian, d-sparse, of norm 1, with known condition number κ. An oracle returns the positions and values of the nonzero entries. A second procedure prepares the input state in time polynomial in log(N). They say explicitly that this is not the traditional problem of solving a linear system. Theorem 5 gives a gate-efficient algorithm making on the order of d κ times a polynomial in log(d κ / ε) queries to both oracles. Gate-efficient, in their definition, means the gate count exceeds the query count by at most a polynomial in the log of the query count and of N. They record that Ambainis had already improved the condition-number dependence of HHL from quadratic to nearly linear, by a method that uses phase estimation, and that copying that method unchanged would bring back a polynomial dependence on 1/ε. Theorem 5 is their bound with both improvements. Theorems 3 and 4 of the same paper are the earlier, quadratic-in-κ algorithms, with the polylogarithmic ε dependence written out: the Chebyshev version uses on the order of d κ² log squared of (d κ / ε) entry queries. The exponential improvement in ε is for preparing the state. They cite HHL's Theorem 6 for the statement that sampling error alone rules out a polylogarithmic dependence on ε for expectation values, unless BQP equals PP.

Applications and limits

Aaronson's 2015 commentary lists four caveats, and each one is already a hypothesis of the theorems above. The vector b has to be loaded into a quantum state. If that loading takes a polynomial in N, the advantage in N is gone at the first step. He notes that a right-hand side with a single 1 and the rest zeros would contradict the Bennett, Bernstein, Brassard, and Vazirani search lower bound if the algorithm still ran in time logarithmic in N. The machine also has to apply e to the power minus i A t. Sparsity, together with a way to locate the nonzero entries, is one sufficient condition he names; if the simulation itself takes polynomial time in N, the advantage is gone. The matrix has to be well conditioned. If κ grows as a power of N, the advantage is gone. His commentary describes the cost as growing nearly linearly with κ. That sentence is his summary, not Theorem 5 and not the 2009 runtime. The proved dependence is quadratic in the 2009 appendix and nearly linear in the later query bound. Finally, the output is a quantum state on the order of log N qubits. A measurement reveals statistical information. Learning one specific entry, he writes, generally means repeating the algorithm on the order of N times.

He also says the algorithm is best read as a template. An application has to supply the state preparation, the simulation, and a measurement that is actually the quantity one wanted, and then compare that restricted problem with the best classical algorithm for the same restricted problem. He discusses one electromagnetic-scattering proposal as an attempt to do that from start to finish, and he reports that the proposal could not rule out a classical algorithm once the restrictions were in place. This page does not re-derive that proposal. The point he draws from HHL's own reduction is the one Theorem 4 already supports: there are linear systems, built by encoding a quantum circuit, on which a fast classical inversion would collapse BQP to BPP, and those systems are as artificial as the circuit they encode.

The dequantization result is a classical algorithm for a different input. Gilyén, Lloyd, and Tang assume length-square sampling access to a low-rank matrix: a data structure from which one can sample a row with probability proportional to its squared Euclidean norm, and query entries, in time polylogarithmic in the dimensions once the structure is built. Theorem 1 says that if the rank is at most k, the operator norm is at most 1, the norm of the pseudoinverse is at most κ, and the projection of b onto the column space has norm on the order of the norm of b, then their algorithm returns a vector within multiplicative error ε of the pseudoinverse applied to b, with probability at least 1 minus a failure parameter. Theorem 2 puts the cost at Õ of κ to the 16, times k to the 6, times the sixth power of the Frobenius norm, divided by ε to the 6. The tilde hides polylogarithmic factors in that cost, in the dimensions, and in the failure probability. The output is an implicit description that supports queries and samples, not a quantum state. A footnote says they are treating matrices whose singular values are not too small, and leaving a well-conditioned-subspace version to later work. None of this is a classical solver for the sparse BQP-complete problem in HHL's Theorem 4. The access model is sampling from a stored low-rank matrix. The rank appears to the sixth power. Logarithmic dependence on the dimension is real inside that model, and it is not a claim that every sparse linear system has a classical polylogarithmic algorithm.

Linear systems is the stub that keeps this comparison next to the phase-estimation page. Phase estimation's own hypotheses, an eigenstate or a decent overlap and a unitary one can actually power, are unchanged by being used as a subroutine.

Reading the source

Read section II of Harrow, Hassidim, and Lloyd for the circuit and the sparsity assumption, the run-time paragraph there for the Õ of log(N) s² κ² / ε expression, and Appendix A.1 together with Theorem 1 for the runtime that includes T_B and for the vector distance. Read Theorems 4 and 6 before improving the dependence on κ or on ε in a quotation. Read Childs, Kothari, and Somma, Problem 1 and Theorem 5, for the state-preparation bound with a logarithmic dependence on the error, and their citation of HHL's Theorem 6 before exporting that logarithm to an expectation value. Read Aaronson for the four caveats. Read Gilyén, Lloyd, and Tang, Theorems 1 and 2, before describing dequantization as a classical HHL. The phase-estimation page is the account of the eigenvalue subroutine and of the eigenstate it assumes.