HHL algorithm

Harrow-Hassidim-Lloyd (HHL) quantum algorithm [1] can be used to solve linear system problems and can provide exponential speedup over the classical conjugate gradient method chosen for this purpose. HHL is the basis of many more advanced algorithms and is important in various applications such as machine learning and modeling of quantum systems.

A linear system problem (LSP) can be represented as the following: Given a matrix $A$ and a vector $b$, find the solution $x$ that satisfies the linear system $$A\vec{x}=\vec{b}.$$

In the quantum case, we consider a special form here for simplicity. We assume $A$ is a Hermitian positive-definite matrix in $\mathbb{C}^{N \times N}$, and $\vec{b} \in \mathbb{C}^{N}$ can be encoded as a quantum state $|b\rangle$. The solution $x$ will also be readout as a quantum state $|x\rangle$. We also assume $N$ be a power of $2$.

In the classical case, the solution now can be represented as $$\vec{x}=A^{−1}\vec{b}.$$ Similarly, in the quantum case, the solution can be represented as $$|x\rangle \propto A^{-1} |b\rangle.$$

Assume the spectral decomposition of $A$ is $$A=\sum_{j=1}^{N}\lambda_j|u_j\rangle\langle u_j|,$$ where $\{\lambda_j\}$ and $\{|u_j\rangle\}$ are the eigenvalues and the eigenvectors of $A$, respectively. Since $\{|u_j\rangle\}$ are orthogonal, $|b\rangle$ can also be decomposed onto this basis: $$|b\rangle=\sum_{j=1}^{N}\beta_j|u_j\rangle.$$ Therefore, $$|x\rangle \propto A^{-1}|b\rangle=\left(\sum_{j=1}^{N}\lambda_j^{-1}|u_j\rangle\langle u_j|\right)\sum_{j=1}^{N}\beta_j|u_j\rangle=\sum_{j=1}^{N}\lambda_j^{-1}\beta_j|u_j\rangle.$$

The idea of the HHL algorithm is to construct a mapping $|u_j\rangle \mapsto \lambda_j^{-1} |u_j\rangle$ for each eigenvector $|u_j\rangle$. Therefore, when the input state is $|b\rangle$, the algorithm will scale each eigenvector $|u_j\rangle$ by a multiplicative factor $\lambda_j^{-1}$, thereby preparing the solution state $|x\rangle$.

The HHL algorithm consists of five major steps, as shown as follows:

• Step 1: State Preparation. We need a unitary $U$ that prepares $|b\rangle$ in the input register from $|0\rangle$.

• Step 2: Phase Estimation. Since $A$ is a Hermitian positive-definite matrix, let us define the unitary $U=e^{iAt}$. Thus $$U=\sum_{j=1}^{N}e^{i\lambda_jt}|u_j\rangle\langle u_j|,$$ i.e., the eigenvalue corresponding to eigenvector $|u_j\rangle$ is $e^{i\lambda_jt}$. Recall that the quantum phase estimation (QPE) algorithm obtains $N\theta$ approximately for eigenvalue $e^{2\pi i\theta}$. Therefore, applying QPE with $U$ will encode $\lambda_j$ in computational basis in the ideal case. The qubits that hold $\lambda_j$ are usually named the clock register. If the clock register contains $m$ qubits, thus it can contain integers from $1$ to $T-1$, where $T=2^m$. Let $U = e^{i A t_0/T}$, where $t_0$ is a constant to be chosen. Thus QPE will yield: $$|b\rangle |0\rangle = \sum_{j=1}^{N} \beta_j |u_j\rangle |0\rangle \mapsto \sum_{j=1}^{N} \beta_j |u_j\rangle |k_j\rangle$$ in the ideal case, where $k_j = \lambda_j t_0/2\pi$ is a integer from $1$ to $T-1$. By the way, the unitary $U=e^{iAt}$ can be implemented via Hamiltonian simulation. A simple way using Trotter decomposition is introduced in [2, Section 4.7].

