Implemented Solovay Kitaev's Algorithm

src/pyqasm/algorithms/solovay_kitaev/basic_approximation.py

@@ -0,0 +1,52 @@
+from math import pi
+import pickle
+import numpy as np
+
+from pyqasm.elements import BasisSet
+import os
+
+
+def basic_approximation(U, target_gate_set, accuracy=0.001, max_tree_depth=10):
+    current_dir = os.path.dirname(os.path.abspath(__file__))
+    gate_set_files = {
+        BasisSet.CLIFFORD_T: os.path.join(current_dir, "cache", "clifford-t_depth-5.pkl"),
+    }
+
+    if target_gate_set not in gate_set_files:
+        raise ValueError(f"Unknown target gate set: {target_gate_set}")
+
+    pkl_file_name = gate_set_files[target_gate_set]
+    try:
+        with open(pkl_file_name, "rb") as file:
+            gate_list = pickle.load(file)
+    except FileNotFoundError:
+        raise FileNotFoundError(f"Pickle file not found: {pkl_file_name}")
+
+    closest_gate = None
+    closest_trace_diff = float("inf")
+
+    for gate in gate_list:
+        gate_matrix = gate["matrix"]
+        tree_depth = gate["depth"]
+
+        # Stop if the maximum depth is exceeded
+        if tree_depth > max_tree_depth:
+            break
+
+        trace_diff = np.abs(np.trace(np.dot(gate_matrix.conj().T, U) - np.identity(2)))
+
+        if trace_diff < accuracy:
+            return gate
+
+        # Update the closest gate if the current one is closer
+        if trace_diff < closest_trace_diff:
+            closest_trace_diff = trace_diff
+            closest_gate = gate
+
+    return closest_gate
+
+
+if __name__ == "__main__":
+    U = np.array([[0.70711,  0.70711j],
+                  [0.70711j, 0.70711]])
+    print(basic_approximation(U, BasisSet.CLIFFORD_T, 0.001, 10))

src/pyqasm/algorithms/solovay_kitaev/generator.py

@@ -0,0 +1,102 @@
+import pickle
+import sys
+from collections import deque
+
+import numpy as np
+
+gate_sets = {
+    "clifford_T": [
+        {
+            "name": "h", 
+            "identity": {
+                "group": "h",
+                "weight": 0.5
+            },
+            "matrix": (1 / np.sqrt(2)) * np.array([[1, 1], [1, -1]])
+        },
+        {
+            "name": "s",
+            "identity": {
+                "group": "s-t",
+                "weight": 0.25
+            },
+            "matrix": np.array([[1, 0], [0, 1j]])
+        },
+        {
+            "name": "t",
+            "identity": {
+                "group": "s-t",
+                "weight": 0.125
+            },
+            "matrix": np.array([[1, 0], [0, np.exp(1j * np.pi / 4)]]),
+        },
+    ]
+}
+
+
+def generate_solovay_kitaev_tree_cache(target_gate_set, max_depth, pkl_file_name):
+    queue = deque([{"name": [], "depth": 0, "matrix": np.eye(2), "identity": {"group": None, "weight": 0}}])
+    result = []
+
+    while queue:
+        node = queue.popleft()
+        if node["depth"] == max_depth:
+            break
+
+        for gate in target_gate_set:
+            new_group = gate["identity"]["group"]
+            new_weight = gate["identity"]["weight"]
+            current_group = node["identity"]["group"]
+            current_weight = node["identity"]["weight"]
+
+            if current_group != new_group:
+                new_node = {
+                    "name": node["name"] + [gate["name"]],
+                    "depth": node["depth"] + 1,
+                    "matrix": np.dot(node["matrix"], gate["matrix"]),
+                    "identity": {"group": new_group, "weight": new_weight}
+                }
+                queue.append(new_node)
+                result.append(new_node)
+            elif current_weight + new_weight < 1:
+                new_node = {
+                    "name": node["name"] + [gate["name"]],
+                    "depth": node["depth"] + 1,
+                    "matrix": np.dot(node["matrix"], gate["matrix"]),
+                    "identity": {"group": current_group, "weight": current_weight + new_weight}
+                }
+                queue.append(new_node)
+                result.append(new_node)
+
+    print(result)
+
+    with open("cache/" + pkl_file_name, "wb") as f:
+        pickle.dump(result, f)
+
+
+if __name__ == "__main__":
+    """
+    How to use:
+
+    Run this file direct, and pass the following command line arguments:
+
+    target_gate_set: The target basis set of which you want to generate cached tree.
+    max_depth: Max depth of the tree which you want to cache.
+    pkl_file_name: Name of the pickel file in with you want to save the generated cache tree.
+
+    Your command will look like this:
+    python generator.py <target_gate_set> <max_depth> <pkl_file_name>
+    eg.:
+    python generator.py clifford_T 10 <pkl_file_name> clifford-t_depth-10.pkl
+
+    The file will be saved in cache dir.
+    """
+
+
+    target_gate_set = sys.argv[1]
+    max_depth = sys.argv[2]
+    pkl_file_name = sys.argv[3]
+
+    t = generate_solovay_kitaev_tree_cache(
+        gate_sets[target_gate_set], int(max_depth), pkl_file_name
+    )

src/pyqasm/algorithms/solovay_kitaev/optimizer.py

