Implement cutting of general 2-qubit unitaries

circuit_knitting/cutting/__init__.py

@@ -57,7 +57,7 @@
     qpd.generate_qpd_weights
     qpd.generate_qpd_samples
     qpd.decompose_qpd_instructions
-    qpd.supported_gates
+    qpd.explicitly_supported_gates
 
 CutQC
 =====

circuit_knitting/cutting/qpd/__init__.py

@@ -18,7 +18,7 @@
     decompose_qpd_instructions,
     WeightType,
     qpdbasis_from_gate,
-    supported_gates,
+    explicitly_supported_gates,
 )
 from .instructions import (
     BaseQPDGate,
@@ -32,7 +32,7 @@
     "generate_qpd_weights",
     "generate_qpd_samples",
     "decompose_qpd_instructions",
-    "supported_gates",
+    "explicitly_supported_gates",
     "QPDBasis",
     "BaseQPDGate",
     "TwoQubitQPDGate",

circuit_knitting/cutting/qpd/qpd.py

@@ -62,6 +62,8 @@
     iSwapGate,
     DCXGate,
 )
+from qiskit.extensions import UnitaryGate
+from qiskit.quantum_info.synthesis.two_qubit_decompose import TwoQubitWeylDecomposition
 from qiskit.utils import deprecate_func
 
 from .qpd_basis import QPDBasis
@@ -557,45 +559,52 @@ def qpdbasis_from_gate(gate: Gate) -> QPDBasis:
     """
     Generate a QPDBasis object, given a supported operation.
 
-    This method currently supports the following operations:
-        - :class:`~qiskit.circuit.library.RXXGate`
-        - :class:`~qiskit.circuit.library.RYYGate`
-        - :class:`~qiskit.circuit.library.RZZGate`
-        - :class:`~qiskit.circuit.library.CRXGate`
-        - :class:`~qiskit.circuit.library.CRYGate`
-        - :class:`~qiskit.circuit.library.CRZGate`
-        - :class:`~qiskit.circuit.library.CXGate`
-        - :class:`~qiskit.circuit.library.CYGate`
-        - :class:`~qiskit.circuit.library.CZGate`
-        - :class:`~qiskit.circuit.library.CHGate`
-        - :class:`~qiskit.circuit.library.CSXGate`
-        - :class:`~qiskit.circuit.library.CSGate`
-        - :class:`~qiskit.circuit.library.CSdgGate`
-        - :class:`~qiskit.circuit.library.CPhaseGate`
-        - :class:`~qiskit.circuit.library.SwapGate`
-        - :class:`~qiskit.circuit.library.iSwapGate`
-        - :class:`~qiskit.circuit.library.DCXGate`
-
-    The above gate names can also be determined by calling
-    :func:`supported_gates`.
+    The operations with explicit support can be obtained by calling
+    :func:`explicitly_supported_gates`.
+
+    Additionally, all two-qubit gates which implement the :meth:`.Gate.to_matrix` method are
+    supported via a KAK decomposition (:class:`.TwoQubitWeylDecomposition`).
 
     Returns:
         The newly-instantiated :class:`QPDBasis` object
 
     Raises:
+        ValueError: Gate not supported.
         ValueError: Cannot decompose gate with unbound parameters.
+        ValueError: ``to_matrix`` conversion of two-qubit gate failed.
     """
     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:
