HHL: Solving Linear Systems of Equations with Quantum Computing
Ronald Vaughn II and Julian Viera
Jan 29, 2021
Introduction
We’ve learned in 6.s089 that quantum computing can potentially provide substantial (i.e. exponential) speedup to many classical algorithms we use today. In class, we’ve discussed algorithms such as Deutsch-Jozsa, Grover’s Search, and Shor’s quantum factoring as prime examples of this speedup. But we’ve left out an important one: HHL. (Fun fact: the first H in HHL is Harrow, one of the MIT professors at our class panel on Wednesday!)
The HHL algorithm, put simply, solves a linear system of equations. This problem is ubiquitous throughout science and engineering and has been an extremely active area of research for decades. In every sector of research, manipulating matrices efficiently has become a massive area of interest. Unlike the classical solutions to the Deutsch-Jozsa and search problems, most of our classical methods for matrix manipulation do work in polynomial time. However, as data analysis becomes more and more powerful (and more and more demanding on today’s computers), the size of these matrices can make even polynomial time too long.
One of the simplest and yet most important problems in linear algebra is solving a system of equations, or calculating the solution to Ax = b. Currently, for reasonably sparse N x N matrices, the best classical algorithms we have take ~O(N) time. Again, this is clearly not as undesirable as exponential time, but we can do better.
Introducing HHL. Developed by Aram Harrow, Avinatan Hassidim, and Seth Lloyd, this is a quantum algorithm that approximates a function of the solution vector of a linear system of equations in ~O(log N) time. This time reduction has significant implications for the speedup of many machine learning algorithms on quantum computers.
Now let’s see how it works.
First, Some Linear Algebra
First, let’s discuss some preliminary linear algebra. We are considering the equations Ax = b where A is Hermitian, which means that we can rewrite A as
In words, A can be written as a the sum of the outer products of its eigenvectors, scaled by its eigenvalues. Therefore, we can write the inverse of A as
Since A is invertible and Hermitian, it must have an orthogonal basis of eigenvectors, and thus we can write b as
The main goal of the HHL algorithm is to solve this system of equations, or get a register in the state
Some Conditions
There are some constraints on the inputs to the HHL algorithm that need mentioning. The matrix A does not have to be Hermitian, but in the case that it is not Hermitian we must transform it into a matrix that is. But the matrix A must be sparse and well-conditioned. A matrix is s-sparse when it has at most s non-zero entries in any row or column, and a matrix is simply sparse if it has at most poly(log N) entries in any row or column. A matrix is well-conditioned when its singular values lie between the reciprocal of its condition number and 1. The condition number 𝜅 of a Hermitian matrix is the ratio of the largest to smallest eigenvalue magnitude and is undefined when the smallest eigenvalue magnitude is 0.
The Algorithm
At a very high level, the algorithm has three big steps to do this: phase estimation, eigenvalue inversion rotation, and inverse phase estimation.
We start with three registers: an ancilla register S, a clock register C, and an input register I. The procedure for the algorithm is as follows:
- Apply Quantum Phase Estimation (QPE) with
Now, the register, expressed in the eigenbasis of A, is now
where here nl represents the number of bits for the binary representation of each eigenvalue and nb is log(N), the number of qubits representing |b⟩.
2. Apply a rotation of the ancilla qubit to get a normalized state of the form
These rotations are conditioned on each eigenvalue of A, and they can be achieved by the application of the Y rotation matrix where
Note also that C here is a scaling factor we use to make sure that all of our quantum states are normalized.
3. The final step of the algorithm is to perform inverse quantum phase estimation. This sets the clock register back to 0’s and leaves the remaining state as
Notice that if the ancilla qubit is 1, then the state of the registers is a normalized equivalent of the solution vector. So we measure this ancilla qubit, and if we do get 1, the new state indeed is
which is component-wise proportional to the solution vector |x⟩.
Runtime Analysis
As we can see, phase estimation is one the main components of the HHL algorithm. And using the runtime we learned in class for QPE, we can see that the runtime for HHL is ~O(log N). However, we must take the conditions on A into mind. In fact, the sparsity and condition number of A do affect the runtime, such that for an s-sparse matrix with condition number 𝜅 the runtime is actually O(log(N)s²𝜅²/ε), where ε is the error parameter.
The best existing classical matrix inversion algorithm involves Gaussian elimination which takes O(N³)time, and for matrices s-sparse with condition number 𝜅, the Conjugate Gradient algorithm can be used to find the solution vector in O(Ns𝜅 log(1/ε)). So clearly, HHL provides an exponential speedup in N. However, it is important to note why that is. Actually reading out the solution vector would take O(N)time, so we can only maintain the logarithmic runtime by sampling the solution vector like ⟨x|M|x⟩, where M is a quantum-mechanical operator. Therefore, HHL is useful mainly in applications where only samples from the solution vector are needed. Additionally, although HHL is exponentially faster than Conjugate Gradient in N, it is polynomially slower in s and 𝜅, so HHL is restricted to only those matrices that are sparse and have low condition numbers.
Our Circuit
Let’s take a closer look at the components of the circuit’s sections previously described above. Below is a circuit that we’ve designed using IBM’s quantum circuit simulator. Here, the ancilla register is represented by qubit q0, the two qubits of the clock register are q1 and q2, and the input register consists only of q3. This 4 qubit circuit is enough to calculate systems of 2 x 2 matrices that follow the restrictions mentioned earlier.
As usual, you can expect all qubits to initialize at |0⟩.
- Before the first barrier, we use the Y arbitrary rotation gate to put q3 in the |b⟩ state. To do this, first normalize the vector b, then alter the probability amplitudes of q3 to match the entries of b by using the Y arbitrary rotation gate. If you have more than one qubit in the input register, the entries of b should match the probability amplitudes of the states of the input register.
- Next, we begin our phase estimation step. The first part of this is to apply hadamard gates on all qubits in the clock register. Then, we apply the conditional unitary gates U on our input register with the control being on the clock register’s qubits. The t in the U operation is usually approximated as being the number of qubits in the clock register (for our case, t=2). Note: The two U gates in the image above are different from each other such that the first U gate actually represents U2, while the second just represents U.
- Moving from 2 to 3, we apply an inverse Quantum Fourier Transform, which is also the last step in the quantum phase estimation.
4. Next, we need to apply what the first image referred to as the “R rotation”. on the ancilla qubit by using multiple conditional Y arbitrary rotation gates. Let r be the number of qubits in the clock register. The first rotation will be conditional on the first clock qubit and will be by π/2r. The next will be conditional on the next clock qubit and will be by π/2r-1. You continue with this pattern until you have a conditional rotation gate for every qubit in the clock register. For the case above it is simple since there are only 2 qubits in the clock register.
5. We now begin the “uncompute” section which will simply be the inverse quantum phase estimation. By barrier 5, we will have applied a Quantum Fourier Transform using all of the clock register and input register qubits.
6. The rest of the inverse quantum phase estimation is taken care of by applying the inverse conditional U gates, then applying hadamard one last time on all clock register qubits. Here, the first U gate will be the inverse of U and the second will be the inverse of U squared.
7. The final step is to take a measurement of the ancilla qubit. This will make it so that the probability amplitudes of the states of our input registry will correspond to the x vector’s entries.
Other Applications
While we’ve only demonstrated here how HHL works and a simple 2 x 2 example, there are many potential applications for the HHL algorithm. Below we list a few areas where there is active research into using the HHL.
Further study has shown that a preconditioner (a transformation that reduces the condition number of a matrix) can be included in HHL. Since HHL is polynomial in the condition number of the matrix, this expanded the class of problems for which using HHL could lead to a potential speedup. This advance was then used to show an algorithm using HHL that determines the radar cross-section (how visible something is to radar) of a complex shape.
An extension of HHL has been developed that efficiently solves the full-time evolution under sparse linear differential equations on a quantum computer. This clearly has many useful applications for a variety of physical systems.
The time required for a least-squares fit (a continuous function to model a set of discrete points) has been found to become very large as the amount of discrete points increases. However, a research group has developed a quantum algorithm that incorporates HHL to solve many cases of least-squares efficiently, eliminating the need for the previously best higher-complexity algorithm.
Of course, HHL is thought to be massively useful in machine learning, where classical algorithms are limited by a polynomial dependence on the volume of data and dimensions of the space. One specific use for HHL in machine learning is in an optimized linear or non-linear binary classifier, or support vector machine. Research has shown that a quantum support vector machine can be used for big data classification and achieve an exponential speedup over classical computers. A quantum algorithm has also been developed for Bayesian training of deep neural networks, with an exponential speedup over classical training due to the use of HHL.
Conclusion
In this blog post, we have provided a general overview of the HHL quantum algorithm for solving a linear system of equations. We’ve gone through the actual procedure of the algorithm, discussing the steps: quantum phase estimation, eigenvalue rotation, and inverse quantum phase estimation. We’ve also analyzed the runtime and some of the limitations of the algorithm. While it provides an exponential speedup in N, it does not provide any speedup for dense matrices, or those with large condition numbers. We also rely on the assumption that the b vector can be constructed efficiently. However, the preparation of arbitrary quantum states is hard and subject to several constraints.
We then provided an example we constructed of a circuit made in IBM Quantum Experience that solves a specific system of two equations, where you can see what each step of the algorithm actually looks like and does.
Lastly, we discussed some other active areas of research where the HHL algorithm has made a significant impact. Clearly, if and when quantum computers become widely accessible and scalable, just as Shor’s algorithm will have tremendous impact on cryptography and security, the HHL algorithm will have a massive impact in just about every area using data analysis.
