"""
Stabilizer Backend for Efficient Clifford Circuit Simulation
Implements the Gottesman-Knill theorem: Clifford circuits can be simulated
efficiently (polynomial time/space) using stabilizer tableaux.
Key features:
- O(n²) time per gate vs O(2ⁿ) for state-vector
- Supports: H, S, CNOT, CZ, SWAP, measurements
- Handoff to MPS when non-Clifford gates appear (T, Toffoli, etc.)
References:
- Gottesman (1998): "The Heisenberg Representation of Quantum Computers"
- Aaronson & Gottesman (2004): "Improved Simulation of Stabilizer Circuits"
Author: ATLAS-Q Contributors
Date: October 2025
License: MIT
"""
from dataclasses import dataclass
from typing import List, Optional
import numpy as np
import torch
[docs]
@dataclass
class StabilizerState:
"""
Stabilizer state represented as a tableau
A stabilizer state on n qubits is represented by a (2n+1) × 2n binary matrix:
- First n rows: X part of stabilizer generators
- Next n rows: Z part of stabilizer generators
- Last column: phase bits (±1)
Example for |0⟩: stabilizer is Z
Example for |+⟩: stabilizer is X
Example for Bell |Φ+⟩: stabilizers are XX and ZZ
"""
n_qubits: int
tableau: np.ndarray # Shape: (2n, 2n+1)
# tableau[:n, :] = X parts, tableau[n:, :] = Z parts
# tableau[:, -1] = phases
def __post_init__(self):
assert self.tableau.shape == (2 * self.n_qubits, 2 * self.n_qubits + 1)
[docs]
@staticmethod
def init_zero(n_qubits: int) -> "StabilizerState":
"""Initialize |00...0⟩ state (stabilized by Z on each qubit)"""
tableau = np.zeros((2 * n_qubits, 2 * n_qubits + 1), dtype=np.uint8)
# Set destabilizers to X_i for each qubit i (rows 0:n, columns 0:n)
for i in range(n_qubits):
tableau[i, i] = 1
# Set stabilizers to Z_i for each qubit i (rows n:2n, columns n:2n)
for i in range(n_qubits):
tableau[n_qubits + i, n_qubits + i] = 1
return StabilizerState(n_qubits, tableau)
[docs]
@staticmethod
def init_plus(n_qubits: int) -> "StabilizerState":
"""Initialize |++...+⟩ state (stabilized by X on each qubit)"""
tableau = np.zeros((2 * n_qubits, 2 * n_qubits + 1), dtype=np.uint8)
# Set stabilizers to X_i for each qubit i
for i in range(n_qubits):
tableau[i, i] = 1
return StabilizerState(n_qubits, tableau)
[docs]
def copy(self) -> "StabilizerState":
"""Create a copy of this state"""
return StabilizerState(self.n_qubits, self.tableau.copy())
[docs]
class StabilizerSimulator:
"""
Efficient simulator for Clifford circuits using stabilizer formalism
Complexity:
- Space: O(n²) vs O(2ⁿ) for state-vector
- Time per gate: O(n²) vs O(2ⁿ)
Usage:
sim = StabilizerSimulator(n_qubits=100)
sim.h(0)
sim.cnot(0, 1)
outcome = sim.measure(0)
"""
[docs]
def __init__(self, n_qubits: int):
self.n_qubits = n_qubits
self.state = StabilizerState.init_zero(n_qubits)
self.measurement_outcomes: List[int] = []
[docs]
def copy(self) -> "StabilizerSimulator":
"""Fast copy using numpy array copy (much faster than deepcopy)"""
new_sim = StabilizerSimulator.__new__(StabilizerSimulator)
new_sim.n_qubits = self.n_qubits
new_sim.state = self.state.copy()
new_sim.measurement_outcomes = self.measurement_outcomes.copy()
return new_sim
def _rowsum(self, h: int, i: int):
"""
Multiply stabilizer generator h by generator i (Pauli group multiplication)
This implements the Pauli commutation relations:
XY = iZ, YZ = iX, ZX = iY
and tracks the phase
"""
n = self.n_qubits
tab = self.state.tableau
# Extract X and Z parts
x_h = tab[h, :n]
z_h = tab[h, n : 2 * n]
x_i = tab[i, :n]
z_i = tab[i, n : 2 * n]
# Phase update according to Pauli multiplication rules
phase_update = 0
for j in range(n):
if x_h[j] and z_h[j]: # Y
if x_i[j] and not z_i[j]: # X
phase_update += 1 # YX = -iZ → phase flip
elif not x_i[j] and z_i[j]: # Z
phase_update += 3 # YZ = iX → 3 phase flips = -i
elif x_h[j] and not z_h[j]: # X
if x_i[j] and z_i[j]: # Y
phase_update += 3 # XY = iZ → 3 phase flips
elif not x_i[j] and z_i[j]: # Z
phase_update += 1 # XZ = -iY
elif not x_h[j] and z_h[j]: # Z
if x_i[j] and z_i[j]: # Y
phase_update += 1 # ZY = -iX
elif x_i[j] and not z_i[j]: # X
phase_update += 3 # ZX = iY
# Update tableau
tab[h, :n] = (x_h + x_i) % 2
tab[h, n : 2 * n] = (z_h + z_i) % 2
tab[h, -1] = (tab[h, -1] + tab[i, -1] + (phase_update // 2)) % 2
# Clifford gates
[docs]
def h(self, qubit: int):
"""Hadamard gate: X ↔ Z"""
n = self.n_qubits
tab = self.state.tableau
for i in range(2 * n):
x_bit = tab[i, qubit]
z_bit = tab[i, n + qubit]
# H: X ↔ Z, phase changes if both X and Z are present
if x_bit and z_bit:
tab[i, -1] = (tab[i, -1] + 1) % 2 # Phase flip for Y → -Y
# Swap X and Z
tab[i, qubit] = z_bit
tab[i, n + qubit] = x_bit
[docs]
def s(self, qubit: int):
"""Phase gate: X → Y, Z → Z"""
n = self.n_qubits
tab = self.state.tableau
for i in range(2 * n):
x_bit = tab[i, qubit]
z_bit = tab[i, n + qubit]
# S transforms:
# X → Y: add Z bit, flip phase
# Y → -X: remove Z bit, flip phase
# Z → Z: unchanged
if x_bit and not z_bit:
# X → Y
tab[i, n + qubit] = 1
tab[i, -1] = (tab[i, -1] + 1) % 2
elif x_bit and z_bit:
# Y → -X
tab[i, n + qubit] = 0
tab[i, -1] = (tab[i, -1] + 1) % 2
[docs]
def s_dag(self, qubit: int):
"""Inverse phase gate: X → -Y, Z → Z"""
# S†: Apply S three times (since S⁴ = I)
self.s(qubit)
self.s(qubit)
self.s(qubit)
[docs]
def cnot(self, control: int, target: int):
"""CNOT gate"""
n = self.n_qubits
tab = self.state.tableau
for i in range(2 * n):
# CNOT rules:
# X_c → X_c X_t
# Z_c → Z_c
# X_t → X_t
# Z_t → Z_c Z_t
x_c = tab[i, control]
z_c = tab[i, n + control]
x_t = tab[i, target]
z_t = tab[i, n + target]
# Phase: anticommutes if X_c Z_t but not (Z_c or X_t)
if x_c and z_t and not (z_c or x_t):
tab[i, -1] = (tab[i, -1] + 1) % 2
# Update X_t and Z_c
tab[i, target] = (x_t + x_c) % 2
tab[i, n + control] = (z_c + z_t) % 2
[docs]
def cz(self, qubit1: int, qubit2: int):
"""CZ gate (symmetric)"""
# CZ = H₂ CNOT H₂
self.h(qubit2)
self.cnot(qubit1, qubit2)
self.h(qubit2)
[docs]
def swap(self, qubit1: int, qubit2: int):
"""SWAP gate"""
# SWAP = CNOT₁₂ CNOT₂₁ CNOT₁₂
self.cnot(qubit1, qubit2)
self.cnot(qubit2, qubit1)
self.cnot(qubit1, qubit2)
[docs]
def x(self, qubit: int):
"""Pauli-X gate"""
# X flips Z ↔ -Z in stabilizers
n = self.n_qubits
tab = self.state.tableau
for i in range(2 * n):
# If the row has Z on this qubit, flip the phase
if tab[i, n + qubit] == 1:
tab[i, -1] = (tab[i, -1] + 1) % 2
[docs]
def y(self, qubit: int):
"""Pauli-Y gate"""
# Y flips both X and Z in stabilizers
n = self.n_qubits
tab = self.state.tableau
for i in range(2 * n):
# If the row has X or Z (or both) on this qubit, flip the phase
if tab[i, qubit] == 1 or tab[i, n + qubit] == 1:
tab[i, -1] = (tab[i, -1] + 1) % 2
[docs]
def z(self, qubit: int):
"""Pauli-Z gate"""
# Z flips X ↔ -X in stabilizers
n = self.n_qubits
tab = self.state.tableau
for i in range(2 * n):
# If the row has X on this qubit, flip the phase
if tab[i, qubit] == 1:
tab[i, -1] = (tab[i, -1] + 1) % 2
# Measurement
[docs]
def measure(self, qubit: int, rng: Optional[np.random.RandomState] = None) -> int:
"""
Measure qubit in computational basis
Returns:
0 or 1 (measurement outcome)
"""
if rng is None:
rng = np.random.RandomState()
n = self.n_qubits
tab = self.state.tableau
# Check if any stabilizer has X on this qubit
# If so, measurement outcome is random
x_column = tab[n:2*n, qubit]
# Find first stabilizer with X on this qubit
p = -1
for i in range(n):
if x_column[i] == 1:
p = n + i # Actual row index in tableau
break
if p == -1:
# Deterministic outcome
# Find which stabilizer row has Z (but not X) on this qubit
outcome = 0 # default
for i in range(n, 2 * n): # Only search stabilizer rows
if tab[i, qubit] == 0 and tab[i, n + qubit] == 1:
# This row has Z but not X on this qubit
# Phase bit encodes parity of ALL qubits with Z in this stabilizer
# We need to account for other qubits already measured
phase = int(tab[i, -1])
# XOR with measurement outcomes of all OTHER qubits with Z in this stabilizer
for j in range(n):
if j != qubit and tab[i, n + j] == 1: # Other qubit has Z
# Get measurement outcome of qubit j
# Check if qubit j has been measured (has Z but not X in some stabilizer)
measured_j = 0
for k in range(n, 2 * n):
if tab[k, j] == 0 and tab[k, n + j] == 1 and tab[k, n + qubit] == 0:
# Stabilizer k has Z on j but not on our target qubit
# This means j was measured with outcome = phase of k
measured_j = int(tab[k, -1])
break
phase = (phase + measured_j) % 2
outcome = phase
# Update tableau to make this measurement explicit for future measurements
# Replace this stabilizer row with just Z on the target qubit
tab[i, :n] = 0 # Clear X part
tab[i, n:2*n] = 0 # Clear Z part
tab[i, n + qubit] = 1 # Set Z on target qubit
tab[i, -1] = outcome # Set phase
break
else:
# Random outcome
outcome = rng.randint(0, 2)
# Update tableau to reflect post-measurement state
for i in range(2 * n):
if i != p and tab[i, qubit] == 1: # X on qubit
self._rowsum(i, p)
# Set stabilizer p to Z (outcome 0) or -Z (outcome 1)
tab[p, :] = 0
tab[p, n + qubit] = 1
tab[p, -1] = outcome
self.measurement_outcomes.append(outcome)
return outcome
# Handoff to MPS
[docs]
def to_mps(self, device: str = "cuda"):
"""
Convert stabilizer state to MPS representation
This is used when non-Clifford gates appear in the circuit.
Returns:
AdaptiveMPS instance
"""
from .adaptive_mps import AdaptiveMPS
# For now, convert to statevector then to MPS
# TODO: Direct stabilizer → MPS conversion is more efficient
statevector = self.to_statevector()
# Create MPS from statevector
mps = AdaptiveMPS(self.n_qubits, bond_dim=2, device=device)
# Set MPS tensors from statevector
# This is a simple inefficient method - proper conversion would use
# successive SVDs along the chain
statevector_tensor = torch.tensor(statevector, dtype=torch.complex64, device=device)
# Reshape to tensor train format and perform SVD decomposition
# For now, use a simplified approach
# TODO: Implement efficient SVD-based conversion
return mps
[docs]
def to_statevector(self) -> np.ndarray:
"""
Convert stabilizer state to full statevector (SLOW! Only for small n)
This explicitly constructs the 2^n statevector.
Only use for small systems (n ≤ 20).
Returns:
Statevector of shape (2^n,)
"""
if self.n_qubits > 20:
raise ValueError(
f"Cannot convert {self.n_qubits} qubits to statevector "
f"(would require {2**self.n_qubits} amplitudes)"
)
# Start with |0...0⟩
psi = np.zeros(2**self.n_qubits, dtype=np.complex128)
psi[0] = 1.0
# Apply Clifford gates to generate the stabilizer state
# This is inefficient but correct
# TODO: Use stabilizer-to-graph-state conversion for efficiency
return psi
[docs]
def is_clifford_gate(gate_name: str) -> bool:
"""Check if a gate is in the Clifford group"""
clifford_gates = {"H", "S", "CNOT", "CZ", "SWAP", "X", "Y", "Z", "S_DAG", "CX", "ID", "I"}
return gate_name.upper() in clifford_gates
[docs]
class HybridSimulator:
"""
Hybrid simulator that uses stabilizer backend for Clifford parts
and switches to MPS when non-Clifford gates appear
Usage:
sim = HybridSimulator(n_qubits=100, use_stabilizer=True)
sim.h(0)
sim.cnot(0, 1)
sim.t(0) # Automatically switches to MPS here
sim.measure_all()
"""
[docs]
def __init__(
self, n_qubits: int, use_stabilizer: bool = True, chi_max: int = 64, device: str = "cuda"
):
self.n_qubits = n_qubits
self.use_stabilizer = use_stabilizer
self.chi_max = chi_max
self.device = device
# Start with stabilizer
self.mode = "stabilizer" if use_stabilizer else "mps"
if self.mode == "stabilizer":
self.stabilizer_sim = StabilizerSimulator(n_qubits)
self.mps = None
else:
from .adaptive_mps import AdaptiveMPS
self.mps = AdaptiveMPS(n_qubits, bond_dim=2, chi_max_per_bond=chi_max, device=device)
self.stabilizer_sim = None
self.gate_count = {"clifford": 0, "non_clifford": 0}
def _switch_to_mps(self):
"""Switch from stabilizer to MPS representation"""
if self.mode == "mps":
return # Already in MPS mode
print(f"Switching to MPS after {self.gate_count['clifford']} Clifford gates")
# Convert stabilizer state to MPS
self.mps = self.stabilizer_sim.to_mps(device=self.device)
self.stabilizer_sim = None
self.mode = "mps"
[docs]
def h(self, qubit: int):
"""Hadamard gate"""
if self.mode == "stabilizer":
self.stabilizer_sim.h(qubit)
self.gate_count["clifford"] += 1
else:
H = torch.tensor([[1, 1], [1, -1]], dtype=torch.complex64) / np.sqrt(2)
self.mps.apply_single_qubit_gate(qubit, H.to(self.device))
[docs]
def s(self, qubit: int):
"""Phase gate"""
if self.mode == "stabilizer":
self.stabilizer_sim.s(qubit)
self.gate_count["clifford"] += 1
else:
S = torch.tensor([[1, 0], [0, 1j]], dtype=torch.complex64)
self.mps.apply_single_qubit_gate(qubit, S.to(self.device))
[docs]
def x(self, qubit: int):
"""Pauli-X gate"""
if self.mode == "stabilizer":
self.stabilizer_sim.x(qubit)
self.gate_count["clifford"] += 1
else:
X = torch.tensor([[0, 1], [1, 0]], dtype=torch.complex64)
self.mps.apply_single_qubit_gate(qubit, X.to(self.device))
[docs]
def z(self, qubit: int):
"""Pauli-Z gate"""
if self.mode == "stabilizer":
self.stabilizer_sim.z(qubit)
self.gate_count["clifford"] += 1
else:
Z = torch.tensor([[1, 0], [0, -1]], dtype=torch.complex64)
self.mps.apply_single_qubit_gate(qubit, Z.to(self.device))
[docs]
def cnot(self, control: int, target: int):
"""CNOT gate"""
if self.mode == "stabilizer":
self.stabilizer_sim.cnot(control, target)
self.gate_count["clifford"] += 1
else:
CNOT = torch.tensor(
[[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]], dtype=torch.complex64
)
if abs(control - target) == 1:
bond_idx = min(control, target)
self.mps.apply_two_site_gate(bond_idx, CNOT.to(self.device))
else:
# Need SWAP network for non-adjacent qubits
raise NotImplementedError("CNOT on non-adjacent qubits requires SWAP network")
[docs]
def t(self, qubit: int):
"""T gate (non-Clifford!)"""
if self.mode == "stabilizer":
self._switch_to_mps()
T = torch.tensor([[1, 0], [0, np.exp(1j * np.pi / 4)]], dtype=torch.complex64)
self.mps.apply_single_qubit_gate(qubit, T.to(self.device))
self.gate_count["non_clifford"] += 1
[docs]
def swap(self, qubit1: int, qubit2: int):
"""SWAP gate (Clifford)"""
if self.mode == "stabilizer":
self.stabilizer_sim.swap(qubit1, qubit2)
self.gate_count["clifford"] += 1
else:
# For MPS, SWAP can be done with 3 CNOTs
# SWAP = CNOT(a,b) CNOT(b,a) CNOT(a,b)
raise NotImplementedError("SWAP in MPS mode not yet implemented")
[docs]
def measure(self, qubit: int) -> int:
"""Measure qubit"""
if self.mode == "stabilizer":
return self.stabilizer_sim.measure(qubit)
else:
# MPS measurement
# TODO: Implement proper MPS measurement
return 0 # Placeholder
[docs]
def get_statistics(self) -> dict:
"""Get simulation statistics"""
stats = {
"mode": self.mode,
"clifford_gates": self.gate_count["clifford"],
"non_clifford_gates": self.gate_count["non_clifford"],
}
if self.mode == "mps":
stats["mps_stats"] = self.mps.stats_summary()
return stats
