HHL algorithm

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.^[1]

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,^[2] 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.^[3][4][5] The first demonstration of a general-purpose version of the algorithm appeared in 2018.^[6]

Overview

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

${\displaystyle \sum _{j\mathop {=} 1}^{N}\beta _{j}|u_{j}\rangle |\lambda _{j}\rangle ,}$

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

${\displaystyle \sum _{i\mathop {=} 1}^{N}\beta _{i}\lambda _{j}^{-1}|u_{j}\rangle =A^{-1}|b\rangle =|x\rangle .}$

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.

Detailed description

Assumptions and initialization

The algorithm requires the following assumptions to hold:

1. The algorithm requires A to be Hermitian so that it can be exponentiated into a unitary operator. If A is not Hermitian, one can define a Hermitian matrix ${\displaystyle \mathbf {C} ={\begin{bmatrix}0&A\\A^{\dagger }&0\end{bmatrix}}}$ and solve ${\displaystyle Cy={\begin{bmatrix}b\\0\end{bmatrix}}}$ to obtain ${\displaystyle y={\begin{bmatrix}0\\x\end{bmatrix}}}$.
2. The algorithm requires an efficient procedure to prepare ${\displaystyle |b\rangle }$. It is assumed that either ${\displaystyle |b\rangle }$ has already been prepared or there exists some B which takes some quantum state ${\displaystyle |\mathrm {initial} \rangle }$ to ${\displaystyle |b\rangle }$ efficiently. Any error in the preparation of ${\displaystyle |b\rangle }$ is ignored.
3. The algorithm assumes that the state ${\displaystyle |\psi _{0}\rangle }$ can be prepared efficiently, where ${\displaystyle |\psi _{0}\rangle$ :={\sqrt {2/T}}\sum _{\tau \mathop {=} 0}^{T-1}\sin \pi \left({\tfrac {\tau +{\tfrac {1}{2}}}{T}}\right)|\tau \rangle } for some large T. The coefficients of ${\displaystyle |\psi _{0}\rangle }$ are chosen to minimize a certain quadratic loss function which induces error in the ${\displaystyle U_{\mathrm {invert} }}$ subroutine described below.
4. The algorithm assumes that the unitary operator ${\displaystyle e^{iAt}}$ can be applied efficiently. This is possible using Hamiltonian simulation if A is s-sparse and efficiently row computable, meaning it has at most s nonzero entries per row which can be computed in time O(s) given a row index. One can then apply ${\displaystyle e^{iAt}}$ in time ${\displaystyle O(\log(N)s^{2}t)}$.

U_invert subroutine

The key subroutine to the algorithm, denoted ${\displaystyle U_{\mathrm {invert} }}$, is defined as follows using phase estimation:

1. Prepare ${\displaystyle |\psi _{0}\rangle ^{C}}$ on register C
2. Apply the conditional Hamiltonian evolution (sum)
3. Apply the Fourier transform to the register C. Denote the resulting basis states with ${\displaystyle |k\rangle }$ for k = 0, ..., T − 1. Define ${\displaystyle \lambda _{k}:=2\pi k/t_{0}}$.
4. Adjoin a three-dimensional register S in the state

${\displaystyle |h(\lambda _{k})\rangle ^{S}:={\sqrt {1-f(\lambda _{k})^{2}-g(\lambda _{k})^{2}}}|\mathrm {nothing} \rangle ^{S}+f(\lambda _{k})|\mathrm {well} \rangle ^{S}+g(\lambda _{k})|\mathrm {ill} \rangle ^{S},}$

1. Reverse steps 1–3, uncomputing any garbage produced along the way.

The phase estimation procedure in steps 1-3 estimates the eigenvalues of A up to error ${\displaystyle \epsilon }$.

The ancilla register in step 4 is needed to construct a state with inverted eigenvalues corresponding to the diagonalized inverse of A. The states 'nothing', 'well' and 'ill' are used to direct the loop body; 'nothing' indicates that the matrix inversion has not yet taken place, 'well' indicates that it has and the loop should halt, and 'ill' indicates that part of ${\displaystyle |b\rangle }$ is in the ill-conditioned subspace of A and the algorithm cannot produce the desired inversion. Producing a state proportional to the inverse of A requires 'well' to be measured, after which the overall state collapses to the desired output.

Main loop

The main loop follows amplitude amplification: starting with ${\displaystyle U_{\mathrm {invert} }B|\mathrm {initial} \rangle }$, repeatedly apply

${\displaystyle U_{\mathrm {invert} }B(I-2|\mathrm {initial} \rangle \langle \mathrm {initial} |)B^{\dagger }U_{\mathrm {invert} }^{\dagger }(I-2|\mathrm {well} \rangle \langle \mathrm {well} |).}$

After each iteration, ${\displaystyle S}$ is measured and will produce a value of 'nothing', 'well', or 'ill.' The loop is repeated until 'well' is measured, which occurs with some probability ${\displaystyle p}$. Using amplitude amplification achieves a given error using ${\displaystyle O(1/{\sqrt {p}})}$ queries, as opposed to ${\displaystyle 1/p}$ using naive repetition.

