The Harrow–Hassidim–Lloyd (HHL) algorithm is a quantum algorithm for obtaining certain limited information about the solution to a system of linear equations, introduced by Aram Harrow, Avinatan Hassidim, and Seth Lloyd. Specifically, the algorithm estimates quadratic functions of the solution vector to a given system.

The algorithm is one of the main fundamental algorithms expected to provide a speedup over their classical counterparts, along with Shor's factoring algorithm and Grover's search algorithm. Assuming the system is sparse, has a low condition number ${\displaystyle \kappa }$, and that the user is only interested in certain information about solution vector and not the entire vector itself, the algorithm has a runtime of ${\displaystyle O(\log(N)\kappa ^{2})}$, where ${\displaystyle N}$ is the number of variables. This offers an exponential speedup over the fastest classical algorithm, which runs in ${\displaystyle O(N\kappa )}$ (or ${\displaystyle O(N{\sqrt {\kappa }})}$ for positive semidefinite matrices).

An implementation of the HHL algorithm was first demonstrated in 2013 by three independent publications, consisting of simple systems on specially designed devices. The first demonstration of a general-purpose version of the algorithm appeared in 2018.

Given an ${\displaystyle N\times N}$ Hermitian matrix ${\displaystyle A}$ and unit vector ${\displaystyle {\vec {b}}\in \mathbb {R} ^{N}}$, the HHL algorithms prepares the quantum state ${\displaystyle |x\rangle }$ whose amplitudes are the entries of the solution ${\displaystyle {\vec {x}}\in \mathbb {R} ^{N}}$ to the linear system ${\displaystyle A{\vec {x}}={\vec {b}}}$. The algorithm cannot efficiently output the solution x itself, but allows one to efficiently estimate ${\displaystyle {\vec {x}}^{T}M{\vec {x}}}$ for a Hermitian matrix ${\displaystyle M}$.

The algorithm first prepares the quantum state ${\displaystyle |b\rangle }$ whose amplitudes are equal to the entries of ${\displaystyle {\vec {b}}}$. Using Hamiltonian simulation, the unitary operator ${\displaystyle e^{iAt}}$ is applied to ${\displaystyle |b\rangle }$ for a superposition of different times t. The algorithm then uses quantum phase estimation to decompose ${\displaystyle |b\rangle }$ in the eigenbasis of ${\displaystyle A}$ and find the corresponding eigenvalues ${\displaystyle \lambda _{j}}$. The state of the system after this step is approximately

where ${\displaystyle u_{j}}$ are the eigenvectors of A and ${\displaystyle \beta _{j}}$ is the j-th coefficient of b in the eigenbasis of A.

We would then like to apply the linear map taking ${\displaystyle |\lambda _{j}\rangle }$ to ${\displaystyle C\lambda _{j}^{-1}|\lambda _{j}\rangle }$ for some constant C. This map is not unitary and must be implemented using a quantum measurement with a nonzero probability of failure. After it succeeds, we have uncomputed the ${\displaystyle |\lambda _{j}\rangle }$ register and have a state proportional to

By performing the quantum measurement corresponding to M, we get an estimate of ${\displaystyle {\vec {x}}^{T}M{\vec {x}}}$. One could use quantum tomography to retrieve all components of x, but this would require repeating the algorithm roughly N times.