Implemented Solovay Kitaev's Algorithm

src/pyqasm/algorithms/__init__.py

@@ -0,0 +1,14 @@
+"""
+Sub module for quantum algorithms.
+
+Functions:
+----------
+    solovay_kitaev: Solovay-Kitaev algorithm for approximating unitary gates.
+
+"""
+
+from pyqasm.algorithms.solovay_kitaev.solovay_kitaev import solovay_kitaev
+
+__all__ = [
+    "solovay_kitaev",
+]

src/pyqasm/algorithms/solovay_kitaev/basic_approximation.py

@@ -0,0 +1,98 @@
+"""
+Definition of the basic approximation algorithm for the Solovay-Kitaev theorem.
+"""
+
+import numpy as np
+
+from pyqasm.algorithms.solovay_kitaev.utils import TU2Matrix, get_tu2matrix_for_basic_approximation
+from pyqasm.elements import BasisSet
+
+
+def rescursive_traversal(
+    target_matrix, approximated_matrix, target_gate_set_list, current_depth, params
+):
+    """Recursively traverse the tree to find the best approximation of the target matrix.
+    Args:
+        target_matrix (np.ndarray): The target matrix to approximate.
+        approximated_matrix (TU2Matrix): The approximated matrix.
+        target_gate_set_list (list): The list of target gates to approximate.
+        current_depth (int): The current depth of the tree.
+        params (tuple): The parameters for the approximation.
+
+    Returns:
+        float: The closest difference between the target and approximated matrix.
+        TU2Matrix: The closest approximated matrix.
+        TU2Matrix: The best approx
+
+    """
+    accuracy, max_tree_depth, best_gate = params
+
+    if current_depth >= max_tree_depth:
+        return best_gate
+
+    for gate in target_gate_set_list:
+        if not approximated_matrix.can_multiple(gate):
+            continue
+        approximated_matrix_copy = approximated_matrix.copy()
+        approximated_matrix = approximated_matrix * gate
+
+        diff = approximated_matrix.distance(target_matrix)
+        if diff < accuracy:
+            best_gate = approximated_matrix.copy()
+            return best_gate
+
+        # Update the closest gate if the current one is closer
+        if diff < best_gate.distance(target_matrix):
+            best_gate = approximated_matrix.copy()
+
+        best_gate = rescursive_traversal(
+            target_matrix,
+            approximated_matrix.copy(),
+            target_gate_set_list,
+            current_depth + 1,
+            (accuracy, max_tree_depth, best_gate),
+        )
+        approximated_matrix = approximated_matrix_copy.copy()
+
+    return best_gate
+
+
+def basic_approximation(target_matrix, target_gate_set, accuracy=0.001, max_tree_depth=3):
+    """Approximate the target matrix using the basic approximation algorithm.
+
+    Args:
+        target_matrix (np.ndarray): The target matrix to approximate.
+        target_gate_set (BasisSet): The target gate set to approximate.
+        accuracy (float): The accuracy of the approximation.
+        max_tree_depth (int): The maximum depth of the tree.
+
+    Returns:
+        TU2Matrix: The approximated matrix.
+    """
+    approximated_matrix = TU2Matrix(np.identity(2), [], None, None)
+    target_gate_set_list = get_tu2matrix_for_basic_approximation(target_gate_set)
+    current_depth = 0
+    best_gate = TU2Matrix(np.identity(2), [], None, None)
+
+    params = (accuracy, max_tree_depth, best_gate)
+
+    best_gate = rescursive_traversal(
+        target_matrix, approximated_matrix.copy(), target_gate_set_list, current_depth, params
+    )
+
+    return best_gate
+
+    # result = None
+
+    # if best_gate:
+    #     result = best_gate.copy()
+    # else:
+    #     result = closest_gate.copy()
+
+    # return result
+
+
+if __name__ == "__main__":
+    target_matrix_U = np.array([[0.70711, 0.70711j], [0.70711j, 0.70711]])
+
+    print(basic_approximation(target_matrix_U, BasisSet.CLIFFORD_T, 0.0001, 3))

src/pyqasm/algorithms/solovay_kitaev/optimizer.py

