HHL Linear Systems Algorithm

The HHL algorithm is the moment quantum computing stopped being about counting oracle queries and started making a claim on the everyday work of numerical science: solving $Ax=b$ . Almost every simulation, every regression, every discretized differential equation bottoms out in a linear system, and classically the cost of solving one grows at least linearly with the number of unknowns $N$ — touch every entry of the answer and you have already paid $\Omega(N)$ . Harrow, Hassidim, and Lloyd's 2009 construction refuses to pay that. It never writes the solution down. Instead it encodes the entire $N$ -component vector $x$ into the amplitudes of an $n=\log_2 N$ qubit state $|x\rangle$ , and it manufactures that state in time polylogarithmic in $N$ . A million-dimensional system lives in twenty qubits; a billion-dimensional one in thirty. The engine is spectral. Written in its own eigenbasis, $A=\sum_j\lambda_j|u_j\rangle\langle u_j|$ , and solving $Ax=b$ means nothing more than dividing each component of $b$ by its eigenvalue $\lambda_j$ . A classical machine does those divisions one at a time; HHL does all $N$ at once. Quantum phase estimation on the unitary $e^{iAt_0}$ tags every eigenvector $|u_j\rangle$ with its eigenvalue $\lambda_j$ in an ancillary register, a single controlled rotation writes an amplitude proportional to $1/\lambda_j$ , and running phase estimation backwards erases the scratch work — leaving a state proportional to $\sum_j\beta_j\lambda_j^{-1}|u_j\rangle=A^{-1}|b\rangle$ . The whole matrix inversion happens in superposition, in one pass. Three caveats keep the claim honest, and understanding them is understanding HHL. The output is a quantum state, so reading all $N$ amplitudes by tomography costs $\Omega(N)$ and erases the speedup — the win is real only for a summary property $\langle x|M|x\rangle$ . The runtime carries the condition number $\kappa$ and the sparsity $s$ , so ill-conditioned or dense systems bleed the advantage away. And the whole edifice assumes coherent oracle access to $A$ and efficient preparation of $|b\rangle$ . Far from weakening the result, these conditions are what make it precise: HHL draws the exact boundary of where quantum computers deliver exponential advantage in linear algebra.

Key takeaways

HHL prepares a quantum state $|x\rangle$ whose amplitudes are the (normalized) entries of the solution to $Ax=b$ , for an $N\times N$ Hermitian, $s$ -sparse, well-conditioned matrix $A$ .
The original runtime is $\tilde O(s^2\kappa^2\log N/\epsilon)$ — exponentially faster in the dimension $N$ than the $O(Ns\sqrt{\kappa}\log(1/\epsilon))$ of classical conjugate gradient.
The mechanism is quantum phase estimation on $e^{iAt_0}$ , which loads each eigenvalue $\lambda_j$ into an ancilla, followed by a controlled rotation that writes an amplitude $\propto 1/\lambda_j$ , then uncomputation.
The speedup is genuine only for extracting a global property $\langle x|M|x\rangle$ ; full read-out by tomography costs $\Omega(N)$ and destroys the advantage.
The condition number $\kappa$ enters because the inversion succeeds with probability $\Theta(1/\kappa^2)$ ; amplitude amplification softens $O(\kappa^2)$ repetitions to $O(\kappa)$ .
Later work (Ambainis; Childs–Kothari–Somma) cut the dependence to $\kappa\,\mathrm{polylog}(1/\epsilon)$ , and HHL is now the seed of the whole field of quantum linear algebra.
Key idea: phase estimation on $e^{iAt_0}$ transfers each eigenvalue $\lambda_j$ of $A$ into an ancilla register, a controlled rotation writes an amplitude $\propto 1/\lambda_j$ , and uncomputing the register leaves a state proportional to $A^{-1}|b\rangle=\sum_j\beta_j\lambda_j^{-1}|u_j\rangle$ — the entire spectral inversion carried out in superposition.

The core insight

Phase estimation on $e^{iAt_0}$ transfers each eigenvalue $\lambda_j$ of $A$ into an ancilla register; a controlled rotation then writes an amplitude $\propto 1/\lambda_j$ , so that uncomputing the eigenvalue register leaves a state proportional to $A^{-1}|b\rangle=\sum_j\beta_j\lambda_j^{-1}|u_j\rangle$ . The entire spectral inversion is done in superposition.

Where this lives today

HHL is small enough to demonstrate and deep enough to seed a field. Its 2×2 and larger instances have been run on photonic processors (Cai and collaborators, 2013), on NMR, and on superconducting devices, each showing the full choreography — state preparation, phase estimation, controlled rotation, uncomputation — in one compact circuit, precisely because a working run is a proof of concept for coherent linear algebra. But the living importance is in what HHL became. It is the archetype of the quantum linear systems problem, now understood through block-encoding and the quantum singular value transformation, where solving $Ax=b$ , simulating Hamiltonians, and estimating spectra are all read off from one unifying framework that traces directly back to the HHL rotation $\lambda\mapsto 1/\lambda$ . HHL-style subroutines drive quantum solvers for linear differential equations (Berry; Childs–Liu–Ostrander), quantum algorithms for least squares and finite-element analysis, and proposed speedups in machine learning — support-vector machines, Gaussian processes, recommendation systems. That last frontier is where the caveats bite hardest: Ewin Tang's dequantization program produced classical sampling algorithms matching several HHL-based proposals, turning the "fine print" into an active research boundary that continually re-tests where quantum advantage in linear algebra genuinely lives. The open problem — for which practically important, well-conditioned, sparse systems the input can be loaded and the output extracted without erasing the speedup — is as central now as it was in 2009.

