Simulation in Cirq

This guide covers quantum circuit simulation, including exact and noisy simulations, parameter sweeps, and the Quantum Virtual Machine (QVM).

Exact Simulation

Basic Simulation

# Create circuit
q0, q1 = cirq.LineQubit.range(2)
circuit = cirq.Circuit(
    cirq.H(q0),
    cirq.CNOT(q0, q1),
    cirq.measure(q0, q1, key='result')
)

# Simulate
simulator = cirq.Simulator()
result = simulator.run(circuit, repetitions=1000)

# Get measurement results
print(result.histogram(key='result'))

State Vector Simulation

# Simulate without measurement to get final state
simulator = cirq.Simulator()
result = simulator.simulate(circuit_without_measurement)

# Access state vector
state_vector = result.final_state_vector
print(f"State vector: {state_vector}")

# Get amplitudes
print(f"Amplitude of |00⟩: {state_vector[0]}")
print(f"Amplitude of |11⟩: {state_vector[3]}")

Density Matrix Simulation

# Use density matrix simulator for mixed states
simulator = cirq.DensityMatrixSimulator()
result = simulator.simulate(circuit)

# Access density matrix
density_matrix = result.final_density_matrix
print(f"Density matrix shape: {density_matrix.shape}")

Step-by-Step Simulation

# Simulate moment-by-moment
simulator = cirq.Simulator()
for step in simulator.simulate_moment_steps(circuit):
    print(f"State after moment {step.moment}: {step.state_vector()}")

Sampling and Measurements

Run Multiple Shots

# Run circuit multiple times
result = simulator.run(circuit, repetitions=10000)

# Access measurement counts
counts = result.histogram(key='result')
print(f"Measurement counts: {counts}")

# Get raw measurements
measurements = result.measurements['result']
print(f"Shape: {measurements.shape}")  # (repetitions, num_qubits)

Expectation Values

# Measure observable expectation value
from cirq import PauliString

observable = PauliString({q0: cirq.Z, q1: cirq.Z})
result = simulator.simulate_expectation_values(
    circuit,
    observables=[observable]
)
print(f"⟨ZZ⟩ = {result[0]}")

Parameter Sweeps

Sweep Over Parameters

# Create parameterized circuit
theta = sympy.Symbol('theta')
q = cirq.LineQubit(0)
circuit = cirq.Circuit(
    cirq.ry(theta)(q),
    cirq.measure(q, key='m')
)

# Define parameter sweep
sweep = cirq.Linspace(key='theta', start=0, stop=2*np.pi, length=50)

# Run sweep
simulator = cirq.Simulator()
results = simulator.run_sweep(circuit, params=sweep, repetitions=1000)

# Process results
for params, result in zip(sweep, results):
    theta_val = params['theta']
    counts = result.histogram(key='m')
    print(f"θ={theta_val:.2f}: {counts}")

Multiple Parameters

# Sweep over multiple parameters
theta = sympy.Symbol('theta')
phi = sympy.Symbol('phi')

circuit = cirq.Circuit(
    cirq.ry(theta)(q0),
    cirq.rz(phi)(q1)
)

# Product sweep (all combinations)
sweep = cirq.Product(
    cirq.Linspace('theta', 0, np.pi, 10),
    cirq.Linspace('phi', 0, 2*np.pi, 10)
)

results = simulator.run_sweep(circuit, params=sweep, repetitions=100)

Zip Sweep (Paired Parameters)

# Sweep parameters together
sweep = cirq.Zip(
    cirq.Linspace('theta', 0, np.pi, 20),
    cirq.Linspace('phi', 0, 2*np.pi, 20)
)

results = simulator.run_sweep(circuit, params=sweep, repetitions=100)

Noisy Simulation

Adding Noise Channels

# Create noisy circuit
noisy_circuit = circuit.with_noise(cirq.depolarize(p=0.01))

# Simulate noisy circuit
simulator = cirq.DensityMatrixSimulator()
result = simulator.run(noisy_circuit, repetitions=1000)

Custom Noise Models

# Apply different noise to different gates
noise_model = cirq.NoiseModel.from_noise_model_like(
    cirq.ConstantQubitNoiseModel(cirq.depolarize(0.01))
)

