A Step-by-Step HHL Algorithm Walkthrough to Enhance Understanding of Critical Quantum Computing Concepts

For readers who have a fresh memory of the basic concepts, they can start reading from Section II, in which the HHL algorithm is discussed step-by-step analytically followed by a numerical example in Section III. The basic concepts mentioned in the Appendix are referred to in the main text and readers are encouraged to review them when needed.

For readers who need reviews on the basic concepts first, they are encouraged to go over the Appendix first before reading the main text.

For readers who have devoted substantial time to learning HHL elsewhere but just need a numerical example to reinforce the understanding, they might start with the numerical example in Section III.

Equations in the Appendix begin with ’V’. If the equations are examples from the main text, the same equation number is used in the Appendix.

We will first give an overview of the problem and the HHL algorithm. Details will be discussed in the following subsections with reference to the Appendix for reviewing basic concepts. A linear system problem (LSP) can be represented as the following

where $A$ is a $N_{b}\times N_{b}$ Hermitian matrix and $\vec{x}$ and $\vec{b}$ are $N_{b}$-dimensional vectors. For simplicity, it is assumed $N_{b}=2^{n_{b}}$, where $n_{b}$ is the number of qubits in the quantum circuit, and $N_{b}$ is the total combinations due number of qubits,$n_{b}$. In matrix representation, qubits are represented by their total combinations($N_{b}\times N_{b}$), or we can say for $N_{b}$ number of unknowns, we need $n_{b}$ qubits to solve unknowns. Dummy equations can be added otherwise to convert the system satisfy this assumption. $A$ and $\vec{b}$ are known and $\vec{x}$ is the unknown to be solved, i.e.

As an example, $A=\begin{pmatrix}1&-\frac{1}{3}\\ -\frac{1}{3}&1\\ \end{pmatrix}$, $\vec{b}=\begin{pmatrix}0\\ 1\\ \end{pmatrix}$, and $\vec{x}=\begin{pmatrix}\frac{3}{8}\\ \frac{9}{8}\\ \end{pmatrix}$ with $n_{b}=1$ and $N_{b}=2^{1}=2$. Readers may refer to Appendix V-13 and Appendix V-14 to review how LSP is solved classically using Gaussian Elimination and Conjugate Gradient Method, respectively.

$A$ is assumed to be Hermitian (See Appendix V-1). If it is not Hermitian, the $A$ can be converted to $\begin{pmatrix}0&A\\ A^{\dagger}&0\\ \end{pmatrix}$, which is Hermitian. Readers may refer to [20] for the more advanced treatment when $A$ is not Hermitian.

Figure 1 shows the schematic of the $HHL$ algorithm and the corresponding circuit to solve LSP. In the HHL quantum algorithm, the $N_{b}$ components of $\vec{b}$ and $\vec{x}$ are encoded as the amplitudes/coefficients (amplitude encoding) of basis states of the $n_{b}$-qubits, $\ket{\;}_{b}$, which form a $\mathbb{C}^{N_{b}}$ Hilbert space. These $n_{b}$ qubits are called b-register. Qubit $n_{b}$ is chosen to be large enough to encode $\vec{b}$, i.e. $2^{n_{b}}$ needs to be the same as the length of the vectors $\vec{b}$ and $\vec{x}$. The matrix $A$ is simulated through Hamiltonian encoding by encoding it as the Hamiltonian of a unitary gate. Appendix V-9 reviews examples of various encoding schemes.

The HHL algorithm has 5 main components, namely state preparation, quantum phase estimation (QPE), ancilla bit rotation, inverse quantum phase estimation (IQPE), and measurement. In this paper, the little-endian convention is used. In a little-endian convention, the rightmost (ending) qubit represents the least significant bit (LSB). For example, in a 4-qubit system, $\ket{0001}$ in binary is $\ket{1}$ in decimal because the $1$ in the basis state $\ket{0001}$ is the LSB, representing $2^{0}$ instead of $2^{3}$ (if it were the most significant bit, MSB). Moreover, in the circuit diagrams, the lowest qubit represents the MSB and the topmost qubit represents the LSB, which is a convention used in qiskit [16] and the IBM-Q platform [17]