The physical picture

Picture the matrix $A$ in its own eigenbasis, where it is simply a list of multipliers: the eigenvector $|u_j\rangle$ gets scaled by $\lambda_j$ . Solving $Ax=b$ is then the most elementary act imaginable — take the input $|b\rangle=\sum_j\beta_j|u_j\rangle$ , and divide each spectral component by its eigenvalue to get $\sum_j\beta_j\lambda_j^{-1}|u_j\rangle$ . The classical difficulty is that the eigenbasis is hidden: you would have to diagonalize $A$ first, an $O(N^3)$ chore. HHL's sleight of hand is that it never diagonalizes anything. Quantum phase estimation, run on the unitary $e^{iAt_0}$ that Hamiltonian simulation builds from $A$ , acts diagonally on each $|u_j\rangle$ automatically — it writes the eigenvalue $\lambda_j$ into a clock register without anyone ever computing an eigenvector. Every eigen-branch is processed in parallel because they all ride the same superposition.

Figure

The HHL circuit worked on a 2×2 example: phase estimation reads each eigenvalue into the clock register, a controlled rotation writes amplitude C/λ onto the flag, phase estimation is undone, and post-selecting the flag on 1 leaves the solution in the input register. For A with eigenvalues λ = 1 and 2 and C = 1 the rotation angles are θ = 180° and 60°, and the flag reads 1 with probability 62.5%

The controlled rotation is the algorithm's one genuinely clever move. Conditioned on the clock register holding $\tilde\lambda_j$ , a single ancilla qubit is rotated by $\theta_j=2\arcsin(C/\tilde\lambda_j)$ , depositing amplitude exactly $C/\tilde\lambda_j$ on its $|1\rangle$ branch. A large eigenvalue gives a small amplitude, a small eigenvalue a large one — precisely the $\lambda\mapsto 1/\lambda$ scaling that inverts $A$ . Post-selecting that flag ancilla on $|1\rangle$ keeps only the inverted branch, but the probability of landing there is $\sim(C/\lambda)^2$ , and the smallest eigenvalue forces $C\lesssim|\lambda_{\min}|\sim 1/\kappa$ . That single constraint is the taproot of every $\kappa$ in the cost. The payoff is the exponential compression of dimension into qubits — and its flip side, the reason you can never cheaply read the answer back out, is drawn in the second figure.

Figure

Runtime versus problem dimension N = 2ⁿ on a logarithmic scale: classical conjugate gradient must touch all N entries and scales as N, while HHL scales as a polynomial in log N. At N = 2²⁰ ≈ 1.05×10⁶ the classical cost is about a million matrix-vector operations while HHL runs in the qubit count n = 20 up to polylog and κ factors — an exponential gap in the dimension that widens without bound

The core equations

The problem is the most familiar one in computational science, with three structural assumptions attached. $A$ is Hermitian so that phase estimation is well-defined; sparse so that $e^{iAt_0}$ can be simulated in polylogarithmic time; and well-conditioned so the inversion does not blow up.

Ax=b,\qquad x=A^{-1}b,\qquad A=A^{\dagger} \qquad (1)

where: $A$ — $N\times N$ Hermitian, $s$ -sparse matrix with condition number $\kappa$ ; $b$ — right-hand side, loaded as $|b\rangle=\sum_i b_i|i\rangle/\lVert b\rVert$ ; $x$ — solution, delivered as $|x\rangle=\sum_i x_i|i\rangle/\lVert x\rVert$ ; $\kappa=|\lambda_{\max}|/|\lambda_{\min}|$ — condition number. Reading: a non-Hermitian $A$ is first embedded in a Hermitian dilation, so Hermiticity costs nothing in generality. Everything happens in the eigenbasis of $A$ , even though it is never computed. The spectral theorem hands us a real spectrum and an orthonormal eigenbasis, and the input decomposes along it.

A=\sum_{j=1}^{N}\lambda_j\,|u_j\rangle\langle u_j|,\qquad |b\rangle=\sum_{j=1}^{N}\beta_j\,|u_j\rangle,\qquad \beta_j=\langle u_j|b\rangle \qquad (2)

where: $\lambda_j$ — $j$ -th real eigenvalue of $A$ ; $|u_j\rangle$ — orthonormal eigenvector; $\beta_j$ — overlap of $|b\rangle$ with $|u_j\rangle$ . Reading: HHL processes each $|u_j\rangle$ branch independently and in parallel, which is why the target below is a term-by-term division. Inverting $A$ is nothing but dividing each spectral coefficient by its eigenvalue — the operation the controlled rotation will realize physically.