# Simulate with noise model
result = cirq.DensityMatrixSimulator(noise=noise_model).run(
    circuit, repetitions=1000
)

See noise.md for comprehensive noise modeling details.

State Histograms

Visualize Results

# Get histogram
result = simulator.run(circuit, repetitions=1000)
counts = result.histogram(key='result')

# Plot
plt.bar(counts.keys(), counts.values())
plt.xlabel('State')
plt.ylabel('Counts')
plt.title('Measurement Results')
plt.show()

State Probability Distribution

# Get state vector
result = simulator.simulate(circuit_without_measurement)
state_vector = result.final_state_vector

# Compute probabilities
probabilities = np.abs(state_vector) ** 2

# Plot
plt.bar(range(len(probabilities)), probabilities)
plt.xlabel('Basis State Index')
plt.ylabel('Probability')
plt.show()

Quantum Virtual Machine (QVM)

QVM simulates realistic quantum hardware with device-specific constraints and noise.

Using Virtual Devices

# Use a virtual Google device

# Get virtual device
device = cirq_google.Sycamore

# Create circuit on device
qubits = device.metadata.qubit_set
circuit = cirq.Circuit(device=device)

# Add operations respecting device constraints
circuit.append(cirq.CZ(qubits[0], qubits[1]))

# Validate circuit against device
device.validate_circuit(circuit)

Noisy Virtual Hardware

# Simulate with device noise
processor = cirq_google.get_engine().get_processor('weber')
noise_props = processor.get_device_specification()

# Create realistic noisy simulator
noisy_sim = cirq.DensityMatrixSimulator(
    noise=cirq_google.NoiseModelFromGoogleNoiseProperties(noise_props)
)

result = noisy_sim.run(circuit, repetitions=1000)

Advanced Simulation Techniques

Custom Initial State

# Start from custom state
initial_state = np.array([1, 0, 0, 1]) / np.sqrt(2)  # |00⟩ + |11⟩

simulator = cirq.Simulator()
result = simulator.simulate(circuit, initial_state=initial_state)

Partial Trace

# Trace out subsystems
result = simulator.simulate(circuit)
full_state = result.final_state_vector

# Compute reduced density matrix for first qubit
from cirq import partial_trace
reduced_dm = partial_trace(result.final_density_matrix, keep_indices=[0])

Intermediate State Access

# Get state at specific moment
simulator = cirq.Simulator()
for i, step in enumerate(simulator.simulate_moment_steps(circuit)):
    if i == 5:  # After 5th moment
        state = step.state_vector()
        print(f"State after moment 5: {state}")
        break

Simulation Performance

Optimizing Large Simulations

1. Use state vector for pure states: Faster than density matrix
2. Avoid density matrix when possible: Exponentially more expensive
3. Batch parameter sweeps: More efficient than individual runs
4. Use appropriate repetitions: Balance accuracy vs computation time

# Efficient: Single sweep
results = simulator.run_sweep(circuit, params=sweep, repetitions=100)

# Inefficient: Multiple individual runs
results = [simulator.run(circuit, param_resolver=p, repetitions=100)
           for p in sweep]

Memory Considerations

# For large systems, monitor state vector size
n_qubits = 20
state_size = 2**n_qubits * 16  # bytes (complex128)
print(f"State vector size: {state_size / 1e9:.2f} GB")

Stabilizer Simulation

For circuits with only Clifford gates, use efficient stabilizer simulation:

# Clifford circuit (H, S, CNOT)
circuit = cirq.Circuit(
    cirq.H(q0),
    cirq.S(q1),
    cirq.CNOT(q0, q1)
)

# Use stabilizer simulator (exponentially faster)
simulator = cirq.CliffordSimulator()
result = simulator.run(circuit, repetitions=1000)

Best Practices

1. Choose appropriate simulator: Use Simulator for pure states, DensityMatrixSimulator for mixed states
2. Use parameter sweeps: More efficient than running individual circuits
3. Validate circuits: Check circuit validity before long simulations
4. Monitor resource usage: Track memory for large-scale simulations
5. Use stabilizer simulation: When circuits contain only Clifford gates
6. Save intermediate results: For long parameter sweeps or optimization runs