3_towards_fault_tolerance.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "165778ef",
   "metadata": {},
   "source": [
    "Towards Fault Tolerant Compiling\n",
    "-----------\n",
    "Recent experiments demonstrating that error-corrected logical qubits can already achieve higher fidelities than physical qubits have us excited about Fault-Tolerant (FT) quantum compilation [1, 2]. In this tutorial, we will highlight some crucial differences between the synthesis of NISQ and FT circuits and discuss how BQSKit can be used to target this new domain.\n",
    "\n",
    "This tutorial requires the `bqskit-ft` extension package and can be install via pip:\n",
    "```bash\n",
    "pip install bqskit-ft\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c306def",
   "metadata": {},
   "source": [
    "Contents\n",
    "-----------\n",
    "* [Fault-Tolerant Synthesis in BQSKit](#fault-tolerant-synthesis-in-bqskit)\n",
    "* [Fault-Tolerant Gate Sets](#fault-tolerant-gate-sets)\n",
    "* [Under the Hood](#under-the-hood)\n",
    "    * [Single Qubit Gates](#single-qubit-gates)\n",
    "    * [Search-based Synthesis](#search-based-synthesis)\n",
    "    * [Synthesis by Diagonalization](#synthesis-by-diagonalization)\n",
    "    * [Synthesis with Constrained Optimization](#synthesis-with-constrained-optimization)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d43c22ca",
   "metadata": {},
   "source": [
    "Fault Tolerant Synthesis in BQSKit\n",
    "----------\n",
    "Let's start with an example of a circuit that we'd like to transpile to the fault-tolerant Clifford+$T$ gate set. Converting to Clifford+$T$ requires that a `FaultTolerantModel` be specified as the `MachineModel`. As we're targeting the Clifford+$T$ gate set, we will use a `CliffordTModel`. Then, we can just use `bqskit.compile`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c482b132",
   "metadata": {},
   "outputs": [],
   "source": [
    "from bqskit.ir import Circuit\n",
    "from bqskit.ir.gates import CNOTGate\n",
    "from bqskit.ir.gates import U3Gate\n",
    "\n",
    "from numpy import pi\n",
    "\n",
    "from utils.printing import print_circuit\n",
    "\n",
    "# A toy circuit to transpile\n",
    "example_circuit = Circuit(2)\n",
    "example_circuit.append_gate(U3Gate(), [0], [pi/2, 0, pi])\n",
    "example_circuit.append_gate(CNOTGate(), (0, 1))\n",
    "example_circuit.append_gate(U3Gate(), [1], [0, pi/4, 0])\n",
    "\n",
    "print('The target transpilation circuit contains:')\n",
    "print_circuit(example_circuit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65d6d056",
   "metadata": {},
   "outputs": [],
   "source": [
    "from bqskit.compiler.compile import compile\n",
    "\n",
    "from bqskit.ft.cliffordt.cliffordtmodel import CliffordTModel\n",
    "\n",
    "# Set up the machine model and compile the circuit\n",
    "machine = CliffordTModel(num_qudits=example_circuit.num_qudits)\n",
    "result = compile(example_circuit, model=machine)\n",
    "\n",
    "# Print out the resulting circuit\n",
    "print_circuit(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "635652ff",
   "metadata": {},
   "source": [
    "Fault-Tolerant Gate Sets\n",
    "----------\n",
    "Compilers are most concerned about transforming circuits when targeting near-term devices so that successful execution is maximized while conforming to the underlying machine's physical constraints. Typically, this means that circuits are optimized to use as few multi-qudit gates as possible, while requiring that circuits be transpiled to a particular gate set and mapped to a restrictive topology. BQSKit provides techniques for addressing each of these aspects of the compilation pipeline.\n",
    "\n",
    "Fault-tolerant (FT) computation poses different challenges. Unlike in the NISQ case, FT compilation assumes that robust error-corrected qubits store information. This feat is accomplished by adding a new layer to the quantum computing stack: a quantum error-correcting code [3, 4, 5]. This new protection does not come for free. Implementing arbitrary unitary operations on logical qubits becomes much more difficult for codes of interest. This is because these codes admit restrictive discrete gate sets, whereas physical qubits can implement continuous rotations. To make matters worse, only a subset of the logical gates needed for universal computation can be easily implemented. These gates are in the Clifford group (generated by $\\{H, S, CNOT\\}$), which contains some familiar faces such as the Pauli matrices ($\\{X, Y, Z\\}$).\n",
    "\n",
    "To make our logical gate set universal, we must add a single non-Clifford gate to the mix [6]. The most popular non-Clifford gate considered is the $T$ gate, which we define as \n",
    "$$\n",
    "T = \\begin{bmatrix}\n",
    "1 & 0 \\\\\n",
    "0 & e^{i\\frac{\\pi}{4}}\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "Even though we cannot implement the $T$ gate directly like we can Clifford gates, teleportation lets us inject special resource states into our circuit that can be used to implement non-Clifford gates. We can prepare these special \"magic\" states using a process called magic state distillation [7]. Because so much must happen behind the scenes to ensure that non-Clifford gates can be executed, they are typically $100$-$1000\\times$ more expensive to implement than their Clifford counterparts. This marks a considerable departure from the optimization objective in NISQ-centric synthesis algorithms. Whereas before, the primary concern was minimizing multi-qubit gates (such as $CNOT$s), now the goal is minimizing non-Clifford gates (such as $T$ gates).\n",
    "\n",
    "Now we can consider the single-qubit Clifford+$T$ gate set, which is simply $\\{H, T\\}$ (note that $T^2=S$ and $T^4=Z$). Equipped with this gate set, the Solovay-Kitaev theorem tells us that we can implement an approximation of any single-qubit unitary up to a precision of $\\epsilon$ in $O(log^{3.97}(\\frac{1}{\\epsilon}))$ gates. Along with the $CNOT$ gate, this allows us to approximate any $n$-qubit operation in a fault-tolerant way!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da44527d",
   "metadata": {},
   "source": [
    "Under the Hood\n",
    "----------\n",
    "Above we saw that compiling to fault-tolerant gate sets is as simple as setting the `MachineModel` to a `FaultTolerantModel`. But BQSKit can perform more complex compilation tasks while targeting fault-tolerant gate sets. In the remaining sections, we will highlight several advanced techniques for fault-tolerant synthesis."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1737c67",
   "metadata": {},
   "source": [
    "### Single Qubit Gates\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "603e47e5",
   "metadata": {},
   "source": [
    "In practice, the algorithm given by the Solovay-Kitaev theorem results in long sequences of gates when approximating single-qubit unitaries. The state of the art for ancilla-free approximate synthesis comes from the Ross-Selinger \"gridsynth\" algorithm [8]. This algorithm produces optimal length sequences of Clifford+$T$ gates approximating $R_Z(\\theta)$ rotations. Any single-qubit unitary can be written as:\n",
    "$$\n",
    "    U = R_Z(\\theta) \\sqrt{X} R_Z(\\phi) \\sqrt{X} R_Z(\\lambda)\n",
    "$$\n",
    "Therefore, gridsynth lets us implement arbitrary single-qubit rotations at most only $3\\times$ away from optimal. Unless $R_Z(\\theta)$ corresponds to gates already in our Clifford+$T$ gate set ($Z$, $S$, or $T$), each $R_Z(\\theta)$ is decomposed into a sequence of around $25\\log\\frac{1}{\\epsilon}$ gates (of which about $\\frac{1}{3}$ are $T$ gates). We will make use of the fact that we can optimally compile $R_Z(\\theta)$ rotations in our approach to implementing multi-qubit unitaries in FT gate sets.\n",
    "\n",
    "For the sake of this tutorial, we will target a gate set of Clifford+$T$+$R_Z(\\theta)$. Post-processing can optionally be used to convert continuous $R_Z$ gates to discrete gates."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee3fa2f6",
   "metadata": {},
   "source": [
    "Synthesis algorithms targeting the Clifford+$T$ gate set suffer a new challenge: how do you implement complex unitaries in discrete gate sets? We will present two methods of addressing this issue: 1) search-based techniques and 2) constrained optimization."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6537527",
   "metadata": {},
   "source": [
    "### Search-Based Synthesis\n",
    "Search-based synthesis algorithms search through the discrete space of candidate circuits. Although the search space grows as $O(|\\text{gateset}|^{\\text{depth}})$, this approach can still be rather effective in certain cases.\n",
    "\n",
    "Say we want to synthesize the following unitary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5f94755c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "easy_target = np.array([\n",
    "    [0.70710678, 0,           0.70710678,  0          ],\n",
    "    [0,          0.5+0.5j,    0,           0.5+0.5j   ],\n",
    "    [0.5+0.5j,   0,          -0.5-0.5j,    0          ],\n",
    "    [0,          0.70710678,  0,          -0.70710678 ],\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2f1fd23",
   "metadata": {},
   "source": [
    "We can perform A* search over discrete gate sets in BQSKit by using a `DiscreteLayerGenerator`. Furthermore, we want to make sure synthesis gives us a very good implementation. We're going to set the `success_threshold=1e-8` to ensure this. The following block should take at most a few minutes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13d89216",
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils.discrete import DiscreteLayerGenerator\n",
    "\n",
    "from bqskit.passes.synthesis import QSearchSynthesisPass\n",
    "from bqskit.compiler import Compiler\n",
    "from bqskit.ir import Circuit\n",
    "\n",
    "success_threshold = 1e-8\n",
    "layer_generator = DiscreteLayerGenerator()\n",
    "synthesis_pass = QSearchSynthesisPass(\n",
    "    success_threshold=success_threshold,\n",
    "    layer_generator=layer_generator,\n",
    ")\n",
    "\n",
    "circuit = Circuit.from_unitary(easy_target)\n",
    "\n",
    "with Compiler() as compiler:\n",
    "    circuit = compiler.compile(circuit, [synthesis_pass])\n",
    "\n",
    "print_circuit(circuit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "911af034",
   "metadata": {},
   "source": [
    "In this case, we were able to find a solution! Now let's say we want to try a slightly different unitary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3a5b2a97",
   "metadata": {},
   "outputs": [],
   "source": [
    "hard_target = np.array([\n",
    "   [0.70370187-0.06930858j, 0, 0.70370187-0.06930858j, 0],\n",
    "   [0, 0.70370187+0.06930858j, 0, 0.70370187+0.06930858j],\n",
    "   [0.70370187+0.06930858j, 0, -0.70370187-0.06930858j, 0],\n",
    "   [0, 0.70370187-0.06930858j, 0, -0.70370187+0.06930858j],\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89843bbb",
   "metadata": {},
   "source": [
    "What happens if you try to use the `DiscreteLayerGenerator` to synthesize this unitary?\n",
    "\n",
    "(Note: you may have to restart the kernel after this part)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "51ac6682",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the DiscreteLayerGenerator to try synthesizing hard_target"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a6355a7",
   "metadata": {},
   "source": [
    "### Synthesis by Diagonalization\n",
    "Pure search-based synthesis has its limits. If the unitary is not easily implementable in the discrete gate set provided to the `DiscreteLayerGenerator`, the synthesis algorithm may spend lots of time checking different candidates before eventually timing out. Many techniques can push the boundaries of what's possible with pure discrete search, but the sheer size of the solution space will always limit these techniques. Finding solutions to a precision of `success_threshold=1e-3` might be possible in certain cases, but beyond that, for high precision, a new technique is needed.\n",
    "\n",
    "To approach this problem, we're going to make an assumption: finding a solution circuit requires $R_Z(\\theta)$ gates. Making this assumption means that our synthesis algorithm is unlikely to find optimal solutions, but a good, suboptimal solution with $R_Z(\\theta)$ gates is better than no solution!\n",
    "\n",
    "How can we modify the search-based synthesis framework to handle unitaries whose implementation requires $R_Z(\\theta)$ gates? Using `QSearch` then converting every `U3Gate` to a sequence of `RZGate`s and `SqrtXGate`s is one option, which will be the topic of the next section, but there's another way.\n",
    "\n",
    "We're going to assume the circuit implementing our `hard_target` unitary is in the form:\n",
    "$$\n",
    "    C(\\bar\\theta) = L \\times D_{\\bar\\theta} \\times R\n",
    "$$\n",
    "where $D_{\\bar\\theta}$ is a diagonal unitary, and the $L$ and $R$ subcircuits can be found by our search-based synthesis algorithm. We relax the need for our search-based synthesis algorithm to find sequences of gates that implement our target to just finding gates that diagonalize our target."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8f6130d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from bqskit.compiler import Compiler\n",
    "from bqskit.ir.circuit import Circuit\n",
    "from bqskit.passes.synthesis import QSearchSynthesisPass\n",
    "\n",
    "from utils.diagonal_synth import replace_pauliz\n",
    "from utils.discrete import DiscreteLayerGenerator\n",
    "from utils.printing import print_circuit\n",
    "\n",
    "success_threshold = 1e-8\n",
    "# double_headed=True means we will try to diagonalize the target\n",
    "layer_generator = DiscreteLayerGenerator(double_headed=True)\n",
    "synthesis_pass = QSearchSynthesisPass(\n",
    "    success_threshold=success_threshold,\n",
    "    layer_generator=layer_generator,\n",
    ")\n",
    "\n",
    "circuit = Circuit.from_unitary(hard_target)\n",
    "\n",
    "# Diagonalize\n",
    "with Compiler() as compiler:\n",
    "    circuit = compiler.compile(circuit, [synthesis_pass])\n",
    "\n",
    "# Decompose diagonal component\n",
    "circuit = replace_pauliz(circuit)\n",
    "\n",
    "print_circuit(circuit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "263e62ce",
   "metadata": {},
   "source": [
    "__Question:__ What similarities do you see between the outputs of synthesis for `easy_target` and `hard_target`?\n",
    "\n",
    "__Question:__ Assume that the best implementation of `hard_target` corresponds to the results from above, but with the $R_Z$ operation converted to ${H,T}$ using gridsynth. Approximately how many candidate circuits might need to be searched to find this circuit without diagonalization? Assume `success_threshold=1e-8`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fb8f283",
   "metadata": {},
   "source": [
    "It turns out that there are many unitaries that appear in real quantum algorithms which are amenable to synthesis by Clifford+T diagonalization. Check out [9] if you'd like to learn about how we trained Reinforcement Learning agents to do exactly this!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "278394b3",
   "metadata": {},
   "source": [
    "### Synthesis with Constrained Optimization\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63df7524",
   "metadata": {},
   "source": [
    "For NISQ style circuit synthesis, BQSKit uses numerical instantiation to solve for the circuit parameters that approximate a target unitary. Using constrained optimization, we can apply that same principle here as well.\n",
    "\n",
    "Using `optimization_level>=2` or providing a unitary as the target in `bqskit.compile` when targeting a `FaultTolerantModel` will automatically select this work flow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "80b51fa6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from bqskit.qis.unitary import UnitaryMatrix\n",
    "from numpy import array\n",
    "\n",
    "hard_target = UnitaryMatrix(array([\n",
    "   [0.70370187-0.06930858j, 0, 0.70370187-0.06930858j, 0],\n",
    "   [0, 0.70370187+0.06930858j, 0, 0.70370187+0.06930858j],\n",
    "   [0.70370187+0.06930858j, 0, -0.70370187-0.06930858j, 0],\n",
    "   [0, 0.70370187-0.06930858j, 0, -0.70370187+0.06930858j],\n",
    "]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86636ee2",
   "metadata": {},
   "outputs": [],
   "source": [
    "from bqskit.compiler.compile import compile\n",
    "\n",
    "from bqskit.ft.cliffordt.cliffordtmodel import CliffordTModel\n",
    "\n",
    "# Set up the machine model and synthesize the target unitary\n",
    "machine = CliffordTModel(num_qudits=hard_target.num_qudits)\n",
    "result = compile(hard_target, model=machine, optimization_level=3)\n",
    "\n",
    "# Print out the resulting circuit\n",
    "print_circuit(result)\n",
    "\n",
    "dist = result.get_unitary().get_distance_from(hard_target)\n",
    "\n",
    "from bqskit.ir.gates.constant.identity import IdentityGate\n",
    "ft_gateset = list(machine.gate_set) + [IdentityGate()]\n",
    "\n",
    "assert dist < 5e-8\n",
    "assert all(op.gate in ft_gateset for op in result)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e76773a9",
   "metadata": {},
   "source": [
    "Using numerical instantiation can result in extra Clifford gates compared to pure search-based approaches. Luckily, optimization of Clifford circuits is a well studied field, and Clifford circuits are relatively cheap to implement. Numerical instantiation based flows may be the only option for finding circuit implementations of complex targets. \n",
    "\n",
    "For details as to how contrained numerical optimization plays into instantiation, see [10]."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f07163a",
   "metadata": {},
   "source": [
    "## References\n",
    "-------------\n",
    "[1] https://arxiv.org/abs/2312.03982, D. Bluvstein, et al., \"Logical quantum processor based on reconfigurable atom arrays\"\n",
    "<br>[2] https://arxiv.org/abs/2404.02280, M. P. da Silva, et al., \"Demonstration of logical qubits and repeated error correction with better-than-physical error rates\"\n",
    "<br>[3] https://journals.aps.org/pra/abstract/10.1103/PhysRevA.52.R2493, P. Shor, \"Scheme for reducing decoherence in quantum computer memory\"\n",
    "<br>[4] https://arxiv.org/abs/quant-ph/9601029, A. Steane, \"Multiple Particle Interference and Quantum Error Correction\"\n",
    "<br>[5] https://arxiv.org/abs/quant-ph/9811052, S. Bravyi, A. Kitaev, \"Quantum codes on a lattice with boundary\"\n",
    "<br>[6] https://iopscience.iop.org/article/10.1070/RM1997v052n06ABEH002155/meta, A. Kitaev, \"Quantum computations: algorithms and error correction\"\n",
    "<br>[7] https://arxiv.org/abs/quant-ph/0403025, S. Bravyi, A. Kitaev, \"Universal Quantum Computation with ideal Clifford gates and noisy ancillas\"\n",
    "<br>[8] https://arxiv.org/abs/1403.2975, N. Ross, P. Selinger, \"Optimal ancilla-free Clifford+T approximations of z-rotations\"\n",
    "<br>[9] https://arxiv.org/abs/2409.00433, M. Weiden, et al., \"Learning to Synthesize Fault-Tolerant Quantum Circuits by Diagonalization\"\n",
    "<br>[10] https://marcdav.is/static/fileshare/Numerical_Synthesis_of_Arbitrary-Multi-Qubit_Unitaries_with_Low_T-Count.pdf, M Grau Davis, \"Numerical Synthesis of Arbitrary Multi-Qubit Unitaries with low T-Count\""
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}