A^{-1}|b\rangle=\sum_{j=1}^{N}\beta_j\,\lambda_j^{-1}\,|u_j\rangle \qquad (3)

where: $\lambda_j^{-1}$ — the inverse eigenvalue, the multiplier that turns $b$ into $x$ along the $j$ -th axis. Reading: the whole algorithm is a physical device for writing the amplitudes $\beta_j\lambda_j^{-1}$ without ever isolating a single $|u_j\rangle$ . Phase estimation of the simulated unitary $e^{iAt_0}$ is the eigenvalue-reading engine: because $e^{iAt_0}|u_j\rangle=e^{i\lambda_j t_0}|u_j\rangle$ , it writes a binary estimate of $\lambda_j$ into a clock register entangled with $|u_j\rangle$ .

|u_j\rangle|0\rangle^{\otimes t}\;\xrightarrow{\ \mathrm{QPE}[e^{iAt_0}]\ }\;|u_j\rangle|\tilde\lambda_j\rangle \qquad (4)

where: $e^{iAt_0}$ — unitary from Hamiltonian simulation of $A$ ; $|\tilde\lambda_j\rangle$ — $t$ -qubit register storing $\lambda_j t_0/2\pi$ to $t$ bits; $t_0$ — evolution time fixed so $\lambda_j t_0\in[0,2\pi)$ , setting the eigenvalue resolution. Reading: $t$ bits of clock give eigenvalue precision $\sim 2^{-t}$ , and resolving $\lambda_{\min}\sim 1/\kappa$ needs $t\gtrsim\log_2(\kappa/\epsilon)$ qubits. The controlled rotation performs $\lambda\mapsto 1/\lambda$ . A flag ancilla is rotated by an angle read off the clock, so its $|1\rangle$ amplitude comes out inversely proportional to the eigenvalue.

|\tilde\lambda_j\rangle|0\rangle\;\to\;|\tilde\lambda_j\rangle\!\left(\sqrt{1-\tfrac{C^2}{\tilde\lambda_j^{2}}}\,|0\rangle+\tfrac{C}{\tilde\lambda_j}\,|1\rangle\right),\qquad \theta_j=2\arcsin\!\left(\tfrac{C}{\tilde\lambda_j}\right) \qquad (5)

where: $C$ — normalization constant obeying $C\le|\lambda_{\min}|\sim 1/\kappa$ so that $C/\tilde\lambda_j\le 1$ ; $|1\rangle$ — flag whose post-selection keeps the inverted branch. Reading: $\sin(\theta_j/2)=C/\tilde\lambda_j$ , so the flag amplitude is exactly proportional to $1/\lambda_j$ ; the cap $C\lesssim 1/\kappa$ is the sole reason the success probability is small. Uncomputing the clock and post-selecting the flag on $|1\rangle$ leaves the input register holding the solution, up to phase-estimation error.

|x\rangle\;\propto\;\sum_{j}\beta_j\,\frac{C}{\lambda_j}\,|u_j\rangle\;=\;C\,A^{-1}|b\rangle \qquad (6)

where: $|x\rangle$ — normalized solution state after uncomputation and post-selection. Reading: the surviving amplitudes are $\beta_j/\lambda_j$ up to the global constant $C$ , which is exactly $A^{-1}|b\rangle$ renormalized. The flag reads $|1\rangle$ with probability set by the branch weights, which is where $\kappa$ enters and where amplitude amplification buys it back.

P_{\mathrm{success}}=\sum_j|\beta_j|^2\frac{C^2}{\lambda_j^2}=\Theta\!\left(\frac{1}{\kappa^2}\right),\qquad N_{\mathrm{amp}}=O\!\left(\frac{1}{\sqrt{P_{\mathrm{success}}}}\right)=O(\kappa) \qquad (7)

where: $P_{\mathrm{success}}$ — probability of measuring the flag in $|1\rangle$ ; $N_{\mathrm{amp}}$ — number of amplitude-amplification rounds. Reading: naive repetition costs $O(\kappa^2)$ runs, but amplifying the flagged branch reduces this quadratically to $O(\kappa)$ — the linear-in- $\kappa$ factor of modern HHL. Assembling Hamiltonian simulation, phase estimation, and amplification gives the runtime, together with its later improvement and the honest classical comparison.

T_{\mathrm{HHL}}=\tilde O\!\left(\frac{s^2\kappa^2\log N}{\epsilon}\right)\xrightarrow{\ \mathrm{CKS}\ }\tilde O\!\big(s\kappa\,\mathrm{polylog}(N/\epsilon)\big),\qquad T_{\mathrm{classical}}=O\!\big(Ns\sqrt{\kappa}\,\log\tfrac1\epsilon\big) \qquad (8)

where: $s$ — row sparsity; $\epsilon$ — target precision; $\tilde O$ — soft-O suppressing log factors; $T_{\mathrm{classical}}$ — conjugate gradient on a sparse positive-definite system. Reading: the decisive change is $N\to\log N$ ; the useful output is a global property $\langle M\rangle=\langle x|M|x\rangle$ , since anything requiring all $N$ amplitudes reinstates the classical cost.

Derivation AFull HHL state evolution end to end