@@ -0,0 +1,94 @@
+"""
+Definition of the optimizer for the Solovay-Kitaev theorem.
+"""
+
+from pyqasm.algorithms.solovay_kitaev.utils import get_identity_weight_group_for_optimizer
+from pyqasm.elements import BasisSet
+
+
+def optimize_gate_sequence(seq: list[str], target_basis_set):
+    """Optimize a gate sequence based on the identity weight group.
+    Args:
+        seq (list[str]): The gate sequence to optimize.
+        target_basis_set (BasisSet): The target basis set.
+
+    Returns:
+        list[str]: The optimized gate sequence.
+    """
+    target_identity_weight_group = get_identity_weight_group_for_optimizer(target_basis_set)
+    for _ in range(int(1e6)):
+        current_group = None
+        current_weight = 0
+        start_index = 0
+        changed = False
+
+        for i, gate_name in enumerate(seq):
+            if gate_name not in target_identity_weight_group:
+                continue
+
+            gate = target_identity_weight_group[gate_name]
+            new_group = gate["group"]
+            new_weight = gate["weight"]
+
+            if current_group is None or new_group != current_group:
+                current_group = new_group
+                current_weight = new_weight
+                start_index = i
+            else:
+                current_weight += new_weight
+
+            if current_weight == 1:
+                seq = seq[:start_index] + seq[i + 1 :]
+                changed = True
+                break
+            if current_weight > 1:
+                remaining_weight = current_weight - 1
+                for key, value in target_identity_weight_group.items():
+                    if value["group"] == current_group and value["weight"] == remaining_weight:
+                        seq = seq[:start_index] + [key] + seq[i + 1 :]
+                        changed = True
+                        break
+                break
+
+        if not changed:
+            return seq
+
+    return seq
+
+
+if __name__ == "__main__":
+    s1 = ["s", "s", "s", "t", "t", "tdg", "sdg", "sdg", "sdg", "tdg", "s", "h", "s"]
+    s2 = [
+        "t",
+        "s",
+        "s",
+        "s",
+        "t",
+        "tdg",
+        "tdg",
+        "sdg",
+        "sdg",
+        "sdg",
+        "t",
+        "s",
+        "s",
+        "s",
+        "t",
+        "tdg",
+        "tdg",
+        "sdg",
+        "sdg",
+        "sdg",
+        "s",
+        "h",
+        "s",
+    ]
+    s3 = ["h", "s", "s", "t", "t", "s", "t"]  # ['h', 't']
+    s4 = ["h", "s", "s", "t", "t", "s", "h"]  # []
+    s5 = ["h", "s", "s", "t", "h", "h", "t", "s", "h", "t"]  # ['t']
+
+    print(optimize_gate_sequence(s1, BasisSet.CLIFFORD_T) == ["s", "h", "s"])
+    print(optimize_gate_sequence(s2, BasisSet.CLIFFORD_T) == ["s", "h", "s"])
+    print(optimize_gate_sequence(s3, BasisSet.CLIFFORD_T) == ["h", "t"])
+    print(optimize_gate_sequence(s4, BasisSet.CLIFFORD_T) == [])
+    print(optimize_gate_sequence(s5, BasisSet.CLIFFORD_T) == ["t"])

src/pyqasm/algorithms/solovay_kitaev/solovay_kitaev.py

