diff --git a/qiskit/synthesis/__init__.py b/qiskit/synthesis/__init__.py
index ff1cd4f123ec..ceed2fb8135e 100644
--- a/qiskit/synthesis/__init__.py
+++ b/qiskit/synthesis/__init__.py
@@ -15,6 +15,8 @@
 Circuit Synthesis (:mod:`qiskit.synthesis`)
 ===========================================
 
+.. automodule:: qiskit.synthesis.discrete_basis
+
 .. currentmodule:: qiskit.synthesis
 
 Evolution Synthesis
diff --git a/qiskit/synthesis/discrete_basis/__init__.py b/qiskit/synthesis/discrete_basis/__init__.py
new file mode 100644
index 000000000000..6f5e21bf22e2
--- /dev/null
+++ b/qiskit/synthesis/discrete_basis/__init__.py
@@ -0,0 +1,91 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 2022.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+r"""
+=======================================================================================
+Approximately decompose 1q gates to a discrete basis using the Solovay-Kitaev algorithm
+=======================================================================================
+
+.. currentmodule:: qiskit.synthesis.discrete_basis
+
+Approximately decompose 1q gates to a discrete basis using the Solovay-Kitaev algorithm.
+
+The Solovay-Kitaev theorem [1] states that any single qubit gate can be approximated to
+arbitrary precision by a set of fixed single-qubit gates, if the set generates a dense
+subset in :math:`SU(2)`. This is an important result, since it means that any single-qubit
+gate can be expressed in terms of a discrete, universal gate set that we know how to implement
+fault-tolerantly. Therefore, the Solovay-Kitaev algorithm allows us to take any
+non-fault tolerant circuit and rephrase it in a fault-tolerant manner.
+
+This implementation of the Solovay-Kitaev algorithm is based on [2].
+
+For example, the following circuit
+
+.. parsed-literal::
+
+         ┌─────────┐
+    q_0: ┤ RX(0.8) ├
+         └─────────┘
+
+can be decomposed into
+
+.. parsed-literal::
+
+    global phase: 7π/8
+         ┌───┐┌───┐┌───┐
+    q_0: ┤ H ├┤ T ├┤ H ├
+         └───┘└───┘└───┘
+
+with an L2-error of approximately 0.01.
+
+
+Examples:
+
+    .. jupyter-execute::
+
+        import numpy as np
+        from qiskit.circuit import QuantumCircuit
+        from qiskit.circuit.library import TGate, HGate, TdgGate
+        from qiskit.converters import circuit_to_dag, dag_to_circuit
+        from qiskit.transpiler.passes.synthesis import SolovayKitaev
+        from qiskit.quantum_info import Operator
+
+        circuit = QuantumCircuit(1)
+        circuit.rx(0.8, 0)
+        dag = circuit_to_dag(circuit)
+
+        print('Original circuit:')
+        print(circuit.draw())
+
+        basis_gates = [TGate(), TdgGate(), HGate()]
+        skd = SolovayKitaev(recursion_degree=2)
+
+        discretized = dag_to_circuit(skd.run(dag))
+
+        print('Discretized circuit:')
+        print(discretized.draw())
+
+        print('Error:', np.linalg.norm(Operator(circuit).data - Operator(discretized).data))
+
+
+References:
+
+    [1]: Kitaev, A Yu (1997). Quantum computations: algorithms and error correction.
+         Russian Mathematical Surveys. 52 (6): 1191–1249.
+         `Online <https://iopscience.iop.org/article/10.1070/RM1997v052n06ABEH002155>`_.
+
+    [2]: Dawson, Christopher M.; Nielsen, Michael A. (2005) The Solovay-Kitaev Algorithm.
+         `arXiv:quant-ph/0505030 <https://arxiv.org/abs/quant-ph/0505030>`_
+
+"""
+
+from .solovay_kitaev import SolovayKitaevDecomposition
diff --git a/qiskit/synthesis/discrete_basis/commutator_decompose.py b/qiskit/synthesis/discrete_basis/commutator_decompose.py
new file mode 100644
index 000000000000..c50a19dfa35e
--- /dev/null
+++ b/qiskit/synthesis/discrete_basis/commutator_decompose.py
@@ -0,0 +1,245 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 2022.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Functions to compute the decomposition of an SO(3) matrix as balanced commutator."""
+
+from __future__ import annotations
+
+import math
+import numpy as np
+from .gate_sequence import _check_is_so3, GateSequence
+
+
+def _compute_trace_so3(matrix: np.ndarray) -> float:
+    """Computes trace of an SO(3)-matrix.
+
+    Args:
+        matrix: an SO(3)-matrix
+
+    Returns:
+        Trace of ``matrix``.
+
+    Raises:
+        ValueError: if ``matrix`` is not an SO(3)-matrix.
+    """
+    _check_is_so3(matrix)
+
+    trace = np.matrix.trace(matrix)
+    trace_rounded = min(trace, 3)
+    return trace_rounded
+
+
+def _compute_rotation_axis(matrix: np.ndarray) -> np.ndarray:
+    """Computes rotation axis of SO(3)-matrix.
+
+    Args:
+        matrix: The SO(3)-matrix for which rotation angle needs to be computed.
+
+    Returns:
+        The rotation axis of the SO(3)-matrix ``matrix``.
+
+    Raises:
+        ValueError: if ``matrix`` is not an SO(3)-matrix.
+    """
+    _check_is_so3(matrix)
+
+    trace = _compute_trace_so3(matrix)
+    theta = math.acos(0.5 * (trace - 1))
+    if math.sin(theta) > 1e-10:
+        x = 1 / (2 * math.sin(theta)) * (matrix[2][1] - matrix[1][2])
+        y = 1 / (2 * math.sin(theta)) * (matrix[0][2] - matrix[2][0])
+        z = 1 / (2 * math.sin(theta)) * (matrix[1][0] - matrix[0][1])
+    else:
+        x = 1.0
+        y = 0.0
+        z = 0.0
+    return np.array([x, y, z])
+
+
+def _solve_decomposition_angle(matrix: np.ndarray) -> float:
+    """Computes angle for balanced commutator of SO(3)-matrix ``matrix``.
+
+    Computes angle a so that the SO(3)-matrix ``matrix`` can be decomposed
+    as commutator [v,w] where v and w are both rotations of a about some axis.
+    The computation is done by solving a trigonometric equation using scipy.optimize.fsolve.
+
+    Args:
+        matrix: The SO(3)-matrix for which the decomposition angle needs to be computed.
+
+    Returns:
+        Angle a so that matrix = [v,w] with v and w rotations of a about some axis.
+
+    Raises:
+        ValueError: if ``matrix`` is not an SO(3)-matrix.
+    """
+    from scipy.optimize import fsolve
+
+    _check_is_so3(matrix)
+
+    trace = _compute_trace_so3(matrix)
+    angle = math.acos((1 / 2) * (trace - 1))
+
+    def objective(phi):
+        rhs = 2 * math.sin(phi / 2) ** 2
+        rhs *= math.sqrt(1 - math.sin(phi / 2) ** 4)
+        lhs = math.sin(angle / 2)
+        return rhs - lhs
+
+    decomposition_angle = fsolve(objective, angle)[0]
+    return decomposition_angle
+
+
+def _compute_rotation_from_angle_and_axis(  # pylint: disable=invalid-name
+    angle: float, axis: np.ndarray
+) -> np.ndarray:
+    """Computes the SO(3)-matrix corresponding to the rotation of ``angle`` about ``axis``.
+
+    Args:
+        angle: The angle of the rotation.
+        axis: The axis of the rotation.
+
+    Returns:
+        SO(3)-matrix that represents a rotation of ``angle`` about ``axis``.
+
+    Raises:
+        ValueError: if ``axis`` is not a 3-dim unit vector.
+    """
+    if axis.shape != (3,):
+        raise ValueError(f"Axis must be a 1d array of length 3, but has shape {axis.shape}.")
+
+    if abs(np.linalg.norm(axis) - 1.0) > 1e-4:
+        raise ValueError(f"Axis must have a norm of 1, but has {np.linalg.norm(axis)}.")
+
+    res = math.cos(angle) * np.identity(3) + math.sin(angle) * _cross_product_matrix(axis)
+    res += (1 - math.cos(angle)) * np.outer(axis, axis)
+    return res
+
+
+def _compute_rotation_between(from_vector: np.ndarray, to_vector: np.ndarray) -> np.ndarray:
+    """Computes the SO(3)-matrix for rotating ``from_vector`` to ``to_vector``.
+
+    Args:
+        from_vector: unit vector of size 3
+        to_vector: unit vector of size 3
+
+    Returns:
+        SO(3)-matrix that brings ``from_vector`` to ``to_vector``.
+
+    Raises:
+        ValueError: if at least one of ``from_vector`` of ``to_vector`` is not a 3-dim unit vector.
+    """
+    from_vector = from_vector / np.linalg.norm(from_vector)
+    to_vector = to_vector / np.linalg.norm(to_vector)
+
+    dot = np.dot(from_vector, to_vector)
+    cross = _cross_product_matrix(np.cross(from_vector, to_vector))
+    rotation_matrix = np.identity(3) + cross + np.dot(cross, cross) / (1 + dot)
+    return rotation_matrix
+
+
+def _cross_product_matrix(v: np.ndarray) -> np.ndarray:
+    """Computes cross product matrix from vector.
+
+    Args:
+        v: Vector for which cross product matrix needs to be computed.
+
+    Returns:
+        The cross product matrix corresponding to vector ``v``.
+    """
+    return np.array([[0, -v[2], v[1]], [v[2], 0, -v[0]], [-v[1], v[0], 0]])
+
+
+def _compute_commutator_so3(a: np.ndarray, b: np.ndarray) -> np.ndarray:
+    """Computes the commutator of the SO(3)-matrices ``a`` and ``b``.
+
+    The computation uses the fact that the inverse of an SO(3)-matrix is equal to its transpose.
+
+    Args:
+        a: SO(3)-matrix
+        b: SO(3)-matrix
+
+    Returns:
+        The commutator [a,b] of ``a`` and ``b`` w
+
+    Raises:
+        ValueError: if at least one of ``a`` or ``b`` is not an SO(3)-matrix.
+    """
+    _check_is_so3(a)
+    _check_is_so3(b)
+
+    a_dagger = np.conj(a).T
+    b_dagger = np.conj(b).T
+
+    return np.dot(np.dot(np.dot(a, b), a_dagger), b_dagger)
+
+
+def commutator_decompose(
+    u_so3: np.ndarray, check_input: bool = True
+) -> tuple[GateSequence, GateSequence]:
+    r"""Decompose an :math:`SO(3)`-matrix, :math:`U` as a balanced commutator.
+
+    This function finds two :math:`SO(3)` matrices :math:`V, W` such that the input matrix
+    equals
+
+    .. math::
+
+        U = V^\dagger W^\dagger V W.
+
+    For this decomposition, the following statement holds
+
+
+    .. math::
+
+        ||V - I||_F, ||W - I||_F \leq \frac{\sqrt{||U - I||_F}}{2},
+
+    where :math:`I` is the identity and :math:`||\cdot ||_F` is the Frobenius norm.
+
+    Args:
+        u_so3: SO(3)-matrix that needs to be decomposed as balanced commutator.
+        check_input: If True, checks whether the input matrix is actually SO(3).
+
+    Returns:
+        Tuple of GateSequences from SO(3)-matrices :math:`V, W`.
+
+    Raises:
+        ValueError: if ``u_so3`` is not an SO(3)-matrix.
+    """
+    if check_input:
+        # assert that the input matrix is really SO(3)
+        _check_is_so3(u_so3)
+
+        identity = np.identity(3)
+        if not (
+            np.allclose(u_so3.dot(u_so3.T), identity) and np.allclose(u_so3.T.dot(u_so3), identity)
+        ):
+            raise ValueError("Input matrix is not orthogonal.")
+
+    angle = _solve_decomposition_angle(u_so3)
+
+    # Compute rotation about x-axis with angle 'angle'
+    vx = _compute_rotation_from_angle_and_axis(angle, np.array([1, 0, 0]))
+
+    # Compute rotation about y-axis with angle 'angle'
+    wy = _compute_rotation_from_angle_and_axis(angle, np.array([0, 1, 0]))
+
+    commutator = _compute_commutator_so3(vx, wy)
+
+    u_so3_axis = _compute_rotation_axis(u_so3)
+    commutator_axis = _compute_rotation_axis(commutator)
+
+    sim_matrix = _compute_rotation_between(commutator_axis, u_so3_axis)
+    sim_matrix_dagger = np.conj(sim_matrix).T
+
+    v = np.dot(np.dot(sim_matrix, vx), sim_matrix_dagger)
+    w = np.dot(np.dot(sim_matrix, wy), sim_matrix_dagger)
+
+    return GateSequence.from_matrix(v), GateSequence.from_matrix(w)
diff --git a/qiskit/synthesis/discrete_basis/gate_sequence.py b/qiskit/synthesis/discrete_basis/gate_sequence.py
new file mode 100644
index 000000000000..751246a24f19
--- /dev/null
+++ b/qiskit/synthesis/discrete_basis/gate_sequence.py
@@ -0,0 +1,414 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 2020.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Algebra utilities and the ``GateSequence`` class."""
+
+from __future__ import annotations
+
+from collections.abc import Sequence
+import math
+import numpy as np
+
+from qiskit.circuit import Gate, QuantumCircuit, Qubit
+from qiskit.dagcircuit import DAGCircuit
+from qiskit.extensions import UnitaryGate
+
+
+class GateSequence:
+    """A class implementing a sequence of gates.
+
+    This class stores the sequence of gates along with the unitary they implement.
+    """
+
+    def __init__(self, gates: Sequence[Gate] = ()) -> None:
+        """Create a new sequence of gates.
+
+        Args:
+            gates: The gates in the sequence. The default is [].
+        """
+        self.gates = list(gates)
+        self.matrices = [np.asarray(gate, dtype=np.complex128) for gate in gates]
+        self.labels = [gate.name for gate in gates]
+
+        # get U(2) representation of the gate sequence
+        u2_matrix = np.identity(2)
+        for matrix in self.matrices:
+            # idea: could this be optimized by a specific numpy operation?
+            u2_matrix = matrix.dot(u2_matrix)
+
+        # convert to SU(2)
+        su2_matrix, global_phase = _convert_u2_to_su2(u2_matrix)
+
+        # convert to SO(3), that's what the Solovay Kitaev algorithm uses
+        so3_matrix = _convert_su2_to_so3(su2_matrix)
+
+        # store the matrix and the global phase
+        self._eulers = None
+        self.name = " ".join(self.labels)
+        self.global_phase = global_phase
+        self.product = so3_matrix
+        self.product_su2 = su2_matrix
+
+    def remove_cancelling_pair(self, indices: Sequence[int]) -> None:
+        """Remove a pair of indices that cancel each other and *do not* change the matrices."""
+        for index in list(indices[::-1]):
+            self.gates.pop(index)
+            self.labels.pop(index)
+
+        # restore name
+        self.name = " ".join(self.labels)
+
+    def __eq__(self, other: "GateSequence") -> bool:
+        """Check if this GateSequence is the same as the other GateSequence.
+
+        Args:
+            other: The GateSequence that will be compared to ``self``.
+
+        Returns:
+            True if ``other`` is equivalent to ``self``, false otherwise.
+
+        """
+        if not len(self.gates) == len(other.gates):
+            return False
+
+        for gate1, gate2 in zip(self.gates, other.gates):
+            if gate1 != gate2:
+                return False
+
+        if self.global_phase != other.global_phase:
+            return False
+
+        return True
+
+    def to_circuit(self):
+        """Convert to a circuit.
+
+        If no gates set but the product is not the identity, returns a circuit with a
+        unitary operation to implement the matrix.
+        """
+        if len(self.gates) == 0 and not np.allclose(self.product, np.identity(3)):
+            circuit = QuantumCircuit(1, global_phase=self.global_phase)
+            su2 = _convert_so3_to_su2(self.product)
+            circuit.unitary(su2, [0])
+            return circuit
+
+        circuit = QuantumCircuit(1, global_phase=self.global_phase)
+        for gate in self.gates:
+            circuit.append(gate, [0])
+
+        return circuit
+
+    def to_dag(self):
+        """Convert to a :class:`.DAGCircuit`.
+
+        If no gates set but the product is not the identity, returns a circuit with a
+        unitary operation to implement the matrix.
+        """
+        qreg = [Qubit()]
+        dag = DAGCircuit()
+        dag.add_qubits(qreg)
+
+        if len(self.gates) == 0 and not np.allclose(self.product, np.identity(3)):
+            su2 = _convert_so3_to_su2(self.product)
+            dag.apply_operation_back(UnitaryGate(su2), qreg)
+            return dag
+
+        dag.global_phase = self.global_phase
+        for gate in self.gates:
+            dag.apply_operation_back(gate, qreg)
+
+        return dag
+
+    def append(self, gate: Gate) -> "GateSequence":
+        """Append gate to the sequence of gates.
+
+        Args:
+            gate: The gate to be appended.
+
+        Returns:
+            GateSequence with ``gate`` appended.
+        """
+        # invalidate euler angles and name
+        self._eulers = None
+
+        # TODO: this recomputes the product whenever we append something, which could be more
+        # efficient by storing the current matrix and just multiplying the input gate to it
+        # self.product = convert_su2_to_so3(self._compute_product(self.gates))
+        matrix = np.array(gate, dtype=np.complex128)
+        su2, phase = _convert_u2_to_su2(matrix)
+        so3 = _convert_su2_to_so3(su2)
+
+        self.product = so3.dot(self.product)
+        self.product_su2 = su2.dot(self.product_su2)
+        self.global_phase = self.global_phase + phase
+
+        self.gates.append(gate)
+        if len(self.labels) > 0:
+            self.name += f" {gate.name}"
+        else:
+            self.name = gate.name
+        self.labels.append(gate.name)
+
+        self.matrices.append(matrix)
+
+        return self
+
+    def adjoint(self) -> "GateSequence":
+        """Get the complex conjugate."""
+        # We're initializing an empty GateSequence and set the state manually, as we can
+        # efficiently infer the adjoint values from the current value instead of recomputing them.
+        adjoint = GateSequence()
+        adjoint.gates = [gate.inverse() for gate in reversed(self.gates)]
+        adjoint.labels = [inv.name for inv in adjoint.gates]
+        adjoint.name = " ".join(adjoint.labels)
+        adjoint.product = np.conj(self.product).T
+        adjoint.product_su2 = np.conj(self.product_su2).T
+        adjoint.global_phase = -self.global_phase
+
+        return adjoint
+
+    def copy(self) -> "GateSequence":
+        """Create copy of the sequence of gates.
+
+        Returns:
+            A new ``GateSequence`` containing copy of list of gates.
+
+        """
+        out = type(self).__new__(type(self))
+        out.labels = self.labels.copy()
+        out.gates = self.gates.copy()
+        out.matrices = self.matrices.copy()
+        out.global_phase = self.global_phase
+        out.product = self.product.copy()
+        out.product_su2 = self.product_su2.copy()
+        out.name = self.name
+        out._eulers = self._eulers
+        return out
+
+    def __len__(self) -> int:
+        """Return length of sequence of gates.
+
+        Returns:
+            Length of list containing gates.
+        """
+        return len(self.gates)
+
+    def __getitem__(self, index: int) -> Gate:
+        """Returns the gate at ``index`` from the list of gates.
+
+        Args
+            index: Index of gate in list that will be returned.
+
+        Returns:
+            The gate at ``index`` in the list of gates.
+        """
+        return self.gates[index]
+
+    def __repr__(self) -> str:
+        """Return string representation of this object.
+
+        Returns:
+            Representation of this sequence of gates.
+        """
+        out = "["
+        for gate in self.gates:
+            out += gate.name
+            out += ", "
+        out += "]"
+        out += ", product: "
+        out += str(self.product)
+        return out
+
+    def __str__(self) -> str:
+        """Return string representation of this object.
+
+        Returns:
+            Representation of this sequence of gates.
+        """
+        out = "["
+        for gate in self.gates:
+            out += gate.name
+            out += ", "
+        out += "]"
+        out += ", product: \n"
+        out += str(self.product)
+        return out
+
+    @classmethod
+    def from_matrix(cls, matrix: np.ndarray) -> "GateSequence":
+        """Initialize the gate sequence from a matrix, without a gate sequence.
+
+        Args:
+            matrix: The matrix, can be SU(2) or SO(3).
+
+        Returns:
+            A ``GateSequence`` initialized from the input matrix.
+
+        Raises:
+            ValueError: If the matrix has an invalid shape.
+        """
+        instance = cls()
+        if matrix.shape == (2, 2):
+            instance.product = _convert_su2_to_so3(matrix)
+        elif matrix.shape == (3, 3):
+            instance.product = matrix
+        else:
+            raise ValueError(f"Matrix must have shape (3, 3) or (2, 2) but has {matrix.shape}.")
+
+        instance.gates = []
+        return instance
+
+    def dot(self, other: "GateSequence") -> "GateSequence":
+        """Compute the dot-product with another gate sequence.
+
+        Args:
+            other: The other gate sequence.
+
+        Returns:
+            The dot-product as gate sequence.
+        """
+        # We're initializing an empty GateSequence and set the state manually, as we can more
+        # efficiently compute the multiplied values from the already constructed matrices.
+        composed = GateSequence()
+        composed.gates = other.gates + self.gates
+        composed.labels = other.labels + self.labels
+        composed.name = " ".join(composed.labels)
+        composed.product = np.dot(self.product, other.product)
+        composed.global_phase = self.global_phase + other.global_phase
+
+        return composed
+
+
+def _convert_u2_to_su2(u2_matrix: np.ndarray) -> tuple[np.ndarray, float]:
+    """Convert a U(2) matrix to SU(2) by adding a global phase."""
+    z = 1 / np.sqrt(np.linalg.det(u2_matrix))
+    su2_matrix = z * u2_matrix
+    phase = np.arctan2(np.imag(z), np.real(z))
+
+    return su2_matrix, phase
+
+
+def _compute_euler_angles_from_so3(matrix: np.ndarray) -> tuple[float, float, float]:
+    """Computes the Euler angles from the SO(3)-matrix u.
+
+    Uses the algorithm from Gregory Slabaugh,
+    see `here <https://www.gregslabaugh.net/publications/euler.pdf>`_.
+
+    Args:
+        matrix: The SO(3)-matrix for which the Euler angles need to be computed.
+
+    Returns:
+        Tuple (phi, theta, psi), where phi is rotation about z-axis, theta rotation about y-axis
+        and psi rotation about x-axis.
+    """
+    matrix = np.round(matrix, decimals=10)
+    if matrix[2][0] != 1 and matrix[2][1] != -1:
+        theta = -math.asin(matrix[2][0])
+        psi = math.atan2(matrix[2][1] / math.cos(theta), matrix[2][2] / math.cos(theta))
+        phi = math.atan2(matrix[1][0] / math.cos(theta), matrix[0][0] / math.cos(theta))
+        return phi, theta, psi
+    else:
+        phi = 0
+        if matrix[2][0] == 1:
+            theta = math.pi / 2
+            psi = phi + math.atan2(matrix[0][1], matrix[0][2])
+        else:
+            theta = -math.pi / 2
+            psi = -phi + math.atan2(-matrix[0][1], -matrix[0][2])
+        return phi, theta, psi
+
+
+def _compute_su2_from_euler_angles(angles: tuple[float, float, float]) -> np.ndarray:
+    """Computes SU(2)-matrix from Euler angles.
+
+    Args:
+        angles: The tuple containing the Euler angles for which the corresponding SU(2)-matrix
+            needs to be computed.
+
+    Returns:
+        The SU(2)-matrix corresponding to the Euler angles in angles.
+    """
+    phi, theta, psi = angles
+    uz_phi = np.array([[np.exp(-0.5j * phi), 0], [0, np.exp(0.5j * phi)]], dtype=complex)
+    uy_theta = np.array(
+        [[math.cos(theta / 2), math.sin(theta / 2)], [-math.sin(theta / 2), math.cos(theta / 2)]],
+        dtype=complex,
+    )
+    ux_psi = np.array(
+        [[math.cos(psi / 2), math.sin(psi / 2) * 1j], [math.sin(psi / 2) * 1j, math.cos(psi / 2)]],
+        dtype=complex,
+    )
+    return np.dot(uz_phi, np.dot(uy_theta, ux_psi))
+
+
+def _convert_su2_to_so3(matrix: np.ndarray) -> np.ndarray:
+    """Computes SO(3)-matrix from input SU(2)-matrix.
+
+    Args:
+        matrix: The SU(2)-matrix for which a corresponding SO(3)-matrix needs to be computed.
+
+    Returns:
+        The SO(3)-matrix corresponding to ``matrix``.
+
+    Raises:
+        ValueError: if ``matrix`` is not an SU(2)-matrix.
+    """
+    _check_is_su2(matrix)
+
+    matrix = matrix.astype(complex)
+    a = np.real(matrix[0][0])
+    b = np.imag(matrix[0][0])
+    c = -np.real(matrix[0][1])
+    d = -np.imag(matrix[0][1])
+    rotation = np.array(
+        [
+            [a**2 - b**2 - c**2 + d**2, 2 * a * b + 2 * c * d, -2 * a * c + 2 * b * d],
+            [-2 * a * b + 2 * c * d, a**2 - b**2 + c**2 - d**2, 2 * a * d + 2 * b * c],
+            [2 * a * c + 2 * b * d, 2 * b * c - 2 * a * d, a**2 + b**2 - c**2 - d**2],
+        ],
+        dtype=float,
+    )
+    return rotation
+
+
+def _convert_so3_to_su2(matrix: np.ndarray) -> np.ndarray:
+    """Converts an SO(3)-matrix to a corresponding SU(2)-matrix.
+
+    Args:
+        matrix: SO(3)-matrix to convert.
+
+    Returns:
+        SU(2)-matrix corresponding to SO(3)-matrix ``matrix``.
+
+    Raises:
+        ValueError: if ``matrix`` is not an SO(3)-matrix.
+    """
+    _check_is_so3(matrix)
+    return _compute_su2_from_euler_angles(_compute_euler_angles_from_so3(matrix))
+
+
+def _check_is_su2(matrix: np.ndarray) -> None:
+    """Check whether ``matrix`` is SU(2), otherwise raise an error."""
+    if matrix.shape != (2, 2):
+        raise ValueError(f"Matrix must have shape (2, 2) but has {matrix.shape}.")
+
+    if abs(np.linalg.det(matrix) - 1) > 1e-4:
+        raise ValueError(f"Determinant of matrix must be 1, but is {np.linalg.det(matrix)}.")
+
+
+def _check_is_so3(matrix: np.ndarray) -> None:
+    """Check whether ``matrix`` is SO(3), otherwise raise an error."""
+    if matrix.shape != (3, 3):
+        raise ValueError(f"Matrix must have shape (3, 3) but has {matrix.shape}.")
+
+    if abs(np.linalg.det(matrix) - 1) > 1e-4:
+        raise ValueError(f"Determinant of matrix must be 1, but is {np.linalg.det(matrix)}.")
diff --git a/qiskit/synthesis/discrete_basis/generate_basis_approximations.py b/qiskit/synthesis/discrete_basis/generate_basis_approximations.py
new file mode 100644
index 000000000000..5fe8302b1ea4
--- /dev/null
+++ b/qiskit/synthesis/discrete_basis/generate_basis_approximations.py
@@ -0,0 +1,163 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 2022.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Functions to generate the basic approximations of single qubit gates for Solovay-Kitaev."""
+
+from __future__ import annotations
+
+import warnings
+import collections
+import numpy as np
+
+import qiskit.circuit.library.standard_gates as gates
+from qiskit.circuit import Gate
+from qiskit.quantum_info.operators.predicates import matrix_equal
+from qiskit.utils import optionals
+
+from .gate_sequence import GateSequence
+
+Node = collections.namedtuple("Node", ("labels", "sequence", "children"))
+
+_1q_inverses = {
+    "i": "i",
+    "x": "x",
+    "y": "y",
+    "z": "z",
+    "h": "h",
+    "t": "tdg",
+    "tdg": "t",
+    "s": "sdg",
+    "sgd": "s",
+}
+
+_1q_gates = {
+    "i": gates.IGate(),
+    "x": gates.XGate(),
+    "y": gates.YGate(),
+    "z": gates.ZGate(),
+    "h": gates.HGate(),
+    "t": gates.TGate(),
+    "tdg": gates.TdgGate(),
+    "s": gates.SGate(),
+    "sdg": gates.SdgGate(),
+    "sx": gates.SXGate(),
+    "sxdg": gates.SXdgGate(),
+}
+
+
+def _check_candidate(candidate, existing_sequences, tol=1e-10):
+    if optionals.HAS_SKLEARN:
+        return _check_candidate_kdtree(candidate, existing_sequences, tol)
+
+    warnings.warn(
+        "The SolovayKitaev algorithm relies on scikit-learn's KDTree for a "
+        "fast search over the basis approximations. Without this, we fallback onto a "
+        "greedy search with is significantly slower. We highly suggest to install "
+        "scikit-learn to use this feature.",
+        category=RuntimeWarning,
+    )
+    return _check_candidate_greedy(candidate, existing_sequences, tol)
+
+
+def _check_candidate_greedy(candidate, existing_sequences, tol=1e-10):
+    # do a quick, string-based check if the same sequence already exists
+    if any(candidate.name == existing.name for existing in existing_sequences):
+        return False
+
+    for existing in existing_sequences:
+        if matrix_equal(existing.product_su2, candidate.product_su2, ignore_phase=True, atol=tol):
+            # is the new sequence less or more efficient?
+            return len(candidate.gates) < len(existing.gates)
+    return True
+
+
+@optionals.HAS_SKLEARN.require_in_call
+def _check_candidate_kdtree(candidate, existing_sequences, tol=1e-10):
+    """Check if there's a candidate implementing the same matrix up to ``tol``.
+
+    This uses a k-d tree search and is much faster than the greedy, list-based search.
+    """
+    from sklearn.neighbors import KDTree
+
+    # do a quick, string-based check if the same sequence already exists
+    if any(candidate.name == existing.name for existing in existing_sequences):
+        return False
+
+    points = np.array([sequence.product.flatten() for sequence in existing_sequences])
+    candidate = np.array([candidate.product.flatten()])
+
+    kdtree = KDTree(points)
+    dist, _ = kdtree.query(candidate)
+
+    return dist[0][0] > tol
+
+
+def _process_node(node: Node, basis: list[str], sequences: list[GateSequence]):
+    inverse_last = _1q_inverses[node.labels[-1]] if node.labels else None
+
+    for label in basis:
+        if label == inverse_last:
+            continue
+
+        sequence = node.sequence.copy()
+        sequence.append(_1q_gates[label])
+
+        if _check_candidate(sequence, sequences):
+            sequences.append(sequence)
+            node.children.append(Node(node.labels + (label,), sequence, []))
+
+    return node.children
+
+
+def generate_basic_approximations(
+    basis_gates: list[str | Gate], depth: int, filename: str | None = None
+) -> list[GateSequence]:
+    """Generates a list of ``GateSequence``s with the gates in ``basic_gates``.
+
+    Args:
+        basis_gates: The gates from which to create the sequences of gates.
+        depth: The maximum depth of the approximations.
+        filename: If provided, the basic approximations are stored in this file.
+
+    Returns:
+        List of ``GateSequences`` using the gates in ``basic_gates``.
+
+    Raises:
+        ValueError: If ``basis_gates`` contains an invalid gate identifier.
+    """
+    basis = []
+    for gate in basis_gates:
+        if isinstance(gate, str):
+            if gate not in _1q_gates.keys():
+                raise ValueError(f"Invalid gate identifier: {gate}")
+            basis.append(gate)
+        else:  # gate is a qiskit.circuit.Gate
+            basis.append(gate.name)
+
+    tree = Node((), GateSequence(), [])
+    cur_level = [tree]
+    sequences = [tree.sequence]
+    for _ in [None] * depth:
+        next_level = []
+        for node in cur_level:
+            next_level.extend(_process_node(node, basis, sequences))
+        cur_level = next_level
+
+    if filename is not None:
+        data = {}
+        for sequence in sequences:
+            gatestring = sequence.name
+            data[gatestring] = sequence.product
+
+        np.save(filename, data)
+
+    return sequences
diff --git a/qiskit/synthesis/discrete_basis/solovay_kitaev.py b/qiskit/synthesis/discrete_basis/solovay_kitaev.py
new file mode 100644
index 000000000000..5657452f702a
--- /dev/null
+++ b/qiskit/synthesis/discrete_basis/solovay_kitaev.py
@@ -0,0 +1,198 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 2020.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Synthesize a single qubit gate to a discrete basis set."""
+
+from __future__ import annotations
+
+import numpy as np
+
+from qiskit.circuit import QuantumCircuit, Gate
+from qiskit.dagcircuit import DAGCircuit
+
+from .gate_sequence import GateSequence
+from .commutator_decompose import commutator_decompose
+from .generate_basis_approximations import generate_basic_approximations, _1q_gates, _1q_inverses
+
+
+class SolovayKitaevDecomposition:
+    """The Solovay Kitaev discrete decomposition algorithm.
+
+    This class is called recursively by the transpiler pass, which is why it is separeted.
+    See :mod:`qiskit.synthesis.discrete_basis` for more information.
+    """
+
+    def __init__(
+        self, basic_approximations: str | dict[str, np.ndarray] | list[GateSequence] | None = None
+    ) -> None:
+        """
+        Args:
+            basic_approximations: A specification of the basic SU(2) approximations in terms
+                of discrete gates. At each iteration this algorithm, the remaining error is
+                approximated with the closest sequence of gates in this set.
+                If a ``str``, this specifies a ``.npy`` filename from which to load the
+                approximation. If a ``dict``, then this contains
+                ``{gates: effective_SO3_matrix}`` pairs,
+                e.g. ``{"h t": np.array([[0, 0.7071, -0.7071], [0, -0.7071, -0.7071], [-1, 0, 0]]}``.
+                If a list, this contains the same information as the dict, but already converted to
+                :class:`.GateSequence` objects, which contain the SO(3) matrix and gates.
+        """
+        if basic_approximations is None:
+            # generate a default basic approximation
+            basic_approximations = generate_basic_approximations(
+                basis_gates=["h", "t", "tdg"], depth=10
+            )
+
+        self.basic_approximations = self.load_basic_approximations(basic_approximations)
+
+    def load_basic_approximations(self, data: list | str | dict) -> list[GateSequence]:
+        """Load basic approximations.
+
+        Args:
+            data: If a string, specifies the path to the file from where to load the data.
+                If a dictionary, directly specifies the decompositions as ``{gates: matrix}``.
+                There ``gates`` are the names of the gates producing the SO(3) matrix ``matrix``,
+                e.g. ``{"h t": np.array([[0, 0.7071, -0.7071], [0, -0.7071, -0.7071], [-1, 0, 0]]}``.
+
+        Returns:
+            A list of basic approximations as type ``GateSequence``.
+
+        Raises:
+            ValueError: If the number of gate combinations and associated matrices does not match.
+        """
+        # is already a list of GateSequences
+        if isinstance(data, list):
+            return data
+
+        # if a file, load the dictionary
+        if isinstance(data, str):
+            data = np.load(data, allow_pickle=True)
+
+        sequences = []
+        for gatestring, matrix in data.items():
+            sequence = GateSequence()
+            sequence.gates = [_1q_gates[element] for element in gatestring.split()]
+            sequence.product = np.asarray(matrix)
+            sequences.append(sequence)
+
+        return sequences
+
+    def run(
+        self, gate_matrix: np.ndarray, recursion_degree: int, return_dag: bool = False
+    ) -> QuantumCircuit | DAGCircuit:
+        r"""Run the algorithm.
+
+        Args:
+            gate_matrix: The 2x2 matrix representing the gate. Does not need to be SU(2).
+            recursion_degree: The recursion degree, called :math:`n` in the paper.
+            return_dag: If ``True`` return a :class:`.DAGCircuit`, else a :class:`.QuantumCircuit`.
+
+        Returns:
+            A one-qubit circuit approximating the ``gate_matrix`` in the specified discrete basis.
+        """
+        # make input matrix SU(2) and get the according global phase
+        z = 1 / np.sqrt(np.linalg.det(gate_matrix))
+        gate_matrix_su2 = GateSequence.from_matrix(z * gate_matrix)
+        global_phase = np.arctan2(np.imag(z), np.real(z))
+
+        # get the decompositon as GateSequence type
+        decomposition = self._recurse(gate_matrix_su2, recursion_degree)
+
+        # simplify
+        _remove_identities(decomposition)
+        _remove_inverse_follows_gate(decomposition)
+
+        # convert to a circuit and attach the right phases
+        if return_dag:
+            out = decomposition.to_dag()
+        else:
+            out = decomposition.to_circuit()
+
+        out.global_phase = decomposition.global_phase - global_phase
+
+        return out
+
+    def _recurse(self, sequence: GateSequence, n: int) -> GateSequence:
+        """Performs ``n`` iterations of the Solovay-Kitaev algorithm on ``sequence``.
+
+        Args:
+            sequence: GateSequence to which the Solovay-Kitaev algorithm is applied.
+            n: number of iterations that the algorithm needs to run.
+
+        Returns:
+            GateSequence that approximates ``sequence``.
+
+        Raises:
+            ValueError: if ``u`` does not represent an SO(3)-matrix.
+        """
+        if sequence.product.shape != (3, 3):
+            raise ValueError("Shape of U must be (3, 3) but is", sequence.shape)
+
+        if n == 0:
+            return self.find_basic_approximation(sequence)
+
+        u_n1 = self._recurse(sequence, n - 1)
+
+        v_n, w_n = commutator_decompose(sequence.dot(u_n1.adjoint()).product)
+
+        v_n1 = self._recurse(v_n, n - 1)
+        w_n1 = self._recurse(w_n, n - 1)
+        return v_n1.dot(w_n1).dot(v_n1.adjoint()).dot(w_n1.adjoint()).dot(u_n1)
+
+    def find_basic_approximation(self, sequence: GateSequence) -> Gate:
+        """Finds gate in ``self._basic_approximations`` that best represents ``sequence``.
+
+        Args:
+            sequence: The gate to find the approximation to.
+
+        Returns:
+            Gate in basic approximations that is closest to ``sequence``.
+        """
+        # TODO explore using a k-d tree here
+
+        def key(x):
+            return np.linalg.norm(np.subtract(x.product, sequence.product))
+
+        best = min(self.basic_approximations, key=key)
+        return best
+
+
+def _remove_inverse_follows_gate(sequence):
+    index = 0
+    while index < len(sequence.gates) - 1:
+        curr_gate = sequence.gates[index]
+        next_gate = sequence.gates[index + 1]
+        if curr_gate.name in _1q_inverses.keys():
+            remove = _1q_inverses[curr_gate.name] == next_gate.name
+        else:
+            remove = curr_gate.inverse() == next_gate
+
+        if remove:
+            # remove gates at index and index + 1
+            sequence.remove_cancelling_pair([index, index + 1])
+            # take a step back to see if we have uncovered a new pair, e.g.
+            # [h, s, sdg, h] at index = 1 removes s, sdg but if we continue at index 1
+            # we miss the uncovered [h, h] pair at indices 0 and 1
+            if index > 0:
+                index -= 1
+        else:
+            # next index
+            index += 1
+
+
+def _remove_identities(sequence):
+    index = 0
+    while index < len(sequence.gates):
+        if sequence.gates[index].name == "id":
+            sequence.gates.pop(index)
+        else:
+            index += 1
diff --git a/qiskit/transpiler/passes/__init__.py b/qiskit/transpiler/passes/__init__.py
index 5148175daff9..c9c5ec33e690 100644
--- a/qiskit/transpiler/passes/__init__.py
+++ b/qiskit/transpiler/passes/__init__.py
@@ -134,7 +134,7 @@
    DAGLongestPath
 
 Synthesis