The goal is to track the complete quantum state — clock, input, and flag registers together — from the initial product $|0\rangle_c^{\otimes t}|b\rangle|0\rangle_f$ through phase estimation, the controlled rotation, uncomputation, and post-selection, and to prove that what survives is proportional to $A^{-1}|b\rangle$ . Working entirely in the eigenbasis of $A$ makes every step a diagonal action on the $|u_j\rangle$ branches, so the argument is exact rather than schematic, and it exposes exactly where phase-estimation error enters.

Step 1. Expand the input in the eigenbasis and set all ancillas to zero: the state is $\sum_j\beta_j\,|0\rangle_c|u_j\rangle|0\rangle_f$ . Why: any vector expands in the orthonormal eigenbasis, and writing it this way lets us follow each eigen-branch $j$ separately, since every gate below acts diagonally on $|u_j\rangle$ .

Step 2. Apply phase estimation of $e^{iAt_0}$ to the clock and input registers, giving $\sum_j\beta_j\,|\tilde\lambda_j\rangle_c|u_j\rangle|0\rangle_f$ . Why: $e^{iAt_0}|u_j\rangle=e^{i\lambda_j t_0}|u_j\rangle$ , so QPE writes the eigenphase into the clock; this is the only place the matrix $A$ enters the amplitudes, and it enters as an eigenvalue label, not as a stored eigenvector.

Step 3. Adjoin the flag and apply the controlled rotation $R_y(2\arcsin(C/\tilde\lambda_j))$ conditioned on the clock, producing $\sum_j\beta_j\,|\tilde\lambda_j\rangle_c|u_j\rangle\big(\sqrt{1-C^2/\tilde\lambda_j^2}\,|0\rangle_f+\tfrac{C}{\tilde\lambda_j}|1\rangle_f\big)$ . Why: the flag amplitude on $|1\rangle$ is now $C/\tilde\lambda_j$ , inversely proportional to the eigenvalue — the physical realization of the term-by-term division in equation (3).

Step 4. Run phase estimation in reverse to disentangle and reset the clock to $|0\rangle_c$ , leaving $\sum_j\beta_j\,|0\rangle_c|u_j\rangle\big(\sqrt{1-C^2/\tilde\lambda_j^2}\,|0\rangle_f+\tfrac{C}{\tilde\lambda_j}|1\rangle_f\big)$ . Why: the clock has done its job of controlling the rotation; uncomputing it un-entangles the eigenvalue register so the final state is a clean product of input and flag, at the cost of the small error incurred when $\tilde\lambda_j\neq\lambda_j$ .

Step 5. Measure the flag and condition on the outcome $|1\rangle$ , projecting onto $\frac{1}{\sqrt{\sum_j|\beta_j|^2 C^2/\tilde\lambda_j^2}}\sum_j\beta_j\tfrac{C}{\tilde\lambda_j}|u_j\rangle$ . Why: post-selection keeps only the inverted-amplitude branch and renormalizes it; the surviving amplitudes are $\propto\beta_j/\lambda_j$ , precisely the spectral coefficients of $A^{-1}|b\rangle$ .

Step 6. Identify the result: in the ideal limit $\tilde\lambda_j=\lambda_j$ the state is $|x\rangle=A^{-1}|b\rangle/\lVert A^{-1}|b\rangle\rVert$ . Why: dividing every eigen-coefficient of $|b\rangle$ by its eigenvalue and renormalizing is the definition of the normalized solution, so the circuit has solved $Ax=b$ in the amplitude encoding. *In the ideal limit $\tilde\lambda_j=\lambda_j$ the derivation reduces to exact inversion; on the 2×2 example it returns $|x\rangle=\sqrt{0.9}\,|0\rangle-\sqrt{0.1}\,|1\rangle$ .*

Go deeper · Derivation B · 1 minThe inversion rotation and the role of the constant C

Derivation BThe inversion rotation and the role of the constant C

The goal is to show precisely how a single controlled single-qubit rotation implements $\lambda\mapsto 1/\lambda$ , and why the normalization constant $C$ cannot be chosen freely but is pinned to $\sim 1/\kappa$ by the smallest eigenvalue. This is the hinge between the mechanism and the cost: the same constraint that makes the rotation valid is the one that caps the success probability, so understanding it is understanding where the $\kappa$ overhead is born.

Step 1. Apply $R_y(\theta_j)$ with $\theta_j=2\arcsin(C/\tilde\lambda_j)$ to a flag qubit in $|0\rangle$ , giving $\cos(\theta_j/2)|0\rangle+\sin(\theta_j/2)|1\rangle$ . Why: the $R_y$ rotation maps $|0\rangle$ to a real superposition whose $|1\rangle$ weight is $\sin(\theta_j/2)$ , so choosing $\theta_j$ lets us dial that weight to any value in $[0,1]$ .

Step 2. Evaluate the half-angle sine: $\sin(\theta_j/2)=\sin(\arcsin(C/\tilde\lambda_j))=C/\tilde\lambda_j$ . Why: this is the crux — the flag amplitude on $|1\rangle$ is exactly $C/\tilde\lambda_j$ , inversely proportional to the eigenvalue, so a large eigenvalue is suppressed and a small one enhanced, which is what inverting $A$ requires.