• Step 3: Controlled Rotation. This step consists of the following sequence of gates which is controlled by an ancilla initialized to $|0\rangle$: $$\sum_{\tau=1}^{T-1} |\tau\rangle \langle \tau| \otimes \mathrm{RY} (2 \arcsin \frac{C}{\tau}),$$ where $C$ is a constant to be chosen, and $\lambda_{\tau} := 2 \pi \tau / t_0$ is the eigenvalue corresponding to the integer $\tau$. Since $$\mathrm{RY} (2 \arcsin \frac{C}{\tau}) |0\rangle = \sqrt{1 - \frac{C^2}{\lambda_{\tau}^2}} |0\rangle + \frac{C}{\lambda_{\tau}} |1\rangle,$$ this yields: $$\sum_{j=1}^{N} \beta_j |u_j\rangle |k_j\rangle |0\rangle \mapsto \sum_{j=1}^{N} \beta_j |u_j\rangle |k_j\rangle \left( \sqrt{1 - \frac{C^2}{\lambda_{j}^2}} |0\rangle + \frac{C}{\lambda_{j}} |1\rangle \right).$$ This state contains all ingredient of $|x\rangle$. However, the clock register and the input register is entangled, thus we cannot obtain $|x\rangle$ if we measure the ancilla register now.

• Step 4: Inverse Phase Estimation. We need to uncompute Phase Estimation to change the clock register back to $|0\rangle$. Therefore, this step is precisely the inverse of Step 2.

• Step 5: Repeat the circuit until the measurement result of the ancilla qubit is 1. In this case, the final state is $$\propto \sum_{j=1}^{N} \frac{\beta_j}{\lambda_j}|u_j\rangle |0\rangle |1\rangle \propto |x\rangle |0\rangle |1\rangle$$ in the idea case. Thus, the algorithm achieves the goal.

Let us write an isQ program to demonstrate the HHL algorithm with the following example of the linear system: $$\frac{1}{6} \begin{pmatrix} 3 & -1\\ -1 & 3 \end{pmatrix}\vec{x}=\begin{pmatrix} 0 \\ 1 \end{pmatrix}.$$

We solve this linear system theoretically first. The spectral decomposition of $A$ is $$A = \lambda_1 |u_1\rangle \langle u_1 | + \lambda_2 |u_2\rangle \langle u_2 | = \frac{1}{3} |+\rangle \langle +| + \frac{2}{3} |-\rangle \langle -|.$$ Thus: $$|b\rangle = \beta_1 |u_1\rangle + \beta_2 |u_2\rangle = \frac{1}{\sqrt{2}} |+\rangle - \frac{1}{\sqrt{2}} |-\rangle.$$ $$A^{-1} |b\rangle = \frac{\beta_1}{\lambda_1} |u_1\rangle + \frac{\beta_2}{\lambda_2} |u_2\rangle = \frac{3}{\sqrt{2}} |+\rangle - \frac{3}{2\sqrt{2}} |-\rangle = \frac{3}{4} \begin{pmatrix} 1 \\ 3 \end{pmatrix}.$$ $$|x\rangle = \frac{A^{-1} |b\rangle}{\| A^{-1} |b\rangle \|} = \frac{2}{\sqrt{5}} |+\rangle - \frac{1}{\sqrt{5}} |-\rangle = \frac{1}{\sqrt{10}} \begin{pmatrix} 1 \\ 3 \end{pmatrix}.$$

To solve this linear system via HHL, we need to bound all the eigenvalues of $A$ in $[1/\kappa, 1]$ first, where $\kappa \ge 1$ is a known parameter. If we know the bound, we can scale $A$ to satisfy the above requirement, thus the resulting $\kappa$ can be seen as an upper bound of the condition number of $A$.

For the above example, we can define $\kappa = 3$. We can choose $t_0 = 2\pi \kappa = 6 \pi$ and $T = \kappa + 1 = 4$, thus the eigenvalues $\lambda = 1/3, 2/3, 3/3$ will be converted to the integers $k=1,2,3$, respectively. The constant $C$ is related to the success probability of the algorithm. The probability of the measurement outcome being $|1\rangle$ is $$\Pr[1] = C^2 \sum_{j=1}^{N} \left| \frac{\beta_j}{\lambda_j} \right|^2 \ge C^2,$$ where $C$ can be chosen as a positive real number without loss of generality. On the other hand, since the quantum state must be a unit vector, $C$ must satisfy $$C \le \min_{\tau = 1,\dots,T-1} \lambda_{\tau} = \frac{2\pi}{t_0}$$ . Thus we can choose $C = 2\pi / t_0 = 1/3$.

