The HHL Algorithm (Introduction)

The goal is to solve $Ax = b$ where $A$ is an $N \times N$ Hermitian matrix and $b \in \mathbb{C}^N$. Encode $b$ as a quantum state $|b\rangle$ on $n = \log_2 N$ qubits, normalized so $\langle b | b \rangle = 1$. The target is the state

Since $A$ is Hermitian it has a spectral decomposition $A = \sum_j \lambda_j |u_j\rangle\langle u_j|$ with real eigenvalues $\lambda_j$ and orthonormal eigenvectors $|u_j\rangle$. Writing $|b\rangle$ in this eigenbasis as $|b\rangle = \sum_j \beta_j |u_j\rangle$, the formal inverse acts eigenvalue-wise:

HHL's central trick is to implement this eigenvalue-wise inversion using Quantum Phase Estimation to read $\lambda_j$ into an ancilla register, a controlled rotation conditioned on that ancilla to multiply the amplitude by something proportional to $1/\lambda_j$, and then uncomputation of the ancilla so the eigenvalue information is discarded while its effect on the amplitude survives. This is a coherent, reversible way of applying a non-unitary operator ($A^{-1}$ is not unitary in general) inside a unitary quantum circuit, using post-selection on an ancilla flag. If $A$ is not Hermitian, a standard embedding restores Hermiticity: define

which is Hermitian, and solve $\tilde{A} \tilde{x} = \tilde{b}$ with $\tilde{b} = (b, 0)^T$; the lower half of the solution recovers $A^{-1}b$ up to sign convention. This means Hermiticity is not really a restriction — it costs only a qubit and a reindexing.

HHL uses three registers: a clock register of $t$ qubits for phase estimation, a single ancilla qubit for the rotation, and the $n$-qubit input register holding $|b\rangle$ and eventually $|x\rangle$. The algorithm assumes access to two black-box subroutines: a way to prepare $|b\rangle$, and a way to implement $e^{iAt}$ for controllable time $t$ — Hamiltonian simulation of $A$. Both are assumed efficient; this assumption is itself a major caveat discussed below.

Step 1 — State preparation. Prepare the input register in state $|b\rangle$. In the eigenbasis of $A$ this is $|b\rangle = \sum_{j} \beta_j |u_j\rangle$, though the algorithm never needs to know the $\beta_j$ or $|u_j\rangle$ explicitly.

Step 2 — Phase estimation. Apply quantum phase estimation using the unitary $U = e^{iAt_0}$ for a fixed simulation time $t_0$. Since $|u_j\rangle$ is an eigenvector of $A$ with eigenvalue $\lambda_j$, it is an eigenvector of $U$ with eigenvalue $e^{i\lambda_j t_0}$. Phase estimation writes an approximation $\tilde\lambda_j$ of $\lambda_j t_0 / (2\pi)$ into the clock register, entangling it with $|u_j\rangle$:

Step 3 — Controlled rotation. Introduce the ancilla qubit initialized to $|0\rangle$, and apply a controlled rotation conditioned on the clock register value $\tilde\lambda_j$:

where $C$ is a normalization constant chosen smaller than the smallest eigenvalue magnitude to keep the square root real. This is implementable by a sequence of controlled $Ry$ rotations, one per possible clock value, each computing the needed angle $\theta(\tilde\lambda_j) = 2\arcsin(C/\tilde\lambda_j)$ classically ahead of time.

Step 4 — Uncompute phase estimation. Apply the inverse of the phase estimation circuit. Because the clock register was only used to control the rotation, and the rotation did not alter which eigenvector branch we are in, inverse phase estimation disentangles the clock register back to $|0\rangle_{\text{clock}}$ (up to phase-estimation approximation error), leaving

Step 5 — Measure the ancilla. Measure the ancilla qubit. Conditioned on outcome $|1\rangle$, which occurs with probability $P_1 = \sum_j |\beta_j|^2 C^2/\lambda_j^2$, the input register collapses (up to normalization) to