-=============
+=========
 
 .. autosummary::
    :toctree: ../stubs/
@@ -143,6 +143,8 @@
    LinearFunctionsSynthesis
    LinearFunctionsToPermutations
    HighLevelSynthesis
+   SolovayKitaev
+   SolovayKitaevSynthesis
 
 Post Layout (Post transpile qubit selection)
 ============================================
@@ -245,6 +247,8 @@
 from .synthesis import LinearFunctionsSynthesis
 from .synthesis import LinearFunctionsToPermutations
 from .synthesis import HighLevelSynthesis
+from .synthesis import SolovayKitaev
+from .synthesis import SolovayKitaevSynthesis
 
 # calibration
 from .calibration import PulseGates
diff --git a/qiskit/transpiler/passes/optimization/collect_1q_runs.py b/qiskit/transpiler/passes/optimization/collect_1q_runs.py
index cc9a28d3efcf..f4f8e1f67223 100644
--- a/qiskit/transpiler/passes/optimization/collect_1q_runs.py
+++ b/qiskit/transpiler/passes/optimization/collect_1q_runs.py
@@ -24,7 +24,7 @@ def run(self, dag):
         The blocks contain "op" nodes in topological order such that all gates
         in a block act on the same qubits and are adjacent in the circuit.
 
-        After the execution, ``property_set['block_list']`` is set to a list of
+        After the execution, ``property_set['run_list']`` is set to a list of
         tuples of "op" node.
         """
         self.property_set["run_list"] = dag.collect_1q_runs()
diff --git a/qiskit/transpiler/passes/optimization/consolidate_blocks.py b/qiskit/transpiler/passes/optimization/consolidate_blocks.py
index 99df3a4d8ee7..3b74a470e533 100644
--- a/qiskit/transpiler/passes/optimization/consolidate_blocks.py
+++ b/qiskit/transpiler/passes/optimization/consolidate_blocks.py
@@ -73,7 +73,7 @@ def run(self, dag):
 
         # compute ordered indices for the global circuit wires
         global_index_map = {wire: idx for idx, wire in enumerate(dag.qubits)}
-        blocks = self.property_set["block_list"]
+        blocks = self.property_set["block_list"] or []
         basis_gate_name = self.decomposer.gate.name
         all_block_gates = set()
         for block in blocks:
diff --git a/qiskit/transpiler/passes/synthesis/__init__.py b/qiskit/transpiler/passes/synthesis/__init__.py
index 2a12be0f70dc..11540d89387a 100644
--- a/qiskit/transpiler/passes/synthesis/__init__.py
+++ b/qiskit/transpiler/passes/synthesis/__init__.py
@@ -16,3 +16,4 @@
 from .plugin import unitary_synthesis_plugin_names
 from .linear_functions_synthesis import LinearFunctionsSynthesis, LinearFunctionsToPermutations
 from .high_level_synthesis import HighLevelSynthesis, HLSConfig
+from .solovay_kitaev_synthesis import SolovayKitaev, SolovayKitaevSynthesis
diff --git a/qiskit/transpiler/passes/synthesis/solovay_kitaev_synthesis.py b/qiskit/transpiler/passes/synthesis/solovay_kitaev_synthesis.py
new file mode 100644
index 000000000000..ab5a530dcb6a
--- /dev/null
+++ b/qiskit/transpiler/passes/synthesis/solovay_kitaev_synthesis.py
@@ -0,0 +1,196 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+"""
+A Solovay-Kitaev synthesis plugin to Qiskit's transpiler.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from qiskit.converters import circuit_to_dag
+from qiskit.dagcircuit import DAGCircuit
+from qiskit.synthesis.discrete_basis.solovay_kitaev import SolovayKitaevDecomposition
+from qiskit.synthesis.discrete_basis.generate_basis_approximations import (
+    generate_basic_approximations,
+)
+from qiskit.transpiler.basepasses import TransformationPass
+from qiskit.transpiler.exceptions import TranspilerError
+
+from .plugin import UnitarySynthesisPlugin
+
+
+class SolovayKitaev(TransformationPass):
+    r"""Approximately decompose 1q gates to a discrete basis using the Solovay-Kitaev algorithm.
+
+    See :mod:`qiskit.synthesis.discrete_basis` for more information.
+    """
+
+    def __init__(
+        self,
+        recursion_degree: int = 3,
+        basic_approximations: str | dict[str, np.ndarray] | None = None,
+    ) -> None:
+        """
+        Args:
+            recursion_degree: The recursion depth for the Solovay-Kitaev algorithm.
+                A larger recursion depth increases the accuracy and length of the
+                decomposition.
+            basic_approximations: The basic approximations for the finding the best discrete
+                decomposition at the root of the recursion. If a string, it specifies the ``.npy``
+                file to load the approximations from. If a dictionary, it contains
+                ``{label: SO(3)-matrix}`` pairs. If None, a default based on the H, T and Tdg gates
+                up to combinations of depth 10 is generated.
+        """
+        super().__init__()
+        self.recursion_degree = recursion_degree
+        self._sk = SolovayKitaevDecomposition(basic_approximations)
+
+    def run(self, dag: DAGCircuit) -> DAGCircuit:
+        """Run the ``SolovayKitaev`` pass on `dag`.
+
+        Args:
+            dag: The input dag.
+
+        Returns:
+            Output dag with 1q gates synthesized in the discrete target basis.
+
+        Raises:
+            TranspilerError: if a gates does not have to_matrix
+        """
+        for node in dag.op_nodes():
+            if not node.op.num_qubits == 1:
+                continue  # ignore all non-single qubit gates
+
+            if not hasattr(node.op, "to_matrix"):
+                raise TranspilerError(
+                    "SolovayKitaev does not support gate without "
+                    f"to_matrix method: {node.op.name}"
+                )
+
+            matrix = node.op.to_matrix()
+
+            # call solovay kitaev
+            approximation = self._sk.run(matrix, self.recursion_degree, return_dag=True)
+
+            # convert to a dag and replace the gate by the approximation
+            dag.substitute_node_with_dag(node, approximation)
+
+        return dag
+
+
+class SolovayKitaevSynthesis(UnitarySynthesisPlugin):
+    """A Solovay-Kitaev Qiskit unitary synthesis plugin.
+
+    This plugin is invoked by :func:`~.compiler.transpile` when the ``unitary_synthesis_method``
+    parameter is set to ``"sk"``.
+
+    This plugin supports customization and additional parameters can be passed to the plugin
+    by passing a dictionary as the ``unitary_synthesis_plugin_config`` parameter of
+    the :func:`~qiskit.compiler.transpile` function.
+
+    Supported parameters in the dictionary:
+
+    basis_approximations (str | dict):
+        The basic approximations for the finding the best discrete decomposition at the root of the
+        recursion. If a string, it specifies the ``.npy`` file to load the approximations from.
+        If a dictionary, it contains ``{label: SO(3)-matrix}`` pairs. If None, a default based on
+        the specified ``basis_gates`` and ``depth`` is generated.
+
+    basis_gates (list):
+        A list of strings specifying the discrete basis gates to decompose to. If None,
+        defaults to ``["h", "t", "tdg"]``.
+
+    depth (int):
+        The gate-depth of the the basic approximations. All possible, unique combinations of the
+        basis gates up to length ``depth`` are considered. If None, defaults to 10.
+
+    recursion_degree (int):
+        The number of times the decomposition is recursively improved. If None, defaults to 3.
+    """
+
+    # we cache an instance of the Solovay-Kitaev class to generate the
+    # computationally expensive basis approximation of single qubit gates only once
+    _sk = None
+
+    @property
+    def max_qubits(self):
+        """Maximum number of supported qubits is ``1``."""
+        return 1
+
+    @property
+    def min_qubits(self):
+        """Minimum number of supported qubits is ``1``."""
+        return 1
+
+    @property
+    def supports_natural_direction(self):
+        """The plugin does not support natural direction, it does not assume
+        bidirectional two qubit gates."""
+        return True
+
+    @property
+    def supports_pulse_optimize(self):
+        """The plugin does not support optimization of pulses."""
+        return False
+
+    @property
+    def supports_gate_lengths(self):
+        """The plugin does not support gate lengths."""
+        return False
+
+    @property
+    def supports_gate_errors(self):
+        """The plugin does not support gate errors."""
+        return False
+
+    @property
+    def supported_bases(self):
+        """The plugin does not support bases for synthesis."""
+        return None
+
+    @property
+    def supports_basis_gates(self):
+        """The plugin does not support basis gates. By default it synthesis to the
+        ``["h", "t", "tdg"]`` gate basis."""
+        return True
+
+    @property
+    def supports_coupling_map(self):
+        """The plugin does not support coupling maps."""
+        return False
+
+    def run(self, unitary, **options):
+
+        # Runtime imports to avoid the overhead of these imports for
+        # plugin discovery and only use them if the plugin is run/used
+        config = options.get("config") or {}
+
+        recursion_degree = config.get("recursion_degree", 3)
+
+        # if we didn't yet construct the Solovay-Kitaev instance, which contains
+        # the basic approximations, do it now
+        if SolovayKitaevSynthesis._sk is None:
+            basic_approximations = config.get("basic_approximations", None)
+            basis_gates = options.get("basis_gates", ["h", "t", "tdg"])
+
+            # if the basic approximations are not generated and not given,
+            # try to generate them if the basis set is specified
+            if basic_approximations is None:
+                depth = config.get("depth", 10)
+                basic_approximations = generate_basic_approximations(basis_gates, depth)
+
+            SolovayKitaevSynthesis._sk = SolovayKitaevDecomposition(basic_approximations)
+
+        approximate_circuit = SolovayKitaevSynthesis._sk.run(unitary, recursion_degree)
+        dag_circuit = circuit_to_dag(approximate_circuit)
+        return dag_circuit
diff --git a/releasenotes/notes/solovay-kitaev-transpiler-pass-bc256c2f3aac28c6.yaml b/releasenotes/notes/solovay-kitaev-transpiler-pass-bc256c2f3aac28c6.yaml
new file mode 100644
index 000000000000..a9bacff16d54
--- /dev/null
+++ b/releasenotes/notes/solovay-kitaev-transpiler-pass-bc256c2f3aac28c6.yaml
@@ -0,0 +1,46 @@
+---
+features:
+  - |
+    Added the :class:`.SolovayKitaev` transpiler pass to run the Solovay-Kitaev algorithm for
+    approximating single-qubit unitaries using a discrete gate set. In combination with the basis
+    translator, this allows to convert any unitary circuit to a universal discrete gate set,
+    which could be implemented fault-tolerantly.
+
+    This pass can e.g. be used after compiling to U and CX gates:
+
+    .. jupyter-execute::
+
+      from qiskit import transpile
+      from qiskit.circuit.library import QFT
+      from qiskit.transpiler.passes.synthesis import SolovayKitaev
+
+      qft = QFT(3)
+
+      # optimize to general 1-qubit unitaries and CX
+      transpiled = transpile(qft, basis_gates=["u", "cx"], optimization_level=1)
+
+      skd = SolovayKitaev()  # uses T Tdg and H as default basis
+      discretized = skd(transpiled)
+
+      print(discretized.count_ops())
+
+    The decomposition can also be used with the unitary synthesis plugin, as
+    the "sk" method on the :class:`~.UnitarySynthesis` transpiler pass:
+
+    .. jupyter-execute::
+
+        from qiskit import QuantumCircuit
+        from qiskit.quantum_info import Operator
+        from qiskit.transpiler.passes import UnitarySynthesis
+
+        circuit = QuantumCircuit(1)
+        circuit.rx(0.8, 0)
+        unitary = Operator(circuit).data
+
+        unitary_circ = QuantumCircuit(1)
+        unitary_circ.unitary(unitary, [0])
+
+        synth = UnitarySynthesis(basis_gates=["h", "s"], method="sk")
+        out = synth(unitary_circ)
+
+        out.draw('mpl')
diff --git a/setup.py b/setup.py
index ebeb6653a3af..4154a56ebb0c 100755
--- a/setup.py
+++ b/setup.py
@@ -101,6 +101,7 @@
         "qiskit.unitary_synthesis": [
             "default = qiskit.transpiler.passes.synthesis.unitary_synthesis:DefaultUnitarySynthesis",
             "aqc = qiskit.transpiler.synthesis.aqc.aqc_plugin:AQCSynthesisPlugin",
+            "sk = qiskit.transpiler.passes.synthesis.solovay_kitaev_synthesis:SolovayKitaevSynthesis",
         ],
         "qiskit.synthesis": [
             "clifford.default = qiskit.transpiler.passes.synthesis.high_level_synthesis:DefaultSynthesisClifford",
diff --git a/test/python/transpiler/test_solovay_kitaev.py b/test/python/transpiler/test_solovay_kitaev.py
new file mode 100644
index 000000000000..95dfd06763d0
--- /dev/null
+++ b/test/python/transpiler/test_solovay_kitaev.py
@@ -0,0 +1,370 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 2020.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Test the Solovay Kitaev transpilation pass."""
+
+import unittest
+import math
+import numpy as np
+import scipy
+
+from ddt import ddt, data
+
+from qiskit.test import QiskitTestCase
+
+from qiskit import transpile
+from qiskit.circuit import QuantumCircuit
+from qiskit.circuit.library import TGate, TdgGate, HGate, SGate, SdgGate, IGate, QFT
+from qiskit.converters import circuit_to_dag, dag_to_circuit
+from qiskit.quantum_info import Operator
+from qiskit.synthesis.discrete_basis.generate_basis_approximations import (
+    generate_basic_approximations,
+)
+from qiskit.synthesis.discrete_basis.commutator_decompose import commutator_decompose
+from qiskit.synthesis.discrete_basis.gate_sequence import GateSequence
+from qiskit.transpiler import PassManager
+from qiskit.transpiler.exceptions import TranspilerError
+from qiskit.transpiler.passes import UnitarySynthesis, Collect1qRuns, ConsolidateBlocks
+from qiskit.transpiler.passes.synthesis import SolovayKitaev, SolovayKitaevSynthesis
+
+
+def _trace_distance(circuit1, circuit2):
+    """Return the trace distance of the two input circuits."""
+    op1, op2 = Operator(circuit1), Operator(circuit2)
+    return 0.5 * np.trace(scipy.linalg.sqrtm(np.conj(op1 - op2).T.dot(op1 - op2))).real
+
+
+def _generate_x_rotation(angle: float) -> np.ndarray:
+    return np.array(
+        [[1, 0, 0], [0, math.cos(angle), -math.sin(angle)], [0, math.sin(angle), math.cos(angle)]]
+    )
+
+
+def _generate_y_rotation(angle: float) -> np.ndarray:
+    return np.array(
+        [[math.cos(angle), 0, math.sin(angle)], [0, 1, 0], [-math.sin(angle), 0, math.cos(angle)]]
+    )
+
+
+def _generate_z_rotation(angle: float) -> np.ndarray:
+    return np.array(
+        [[math.cos(angle), -math.sin(angle), 0], [math.sin(angle), math.cos(angle), 0], [0, 0, 1]]
+    )
+
+
+def is_so3_matrix(array: np.ndarray) -> bool:
+    """Check if the input array is a SO(3) matrix."""
+    if array.shape != (3, 3):
+        return False
+
+    if abs(np.linalg.det(array) - 1.0) > 1e-10:
+        return False
+
+    if False in np.isreal(array):
+        return False
+
+    return True
+
+
+@ddt
+class TestSolovayKitaev(QiskitTestCase):
+    """Test the Solovay Kitaev algorithm and transformation pass."""
+
+    def setUp(self):
+        super().setUp()
+
+        self.basic_approx = generate_basic_approximations([HGate(), TGate(), TdgGate()], 3)
+
+    def test_unitary_synthesis(self):
+        """Test the unitary synthesis transpiler pass with Solovay-Kitaev."""
+        circuit = QuantumCircuit(2)
+        circuit.rx(0.8, 0)
+        circuit.cx(0, 1)
+        circuit.x(1)
+
+        _1q = Collect1qRuns()
+        _cons = ConsolidateBlocks()
+        _synth = UnitarySynthesis(["h", "s"], method="sk")
+        passes = PassManager([_1q, _cons, _synth])
+        compiled = passes.run(circuit)
+
+        diff = np.linalg.norm(Operator(compiled) - Operator(circuit))
+        self.assertLess(diff, 1)
+        self.assertEqual(set(compiled.count_ops().keys()), {"h", "s", "cx"})
+
+    def test_plugin(self):
+        """Test calling the plugin directly."""
+        circuit = QuantumCircuit(1)
+        circuit.rx(0.8, 0)
+
+        unitary = Operator(circuit).data
+
+        plugin = SolovayKitaevSynthesis()
+        out = plugin.run(unitary, basis_gates=["h", "s"])
+
+        reference = QuantumCircuit(1, global_phase=3 * np.pi / 4)
+        reference.h(0)
+        reference.s(0)
+        reference.h(0)
+
+        self.assertEqual(dag_to_circuit(out), reference)
+
+    def test_generating_default_approximation(self):
+        """Test the approximation set is generated by default."""
+        skd = SolovayKitaev()
+        circuit = QuantumCircuit(1)
+        dummy = skd(circuit)
+
+        self.assertIsNotNone(skd._sk.basic_approximations)
+
+    def test_i_returns_empty_circuit(self):
+        """Test that ``SolovayKitaev`` returns an empty circuit when
+        it approximates the I-gate."""
+        circuit = QuantumCircuit(1)
+        circuit.i(0)
+
+        skd = SolovayKitaev(3, self.basic_approx)
+
+        decomposed_circuit = skd(circuit)
+        self.assertEqual(QuantumCircuit(1), decomposed_circuit)
+
+    def test_exact_decomposition_acts_trivially(self):
+        """Test that the a circuit that can be represented exactly is represented exactly."""
+        circuit = QuantumCircuit(1)
+        circuit.t(0)
+        circuit.h(0)
+        circuit.tdg(0)
+
+        synth = SolovayKitaev(3, self.basic_approx)
+
+        dag = circuit_to_dag(circuit)
+        decomposed_dag = synth.run(dag)
+        decomposed_circuit = dag_to_circuit(decomposed_dag)
+        self.assertEqual(circuit, decomposed_circuit)
+
+    def test_fails_with_no_to_matrix(self):
+        """Test failer if gate does not have to_matrix."""
+        circuit = QuantumCircuit(1)
+        circuit.initialize("0")
+
+        synth = SolovayKitaev(3, self.basic_approx)
+
+        dag = circuit_to_dag(circuit)
+
+        with self.assertRaises(TranspilerError) as cm:
+            _ = synth.run(dag)
+
+        self.assertEqual(
+            "SolovayKitaev does not support gate without to_matrix method: initialize",
+            cm.exception.message,
+        )
+
+    def test_str_basis_gates(self):
+        """Test specifying the basis gates by string works."""
+        circuit = QuantumCircuit(1)
+        circuit.rx(0.8, 0)
+
+        basic_approx = generate_basic_approximations(["h", "t", "s"], 3)
+        synth = SolovayKitaev(2, basic_approx)
+
+        dag = circuit_to_dag(circuit)
+        discretized = dag_to_circuit(synth.run(dag))
+
+        reference = QuantumCircuit(1, global_phase=7 * np.pi / 8)
+        reference.h(0)
+        reference.t(0)
+        reference.h(0)
+
+        self.assertEqual(discretized, reference)
+
+    def test_approximation_on_qft(self):
+        """Test the Solovay-Kitaev decomposition on the QFT circuit."""
+        qft = QFT(3)
+        transpiled = transpile(qft, basis_gates=["u", "cx"], optimization_level=1)
+
+        skd = SolovayKitaev(1)
+
+        with self.subTest("1 recursion"):
+            discretized = skd(transpiled)
+            self.assertLess(_trace_distance(transpiled, discretized), 15)
+
+        skd.recursion_degree = 2
+        with self.subTest("2 recursions"):
+            discretized = skd(transpiled)
+            self.assertLess(_trace_distance(transpiled, discretized), 7)
+
+
+@ddt
+class TestGateSequence(QiskitTestCase):
+    """Test the ``GateSequence`` class."""
+
+    def test_append(self):
+        """Test append."""
+        seq = GateSequence([IGate()])
+        seq.append(HGate())
+
+        ref = GateSequence([IGate(), HGate()])
+        self.assertEqual(seq, ref)
+
+    def test_eq(self):
+        """Test equality."""
+        base = GateSequence([HGate(), HGate()])
+        seq1 = GateSequence([HGate(), HGate()])
+        seq2 = GateSequence([IGate()])
+        seq3 = GateSequence([HGate(), HGate()])
+        seq3.global_phase = 0.12
+        seq4 = GateSequence([IGate(), HGate()])
+
+        with self.subTest("equal"):
+            self.assertEqual(base, seq1)
+
+        with self.subTest("same product, but different repr (-> false)"):
+            self.assertNotEqual(base, seq2)
+
+        with self.subTest("differing global phase (-> false)"):
+            self.assertNotEqual(base, seq3)
+
+        with self.subTest("same num gates, but different gates (-> false)"):
+            self.assertNotEqual(base, seq4)
+
+    def test_to_circuit(self):
+        """Test converting a gate sequence to a circuit."""
+        seq = GateSequence([HGate(), HGate(), TGate(), SGate(), SdgGate()])
+        ref = QuantumCircuit(1)
+        ref.h(0)
+        ref.h(0)
+        ref.t(0)
+        ref.s(0)
+        ref.sdg(0)
+        # a GateSequence is SU(2), so add the right phase
+        z = 1 / np.sqrt(np.linalg.det(Operator(ref)))
+        ref.global_phase = np.arctan2(np.imag(z), np.real(z))
+
+        self.assertEqual(seq.to_circuit(), ref)
+
+    def test_adjoint(self):
+        """Test adjoint."""
+        seq = GateSequence([TGate(), SGate(), HGate(), IGate()])
+        inv = GateSequence([IGate(), HGate(), SdgGate(), TdgGate()])
+
+        self.assertEqual(seq.adjoint(), inv)
+
+    def test_copy(self):
+        """Test copy."""
+        seq = GateSequence([IGate()])
+        copied = seq.copy()
+        seq.gates.append(HGate())
+
+        self.assertEqual(len(seq.gates), 2)
+        self.assertEqual(len(copied.gates), 1)
+
+    @data(0, 1, 10)
+    def test_len(self, n):
+        """Test __len__."""
+        seq = GateSequence([IGate()] * n)
+        self.assertEqual(len(seq), n)
+
+    def test_getitem(self):
+        """Test __getitem__."""
+        seq = GateSequence([IGate(), HGate(), IGate()])
+        self.assertEqual(seq[0], IGate())
+        self.assertEqual(seq[1], HGate())
+        self.assertEqual(seq[2], IGate())
+
+        self.assertEqual(seq[-2], HGate())
+
+    def test_from_su2_matrix(self):
+        """Test from_matrix with an SU2 matrix."""
+        matrix = np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2)
+        matrix /= np.sqrt(np.linalg.det(matrix))
+        seq = GateSequence.from_matrix(matrix)
+
+        ref = GateSequence([HGate()])
+
+        self.assertEqual(seq.gates, list())
+        self.assertTrue(np.allclose(seq.product, ref.product))
+        self.assertEqual(seq.global_phase, 0)
+
+    def test_from_so3_matrix(self):
+        """Test from_matrix with an SO3 matrix."""
+        matrix = np.array([[0, 0, -1], [0, -1, 0], [-1, 0, 0]])
+        seq = GateSequence.from_matrix(matrix)
+
+        ref = GateSequence([HGate()])
+
+        self.assertEqual(seq.gates, list())
+        self.assertTrue(np.allclose(seq.product, ref.product))
+        self.assertEqual(seq.global_phase, 0)
+
+    def test_from_invalid_matrix(self):
+        """Test from_matrix with invalid matrices."""
+        with self.subTest("2x2 but not SU2"):
+            matrix = np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2)
+            with self.assertRaises(ValueError):
+                _ = GateSequence.from_matrix(matrix)
+
+        with self.subTest("not 2x2 or 3x3"):
+            with self.assertRaises(ValueError):
+                _ = GateSequence.from_matrix(np.array([[1]]))
+
+    def test_dot(self):
+        """Test dot."""
+        seq1 = GateSequence([HGate()])
+        seq2 = GateSequence([TGate(), SGate()])
+        composed = seq1.dot(seq2)
+
+        ref = GateSequence([TGate(), SGate(), HGate()])
+
+        # check the product matches
+        self.assertTrue(np.allclose(ref.product, composed.product))
+
+        # check the circuit & phases are equivalent
+        self.assertTrue(Operator(ref.to_circuit()).equiv(composed.to_circuit()))
+
+
+@ddt
+class TestSolovayKitaevUtils(QiskitTestCase):
+    """Test the public functions in the Solovay Kitaev utils."""
+
+    @data(
+        _generate_x_rotation(0.1),
+        _generate_y_rotation(0.2),
+        _generate_z_rotation(0.3),
+        np.dot(_generate_z_rotation(0.5), _generate_y_rotation(0.4)),
+        np.dot(_generate_y_rotation(0.5), _generate_x_rotation(0.4)),
+    )
+    def test_commutator_decompose_return_type(self, u_so3: np.ndarray):
+        """Test that ``commutator_decompose`` returns two SO(3) gate sequences."""
+        v, w = commutator_decompose(u_so3)
+        self.assertTrue(is_so3_matrix(v.product))
+        self.assertTrue(is_so3_matrix(w.product))
+        self.assertIsInstance(v, GateSequence)
+        self.assertIsInstance(w, GateSequence)
+
+    @data(
+        _generate_x_rotation(0.1),
+        _generate_y_rotation(0.2),
+        _generate_z_rotation(0.3),
+        np.dot(_generate_z_rotation(0.5), _generate_y_rotation(0.4)),
+        np.dot(_generate_y_rotation(0.5), _generate_x_rotation(0.4)),
+    )
+    def test_commutator_decompose_decomposes_correctly(self, u_so3):
+        """Test that ``commutator_decompose`` exactly decomposes the input."""
+        v, w = commutator_decompose(u_so3)
+        v_so3 = v.product
+        w_so3 = w.product
+        actual_commutator = np.dot(v_so3, np.dot(w_so3, np.dot(np.conj(v_so3).T, np.conj(w_so3).T)))
+        self.assertTrue(np.allclose(actual_commutator, u_so3))
+
+
+if __name__ == "__main__":
+    unittest.main()