Repository: PennyLaneAI/pennylane
Pull request #2225 (2022-02-23 23:23:01)

add `pow` method to operators
This PR adds the `Operator.pow` method.  This method will provide easier ways to decompose an operator raised to a power. It will be purely developer-facing.  Users will raise operators to powers through either a `qml.pow` constructor or `op ** n`.

Later PR's will:
* Add a `Pow` wrapper arithmetic class
* Add a `qml.pow` contructor transform
* Add a `__pow__` method.

When the power does not exist for the provided exponent, it raises a `qml.operation.PowUndefinedError`.

github-actions[bot]: Hello. You may have forgotten to update the changelog!
Please edit [`doc/releases/changelog-dev.md`](/PennyLaneAI/pennylane/blob/master/doc/releases/changelog-dev.md) with:
* A one-to-two sentence description of the change. You may include a small working example for new features.
* A link back to this PR.
* Your name (or GitHub username) in the contributors section.
codecov[bot]: # [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2225?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#2225](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2225?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (2f04401) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/1fdbcb56e0695668f4e2987550e5999e06fbf640?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (1fdbcb5) will **decrease** coverage by `3.41%`.
> The diff coverage is `36.09%`.

[![Impacted file tree graph](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2225/graphs/tree.svg?width=650&height=150&src=pr&token=9S3YzlUgGk&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2225?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)

```diff
@@            Coverage Diff             @@
##           master    #2225      +/-   ##
==========================================
- Coverage   99.26%   95.84%   -3.42%     
==========================================
  Files         231      231              
  Lines       18221    18405     +184     
==========================================
- Hits        18087    17641     -446     
- Misses        134      764     +630     
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2225?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Coverage Δ | |
|---|---|---|
| [pennylane/ops/qubit/non\_parametric\_ops.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2225/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9xdWJpdC9ub25fcGFyYW1ldHJpY19vcHMucHk=) | `84.86% <17.64%> (-15.14%)` | :arrow_down: |
| [pennylane/ops/qubit/observables.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2225/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9xdWJpdC9vYnNlcnZhYmxlcy5weQ==) | `92.30% <20.00%> (-7.70%)` | :arrow_down: |
| [pennylane/ops/qubit/matrix\_ops.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2225/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9xdWJpdC9tYXRyaXhfb3BzLnB5) | `91.86% <25.00%> (-7.23%)` | :arrow_down: |
| [pennylane/ops/qubit/parametric\_ops.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2225/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9xdWJpdC9wYXJhbWV0cmljX29wcy5weQ==) | `97.33% <50.00%> (-2.67%)` | :arrow_down: |
| [pennylane/ops/qubit/qchem\_ops.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2225/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9xdWJpdC9xY2hlbV9vcHMucHk=) | `99.61% <50.00%> (-0.39%)` | :arrow_down: |
| [pennylane/operation.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2225/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wZXJhdGlvbi5weQ==) | `95.80% <66.66%> (-0.15%)` | :arrow_down: |
| [pennylane/ops/cv.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2225/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9jdi5weQ==) | `97.24% <73.80%> (-2.76%)` | :arrow_down: |
| [pennylane/interfaces/batch/jax.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2225/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2ludGVyZmFjZXMvYmF0Y2gvamF4LnB5) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| [pennylane/interfaces/batch/jax\_jit.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2225/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2ludGVyZmFjZXMvYmF0Y2gvamF4X2ppdC5weQ==) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| [pennylane/interfaces/jax.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2225/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2ludGVyZmFjZXMvamF4LnB5) | `0.00% <0.00%> (-97.44%)` | :arrow_down: |
| ... and [20 more](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2225/diff?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2225?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2225?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Last update [1fdbcb5...2f04401](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2225?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).

josh146: wow @albi3ro this is one powerhouse of a PR - I didn't realize you had added `__pow__` to every operator class! Even if this doesn't get merged in its current state (we still need to coordinate with @mariaschuld), this will be a good source of exponentiation relations between different operations.
albi3ro: @mariaschuld Thanks for all the good comments. I'll get to responding to the rest in a moment, but I first wanted to address why I decided to return a list of operations instead of a single operation.

First, this will be used in `Pow.compute_decomposition`.  Decompositions return lists, so even if `pow` returns a single operation, it will get promptly wrapped into a list in its one use case.

Second, how do we handle the identity use case? I was initially having the method return an `Identity`, but that identity is basically just a placeholder that we will have to deal with down the road somewhere.  Plenty of operations, like `Identity`, `Barrier`, `WireCut`, etc. decompose to an empty list, so there's nothing wrong with having a power operator decompose to an empty list.

What happens if the base gate is across multiple wires, like `CNOT`? Do we return an `Identity` on every wire, or just the first one? If it's on every wire, we need to return a list of operations anyway.  If it's just the first one, what's so special about that wire? It's no longer indicating all the places something used to be, so it's not really serving it's purpose as longer-term memory.

And lastly, what if we want to support more intricate rules later? For example, we could decide to decompose `SX **-1` as both `PauliX` and `SX` together. In that situation, we would need to return a list instead of just a single operation.
albi3ro: @josh146 So maybe my last point about longer decomposition might not hold as well in the future when we have a `Prod` class, but I still think the other arguments hold water.

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 703f81390..417ae63e4 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -39,6 +39,10 @@
 
 <h3>Improvements</h3>
 
+* The developer-facing `pow` method has been added to `Operator` with concrete implementations
+  for many classes.
+  [(#2225)](https://github.com/PennyLaneAI/pennylane/pull/2225)
+
 * Test classes are created in qchem test modules to group the integrals and matrices unittests.
   [(#2545)](https://github.com/PennyLaneAI/pennylane/pull/2545)
 
diff --git a/pennylane/operation.py b/pennylane/operation.py
index 561dad215..82389e930 100644
--- a/pennylane/operation.py
+++ b/pennylane/operation.py
@@ -255,6 +255,10 @@ class AdjointUndefinedError(OperatorPropertyUndefined):
     """Raised when an Operator's adjoint version is undefined."""
 
 
+class PowUndefinedError(OperatorPropertyUndefined):
+    """Raised when an Operator's power is undefined."""
+
+
 class GeneratorUndefinedError(OperatorPropertyUndefined):
     """Exception used to indicate that an operator
     does not have a generator"""
@@ -1025,6 +1029,25 @@ class Operator(abc.ABC):
         """
         raise GeneratorUndefinedError(f"Operation {self.name} does not have a generator")
 
+    def pow(self, z):
+        """A list of new operators equal to this one raised to the given power.
+
+        Args:
+            z (float): exponent for the operator
+
+        Returns:
+            list[:class:`~.operation.Operator`]
+
+        """
+        # Child methods may call super().pow(z%period) where op**period = I
+        # For example, PauliX**2 = I, SX**4 = I
+        # Hence we define 0 and 1 special cases here.
+        if z == 0:
+            return []
+        if z == 1:
+            return [self.__copy__()]
+        raise PowUndefinedError
+
     def queue(self, context=qml.QueuingContext):
         """Append the operator to the Operator queue."""
         context.append(self)
diff --git a/pennylane/ops/identity.py b/pennylane/ops/identity.py
index a5483b70e..12d895cb2 100644
--- a/pennylane/ops/identity.py
+++ b/pennylane/ops/identity.py
@@ -148,3 +148,6 @@ class Identity(CVObservable, Operation):
 
     def adjoint(self):  # pylint:disable=arguments-differ
         return Identity(wires=self.wires)
+
+    def pow(self, n):
+        return [Identity(wires=self.wires)]
diff --git a/pennylane/ops/qubit/matrix_ops.py b/pennylane/ops/qubit/matrix_ops.py
index da4508c9a..b67ed2455 100644
--- a/pennylane/ops/qubit/matrix_ops.py
+++ b/pennylane/ops/qubit/matrix_ops.py
@@ -153,6 +153,11 @@ class QubitUnitary(Operation):
     def adjoint(self):
         return QubitUnitary(qml.math.T(qml.math.conj(self.matrix())), wires=self.wires)
 
+    def pow(self, z):
+        if isinstance(z, int):
+            return [QubitUnitary(qml.math.linalg.matrix_power(self.matrix(), z), wires=self.wires)]
+        return super().pow(z)
+
     def _controlled(self, wire):
         ControlledQubitUnitary(*self.parameters, control_wires=wire, wires=self.wires)
 
@@ -317,6 +322,17 @@ class ControlledQubitUnitary(QubitUnitary):
     def control_wires(self):
         return self.hyperparameters["control_wires"]
 
+    def pow(self, z):
+        if isinstance(z, int):
+            return [
+                ControlledQubitUnitary(
+                    qml.math.linalg.matrix_power(self.data[0], z),
+                    control_wires=self.control_wires,
+                    wires=self.hyperparameters["u_wires"],
+                )
+            ]
+        return super().pow(z)
+
     def _controlled(self, wire):
         ctrl_wires = sorted(self.control_wires + wire)
         ControlledQubitUnitary(
@@ -439,6 +455,12 @@ class DiagonalQubitUnitary(Operation):
     def adjoint(self):
         return DiagonalQubitUnitary(qml.math.conj(self.parameters[0]), wires=self.wires)
 
+    def pow(self, z):
+        if isinstance(self.data[0], list):
+            return [DiagonalQubitUnitary([(x + 0.0j) ** z for x in self.data[0]], wires=self.wires)]
+        casted_data = qml.math.cast(self.data[0], np.complex128)
+        return [DiagonalQubitUnitary(casted_data**z, wires=self.wires)]
+
     def _controlled(self, control):
         DiagonalQubitUnitary(
             qml.math.concatenate([np.ones_like(self.parameters[0]), self.parameters[0]]),
diff --git a/pennylane/ops/qubit/non_parametric_ops.py b/pennylane/ops/qubit/non_parametric_ops.py
index 3c6d542e1..fa96baab1 100644
--- a/pennylane/ops/qubit/non_parametric_ops.py
+++ b/pennylane/ops/qubit/non_parametric_ops.py
@@ -158,6 +158,9 @@ class Hadamard(Observable, Operation):
         # H = RZ(\pi) RY(\pi/2) RZ(0)
         return [np.pi, np.pi / 2, 0.0]
 
+    def pow(self, z):
+        return super().pow(z % 2)
+
 
 class PauliX(Observable, Operation):
     r"""PauliX(wires)
@@ -288,6 +291,12 @@ class PauliX(Observable, Operation):
     def adjoint(self):
         return PauliX(wires=self.wires)
 
+    def pow(self, z):
+        z_mod2 = z % 2
+        if abs(z_mod2 - 0.5) < 1e-6:
+            return [SX(wires=self.wires)]
+        return super().pow(z_mod2)
+
     def _controlled(self, wire):
         CNOT(wires=Wires(wire) + self.wires)
 
@@ -427,6 +436,9 @@ class PauliY(Observable, Operation):
     def adjoint(self):
         return PauliY(wires=self.wires)
 
+    def pow(self, z):
+        return super().pow(z % 2)
+
     def _controlled(self, wire):
         CY(wires=Wires(wire) + self.wires)
 
@@ -554,6 +566,20 @@ class PauliZ(Observable, Operation):
     def adjoint(self):
         return PauliZ(wires=self.wires)
 
+    def pow(self, z):
+        z_mod2 = z % 2
+        if z_mod2 == 0:
+            return []
+        if z_mod2 == 1:
+            return [self.__copy__()]
+
+        if abs(z_mod2 - 0.5) < 1e-6:
+            return [S(wires=self.wires)]
+        if abs(z_mod2 - 0.25) < 1e-6:
+            return [T(wires=self.wires)]
+
+        return [qml.PhaseShift(np.pi * z_mod2, wires=self.wires)]
+
     def _controlled(self, wire):
         CZ(wires=Wires(wire) + self.wires)
 
@@ -658,6 +684,18 @@ class S(Operation):
         op.inverse = not self.inverse
         return op
 
+    def pow(self, z):
+        z_mod4 = z % 4
+        pow_map = {
+            0: lambda op: [],
+            0.5: lambda op: [T(wires=op.wires)],
+            1: lambda op: [op.__copy__()],
+            2: lambda op: [PauliZ(wires=op.wires)],
+        }
+        return pow_map.get(z_mod4, lambda op: [qml.PhaseShift(np.pi * z_mod4 / 2, wires=op.wires)])(
+            self
+        )
+
     def single_qubit_rot_angles(self):
         # S = RZ(\pi/2) RY(0) RZ(0)
         return [np.pi / 2, 0.0, 0.0]
@@ -754,6 +792,18 @@ class T(Operation):
         """
         return [qml.PhaseShift(np.pi / 4, wires=wires)]
 
+    def pow(self, z):
+        z_mod8 = z % 8
+        pow_map = {
+            0: lambda op: [],
+            1: lambda op: [op.__copy__()],
+            2: lambda op: [S(wires=op.wires)],
+            4: lambda op: [PauliZ(wires=op.wires)],
+        }
+        return pow_map.get(z_mod8, lambda op: [qml.PhaseShift(np.pi * z_mod8 / 4, wires=op.wires)])(
+            self
+        )
+
     def adjoint(self):
         op = T(wires=self.wires)
         op.inverse = not self.inverse
@@ -865,6 +915,12 @@ class SX(Operation):
         ]
         return decomp_ops
 
+    def pow(self, z):
+        z_mod4 = z % 4
+        if z_mod4 == 2:
+            return [PauliX(wires=self.wires)]
+        return super().pow(z_mod4)
+
     def adjoint(self):
         op = SX(wires=self.wires)
         op.inverse = not self.inverse
@@ -931,6 +987,9 @@ class CNOT(Operation):
     def adjoint(self):
         return CNOT(wires=self.wires)
 
+    def pow(self, z):
+        return super().pow(z % 2)
+
     def _controlled(self, wire):
         Toffoli(wires=Wires(wire) + self.wires)
 
@@ -1020,6 +1079,9 @@ class CZ(Operation):
     def adjoint(self):
         return CZ(wires=self.wires)
 
+    def pow(self, z):
+        return super().pow(z % 2)
+
     @property
     def control_wires(self):
         return Wires(self.wires[0])
@@ -1112,6 +1174,9 @@ class CY(Operation):
     def adjoint(self):
         return CY(wires=self.wires)
 
+    def pow(self, z):
+        return super().pow(z % 2)
+
     @property
     def control_wires(self):
         return Wires(self.wires[0])
@@ -1190,6 +1255,9 @@ class SWAP(Operation):
         ]
         return decomp_ops
 
+    def pow(self, z):
+        return super().pow(z % 2)
+
     def adjoint(self):
         return SWAP(wires=self.wires)
 
@@ -1309,6 +1377,12 @@ class ISWAP(Operation):
         op.inverse = not self.inverse
         return op
 
+    def pow(self, z):
+        z_mod2 = z % 2
+        if abs(z_mod2 - 0.5) < 1e-6:
+            return [SISWAP(wires=self.wires)]
+        return super().pow(z_mod2)
+
 
 class SISWAP(Operation):
     r"""SISWAP(wires)
@@ -1437,6 +1511,10 @@ class SISWAP(Operation):
         ]
         return decomp_ops
 
+    def pow(self, z):
+        z_mod4 = z % 4
+        return [ISWAP(wires=self.wires)] if z_mod4 == 2 else super().pow(z_mod4)
+
     def adjoint(self):
         op = SISWAP(wires=self.wires)
         op.inverse = not self.inverse
@@ -1544,6 +1622,9 @@ class CSWAP(Operation):
         ]
         return decomp_ops
 
+    def pow(self, z):
+        return super().pow(z % 2)
+
     def adjoint(self):
         return CSWAP(wires=self.wires)
 
@@ -1682,6 +1763,9 @@ class Toffoli(Operation):
     def adjoint(self):
         return Toffoli(wires=self.wires)
 
+    def pow(self, z):
+        return super().pow(z % 2)
+
     @property
     def control_wires(self):
         return Wires(self.wires[:2])
@@ -1854,6 +1938,9 @@ class MultiControlledX(Operation):
             control_values=self.hyperparameters["control_values"],
         )
 
+    def pow(self, z):
+        return super().pow(z % 2)
+
     @staticmethod
     def compute_decomposition(wires=None, work_wires=None, control_values=None, **kwargs):
         r"""Representation of the operator as a product of other operators (static method).
@@ -2056,6 +2143,9 @@ class Barrier(Operation):
     def adjoint(self, do_queue=True):
         return Barrier(wires=self.wires, do_queue=do_queue)
 
+    def pow(self, z):
+        return [self.__copy__()]
+
 
 class WireCut(Operation):
     r"""WireCut(wires)
@@ -2103,3 +2193,6 @@ class WireCut(Operation):
 
     def adjoint(self):
         return WireCut(wires=self.wires)
+
+    def pow(self, z):
+        return [self.__copy__()]
diff --git a/pennylane/ops/qubit/observables.py b/pennylane/ops/qubit/observables.py
index 1d3eecefe..0e70e2d3d 100644
--- a/pennylane/ops/qubit/observables.py
+++ b/pennylane/ops/qubit/observables.py
@@ -441,3 +441,6 @@ class Projector(Observable):
         []
         """
         return []
+
+    def pow(self, z):
+        return [self.__copy__()] if (isinstance(z, int) and z > 0) else super().pow(z)
diff --git a/pennylane/ops/qubit/parametric_ops.py b/pennylane/ops/qubit/parametric_ops.py
index 11c073870..09c39fbd2 100644
--- a/pennylane/ops/qubit/parametric_ops.py
+++ b/pennylane/ops/qubit/parametric_ops.py
@@ -106,6 +106,9 @@ class RX(Operation):
     def adjoint(self):
         return RX(-self.data[0], wires=self.wires)
 
+    def pow(self, z):
+        return [RX(self.data[0] * z, wires=self.wires)]
+
     def _controlled(self, wire):
         CRX(*self.parameters, wires=wire + self.wires)
 
@@ -184,6 +187,9 @@ class RY(Operation):
     def adjoint(self):
         return RY(-self.data[0], wires=self.wires)
 
+    def pow(self, z):
+        return [RY(self.data[0] * z, wires=self.wires)]
+
     def _controlled(self, wire):
         CRY(*self.parameters, wires=wire + self.wires)
 
@@ -295,6 +301,9 @@ class RZ(Operation):
     def adjoint(self):
         return RZ(-self.data[0], wires=self.wires)
 
+    def pow(self, z):
+        return [RZ(self.data[0] * z, wires=self.wires)]
+
     def _controlled(self, wire):
         CRZ(*self.parameters, wires=wire + self.wires)
 
@@ -433,6 +442,9 @@ class PhaseShift(Operation):
     def adjoint(self):
         return PhaseShift(-self.data[0], wires=self.wires)
 
+    def pow(self, z):
+        return [PhaseShift(self.data[0] * z, wires=self.wires)]
+
     def _controlled(self, wire):
         ControlledPhaseShift(*self.parameters, wires=wire + self.wires)
 
@@ -588,6 +600,9 @@ class ControlledPhaseShift(Operation):
     def adjoint(self):
         return ControlledPhaseShift(-self.data[0], wires=self.wires)
 
+    def pow(self, z):
+        return [ControlledPhaseShift(self.data[0] * z, wires=self.wires)]
+
     @property
     def control_wires(self):
         return Wires(self.wires[0])
@@ -879,6 +894,9 @@ class MultiRZ(Operation):
     def adjoint(self):
         return MultiRZ(-self.parameters[0], wires=self.wires)
 
+    def pow(self, z):
+        return [MultiRZ(self.data[0] * z, wires=self.wires)]
+
 
 class PauliRot(Operation):
     r"""
@@ -1165,6 +1183,9 @@ class PauliRot(Operation):
     def adjoint(self):
         return PauliRot(-self.parameters[0], self.hyperparameters["pauli_word"], wires=self.wires)
 
+    def pow(self, z):
+        return [PauliRot(self.data[0] * z, self.hyperparameters["pauli_word"], wires=self.wires)]
+
 
 class CRX(Operation):
     r"""
@@ -1310,6 +1331,9 @@ class CRX(Operation):
     def adjoint(self):
         return CRX(-self.data[0], wires=self.wires)
 
+    def pow(self, z):
+        return [CRX(self.data[0] * z, wires=self.wires)]
+
     @property
     def control_wires(self):
         return Wires(self.wires[0])
@@ -1449,6 +1473,9 @@ class CRY(Operation):
     def adjoint(self):
         return CRY(-self.data[0], wires=self.wires)
 
+    def pow(self, z):
+        return [CRY(self.data[0] * z, wires=self.wires)]
+
     @property
     def control_wires(self):
         return Wires(self.wires[0])
@@ -1616,6 +1643,9 @@ class CRZ(Operation):
     def adjoint(self):
         return CRZ(-self.data[0], wires=self.wires)
 
+    def pow(self, z):
+        return [CRZ(self.data[0] * z, wires=self.wires)]
+
     @property
     def control_wires(self):
         return Wires(self.wires[0])
@@ -1879,6 +1909,9 @@ class U1(Operation):
     def adjoint(self):
         return U1(-self.data[0], wires=self.wires)
 
+    def pow(self, z):
+        return [U1(self.data[0] * z, wires=self.wires)]
+
 
 class U2(Operation):
     r"""
@@ -2242,6 +2275,9 @@ class IsingXX(Operation):
         (phi,) = self.parameters
         return IsingXX(-phi, wires=self.wires)
 
+    def pow(self, z):
+        return [IsingXX(self.data[0] * z, wires=self.wires)]
+
 
 class IsingYY(Operation):
     r"""
@@ -2348,6 +2384,9 @@ class IsingYY(Operation):
         (phi,) = self.parameters
         return IsingYY(-phi, wires=self.wires)
 
+    def pow(self, z):
+        return [IsingYY(self.data[0] * z, wires=self.wires)]
+
 
 class IsingZZ(Operation):
     r"""
@@ -2485,3 +2524,6 @@ class IsingZZ(Operation):
     def adjoint(self):
         (phi,) = self.parameters
         return IsingZZ(-phi, wires=self.wires)
+
+    def pow(self, z):
+        return [IsingZZ(self.data[0] * z, wires=self.wires)]
diff --git a/pennylane/ops/qubit/qchem_ops.py b/pennylane/ops/qubit/qchem_ops.py
index 0d0513494..018959c00 100644
--- a/pennylane/ops/qubit/qchem_ops.py
+++ b/pennylane/ops/qubit/qchem_ops.py
@@ -154,6 +154,9 @@ class SingleExcitation(Operation):
         (phi,) = self.parameters
         return SingleExcitation(-phi, wires=self.wires)
 
+    def pow(self, z):
+        return [SingleExcitation(self.data[0] * z, wires=self.wires)]
+
     def label(self, decimals=None, base_label=None, cache=None):
         return super().label(decimals=decimals, base_label=base_label or "G", cache=cache)
 
@@ -515,6 +518,9 @@ class DoubleExcitation(Operation):
         ]
         return qml.Hamiltonian(coeffs, obs)
 
+    def pow(self, z):
+        return [DoubleExcitation(self.data[0] * z, wires=self.wires)]
+
     def __init__(self, phi, wires, do_queue=True, id=None):
         super().__init__(phi, wires=wires, do_queue=do_queue, id=id)
 


Test patch:
diff --git a/tests/ops/qubit/test_matrix_ops.py b/tests/ops/qubit/test_matrix_ops.py
index 3366e7809..1a03f97dc 100644
--- a/tests/ops/qubit/test_matrix_ops.py
+++ b/tests/ops/qubit/test_matrix_ops.py
@@ -34,6 +34,35 @@ from gate_data import (
 class TestQubitUnitary:
     """Tests for the QubitUnitary class."""
 
+    def test_qubit_unitary_noninteger_pow(self):
+        """Test QubitUnitary raised to a non-integer power raises an error."""
+        U = np.array(
+            [[0.98877108 + 0.0j, 0.0 - 0.14943813j], [0.0 - 0.14943813j, 0.98877108 + 0.0j]]
+        )
+
+        op = qml.QubitUnitary(U, wires="a")
+
+        with pytest.raises(qml.operation.PowUndefinedError):
+            op.pow(0.123)
+
+    @pytest.mark.parametrize("n", (1, 3, -1, -3))
+    def test_qubit_unitary_pow(self, n):
+        """Test qubit unitary raised to an integer power."""
+        U = np.array(
+            [[0.98877108 + 0.0j, 0.0 - 0.14943813j], [0.0 - 0.14943813j, 0.98877108 + 0.0j]]
+        )
+
+        op = qml.QubitUnitary(U, wires="a")
+        new_ops = op.pow(n)
+
+        assert len(new_ops) == 1
+        assert new_ops[0].wires == op.wires
+
+        mat_to_pow = qml.math.linalg.matrix_power(qml.matrix(op), n)
+        new_mat = qml.matrix(new_ops[0])
+
+        assert qml.math.allclose(mat_to_pow, new_mat)
+
     @pytest.mark.autograd
     @pytest.mark.parametrize("U,num_wires", [(H, 1), (np.kron(H, H), 2)])
     def test_qubit_unitary_autograd(self, U, num_wires):
@@ -256,6 +285,17 @@ class TestDiagonalQubitUnitary:
         assert np.allclose(res_static, expected, atol=tol)
         assert np.allclose(res_dynamic, expected, atol=tol)
 
+    @pytest.mark.parametrize("n", (2, -1, 0.12345))
+    @pytest.mark.parametrize("diag", ([1.0, -1.0], np.array([1.0, -1.0])))
+    def test_pow(self, n, diag):
+        """Test pow method returns expected results."""
+        op = qml.DiagonalQubitUnitary(diag, wires="b")
+        pow_ops = op.pow(n)
+        assert len(pow_ops) == 1
+
+        for x_op, x_pow in zip(op.data[0], pow_ops[0].data[0]):
+            assert (x_op + 0.0j) ** n == x_pow
+
     def test_error_matrix_not_unitary(self):
         """Tests that error is raised if diagonal by `compute_matrix` does not lead to a unitary"""
         with pytest.raises(ValueError, match="Operator must be unitary"):
@@ -539,6 +579,41 @@ class TestControlledQubitUnitary:
         with pytest.raises(qml.operation.DecompositionUndefinedError):
             qml.ControlledQubitUnitary(mat, wires=0, control_wires=1).decomposition()
 
+    @pytest.mark.parametrize("n", (2, -1, -2))
+    def test_pow(self, n):
+        """Tests the metadata and unitary for a controlledQubitUnitary raised to a power."""
+        U1 = np.array(
+            [
+                [0.73708696 + 0.61324932j, 0.27034258 + 0.08685028j],
+                [-0.24979544 - 0.1350197j, 0.95278437 + 0.1075819j],
+            ]
+        )
+
+        op = qml.ControlledQubitUnitary(U1, control_wires=("b", "c"), wires="a")
+
+        pow_ops = op.pow(n)
+        assert len(pow_ops) == 1
+
+        assert pow_ops[0].hyperparameters["u_wires"] == op.hyperparameters["u_wires"]
+        assert pow_ops[0].control_wires == op.control_wires
+
+        op_mat_to_pow = qml.math.linalg.matrix_power(op.data[0], n)
+        assert qml.math.allclose(pow_ops[0].data[0], op_mat_to_pow)
+
+    def test_noninteger_pow(self):
+        """Test that a ControlledQubitUnitary raised to a non-integer power raises an error."""
+        U1 = np.array(
+            [
+                [0.73708696 + 0.61324932j, 0.27034258 + 0.08685028j],
+                [-0.24979544 - 0.1350197j, 0.95278437 + 0.1075819j],
+            ]
+        )
+
+        op = qml.ControlledQubitUnitary(U1, control_wires=("b", "c"), wires="a")
+
+        with pytest.raises(qml.operation.PowUndefinedError):
+            op.pow(0.12)
+
 
 label_data = [
     (X, qml.QubitUnitary(X, wires=0)),
diff --git a/tests/ops/qubit/test_non_parametric_ops.py b/tests/ops/qubit/test_non_parametric_ops.py
index 2dfd3424d..7201ce028 100644
--- a/tests/ops/qubit/test_non_parametric_ops.py
+++ b/tests/ops/qubit/test_non_parametric_ops.py
@@ -981,6 +981,162 @@ class TestMultiControlledX:
         assert np.allclose(mat1, mat2)
 
 
+period_two_ops = (
+    qml.PauliX(0),
+    qml.PauliY(0),
+    qml.PauliZ(0),
+    qml.Hadamard("a"),
+    qml.CNOT(wires=(0, 1)),
+    qml.CZ(wires=(0, 1)),
+    qml.CY(wires=(0, 1)),
+    qml.SWAP(wires=(0, 1)),
+    qml.ISWAP(wires=(0, 1)),
+    qml.CSWAP(wires=(0, 1, 2)),
+    qml.Toffoli(wires=(0, 1, 2)),
+    qml.MultiControlledX(wires=(0, 1, 2, 3)),
+)
+
+
+class TestPowMethod:
+    @pytest.mark.parametrize("op", period_two_ops)
+    @pytest.mark.parametrize("n", (1, 5, -1, -5))
+    def test_period_two_pow_odd(self, op, n):
+        """Test that ops with a period of 2 raised to an odd power are the same as the original op."""
+        assert op.pow(n)[0].__class__ is op.__class__
+
+    @pytest.mark.parametrize("op", period_two_ops)
+    @pytest.mark.parametrize("n", (2, 6, 0, -2))
+    def test_period_two_pow_even(self, op, n):
+        """Test that ops with a period of 2 raised to an even power are empty lists."""
+        assert len(op.pow(n)) == 0
+
+    @pytest.mark.parametrize("op", period_two_ops)
+    def test_period_two_noninteger_power(self, op):
+        """Test that ops with a period of 2 raised to a non-integer power raise an error."""
+        if op.__class__ is qml.PauliZ:
+            pytest.skip("PauliZ can be raised to any power.")
+        with pytest.raises(qml.operation.PowUndefinedError):
+            op.pow(1.234)
+
+    @pytest.mark.parametrize("n", (0.5, 2.5, -1.5))
+    def test_paulix_squareroot(self, n):
+        """Check that the square root of PauliX is SX"""
+        op = qml.PauliX(0)
+
+        pow_ops = op.pow(n)
+        assert len(pow_ops) == 1
+        assert pow_ops[0].__class__ is qml.SX
+
+        sqrt_mat = qml.matrix(op.pow)(n)
+        sqrt_mat_squared = qml.math.linalg.matrix_power(sqrt_mat, 2)
+
+        assert qml.math.allclose(sqrt_mat_squared, qml.matrix(op))
+
+    @pytest.mark.parametrize("n", (0.5, 2.5, -1.5))
+    def test_pauliz_squareroot(self, n):
+        """Check that the square root of PauliZ is S"""
+        assert qml.PauliZ(0).pow(n)[0].__class__ is qml.S
+
+        op = qml.PauliZ(0)
+        sqrt_mat = qml.matrix(op.pow)(n)
+        sqrt_mat_squared = qml.math.linalg.matrix_power(sqrt_mat, 2)
+
+        assert qml.math.allclose(sqrt_mat_squared, qml.matrix(op))
+
+    @pytest.mark.parametrize("n", (0.25, 2.25, -1.75))
+    def test_pauliz_fourthroot(self, n):
+        """Check that the fourth root of PauliZ is T."""
+        assert qml.PauliZ(0).pow(n)[0].__class__ is qml.T
+
+        op = qml.PauliZ(0)
+        quad_mat = qml.matrix(op.pow)(n)
+        quad_mat_pow = qml.math.linalg.matrix_power(quad_mat, 4)
+
+        assert qml.math.allclose(quad_mat_pow, qml.matrix(op))
+
+    @pytest.mark.parametrize("n", (0.12, -3.462, 3.693))
+    def test_pauliz_general_power(self, n):
+        """Check that PauliZ raised to an non-integer power that's not the square root
+        results in a PhaseShift."""
+        op_pow = qml.PauliZ(0).pow(n)
+
+        assert len(op_pow) == 1
+        assert op_pow[0].__class__ is qml.PhaseShift
+        assert qml.math.allclose(op_pow[0].data[0], np.pi * (n % 2))
+
+    @pytest.mark.parametrize("n", (0.5, 2.5, -1.5))
+    def test_ISWAP_sqaure_root(self, n):
+        """Test that SISWAP is the square root of ISWAP."""
+        op = qml.ISWAP(wires=(0, 1))
+
+        assert op.pow(n)[0].__class__ is qml.SISWAP
+
+        sqrt_mat = qml.matrix(op.pow)(n)
+        sqrt_mat_squared = qml.math.linalg.matrix_power(sqrt_mat, 2)
+        assert qml.math.allclose(sqrt_mat_squared, qml.matrix(op))
+
+    @pytest.mark.parametrize("offset", (0, 4, -4))
+    def test_S_pow(self, offset):
+        op = qml.S("a")
+
+        assert len(op.pow(0 + offset)) == 0
+
+        assert op.pow(0.5 + offset)[0].__class__ is qml.T
+        assert op.pow(1 + offset)[0].__class__ is qml.S
+        assert op.pow(2 + offset)[0].__class__ is qml.PauliZ
+
+        n = 1.234
+        op_pow = op.pow(n + offset)
+        assert op_pow[0].__class__ is qml.PhaseShift
+        assert qml.math.allclose(op_pow[0].data[0], np.pi * n / 2)
+
+    @pytest.mark.parametrize("offset", (0, 8, -8))
+    def test_T_pow(self, offset):
+        """Test the powers of the T gate."""
+        op = qml.T("b")
+
+        assert len(op.pow(0 + offset)) == 0
+        assert op.pow(1 + offset)[0].__class__ is qml.T
+        assert op.pow(2 + offset)[0].__class__ is qml.S
+        assert op.pow(4 + offset)[0].__class__ is qml.PauliZ
+
+        n = 1.234
+        op_pow = op.pow(n + offset)
+        assert op_pow[0].__class__ is qml.PhaseShift
+        assert qml.math.allclose(op_pow[0].data[0], np.pi * n / 4)
+
+    @pytest.mark.parametrize("offset", (0, 4, -4))
+    def test_SX_pow(self, offset):
+        op = qml.SX("d")
+
+        assert len(op.pow(0 + offset)) == 0
+
+        assert op.pow(1 + offset)[0].__class__ is qml.SX
+        assert op.pow(2 + offset)[0].__class__ is qml.PauliX
+
+        with pytest.raises(qml.operation.PowUndefinedError):
+            op.pow(2.43 + offset)
+
+    @pytest.mark.parametrize("offset", (0, 4, -4))
+    def test_SISWAP_pow(self, offset):
+        """Test powers of the SISWAP operator"""
+        op = qml.SISWAP(wires=("b", "c"))
+
+        assert len(op.pow(0 + offset)) == 0
+        assert op.pow(1 + offset)[0].__class__ is qml.SISWAP
+        assert op.pow(2 + offset)[0].__class__ is qml.ISWAP
+
+        with pytest.raises(qml.operation.PowUndefinedError):
+            op.pow(2.34 + offset)
+
+    @pytest.mark.parametrize("op", (qml.WireCut(0), qml.Barrier(0)))
+    @pytest.mark.parametrize("n", (2, 0.123, -2.3))
+    def test_pow_independent_ops(self, op, n):
+        """Assert that the pow-independent ops WireCut and Barrier can be raised
+        to any power and just return a copy."""
+        assert op.pow(n)[0].__class__ is op.__class__
+
+
 label_data = [
     (qml.Identity(0), "I", "I"),
     (qml.Hadamard(0), "H", "H"),
diff --git a/tests/ops/qubit/test_observables.py b/tests/ops/qubit/test_observables.py
index cbfb522a0..ee6c6a987 100644
--- a/tests/ops/qubit/test_observables.py
+++ b/tests/ops/qubit/test_observables.py
@@ -472,6 +472,22 @@ class TestProjector:
         assert np.allclose(res_dynamic, expected, atol=tol)
         assert np.allclose(res_static, expected, atol=tol)
 
+    def test_pow_zero(self):
+        """Assert that the projector raised to zero is an empty list."""
+
+        basis_state = np.array([0, 1])
+        op = qml.Projector(basis_state, wires=(0, 1))
+        assert len(op.pow(0)) == 0
+
+    @pytest.mark.parametrize("n", (1, 3))
+    def test_pow_non_zero_positive_int(self, n):
+        """Test that the projector raised to a positive integer is just a copy."""
+        basis_state = np.array([0, 1])
+        op = qml.Projector(basis_state, wires=(0, 1))
+        pow_op = op.pow(n)[0]
+        assert pow_op.__class__ is qml.Projector
+        assert qml.math.allclose(pow_op.data[0], op.data[0])
+
 
 label_data = [
     (qml.Hermitian(np.eye(2), wires=1), "𝓗"),
diff --git a/tests/ops/qubit/test_parametric_ops.py b/tests/ops/qubit/test_parametric_ops.py
index 4f6b2dff0..0fbda53a6 100644
--- a/tests/ops/qubit/test_parametric_ops.py
+++ b/tests/ops/qubit/test_parametric_ops.py
@@ -1787,6 +1787,46 @@ class TestLabel:
         assert op3.label(decimals=0) == "Rot\n(x,\ny,\nz)"
 
 
+pow_parametric_ops = (
+    qml.RX(1.234, wires=0),
+    qml.RY(2.345, wires=0),
+    qml.RZ(3.456, wires=0),
+    qml.PhaseShift(6.78, wires=0),
+    qml.ControlledPhaseShift(0.234, wires=(0, 1)),
+    qml.MultiRZ(-0.4432, wires=(0, 1, 2)),
+    qml.PauliRot(0.5, "X", wires=0),
+    qml.CRX(-6.5432, wires=(0, 1)),
+    qml.CRY(-0.543, wires=(0, 1)),
+    qml.CRZ(1.234, wires=(0, 1)),
+    qml.U1(1.23, wires=0),
+    qml.IsingXX(-2.345, wires=(0, 1)),
+    qml.IsingYY(3.1652, wires=(0, 1)),
+    qml.IsingZZ(1.789, wires=("a", "b")),
+)
+
+
+class TestParametricPow:
+    @pytest.mark.parametrize("op", pow_parametric_ops)
+    @pytest.mark.parametrize("n", (2, -1, 0.2631, -0.987))
+    def test_pow_method_parametric_ops(self, op, n):
+        """Assert that a matrix raised to a power is the same as multiplying the data by n for relevant ops."""
+        pow_op = op.pow(n)
+
+        assert len(pow_op) == 1
+        assert pow_op[0].__class__ is op.__class__
+        assert all((d1 == d2 * n for d1, d2 in zip(pow_op[0].data, op.data)))
+
+    @pytest.mark.parametrize("op", pow_parametric_ops)
+    @pytest.mark.parametrize("n", (3, -2))
+    def test_pow_matrix(self, op, n):
+        """Test that the matrix of an op first raised to a power is the same as the
+        matrix raised to the power.  This test only can work for integer powers."""
+        op_mat = qml.matrix(op)
+        pow_mat = qml.matrix(op.pow)(n)
+
+        assert qml.math.allclose(qml.math.linalg.matrix_power(op_mat, n), pow_mat)
+
+
 control_data = [
     (qml.Rot(1, 2, 3, wires=0), Wires([])),
     (qml.RX(1.23, wires=0), Wires([])),
diff --git a/tests/ops/qubit/test_qchem_ops.py b/tests/ops/qubit/test_qchem_ops.py
index 0ea35f1ba..0a7d6608f 100644
--- a/tests/ops/qubit/test_qchem_ops.py
+++ b/tests/ops/qubit/test_qchem_ops.py
@@ -16,7 +16,7 @@ Unit tests for the available qubit operations for quantum chemistry purposes.
 """
 import pytest
 import numpy as np
-from scipy.linalg import expm
+from scipy.linalg import expm, fractional_matrix_power
 
 import pennylane as qml
 from pennylane import numpy as pnp
@@ -199,6 +199,21 @@ class TestSingleExcitation:
         exp = SingleExcitation(phi)
         assert np.allclose(res, exp)
 
+    @pytest.mark.parametrize("n", (2, -2, 1.3, -0.6))
+    def test_single_excitatation_pow(self, n):
+
+        op = qml.SingleExcitation(1.234, wires=(0, 1))
+
+        pow_ops = op.pow(n)
+        assert len(pow_ops) == 1
+        assert pow_ops[0].__class__ is qml.SingleExcitation
+
+        mat = qml.matrix(op)
+        pow_mat = qml.matrix(op.pow)(n)
+
+        mat_then_pow = fractional_matrix_power(mat, n)
+        assert qml.math.allclose(pow_mat, mat_then_pow)
+
     @pytest.mark.parametrize("phi", [-0.1, 0.2, np.pi / 4])
     def test_single_excitation_plus_matrix(self, phi):
         """Tests that the SingleExcitationPlus operation calculates the correct matrix"""
@@ -373,6 +388,21 @@ class TestDoubleExcitation:
 
         assert np.allclose(res, exp)
 
+    @pytest.mark.parametrize("n", (2, -2, 1.3, -0.6))
+    def test_double_excitatation_pow(self, n):
+        """Test the double excitation pow method."""
+        op = qml.DoubleExcitation(1.234, wires=(0, 1, 2, 3))
+
+        pow_ops = op.pow(n)
+        assert len(pow_ops) == 1
+        assert pow_ops[0].__class__ is qml.DoubleExcitation
+
+        mat = qml.matrix(op)
+        pow_mat = qml.matrix(op.pow)(n)
+
+        mat_then_pow = fractional_matrix_power(mat, n)
+        assert qml.math.allclose(pow_mat, mat_then_pow)
+
     @pytest.mark.parametrize("phi", [-0.1, 0.2, 0.5])
     def test_double_excitation_decomp(self, phi):
         """Tests that the DoubleExcitation operation calculates the correct decomposition.
diff --git a/tests/ops/test_identity.py b/tests/ops/test_identity.py
index 21598e87c..f6aea18ad 100644
--- a/tests/ops/test_identity.py
+++ b/tests/ops/test_identity.py
@@ -12,6 +12,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 """Unit tests for the Identity Operator."""
+import pytest
 
 from pennylane import Identity
 import numpy as np
@@ -36,6 +37,16 @@ def test_label_method():
     assert Identity(wires=0).label() == "I"
 
 
+@pytest.mark.parametrize("n", (2, -3, 3.455, -1.29))
+def test_identity_pow(n):
+    """Test that the identity raised to any power is simply a single copy."""
+    op = Identity("b")
+    pow_ops = op.pow(n)
+    assert len(pow_ops) == 1
+    assert pow_ops[0].__class__ is Identity
+    assert pow_ops[0].wires == op.wires
+
+
 def test_matrix_representation(tol):
     """Test the matrix representation"""
     res_static = Identity.compute_matrix()
diff --git a/tests/test_operation.py b/tests/test_operation.py
index 0e23afb81..bee4221c6 100644
--- a/tests/test_operation.py
+++ b/tests/test_operation.py
@@ -1389,6 +1389,21 @@ class TestDefaultRepresentations:
         with pytest.raises(qml.operation.GeneratorUndefinedError):
             gate.generator()
 
+    def test_pow_zero(self):
+        """Test that the default of an operation raised to a zero power is an empty array."""
+        assert len(gate.pow(0)) == 0
+
+    def test_pow_one(self):
+        """Test that the default of an operation raised to the power of one is a copy."""
+        pow_gate = gate.pow(1)
+        assert len(pow_gate) == 1
+        assert pow_gate[0].__class__ is gate.__class__
+
+    def test_pow_undefined(self):
+        """Tests that custom error is raised in the default pow decomposition."""
+        with pytest.raises(qml.operation.PowUndefinedError):
+            gate.pow(1.234)
+
 
 class TestChannel:
     """Unit tests for the Channel class"""