which is exactly $|x\rangle$ up to normalization. This is the crux: postselecting on the ancilla is what implements the non-unitary inversion inside an otherwise unitary circuit. ``` Algorithm HHL(A, b, t): Input: A : N x N Hermitian matrix, N = 2^n, sparsity s, condition number kappa b : classical description sufficient to prepare state |b> on n qubits t : number of clock qubits for phase estimation (precision control) Output: Quantum state |x> proportional to A^{-1}|b> (NOT a classical vector) 1. Allocate registers: clock <- t qubits, initialized |0...0> anc <- 1 qubit, initialized |0> input <- n qubits, initialized |0...0> 2. Prepare input register: input <- |b> // via a state-preparation oracle 3. Phase estimation of U = exp(i A t0): clock, input <- PhaseEstimation(U, input, clock) // entangles clock with eigenbasis of A: // sum_j beta_j |u_j> |lambda_tilde_j> 4. Controlled rotation on ancilla, conditioned on clock register value L: for each computational basis value L of clock register: theta(L) <- 2 * arcsin( C / lambda_from(L) ) // precomputed classically apply controlled-Ry(theta(L)) on anc, controlled by clock == L 5. Uncompute phase estimation: clock, input <- InversePhaseEstimation(U, input, clock) // clock returns (approximately) to |0...0> 6. Measure anc in computational basis: if outcome == 1: input register now holds state proportional to A^{-1}|b> = |x> ACCEPT else: discard and retry from step 2 // postselection failure 7. Return input register (as a quantum state, to be consumed by a subsequent expectation-value estimation circuit, not read out directly) ```

Correctness argument. The proof rests on linearity and the eigenbasis decomposition. Because $A$ is Hermitian, $\{|u_j\rangle\}$ forms an orthonormal basis, so every step of the circuit can be analyzed eigenvector-by-eigenvector and then summed by linearity of the unitary operations. Phase estimation is exact in the limit of infinite clock qubits: it maps $|u_j\rangle|0\rangle \to |u_j\rangle|\lambda_j t_0/2\pi\rangle$ exactly when $\lambda_j t_0/2\pi$ is exactly representable in $t$ bits, and approximately otherwise, with error analyzed in the Quantum Phase Estimation article. The controlled rotation is a deterministic, coherent function of the clock register value, so it commutes correctly with the entangled superposition — each branch $j$ receives its own correctly computed angle $\theta(\tilde\lambda_j)$. Uncomputing phase estimation is valid because the rotation did not touch the clock or input registers, only entangled a new ancilla with the already-existing clock value; inverting the exact same unitary used to create the clock entanglement returns the clock to $|0\rangle$ in each branch, leaving the eigenvector label undisturbed. Finally, projecting the ancilla onto $|1\rangle$ selects exactly the $1/\lambda_j$ amplitude component from each branch, which is by construction the action of $A^{-1}$ in the eigenbasis. No step requires knowing the eigenvalues or eigenvectors classically; they only ever appear encoded in superposition.

Success probability. The probability of measuring ancilla $=1$ is

Choosing $C = O(1/\kappa)$ where $\kappa = \lambda_{\max}/\lambda_{\min}$ is the condition number keeps the rotation argument valid for the smallest eigenvalue, but it also makes $P_1 = O(1/\kappa^2)$ in the worst case. Amplitude amplification, using the same reasoning as Grover's Algorithm, boosts this to $O(1/\kappa)$ repetitions instead of $O(1/\kappa^2)$ by amplifying the amplitude linearly instead of relying on repeated sampling of the squared probability.

Time complexity derivation. Three costs stack up. 1. Hamiltonian simulation. Implementing $e^{iAt_0}$ for an $s$-sparse $N \times N$ Hermitian matrix costs $O(s^2 t_0 \, \text{poly}(\log N)/\epsilon_{\text{sim}})$ using standard sparse Hamiltonian simulation techniques (product formulas or more advanced methods), where $\epsilon_{\text{sim}}$ is the simulation error tolerance. 2. Phase estimation precision. To resolve eigenvalues to relative precision sufficient for the controlled rotation to be accurate, the clock register needs $t = O(\log(\kappa/\epsilon))$ qubits, and phase estimation must run for evolution time $t_0 = O(\kappa/\epsilon)$ to distinguish eigenvalues spanning the range down to $\lambda_{\min} \sim \lambda_{\max}/\kappa$. 3. Amplitude amplification. Boosting success probability from $O(1/\kappa^2)$ to $O(1)$ costs a multiplicative factor of $O(\kappa)$ in repetitions of the whole circuit. Multiplying these factors, the total runtime is