Step 3. Demand validity for every branch: $\arcsin$ needs $C/\tilde\lambda_j\le 1$ , hence $C\le\min_j\tilde\lambda_j\approx|\lambda_{\min}|$ . Normalizing the spectrum so $|\lambda_{\max}|\le 1$ gives $|\lambda_{\min}|=|\lambda_{\max}|/\kappa\sim 1/\kappa$ , so $C\sim 1/\kappa$ . Why: the smallest eigenvalue is the binding constraint; the whole inverse-condition-number factor enters here, purely to keep the rotation angle real for the hardest-to-invert eigen-direction.

Step 4. Read off the consequence: the flag amplitude for the smallest eigenvalue saturates near $1$ , but for typical eigenvalues it is $O(1/\kappa)$ , so the flagged branch carries weight $O(1/\kappa^2)$ . Why: the constant that makes the rotation legal simultaneously shrinks the post-selected probability, which is exactly why HHL's cost tracks $\kappa$ — the mechanism and its price are two faces of the choice of $C$ . *For the 2×2 example with $C=1$ : $\lambda=1$ gives $\theta=2\arcsin(1)=180°$ and flag amplitude $1$ , while $\lambda=2$ gives $\theta=2\arcsin(0.5)=60°$ and flag amplitude $0.5$ — each matching $C/\lambda$ exactly.*

Go deeper · Derivation C · 1 minCondition-number cost and amplitude amplification

Derivation CCondition-number cost and amplitude amplification

The goal is to derive the success probability of the flag measurement from the Born rule and then show how amplitude amplification converts the naive $O(\kappa^2)$ repetition count into $O(\kappa)$ , which is the single most important resource statement in modern analyses of HHL. The argument is short but it is where the algorithm's practical viability is decided, so each factor of $\kappa$ is tracked explicitly rather than absorbed into soft-O notation.

Step 1. Apply the Born rule to the flagged branch: $P_{\mathrm{success}}=\sum_j|\beta_j|^2(C/\lambda_j)^2$ , the summed squared amplitude on $|1\rangle$ . Why: post-selection succeeds with probability equal to the total weight of the $|1\rangle$ component, and because the eigen-branches are orthogonal these weights simply add.

Step 2. Insert $C\sim 1/\kappa$ and the normalization $|\lambda_j|\le 1$ : each term is at least $(C/\lambda_{\max})^2=\Theta(1/\kappa^2)$ , and since $\sum_j|\beta_j|^2=1$ the sum is $\Theta(1/\kappa^2)$ . Why: the constant $C$ from the previous derivation directly sets the scale; the probability is small precisely because $C$ had to be small to invert the smallest eigenvalue.

Step 3. Count naive repetitions: succeeding once from probability $p=\Theta(1/\kappa^2)$ takes $O(1/p)=O(\kappa^2)$ independent runs. Why: without amplification, post-selection is just rejection sampling, and the expected number of trials to a first success is the reciprocal of the per-trial probability.

Step 4. Amplify instead: amplitude amplification boosts a success amplitude $\sqrt p$ to $\Theta(1)$ in $O(1/\sqrt p)=O(\kappa)$ rounds, each one call to the HHL circuit and its inverse. Why: amplitude amplification is the Grover-style quadratic speedup applied to the flag, turning $O(\kappa^2)$ probabilistic repetitions into $O(\kappa)$ coherent iterations — the reason the headline runtime is linear, not quadratic, in $\kappa$ . *For the 2×2 example, $P_{\mathrm{success}}=\tfrac12(1)^2+\tfrac12(0.5)^2=0.625$ , equal to $C^2\lVert A^{-1}b\rVert^2$ with $C=1$ .*

Limiting cases

How HHL's cost and meaning shift at the edges of its assumptions:

Case	Result	Meaning
Well-conditioned	$\kappa=O(1)$	Success probability $\Theta(1)$ , no amplification needed; the $\log N$ dependence dominates and the exponential-in-dimension speedup is at its cleanest.
Ill-conditioned	$\kappa\to\infty$	Runtime $\propto\kappa^2$ (original) or $\kappa$ (improved); as $\lambda_{\min}\to 0$ the rotation angle saturates and $P_{\mathrm{success}}\to 0$ , so the speedup vanishes for near-singular systems.
High precision	$\epsilon\to 0$	Original HHL scales as $1/\epsilon$ from phase-estimation resolution; Childs–Kothari–Somma reach $\mathrm{polylog}(1/\epsilon)$ via a linear-combination-of-unitaries approximation of $1/x$ .
Dense matrix	$s\to N$	Sparsity enters the simulation of $e^{iAt_0}$ ; when $A$ is dense, $s=\Theta(N)$ makes the simulation cost $\mathrm{poly}(N)$ and the exponential advantage over classical solvers disappears.
Full read-out demanded	tomography of all $x_i$	Recovering all $N$ amplitudes needs $\Omega(N)$ measurements, restoring linear-in- $N$ cost; the advantage holds only for an observable $\langle x\|M\|x\rangle$ with an efficiently implementable $M$ .

