Operator apply permutation

qiskit/quantum_info/operators/operator.py

@@ -198,6 +198,85 @@ def from_label(cls, label):
                 op = op.compose(label_mats[char], qargs=[qubit])
         return op
 
+    def apply_permutation(self, perm: list, front: bool = False):
+        """Modifies operator's data by composing it with a permutation.
+
+        Args:
+            perm (list): permutation pattern, describing which qubits
+                occupy the positions 0, 1, 2, etc. after applying the permutation.
+            front (bool): When set to ``True`` the permutation is applied before the
+                operator, when set to ``False`` the permutation is applied after the
+                operator.
+        Returns:
+            Operator: The modified operator.
+
+        Raises:
+            QiskitError: if the size of the permutation pattern does not match the
+                dimensions of the operator.
+        """
+
+        # See https://github.com/Qiskit/qiskit-terra/pull/9403 for the math
+        # behind the following code.
+
+        inv_perm = np.argsort(perm)
+        raw_shape_l = self._op_shape.dims_l()
+        n_dims_l = len(raw_shape_l)
+        raw_shape_r = self._op_shape.dims_r()
+        n_dims_r = len(raw_shape_r)
+
+        if front:
+            # The permutation is applied first, the operator is applied after;
+            # however, in terms of matrices, we compute [O][P].
+
+            if len(perm) != n_dims_r:
+                raise QiskitError(
+                    "The size of the permutation pattern does not match dimensions of the operator."
+                )
+
+            # shape: original on left, permuted on right
+            shape_l = self._op_shape.dims_l()
+            shape_r = tuple(raw_shape_r[n_dims_r - n - 1] for n in reversed(perm))
+
+            # axes order: id on left, inv-permuted on right
+            axes_l = tuple(x for x in range(self._op_shape._num_qargs_l))
+            axes_r = tuple(self._op_shape._num_qargs_l + x for x in (np.argsort(perm[::-1]))[::-1])
+
+            # updated shape: original on left, permuted on right
+            new_shape_l = self._op_shape.dims_l()
+            new_shape_r = tuple(raw_shape_r[n_dims_r - n - 1] for n in reversed(inv_perm))
+
+        else:
+            # The operator is applied first, the permutation is applied after;
+            # however, in terms of matrices, we compute [P][O].
+
+            if len(perm) != n_dims_l:
+                raise QiskitError(
+                    "The size of the permutation pattern does not match dimensions of the operator."
+                )
+
+            # shape: inv-permuted on left, original on right
+            shape_l = tuple(raw_shape_l[n_dims_l - n - 1] for n in reversed(inv_perm))
+            shape_r = self._op_shape.dims_r()
+
+            # axes order: permuted on left, id on right
+            axes_l = tuple((np.argsort(inv_perm[::-1]))[::-1])
+            axes_r = tuple(
+                self._op_shape._num_qargs_l + x for x in range(self._op_shape._num_qargs_r)
+            )
+
+            # updated shape: permuted on left, original on right
+            new_shape_l = tuple(raw_shape_l[n_dims_l - n - 1] for n in reversed(perm))
+            new_shape_r = self._op_shape.dims_r()
+
+        # Computing the new operator
+        split_shape = shape_l + shape_r
+        axes_order = axes_l + axes_r
+        new_mat = (
+            self._data.reshape(split_shape).transpose(axes_order).reshape(self._op_shape.shape)
+        )
+        new_op = Operator(new_mat, input_dims=new_shape_r, output_dims=new_shape_l)
+        return new_op
+
     @classmethod
     def from_circuit(cls, circuit, ignore_set_layout=False, layout=None, final_layout=None):
         """Create a new Operator object from a :class:`.QuantumCircuit`
@@ -261,14 +340,8 @@ def from_circuit(cls, circuit, ignore_set_layout=False, layout=None, final_layou
         op._append_instruction(instruction, qargs=qargs)
         # If final layout is set permute output indices based on layout
         if final_layout is not None:
-            # TODO: Do this without the intermediate Permutation object by just
-            # operating directly on the array directly
-            from qiskit.circuit.library import Permutation  # pylint: disable=cyclic-import
-
-            final_physical_to_virtual = final_layout.get_physical_bits()
             perm_pattern = [final_layout._v2p[v] for v in circuit.qubits]
-            perm_op = Operator(Permutation(len(final_physical_to_virtual), perm_pattern))
-            op &= perm_op
+            op = op.apply_permutation(perm_pattern, front=False)
         return op
 
     def is_unitary(self, atol=None, rtol=None):

releasenotes/notes/operator-apply-permutation-c113c05513cb7881.yaml

@@ -0,0 +1,19 @@
+---
+features:
+  - |
+    Added a method :meth:`qiskit.quantum_info.Operator.apply_permutation` that
+    pre-composes or post-composes an Operator with a Permutation. This method
+    works for general qudits.
+
+    Here is an example to calculate :math:`P^\dagger.O.P` which reorders Operator's bits::
+
+        import numpy as np
+        from qiskit.quantum_info.operators import Operator
+
+        op = Operator(np.array(range(576)).reshape((24, 24)), input_dims=(2, 3, 4), output_dims=(2, 3, 4))
+        perm = [1, 2, 0]
+        inv_perm = [2, 0, 1]
+        conjugate_op = op.apply_permutation(inv_perm, front=True).apply_permutation(perm, front=False)
+
+    The conjugate operator has dimensions `(4, 2, 3) x (4, 2, 3)`, which is consistent with permutation
+    moving qutrit to position 0, qubit to position 1, and the 4-qudit to position 2.

test/python/quantum_info/operators/test_operator.py

@@ -17,7 +17,9 @@
 import unittest
 import logging
 import copy
+from test import combine
 import numpy as np
+from ddt import ddt
 from numpy.testing import assert_allclose
 import scipy.linalg as la
 
@@ -30,6 +32,7 @@
 from qiskit.quantum_info.operators.predicates import matrix_equal
 from qiskit.compiler.transpiler import transpile
 from qiskit.circuit import Qubit
+from qiskit.circuit.library import Permutation, PermutationGate
 
 logger = logging.getLogger(__name__)
 
@@ -91,6 +94,7 @@ def simple_circuit_with_measure(self):
         return circ
 
 
+@ddt
 class TestOperator(OperatorTestCase):
     """Tests for Operator linear operator class."""
 
@@ -1050,6 +1054,138 @@ def test_from_circuit_mixed_reg_loose_bits_transpiled(self):
         result = Operator.from_circuit(tqc)
         self.assertTrue(Operator(circuit).equiv(result))
 
+    def test_apply_permutation_back(self):
+        """Test applying permutation to the operator,
+        where the operator is applied first and the permutation second."""
+        op = Operator(self.rand_matrix(64, 64))
+        pattern = [1, 2, 0, 3, 5, 4]
+
+        # Consider several methods of computing this operator and show
+        # they all lead to the same result.
+
+        # Compose the operator with the operator constructed from the
+        # permutation circuit.
+        op2 = op.copy()
+        perm_op = Operator(Permutation(6, pattern))
+        op2 &= perm_op
+
+        # Compose the operator with the operator constructed from the
+        # permutation gate.
+        op3 = op.copy()
+        perm_op = Operator(PermutationGate(pattern))
+        op3 &= perm_op
+
+        # Modify the operator using apply_permutation method.
+        op4 = op.copy()
+        op4 = op4.apply_permutation(pattern, front=False)
+
+        self.assertEqual(op2, op3)
+        self.assertEqual(op2, op4)
+
+    def test_apply_permutation_front(self):
+        """Test applying permutation to the operator,
+        where the permutation is applied first and the operator second"""
+        op = Operator(self.rand_matrix(64, 64))
+        pattern = [1, 2, 0, 3, 5, 4]
+
+        # Consider several methods of computing this operator and show
+        # they all lead to the same result.
+
+        # Compose the operator with the operator constructed from the
+        # permutation circuit.
+        op2 = op.copy()
+        perm_op = Operator(Permutation(6, pattern))
+        op2 = perm_op & op2
+
+        # Compose the operator with the operator constructed from the
+        # permutation gate.
+        op3 = op.copy()
+        perm_op = Operator(PermutationGate(pattern))
+        op3 = perm_op & op3
+
+        # Modify the operator using apply_permutation method.
+        op4 = op.copy()
+        op4 = op4.apply_permutation(pattern, front=True)
+
+        self.assertEqual(op2, op3)
+        self.assertEqual(op2, op4)
+
+    def test_apply_permutation_qudits_back(self):
+        """Test applying permutation to the operator with heterogeneous qudit spaces,
+        where the operator O is applied first and the permutation P second.
+        The matrix of the resulting operator is the product [P][O] and
+        corresponds to suitably permuting the rows of O's matrix.
+        """
+        mat = np.array(range(6 * 6)).reshape((6, 6))
+        op = Operator(mat, input_dims=(2, 3), output_dims=(2, 3))
+        perm = [1, 0]
+        actual = op.apply_permutation(perm, front=False)
+
+        # Rows of mat are ordered to 00, 01, 02, 10, 11, 12;
+        # perm maps these to 00, 10, 20, 01, 11, 21,
+        # while the default ordering is 00, 01, 10, 11, 20, 21.
+        permuted_mat = mat.copy()[[0, 2, 4, 1, 3, 5]]
+        expected = Operator(permuted_mat, input_dims=(2, 3), output_dims=(3, 2))
+        self.assertEqual(actual, expected)
+
+    def test_apply_permutation_qudits_front(self):
+        """Test applying permutation to the operator with heterogeneous qudit spaces,
+        where the permutation P is applied first and the operator O is applied second.
+        The matrix of the resulting operator is the product [O][P] and
+        corresponds to suitably permuting the columns of O's matrix.
+        """
+        mat = np.array(range(6 * 6)).reshape((6, 6))
+        op = Operator(mat, input_dims=(2, 3), output_dims=(2, 3))
+        perm = [1, 0]
+        actual = op.apply_permutation(perm, front=True)
+
+        # Columns of mat are ordered to 00, 01, 02, 10, 11, 12;
+        # perm maps these to 00, 10, 20, 01, 11, 21,
+        # while the default ordering is 00, 01, 10, 11, 20, 21.
+        permuted_mat = mat.copy()[:, [0, 2, 4, 1, 3, 5]]
+        expected = Operator(permuted_mat, input_dims=(3, 2), output_dims=(2, 3))
+        self.assertEqual(actual, expected)
+
+    @combine(
+        dims=((2, 3, 4, 5), (5, 2, 4, 3), (3, 5, 2, 4), (5, 3, 4, 2), (4, 5, 2, 3), (4, 3, 2, 5))
+    )
+    def test_reverse_qargs_as_apply_permutation(self, dims):
+        """Test reversing qargs by pre- and post-composing with reversal
+        permutation.
+        """
+        perm = [3, 2, 1, 0]
+        op = Operator(
+            np.array(range(120 * 120)).reshape((120, 120)), input_dims=dims, output_dims=dims
+        )
+        op2 = op.reverse_qargs()
+        op3 = op.apply_permutation(perm, front=True).apply_permutation(perm, front=False)
+        self.assertEqual(op2, op3)
+
+    def test_apply_permutation_exceptions(self):
+        """Checks that applying permutation raises an error when dimensions do not match."""
+        op = Operator(
+            np.array(range(24 * 30)).reshape((24, 30)), input_dims=(6, 5), output_dims=(2, 3, 4)
+        )
+
+        with self.assertRaises(QiskitError):
+            op.apply_permutation([1, 0], front=False)
+        with self.assertRaises(QiskitError):
+            op.apply_permutation([2, 1, 0], front=True)
+
+    def test_apply_permutation_dimensions(self):
+        """Checks the dimensions of the operator after applying permutation."""
+        op = Operator(
+            np.array(range(24 * 30)).reshape((24, 30)), input_dims=(6, 5), output_dims=(2, 3, 4)
+        )
+        op2 = op.apply_permutation([1, 2, 0], front=False)
+        self.assertEqual(op2.output_dims(), (4, 2, 3))
+
+        op = Operator(
+            np.array(range(24 * 30)).reshape((30, 24)), input_dims=(2, 3, 4), output_dims=(6, 5)
+        )
+        op2 = op.apply_permutation([2, 0, 1], front=True)
+        self.assertEqual(op2.input_dims(), (4, 2, 3))
+
 
 if __name__ == "__main__":
     unittest.main()