home/blog/efficient-sim

Efficient simulation of quantum circuits

2024.05.17

Here’s an interesting tidbit from my research, concerning different approaches toward simulating systems of qubits in the context of measurement-induced phase transitions. There’s this idea that for a certain class of quantum states and a certain set of transformations one can apply to those states, classical simulation of the system becomes feasible, even for relatively large system sizes. All this stuff has become well-known over the past couple of decades or so, but like a lot of other things I do, I thought I’d take a stab at creating my own simulator for this kind of system. My approach turned out not to be efficient (although it has some nice pedagogical properties), so here’s a high-level overview of three approaches to classical simulation: one which is obviously and definitely not efficient, one which seems like it could be efficient but isn’t (my approach), and the one that everyone uses because it actually is efficient.

Introduction/Motivation

Measurement-induced phase transitions have been a recent focus of computational quantum information science. By applying alternating layers of unitary rotation and projective measurement at variable rates to a system of qubits, one can observe a phase transition in the rate at which the mount of entanglement grows with the size of the system, versus that rate at which projective measurements occur. We consider circuits that look like this:

where each can be any -qubit unitary and measurements are single-qubit only, occurring with fixed probability . Commonly the simulated circuits look more like this:

where each can be any single-qubit unitary, and each can be any two-qubit unitary applied to alternating left- and right-nearest neighbors on each layer of the circuit. (Usually these systems are treated in 1D for simplicity.)

Intuitively, this kind of circuit is like a Trotterized version of what one can imagine happening in real-world quantum systems: We take some initial state, apply random dynamics and some probability of collapsing superpositions due to interaction with some background. Randomized dynamics tend to increase the amount of entanglement in the system, while the measurements decrease it. In particular, the entropy of entanglement in the system is an extrinsic quantity with some relationship to the system size . The idea is that there exists some critical measurement probability , below which the entanglement entropy grows as a volume law and above which it grows as an area law. For 1D, gives , gives , and gives . This is the phase transition.

One of the reasons this is an attractive system to simulate is that there are few constraints on what exactly needs to happen in these circuits in order to see a phase transition, only that there needs to be alternation between unitary dynamics and probabilistic projective measurements. The Gottesman-Knill theorem¹ guarantees that if you limit your system to stabilizer states and Clifford gates (i.e. the six cardinal directions on the Bloch sphere and only gates generated by Hadamard, Z(π/2), and CNOT), the system can be efficiently simulated (sub-exponential scaling in the number of qubits in both time and space) on a classical computer.

Approach 1: State vectors and matrices

So here’s the naive approach to this: Simply represent the states and gates in the standard (Z) basis as complex-valued vectors and matrices. The application of gates and measurements is just (possibly sums over) matrix-vector products, which are well-understood and ubiquitous, meaning there are plenty of libraries one could use for the actual computation. Since all of the math describing this stuff is already phrased using linear algebra anyway, the conceptual overhead in this approach is quite low.

The problem, of course, is that the number of complex values needed to describe an n-qubit state is precisely 2ⁿ. Hence, gates are represented by 2ⁿ×2ⁿ matrices, and the simulation quickly explodes in both time and space.

Approach 2: Trees

The zero-th order correction to the vector/matrix approach is to not represent the states and gates with complex numbers. Even disregarding the sizes of the vectors and matrices, complex numbers model a continuous space and if we’re limited to only six positions on the Bloch sphere, then we definitely don’t need all those bits. The basic motivation for this representation comes from these observations:

1. Single-qubit gates only apply rotations to single qubits (pretty obvious), so an n-qubit register should be represented as a list of n Bloch sphere positions to obviate the need to represent, e.g. the state with two amplitudes. By treating all Bloch sphere positions as base objects, we could, for example, have the Hadamard gate just convert a to a .
2. Two-qubit gates other than SWAP will in some cases also only apply rotations to one qubit, but in others they’ll create irreducible superpositions of n-qubit states (these are the cases where they create entanglement in the form of Bell states). However, it can be seen that in these cases, a basis (X, Y, or Z) can always be chosen such that the resulting state can be written down as the superposition of at most two basis states.