Common trapHHL hands you the classical solution vector x

It does not. HHL outputs a quantum state $|x\rangle$ whose amplitudes are the normalized entries $x_i$ of the solution — not a list of numbers you can print. Reading out all $N$ amplitudes requires quantum state tomography, which costs $\Omega(N)$ measurements and immediately erases the exponential speedup, dragging you back to classical scaling. The advantage is genuine only when the question you ask about $x$ is itself "small" — a single expectation value $\langle x|M|x\rangle$ , a few moments, an overlap with a fixed state, or a binary decision — and only when $\kappa$ , $s$ , and $1/\epsilon$ are modest and $|b\rangle$ can be prepared efficiently and $A$ is accessible by a coherent oracle. Strip away any of these and the "exponential speedup" is a mirage. HHL is not a faster way to get the vector $x$ ; it is a way to learn one global fact about $x$ without ever forming it.

Worked examples

1 · A 2×2 system solved by hand, end to end (fully worked)

Solve $Ax=b$ with $A=\begin{pmatrix}1.5&0.5\\0.5&1.5\end{pmatrix}$ and $b=(1,0)^{\top}$ on the HHL circuit: find the eigenbasis, the rotation angles, the output state, and the success probability. Take $C=1$ .

Step 1. Diagonalize $A$ . Its eigenvalues are $\lambda=1.5\pm0.5=\{2,1\}$ with eigenvectors $|u_1\rangle=|+\rangle=(|0\rangle+|1\rangle)/\sqrt2$ (for $\lambda_1=2$ ) and $|u_2\rangle=|-\rangle=(|0\rangle-|1\rangle)/\sqrt2$ (for $\lambda_2=1$ ). The condition number is $\kappa=\lambda_{\max}/\lambda_{\min}=2/1=\mathbf{2}$ — very well-conditioned.

Step 2. Load $b$ . Normalized, $|b\rangle=|0\rangle=\tfrac{1}{\sqrt2}(|+\rangle+|-\rangle)$ , so the eigen-overlaps are $\beta_1=\beta_2=\mathbf{1/\sqrt2}$ . Both eigen-branches start with equal weight.

Step 3. Rotate to invert. The angles are $\theta_j=2\arcsin(C/\lambda_j)$ : for $\lambda_2=1$ , $\theta=2\arcsin(1)=\mathbf{180°}$ (flag amplitude $C/\lambda=1$ ); for $\lambda_1=2$ , $\theta=2\arcsin(0.5)=\mathbf{60°}$ (flag amplitude $0.5$ ). The smaller eigenvalue gets the larger amplitude, as inversion demands.

Step 4. Assemble the post-selected output. The flagged branch is $\propto\sum_j\beta_j(C/\lambda_j)|u_j\rangle=\tfrac{1}{\sqrt2}\big(\tfrac12|+\rangle+1\cdot|-\rangle\big)$ . Converting to the computational basis gives $0.75|0\rangle-0.25|1\rangle$ , which is exactly $A^{-1}b=(0.75,-0.25)^{\top}$ . Normalizing ( $\lVert x\rVert=\sqrt{0.625}=0.7906$ ) yields $|x\rangle=\mathbf{\sqrt{0.9}\,|0\rangle-\sqrt{0.1}\,|1\rangle}$ , i.e. amplitudes $0.9487$ and $-0.3162$ .

Step 5. Success probability. $P_{\mathrm{success}}=\sum_j|\beta_j|^2(C/\lambda_j)^2=\tfrac12(1)^2+\tfrac12(0.5)^2=0.5+0.125=\mathbf{0.625}$ . Because $\kappa=2$ is tiny, the flag reads $1$ on 62.5% of runs — no amplification required. Result: $|x\rangle=\sqrt{0.9}\,|0\rangle-\sqrt{0.1}\,|1\rangle$ , recovered with probability 62.5%, exactly matching $A^{-1}b=(0.75,-0.25)$ .

2 · Sizing the exponential speedup (N = 2³⁰)

A sparse system has $N=2^{30}\approx1.07\times10^9$ unknowns, row sparsity $s=10$ , condition number $\kappa=100$ , and you want one expectation value $\langle x|M|x\rangle$ to precision $\epsilon=10^{-2}$ . Estimate the register sizes and compare the quantum and classical costs.

Step 1. Qubit budget. The solution register needs $n=\log_2 N=\mathbf{30}$ qubits — the entire billion-dimensional vector lives here. The clock needs $t\approx\log_2(\kappa/\epsilon)=\log_2(10^4)\approx\mathbf{14}$ qubits to resolve eigenvalues down to $\lambda_{\min}\sim1/\kappa$ , plus one flag ancilla: about $\mathbf{45}$ qubits total to represent a $10^9$ -dimensional linear solve.

Step 2. Classical cost. Conjugate gradient must perform matrix-vector products over all $N$ entries: $T_{\mathrm{classical}}\approx Ns\sqrt{\kappa}\approx(1.07\times10^9)(10)(10)\approx\mathbf{1.1\times10^{11}}$ operations — every one of the billion components must be touched at least once.

