From 37cf5813845f2feb590e4459a369d0798ab851f7 Mon Sep 17 00:00:00 2001 From: Ashwin Kaliyaperumal Date: Wed, 12 Feb 2025 17:01:40 -0800 Subject: [PATCH 1/3] created solovay kitaev implementation created SV algorith implementation with temporary testing scaffolding at the bottom, doesn't work right now -> approx matrix is considerably off from target --- src/solovay_kitaev.py | 158 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 158 insertions(+) create mode 100644 src/solovay_kitaev.py diff --git a/src/solovay_kitaev.py b/src/solovay_kitaev.py new file mode 100644 index 0000000..73cc6f3 --- /dev/null +++ b/src/solovay_kitaev.py @@ -0,0 +1,158 @@ +import numpy as np +from typing import List, Tuple + +class SU2Matrix: + """Class representing a 2x2 Special Unitary matrix.""" + def __init__(self, matrix: np.ndarray): + self.matrix = matrix + + def __mul__(self, other: 'SU2Matrix') -> 'SU2Matrix': + return SU2Matrix(np.dot(self.matrix, other.matrix)) + + def dagger(self) -> 'SU2Matrix': + """Returns the conjugate transpose.""" + return SU2Matrix(self.matrix.conj().T) + + def distance(self, other: 'SU2Matrix') -> float: + """Calculates the operator norm distance between two matrices.""" + diff = self.matrix - other.matrix + return np.linalg.norm(diff) + +def group_commutator(a: SU2Matrix, b: SU2Matrix) -> SU2Matrix: + """Compute the group commutator [a,b] = aba^{-1}b^{-1}.""" + return a * b * a.dagger() * b.dagger() + +def find_basic_approximation(target: SU2Matrix, basic_gates: List[SU2Matrix]) -> SU2Matrix: + """ + Find the best approximation of the target unitary from the basic gate set. + + Args: + target: Target unitary matrix to approximate + basic_gates: List of available basic gates + + Returns: + Best approximating unitary from the basic gate set + """ + best_dist = float('inf') + best_gate = None + + for gate in basic_gates: + dist = target.distance(gate) + if dist < best_dist: + best_dist = dist + best_gate = gate + + return best_gate + +def decompose_group_element( + target: SU2Matrix, + basic_gates: List[SU2Matrix], + depth: int +) -> Tuple[List[SU2Matrix], float]: + """ + Decompose a target unitary into a sequence of basic gates using Solovay-Kitaev. + + Args: + target: Target unitary matrix to decompose + basic_gates: List of available basic gates + depth: Recursion depth for the algorithm + + Returns: + Tuple of (sequence of gates, approximation error) + """ + if depth == 0: + best_approx = find_basic_approximation(target, basic_gates) + return [best_approx], target.distance(best_approx) + + # Recursive approximation + prev_sequence, prev_error = decompose_group_element(target, basic_gates, depth - 1) + + # If previous approximation is good enough, return it + if prev_error < 1e-6: + return prev_sequence, prev_error + + # Compute the error unitary + approx = prev_sequence[0] + for gate in prev_sequence[1:]: + approx = approx * gate + + error = target * approx.dagger() + + # Find Va and Vb such that their group commutator approximates the error + # This is a simplified version - in practice, you'd need a more sophisticated search + best_va = None + best_vb = None + best_error = float('inf') + + for va in basic_gates: + for vb in basic_gates: + comm = group_commutator(va, vb) + curr_error = error.distance(comm) + if curr_error < best_error: + best_error = curr_error + best_va = va + best_vb = vb + + # Construct the final sequence + result_sequence = [] + result_sequence.extend(prev_sequence) + + # Add correction terms + if best_va is not None and best_vb is not None: + result_sequence.extend([best_va, best_vb, best_va.dagger(), best_vb.dagger()]) + + # Calculate final error + final_approx = result_sequence[0] + for gate in result_sequence[1:]: + final_approx = final_approx * gate + final_error = target.distance(final_approx) + + return result_sequence, final_error + +def solovay_kitaev( + target: np.ndarray, + basic_gates: List[np.ndarray], + depth: int = 3 +) -> List[np.ndarray]: + """ + Main function to run the Solovay-Kitaev algorithm. + + Args: + target: Target unitary matrix as numpy array + basic_gates: List of basic gates as numpy arrays + depth: Recursion depth + + Returns: + List of gates that approximate the target unitary + """ + # Convert inputs to SU2Matrix objects + target_su2 = SU2Matrix(target) + basic_gates_su2 = [SU2Matrix(gate) for gate in basic_gates] + + # Run the decomposition + sequence, error = decompose_group_element(target_su2, basic_gates_su2, depth) + + # Convert back to numpy arrays + return [gate.matrix for gate in sequence] + +# Example usage: +if __name__ == "__main__": + # Define some basic gates (Pauli matrices and their combinations) + I = np.array([[1, 0], [0, 1]], dtype=complex) + X = np.array([[0, 1], [1, 0]], dtype=complex) + Y = np.array([[0, -1j], [1j, 0]], dtype=complex) + Z = np.array([[1, 0], [0, -1]], dtype=complex) + H = np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2) + + basic_gates = [I, X, Y, Z, H] + + # Define a target unitary + theta = np.pi / 8 + target = np.array([ + [np.cos(theta), -np.sin(theta)], + [np.sin(theta), np.cos(theta)] + ], dtype=complex) + + # Run the algorithm + sequence = solovay_kitaev(target, basic_gates, depth=3) + print(f"Found sequence of {len(sequence)} gates") \ No newline at end of file From e1d7a606beb5c1ba58cdd0033fa2729d966b2e4b Mon Sep 17 00:00:00 2001 From: Ashwin Kaliyaperumal Date: Wed, 12 Feb 2025 17:07:22 -0800 Subject: [PATCH 2/3] created solovay kitaev implementation created SV algorithm implementation with temporary testing scaffolding at the bottom, doesn't work right now -> approx matrix is considerably off from target --- src/pyqasm/maps/decomposition_rules.py | 133 ++++++++++++++++ src/pyqasm/maps/solovay-kitaev.py | 202 +++++++++++++++++++++++++ 2 files changed, 335 insertions(+) create mode 100644 src/pyqasm/maps/decomposition_rules.py create mode 100644 src/pyqasm/maps/solovay-kitaev.py diff --git a/src/pyqasm/maps/decomposition_rules.py b/src/pyqasm/maps/decomposition_rules.py new file mode 100644 index 0000000..57b2bf6 --- /dev/null +++ b/src/pyqasm/maps/decomposition_rules.py @@ -0,0 +1,133 @@ +# Copyright (C) 2024 qBraid +# +# This file is part of pyqasm +# +# Pyqasm is free software released under the GNU General Public License v3 +# or later. You can redistribute and/or modify it under the terms of the GPL v3. +# See the LICENSE file in the project root or . +# +# TH_summary_ERE IS NO WARRANTY for pyqasm, as per Section 15 of the GPL v3. + +""" +Definition of the decomposition rules for the gates in the basis sets. +""" +from enum import Enum +from pyqasm.elements import BasisSet +from pyqasm.maps.expressions import CONSTANTS_MAP + +class AppliedQubit(Enum): + """Enum to represent the qubits that are involved in the decomposition of a gate. + """ + QUBIT1 = 0 # Control qubit + QUBIT2 = 1 # Target qubit + +# Decomposition rules for the gates in the basis sets. +# The rules are defined as a dictionary where the key is the gate name. +# The value is a list of dictionaries. +# Each dictionary in the list represents a step in the decomposition of the gate. +# Each step is defined by a dictionary with the following keys: +# - gate: The name of the gate to be applied. +# - param: The parameter of the gate. +# - target_bit: The target qubit of the gate. +# This key is only used for decomposition of gates that operate on two qubits. +# - controll_bit: The control qubit of the gate. +# This key is only used for gates that operate on two qubits. +# +DECOMPOSITION_RULES = { + + BasisSet.ROTATIONAL_CX: { + "x": [ + {"gate": "rx", "param": CONSTANTS_MAP["pi"]} + ], + "y": [ + {"gate": "ry", "param": CONSTANTS_MAP["pi"]} + ], + "z": [ + {"gate": "rz", "param": CONSTANTS_MAP["pi"]} + ], + "h": [ + {"gate": "ry", "param": CONSTANTS_MAP["pi"]/2}, + {"gate": "rx", "param": CONSTANTS_MAP["pi"]} + ], + "s": [ + {"gate": "rz", "param": CONSTANTS_MAP["pi"]/2} + ], + "t": [ + {"gate": "rz", "param": CONSTANTS_MAP["pi"]/4} + ], + "sx": [ + {"gate": "rx", "param": CONSTANTS_MAP["pi"]/2} + ], + "sdg": [ + {"gate": "rz", "param": -CONSTANTS_MAP["pi"]/2} + ], + "tgd": [ + {"gate": "rz", "param": -CONSTANTS_MAP["pi"]/4} + ], + "cz": [ + {"gate": "ry", "param": CONSTANTS_MAP["pi"]/2, "target_bit": AppliedQubit.QUBIT2}, + {"gate": "rx", "param": CONSTANTS_MAP["pi"], "target_bit": AppliedQubit.QUBIT2}, + {"gate": "cx", "controll_bit": AppliedQubit.QUBIT1, "target_bit": AppliedQubit.QUBIT2}, + {"gate": "ry", "param": CONSTANTS_MAP["pi"]/2, "target_bit": AppliedQubit.QUBIT2}, + {"gate": "rx", "param": CONSTANTS_MAP["pi"], "target_bit": AppliedQubit.QUBIT2} + ], + "swap": [ + {"gate": "cx", "controll_bit": AppliedQubit.QUBIT1, "target_bit": AppliedQubit.QUBIT2}, + {"gate": "cx", "controll_bit": AppliedQubit.QUBIT2, "target_bit": AppliedQubit.QUBIT1}, + {"gate": "cx", "controll_bit": AppliedQubit.QUBIT1, "target_bit": AppliedQubit.QUBIT2} + ] + }, + BasisSet.CLIFFORD_T: { + "x": [ + {"gate": "h"}, + {"gate": "s"}, + {"gate": "s"}, + {"gate": "h"} + ], + "y": [ + {"gate": "s"}, + {"gate": "s"}, + {"gate": "h"}, + {"gate": "s"}, + {"gate": "s"}, + {"gate": "h"} + ], + "z": [ + {"gate": "s"}, + {"gate": "s"} + ], + "sx": [ + {"gate": "s"}, + {"gate": "s"}, + {"gate": "s"}, + {"gate": "h"}, + {"gate": "s"}, + {"gate": "s"}, + {"gate": "s"} + ], + "cz": [ + {"gate": "h", "target_bit": AppliedQubit.QUBIT2}, + {"gate": "cx", "controll_bit": AppliedQubit.QUBIT1, "target_bit": AppliedQubit.QUBIT2}, + {"gate": "h", "target_bit": AppliedQubit.QUBIT2}, + ], + "swap": [ + {"gate": "cx", "controll_bit": AppliedQubit.QUBIT1, "target_bit": AppliedQubit.QUBIT2}, + {"gate": "cx", "controll_bit": AppliedQubit.QUBIT2, "target_bit": AppliedQubit.QUBIT1}, + {"gate": "cx", "controll_bit": AppliedQubit.QUBIT1, "target_bit": AppliedQubit.QUBIT2} + ] + } +} + +# """TODO: Implement the Solovay-Kitaev algorithm""" +# +def solovay_kitaev_algo(gate_name, param, accuracy): + pass + + +def solovay_kitaev_helper(gate, depth): + if (depth == 0): + #return BasicApproximation(gate_name) + pass + else: + gate_stepDown = solovay_kitaev_helper(gate, depth - 1) + \ No newline at end of file diff --git a/src/pyqasm/maps/solovay-kitaev.py b/src/pyqasm/maps/solovay-kitaev.py new file mode 100644 index 0000000..7f43723 --- /dev/null +++ b/src/pyqasm/maps/solovay-kitaev.py @@ -0,0 +1,202 @@ +import numpy as np +from typing import List, Tuple + +class SU2Matrix: + """Class representing a 2x2 Special Unitary matrix.""" + def __init__(self, matrix: np.ndarray): + self.matrix = matrix + + def __mul__(self, other: 'SU2Matrix') -> 'SU2Matrix': + return SU2Matrix(np.dot(self.matrix, other.matrix)) + + def dagger(self) -> 'SU2Matrix': + """Returns the conjugate transpose.""" + return SU2Matrix(self.matrix.conj().T) + + def distance(self, other: 'SU2Matrix') -> float: + """Calculates the operator norm distance between two matrices.""" + diff = self.matrix - other.matrix + return np.linalg.norm(diff) + +def group_commutator(a: SU2Matrix, b: SU2Matrix) -> SU2Matrix: + """Compute the group commutator [a,b] = aba^{-1}b^{-1}.""" + return a * b * a.dagger() * b.dagger() + +def find_basic_approximation(target: SU2Matrix, basic_gates: List[SU2Matrix]) -> SU2Matrix: + """ + Find the best approximation of the target unitary from the basic gate set. + + Args: + target: Target unitary matrix to approximate + basic_gates: List of available basic gates + + Returns: + Best approximating unitary from the basic gate set + """ + best_dist = float('inf') + best_gate = None + + for gate in basic_gates: + dist = target.distance(gate) + if dist < best_dist: + best_dist = dist + best_gate = gate + + return best_gate + +def decompose_group_element(target: SU2Matrix, basic_gates: List[SU2Matrix], depth: int) -> Tuple[List[SU2Matrix], float]: + """ + Decompose a target unitary into a sequence of basic gates using Solovay-Kitaev. + + Args: + target: Target unitary matrix to decompose + basic_gates: List of available basic gates + depth: Recursion depth for the algorithm + + Returns: + Tuple of (sequence of gates, approximation error) + """ + if depth == 0: + best_approx = find_basic_approximation(target, basic_gates) + return [best_approx], target.distance(best_approx) + + # Recursive approximation + prev_sequence, prev_error = decompose_group_element(target, basic_gates, depth - 1) + + # If previous approximation is good enough, return it + # ERROR IS HARD CODED RIGHT NOW -> CHANGE THIS TO FIT USER-INPUT + if prev_error < 1e-6: + return prev_sequence, prev_error + + # Compute the error unitary + approx = prev_sequence[0] + for gate in prev_sequence[1:]: + approx = approx * gate + + error = target * approx.dagger() + + # Find Va and Vb such that their group commutator approximates the error + best_v = None + best_w = None + best_error = float('inf') + + for v in basic_gates: + for w in basic_gates: + comm = group_commutator(v, w) + curr_error = error.distance(comm) + if curr_error < best_error: + best_error = curr_error + best_v = v + best_w = w + + # Construct the final sequence + result_sequence = [] + result_sequence.extend(prev_sequence) + + # Add correction terms + if best_v is not None and best_w is not None: + v_sequence, error = decompose_group_element(best_v, basic_gates, depth - 1) + w_sequence, error = decompose_group_element(best_w, basic_gates, depth - 1) + + v_approx = v_sequence[0] + for gate in v_sequence[1:]: + v_approx = v_approx * gate + + w_approx = w_sequence[0] + for gate in w_sequence[1:]: + w_approx = w_approx * gate + + result_sequence.extend([best_v, best_w, v_approx.dagger(), w_approx.dagger()]) + + # Calculate final error + final_approx = result_sequence[0] + for gate in result_sequence[1:]: + final_approx = final_approx * gate + final_error = target.distance(final_approx) + + return result_sequence, final_error + +def solovay_kitaev(target: np.ndarray, basic_gates: List[np.ndarray], depth: int = 3) -> List[np.ndarray]: + """ + Main function to run the Solovay-Kitaev algorithm. + + Args: + target: Target unitary matrix as numpy array + basic_gates: List of basic gates as numpy arrays + depth: Recursion depth + + Returns: + List of gates that approximate the target unitary + """ + # Convert inputs to SU2Matrix objects + target_su2 = SU2Matrix(target) + basic_gates_su2 = [SU2Matrix(gate) for gate in basic_gates] + + # Run the decomposition + sequence, error = decompose_group_element(target_su2, basic_gates_su2, depth) + + # Convert back to numpy arrays + return [gate.matrix for gate in sequence] + +def gate_to_name(gate: np.ndarray) -> str: + """Convert a gate (numpy array) to its standard name.""" + # Helper function to check if matrices are equal within numerical precision + # I just made up a tolerance, does it matter too much? I thought this would make my life easier + def is_equal(A, B, tol=1e-10): + return np.allclose(A, B, rtol=tol, atol=tol) + + # Define gates and their names + gate_map = { + 'X': X, + 'Y': Y, + 'Z': Z, + 'H': H, + 'S': S, + 'S_dagger': S.conj().T, + 'T': T, + 'T_dagger': T.conj().T, + } + + # Check each known gate + for name, reference_gate in gate_map.items(): + if is_equal(gate, reference_gate): + return name + + # If no match found, return a generic descriptor + return "U" # Unknown gate + +# FOR TESTING - DELETE BEFORE MERGING +if __name__ == "__main__": + X = np.array([[0, 1], [1, 0]], dtype=complex) + Y = np.array([[0, -1j], [1j, 0]], dtype=complex) + Z = np.array([[1, 0], [0, -1]], dtype=complex) + H = np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2) + S = np.array([[1, 0], [0, 1j]], dtype=complex) + T = np.array([[1, 0], [0, np.exp((1j*np.pi)/4)]], dtype=complex) + + basis_gates = [H, T, S, T.conj().T, S.conj().T] + + # Define a target unitary + theta = np.pi / 4 + target = np.array([ + [np.cos(theta/2), -np.sin(theta/2)], + [np.sin(theta/2), np.cos(theta/2)] + ], dtype=complex) + + # Run the algorithm + sequence = solovay_kitaev(target, basis_gates, depth=3) + + sequence_to_string = "" + + for matrix in sequence: + sequence_to_string += gate_to_name(matrix) + " " + + # Compare target matrix and approx matrix + print(target) + + approx = sequence[0] + for gate in sequence[1:]: + approx = approx * gate + print(approx) + + print(sequence_to_string) \ No newline at end of file From 4a5827bb04ae671d7a2cde1cd2da1c232a9c8503 Mon Sep 17 00:00:00 2001 From: Ashwin Kaliyaperumal Date: Wed, 12 Feb 2025 17:13:42 -0800 Subject: [PATCH 3/3] updated solovay kitaev.py using updated code with proper testing scaffolding (prints out approx and target array) --- src/solovay_kitaev.py | 104 ++++++++++++++++++++++++++++++------------ 1 file changed, 74 insertions(+), 30 deletions(-) diff --git a/src/solovay_kitaev.py b/src/solovay_kitaev.py index 73cc6f3..f9f4c93 100644 --- a/src/solovay_kitaev.py +++ b/src/solovay_kitaev.py @@ -44,11 +44,7 @@ def find_basic_approximation(target: SU2Matrix, basic_gates: List[SU2Matrix]) -> return best_gate -def decompose_group_element( - target: SU2Matrix, - basic_gates: List[SU2Matrix], - depth: int -) -> Tuple[List[SU2Matrix], float]: +def decompose_group_element(target: SU2Matrix, basic_gates: List[SU2Matrix], depth: int) -> Tuple[List[SU2Matrix], float]: """ Decompose a target unitary into a sequence of basic gates using Solovay-Kitaev. @@ -68,6 +64,7 @@ def decompose_group_element( prev_sequence, prev_error = decompose_group_element(target, basic_gates, depth - 1) # If previous approximation is good enough, return it + # ERROR IS HARD CODED RIGHT NOW -> CHANGE THIS TO FIT USER-INPUT if prev_error < 1e-6: return prev_sequence, prev_error @@ -79,27 +76,37 @@ def decompose_group_element( error = target * approx.dagger() # Find Va and Vb such that their group commutator approximates the error - # This is a simplified version - in practice, you'd need a more sophisticated search - best_va = None - best_vb = None + best_v = None + best_w = None best_error = float('inf') - for va in basic_gates: - for vb in basic_gates: - comm = group_commutator(va, vb) + for v in basic_gates: + for w in basic_gates: + comm = group_commutator(v, w) curr_error = error.distance(comm) if curr_error < best_error: best_error = curr_error - best_va = va - best_vb = vb - + best_v = v + best_w = w + # Construct the final sequence result_sequence = [] result_sequence.extend(prev_sequence) # Add correction terms - if best_va is not None and best_vb is not None: - result_sequence.extend([best_va, best_vb, best_va.dagger(), best_vb.dagger()]) + if best_v is not None and best_w is not None: + v_sequence, error = decompose_group_element(best_v, basic_gates, depth - 1) + w_sequence, error = decompose_group_element(best_w, basic_gates, depth - 1) + + v_approx = v_sequence[0] + for gate in v_sequence[1:]: + v_approx = v_approx * gate + + w_approx = w_sequence[0] + for gate in w_sequence[1:]: + w_approx = w_approx * gate + + result_sequence.extend([best_v, best_w, v_approx.dagger(), w_approx.dagger()]) # Calculate final error final_approx = result_sequence[0] @@ -109,11 +116,7 @@ def decompose_group_element( return result_sequence, final_error -def solovay_kitaev( - target: np.ndarray, - basic_gates: List[np.ndarray], - depth: int = 3 -) -> List[np.ndarray]: +def solovay_kitaev(target: np.ndarray, basic_gates: List[np.ndarray], depth: int = 3) -> List[np.ndarray]: """ Main function to run the Solovay-Kitaev algorithm. @@ -135,24 +138,65 @@ def solovay_kitaev( # Convert back to numpy arrays return [gate.matrix for gate in sequence] -# Example usage: +def gate_to_name(gate: np.ndarray) -> str: + """Convert a gate (numpy array) to its standard name.""" + # Helper function to check if matrices are equal within numerical precision + # I just made up a tolerance, does it matter too much? I thought this would make my life easier + def is_equal(A, B, tol=1e-10): + return np.allclose(A, B, rtol=tol, atol=tol) + + # Define gates and their names + gate_map = { + 'X': X, + 'Y': Y, + 'Z': Z, + 'H': H, + 'S': S, + 'S_dagger': S.conj().T, + 'T': T, + 'T_dagger': T.conj().T, + } + + # Check each known gate + for name, reference_gate in gate_map.items(): + if is_equal(gate, reference_gate): + return name + + # If no match found, return a generic descriptor + return "U" # Unknown gate + +# FOR TESTING - DELETE BEFORE MERGING if __name__ == "__main__": - # Define some basic gates (Pauli matrices and their combinations) - I = np.array([[1, 0], [0, 1]], dtype=complex) X = np.array([[0, 1], [1, 0]], dtype=complex) Y = np.array([[0, -1j], [1j, 0]], dtype=complex) Z = np.array([[1, 0], [0, -1]], dtype=complex) H = np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2) + S = np.array([[1, 0], [0, 1j]], dtype=complex) + T = np.array([[1, 0], [0, np.exp((1j*np.pi)/4)]], dtype=complex) - basic_gates = [I, X, Y, Z, H] + basis_gates = [H, T, S, T.conj().T, S.conj().T] # Define a target unitary - theta = np.pi / 8 + theta = np.pi / 4 target = np.array([ - [np.cos(theta), -np.sin(theta)], - [np.sin(theta), np.cos(theta)] + [np.cos(theta/2), -np.sin(theta/2)], + [np.sin(theta/2), np.cos(theta/2)] ], dtype=complex) # Run the algorithm - sequence = solovay_kitaev(target, basic_gates, depth=3) - print(f"Found sequence of {len(sequence)} gates") \ No newline at end of file + sequence = solovay_kitaev(target, basis_gates, depth=3) + + sequence_to_string = "" + + for matrix in sequence: + sequence_to_string += gate_to_name(matrix) + " " + + # Compare target matrix and approx matrix + print(target) + + approx = sequence[0] + for gate in sequence[1:]: + approx = approx * gate + print(approx) + + print(sequence_to_string)