After successfully measuring 'well' on ${\displaystyle S}$ the system will be in a state proportional to

${\displaystyle \sum _{i\mathop {=} 1}^{N}\beta _{i}\lambda _{j}^{-1}|u_{j}\rangle =A^{-1}|b\rangle =|x\rangle .}$

The quantum measurement corresponding to M then gives an estimate of ${\displaystyle {\vec {x}}^{T}M{\vec {x}}}$.

Analysis

Classical efficiency

The best classical algorithm which produces the actual solution vector ${\displaystyle {\vec {x}}}$ is Gaussian elimination, which runs in ${\displaystyle O(N^{3})}$ time.

If A is s-sparse and positive semi-definite, then the Conjugate Gradient method can be used to find the solution vector ${\displaystyle {\vec {x}}}$, which can be found in ${\displaystyle O(Ns\kappa )}$ time by minimizing the quadratic function ${\displaystyle \lVert A{\vec {x}}-{\vec {b}}\rVert ^{2}}$.

When only a summary statistic of the solution vector ${\displaystyle {\vec {x}}}$ is needed, as is the case for the quantum algorithm for linear systems of equations, a classical computer can find an estimate of ${\displaystyle {\vec {x}}^{\dagger }M{\vec {x}}}$ in ${\displaystyle O(N{\sqrt {\kappa }})}$.

Quantum efficiency

The runtime of the quantum algorithm for solving systems of linear equations originally proposed by Harrow et al. was shown to be ${\displaystyle O(\kappa ^{2}\log N/\varepsilon )}$, where ${\displaystyle \varepsilon >0}$ is the error parameter and ${\displaystyle \kappa }$ is the condition number of ${\displaystyle A}$. This was subsequently improved to ${\displaystyle O(\kappa \log ^{3}\kappa \log N/\varepsilon ^{3})}$ by Andris Ambainis^[7] and to ${\displaystyle O(\kappa \log N/\varepsilon )}$ for large condition number cases by Peniel Tsemo et al^[8], and a quantum algorithm with runtime polynomial in ${\displaystyle \log(1/\varepsilon )}$ was developed by Childs et al.^[9] Since the HHL algorithm maintains its logarithmic scaling in ${\displaystyle N}$ only for sparse or low rank matrices, Wossnig et al.^[10] extended the HHL algorithm based on a quantum singular value estimation technique and provided a linear system algorithm for dense matrices which runs in ${\displaystyle O({\sqrt {N}}\log N\kappa ^{2})}$ time compared to the ${\displaystyle O(N\log N\kappa ^{2})}$ of the standard HHL algorithm.

Optimality

An important factor in the performance of the matrix inversion algorithm is the condition number ${\displaystyle \kappa }$, which represents the ratio of ${\displaystyle A}$'s largest and smallest eigenvalues. As the condition number increases, the ease with which the solution vector can be found using gradient descent methods such as the conjugate gradient method decreases, as ${\displaystyle A}$ becomes closer to a matrix which cannot be inverted and the solution vector becomes less stable. This algorithm assumes that all singular values of the matrix ${\displaystyle A}$ lie between ${\displaystyle {\frac {1}{\kappa }}}$ and 1, in which case the claimed run-time proportional to ${\displaystyle \kappa ^{2}}$ will be achieved. Therefore, the speedup over classical algorithms is increased further when ${\displaystyle \kappa }$ is a ${\displaystyle \mathrm {poly} (\log(N))}$.^[1]

If the run-time of the algorithm were made poly-logarithmic in ${\displaystyle \kappa }$ then problems solvable on n qubits could be solved in poly(n) time, causing the complexity class BQP to be equal to PSPACE.^[1]

Error analysis

Performing the Hamiltonian simulation, which is the dominant source of error, is done by simulating ${\displaystyle e^{iAt}}$. Assuming that ${\displaystyle A}$ is s-sparse, this can be done with an error bounded by a constant ${\displaystyle \varepsilon }$, which will translate to the additive error achieved in the output state ${\displaystyle |x\rangle }$.

The phase estimation step errs by ${\displaystyle O\left({\frac {1}{t_{0}}}\right)}$ in estimating ${\displaystyle \lambda }$, which translates into a relative error of ${\displaystyle O\left({\frac {1}{\lambda t_{0}}}\right)}$ in ${\displaystyle \lambda ^{-1}}$. If ${\displaystyle \lambda \geq 1/\kappa }$, taking ${\displaystyle t_{0}=O(\kappa \varepsilon )}$ induces a final error of ${\displaystyle \varepsilon }$. This requires that the overall run-time efficiency be increased proportional to ${\displaystyle O\left({\frac {1}{\varepsilon }}\right)}$ to minimize error.

Experimental realization

While a general-purpose quantum computer does not yet exist, one can still try to execute a proof of concept implementation of the HHL algorithm. This remained a challenge for years, until three groups independently did so in 2013.

On February 5, 2013, a group led by Stefanie Barz reported an implementation of the HHL algorithm on a photonic quantum computer. The implementation used two consecutive entangling gates on the same pair of polarization-encoded qubits. Two separately controlled NOT gates were realized where the successful operation of the first was heralded by a measurement of two ancillary photons. Experimental measurements of the fidelity in the obtained output state ranged from 64.7% to 98.1% due to the influence of higher-order emissions from spontaneous parametric down-conversion.^[4]

