[Feature]: two classical simulation algorithms for Clifford+T circuits

[Fe-r-oz]

Is your feature request related to a problem? Please describe.

The paper titled Improved classical simulation of quantum circuits dominated by Clifford gates[1] by Sergey and David extends the Gottesman-Knill theorem by providing two classical simulation algorithms for $Clifford+T$ quantum circuits by introducing new algorithmic tools for stabilizer simulators in the following three verticals: gadgetized circuits (for decomposing magic states), stabilizer rank of magic states, and fast norm estimation algorithm.

A $Clifford+T$ quantum circuit of length $m$ operating on $n$ qubits is described by the unitary operator $U = U_m \cdots U_2 U_1$, where each $U_j$ is a one- or two-qubit gate from the set ${H, S, T, \text{CNOT}}$. In this set, $S = T^2$. Let $m = c + t$, where $c$ denotes the number of Clifford gates ($H, S, \text{CNOT}$) and $t$ denotes the number of T-gates, referred to as the T-count. When $U$ is applied to the initial state $|0^n\rangle$ and a specific output register $Q_{\text{out}} \subseteq [n]$ is measured in the 0,1-basis, a random bit string $x$ of length $w = |Q_{\text{out}}|$ is produced [1].

CREDITS: Sergey Bravyi: Improved classical simulation of quantum circuits dominated by Clifford gates

The probability of obtaining the bit string $x$ is given by $P_{\text{out}}(x) = \langle 0^n | U^\dagger \Pi(x) U | 0^n \rangle$, where $\Pi(x)$ projects $Q_{\text{out}}$ onto the basis state $|x\rangle$ and acts trivially on the other qubits [1]. The magic state is $|A\rangle = 2^{-1/2} (|0\rangle + e^{i\pi/4} |1\rangle)$.

Describe the solution you’d like

The first step is to check out the talk by Sergey on his paper: Sergey Bravyi: Improved classical simulation of quantum circuits dominated by Clifford gates

They introduce two simulation tasks, each aligned with a classical simulation algorithm:

1. Strong Simulation: Computes the probability of a given output $x$.

• Task: Approximate the probability $P_{\text{out}}(x)$ for a given string $x \in [0, 1]^w$ with a specified relative error $\epsilon$ and a failure probability $p_f$. Context: Consider $Clifford+T$ circuit with an output register of $w$ qubits. Measuring each output qubit produces a string $x$ drawn from a probability distribution $P_{\text{out}}$. Note: $n$: number of qubits, $c$ is number of clifford gates, $t$ is number of $T$ gates, $w$ is the output width.
• Runtime: $\text{poly}(n, c, t) + 2^{0.47t} \cdot t^3 \epsilon^{-2}$.
• Remark by Sergey: Strong simulation is difficult (and unreasonable) even for quantum computers, as they can only sample from the probability distribution, not approximate output probabilities.

2. Weak Simulation: Samples $x$ from the output probability distribution $P_{\text{out}}(x, y)$.

• Task: Sample the output string $x$ from a distribution that is $\epsilon$-close to $P_{\text{out}}$ with respect to the $L_1$-norm.
• Runtime: $\text{poly}(n, c, t) + 2^{0.23t} \cdot t^3 w^4 \epsilon^{-5}$.
• Parallelism and Implementation: Both algorithms are divided into independent subroutines as presented in their paper with a runtime of $O(t^3)$ each, allowing for significant parallelism. Pseudocode for the main subroutines (5 in total) is provided in the paper. Specifically, the simulation addressed a quantum algorithm for solving the hidden shift problem for non-linear Boolean functions.
• The classical sampling algorithm was implemented in MATLAB and used to simulate benchmark quantum circuits with:
  • $n = 40$ qubits,
  • Several hundred Clifford gates,
  • T-count $t \leq 48$.
• Sergey notes in his talk that it took them 5 CPU hours to simulate $t \approx 50$, $n \approx 50$, $c \approx 1000$, $w = 1$, $\epsilon \approx 10%$ for weak simulation.

There are basically three steps and one preprocessing step (each colored differently) for classicaly simulating the $Clifford+T$ circuit. The Gottesman-Knill preprocessing step that Sergey discusses (at 25:07 minute mark) in his talk is in blue.

The following implementation flow is consistent with both Sergey’s paper and his talk.

There are a lot of small subtasks/low hanging fruits that spawn out from this paper on the road the implementing these two classical simulation algorithms for $Clifford+T$ circuits.

For instance, some of these subtasks are:

• functionality to implement T-gadget so that we can replace T gate with T-gadget.
  
  CREDITS: Improved classical simulation of quantum circuits dominated by Clifford gates
• functionality to calculate the stabilizer rank $\chi$, given a stabilizer state. The stabilizer rank of a quantum state $|\psi\rangle$ is the minimal $\chi$ such that $\psi$ can be written as a linear combination of $\chi$ stabilizer states:

$$|\psi\rangle = \sum_{a=1}^{\chi}c_{a}|\psi_{a}\rangle$$

for $c_{a} \in \mathbb{C}$ and stabilizer states $|\psi_{a}\rangle$.  

• In the Appendix A, Sergey makes connection to the quadratic forms and how it is used to represent the stabilizer state uniquely.

According to Lemma 4, equation (47) of Sergey's paper, Any $n$-qubit stabilizer state $|\psi\rangle$ can be uniquely expressed as:

$$|\psi\rangle = 2^{-\frac{k}{2}} \sum_{x \in K} e^{\frac{i\pi}{4} q(x)} |x\rangle$$

where $K \subseteq \mathbb{F}_2^n$ is an affine subspace of dimension $0 \leq k \leq n$ and $q: K \rightarrow \mathbb{Z}_8$ is a quadratic form. In the subroutines presented in appendix of his paper, he relies on this representation. Exploring how this can be implemented with help of Nemo would be nice.

Describe alternatives you’ve considered

In progress.

Additional context

The authors introduced two algorithms for the classical simulation of quantum circuits that utilize the Clifford+T gate set. The runtime for these algorithms scales polynomially with the number of qubits and the number of Clifford gates in the circuit but exponentially with the number of T gates. Despite the exponential factor, this scaling remains 'mild', suggesting that classical simulations of Clifford+T circuits with several hundred qubits and a T-count $t \leq 50$ can be executed on a medium-sized computer cluster [1].

Please consider reading through Sergey's second paper on $Clifford+T$ which he refers when talking about "sketch of techniques section" in the aforementioned paper. The second paper is Trading classical and quantum computational resources. Additionally, please examine the difference between the RandomStabilizerState subroutine and the corresponding implementation (if) provided here. Specifically, is the primary difference related to the use of quadratic forms in representing the stabilizer state use by the subroutine RandomStabilizerState ?