Hence, states in this approach are modeled as quasi-binary trees, where the leaves are single basis states, and branching denotes a superposition. I think in the case where you always start from a fixed single basis sate, the tree will always be completely balanced, but what follows below will handle the general case.

We can start off with some basic definitions:

-- A qubit is one of the six cardinal directions.
data Qubit = Xp -- ∣+x⟩
           | Xm -- ∣-x⟩
           | Yp -- ∣+y⟩
           | Ym -- ∣-y⟩
           | Zp -- ∣+z⟩
           | Zm -- ∣-z⟩
           deriving (Eq)

-- Phase angles can be limited to integer multiples of π/4.
data Phase = Pi0  -- 0π
           | Pi1q -- π/4
           | Pi1h -- π/2
           | Pi3q -- 3π/4
           | Pi   -- π
           | Pi5q -- 5π/4
           | Pi3h -- 3π/2
           | Pi7q -- 7π/4

-- Define the gate set: each possible gate carries one or two qubit indices.
data Gate = H Int -- Hadamard
          | X Int -- π rotation about X
          | Y Int -- π rotation about Y
          | Z Int -- π rotation about Z
          | S Int -- π/2 rotation about Z
          | CX Int Int -- Z-controlled π rotation about X
          | Swap Int Int -- swap two qubits

-- A register (basis) state is just a list of `Qubit`s.
type Register = [Qubit]

-- Create a new register of qubits initialized to all `Zp`.
regNew :: Int -> Register
regNew = flip replicate Zp

This definition of a qubit is a variable-basis representation. That is, the basis of the representation can change to suit the state. In principle, we only really need two of the six to describe a state, but adding in the four others allows us to reduce a two-term, single-qubit superposition to just one term.

Recursive data structures are a pleasure in high-level languages, so let’s go with that to define a general pure state:

-- A pure state of an n-qubit register.
data Pure = Null -- a physically impossible null state, useful for measurements
          | Single Register -- only a single basis state
          | Superpos Pure Pure Phase -- an even superposition of two basis
                                     -- states with a relative phase

pureNew :: Int -> Pure
pureNew = Single . regNew

A Pure is a quasi-binary tree in the sense that every node in the tree is either Null/Single or a Superpos with two children, and there are no restrictions on the locations of either kind of node. The Phase attached to a Superpos is the phase of the Pure on the right relative to the Pure on the left, and the magnitude of the amplitude of a given Single is , where is its depth in the tree. So the pure state

has a tree that looks like this:

and the corresponding structure in Haskell is:

Superpos
  ( Superpos
    (Single [Zp, Zp, Zp])
    (Single [Zm, Zp, Zm])
    Pi0
  )
  ( Superpos
    (Single [Zp, Zm, Zp])
    (Single [Zm, Zm, Zm])
    Pi
  )
  Pi0

From here, operations can be performed on the tree by traversing it. Single-qubit rotations can be applied to modify the appropriate qubit in a Single, do nothing on a Null, or else be applied to the left and right branches of a Superpos:

