# This original code is sourced from:
# https://github.com/SRI-International/QC-App-Oriented-Benchmarks/blob/master/hhl/qiskit/hhl_benchmark.py
# https://github.com/SRI-International/QC-App-Oriented-Benchmarks/blob/master/hhl/qiskit/sparse_Ham_sim.py
#
# Original repository link:
# https://github.com/SRI-International/QC-App-Oriented-Benchmarks?tab=Apache-2.0-1-ov-file
#
#
# This program is adapted for use as a benchmark in controlled software testing experiments.
# Modifications made to the original code include:
# - Stucture reorganization:
# - Modify the code in the `QuantumCircuit` class
# - Add type hints for better code readability
# - Remove measurements from subcircuits, as they can be added in the main circuit if needed
# - Dependency decoupling:
# - Copy `_alpha2theta` along with its dependent functions within this module
# - Replace the Hamiltonian simulation code with the `SparseHamiltonian` class from the `hamiltonian` module
#
#
# Copyright [yyyy] [name of copyright owner]
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""
HHL Benchmark Program - Qiskit
"""
import numpy as np
pi = np.pi
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.circuit.library import U1Gate
from qolumbina.programs.hamiltonian import SparseHamiltonian
from qolumbina.programs.pauli_rotations import UniformlyControlledRotations
from sympy.combinatorics.graycode import GrayCode
# saved circuits for display
QC_ = None
U_ = None
UI_ = None
QFT_ = None
QFTI_ = None
HP_ = None
INVROT_ = None
# ---------- benchmark registration ----------
from ..benchmark_registry import register_benchmark
from pathlib import Path
@register_benchmark(
Path(__file__).stem,
family=Path(__file__).resolve().parent.name,
description="Harrow-Hassidim-Lloyd (HHL) algorithm for solving linear systems of equations",
class_name="HarrowHassidimLloyd",
source={
"repo": "https://github.com/SRI-International/QC-App-Oriented-Benchmarks?tab=Apache-2.0-1-ov-file",
"file": "hhl/qiskit/hhl_benchmark.py",
"sdk": "Qiskit",
"available_doc": True
},
testability_refactoring=[
"Dependency decoupling",
"Structure reorganization",
"Input validation"
]
)
def create_harrow_hassidim_lloyd(
A: np.ndarray,
b: int,
num_clock_qubits: int,
C: float = 1/4,
if_barrier: bool = False
):
return HarrowHassidimLloyd(
A=A,
b=b,
num_clock_qubits=num_clock_qubits,
C=C,
if_barrier=if_barrier
)
############### Circuit Definitions
''' replaced with code below ...
def qft_dagger(qc, clock, n):
qc.h(clock[1]);
for j in reversed(range(n)):
for k in reversed(range(j+1,n)):
qc.cu1(-np.pi/float(2**(k-j)), clock[k], clock[j]);
qc.h(clock[0]);
def qft(qc, clock, n):
qc.h(clock[0]);
for j in reversed(range(n)):
for k in reversed(range(j+1,n)):
qc.cu1(np.pi/float(2**(k-j)), clock[k], clock[j]);
qc.h(clock[1]);
'''
'''
DEVNOTE: the QFT and IQFT are defined here as they are in the QFT benchmark - almost;
Here, the sign of the angles is reversed and the QFT is actually used as the inverse QFT.
This is an inconsistency that needs to be resolved later.
The QPE part of the algorithm should be using the inverse QFT, but the qubit order is not correct.
The QFT as defined in the QFT benchmark operates on qubits in the opposite order from the HHL pattern.
'''
[docs]
class HarrowHassidimLloyd(QuantumCircuit):
r""" Generate top-level circuit for HHL algorithm :math:`A\vec{x} = \vec{b}`
"""
def __init__(
self,
A: np.ndarray,
b: int,
num_clock_qubits: int,
C: float = 1/4,
if_barrier: bool = False,
name: str | None = None
):
r"""
Args:
A: The Hermitian matrix $A$ in the linear equation :math:`A\vec{x} = \vec{b}`.
b: The integer to initialize the b-register to the basis state :math:`\ket{b}_n`.
num_clock_qubits: The number of clock qubits :math:`n_t` used in the QPE.
C: The normalization constant used in the controlled rotation step. Default is ``1/4``.
if_barrier: Whether to add barriers in the circuit.
name: Optional name for the circuit.
Raises:
ValueError: If ``b`` is not between 0 and the dimension of ``A`` - 1.
ValueError: If ``C`` is not between 0 and 1.
ValueError: If ``A`` is not a valid input matrix
(e.g., not Hermitian, not sparse, or singular).
"""
self._A = A
self._b = b
self._num_clock_qubits = num_clock_qubits
self._C = C
self._if_barrier = if_barrier
# b is an integer representing the initial basis state |b>
if b < 0 or b >= A.shape[0]:
raise ValueError(f"b must be between 0 and {A.shape[0]-1}, but got b={b}")
# C must be positive and less than 1
if C <= 0 or C > 1:
raise ValueError("C must be between (0, 1]")
super().__init__(2*int(np.log2(A.shape[0]))+num_clock_qubits+1, name=name or "HHL")
self._build()
def _validate_matrix_A(self, A: np.ndarray) -> None:
""" Validate the input matrix A for HHL algorithm """
# Hermitian check
if not np.allclose(A, A.conj().T):
raise ValueError("Matrix A must be Hermitian (A != A.conj().T)")
# Sparse check (here assume "s-sparse": max nonzeros per row/col <= s)
s = np.max(np.sum(A != 0, axis=0)) # Non-zero element numbers in each column
if s >= A.shape[0]: # For n x n matrix, s should be < n
raise ValueError(f"Matrix A must be sparse. Max nonzeros per row/col = {s}")
# Non-singular / minimal eigenvalue
eigvals = np.linalg.eigvalsh(A) # Hermitian matrix eigenvalues
if np.any(np.isclose(eigvals, 0)):
raise ValueError("Matrix A is singular (has zero eigenvalues)")
[docs]
def dot_product(self, str1: str, str2: str) -> int:
r"""
Compute the dot product between 2 binary strings
"""
prod = 0
for j in range(len(str1)):
if str1[j] == '1' and str2[j] == '1':
prod = (prod + 1)%2
return prod
[docs]
def conversion_matrix(self, N: int) -> np.ndarray:
r"""
Compute the conversion matrix from alpha to theta
"""
M = np.zeros((N,N))
n = int(np.log2(N))
gc_list = list(GrayCode(n).generate_gray()) # list of gray code strings
for i in range(N):
g_i = gc_list[i]
for j in range(N):
b_j = np.binary_repr(j, width=n)[::-1]
M[i,j] = (-1)**self.dot_product(g_i, b_j)/(2**n)
return M
[docs]
def alpha2theta(self, alpha: list[float]) -> np.ndarray:
r"""
Args:
alpha : list of angles that get applied controlled on :math:`0,...,2^n-1`
theta : list of angles occuring in circuit construction
"""
N = len(alpha)
M = self.conversion_matrix(N)
theta = M @ np.array(alpha)
return theta
[docs]
def initialize_state(self, qc: QuantumCircuit, qreg: QuantumRegister, b: int) -> QuantumCircuit:
r"""
Args:
b: initial basis state :math:`\ket{b}_n` to prepare on qreg
"""
n = qreg.size
b_bin = np.binary_repr(b, width=n)
for q in range(n):
if b_bin[n-1-q] == '1':
qc.x(qreg[q])
return qc
[docs]
def IQFT(self, qreg: QuantumRegister) -> None:
r"""
inverse QFT
Args:
qreg : QuantumRegister belonging to qc
Note:
This method does not include SWAP at end of the circuit
"""
n = int(qreg.size)
for i in reversed(range(n)):
for j in range(i+1,n):
phase = -pi/2**(j-i)
self.cp(phase, qreg[i], qreg[j])
[docs]
def QFT(self, qreg: QuantumRegister) -> None:
r"""
QFT
Args:
qreg : QuantumRegister belonging to qc
Note:
This method does not include SWAP at end of circuit
"""
n = int(qreg.size)
for i in range(n):
self.h(qreg[i])
for j in reversed(range(i+1,n)):
phase = pi/2**(j-i)
self.cp(phase, qreg[i], qreg[j])
[docs]
def inv_qft_gate(self, input_size: int, method: int = 1) -> QuantumCircuit:
#global QFT_
qr = QuantumRegister(input_size)
#qc = QuantumCircuit(qr, name="qft")
qc = QuantumCircuit(qr, name="IQFT")
if method == 1:
# Generate multiple groups of diminishing angle CRZs and H gate
for i_qubit in range(0, input_size):
# start laying out gates from highest order qubit (the hidx)
hidx = input_size - i_qubit - 1
# if not the highest order qubit, add multiple controlled RZs of decreasing angle
if hidx < input_size - 1:
num_crzs = i_qubit
for j in range(0, num_crzs):
divisor = 2 ** (num_crzs - j)
#qc.crz( math.pi / divisor , qr[hidx], qr[input_size - j - 1])
##qc.crz( -np.pi / divisor , qr[hidx], qr[input_size - j - 1])
qc.cp(-np.pi / divisor, qr[hidx], qr[input_size - j - 1])
# followed by an H gate (applied to all qubits)
qc.h(qr[hidx])
elif method == 2:
# apply IQFT to register
for i in range(input_size)[::-1]:
for j in range(i+1,input_size):
phase = -np.pi/2**(j-i)
qc.cp(phase, qr[i], qr[j])
qc.h(qr[i])
if self._if_barrier:
qc.barrier()
return qc
[docs]
def qft_gate(self, input_size: int, method: int = 1):
#global QFTI_
qr = QuantumRegister(input_size)
#qc = QuantumCircuit(qr, name="inv_qft")
qc = QuantumCircuit(qr, name="QFT")
if method == 1:
# Generate multiple groups of diminishing angle CRZs and H gate
for i_qubit in reversed(range(0, input_size)):
# start laying out gates from highest order qubit (the hidx)
hidx = input_size - i_qubit - 1
# precede with an H gate (applied to all qubits)
qc.h(qr[hidx])
# if not the highest order qubit, add multiple controlled RZs of decreasing angle
if hidx < input_size - 1:
num_crzs = i_qubit
for j in reversed(range(0, num_crzs)):
divisor = 2 ** (num_crzs - j)
#qc.crz( -math.pi / divisor , qr[hidx], qr[input_size - j - 1])
##qc.crz( np.pi / divisor , qr[hidx], qr[input_size - j - 1])
qc.cp( np.pi / divisor , qr[hidx], qr[input_size - j - 1])
elif method == 2:
# apply QFT to register
for i in range(input_size):
qc.h(qr[i])
for j in range(i+1, input_size):
phase = np.pi/2**(j-i)
qc.cp(phase, qr[i], qr[j])
if self._if_barrier:
qc.barrier()
return qc
[docs]
def ctrl_u(self, exponent: int) -> tuple:
r"""
Construct the U gates for A
Args:
exponent: exponent of U
"""
qc = QuantumCircuit(1, name=f"U^{exponent}")
for i in range(exponent):
#qc.u(np.pi/2, -np.pi/2, np.pi/2, 3*np.pi/4, target);
#qc.cu(np.pi/2, -np.pi/2, np.pi/2, 3*np.pi/4, control, target);
qc.u(np.pi/2, -np.pi/2, np.pi/2, 0)
cu_gate = qc.to_gate().control(1)
return cu_gate, qc
[docs]
def ctrl_ui(self, exponent: int) -> tuple:
"""
Construct the U^-1 gates for A
Args:
exponent: exponent of U
"""
qc = QuantumCircuit(1, name=f"U^-{exponent}")
for i in range(exponent):
#qc.u(np.pi/2, -np.pi/2, np.pi/2, 3*np.pi/4, target);
#qc.cu(np.pi/2, -np.pi/2, np.pi/2, 3*np.pi/4, control, target);
qc.u(np.pi/2, np.pi/2, -np.pi/2, 0)
cu_gate = qc.to_gate().control(1)
return cu_gate, qc
[docs]
def qpe(
self,
clock: QuantumRegister,
target: QuantumRegister,
extra_qubits: QuantumRegister,
ancilla: QuantumRegister,
A: np.ndarray,
method: int = 1
):
if self._if_barrier:
self.barrier()
r''' original code from Hector_Wong
# e^{i*A*t}
#qc.cu(np.pi/2, -np.pi/2, np.pi/2, 3*np.pi/4, clock[0], target, label='U');
# The CU gate is equivalent to a CU1 on the control bit followed by a CU3
qc.u1(3*np.pi/4, clock[0]);
qc.cu3(np.pi/2, -np.pi/2, np.pi/2, clock[0], target);
# e^{i*A*t*2}
#qc.cu(np.pi, np.pi, 0, 0, clock[1], target, label='U2');
qc.cu3(np.pi, np.pi, 0, clock[1], target);
qc.barrier();
'''
# apply series of controlled U operations to the state |1>
# does nothing to state |0>
# DEVNOTE: have not found a way to create a controlled operation that contains a U gate
# with the global phase; instead do it piecemeal for now
if method == 1:
repeat = 1
#for j in reversed(range(len(clock))):
for j in (range(len(clock))):
# create U with exponent of 1, but in a loop repeating N times
for k in range(repeat):
# this global phase is applied to clock qubit
# qc.u1() is replaced by the following code
qc_u1 = U1Gate(3*np.pi/4)
self.append(qc_u1, [clock[j]])
# apply the rest of U controlled by clock qubit
#cp, _ = ctrl_u(repeat)
cp, _ = self.ctrl_u(1)
self.append(cp, [clock[j], target])
repeat *= 2
if self._if_barrier:
self.barrier()
#Define global U operator as the phase operator (for printing later)
_, U_ = self.ctrl_u(1)
if method == 2:
for j in range(len(clock)):
control = clock[j]
phase = -(2*np.pi)*2**j
# con_H_sim = shs.control_Ham_sim(A, phase)
con_H_sim = SparseHamiltonian(A, phase, controlled_ver=True)
qubits = [control] + [q for q in target] + [q for q in extra_qubits] + [ancilla[0]]
self.append(con_H_sim, qubits)
# Perform an inverse QFT on the register holding the eigenvalues
self.append(self.inv_qft_gate(len(clock), method), clock)
[docs]
def inv_qpe(
self,
clock: QuantumRegister,
target: QuantumRegister,
extra_qubits: QuantumRegister,
ancilla: QuantumRegister,
A: np.ndarray,
method: int = 1
):
"""
Append a series of Inverse Quantum Phase Estimation gates to the circuit
Args:
clock : clock register qubits
target : target register qubits
extra_qubits : extra qubits for Hamiltonian simulation
ancilla : ancilla qubit
A : sparse Hermitian matrix
method : method for QFT implementation
"""
# Perform a QFT on the register holding the eigenvalues
self.append(self.qft_gate(len(clock), method), clock)
if self._if_barrier:
self.barrier()
if method == 1:
''' original code from Hector_Wong
# e^{i*A*t*2}
#qc.cu(np.pi, np.pi, 0, 0, clock[1], target, label='U2');
qc.cu3(np.pi, np.pi, 0, clock[1], target);
# e^{i*A*t}
#qc.cu(np.pi/2, np.pi/2, -np.pi/2, -3*np.pi/4, clock[0], target, label='U');
# The CU gate is equivalent to a CU1 on the control bit followed by a CU3
qc.u1(-3*np.pi/4, clock[0]);
qc.cu3(np.pi/2, np.pi/2, -np.pi/2, clock[0], target);
qc.barrier()
'''
# apply inverse series of controlled U operations to the state |1>
# does nothing to state |0>
# DEVNOTE: have not found a way to create a controlled operation that contains a U gate
# with the global phase; instead do it piecemeal for now
repeat = 2 ** (len(clock) - 1)
for j in reversed(range(len(clock))):
#for j in (range(len(clock))):
# create U with exponent of 1, but in a loop repeating N times
for k in range(repeat):
# this global phase is applied to clock qubit
qc_u1 = U1Gate(-3*np.pi/4)
self.append(qc_u1, [clock[j]])
# apply the rest of U controlled by clock qubit
#cp, _ = ctrl_u(repeat)
cp, _ = self.ctrl_ui(1)
self.append(cp, [clock[j], target])
repeat = int(repeat / 2)
self.barrier()
#Define global U operator as the phase operator (for printing later)
_, UI_ = self.ctrl_ui(1)
if method == 2:
for j in reversed(range(len(clock))):
control = clock[j]
phase = (2*np.pi)*2**j
# con_H_sim = shs.control_Ham_sim(A, phase)
con_H_sim = SparseHamiltonian(A, phase, controlled_ver=True)
qubits = [control] + [q for q in target] + [q for q in extra_qubits] + [ancilla[0]]
self.append(con_H_sim, qubits)
def _build(self):
""" Generate top-level circuit for HHL algo A|x>=|b>
A : sparse Hermitian matrix
b (int): between 0,...,2^n-1. Initial basis state |b>
"""
# save smaller circuit example for display
global QC_, U_, UI_, QFT_, QFTI_, HP_, INVROT_
# read in number of qubits
n_t = self._num_clock_qubits # number of qubits in clock register
b = self._b
C = self._C
A = self._A
N = len(A) # dimension of A
n = int(np.log2(N)) # number of qubits in input register
# num_qubits = 2*n + n_t + 1
''' Define sets of qubits for this algorithm '''
# Create 'input' quantum and classical measurement register
qr = QuantumRegister(n, name='input')
qr_b = QuantumRegister(n, name='in_anc') # ancillas for Hamiltonian simulation (?)
# Create 'clock' quantum register
qr_t = QuantumRegister(n_t, name='clock') # for phase estimation
# Create 'ancilla' quantum and classical measurement register
qr_a = QuantumRegister(1, name='ancilla') # ancilla qubit
# Create the top-level HHL circuit, with all the registers
qc = QuantumCircuit(qr, qr_b, qr_t, qr_a)
# Initialize the |b> state - the 'input'
qc = self.initialize_state(qc, qr, b)
if self._if_barrier:
qc.barrier()
# Hadamard the phase estimation register - the 'clock'
for q in range(n_t):
qc.h(qr_t[q])
if self._if_barrier:
qc.barrier()
''' Perform Quantum Phase Estimation on input (b), clock, and ancilla '''
# perform controlled e^(i*A*t)
for q in range(n_t):
control = qr_t[q]
anc = qr_a[0]
phase = -(2*pi)*2**q
# con_H_sim = shs.control_Ham_sim(n, A, phase)
qc_u = SparseHamiltonian(A, phase, controlled_ver=True)
if phase <= 0:
qc_u.name = "e^{" + str(q) + "iAt}"
else:
qc_u.name = "e^{-" + str(q) + "iAt}"
if U_ == None:
U_ = qc_u
qc.append(qc_u, qr[0:len(qr)] + qr_b[0:len(qr_b)] + [control] + [anc])
if self._if_barrier:
qc.barrier()
''' Perform Inverse Quantum Fourier Transform on clock qubits '''
#qc = IQFT(qc, qr_t)
qc_qfti = self.inv_qft_gate(n_t, method=2)
qc.append(qc_qfti, qr_t)
if QFTI_ == None:
QFTI_ = qc_qfti
if self._if_barrier:
qc.barrier()
''' Perform inverse rotation with ancilla '''
# reset ancilla
qc.reset(qr_a[0])
# compute angles for inversion rotations
alpha = [2*np.arcsin(C)]
for x in range(1, 2**n_t):
x_bin_rev = np.binary_repr(x, width=n_t)[::-1]
lam = int(x_bin_rev,2)/(2**n_t)
if lam < C:
alpha.append(0)
elif lam >= C:
alpha.append(2*np.arcsin(C/lam))
# create uniformly controlled rotation
theta = self.alpha2theta(alpha)
# do inversion step
# from `qc_invrot = ucr.uniformly_controlled_rot(n_t, theta)`
qc_invrot = UniformlyControlledRotations(theta_list=theta.tolist())
qc.append(qc_invrot, qr_t[0:len(qr_t)] + [qr_a[0]])
if INVROT_ == None:
INVROT_ = qc_invrot
# and measure ancilla
# qc.measure(qr_a[0], cr_a[0])
# qc.reset(qr_a[0])
if self._if_barrier:
qc.barrier()
''' Perform Quantum Fourier Transform on clock qubits '''
#qc = QFT(qc, qr_t)
qc_qft = self.qft_gate(n_t, method=2)
qc.append(qc_qft, qr_t)
if QFT_ == None:
QFT_ = qc_qft
if self._if_barrier:
qc.barrier()
''' Perform Inverse Quantum Phase Estimation on input (b), clock, and ancilla '''
# uncompute phase estimation
# perform controlled e^(-i*A*t)
for q in reversed(range(n_t)):
control = qr_t[q]
phase = (2*pi)*2**q
# qc_ui = shs.control_Ham_sim(n, A, phase)
qc_ui = SparseHamiltonian(A, phase, controlled_ver=True)
if phase <= 0:
qc_ui.name = "e^{" + str(q) + "iAt}"
else:
qc_ui.name = "e^{-" + str(q) + "iAt}"
if UI_ == None:
UI_ = qc_ui
qc.append(qc_ui, qr[0:len(qr)] + qr_b[0:len(qr_b)] + [control] + [anc])
if self._if_barrier:
qc.barrier()
# Hadamard (again) the phase estimation register - the 'clock'
for q in range(n_t):
qc.h(qr_t[q])
if self._if_barrier:
qc.barrier()
''' Perform final measurements '''
# measure ancilla and main register
# qc.measure(qr[0:], cr[0:])
if QC_ == None:
QC_ = qc
#print(f"... made circuit = \n{QC_}")
self.compose(qc, inplace=True)