Noise simulation in qsimcirq

Noisy gates in Cirq are represented by Channels, which can act as one of a set of gates depending on the state of the circuit. The Cirq tutorial on noise explains how to construct these objects and add them to your circuits.

Setup

Install the Cirq and qsimcirq packages:

try:
    import cirq
except ImportError:
    !pip install cirq --quiet
    import cirq

try:
    import qsimcirq
except ImportError:
    !pip install qsimcirq --quiet
    import qsimcirq

It is possible to simulate channels with density matrices, which combine all possible channel behaviors, but the overhead is steep: a density matrix requires O(4^N) storage for N qubits.

In qsimcirq, noisy circuits are instead simulated as "trajectories": the behavior of each Channel is determined probabilistically at runtime. This permits much larger simulations at the cost of only capturing one such "trajectory" per execution.

Performance

Noisy circuits tend to be more expensive to simulate than their noiseless equivalents, but qsim is optimized to avoid these overheads when possible. In particular, the less incoherent noise (i.e. non-unitary effects) that a Channel has, the closer its performance will be to the noiseless case for a single repetition.

Simulating many repetitions of a noisy circuit requires executing the entire circuit once for each repetition due to the nondeterministic nature of noisy operations.

Constructing noisy circuits

Cirq provides a number of tools for constructing noisy circuits. For the purpose of this tutorial, we will focus on two common types of noise: T1 ("amplitude damping") and T2 ("phase damping"). These can be created in Cirq with cirq.amplitude_damp and cirq.phase_damp, as shown below:

q0, q1 = cirq.LineQubit.range(2)

circuit = cirq.Circuit(
    # Perform a Hadamard on both qubits
    cirq.H(q0), cirq.H(q1),
    # Apply amplitude damping to q0 with probability 0.1
    cirq.amplitude_damp(gamma=0.1).on(q0),
    # Apply phase damping to q1 with probability 0.1
    cirq.phase_damp(gamma=0.1).on(q1),
)

Simulating noisy circuits

Simulating this circuit works exactly the same as simulating a noiseless circuit: simply construct a simulator object and simulate. QSimSimulator will automatically switch over to the noisy simulator if it detects noise (i.e. Channels) in your circuit.

qsim_simulator = qsimcirq.QSimSimulator()
results = qsim_simulator.simulate(circuit)
print(results.final_state_vector)

[0.52631575+0.j 0.49930704+0.j 0.49930704+0.j 0.47368425+0.j]

It's important to note that unlike density-matrix simulations, this result (from a single repetition) is stochastic in nature. Running the circuit multiple times may yield different results, but each result generated is a possible outcome of the provided circuit.

Other simulation modes

Noisy circuit simulation in qsimcirq supports all of the same simulation modes as the noiseless simulator, including:

Measurement Sampling

# Simulate measuring at the end of the circuit.
measured_circuit = circuit + cirq.measure(q0, q1, key='m')
measure_results = qsim_simulator.run(measured_circuit, repetitions=5)
print(measure_results)

m=01101, 00100

Amplitude evaluation

# Calculate only the amplitudes of the |00) and |01) states.
amp_results = qsim_simulator.compute_amplitudes(
    circuit, bitstrings=[0b00, 0b01])
print(amp_results)

[(0.5263157486915588+0j), (0.4993070363998413+0j)]

Expectation values

Expectation values can only be estimated from trajectories, but the accuracy of these estimates can be increased by simulating the circuit additional times. This is demonstrated below.

# Set the "noisy repetitions" to 100.
# This parameter only affects expectation value calculations.
options = {'r': 100}
# Also set the random seed to get reproducible results.
ev_simulator = qsimcirq.QSimSimulator(qsim_options=options, seed=1)
# Define observables to measure: <Z> for q0 and <X> for q1.
pauli_sum1 = cirq.Z(q0)
pauli_sum2 = cirq.X(q1)
# Calculate expectation values for the given observables.
ev_results = ev_simulator.simulate_expectation_values(
    circuit,
    observables=[pauli_sum1, pauli_sum2],
)
print(ev_results)

[(0.13789467364549637+0j), (0.9386972016096116+0j)]

The output is a list of expectation values, one for each observable.