+        try:
+            mat = gate.to_matrix()
+        except Exception as ex:
+            raise ValueError("`to_matrix` conversion of two-qubit gate failed") from ex
+        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]:
+
+def explicitly_supported_gates() -> set[str]:
     """
-    Return a set of gate names supported for automatic decomposition.
+    Return a set of instruction names with explicit support for automatic decomposition.
+
+    These instructions are *explicitly* supported by :func:`qpdbasis_from_gate`.
+    Other instructions may be supported too, via a KAK decomposition.
 
     Returns:
         Set of gate names supported for automatic decomposition.
@@ -619,6 +628,54 @@ def _copy_unique_sublists(lsts: tuple[list, ...], /) -> tuple[list, ...]:
     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]:
+    r"""
+    Exponentiate the non-local portion of a KAK decomposition.
+
+    This implements Eq. (6) of https://arxiv.org/abs/2006.11174v2:
+
+    .. math::
+
+       \exp [ i ( \sum_\alpha^3 \theta_\alpha \, \sigma_\alpha \otimes \sigma_\alpha ) ]
+       =
+       \sum_{\alpha=0}^3 u_\alpha \, \sigma_\alpha \otimes \sigma_\alpha
+
+    where each :math:`\theta_\alpha` is assumed to be real, and
+    :math:`u_\alpha` is complex in general.
+    """
+    theta = np.asarray(theta)
+    if theta.shape != (3,):
+        raise ValueError(
+            f"theta vector has wrong shape: {theta.shape} (1D vector of length 3 expected)"
+        )
+    # First, we note that if we choose the basis vectors II, XX, YY, and ZZ,
+    # then the following matrix represents one application of the summation in
+    # the exponential:
+    #
+    #   0   θx  θy  θz
+    #   θx  0  -θz -θy
+    #   θy -θz  0  -θx
+    #   θz -θy -θx  0
+    #
+    # This matrix is symmetric and can be exponentiated by diagonalizing it.
+    # Its eigendecomposition is given by:
+    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)
+    # Finally, we exponentiate the eigenvalues of the matrix in diagonal form.
+    # We also project to the vector [1,0,0,0] on the right, since the
+    # multiplicative identity is given by II.
+    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:

releasenotes/notes/supported-gates-d2156f58bc07fc7a.yaml

@@ -1,4 +1,4 @@
 ---
 upgrade:
   - |
-    Addition of :func:`~circuit_knitting.cutting.qpd.supported_gates` function, which returns the names of all gates which may be automatically decomposed using :func:`~circuit_knitting.cutting.qpd.qpdbasis_from_gate`.
+    Addition of :func:`~circuit_knitting.cutting.qpd.explicitly_supported_gates` function, which returns the names of all gates which are explicitly supported by :func:`~circuit_knitting.cutting.qpd.qpdbasis_from_gate`.

test/cutting/qpd/test_qpd.py

@@ -31,6 +31,7 @@
     RXXGate,
     RYYGate,
     RZZGate,
+    RZXGate,
 )
 
 from circuit_knitting.utils.iteration import unique_by_eq
@@ -45,6 +46,7 @@
     _generate_qpd_weights,
     _generate_exact_weights_and_conditional_probabilities,
     _nonlocal_qpd_basis_from_u,
+    _u_from_thetavec,
 )
 
 
@@ -256,6 +258,7 @@ 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))),
         (CPhaseGate(np.pi / 7), 1 + 2 * np.abs(np.sin(np.pi / 14))),
         (CSGate(), 1 + np.sqrt(2)),
         (CSdgGate(), 1 + np.sqrt(2)),
@@ -422,8 +425,8 @@ def from_theta(theta):
             )
             assert weights[map_ids][1] == WeightType.SAMPLED
 
-    def test_supported_gates(self):
-        gates = supported_gates()
+    def test_explicitly_supported_gates(self):
+        gates = explicitly_supported_gates()
         self.assertEqual(
             {
                 "rxx",
@@ -456,3 +459,30 @@ def test_nonlocal_qpd_basis_from_u(self):
                 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)"
+            )
+
+    def test_qpdbasis_from_gate_errors(self):
+        with self.subTest("to_matrix fails"):
+            with pytest.raises(ValueError) as e_info:
+                # https://github.com/Qiskit/qiskit-terra/issues/10396
+                qpdbasis_from_gate(CSXGate().inverse())
+            assert (
+                e_info.value.args[0]
+                == "`to_matrix` conversion of two-qubit gate failed"
+            )

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(XXMinusYYGate(0.1))
-        assert e_info.value.args[0] == "Gate not supported: xx_minus_yy"
+            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,6 +21,9 @@
     RXXGate,
     RYYGate,
     RZZGate,
+    RZXGate,
+    XXPlusYYGate,
+    XXMinusYYGate,
     CHGate,
     CXGate,
     CYGate,
@@ -83,6 +86,10 @@ def append_random_unitary(circuit: QuantumCircuit, qubits):
         [RXXGate(np.pi / 3), CRYGate(np.pi / 7)],
         [CPhaseGate(np.pi / 3)],
         [RXXGate(np.pi / 3), CPhaseGate(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(