Implement cutting of general 2-qubit unitaries

circuit_knitting/cutting/qpd/qpd.py

@@ -32,9 +32,10 @@
     ZGate,
     HGate,
     SGate,
-    TGate,
     SdgGate,
     SXGate,
+    SXdgGate,
+    TGate,
     RXGate,
     RYGate,
     RZGate,
@@ -50,7 +51,12 @@
     CRZGate,
     ECRGate,
     CSXGate,
+    SwapGate,
+    iSwapGate,
+    DCXGate,
 )
+from qiskit.extensions import UnitaryGate
+from qiskit.quantum_info.synthesis.two_qubit_decompose import TwoQubitWeylDecomposition
 
 from .qpd_basis import QPDBasis
 from .instructions import BaseQPDGate, TwoQubitQPDGate, QPDMeasure
@@ -221,10 +227,25 @@ def qpdbasis_from_gate(gate: Gate) -> QPDBasis:
     try:
         f = _qpdbasis_from_gate_funcs[gate.name]
     except KeyError:
-        raise ValueError(f"Gate not supported: {gate.name}") from None
+        pass
     else:
         return f(gate)
 
+    if isinstance(gate, Gate) and gate.num_qubits == 2:
+        mat = gate.to_matrix()
+        d = TwoQubitWeylDecomposition(mat)
+        u = _u_from_thetavec([d.a, d.b, d.c])
+        retval = _nonlocal_qpd_basis_from_u(u)
+        for operations in unique_by_id(m[0] for m in retval.maps):
+            operations.insert(0, UnitaryGate(d.K2r))
+            operations.append(UnitaryGate(d.K1r))
+        for operations in unique_by_id(m[1] for m in retval.maps):
+            operations.insert(0, UnitaryGate(d.K2l))
+            operations.append(UnitaryGate(d.K1l))
+        return retval
+
+    raise ValueError(f"Gate not supported: {gate.name}") from None
+
 
 def supported_gates() -> set[str]:
     """
@@ -236,6 +257,176 @@ def supported_gates() -> set[str]:
     return set(_qpdbasis_from_gate_funcs)
 
 
+def _copy_unique_sublists(lsts: tuple[list, ...], /) -> tuple[list, ...]:
+    """
+    Copy each list in a sequence of lists while preserving uniqueness.
+
+    This is useful to ensure that the two sets of ``maps`` in a
+    :class:`QPDBasis` will be independent of each other.  This enables one to
+    subsequently edit the ``maps`` independently of each other (e.g., to apply
+    single-qubit pre- or post-rotations.
+    """
+    copy_by_id: dict[int, list] = {}
+    for lst in lsts:
+        if id(lst) not in copy_by_id:
+            copy_by_id[id(lst)] = lst.copy()
+    return tuple(copy_by_id[id(lst)] for lst in lsts)
+
+
+def _u_from_thetavec(
+    theta: np.typing.NDArray[np.float64] | Sequence[float], /
+) -> np.typing.NDArray[np.complex128]:
+    theta = np.asarray(theta)
+    if theta.shape != (3,):
+        raise ValueError(
+            f"theta vector has wrong shape: {theta.shape} (1D vector of length 3 expected)"
+        )
+    eigvals = np.array(
+        [
+            -np.sum(theta),
+            -theta[0] + theta[1] + theta[2],
+            -theta[1] + theta[2] + theta[0],
+            -theta[2] + theta[0] + theta[1],
+        ]
+    )
+    eigvecs = np.ones([1, 1]) / 2 - np.eye(4)
+    return np.transpose(eigvecs) @ (np.exp(1j * eigvals) * eigvecs[:, 0])
+
+
+def _nonlocal_qpd_basis_from_u(
+    u: np.typing.NDArray[np.complex128] | Sequence[complex], /
+) -> QPDBasis:
+    u = np.asarray(u)
+    if u.shape != (4,):
+        raise ValueError(
+            f"u vector has wrong shape: {u.shape} (1D vector of length 4 expected)"
+        )
+    # The following operations are described in Sec. 2.3 of
+    # https://quantum-journal.org/papers/q-2021-01-28-388/
+    #
+    # Projective measurements in each basis
+    A0x = [HGate(), QPDMeasure(), HGate()]
+    A0y = [SdgGate(), HGate(), QPDMeasure(), HGate(), SGate()]
+    A0z = [QPDMeasure()]
+    # Single qubit rotations that swap two axes.  There are "plus" and "minus"
+    # versions of these rotations.  The "minus" rotations also flip the sign
+    # along that axis.
+    Axyp = [SGate(), YGate()]
+    Axym = [ZGate()] + Axyp
+    Ayzp = [SXGate(), ZGate()]
+    Ayzm = [XGate()] + Ayzp
+    Azxp = [HGate()]
+    Azxm = [YGate()] + Azxp
+    # Single qubit rotations by ±pi/4 about each axis.
+    B0xp = [SXGate()]
+    B0xm = [SXdgGate()]
+    B0yp = [RYGate(0.5 * np.pi)]
+    B0ym = [RYGate(-0.5 * np.pi)]
+    B0zp = [SGate()]
+    B0zm = [SdgGate()]
+    # Projective measurements, each followed by the proper flip.
+    Bxy = A0z + [XGate()]
+    Byz = A0x + [YGate()]
+    Bzx = A0y + [ZGate()]
+    # The following values occur repeatedly in the coefficients
+    uu01 = u[0] * np.conj(u[1])
+    uu02 = u[0] * np.conj(u[2])
+    uu03 = u[0] * np.conj(u[3])
+    uu12 = u[1] * np.conj(u[2])
+    uu23 = u[2] * np.conj(u[3])
+    uu31 = u[3] * np.conj(u[1])
+    coeffs, maps1, maps2 = zip(
+        # First line of Eq. (19) in
+        # https://quantum-journal.org/papers/q-2021-01-28-388/
+        (np.abs(u[0]) ** 2, [], []),  # Identity
+        (np.abs(u[1]) ** 2, [XGate()], [XGate()]),
+        (np.abs(u[2]) ** 2, [YGate()], [YGate()]),
+        (np.abs(u[3]) ** 2, [ZGate()], [ZGate()]),
+        # Second line
+        (2 * np.real(uu01), A0x, A0x),
+        (2 * np.real(uu02), A0y, A0y),
+        (2 * np.real(uu03), A0z, A0z),
+        (0.5 * np.real(uu12), Axyp, Axyp),
+        (-0.5 * np.real(uu12), Axyp, Axym),
+        (-0.5 * np.real(uu12), Axym, Axyp),
+        (0.5 * np.real(uu12), Axym, Axym),
+        (0.5 * np.real(uu23), Ayzp, Ayzp),
+        (-0.5 * np.real(uu23), Ayzp, Ayzm),
+        (-0.5 * np.real(uu23), Ayzm, Ayzp),
+        (0.5 * np.real(uu23), Ayzm, Ayzm),
+        (0.5 * np.real(uu31), Azxp, Azxp),
+        (-0.5 * np.real(uu31), Azxp, Azxm),
+        (-0.5 * np.real(uu31), Azxm, Azxp),
+        (0.5 * np.real(uu31), Azxm, Azxm),
+        (-0.5 * np.real(uu01), B0xp, B0xp),
+        (0.5 * np.real(uu01), B0xp, B0xm),
+        (0.5 * np.real(uu01), B0xm, B0xp),
+        (-0.5 * np.real(uu01), B0xm, B0xm),
+        (-0.5 * np.real(uu02), B0yp, B0yp),
+        (0.5 * np.real(uu02), B0yp, B0ym),
+        (0.5 * np.real(uu02), B0ym, B0yp),
+        (-0.5 * np.real(uu02), B0ym, B0ym),
+        (-0.5 * np.real(uu03), B0zp, B0zp),
+        (0.5 * np.real(uu03), B0zp, B0zm),
+        (0.5 * np.real(uu03), B0zm, B0zp),
+        (-0.5 * np.real(uu03), B0zm, B0zm),
+        (-2 * np.real(uu12), Bxy, Bxy),
+        (-2 * np.real(uu23), Byz, Byz),
+        (-2 * np.real(uu31), Bzx, Bzx),
+        # Third line
+        (np.imag(uu01), A0x, B0xp),
+        (-np.imag(uu01), A0x, B0xm),
+        (np.imag(uu01), B0xp, A0x),
+        (-np.imag(uu01), B0xm, A0x),
+        (np.imag(uu02), A0y, B0yp),
+        (-np.imag(uu02), A0y, B0ym),
+        (np.imag(uu02), B0yp, A0y),
+        (-np.imag(uu02), B0ym, A0y),
+        (np.imag(uu03), A0z, B0zp),
+        (-np.imag(uu03), A0z, B0zm),
+        (np.imag(uu03), B0zp, A0z),
+        (-np.imag(uu03), B0zm, A0z),
+        (np.imag(uu12), Axyp, Bxy),
+        (-np.imag(uu12), Axym, Bxy),
+        (np.imag(uu12), Bxy, Axyp),
+        (-np.imag(uu12), Bxy, Axym),
+        (np.imag(uu23), Ayzp, Byz),
+        (-np.imag(uu23), Ayzm, Byz),
+        (np.imag(uu23), Byz, Ayzp),
+        (-np.imag(uu23), Byz, Ayzm),
+        (np.imag(uu31), Azxp, Bzx),
+        (-np.imag(uu31), Azxm, Bzx),
+        (np.imag(uu31), Bzx, Azxp),
+        (-np.imag(uu31), Bzx, Azxm),
+    )
+    maps = list(zip(maps1, _copy_unique_sublists(maps2)))
+    return QPDBasis(maps, coeffs)
+
+
+@_register_qpdbasis_from_gate("swap")
+def _(gate: SwapGate):
+    return _nonlocal_qpd_basis_from_u([(1 + 1j) / np.sqrt(8)] * 4)
+
+
+@_register_qpdbasis_from_gate("iswap")
+def _(gate: iSwapGate):
+    return _nonlocal_qpd_basis_from_u([0.5, 0.5j, 0.5j, 0.5])
+
+
+@_register_qpdbasis_from_gate("dcx")
+def _(gate: DCXGate):
+    retval = qpdbasis_from_gate(iSwapGate())
+    # Modify basis according to DCXGate definition in Qiskit circuit library
+    # https://github.com/Qiskit/qiskit-terra/blob/e9f8b7c50968501e019d0cb426676ac606eb5a10/qiskit/circuit/library/standard_gates/equivalence_library.py#L938-L944
+    for operations in unique_by_id(m[0] for m in retval.maps):
+        operations.insert(0, SdgGate())
+        operations.insert(0, HGate())
+    for operations in unique_by_id(m[1] for m in retval.maps):
+        operations.insert(0, SdgGate())
+        operations.append(HGate())
+    return retval
+
+
 @_register_qpdbasis_from_gate("rxx", "ryy", "rzz", "crx", "cry", "crz")
 def _(gate: RXXGate | RYYGate | RZZGate | CRXGate | CRYGate | CRZGate):
     # Constructing a virtual two-qubit gate by sampling single-qubit operations - Mitarai et al

releasenotes/notes/additional-gates-f4ed6c0e8dc3a9be.yaml

@@ -6,3 +6,5 @@ features:
     - :class:`~qiskit.circuit.library.CYGate`
     - :class:`~qiskit.circuit.library.SXGate`
     - :class:`~qiskit.circuit.library.ECRGate`
+    - :class:`~qiskit.circuit.library.SwapGate`
+    - :class:`~qiskit.circuit.library.iSwapGate`

test/cutting/qpd/test_qpd.py

@@ -28,6 +28,7 @@
     RXXGate,
     RYYGate,
     RZZGate,
+    RZXGate,
 )
 
 from circuit_knitting.utils.iteration import unique_by_eq
@@ -38,6 +39,10 @@
     generate_qpd_samples,
 )
 from circuit_knitting.cutting.qpd.qpd import *
+from circuit_knitting.cutting.qpd.qpd import (
+    _nonlocal_qpd_basis_from_u,
+    _u_from_thetavec,
+)
 
 
 @ddt
@@ -248,6 +253,10 @@ def test_decompose_qpd_instructions(self):
         (RXXGate(np.pi / 7), 1 + 2 * np.abs(np.sin(np.pi / 7))),
         (RYYGate(np.pi / 7), 1 + 2 * np.abs(np.sin(np.pi / 7))),
         (RZZGate(np.pi / 7), 1 + 2 * np.abs(np.sin(np.pi / 7))),
+        (RZXGate(np.pi / 7), 1 + 2 * np.abs(np.sin(np.pi / 7))),
+        (SwapGate(), 7),
+        (iSwapGate(), 7),
+        (DCXGate(), 7),
     )
     @unpack
     def test_optimal_kappa_for_known_gates(self, instruction, gamma):
@@ -299,6 +308,35 @@ def test_supported_gates(self):
                 "ch",
                 "csx",
                 "ecr",
+                "swap",
+                "iswap",
+                "dcx",
             },
             gates,
         )
+
+    def test_nonlocal_qpd_basis_from_u(self):
+        with self.subTest("Invalid shape"):
+            with pytest.raises(ValueError) as e_info:
+                _nonlocal_qpd_basis_from_u([1, 2, 3])
+            assert (
+                e_info.value.args[0]
+                == "u vector has wrong shape: (3,) (1D vector of length 4 expected)"
+            )
+
+    @data(
+        ([np.pi / 4] * 3, [(1 + 1j) / np.sqrt(8)] * 4),
+        ([np.pi / 4, np.pi / 4, 0], [0.5, 0.5j, 0.5j, 0.5]),
+    )
+    @unpack
+    def test_u_from_thetavec(self, theta, expected):
+        assert _u_from_thetavec(theta) == pytest.approx(expected)
+
+    def test_u_from_thetavec_exceptions(self):
+        with self.subTest("Invalid shape"):
+            with pytest.raises(ValueError) as e_info:
+                _u_from_thetavec([0, 1, 2, 3])
+            assert (
+                e_info.value.args[0]
+                == "theta vector has wrong shape: (4,) (1D vector of length 3 expected)"
+            )

test/cutting/qpd/test_qpd_basis.py

@@ -120,8 +120,8 @@ def test_eq(self):
 
     def test_unsupported_gate(self):
         with pytest.raises(ValueError) as e_info:
-            QPDBasis.from_gate(SwapGate())
-        assert e_info.value.args[0] == "Gate not supported: swap"
+            QPDBasis.from_gate(C3XGate())
+        assert e_info.value.args[0] == "Gate not supported: mcx"
 
     def test_unbound_parameter(self):
         with pytest.raises(ValueError) as e_info:

test/cutting/test_cutting_roundtrip.py

@@ -21,15 +21,21 @@
     RXXGate,
     RYYGate,
     RZZGate,
+    RZXGate,
+    XXPlusYYGate,
+    XXMinusYYGate,
     CHGate,
     CXGate,
     CYGate,
     CZGate,
     CRXGate,
     CRYGate,
     CRZGate,
-    ECRGate,
     CSXGate,
+    ECRGate,
+    SwapGate,
+    iSwapGate,
+    DCXGate,
 )
 from qiskit.extensions import UnitaryGate
 from qiskit.quantum_info import PauliList, random_unitary
@@ -52,6 +58,9 @@ def append_random_unitary(circuit: QuantumCircuit, qubits):
 
 @pytest.fixture(
     params=[
+        [SwapGate()],
+        [iSwapGate()],
+        [DCXGate()],
         [CXGate()],
         [CYGate()],
         [CZGate()],
@@ -70,6 +79,10 @@ def append_random_unitary(circuit: QuantumCircuit, qubits):
         [CRYGate(np.pi / 7)],
         [CRZGate(np.pi / 11)],
         [RXXGate(np.pi / 3), CRYGate(np.pi / 7)],
+        [UnitaryGate(random_unitary(2**2))],
+        [RZXGate(np.pi / 5)],
+        [XXPlusYYGate(7 * np.pi / 11)],
+        [XXMinusYYGate(11 * np.pi / 17)],
     ]
 )
 def example_circuit(
@@ -126,7 +139,7 @@ def test_cutting_exact_reconstruction(example_circuit):
     quasi_dists, coefficients = execute_experiments(
         circuits=subcircuits,
         subobservables=subobservables,
-        num_samples=1500,
+        num_samples=np.inf,
         samplers=sampler,
     )
     simulated_expvals = reconstruct_expectation_values(