in the original 2009 analysis (later refinements using improved Hamiltonian simulation and variable-time amplitude amplification reduce the $\kappa$ dependence toward linear, and the $\epsilon$ dependence toward logarithmic, but the qualitative shape — polylog in $N$, polynomial in $\kappa$ and $1/\epsilon$ — persists). The key point is the $\log N$ scaling: for classical general-purpose Gaussian elimination, cost is $O(N^3)$; even fast sparse conjugate-gradient methods cost $O(N s \kappa \log(1/\epsilon))$. HHL trades polynomial dependence on $N$ for polynomial dependence on $\kappa$ and $1/\epsilon$, which is a genuine exponential speedup in $N$ when $\kappa$ and $1/\epsilon$ are modest.

Space complexity. The circuit uses $n + t + 1$ qubits: $n = \log_2 N$ for the input register, $t = O(\log(\kappa/\epsilon))$ for the clock, and one ancilla. This is exponentially smaller than the $N$ classical numbers needed to write down $x$ explicitly — which is exactly why HHL cannot simply hand you $x$ without destroying the speedup.

Honest tradeoffs.

• The speedup is conditional on $A$ being efficiently simulable — dense, unstructured matrices with no exploitable sparsity or structure do not admit efficient Hamiltonian simulation, and the whole algorithm stalls at the first step.
• The speedup degrades quadratically (pre-amplification) or linearly (post-amplification) in the condition number $\kappa$; for ill-conditioned systems ($\kappa$ growing with $N$), the advantage can vanish entirely.
• State preparation of $|b\rangle$ is assumed efficient; if $b$ is an arbitrary classical vector with no succinct description, preparing $|b\rangle$ can itself cost $O(N)$, erasing the speedup before the algorithm even starts.
• Reading out $x$ fully costs $O(N)$ via state tomography, again erasing the advantage — usable output is restricted to expectation values $\langle x|M|x\rangle$ for some observable $M$, or sampling-based summary statistics.

Take a small, fully traceable instance. Let

$A$ is Hermitian (real symmetric) with eigenvalues and eigenvectors:

Condition number $\kappa = \lambda_1/\lambda_2 = 2$ — a friendly, well-conditioned case, chosen deliberately non-symmetric in its input so the eigen-decomposition step is non-trivial. Expand $b$ in the eigenbasis:

Convert back to the computational basis:

Check against classical inversion: $A^{-1} = \frac{1}{8}\begin{pmatrix}3 & -1\\ -1 & 3\end{pmatrix}$, so

which normalizes to exactly $\frac{1}{\sqrt{10}}(3,-1)^T$ — matching the HHL output state. This confirms the algorithm reproduces the correctly normalized solution direction; the overall scale (the $1/8$ classical prefactor) is discarded, because quantum states are normalized by construction. Any application of HHL must separately extract this scale (e.g. via the ancilla acceptance probability $P_1$, which is proportional to $\|A^{-1}b\|^2$) if the magnitude of $x$ matters and not just its direction.

HHL is foundational to the broader research program of quantum linear algebra, which underlies many proposed quantum machine learning algorithms: quantum support vector machines, quantum principal component analysis, and quantum recommendation systems all use HHL-style matrix inversion as a subroutine, extracting an expectation value or classification decision rather than the raw solution vector. It also motivates ongoing research into Hamiltonian simulation methods, since HHL's practical runtime is bottlenecked by how efficiently $e^{iAt}$ can be implemented for the matrices that arise in applications — a topic tightly linked to Quantum Simulation more broadly. In finance and engineering, HHL-inspired algorithms have been proposed for solving discretized partial differential equations (e.g. in computational fluid dynamics) where the matrix $A$ comes from a sparse finite-difference stencil with modest condition number, precisely the regime where the algorithm's assumptions are most plausible. Equally important is HHL's role as a cautionary case study: the 2018 "dequantization" results, showing certain claimed speedups can be matched classically under comparable sampling assumptions, taught the field to scrutinize input/output models as carefully as the quantum circuit itself — a methodological lesson that now shapes how every new quantum algorithm, including those built on Quantum Phase Estimation, is evaluated for genuine advantage.

Ask the tutor: How does amplitude amplification concretely boost HHL's success probability from $O(1/\kappa^2)$ to $O(1/\kappa)$?