Step 3. Quantum cost. HHL scales as $\mathrm{poly}(\log N)$ in the dimension: with the improved $\tilde O(s\kappa\,\mathrm{polylog}(N/\epsilon))$ bound and $\log_2 N=30$ , the dimension contributes only polynomially in $30$ , not in $10^9$ . The $N\to\log N$ collapse is the whole point: a factor of $\sim N/\mathrm{poly}(\log N)\approx 10^9/10^3\sim\mathbf{10^6}$ in the dimension dependence alone.

Step 4. The catch, quantified. That $10^6$ -fold dimension advantage survives only if $M$ is efficiently measurable, $|b\rangle$ is efficiently preparable, and the oracle for $A$ is cheap; and it is throttled by $\kappa=100$ (a factor of $100$ overhead) and by $s=10$ . Ask for the full vector instead of $\langle x|M|x\rangle$ and the $\Omega(N)\approx10^9$ read-out cost erases everything. The dimension speedup is exponential ( $N\to\log N$ ), but it is realized only through a global observable and only when the input/output and conditioning are benign.

Physics of the portrait

HHL is a theorem in computational complexity, yet it sits in every quantum-information course because it drew the boundary of quantum advantage for the most practical task in all of numerical mathematics — and then a decade of "fine print" mapped exactly how sharp that boundary is.

Origins · 1 minAram Harrow and how this was discovered

Aram Harrow

b. 1980 · MIT · quantum algorithms and the theory of quantum information

Aram W. Harrow, together with Avinatan Hassidim and Seth Lloyd, published Quantum Algorithm for Linear Systems of Equations in Physical Review Letters in 2009, and the three initials gave the algorithm its name. The construction was a conceptual jolt: for two decades the flagship quantum algorithms — Deutsch–Jozsa, Simon, Shor, Grover — attacked number-theoretic or search problems, but HHL reached into the daily machinery of engineering and science, the linear solve at the heart of every finite-element model and every least-squares fit, and claimed an exponential speedup in the dimension. Harrow, a physicist trained at MIT under Ike Chuang and now a professor there, has spent his career on precisely the questions HHL forces: what quantum computers can and cannot do faster, and how to prove it. The algorithm's reception drove its own literature. Scott Aaronson's 2015 essay Read the fine print laid out the four caveats — state preparation, condition number, sparsity, and read-out — that any honest speedup claim must clear, and Ewin Tang's 2018 "dequantization" showed that some HHL-based machine-learning proposals had classical analogues after all, sharpening rather than shrinking the result. Meanwhile Ambainis, and later Childs, Kothari, and Somma, improved the $\kappa^2$ and $1/\epsilon$ dependences, and the modern framework of block-encoding and quantum singular value transformation absorbed HHL as a special case. What began as one surprising circuit became the founding document of quantum linear algebra.

History

1996

— Kitaev frames eigenvalue estimation as quantum phase estimation, the engine HHL will use

2009

— Harrow, Hassidim, and Lloyd publish the quantum algorithm for linear systems of equations

2015

— Aaronson's essay Read the fine print catalogues the input, conditioning, and read-out caveats

2017

— Childs, Kothari, and Somma improve the precision dependence to polylog(1/ε)

2018

— Tang dequantizes several HHL-based recommendation and machine-learning speedups

Going further

Texts. Harrow, Hassidim, Lloyd, Quantum Algorithm for Linear Systems of Equations, Phys. Rev. Lett. 103, 150502 (2009) · Nielsen & Chuang, Quantum Computation and Quantum Information (10th anniversary ed.), Ch. 4–5 (phase estimation and Hamiltonian simulation) · Preskill, Lecture Notes for Physics 219, Ch. 6 · Childs, Kothari, Somma, Quantum algorithm for systems of linear equations with exponentially improved dependence on precision, SIAM J. Comput. 46, 1920 (2017) · Aaronson, Read the fine print, Nature Physics 11, 291 (2015) · Dervovic et al., Quantum linear systems algorithms: a primer, arXiv:1802.08227 (2018). The rotation identity $\sin(\theta_j/2)=C/\tilde\lambda_j$ with $\theta_j=2\arcsin(C/\tilde\lambda_j)$ is elementary trigonometry. The success probability $\Theta(1/\kappa^2)$ and the $O(\kappa)$ amplification count follow from the Brassard–Høyer–Mosca–Tapp amplitude-amplification bound. The 2×2 worked example ( $A$ with eigenvalues $1,2$ ; $\kappa=2$ ; angles $180°,60°$ ; $|x\rangle=\sqrt{0.9}\,|0\rangle-\sqrt{0.1}\,|1\rangle$ ; $P_{\mathrm{success}}=0.625$ ) and the $N=2^{30}$ resource estimate ( $n=30$ , $t\approx14$ clock qubits, classical $\approx1.1\times10^{11}$ operations) were checked numerically. The improved runtime $\tilde O(s\kappa\,\mathrm{polylog}(N/\epsilon))$ is from Childs–Kothari–Somma (2017).

Prerequisites. Quantum Phase Estimation · Quantum Fourier Transform · Amplitude Amplification. Leads to. Grover's Search Algorithm · Amplitude Amplification.