pureApplyGate :: Gate -> Pure -> Pure
pureApplyGate gate state = fst $ pureApplyGateInner gate state
  where pureApplyGateInner :: Gate -> Pure -> (Pure, Maybe Phase)
        pureApplyGateInner _ Null = (Null, Nothing)
        pureApplyGateInner g (Single reg) = case g of -- gate action...
        pureApplyGateInner g (Superpos l r ph) =
          let (l', maybePhL) = pureApplyGateInner g l
              (r', maybePhR) = pureApplyGateInner g r
           in case (maybePhL, maybePhR) of
                -- handle null cases and propagate phases...

If a phase is incurred on a Single, it should be propagated back up the tree to modify the relative phases of relevant superpositions. For the recursive call in the Superpos case, phases coming from the right branch should be absorbed into the relative phase of the superposition, and phases coming from the left branch should be propagated. The total update to the relative phase is . Any phases that make it up to the root of the tree are discarded.

CNOT gates can be handled by either applying the appropriate single-qubit rotations or, in the case of a Bell state output, converting a Single to a Superpos. Here, we want to think about what the CNOT gate does in different bases, which can be seen by considering the matrix

for all with being the appropriate basis transformation. Here is the usual matrix in the Z-basis,

and changes two-qubit bases from to ,

Of the nine possible combinations for the input state bases (), the bases , , , and create Bell states of various symmetries and phases; the rest apply a rotation to either the control or target qubits and/or a phase.

Measurement here is most easily performed with the understanding that it’s a post-selected measurement (i.e. the measurement was applied, and the output state is what you’d get by filtering for a desired outcome; this will make more sense in the context of mixed states below). For a measurement on a single qubit in a given basis, the qubit state can either align with the measurement basis or not. In the former case, nothing needs to be done unless the post-selection is incompatible with the state of the targeted qubit, in which case nodes can be replaced with Null as needed. In the latter case, we just have to set the qubit to agree with the post-selection.

-- Perform a projective measurement on a single qubit, post-selected to a
-- particular outcome.
pureMeasure :: Int -> Qubit -> Pure -> Pure
pureMeasure idx outcome state = fst $ pureMeasureInner idx outcome state
  where pureMeasureInner :: Int -> Qubit -> Pure -> (Pure, Maybe Phase)
        pureMeasureInner _ _ Null = (Null, Nothing)
        pureMeasureInner idx outcome (Single reg) = -- set the qubit state...
        pureMeasureInner idx outcome (Superpos l r ph) =
          let (l', maybePhL) = pureMeasureInner idx outcome l
              (r', maybePhr) = pureMeasureInner idx outcome r
           in case (maybePhL, maybePhR) of
                -- handle null cases and propagate phases...

This handling of measurement leads naturally into mixed states, where the non-aligning case produces an even classical mixture of two pureMeasure outputs, one post-selected to the plus state of the bases, and one to the minus state. Thus we have another tree representation:

-- A classical mixture of pure states.
data State = Empty -- a physically impossible null state
           | Pure Pure -- a single `Pure` state
           | Mixed State State -- a mixed state comprising the even classical
                               -- mixture of two pure states

-- Create a new state consisting of only the default pure state
stateNew :: Int -> State
stateNew = Pure . pureNew

The probability of a given pure state in the classical mixture is again related to its depth in the tree, but here we obviously don’t have to keep track of phases. Gates are applied by recursively traversing the tree,

stateApplyGate :: Gate -> State -> State
stateApplyGate _ Empty = Empty
stateApplyGate g (Pure pure) = Pure $ pureApplyGate g pure
stateApplyGate g (Mixed l r) = Mixed l' r'
  where l' = stateApplyGate g l
        r' = stateApplyGate g r

and measurements behave similarly to CNOTs on Pures:

-- The basis for a given measurement
data Basis = MeasX | MeasY | MeasZ

-- Return the (plus, minus) states of a basis.
basisOutcomes :: Basis -> (Qubit, Qubit)
basisOutcomes MeasX = (Xp, Xm)
basisOutcomes MeasY = (Yp, Ym)
basisOutcomes MeasZ = (Zp, Zm)

-- Perform a projective measurement on a single qubit in a given basis.
stateMeasure :: Int -> Basis -> State -> State
stateMeasure _ _ Empty = Empty
stateMeasure idx basis (Pure pure) =
  let (op, om) = basisOutcomes basis
      pureP = pureMeasure idx op pure
      pureM = pureMeasure idx om pure
   in -- handle null cases...
stateMeasure idx basis (Mixed l r) =
  let l' = stateMeasure idx basis l
      r' = stateMeasure idx basis r
   in -- handle null cases...

The reason why these trees are not efficient lies in their recursive nature. Although register states are not properly state vectors and gates are not applied via matrix multiplication, the overall tree will still grow on average when gates are applied to a Pure. (More precisely, a given CNOT gate will have, from a coarse average, about 4/9 probability of increasing the height of the tree uniformly.) Only post-selected measurements can shrink the tree, so for cases where the measurement rate is less than the rate at which CNOTs are applied, the size of a given Pure will grow basically exponentially in the depth of the circuit. In principle, there should be many cases where the phases in the tree will be such that some terms should cancel, but it’s rather hard in this representation to capture that behavior. The exponential growth of the tree then also exponentially increases the number of operations needed to apply gates. When States are used, whole Pures are then effectively duplicated, leading to further exponential growth.

(CS people are probably rolling their eyes or laughing up their sleeves right now at the notion that a tree could be efficient in this context…)

Approach 3: Gottesman-Knill tableaus

So here’s the approach that works. The Gottesman-Knill (GK) tableau is a representation of stabilizer states that is directly derived from the GK theorem. To be more precise about it compared to its brief mention above, the theorem says that stabilizer states can be identified by the amplitudes of basis states, but by their stabilizer group (here, “group” is defined in the mathematical sense) intersected with the group formed by all n-qubit Pauli operators. That is, the set of all operators formed as the tensor product of n Pauli operators (with a sign) of which a given state is a +1 eigenstate. This set for a given stabilizer state uniquely identifies , meaning that (given another stabilizer state ) .

The reason why this allows for efficient simulation goes like this: Although the size of this group, , is exponential in the number of qubits (in fact, ), there is a well-known theorem from group theory that says any group is generated by at most elements. Therefore the generating group of -qubit Pauli stabilizers for a given state is of size . Since there are four Pauli operators, each operator can be represented using only two bits, so each -qubit Pauli stabilizer in this generating group (with its sign) requires bits. Hence, any given stabilizer state can be represented with only bits.

GK shows that the application of the usual Clifford generator gates can be understood as update rules on these bits, each having only time complexity. Measurements have similar update rules (with the understanding that in this scheme, measurement is only in the Z-basis and is allowed to be probabilistic), but these are due to what is essentially matrix inversion. In a later paper with Scott Aaronson², this was reduced to at the cost of storing an additional bits to identify a generating set of what they call the “destabilizer” group for the state (which, together with the generating set for the stabilizer group, generate the complete -qubit Pauli group). Most recently, Craig Gidney from Google Quantum further reduced this to by additionally storing the inverse of the resulting tableau (described below)³.

In total, the bits describing the (de)stabilizer generators are assembled into a binary matrix conventionally called a tableau:

where the top half encodes the destabilizer generators and the bottom half encodes the stabilizer generators. The -th Pauli of the -th (de)stabilizer is identified by the and entries as

and the value of (for destabilizers, ) gives the sign. This scheme produces time complexity, and allows stabilizer circuits for even qubits to be simulated on just my laptop!

As an aside, there’s also a nice connection to graph states here. Graph states form a subset of stabilizer states that can be drawn as graphs, where nodes represent qubits in and edges represent the application of CZ gates. When a CZ gate is applied between two qubits, it adds an edge between their nodes if it doesn’t exist, or removes it if it does. The connection is that every graph state is stabilized by a set of stabilizers where the -th stabilizer is constructed from the following rules:

1. The Pauli on the -th qubit is
2. The Paulis on each qubit with an edge to the -th qubit are all
3. The Paulis applied to all other qubits are all

In the tableau representation, this corresponds to a stabilizer block that is the identity in the portion and the adjacency matrix of the graph in the portion. And since every graph state can be formed by taking an initial state (whose tableau is the identity matrix in the destabilizer X and stabilizer Z blocks and zeros everywhere else), applying Hadamards to every qubit, and then CZs according to the edges in the graph, the destabilizer portion of the tableau is just all zeros in the X portion and the identity in the Z portion. Hence, there’s an easy algorithm to convert a graph state to a full tableau which is just traversing the graph. Additionally, Vijayan et al.⁴ have a algorithm (based on Gaussian elimination and other operations equivalent to applying Clifford gates) to convert an arbitrary stabilizer state into a graph state. The fact that such an algorithm exists is what’s known as the local equivalence between stabilizer states and graph states.

Footnotes

1. ↵ D. Gottesman, “The Heisenberg representation of quantum computers,” arXiv:9807006.
2. ↵ S. Aaronson and D. Gottesman, “Improved simulation of stabilizer circuits,” arXiv:0406196.
3. ↵ C. Gidney, “Stim: a fast stabilizer circuit simulator,” Quantum 5, 497 (2021).
4. ↵ M. K. Vijayan et al., “Compilation of algorithm-specific graph states for quantum circuits,” arXiv:2209.07345.

Changelog

• Originally posted 2024.05.17
• Updated in Typst reformat 2026.03.17 with minor wording adjustments