✨ Implement HHL algorithm (#582)

Quite some time ago, the HHL implementation based on `qiskit.algorithms`
has been removed. Now, the time has come to re-add this benchmark.

In this current implementation, a simplified version is implemented that
solves a specific linear equation system and the `num_qubits` scale the
number of QPE qubits. This gives the advantages of getting a circuit
that actually has the specified number of qubits in comparison to
before.

However, while I think this implementation could and should be improved
in the future, this might be a fair starting point.

This PR resolves #216.

---------

Signed-off-by: Nils Quetschlich <[email protected]>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

Author: nquetschlich

src/mqt/bench/benchmark_generation.py

@@ -46,6 +46,7 @@ def get_supported_benchmarks() -> list[str]:
         "ghz",
         "graphstate",
         "grover",
+        "hhl",
         "qaoa",
         "qft",
         "qftentangled",

src/mqt/bench/benchmarks/hhl.py

@@ -0,0 +1,78 @@
+# Copyright (c) 2023 - 2025 Chair for Design Automation, TUM
+# Copyright (c) 2025 Munich Quantum Software Company GmbH
+# All rights reserved.
+#
+# SPDX-License-Identifier: MIT
+#
+# Licensed under the MIT License
+
+"""HHL Benchmark Circuit Generation."""
+
+from __future__ import annotations
+
+import numpy as np
+from qiskit import ClassicalRegister, QuantumCircuit, QuantumRegister
+from qiskit.circuit.library import QFTGate
+
+
+def create_circuit(num_qubits: int) -> QuantumCircuit:
+    """HHL algorithm for a fixed 2x2 Hermitian matrix A using scalable QPE precision.
+
+    This implementation simulates a more accurate model of the HHL algorithm
+    for the matrix A = [[1, 1], [1, 3]] with known eigenvalues.
+
+    Args:
+        num_qubits: Number of qubits in the phase estimation register (precision control)
+
+    Returns:
+        QuantumCircuit: Qiskit circuit implementing HHL
+    """
+    assert num_qubits >= 3, "Number of qubits must be at least 3 for HHL."
+    num_qpe_qubits = num_qubits - 2
+    qr_sys = QuantumRegister(1, name="sys")  # System qubit (|b⟩)
+    qr_eig = QuantumRegister(num_qpe_qubits, name="phase")  # Eigenvalue estimation
+    qr_anc = QuantumRegister(1, name="ancilla")  # Ancilla for rotation
+    cr = ClassicalRegister(1, name="c")  # Classical for system measurement
+    qc = QuantumCircuit(qr_sys, qr_eig, qr_anc, cr)
+
+    # Step 1: Prepare |b⟩ = |1⟩
+    qc.x(qr_sys[0])
+
+    # Step 2: Apply Hadamards to phase register (QPE start)
+    qc.h(qr_eig)
+
+    # Step 3: Controlled-unitary approximation simulating controlled-e^{iAt}
+    # A = [[1, 1], [1, 3]], eigenvalues: lam1 ≈ 0.382, lam2 ≈ 3.618
+    # We simulate one eigenvalue's evolution as approximation
+    t = 2 * np.pi
+    lam = 3.618  # simulate dominant eigenvalue for better realism
+
+    for j in range(num_qpe_qubits):
+        angle = t * lam / (2 ** (j + 1))
+        qc.cp(angle, qr_eig[j], qr_sys[0])
+
+    # Step 4: Apply inverse QFT
+    qc.append(QFTGate(num_qpe_qubits).inverse(), qr_eig)
+
+    # Step 5: Controlled Ry rotations on ancilla (based on actual eigenvalue)
+    for j in range(num_qpe_qubits):
+        # Approximate eigenvalue encoded in basis state |j⟩
+        # Use inverse λ (scaled appropriately)
+        estimated_lambda = lam / (2 ** (num_qpe_qubits - j))
+        if estimated_lambda > 0:
+            inv_lambda = 1.0 / estimated_lambda
+            inv_lambda = np.clip(inv_lambda, -1, 1)  # valid for arcsin
+            theta = 2 * np.arcsin(inv_lambda)
+            qc.cry(theta, qr_eig[j], qr_anc[0])
+
+    # Step 6: QPE uncomputation (apply QFT + reverse controlled unitary)
+    qc.append(QFTGate(num_qpe_qubits), qr_eig)
+    for j in reversed(range(num_qpe_qubits)):
+        angle = -t * lam / (2 ** (j + 1))
+        qc.cp(angle, qr_eig[j], qr_sys[0])
+
+    # Step 7: Final Hadamards and measurement
+    qc.h(qr_eig)
+    qc.measure(qr_sys[0], cr[0])
+    qc.name = "hhl"
+    return qc

tests/test_bench.py

@@ -52,6 +52,7 @@
     ghz,
     graphstate,
     grover,
+    hhl,
     qaoa,
     qft,
     qftentangled,
@@ -97,6 +98,7 @@ def output_path() -> str:
         (dj, 3),
         (graphstate, 3),
         (grover, 3),
+        (hhl, 3),
         (qaoa, 3),
         (qft, 3),
         (qftentangled, 3),
@@ -186,6 +188,7 @@ def test_dj_constant_oracle() -> None:
         # Algorithm-level tests
         ("dj", BenchmarkLevel.ALG, 3, None, None),
         ("wstate", BenchmarkLevel.ALG, 3, None, None),
+        ("hhl", BenchmarkLevel.ALG, 3, None, None),
         ("shor", BenchmarkLevel.ALG, 18, None, None),
         ("grover", BenchmarkLevel.ALG, 3, None, None),
         ("qwalk", BenchmarkLevel.ALG, 3, None, None),