Use pauli_rep to compute qml.matrix() whenever possible

doc/releases/changelog-dev.md

@@ -223,13 +223,14 @@
   [stim](https://github.com/quantumlib/Stim) `v1.13.0`.
   [(#5409)](https://github.com/PennyLaneAI/pennylane/pull/5409)
 
+* `Sum.matrix` and `Prod.matrix` use the pauli representation for matrix computation,
+  which results in faster computation of the matrix.
+  [(5392)](https://github.com/PennyLaneAI/pennylane/pull/5392)
+
 * `qml.specs` and `qml.Tracker` now return information about algorithmic errors for the qnode as well.
   [(#5464)](https://github.com/PennyLaneAI/pennylane/pull/5464)
   [(#5465)](https://github.com/PennyLaneAI/pennylane/pull/5465)
 
-* `qml.specs` now returns information regarding algorithmic errors for the qnode as well.
-  [(#5464)](https://github.com/PennyLaneAI/pennylane/pull/5464)
-
 * `qml.transforms.hamiltonian_expand` can now handle multi-term observables with a constant offset.
   [(#5414)](https://github.com/PennyLaneAI/pennylane/pull/5414)
 

pennylane/ops/functions/matrix.py

@@ -211,8 +211,10 @@ def circuit():
     if isinstance(op, qml.operation.Tensor) and wire_order is not None:
         op = 1.0 * op  # convert to a Hamiltonian
 
-    if isinstance(op, qml.ops.Hamiltonian):
+    if (pr := op.pauli_rep) is not None:
+        return pr.to_mat(wire_order=wire_order or op.wires, format="csr").toarray()
 
+    if isinstance(op, qml.ops.Hamiltonian):
         return op.sparse_matrix(wire_order=wire_order).toarray()
 
     try:

pennylane/ops/op_math/prod.py

@@ -295,6 +295,9 @@ def decomposition(self):
     def matrix(self, wire_order=None):
         """Representation of the operator as a matrix in the computational basis."""
 
+        # if self.pauli_rep:  # Get the matrix from the PauliSentence representation
+        #     return self.pauli_rep.to_mat(wire_order=wire_order or self.wires, format="csr").toarray()
+
         mats: List[TensorLike] = []
         batched: List[bool] = []  # batched[i] tells if mats[i] is batched or not
         for ops in self.overlapping_ops:

pennylane/ops/op_math/sum.py

@@ -318,6 +318,9 @@ def matrix(self, wire_order=None):
         Returns:
             tensor_like: matrix representation
         """
+        # if self.pauli_rep:  # Get the matrix from the PauliSentence representation
+        #     return self.pauli_rep.to_mat(wire_order=wire_order or self.wires, format="csr").toarray()
+
         gen = (
             (qml.matrix(op) if isinstance(op, qml.ops.Hamiltonian) else op.matrix(), op.wires)
             for op in self

pennylane/pauli/pauli_arithmetic.py

@@ -226,7 +226,7 @@ def __hash__(self):
     def _matmul(self, other):
         """Private matrix multiplication that returns (pauli_word, coeff) tuple for more lightweight processing"""
         base, iterator, swapped = (
-            (self, other, False) if len(self) > len(other) else (other, self, True)
+            (self, other, False) if len(self) >= len(other) else (other, self, True)
         )
         result = copy(dict(base))
         coeff = 1

tests/circuit_graph/test_qasm.py

@@ -117,9 +117,9 @@ def test_to_ApproxTimeEvolution(self):
             include "qelib1.inc";
             qreg q[2];
             creg c[2];
-            cx q[0],q[1];
-            rz(2.0) q[1];
-            cx q[0],q[1];
+            cx q[1],q[0];
+            rz(2.0) q[0];
+            cx q[1],q[0];
             measure q[0] -> c[0];
             measure q[1] -> c[1];
             """

tests/ops/functions/test_matrix.py

@@ -569,8 +569,8 @@ def circuit():
             qml.matrix(circuit, wire_order=wires)()
 
         with pytest.raises(
-            TransformError,
-            match=r"Wires in circuit \[0\] are inconsistent with those in wire_order \[1\]",
+            ValueError,
+            match=r"the matrix for the specified wire order because it does not contain all the Pauli sentence",
         ):
             qml.matrix(qml.PauliX(0), wire_order=[1])
 

tests/ops/op_math/test_prod.py

@@ -255,6 +255,15 @@ def test_terms_pauli_rep(self, op, coeffs_true, ops_true):
         assert coeffs == coeffs_true
         assert ops1 == ops_true
 
+    def test_terms_pauli_rep_wire_order(self):
+        """Test that the wire order of the terms is the same as the wire order of the original
+        operands when the Prod has a valid pauli_rep"""
+        H = qml.prod(X(0), X(1), X(2))
+        _, H_ops = H.terms()
+
+        assert len(H_ops) == 1
+        assert H_ops[0].wires == H.wires
+
     def test_batch_size(self):
         """Test that batch size returns the batch size of a base operation if it is batched."""
         x = qml.numpy.array([1.0, 2.0, 3.0])
@@ -1035,8 +1044,8 @@ def test_pauli_rep_order(self):
         """
         op = qml.prod(qml.PauliX(0), qml.PauliY(1), qml.PauliZ(2))
         pw = list(op.pauli_rep.keys())[0]
-        assert list(pw.keys()) == [1, 0, 2]
-        assert list(pw.values()) == ["Y", "X", "Z"]
+        assert list(pw.keys()) == [0, 1, 2]
+        assert list(pw.values()) == ["X", "Y", "Z"]
 
     @pytest.mark.parametrize("op, rep", op_pauli_reps)
     def test_pauli_rep(self, op, rep):
@@ -1119,7 +1128,7 @@ def test_simplify_method_product_of_sums(self):
         """Test the simplify method with a product of sums."""
         prod_op = Prod(qml.PauliX(0) + qml.RX(1, 0), qml.PauliX(1) + qml.RX(1, 1), qml.Identity(3))
         final_op = qml.sum(
-            Prod(qml.PauliX(1), qml.PauliX(0)),
+            Prod(qml.PauliX(0), qml.PauliX(1)),
             qml.PauliX(0) @ qml.RX(1, 1),
             qml.PauliX(1) @ qml.RX(1, 0),
             qml.RX(1, 0) @ qml.RX(1, 1),

tests/ops/op_math/test_sum.py

@@ -212,6 +212,24 @@ def test_terms_pauli_rep(self, op, coeffs_true, ops_true):
         assert coeffs == coeffs_true
         assert ops1 == ops_true
 
+    def test_terms_pauli_rep_wire_order(self):
+        """Test that the wire order of the terms is the same as the wire order of the original
+        operands when the Sum has a valid pauli_rep"""
+        w0, w1, w2, w3 = [0, 1, 2, 3]
+        coeffs = [0.5, -0.5]
+
+        obs = [
+            qml.X(w0) @ qml.Y(w1) @ qml.X(w2) @ qml.Z(w3),
+            qml.X(w0) @ qml.X(w1) @ qml.Y(w2) @ qml.Z(w3),
+        ]
+
+        H = qml.dot(coeffs, obs)
+        _, H_ops = H.terms()
+
+        assert all(o1.wires == o2.wires for o1, o2 in zip(obs, H_ops))
+        assert H_ops[0] == qml.prod(qml.X(w0), qml.Y(w1), qml.X(w2), qml.Z(w3))
+        assert H_ops[1] == qml.prod(qml.X(w0), qml.X(w1), qml.Y(w2), qml.Z(w3))
+
     coeffs_ = [1.0, 1.0, 1.0, 3.0, 4.0, 4.0, 5.0]
     h6 = qml.sum(
         qml.s_prod(2.0, qml.prod(qml.Hadamard(0), qml.PauliZ(10))),

tests/pauli/test_measurement_transformations.py

@@ -120,7 +120,7 @@ def test_diagonalize_pauli_word_catch_non_pauli_word(self, non_pauli_word):
         (
             [PauliX(0) @ PauliY(1), PauliX(0) @ PauliZ(2)],
             (
-                [RX(np.pi / 2, wires=[1]), RY(-np.pi / 2, wires=[0])],
+                [RY(-np.pi / 2, wires=[0]), RX(np.pi / 2, wires=[1])],
                 [PauliZ(wires=[0]) @ PauliZ(wires=[1]), PauliZ(wires=[0]) @ PauliZ(wires=[2])],
             ),
         ),

tests/pauli/test_pauli_utils.py

@@ -926,10 +926,7 @@ def test_diagonalize_pauli_word_catch_non_pauli_word(self, non_pauli_word):
         (
             [PauliX(0) @ PauliY(1), PauliX(0) @ PauliZ(2)],
             (
-                [
-                    RX(np.pi / 2, wires=[1]),
-                    RY(-np.pi / 2, wires=[0]),
-                ],
+                [RY(-np.pi / 2, wires=[0]), RX(np.pi / 2, wires=[1])],
                 [PauliZ(wires=[0]) @ PauliZ(wires=[1]), PauliZ(wires=[0]) @ PauliZ(wires=[2])],
             ),
         ),

tests/templates/test_subroutines/test_approx_time_evolution.py

@@ -66,7 +66,7 @@ class TestDecomposition:
             ),
             (
                 2,
-                qml.Hamiltonian([2, 0.5], [qml.PauliX("a"), qml.PauliZ("b") @ qml.PauliX("a")]),
+                qml.Hamiltonian([2, 0.5], [qml.PauliX("a"), qml.PauliX("a") @ qml.PauliZ("b")]),
                 2,
                 [
                     qml.PauliRot(4.0, "X", wires=["a"]),
@@ -87,7 +87,7 @@ class TestDecomposition:
                     [2, 0.5, 0.5],
                     [
                         qml.PauliX("a"),
-                        qml.PauliZ(-15) @ qml.PauliX("a"),
+                        qml.PauliX("a") @ qml.PauliZ(-15),
                         qml.Identity(0) @ qml.PauliY(-15),
                     ],
                 ),

tests/test_qaoa.py

@@ -1166,12 +1166,12 @@ def make_mixer_layer_test_cases():
         [
             qaoa.xy_mixer(Graph([(0, 1), (1, 2), (2, 0)])),
             [
-                qml.PauliRot(1.0, "XX", wires=[1, 0]),
-                qml.PauliRot(1.0, "YY", wires=[1, 0]),
-                qml.PauliRot(1.0, "XX", wires=[2, 0]),
-                qml.PauliRot(1.0, "YY", wires=[2, 0]),
-                qml.PauliRot(1.0, "XX", wires=[2, 1]),
-                qml.PauliRot(1.0, "YY", wires=[2, 1]),
+                qml.PauliRot(1.0, "XX", wires=[0, 1]),
+                qml.PauliRot(1.0, "YY", wires=[0, 1]),
+                qml.PauliRot(1.0, "XX", wires=[0, 2]),
+                qml.PauliRot(1.0, "YY", wires=[0, 2]),
+                qml.PauliRot(1.0, "XX", wires=[1, 2]),
+                qml.PauliRot(1.0, "YY", wires=[1, 2]),
             ],
         ],
     ]
@@ -1186,9 +1186,9 @@ def make_cost_layer_test_cases():
         [
             qaoa.maxcut(Graph([(0, 1), (1, 2), (2, 0)]))[0],
             [
-                qml.PauliRot(1.0, "ZZ", wires=[1, 0]),
-                qml.PauliRot(1.0, "ZZ", wires=[2, 0]),
-                qml.PauliRot(1.0, "ZZ", wires=[2, 1]),
+                qml.PauliRot(1.0, "ZZ", wires=[0, 1]),
+                qml.PauliRot(1.0, "ZZ", wires=[0, 2]),
+                qml.PauliRot(1.0, "ZZ", wires=[1, 2]),
             ],
         ],
     ]