As shown in Figure 1, besides the b-register, which belongs to the more significant bits, there are two more sets of inputs to the algorithm. The first set is sometimes called the c-register because it is related to the time (clock) in the controlled rotation in the QPE part. Therefore, they are also called the clock qubits.The c-register stores the values of the phase of the eigenvalues of the $A$ matrix after the QPE. There are $n$ qubits in the c-register. Since basis encoding is used (i.e. the phase value is encoded as the basis number (See Appendix V-9), the value of $n$ determines how accurately the phase can be stored. A larger $n$ results in higher accuracy when the encoding is not exact. We set $N=2^{n}$.

The last set of qubits is the ancilla qubit $\ket{\;}_{a}$ which is the LSB. The ancilla qubit, as its name implies, is important to help achieve the goal although it will be discarded at the end, as will be detailed later.

The matrix $A$, which is a Hamiltonian, may be written as a linear combination of the outer products of its eigenvectors, $\ket{u_{i}}\bra{u_{i}}$ weighted by its eigenvalues, $\lambda_{i}$,in Eq.(3), (See Appendix V-8).

Since A is diagonal in its eigenvector basis, its inverse is simply, $A^{-1}=\sum_{i=0}^{2^{n_{b}}-1}\lambda_{i}^{-1}\ket{u_{i}}\bra{u_{i}}$. $\vec{b}$ can be also expressed in the basis formed by the eigenvectors of A, such that

Therefore, Eq. (2) can be encoded as,

The equation needs also to be prepared so that the eigenvectors, $\ket{u_{i}}$, and $\ket{b}$ are normalized so they can be properly represented as a unit vector in quantum computing. Therefore, (4) and Eq. (5) require

There are total $n_{b}+n+1$ qubits and they are initialized as

In the state preparation, $\ket{0\cdots 0}_{b}$ in the b-register needs to be rotated to have the amplitudes correspond to the coefficients of $\vec{b}$. That is

The vector $\vec{b}$ is represented in a column form on the left with coefficients $\beta^{\prime}s$. This is also a valid representation of $\ket{b}$. On the right, the corresponding basis of the Hilbert space formed by the $n_{b}$ qubits is written explicitly. Therefore,

From now on, some of the subscripts of the kets will be omitted when there is no ambiguity. Since the state preparation depends on the actual value of $\vec{b}$, it will be discussed in more detail in the numerical example.

Quantum phase estimation (QPE) is also an eigenvalue estimation algorithm. QPE has three components, namely the superposition of the clock qubits through Hadamard gates, controlled rotation, and inverse quantum Fourier transform (IQFT). The goal of QPE is to estimate the phase of the eigenvalues of the unitary rotation matrix, $U=e^{iAt}$, in the controlled gate, $C-U$, (Fig. 1) used in the QPE. Again, this gate encodes the matrix $A$ as its Hamiltonian. It is also instructive to note that the eigenvalues of $U$ must be roots of unity (i.e. in the form of $e^{i\theta}$) as $U$ is unitary. Therefore, the phase of the eigenvalue of the gate is proportional to the eigenvalue of $A$. As a result, by using QPE in the HHL algorithm, it is expected the eigenvalues of $A$ will be encoded in the c-register after the QPE, i.e. at $\ket{\Psi_{4}}$. As it will be clear later, the eigenvalues are only encoded through basis encoding. The c-register does not store the exact eigenvalues.

Here we assume the readers are already familiar with IQFT and it will not be explained in detail. Readers may review the basic concepts in Appendices V-11 and V-12.

In the first step of QPE, Hadamard gates are applied to the clock qubits to create a superposition of the clock qubits,

In the controlled rotation part, controlled gates are applied to $\ket{b}$ with the clock qubits as the controlling qubits (Figure 2). The number of the qubit, $n$, of the c-register determines the number of the controlled gates. The gates are in the form of $U^{2^{r}}$, where $r$ is the index of the clock qubit. Also, $U=e^{iAt}$. For the most significant bit in the c-register, $\ket{c_{n-1}}$, it controls the gate $U^{2^{n-1}}$ on the b-register while the least significant one, $\ket{c_{0}}$, controls the gate $U^{2^{0}}=U$ on the b-register.

To understand how it works, we begin by assuming that $\ket{b}$ is an eigenvector of $U$ with eigenvalue $e^{2\pi i\phi}$. The eigenvalue is written in this form so that the phase, $\phi$, will be encoded as the basis state in the c-register (See Eq. (II-C)). Therefore, based on the definition of eigenvalues and eigenvectors (See Section V-8),

When the control clock qubit is $\ket{0}$, $\ket{b}$ will not be affected. If the clock bit is $\ket{1}$, $U$ will be applied to $\ket{b}$. This is equivalent to multiplying $e^{2\pi i\phi 2^{j}}$ in front of the $\ket{1}$ of the $j$th clock qubit, $\ket{c_{j}}$, as one can assign the prefactor to the controlling qubit. Therefore, after the controlled-$U$ operation, we have

In the IQFT part,(II-C), only the clock qubits are affected. Note that in certain literature, this is called Quantum Fourier Transform (QFT) (Appendix V-11).

Due to interference, only $\ket{y}$ satisfying the condition $\phi-y/N=0$ will have a finite amplitude of $\sum_{k=0}^{2^{n}-1}e^{0}=2^{n}$. Otherwise, the amplitude is $\sum_{k=0}^{2^{n}-1}e^{2\pi ik(\phi-y/N)}=0$ due to destructive interference. By ignoring the states with zero amplitude, we may rewrite $\ket{\Psi_{4}}$ as

Therefore, in QPE, the clock qubits are used to represent the phase information of $U$, which is $\phi$, and the accuracy depends on the number of qubits, $n$.

Since in Hamiltonian encoding, $U$ is related to $A$ through

where $t$ is the evolution time for that Hamiltonian. $U$ is also diagonal in $A^{\prime}s$ eigenvector, $\ket{u_{i}}$, basis. If $\ket{b}=\ket{u_{j}}$,

By equating $i\lambda_{j}t$ to $2\pi i\phi$ in Eq. (13), we get $\phi=\lambda_{j}t/2\pi$ and Eq. (II-C) becomes

Thus the eigenvalues of $A$ have been encodeded in the clock qubits (basis encoding). However, in general, given in (4), by superposition,

The $\lambda_{j}$ are usually not integers. We will choose $t$ so that $\tilde{\lambda_{j}}=N\lambda_{j}t/2\pi$ are integers. Therefore, $\tilde{\lambda_{j}}$ are usually scaled version of $\lambda_{j}$.

$\Psi_{4}$ can be rewritten as

The next step is to rotate the ancilla qubit, $\ket{0}_{a}$, based on the encoded eigenvalues in the c-register, such that,

where $C$ is a constant. The goal is to show why this is useful.

When the ancilla qubit is measured, the ancilla qubit wavefunction will collapse to either $\ket{0}$ or $\ket{1}$. If it is $\ket{0}$, the result will be discarded and the computation will be repeated until the measurement is $\ket{1}$. Therefore, the final wavefunction of interest is

where the prefactor is due to normalization after measurement. Since $\lvert\frac{C}{\tilde{\lambda_{j}}}\rvert^{2}$ is the probabily of obtaining $\ket{1}$ when the ancilla bit is measured, $C$ should be chosen to be as large as possible. Compared to (5), the result resembles the answer $\ket{x}$ that we are looking for. However, we can only obtain the correct result if the b-register is measured in the eigenvector basis (i.e. $\ket{u_{j}}$ instead of $\ket{0}/\ket{1}$). However, the b-register is entangled with the clock qubits, $\ket{\tilde{\lambda_{j}}}$. This means that we cannot factorize the result into a tensor product of the c-register and b-register (See the discussion after (5) and Appendix V-6). As a result, we cannot convert the b-register into the $\ket{0}/\ket{1}$ measurement basis with the desired amplitudes. We will need to uncompute the state so that it gives the right results in the $\ket{0}/\ket{1}$ measurement during which the b-register and c-register will be unentangled.

The measurement of the ancilla qubit can be and is usually performed after uncomputation. However, since the ancilla bit is not involved in any operations after the controlled rotation, measuring the ancilla bit before the uncomputation gives the same result. For simplicity in the derivation, it is thus performed before the uncomputation.

The uncomputation is done by using inverse QPE. Firstly, QFT is applied to the clock qubits as,

Then inverse controlled-rotations of the b-register by the clock qubits are applied with $U^{-1}=e^{-iAt}$. Similar to the forward process, when the controlling $r$th clock qubit is $\ket{0}$, $\ket{u_{j}}$ will not be affected. If the $r$th clock qubit is $\ket{1}$, $(U^{-1})^{2^{r}}$ will be applied to $\ket{u_{j}}$. This is equivalent to multiplying $e^{-i\lambda_{j}ty}$ if the c-register is $\ket{y}$. This is due to the similar argument in Eq. (II-C) and the fact that $2\pi i\phi=i\lambda_{j}t$. Therefore,

Since we earlier chose to set $\tilde{\lambda_{j}}=N\lambda_{j}t/2\pi$, therefore, the two exponential terms cancel each other and

The clock qubits and the b-register are now unentangled and the b-register stores $\ket{x}$. By applying the Hadamard gate on the clock qubits, finally, we have

If $C$ is real and by using (7),

Here, the solution $\ket{x}$ (Eq. 5) is stored in the b-register successfully.

In this example, the matrix $A$ and vector $\vec{b}$ are set to be

The eigenvectors of $A$ are $\vec{u_{0}}=\begin{pmatrix}\frac{-1}{\sqrt{2}}\\ \frac{-1}{\sqrt{2}}\end{pmatrix}$, $\vec{u_{1}}=\begin{pmatrix}\frac{-1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}\end{pmatrix}$ with eigenvalues $\lambda_{0}=\frac{2}{3}$ and $\lambda_{1}=\frac{4}{3}$ respectively. We need to using basis encoding to encode the eigenvalues in the basis formed by the clock qubit and 2 qubits are needed by encoding $\lambda_{0}$ as $\ket{01}$ and $\lambda_{1}$ as $\ket{10}$ so that it maintains the ratio of $\lambda_{1}/\lambda_{0}=2$. This means $\tilde{\lambda_{0}}=1$ and $\tilde{\lambda_{1}}=2$ or in other words, $\tilde{\ket{\lambda_{0}}}=\ket{01}$ and $\tilde{\ket{\lambda_{1}}}=\ket{10}$. This gives a perfect encoding with $n=2$ (i.e. $N=4$). Therefore, $t$ is chosen to be $\frac{3\pi}{4}$ to achieve the encoding scheme since $\tilde{\lambda_{j}}=N\lambda_{j}t/2\pi$.

Since $\vec{b}$ is a 2-dimensional complex vector, it can be encoded using 1 qubit and, thus, $n_{b}=1$.

The solution to the LSP is found to be

whereby, the ratio of $\lvert x_{0}\rvert^{2}$ to $\lvert x_{1}\rvert^{2}$ is $1:9$.

In reality, we expect the controlled-$U$ operation to be implemented by Hamiltonian simulation [18]. However, to understand the algorithm, we will derive the matrix for $U$ and then map this to the $Controlled-U$ gate used in IBM-Q directly. Since $n=2$, there are two operations needed, namely $U^{2^{1}}=U^{2}$ and $U^{2^{0}}=U$, controlled by $c_{1}$ and $c_{0}$, respectively.

In order to find the corresponding matrix for $U^{2}=e^{i2At}$ and $U=e^{iAt}$, we need to perform similarity transformation on $i2At$ and $iAt$, exponentiate then, and transforms back to the original basis.

The transformation matrix from the original basis to the eigenvector basis is

Since $V$ is real and symmetric,its conjugate, $V^{\dagger}$ equals itself.

The diagonalized $A$, i.e. expressed in the basis formed by $\vec{u_{0}}$ and $\vec{u_{0}}$, is

As it is diagonal, $U$ can be obtained by exponentiation of the elements accordingly.

where $t=3\pi/4$ as mentioned earlier is used. Also,

It is worth noting that both are naturally unitary which is a requirement for a quantum operation.

To obtain $U$ and $U^{2}$ in the original basis, we need to apply similarity transformation again in the reverse direction,

To implement $U$ and $U^{2}$, a 4-parameter arbitrary unitary gate with global phase [19],

By choosing $\theta=\pi,\phi=\pi,\lambda=0,\gamma=0$, $U^{2}$ is implemented.

By choosing $\theta=\pi/2,\phi=-\pi/2,\lambda=\pi/2,\gamma=3\pi/4$, $U$ is implemented.

For the IQPE part, we also need to implement $U^{-1}$ and $U^{-2}$. Since in this example, ${(U^{2})}^{-1}=U^{2}$, one can use the same set of parameters to implement $(U^{2})^{-1}$.

However,

We need to choose $\theta=\pi/2,\phi=\pi/2,\lambda=-\pi/2,\gamma=-3\pi/4$ to implement $U^{-1}$.

The controlled version of matrix $U^{\prime}$ can then be constructed using

Note that in this equation, only the controlling clock bit and the b-register are included for simplicity. For example,

The coefficients of $\ket{0}$ and $\ket{1}$ of the ancilla bit after rotation in Eq. (22) are $\sqrt{1-\frac{C^{2}}{\tilde{\lambda_{j}}^{2}}}$ and $\frac{C}{\tilde{\lambda_{j}}}$, respectively. The sum of the square of the magnitudes of the coefficients is $1$ as required. This means also $C\leq\tilde{\lambda_{j}}$. Since the minimal $\tilde{\lambda_{j}}$ is $1$, we will set $C=1$ to maximize the probability of measuring $\ket{1}$ during the ancilla bit measurement.

The transformation of $\ket{0}_{a}$ to $\sqrt{1-\frac{1}{\tilde{\lambda_{j}}^{2}}}\ket{0}_{a}+\frac{1}{\tilde{\lambda_{j}}}\ket{1}_{a}$ is known to be equivalent to $RY(\theta)$ rotation,

with $\theta=2\arcsin(\frac{1}{\tilde{\lambda_{j}}})$. One can check this by multiplying $RY(\theta)$ to $\ket{0}_{a}$.

Therefore, we will establish a function to implement this rotation and this function only need to be valid when the input are the encoded eigenvalues because only encoded eigenvalues have zero magnitudes in the c-register as shown in (21). The function is defined as

where $c$ is the value of the clock qubits and $c_{1}c_{0}$ is its binary form.

Since only $\ket{\tilde{\lambda_{j}}}$ has non-zero amplitude in (21), we only need to set up (44) such that it is correct for $\ket{c}$=$\ket{01}$ and $\ket{10}$, namely

The following function can achieve the goal,

Therefore, the controlled rotation can be implemented as

where the operators operate on qubits $\ket{c_{1}}$, $\ket{c_{0}}$, and $\ket{a}$ from left to right, respectively.

An HHL circuit for the numerical example is then built and shown in Figure 3. We will then walk through the circuit using numerical substitution.

The algorithm begins with

X-gate is then applied to convert $\ket{0}_{b}$ to $\ket{1}_{b}$ with

After applying the Hadamard gates to create a superposition among the clock qubits,

Before applying the CU3(controlled rotation of ancillary qubit) gates in the bra-ket notation, it will be convenient to perform a basis change to the eigenvector basis of $A$. Since $\ket{1}=\frac{1}{\sqrt{2}}(-\ket{u_{0}}+\ket{u_{1}})$, we have $b_{0}=\frac{-1}{\sqrt{2}}$ and $b_{1}=\frac{1}{\sqrt{2}}$. Therefore,

In the controlled rotation operations, when the corresponding c-register is $\ket{k}_{c}$, a phase change of $\phi_{j}=k\lambda_{j}t/2\pi$ is added (i.e. multiplied by $e^{2\pi i\phi_{j}}$) for $\ket{u_{j}}$. Since $t=\frac{3\pi}{4}$, $\lambda_{0}=\frac{2}{3}$ and $\lambda_{1}=\frac{4}{3}$, we have