Solving linear system problems with HHL#
The Harrow-Hassidim-Lloyd (HHL) quantum algorithm offers an exponential speed-up over classical methods for solving linear system problems \(Ax=b\) for certain sparse matrices \(A\).
The implementation features hybrid quantum-classical workflows and is compiled using Catalyst, a quantum just-in-time (QJIT) compiler framework.
The goal of this demo is to showcase how Qrisp and Catalyst complement each other for implemententing advanced quantum algorithms and compling them for practically relevant problem sizes through the Jasp compilation pipeline.
HHL algorithm in theory#
Given an \(N\)-by-\(N\) Hermitian matrix \(A\) and an \(N\)-dimensional vector \(b\), the Quantum Linear Systems Problem (QSLP) consists of preparing a quantum state \(\ket{x}\) with amplitudes proportional to the solution \(x\) of the linear system of equations \(Ax=b\). Thereby, it can exhibit an exponential speedup over classical methods for certain sparse matrices \(A\). The HHL quantum algorithm and, more generally, quantum linear systems algorithms, hold significant promise for accelerating computations in fields that rely heavily on solving linear systems of equations, such as solving differential equations, or accelerating machine learning.
In its eigenbasis, the matrix \(A\) can be written as
where \(\ket{u_i}\) is an eigenvector of \(A\) corresponding to the eigenvalue \(\lambda_i\).
We define the quantum states \(\ket{b}\) and \(\ket{x}\) as
where \(\ket{b}\) and \(\ket{x}\) are expressed in the eigenbasis of \(A\).
Solving the linerar system amounts to
You might wonder why we can’t just apply \(A^{-1}\) directly to \(\ket{b}\)? This is because, in general, the matix \(A\) is not unitary. However, we will circumnavigate this by exploiting that the Hamiltonian evolution \(U=e^{itA}\) is unitary for a Hemitian matrix \(A\). And this brings us to the HHL algorithm.
In theory, the HHL algorithm can be described as follows:
Step 1: We start by preparing the state
\[\ket{\Psi_1} = \ket{b} = \sum_i \beta_i\ket{u_i}\]Step 2: Applying quantum phase estimation with respect to the Hamiltonian evolution \(U=e^{itA}\) yields the state
\[\ket{\Psi_2} = \sum_i \beta_i\ket{u_i}\ket{\lambda_it/2\pi} = \sum_i \beta_i\ket{u_i}\ket{\widetilde{\lambda}_i}\]To simplify notation, we write \(\widetilde{\lambda}_i=\lambda_it/2\pi\).
Step 3: Performing the inversion of the eigenvalues \(\widetilde{\lambda}_i\rightarrow\widetilde{\lambda}_i^{-1}\) yields the state
\[\ket{\Psi_3} = \sum_i \beta_i\ket{u_i}\ket{\widetilde{\lambda}_i}\ket{\widetilde{\lambda}_i^{-1}}\]Step 4: The amplitudes are multiplied by the inverse eigenvalues \(\widetilde{\lambda}_i^{-1}\) to obtain the state
\[\ket{\Psi_4} = \sum_i \lambda_i^{-1}\beta_i\ket{u_i}\ket{\widetilde{\lambda}_i}\ket{\widetilde{\lambda}_i^{-1}}\]This is achieved by means of a repeat-until-success procedure that applies Steps 1-3 as a subroutine. Stay tuned for more details below!
Step 5: As a final step, we uncompute the variables \(\ket{\widetilde{\lambda}^{-1}}\) and \(\ket{\widetilde{\lambda}}\), and obtain the state
\[\ket{\Psi_5} = \sum_i \lambda_i^{-1}\beta_i\ket{u_i} = \ket{x}\]
And that’s the HHL algorithm. The variable initialized in state \(\ket{b}\) is now found in state \(\ket{x}\)!
As shown in the original paper, the runtime of this algorithm is \(\mathcal{O}(\log(N)s^2\kappa^2/\epsilon)\) where \(s\) and \(\kappa\) are the sparsity and condition number of the matrix \(A\), respectively, and \(\epsilon\) is the precison of the solution. The logarithmic dependence on the dimension \(N\) is the source of an exponential advantage over classical methods.
HHL implementation in practice#
Let’s put theory into practice and dive into an implementation of the HHL algorithm.
As a first step, we define a function fake_inversion that performs the inversion \(\lambda\mapsto\lambda^{-1}\). In this example, we restrict ourselves to an implementation that works for values \(\lambda=2^{-k}\) for \(k\in\mathbb{N}\). (A general inversion is available in Qrisp and will soon be updated to be compatible with QJIT compilation!)
[1]:
from qrisp import *
def fake_inversion(qf, res=None):
if res is None:
res = QuantumFloat(qf.size + 1)
for i in jrange(qf.size):
cx(qf[i], res[qf.size - i])
return res
Note: Essentially, the controlled-NOT operations in the loop reverse the positions of the bits in input variable and place them in the result variable in the opposite order. For example, for \(\lambda=2^{-3}\), which is \(0.001\) in binary, the function would produce \(\lambda^{-1}=2^3\), which in binary is 1000.
Let’s see if it works as intended!
[2]:
qf = QuantumFloat(3, -3)
x(qf[2])
dicke_state(qf, 1)
res = fake_inversion(qf)
print(multi_measurement([qf, res]))
{(0.125, 8): 0.3333333333333333, (0.25, 4): 0.3333333333333333, (0.5, 2): 0.3333333333333333}
Next, we define the function HHL_encoding that performs Steps 1-4 and prepares the state \(\ket{\Psi_4}\). But, how do we get the values \(\widetilde{\lambda}^{-1}_i\) into the amplitudes of the states, i.e. how do we go from \(\ket{\Psi_3}\) to \(\ket{\Psi_4}\)?
Recently, efficient methods for black-box quantum state preparation that avoid arithmetic were proposed, see Bausch, Sanders et al. In this demo, we use a routine proposed in the latter reference which is based on a comparison between integers. This is implemented via the aforementioned comparisons of QuantumFloats.
To simplify the notation, we write \(y^{(i)}=\widetilde{\lambda}^{-1}_i\). Recall that the values \(y^{(i)}\) represent unsigned integers between \(0\) and \(2^n-1\).
Starting from the state
we prepare a uniform superposition of \(2^n\) states in a case_indicator QuantumFloat.
Next we calculate the comparison \(a\geq b\) between the res and the case_indicator into a QuantumBool qbl.
Finally, the case_indicator is unprepared with \(n\) Hadamards and we obtain the state
where \(\ket{\Phi}\) is an orthogonal state with the last variables not in \(\ket{0}_{\text{case}}\ket{0}_{\text{qbl}}\).
Hence, upon measuring the case_indicator in state \(\ket{0}\) and the target qbl in state \(\ket{0}\), the desired state is prepared.
Steps 1-4 are preformed as a repeat-until-success (RUS) routine. This decorator converts the function to be executed within a repeat-until-success (RUS) procedure. The function must return a boolean value as first return value and is repeatedly executed until the first return value is True.
[ ]:
@RUS(static_argnums=[0, 1])
def HHL_encoding(b, hamiltonian_evolution, n, precision):
# Prepare the state |b>. Step 1
qf = QuantumFloat(n)
# Reverse the endianness for compatibility with Hamiltonian simulation.
prepare(qf, b, reversed=True)
qpe_res = QPE(qf, hamiltonian_evolution, precision=precision) # Step 2
inv_res = fake_inversion(qpe_res) # Step 3
case_indicator = QuantumFloat(inv_res.size)
with conjugate(h)(case_indicator):
qbl = case_indicator >= inv_res
cancellation_bool = (measure(case_indicator) == 0) & (measure(qbl) == 0)
# The first return value is a boolean.
# Additional return values are QuantumVariables.
return cancellation_bool, qf, qpe_res, inv_res
The probability of success could be further increased by oblivious amplitude amplification in order to obain an optimal asymptotic scaling.
Finally, we put all things together into the HHL function.
This function takes the follwoing arguments:
bThe vector \(b\).hamiltonian_evolutionA function performing hamiltonian_evolution \(e^{itA}\).nThe number of qubits encoding the state \(\ket{b}\) (\(N=2^n\)).precisionThe precison of the quantum phase estimation.
The HHL function uses the previously defined subroutine to prepare the state \(\ket{\Psi_4}\) and subsequently uncomputes the \(\ket{\widetilde{\lambda}}\) and \(\ket{\lambda}\) quantum variables leaving the first variable, that was initialized in state \(\ket{b}\), in the target state \(\ket{x}\).
[4]:
def HHL(b, hamiltonian_evolution, n, precision):
qf, qpe_res, inv_res = HHL_encoding(b, hamiltonian_evolution, n, precision)
with invert():
QPE(qf, hamiltonian_evolution, target=qpe_res)
fake_inversion(qpe_res, res=inv_res)
# Reverse the endianness for compatibility with Hamiltonian simulation.
for i in jrange(qf.size // 2):
swap(qf[i], qf[n - i - 1])
return qf
Applying HHL to solve systems of linear equations#
Let’s try a first simple example. First, the matrix \(A\) is repesented as a Pauli operator \(H\) and the Hamiltonian evolution unitary \(U=e^{itH}\) is obtained by trotterization with 1 step (as the Pauli terms commute in this case). We choose \(t=\pi\) to ensure that \(\widetilde{\lambda}_i=\lambda_i t/2\pi\) are of the form \(2^{-k}\) for a positive integer \(k\).
This is enabled by the Qrisp’s QubitOperator class providing the tools to describe, optimize and efficiently simulate quantum Hamiltonians.
[ ]:
from qrisp.operators import QubitOperator
import numpy as np
A = np.array([[3 / 8, 1 / 8], [1 / 8, 3 / 8]])
b = np.array([1, 1])
H = QubitOperator.from_matrix(A).to_pauli()
def U(qf):
# By default e^{-itH} is performed. Therefore, we set t=-pi.
H.trotterization()(qf, t=-np.pi, steps=1)
The terminal_sampling decorator performs a hybrid simulation and afterwards samples from the resulting quantum state. We convert the resulting measurement probabilities to amplitudes by appling the square root. Note that, minus signs of amplitudes cannot be recovered from measurement probabilities.
[6]:
@terminal_sampling
def main():
x = HHL(tuple(b), U, 1, 3)
return x
res_dict = main()
for k, v in res_dict.items():
res_dict[k] = v**0.5
print(res_dict)
{0.0: 0.7071067811865476, 1.0: 0.7071067811865476}
Finally, let’s compare to the classical result.
[7]:
x = (np.linalg.inv(A) @ b) / np.linalg.norm(np.linalg.inv(A) @ b)
print(x)
[0.70710678 0.70710678]
And voila! Now, let’s tackle some more complicated examples! Next, we try some randomly generated matrices whose eigenvalues are inverse powers of 2, i.e. of the form \(2^{-k}\) for \(k<K\).
To facilitate fast simulations, we restrict ourselves to \(K=4\) (required precision of QPE) as the runtime of the HHL algorithm scales linearly in the inverse precision \(\epsilon=2^{-K}\) (and therefore exponentially in \(K\)).
[8]:
def hermitian_matrix_with_power_of_2_eigenvalues(n):
# Generate eigenvalues as inverse powers of 2.
eigenvalues = 1 / np.exp2(np.random.randint(1, 4, size=n))
# Generate a random unitary matrix.
Q, _ = np.linalg.qr(np.random.randn(n, n))
# Construct the Hermitian matrix.
A = Q @ np.diag(eigenvalues) @ Q.conj().T
return A
# Example
n = 3
A = hermitian_matrix_with_power_of_2_eigenvalues(2**n)
H = QubitOperator.from_matrix(A).to_pauli()
def U(qf):
H.trotterization()(qf, t=-np.pi, steps=5)
b = np.random.randint(0, 2, size=2**n)
print("Hermitian matrix A:")
print(A)
print("Eigenvalues:")
print(np.linalg.eigvals(A))
print("b:")
print(b)
Hermitian matrix A:
[[ 0.27161849 0.04084016 -0.02994612 -0.01343795 0.00456208 -0.02389833
-0.03865614 -0.01028375]
[ 0.04084016 0.35006745 -0.06348883 0.04375543 0.00892184 -0.05256036
-0.05969854 -0.03136835]
[-0.02994612 -0.06348883 0.39879834 0.02868921 -0.05439778 0.0555282
-0.05704449 0.00337529]
[-0.01343795 0.04375543 0.02868921 0.47607339 -0.01603151 -0.00157681
0.03290459 -0.0338936 ]
[ 0.00456208 0.00892184 -0.05439778 -0.01603151 0.27284976 -0.01434681
0.04061657 0.00427224]
[-0.02389833 -0.05256036 0.0555282 -0.00157681 -0.01434681 0.2826893
0.01797763 0.01245469]
[-0.03865614 -0.05969854 -0.05704449 0.03290459 0.04061657 0.01797763
0.43479804 0.02610161]
[-0.01028375 -0.03136835 0.00337529 -0.0338936 0.00427224 0.01245469
0.02610161 0.26310524]]
Eigenvalues:
[0.25 0.25 0.5 0.5 0.5 0.25 0.25 0.25]
b:
[0 0 0 1 0 0 1 1]
[9]:
@terminal_sampling
def main():
x = HHL(tuple(b), U, n, 4)
return x
res_dict = main()
for k, v in res_dict.items():
res_dict[k] = v**0.5
np.array([res_dict[key] for key in sorted(res_dict)])
[9]:
array([0.09903807, 0.08063741, 0.03992062, 0.43921894, 0.04622294,
0.04614081, 0.40873337, 0.78603666])
Let’s compare to the classical solution:
[10]:
x = (np.linalg.inv(A) @ b) / np.linalg.norm(np.linalg.inv(A) @ b)
print(x)
[ 0.09933223 0.0753401 0.03977885 0.43778338 -0.04595318 -0.04595032
0.40797335 0.7877531 ]
Yup, close enough… That’s all folks!
Step-by-step recap#
Let’s rewind for a second, take a deep breath, and go through the steps and concepts you learned so far.
Equipped with a theoretical introduction to HHL and outlining the steps required to perform this algorithm, you got to see how to first encode the first 4 steps and making use of the repeat until success feature of Jasp. Then, putting everything together, we combined the previously defined building blocks (read: Python functions) - the HHL_encoding and QPE - into a simple function. With a brief feature apperance of Hamiltonian simulation you then successfully managed to solve two systems of linear equations.
In conclusion, let’s take a moment to appreciate one last time how elegantly we can call the HHL algorithm:
x = HHL(b, hamiltonian_evolution, n, precision)
As qrispy as always!
QJIT compilation with Catalyst#
The solution that we saw so far only ran within the Python-based Qrisp internal simulator, which doesn’t have to care about silly things such as coherence time. Unfortunately, this is not (yet) the case for real quantum hardware so it is of paramount importance that the classical real-time computations (such as the while loop within the repeat-until-success steps) are executed as fast as possible. For this reason the QIR specification was created.
QIR essentially embeds quantum aspects into LLVM, which is the foundation of a lot of modern compiler infrastructure. This implies QIR based software stacks are able to integrate a large part of established classical software and also express real-time control structures. Read more about QIR at the QIR Alliance website and the Advancing hybrid quantum–classical computation with real-time execution paper.
Jasp has been built with a direct Catalyst integration implying Jasp programs can be converted to QIR via the Catalyst pipeline. This conversion is straightforward: You simply capture the Qrisp computation using the make_jaspr function and then call to_qir.
[11]:
def main():
x = HHL(tuple(b), U, n, 4)
# Note that we have to return a classical value
# (in this case the measurement result of the
# quantum variable returned by the HHL algorithm)
# Within the above examples, we used the terminal_sampling
# decorator, which is a convenience feature and allows
# a much faster sampling procedure.
# The terminal_sampling decorator expects a function returning
# quantum variables, while most other evaluation modes require
# classical return values.
return measure(x)
jaspr = make_jaspr(main)()
qir_str = jaspr.to_qir()
# Print only the first few lines - the whole string is very long.
print(qir_str[:2000])
; ModuleID = 'LLVMDialectModule'
source_filename = "LLVMDialectModule"
target datalayout = "e-m:o-i64:64-i128:128-n32:64-S128-Fn32"
target triple = "arm64-apple-darwin24.5.0"
@"{'shots': 0, 'mcmc': False, 'num_burnin': 0, 'kernel_name': None}" = internal constant [66 x i8] c"{'shots': 0, 'mcmc': False, 'num_burnin': 0, 'kernel_name': None}\00"
@LightningSimulator = internal constant [19 x i8] c"LightningSimulator\00"
@"/opt/anaconda3/envs/qrisp/lib/python3.10/site-packages/pennylane_lightning/liblightning_qubit_catalyst.dylib" = internal constant [109 x i8] c"/opt/anaconda3/envs/qrisp/lib/python3.10/site-packages/pennylane_lightning/liblightning_qubit_catalyst.dylib\00"
@__constant_xi64_3 = private constant i64 1
@__constant_xi1 = private constant i1 false
@__constant_xi64_2 = private constant i64 0
@__constant_1024xi64 = private constant [1024 x i64] zeroinitializer
@__constant_xi64_1 = private constant i64 30
@__constant_xi64_0 = private constant i64 4
@__constant_xi64 = private constant i64 3
@__constant_30xi64 = private constant [30 x i64] [i64 30, i64 29, i64 28, i64 27, i64 26, i64 25, i64 24, i64 23, i64 22, i64 21, i64 20, i64 19, i64 18, i64 17, i64 16, i64 15, i64 14, i64 13, i64 12, i64 11, i64 10, i64 9, i64 8, i64 7, i64 6, i64 5, i64 4, i64 3, i64 2, i64 1]
@__constant_xf64_34 = private constant double 0xBF8D02B4A78574FA
@__constant_xf64_33 = private constant double 0xBFCBA5614317CB35
@__constant_xf64_32 = private constant double 0x3F998A8430AEF7BA
@__constant_xf64_31 = private constant double 0x3F9FCFE0804A8850
@__constant_xf64_30 = private constant double 0x3F807D3742E47A2C
@__constant_xf64_29 = private constant double 0xBF6B69C1DC783640
@__constant_xf64_28 = private constant double 0xBF40F04DDF80EBDC
@__constant_xf64_27 = private constant double 0xBF840EC9BEB08EA8
@__constant_xf64_26 = private constant double 0xBF80EBD1A8EC99B0
@__constant_xf64_25 = private constant double 0xBF8296060CFC377A
@__constant_xf64_24 = private constant double 0x3F7234DFA
The Catalyst runtime#
Jasp is also capable of targeting the Catalyst execution runtime (i.e., the Lightning simulator). There are, however, still some simulator features to be implemented on Jasp side, which prevents efficient sampling. We restrict the demonstration to the smaller examples from above to limit the overall execution time required. (Warning: The execution may take more than 15 minutes.)
[12]:
A = np.array([[3 / 8, 1 / 8], [1 / 8, 3 / 8]])
b = np.array([1, 1])
H = QubitOperator.from_matrix(A).to_pauli()
# By default e^{-itH} is performed. Therefore, we set t=-pi.
def U(qf):
H.trotterization()(qf, t=-np.pi, steps=1)
@qjit
def main():
x = HHL(tuple(b), U, 1, 3)
return measure(x)
samples = []
for i in range(5):
samples.append(float(main()))
print(samples)
[0.0, 0.0, 1.0, 0.0, 0.0]
Scrolling back to the terminal_sampling cell, we see that the expected distribution is 50/50 between one and zero, which roughly agrees to the result of the previous cell and concludes this demo.
Conclusion#
In this demo, we have shown how to implement the HHL algorithm featuring classical real-time computations in the high-level language Qrisp. This algorithm is important in a variety of use cases such as solving differential equations, accelerating machine learning, and more generally, any task that involves solving linear systems of equations. Along the way, we have dipped into advanced concepts such as Linear Combination of Unitaries and Hamiltonian simulation. Moreover, we have demonstrated how Qrisp and Catalyst complement each other for translating a high-level implementation into low-level QIR.