Circuit superoperator (#4550)

cirq-core/cirq/circuits/circuit.py

@@ -22,9 +22,9 @@
 import abc
 import enum
 import html
+import itertools
 import math
 from collections import defaultdict
-from itertools import groupby
 from typing import (
     AbstractSet,
     Any,
@@ -999,6 +999,20 @@ def unitary(
         result = _apply_unitary_circuit(self, state, qs, dtype)
         return result.reshape((side_len, side_len))
 
+    def _superoperator_(self) -> np.ndarray:
+        """Compute superoperator matrix for quantum channel specified by this circuit."""
+        all_qubits = self.all_qubits()
+        n = len(all_qubits)
+        if n > 10:
+            raise ValueError(f"{n} > 10 qubits is too many to compute superoperator")
+
+        circuit_superoperator = np.eye(4 ** n)
+        for moment in self:
+            full_moment = moment.expand_to(all_qubits)
+            moment_superoperator = full_moment._superoperator_()
+            circuit_superoperator = moment_superoperator @ circuit_superoperator
+        return circuit_superoperator
+
     def final_state_vector(
         self,
         initial_state: 'cirq.STATE_VECTOR_LIKE' = 0,
@@ -2638,4 +2652,4 @@ def _group_until_different(items: Iterable[TIn], key: Callable[[TIn], TKey], val
     Yields:
         Tuples containing the group key and item values.
     """
-    return ((k, [val(i) for i in v]) for (k, v) in groupby(items, key))
+    return ((k, [val(i) for i in v]) for (k, v) in itertools.groupby(items, key))

cirq-core/cirq/circuits/circuit_test.py

@@ -11,10 +11,11 @@
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 # See the License for the specific language governing permissions and
 # limitations under the License.
+import itertools
 import os
 from collections import defaultdict
 from random import randint, random, sample, randrange
-from typing import Optional, Tuple, TYPE_CHECKING
+from typing import Iterator, Optional, Tuple, TYPE_CHECKING
 
 import numpy as np
 import pytest
@@ -69,6 +70,8 @@ def can_add_operation_into_moment(
 if TYPE_CHECKING:
     import cirq
 
+q0, q1, q2, q3 = cirq.LineQubit.range(4)
+
 
 class _MomentAndOpTypeValidatingDeviceType(cirq.Device):
     def validate_operation(self, operation):
@@ -2851,6 +2854,159 @@ def _decompose_(self, qubits):
     cirq.testing.assert_allclose_up_to_global_phase(mat, mat_expected, atol=1e-8)
 
 
+def test_circuit_superoperator_too_many_qubits():
+    circuit = cirq.Circuit(cirq.IdentityGate(num_qubits=11).on(*cirq.LineQubit.range(11)))
+    with pytest.raises(ValueError, match="too many"):
+        _ = circuit._superoperator_()
+
+
+@pytest.mark.parametrize(
+    'circuit, expected_superoperator',
+    (
+        (cirq.Circuit(cirq.I(q0)), np.eye(4)),
+        (cirq.Circuit(cirq.IdentityGate(2).on(q0, q1)), np.eye(16)),
+        (
+            cirq.Circuit(cirq.H(q0)),
+            np.array([[1, 1, 1, 1], [1, -1, 1, -1], [1, 1, -1, -1], [1, -1, -1, 1]]) / 2,
+        ),
+        (cirq.Circuit(cirq.S(q0)), np.diag([1, -1j, 1j, 1])),
+        (cirq.Circuit(cirq.depolarize(0.75).on(q0)), np.outer([1, 0, 0, 1], [1, 0, 0, 1]) / 2),
+        (
+            cirq.Circuit(cirq.X(q0), cirq.depolarize(0.75).on(q0)),
+            np.outer([1, 0, 0, 1], [1, 0, 0, 1]) / 2,
+        ),
+        (
+            cirq.Circuit(cirq.Y(q0), cirq.depolarize(0.75).on(q0)),
+            np.outer([1, 0, 0, 1], [1, 0, 0, 1]) / 2,
+        ),
+        (
+            cirq.Circuit(cirq.Z(q0), cirq.depolarize(0.75).on(q0)),
+            np.outer([1, 0, 0, 1], [1, 0, 0, 1]) / 2,
+        ),
+        (
+            cirq.Circuit(cirq.H(q0), cirq.depolarize(0.75).on(q0)),
+            np.outer([1, 0, 0, 1], [1, 0, 0, 1]) / 2,
+        ),
+        (cirq.Circuit(cirq.H(q0), cirq.H(q0)), np.eye(4)),
+        (
+            cirq.Circuit(cirq.H(q0), cirq.CNOT(q1, q0), cirq.H(q0)),
+            np.diag([1, 1, 1, -1, 1, 1, 1, -1, 1, 1, 1, -1, -1, -1, -1, 1]),
+        ),
+    ),
+)
+def test_circuit_superoperator_fixed_values(circuit, expected_superoperator):
+    """Tests Circuit._superoperator_() on a few simple circuits."""
+    assert np.allclose(circuit._superoperator_(), expected_superoperator)
+
+
+@pytest.mark.parametrize(
+    'rs, n_qubits',
+    (
+        ([0.1, 0.2], 1),
+        ([0.1, 0.2], 2),
+        ([0.8, 0.9], 1),
+        ([0.8, 0.9], 2),
+        ([0.1, 0.2, 0.3], 1),
+        ([0.1, 0.2, 0.3], 2),
+        ([0.1, 0.2, 0.3], 3),
+    ),
+)
+def test_circuit_superoperator_depolarizing_channel_compositions(rs, n_qubits):
+    """Tests Circuit._superoperator_() on compositions of depolarizing channels."""
+
+    def pauli_error_probability(r: float, n_qubits: int) -> float:
+        """Computes Pauli error probability for given depolarization parameter.
+
+        Pauli error is what cirq.depolarize takes as argument. Depolarization parameter
+        makes it simple to compute the serial composition of depolarizing channels. It
+        is multiplicative under channel composition.
+        """
+        d2 = 4 ** n_qubits
+        return (1 - r) * (d2 - 1) / d2
+
+    def depolarize(r: float, n_qubits: int) -> cirq.DepolarizingChannel:
+        """Returns depolarization channel with given depolarization parameter."""
+        return cirq.depolarize(pauli_error_probability(r, n_qubits=n_qubits), n_qubits=n_qubits)
+
+    qubits = cirq.LineQubit.range(n_qubits)
+    circuit1 = cirq.Circuit(depolarize(r, n_qubits).on(*qubits) for r in rs)
+    circuit2 = cirq.Circuit(depolarize(np.prod(rs), n_qubits).on(*qubits))
+
+    cm1 = circuit1._superoperator_()
+    cm2 = circuit2._superoperator_()
+    assert np.allclose(cm1, cm2)
+
+
+def density_operator_basis(n_qubits: int) -> Iterator[np.ndarray]:
+    """Yields operator basis consisting of density operators."""
+    RHO_0 = np.array([[1, 0], [0, 0]], dtype=np.complex64)
+    RHO_1 = np.array([[0, 0], [0, 1]], dtype=np.complex64)
+    RHO_2 = np.array([[1, 1], [1, 1]], dtype=np.complex64) / 2
+    RHO_3 = np.array([[1, -1j], [1j, 1]], dtype=np.complex64) / 2
+    RHO_BASIS = (RHO_0, RHO_1, RHO_2, RHO_3)
+
+    if n_qubits < 1:
+        yield np.array(1)
+        return
+    for rho1 in RHO_BASIS:
+        for rho2 in density_operator_basis(n_qubits - 1):
+            yield np.kron(rho1, rho2)
+
+
+@pytest.mark.parametrize(
+    'circuit, initial_state',
+    itertools.chain(
+        itertools.product(
+            [
+                cirq.Circuit(cirq.I(q0)),
+                cirq.Circuit(cirq.X(q0)),
+                cirq.Circuit(cirq.Y(q0)),
+                cirq.Circuit(cirq.Z(q0)),
+                cirq.Circuit(cirq.S(q0)),
+                cirq.Circuit(cirq.T(q0)),
+            ],
+            density_operator_basis(n_qubits=1),
+        ),
+        itertools.product(
+            [
+                cirq.Circuit(cirq.H(q0), cirq.CNOT(q0, q1)),
+                cirq.Circuit(cirq.depolarize(0.2).on(q0), cirq.CNOT(q0, q1)),
+                cirq.Circuit(
+                    cirq.X(q0),
+                    cirq.amplitude_damp(0.2).on(q0),
+                    cirq.depolarize(0.1).on(q1),
+                    cirq.CNOT(q0, q1),
+                ),
+            ],
+            density_operator_basis(n_qubits=2),
+        ),
+        itertools.product(
+            [
+                cirq.Circuit(
+                    cirq.depolarize(0.1, n_qubits=2).on(q0, q1),
+                    cirq.H(q2),
+                    cirq.CNOT(q1, q2),
+                    cirq.phase_damp(0.1).on(q0),
+                ),
+                cirq.Circuit(cirq.H(q0), cirq.H(q1), cirq.TOFFOLI(q0, q1, q2)),
+            ],
+            density_operator_basis(n_qubits=3),
+        ),
+    ),
+)
+def test_compare_circuits_superoperator_to_simulation(circuit, initial_state):
+    """Compares action of circuit superoperator and circuit simulation."""
+    superoperator = circuit._superoperator_()
+    vectorized_initial_state = np.reshape(initial_state, np.prod(initial_state.shape))
+    vectorized_final_state = superoperator @ vectorized_initial_state
+    actual_state = np.reshape(vectorized_final_state, initial_state.shape)
+
+    sim = cirq.DensityMatrixSimulator()
+    expected_state = sim.simulate(circuit, initial_state=initial_state).final_density_matrix
+
+    assert np.allclose(actual_state, expected_state)
+
+
 @pytest.mark.parametrize('circuit_cls', [cirq.Circuit, cirq.FrozenCircuit])
 def test_expanding_gate_symbols(circuit_cls):
     class MultiTargetCZ(cirq.Gate):

cirq-core/cirq/ops/moment.py

@@ -24,13 +24,18 @@
     Iterator,
     overload,
     Optional,
+    Sequence,
     Tuple,
     TYPE_CHECKING,
     TypeVar,
     Union,
 )
 
-from cirq import protocols, ops, value
+import itertools
+
+import numpy as np
+
+from cirq import protocols, ops, qis, value
 from cirq._import import LazyLoader
 from cirq.ops import raw_types
 from cirq.protocols import circuit_diagram_info_protocol
@@ -327,6 +332,80 @@ def transform_qubits(
         """
         return self.__class__(op.transform_qubits(qubit_map) for op in self.operations)
 
+    def expand_to(self, qubits: Iterable['cirq.Qid']) -> 'cirq.Moment':
+        """Returns self expanded to given superset of qubits by making identities explicit."""
+        if not self.qubits.issubset(qubits):
+            raise ValueError(f'{qubits} is not a superset of {self.qubits}')
+
+        operations = list(self.operations)
+        for q in set(qubits) - self.qubits:
+            operations.append(ops.I(q))
+        return Moment(*operations)
+
+    def _kraus_(self) -> Sequence[np.ndarray]:
+        r"""Returns Kraus representation of self.
+
+        The method computes a Kraus representation of self from Kraus representations of its
+        constituent operations by taking the tensor product of Kraus operators along all paths
+        corresponding to the possible choices of the Kraus operators of the operations. More
+        precisely, it computes all terms in the expression
+
+        $$
+        \sum_{i_1} \sum_{i_2} \dots \sum_{i_m} \bigotimes_{k=1}^m K_{k,i_k}
+        $$
+
+        where $K_{k,j}$ is the jth Kraus operator of the kth operation in self. Each term becomes
+        an element in the sequence returned by this method.
+
+        Args:
+            self: This Moment.
+        Returns:
+            A Kraus representation of self.
+        Raises:
+            ValueError: If self uses more than ten qubits as the length of the resulting sequence
+            is the product of the lengths of the Kraus representations returned by _kraus_ for
+            each constituent operation.
+        """
+        qubits = sorted(self.qubits)
+        n = len(qubits)
+        if n < 1:
+            return (np.array([[1 + 0j]]),)
+        if n > 10:
+            raise ValueError(f'Cannot compute Kraus representation of moment with {n} > 10 qubits')
+
+        qubit_to_row_subscript = dict(zip(qubits, 'abcdefghij'))
+        qubit_to_col_subscript = dict(zip(qubits, 'ABCDEFGHIJ'))
+
+        def row_subscripts(qs: Sequence['cirq.Qid']) -> str:
+            return ''.join(qubit_to_row_subscript[q] for q in qs)
+
+        def col_subscripts(qs: Sequence['cirq.Qid']) -> str:
+            return ''.join(qubit_to_col_subscript[q] for q in qs)
+
+        def kraus_tensors(op: 'cirq.Operation') -> Sequence[np.ndarray]:
+            return tuple(np.reshape(k, (2, 2) * len(op.qubits)) for k in protocols.kraus(op))
+
+        input_subscripts = ','.join(
+            row_subscripts(op.qubits) + col_subscripts(op.qubits) for op in self.operations
+        )
+        output_subscripts = row_subscripts(qubits) + col_subscripts(qubits)
+        assert len(input_subscripts) == 2 * n + len(self.operations) - 1
+        assert len(output_subscripts) == 2 * n
+
+        transpose = input_subscripts + '->' + output_subscripts
+
+        r = []
+        d = 2 ** n
+        kss = [kraus_tensors(op) for op in self.operations]
+        for ks in itertools.product(*kss):
+            k = np.einsum(transpose, *ks)
+            r.append(np.reshape(k, (d, d)))
+        return r
+
+    def _superoperator_(self) -> np.ndarray:
+        """Returns superoperator representation of self."""
+        return qis.kraus_to_superoperator(self._kraus_())
+
     def _json_dict_(self) -> Dict[str, Any]:
         return protocols.obj_to_dict_helper(self, ['operations'])