@@ -0,0 +1,76 @@
+import numpy as np
+from pyqasm.elements import BasisSet
+
+IDENTITY_WEIGHT_GROUP = {
+    BasisSet.CLIFFORD_T: {
+        "h": {
+            "group": "h",
+            "weight": 0.5
+        },
+        "s": {
+            "group": "s-t",
+            "weight": 0.25
+        },
+        "t": {
+            "group": "s-t",
+            "weight": 0.125
+        },
+        "sdg": {
+            "group": "sdg-tdg",
+            "weight": 0.25
+        },
+        "tdg": {
+            "group": "sdg-tdg",
+            "weight": 0.125
+        },
+    }
+}
+
+def optimize_gate_sequnce(seq: list[str], target_basis_set):
+    target_identity_weight_group = IDENTITY_WEIGHT_GROUP[target_basis_set]
+    while True:
+        current_group = None
+        current_weight = 0
+        start_index = 0
+        changed = False
+
+        for i, gate_name in enumerate(seq):
+            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
+            elif 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    
+
+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_sequnce(s1, BasisSet.CLIFFORD_T) == ['s', 'h', 's'])
+    print(optimize_gate_sequnce(s2, BasisSet.CLIFFORD_T) == ['s', 'h', 's'])
+    print(optimize_gate_sequnce(s3, BasisSet.CLIFFORD_T) == ['h', 't'])
+    print(optimize_gate_sequnce(s4, BasisSet.CLIFFORD_T) == [])
+    print(optimize_gate_sequnce(s5, BasisSet.CLIFFORD_T) == ['t'])

src/pyqasm/algorithms/solovay_kitaev/solovay_kitaev.py

@@ -0,0 +1,147 @@
+import numpy as np
+from typing import List, Tuple
+
+from pyqasm.algorithms.solovay_kitaev.generator import gate_sets
+from pyqasm.algorithms.solovay_kitaev.optimizer import optimize_gate_sequnce
+from pyqasm.maps.gates import BASIS_GATE_MAP, SELF_INVERTING_ONE_QUBIT_OP_SET, ST_GATE_INV_MAP
+from pyqasm.algorithms.solovay_kitaev.basic_approximation import basic_approximation
+from pyqasm.elements import BasisSet
+
+gate_matrix = {
+    "h": (1 / np.sqrt(2)) * np.array([[1, 1], [1, -1]]),
+    "s": np.array([[1, 0], [0, 1j]]),
+    "t": np.array([[1, 0], [0, np.exp(1j * np.pi / 4)]]),
+    "sdg": np.array([[1, 0], [0, 1j]]).conj().T,
+    "tdg": np.array([[1, 0], [0, np.exp(1j * np.pi / 4)]]).conj().T,
+}
+
+class SU2Matrix:
+    """Class representing a 2x2 Special Unitary matrix."""
+    def __init__(self, matrix: np.ndarray, name: List[str]):
+        self.matrix = matrix
+        self.name = name
+
+    def __mul__(self, other: 'SU2Matrix') -> 'SU2Matrix':
+        matrix = np.dot(self.matrix, other.matrix)
+        name = self.name.copy()
+        name.extend(other.name)
+        return SU2Matrix(matrix, name)
+
+    def dagger(self) -> 'SU2Matrix':
+        """Returns the conjugate transpose."""
+        matrix = self.matrix.conj().T
+        name = []
+        for n in self.name[::-1]:
+            name.append(self._get_dagger_gate_name(n))
+
+        return SU2Matrix(matrix, name)
+
+    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 _get_dagger_gate_name(self,name: str):
+        if name in SELF_INVERTING_ONE_QUBIT_OP_SET:
+            return name
+        else:
+            return ST_GATE_INV_MAP[name]
+
+    def __str__(self):
+        return f"name: {self.name}, matrix: {self.matrix}"
+
+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(U: SU2Matrix, target_basis_set, accuracy=0.001, max_tree_depth=10) -> SU2Matrix:
+    gates = basic_approximation(U, target_basis_set, accuracy, max_tree_depth)
+    return SU2Matrix(gates["matrix"], gates["name"])
+
+def decompose_group_element(target: SU2Matrix, target_gate_set, basic_gates: List[SU2Matrix], depth: int) -> Tuple[List[SU2Matrix], float]:
+
+    if depth == 0:
+        best_approx = find_basic_approximation(target.matrix, target_gate_set)
+        return best_approx, target.distance(best_approx)
+
+    # Recursive approximation
+    prev_sequence, prev_error = decompose_group_element(target, target_gate_set, 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
+
+    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)
+        w_sequence, error = decompose_group_element(best_w, target_gate_set, basic_gates, depth - 1)
+
+        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) -> List[np.ndarray]:
+    """
+    Main function to run the Solovay-Kitaev algorithm.
+
+    Args:
+        target: Target unitary matrix as numpy array
+        target_basis_set: The target basis set to rebase the module to.
+        depth: Recursion depth
+
+    Returns:
+        List of gates that approximate the target unitary
+    """
+    # Convert inputs to SU2Matrix objects
+    target_su2 = SU2Matrix(target, [])
+
+    target_basis_gate_list = BASIS_GATE_MAP[target_basis_set]
+    basic_gates_su2 = [SU2Matrix(gate_matrix[gate], [gate]) for gate in target_basis_gate_list if gate != "cx"]
+
+
+    # Run the decomposition
+    sequence, error = decompose_group_element(target_su2, target_basis_set, basic_gates_su2, depth)
+
+    return sequence
+    # return optimize_gate_sequnce(sequence, target_basis_set)
+
+if __name__ == '__main__':  
+    U = np.array([[0.70711,  0.70711j],
+                  [0.70711j, 0.70711]])
+
+    r0 = solovay_kitaev(U, BasisSet.CLIFFORD_T, depth=0)
+    print(r0.name)  # Output: ['s', 'h', 's']
+
+    r1 = solovay_kitaev(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(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
+
+    # Test optimizer
+    print(optimize_gate_sequnce(r2.name, BasisSet.CLIFFORD_T))