@@ -0,0 +1,169 @@
+"""
+Definition of the Solovay-Kitaev algorithm.
+"""
+
+from typing import List, Tuple
+
+import numpy as np
+
+from pyqasm.algorithms.solovay_kitaev.basic_approximation import basic_approximation
+from pyqasm.algorithms.solovay_kitaev.optimizer import optimize_gate_sequence
+from pyqasm.algorithms.solovay_kitaev.utils import (
+    SU2Matrix,
+    get_su2matrix_for_solovay_kitaev_algorithm,
+)
+from pyqasm.elements import BasisSet
+
+
+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_matrix: SU2Matrix, target_basis_set, use_optimization, accuracy=1e-6, max_tree_depth=10
+) -> SU2Matrix:
+    """Find the basic approximation of a target matrix.
+
+    Args:
+        target_matrix (SU2Matrix): The target matrix to approximate.
+        target_basis_set (BasisSet): The basis set to use for the approximation.
+        use_optimization (bool): Whether to use optimization to reduce the number of gates.
+        accuracy (float): The accuracy of the approximation.
+        max_tree_depth (int): The maximum depth of the tree.
+
+    Returns:
+        SU2Matrix: The basic approximation of the target matrix.
+    """
+    gates = basic_approximation(target_matrix, target_basis_set, accuracy, max_tree_depth)
+    if use_optimization:
+        gates.name = optimize_gate_sequence(gates.name, target_basis_set)
+    return SU2Matrix(gates.matrix, gates.name)
+
+
+def decompose_group_element(
+    target: SU2Matrix,
+    target_gate_set,
+    basic_gates: List[SU2Matrix],
+    depth,
+    accuracy,
+    use_optimization,
+) -> Tuple[List[SU2Matrix], float]:
+    """Decompose a group element into a sequence of basic gates.
+
+    Args:
+        target (SU2Matrix): The target group element.
+        target_gate_set (BasisSet): The target gate set.
+        basic_gates (List[SU2Matrix]): The basic gates to use for the approximation.
+        depth (int): The depth of the approximation.
+        accuracy (float): The accuracy of the approximation.
+        use_optimization (bool): Whether to use optimization to reduce the number of gates.
+
+    Returns:
+        Tuple[List[SU2Matrix], float]: The sequence of basic gates and the error.
+    """
+
+    if depth == 0:
+        best_approx = find_basic_approximation(
+            target.matrix, target_gate_set, use_optimization=use_optimization
+        )
+        return best_approx, target.distance(best_approx)
+
+    # Recursive approximation
+    prev_sequence, prev_error = decompose_group_element(
+        target, target_gate_set, basic_gates, depth - 1, accuracy, use_optimization
+    )
+
+    # If previous approximation is good enough, return it
+    # ERROR IS HARD CODED RIGHT NOW -> CHANGE THIS TO FIT USER-INPUT
+    if prev_error < accuracy:
+        return prev_sequence, prev_error
+
+    error = target * prev_sequence.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
+
+    result = 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, target_gate_set, basic_gates, depth - 1, accuracy, use_optimization
+        )
+        w_sequence, error = decompose_group_element(
+            best_w, target_gate_set, basic_gates, depth - 1, accuracy, use_optimization
+        )
+
+        result = group_commutator(v_sequence, w_sequence) * prev_sequence
+
+    final_error = target.distance(result)
+
+    return result, final_error
+
+
+def solovay_kitaev(
+    target: np.ndarray, target_basis_set, depth: int = 3, accuracy=1e-6, use_optimization=True
+) -> List[np.ndarray]:
+    """
+    Main function to run the Solovay-Kitaev algorithm.
+
+    Args:
+        target: The target unitary matrix to approximate
+        target_basis_set: The basis set to use for the approximation
+        depth: The depth of the approximation
+        accuracy: The accuracy of the approximation
+        use_optimization: Whether to use optimization to reduce the number
+
+    Returns:
+        A list of gates that approximate the target unitary matrix
+    """
+    # Convert inputs to SU2Matrix objects
+    target_su2 = SU2Matrix(target, [])
+
+    basic_gates_su2 = get_su2matrix_for_solovay_kitaev_algorithm(target_basis_set)
+
+    # Run the decomposition
+    sequence, _ = decompose_group_element(
+        target_su2, target_basis_set, basic_gates_su2, depth, accuracy, use_optimization
+    )
+
+    if use_optimization:
+        sequence.name = optimize_gate_sequence(sequence.name, target_basis_set)
+        return sequence
+
+    return sequence
+
+
+if __name__ == "__main__":
+    target_matrix_U = np.array([[0.70711, 0.70711j], [0.70711j, 0.70711]])
+
+    r0 = solovay_kitaev(target_matrix_U, BasisSet.CLIFFORD_T, depth=0)
+    print(r0.name)  # Output: ['s', 'h', 's']
+
+    r1 = solovay_kitaev(target_matrix_U, BasisSet.CLIFFORD_T, depth=1)
+    print(
+        r1.name
+    )  # Output: ['s', 's', 's', 't', 't', 'tdg', 'sdg', 'sdg', 'sdg', 'tdg', 's', 'h', 's']
+
+    r2 = solovay_kitaev(target_matrix_U, BasisSet.CLIFFORD_T, depth=2)
+    print(r2.name)  # Output: ['t', 's', 's', 's', 't',
+    #             'tdg', 'tdg', 'sdg', 'sdg', 'sdg',
+    #             't', 's', 's', 's', 't',
+    #             'tdg', 'tdg', 'sdg', 'sdg', 'sdg',
+    #             's', 'h', 's']
+
+    print(np.allclose(r0.matrix, r1.matrix))  # Output: True
+    print(np.allclose(r1.matrix, r2.matrix))  # Output: True
+    print(np.allclose(r2.matrix, r0.matrix))  # Output: True