Ask the tutor: Ask me to run the 2×2 example through every register — phase estimation, the $180°$ and $60°$ rotations, and the post-selection that returns $\sqrt{0.9}\,|0\rangle-\sqrt{0.1}\,|1\rangle$ — or to show how the fine-print caveats decide whether a real-world linear solve keeps its exponential speedup.

Does HHL solve linear systems faster than my laptop?

For getting the actual answer vector, no. HHL never produces the list of numbers $x_1,\dots,x_N$ ; it produces a quantum state $|x\rangle$ whose amplitudes are those numbers. If you need all of them, you must measure the state many times — tomography costing $\Omega(N)$ — which is no better than a classical solve. HHL is faster only when you want a single global property of the solution, such as an expectation value $\langle x|M|x\rangle$ , an overlap with a reference state, or a yes/no answer, and only when the matrix is sparse and well-conditioned and the input is loadable efficiently. It is best thought of as a way to learn one fact about the solution to a gigantic system without ever forming the solution.

Why must A be Hermitian, and what if it is not?

Hermiticity guarantees a real spectrum and a complete orthonormal eigenbasis, which is exactly what makes phase estimation of $e^{iAt_0}$ well-defined — the unitary has eigenphases $e^{i\lambda_j t_0}$ that the clock register can read. If $A$ is not Hermitian, you do not give up: you embed it in the Hermitian dilation $\tilde A=\begin{pmatrix}0&A\\A^{\dagger}&0\end{pmatrix}$ and solve $\tilde A\,y=\binom{b}{0}$ , whose solution has $A^{-1}b$ sitting in one block. The dilation doubles the dimension (one extra qubit) and leaves the sparsity and condition number essentially unchanged, so it is a free move in the asymptotics. Hermiticity is therefore a convenience of exposition, not a real restriction on the problems HHL can address.

What exactly is the condition number, and why does it hurt so much?

The condition number $\kappa=|\lambda_{\max}|/|\lambda_{\min}|$ measures how close $A$ is to singular: a large $\kappa$ means some eigen-direction is barely stretched, so inverting it demands a huge amplitude. In HHL that shows up twice. First, the inversion constant must satisfy $C\lesssim|\lambda_{\min}|\sim1/\kappa$ so the rotation angle stays valid for the smallest eigenvalue, which shrinks the success probability to $\Theta(1/\kappa^2)$ . Second, phase estimation must resolve eigenvalues down to $\sim1/\kappa$ , needing $\sim\log_2\kappa$ extra clock qubits and longer evolution times. Amplitude amplification softens the first cost from $\kappa^2$ to $\kappa$ , but nothing removes the $\kappa$ entirely — an ill-conditioned system is genuinely hard, quantum or classical, because near-singularity is a real obstruction, not an artifact of the algorithm.

Where does the exponential speedup actually come from?

From compressing dimension into qubits. An $N$ -component vector lives in the amplitudes of $\log_2 N$ qubits, so a million-dimensional problem uses twenty qubits and a billion-dimensional one uses thirty. Phase estimation then acts on all $N$ eigen-components in a single superposition, and the controlled rotation inverts every eigenvalue at once — the operations a classical machine would do one eigenvalue at a time happen in parallel across the whole spectrum. That is the source of the $N\to\log N$ collapse. The catch is symmetric: the same compression that lets you hold $N$ numbers cheaply prevents you from reading them cheaply. So the speedup is real for the compute step and forbidden for the read-out step, which is why HHL is useful only for questions whose answers are themselves small.

How does HHL improve on the original κ² and 1/ε scaling?

Two separate refinements. The $\kappa^2$ in the original success probability becomes $\kappa$ through amplitude amplification — Ambainis' variable-time amplitude amplification handles the fact that different eigen-branches take different times to invert, squeezing the repetition count down without spoiling the branches that finish early. The $1/\epsilon$ from phase-estimation resolution and Hamiltonian-simulation error becomes $\mathrm{polylog}(1/\epsilon)$ in the Childs–Kothari–Somma algorithm, which sidesteps phase estimation entirely and instead approximates the function $1/x$ directly by a linear combination of unitaries built from a Chebyshev or Fourier expansion. The modern statement, $\tilde O(s\kappa\,\mathrm{polylog}(N/\epsilon))$ , folds both improvements together and is what block-encoding and quantum singular value transformation reproduce as a special case.

Which quantum subroutines does HHL call, and in what order?

Five stages, in sequence. First, Hamiltonian simulation turns the matrix $A$ into the unitary $e^{iAt_0}$ — sparsity and oracle access are what make this efficient. Second, quantum phase estimation (which ends in an inverse quantum Fourier transform) writes each eigenvalue into the clock register. Third, a controlled $R_y$ rotation on a flag ancilla imprints the amplitude $C/\lambda_j$ , performing the inversion. Fourth, the inverse of phase estimation uncomputes and resets the clock so the eigenvalue register disentangles. Fifth, amplitude amplification on the flag boosts the post-selection success probability from $\Theta(1/\kappa^2)$ toward $\Theta(1)$ . This is why HHL sits downstream of phase estimation, the Fourier transform, and amplitude amplification in the dependency graph — it is a composite of all three.