On February 8, 2013, Pan et al. reported a proof-of-concept experimental demonstration of the quantum algorithm using a 4-qubit NMR quantum computer. The implementation was tested using linear systems of 2 variables. Across three experiments, the solution vector was obtained with over 96% fidelity.^[5]

On February 18, 2013, Cai et al. reported an experimental demonstration solving 2-by-2 linear systems. The quantum circuit was optimized and compiled into a linear optical network with four photonic qubits and four controlled logic gates, which were used to coherently implement the subroutines of the HHL algorithm. For various input vectors, the realization gave solutions with fidelities ranging from 0.825 to 0.993.^[11]

Another experimental demonstration using NMR for solving an 8*8 system was reported by Wen et al.^[12] in 2018 using the algorithm developed by Subaşı et al.^[13]

Proposed applications

Several concrete applications of the HHL algorithm have been proposed, which analyze the algorithm's input assumptions and output guarantees for particular problems.

Electromagnetic scattering

Clader et al. gave a version of the HHL algorithm which allows a preconditioner to be included, which can be used improve the dependence on the condition number. The algorithm was applied to solving for the radar cross-section of a complex shape, which was one of the first examples of an application of the HHL algorithm to a concrete problem.^[14]

Solving linear differential equations

Berry proposed an algorithm for solving linear, time-dependent initial value problems using the HHL algorithm.^[15]

Solving nonlinear differential equations

Two groups proposed^[16] efficient algorithms for numerically integrating dissipative nonlinear ordinary differential equations. Liu et al.^[17] utilized Carleman linearization for second order equations and Lloyd et al.^[18] used a mean field linearization method inspired by the nonlinear Schrödinger equation for general order nonlinearities. The resulting linear equations are solved using quantum algorithms for linear differential equations.

Finite element method

The finite element method approximates linear partial differential equations using large systems of linear equations. Montanaro and Pallister demonstrate that the HHL algorithm can achieve a polynomial quantum speedup for the resulting linear systems. Exponential speedups are not expected for problems in a fixed dimension or for which the solution meets certain smoothness conditions, such as certain high-order problems in many-body dynamics, or some problems in computational finance.^[19]

Least-squares fitting

Wiebe et al. gave a quantum algorithm to determine the quality of a least-squares fit. The optimal coefficients cannot be calculated directly from the output of the quantum algorithm, but the algorithm still outputs the optimal least-squares error.^[20]

Machine learning

Main article: Quantum machine learning

Many quantum machine learning algorithms have been developed, a large number of which use the HHL algorithm as a subroutine. The runtime of certain classical algorithms is often polynomial in the size and dimension of a dataset, while the HHL algorithm can give an exponential speedup in some cases. However, a line work initiated by Ewin Tang has found that for most quantum machine learning algorithms, there are classical algorithms giving the same exponential speedups with similar input assumptions.

Finance

Proposals for using HHL in finance include solving partial differential equations for the Black–Scholes equation and determining portfolio optimization via a Markowitz solution.^[21]

Quantum chemistry

The linearized coupled cluster method in quantum chemistry can be recast as a system of linear equations. In 2023, Baskaran et al. proposed the use of HHL algorithm to solve the resulting linear systems.^[22] The number of state register qubits in the quantum algorithm is the logarithm of the number of excitations, offering an exponential reduction in the number of required qubits when compared to using the variational quantum eigensolver or quantum phase estimation.

Implementation difficulties

Recognizing the importance of the HHL algorithm in the field of quantum machine learning, Scott Aaronson^[23] analyzes the caveats and factors that could limit the actual quantum advantage of the algorithm.

1. the solution vector, ${\displaystyle |b\rangle }$, has to be efficiently prepared in the quantum state. If the vector is not close to uniform, the state preparation is likely to be costly, and if it takes ${\displaystyle O(n^{c})}$ steps the exponential advantage of HHL would vanish.
2. the QPE phases calls for the generation of the unitary ${\displaystyle e^{iAt}}$, and its controlled application. The efficiency of this step depends on the ${\displaystyle A}$ matrix being sparse and 'well conditioned' (low ${\displaystyle \kappa }$). Otherwise, the application of ${\displaystyle e^{iAt}}$ would grow as ${\displaystyle O(n^{c})}$ and once again, the algorithm's quantum advantage would vanish.
3. lastly, the vector, ${\displaystyle |x\rangle }$, is not readily accessible. The HHL algorithm enables learning a 'summary' of the vector, namely the result of measuring the expectation of an operator ${\displaystyle \langle x|M|x\rangle }$. If actual values of ${\displaystyle {\vec {x}}}$ are needed, then HHL would need to be repeated ${\displaystyle O(n)}$ times, killing the exponential speed-up. However, three ways of avoiding getting the actual values have been proposed: first, if only some properties of the solution are needed;^[24] second, if the results are needed only to feed downstream matrix operations; third, if only a sample of the solution is needed.^[25]

See also

• Differentiable programming