Now we have determined all the parameters in the HHL algorithm. The complete source code is as follows.

import std;
import qft;


/* 
State Preparation for |b>. 
*/
unit state_preparation(qbit q) {
    X(q);
}


/* 
Hamiltonian Simulation.
Here A = (1/2) I - (1/6) X, 
Thus e^{iAt} = e^{i (1/2) I t} e^{i (-1/6) X t}.
Here t = 1.5 * pi.
For more complicated case, Trotter decomposion is a straightforward way.
*/
unit hamiltonian_simulation(qbit q) {
    double alpha = 1.0/2;
    double beta = -1.0/6;
    double t = 1.5 * pi;
    GPhase(alpha * t);
    H(q);
    Rz(-2 * beta * t, q);
    H(q);
} deriving gate


/* 
Phase Estimation.
Here we set 2 qubits for clock register, where q_clock[0] is LSB.
*/
unit phase_estimation(qbit q_input, qbit q_clock[2]) {
    H(q_clock);
    ctrl hamiltonian_simulation(q_clock[0], q_input);
    ctrl hamiltonian_simulation(q_clock[1], q_input);
    ctrl hamiltonian_simulation(q_clock[1], q_input);
    qft_inv(q_clock);
}


/* 
Controlled Rotation.
Here theta_tau = 2 * arcsin (1/tau) for tau = 1,2,3, respectively.
The gate sequence is sum_{tau=1,2,3} |tau><tau| otimes RY(theta_tau).
*/
unit controlled_rotation(qbit q_clock[], qbit q_ancilla) {
    double theta_1 = pi;
    double theta_2 = pi / 3;
    double theta_3 = 0.679673818908;
    nctrl ctrl Ry(theta_1, q_clock[1], q_clock[0], q_ancilla);
    ctrl nctrl Ry(theta_2, q_clock[1], q_clock[0], q_ancilla);
    ctrl ctrl Ry(theta_3, q_clock[1], q_clock[0], q_ancilla);
}


/* 
Inversed Phase Estimation.
*/
unit inv_phase_estimation(qbit q_input, qbit q_clock[2]) {
    qft(q_clock);
    ctrl inv hamiltonian_simulation(q_clock[0], q_input);
    ctrl inv hamiltonian_simulation(q_clock[1], q_input);
    ctrl inv hamiltonian_simulation(q_clock[1], q_input);
    H(q_clock);
}


/* 
The HHL process.
*/
int hhl_process() {
    while (true) {
        qbit q_input;
        qbit q_clock[2];
        qbit q_ancilla;
        state_preparation(q_input);
        phase_estimation(q_input, q_clock);
        controlled_rotation(q_clock, q_ancilla);
        inv_phase_estimation(q_input, q_clock);
        // q_input will be ~ |x> if the measurement outcome of q_ancilla is |1> 
        if (M(q_ancilla)) {
            return M(q_input);
        }
    }
}


/* 
Repeat the HHL process for statistic outcome.
*/
unit statistic(int count) {
    int outcome_array[] = [0, 0];
    for i in 0 : count {
        int outcome = hhl_process();
        outcome_array[outcome] += 1;
    }
    for outcome in outcome_array {
        print outcome;
    }
}


unit main(){
    int count = 100;
    statistic(count);
}

Since the solution is $|x\rangle = (1,3)^T$ without normalization, the printed statistics should be approximately $1:9$ for 0s and 1s.

Reference

1. A. W. Harrow, A. Hassidim, and S. Lloyd, "Quantum Algorithm for Linear Systems of Equations," Phys. Rev. Lett. 103, 150502, 2009.
2. M. A. Nielsen, and I. L. Chuang. "Quantum Computation and Quantum Information: 10th Anniversary Edition." Cambridge University Press, 2010.