Cirq: Getting The Instances (samples) of a Quantum Circuit with Probabilistic Unitaries or Mixtures

Question

The wavefunction simulator in Cirq uses a Monte Carlo approach to simulate a certain subset of quantum noise channels, namely through probabilistic/stochastic application of unitary gates. These are called MixedUnitaryGates, mixtures, etc. in Cirq.

Instead of creating the Kraus operators of a channel to deterministically apply to the NxN dimensional density matrix, the simulator here samples a unitary gate for given probabilities, gets a noisy circuit with all unitary operations, and apply this circuit/unitary on the initial wavefunction to acquire the final wavefunction. Given a noisy circuit, the cirq.Simulator() does exactly this way to output a final state vector. Here, each instance gives a different state vector.

I could not exactly verify Cirq's pipeline. But I need an array/list of these sampled instances of noisy circuits. Do you know if Cirq has a function to sample unitary circuits out of noisy circuits with mixtures (i.e, probabilistic unitaries)?

Let's say I want to have M=100 circuit samples from a circuit with depolarizing channel with p=0.03:

circuit=cirq.Circuit()
noisy_circuit=circuit.with_noise(cirq.depolarize(0.1))

#After calling the mystery function we sample M=100 instances of unitary circuits

sampled_circuits=MysteryFunction(noisy_circuit,M=100)

Is there a MysteryFunction of the above form? Is there a way to sample these unitary circuits out of circuits with channels that support Mixture format?

Thank you

Answer 1

There's no such function built in, but you can use map_operations to iterate through the circuit and replace mixture ops with a selected unitary gate as follows.

def sample_mixture_as_unitary_circuit(circuit: cirq.Circuit, prng: np.random.RandomState):
    def sample_op(op: cirq.Operation, _: int):
        if cirq.has_mixture(op) and not cirq.has_unitary(op):
            mixture = cirq.mixture(op)
            probabilities, unitaries = zip(*mixture)
            index = prng.choice(len(unitaries), p=probabilities)
            op = cirq.MatrixGate(unitaries[index]).on(*op.qubits)
        return op

    return cirq.map_operations(circuit, sample_op)

To use that, just apply to the noisy circuit in a loop:

q = cirq.LineQubit(0)
circuit = cirq.Circuit(cirq.X(q))
noisy_circuit = circuit.with_noise(cirq.depolarize(0.3))
prng = np.random.RandomState()
print()
print('original')
print(noisy_circuit)
print()
print('samples')
circuits = [sample_mixture_as_unitary_circuit(noisy_circuit, prng) for _ in range(10)]
for c in circuits:
    print(c)

Example output:

original
0: ───X───D(0.3)[cirq.VirtualTag()]───

samples
          ┌             ┐
0: ───X───│0.+0.j 0.-1.j│───
          │0.+1.j 0.+0.j│
          └             ┘
          ┌     ┐
0: ───X───│1. 0.│───
          │0. 1.│
          └     ┘
          ┌     ┐
0: ───X───│1. 0.│───
          │0. 1.│
          └     ┘
          ┌     ┐
0: ───X───│1. 0.│───
          │0. 1.│
          └     ┘
          ┌     ┐
0: ───X───│1. 0.│───
          │0. 1.│
          └     ┘
          ┌     ┐
0: ───X───│1. 0.│───
          │0. 1.│
          └     ┘
          ┌     ┐
0: ───X───│1. 0.│───
          │0. 1.│
          └     ┘
          ┌             ┐
0: ───X───│0.+0.j 1.+0.j│───
          │1.+0.j 0.+0.j│
          └             ┘
          ┌     ┐
0: ───X───│1. 0.│───
          │0. 1.│
          └     ┘
          ┌     ┐
0: ───X───│1. 0.│───
          │0. 1.│
          └     ┘