Repository: PennyLaneAI/pennylane
Pull request #2198 (2022-02-11 19:34:02)

Denote matrices in circuit drawing
Currently, the new text and mpl drawer omit matrix-valued parameters from operator labels, while the old text drawer both includes them and prints them below the circuit diagram.

This PR adds a `cache` keyword to the `Operator.label` method.  This cache keyword is an optional dictionary that can store a list of encountered matrices.

The `cache` dictionary keyword of `qml.drawer.tape_text`  now stores a list of `'matrices'`

```python
with qml.tape.QuantumTape() as tape:
    qml.StronglyEntanglingLayers(params, wires=range(4))
    qml.QubitUnitary(U, wires=(0,1))
    qml.QubitUnitary(U, wires=(3,4))
print(qml.drawer.tape_text(tape, decimals=2))
```
```
0: ─╭StronglyEntanglingLayers(M0)─╭U(M1)─┤  
1: ─├StronglyEntanglingLayers(M0)─╰U(M1)─┤  
2: ─├StronglyEntanglingLayers(M0)────────┤  
3: ─╰StronglyEntanglingLayers(M0)─╭U(M1)─┤  
4: ───────────────────────────────╰U(M1)─┤  
```

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/2198?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#2198](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2198?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (7a5f1ec) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/d2b1077963053eecfe257270319f91b40d42d2a2?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (d2b1077) will **decrease** coverage by `2.78%`.
> The diff coverage is `99.21%`.

[![Impacted file tree graph](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2198/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/2198?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)

```diff
@@            Coverage Diff             @@
##           master    #2198      +/-   ##
==========================================
- Coverage   99.21%   96.43%   -2.79%     
==========================================
  Files         231      231              
  Lines       17916    17937      +21     
==========================================
- Hits        17775    17297     -478     
- Misses        141      640     +499     
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2198?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Coverage Δ | |
|---|---|---|
| [pennylane/operation.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2198/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wZXJhdGlvbi5weQ==) | `96.11% <95.83%> (+0.06%)` | :arrow_up: |
| [pennylane/drawer/tape\_text.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2198/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2RyYXdlci90YXBlX3RleHQucHk=) | `100.00% <100.00%> (ø)` | |
| [pennylane/ops/cv.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2198/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9jdi5weQ==) | `100.00% <100.00%> (ø)` | |
| [pennylane/ops/identity.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2198/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9pZGVudGl0eS5weQ==) | `100.00% <100.00%> (ø)` | |
| [pennylane/ops/qubit/arithmetic\_ops.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2198/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9xdWJpdC9hcml0aG1ldGljX29wcy5weQ==) | `100.00% <100.00%> (ø)` | |
| [pennylane/ops/qubit/hamiltonian.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2198/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9xdWJpdC9oYW1pbHRvbmlhbi5weQ==) | `99.04% <100.00%> (ø)` | |
| [pennylane/ops/qubit/matrix\_ops.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2198/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9xdWJpdC9tYXRyaXhfb3BzLnB5) | `99.12% <100.00%> (ø)` | |
| [pennylane/ops/qubit/non\_parametric\_ops.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2198/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9xdWJpdC9ub25fcGFyYW1ldHJpY19vcHMucHk=) | `97.35% <100.00%> (ø)` | |
| [pennylane/ops/qubit/observables.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2198/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9xdWJpdC9vYnNlcnZhYmxlcy5weQ==) | `100.00% <100.00%> (ø)` | |
| [pennylane/ops/qubit/parametric\_ops.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2198/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9xdWJpdC9wYXJhbWV0cmljX29wcy5weQ==) | `99.82% <100.00%> (ø)` | |
| ... and [25 more](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2198/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/2198?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/2198?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Last update [d2b1077...7a5f1ec](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2198?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).

albi3ro: [sc-14746]

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 78a419eac..60c36ef58 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -4,10 +4,17 @@
 
 <h3>New features since last release</h3>
 
-* The text based drawer accessed via `qml.draw` has been overhauled. The new drawer has 
-  a `decimals` keyword for controlling parameter rounding, a different algorithm for determining positions, 
-  deprecation of the `charset` keyword, and minor cosmetic changes.
+* The text based drawer accessed via `qml.draw` has been overhauled.
   [(#2128)](https://github.com/PennyLaneAI/pennylane/pull/2128)
+  [(#2198)](https://github.com/PennyLaneAI/pennylane/pull/2198)
+
+  The new drawer has:
+
+  * a `decimals` keyword for controlling parameter rounding
+  * a `show_matrices` keyword for controlling display of matrices
+  * a different algorithm for determining positions
+  * deprecation of the `charset` keyword
+  * additional minor cosmetic changes
 
   ```
   @qml.qnode(qml.device('lightning.qubit', wires=2))
diff --git a/pennylane/drawer/tape_text.py b/pennylane/drawer/tape_text.py
index 735111146..26d0f3fa4 100644
--- a/pennylane/drawer/tape_text.py
+++ b/pennylane/drawer/tape_text.py
@@ -36,7 +36,7 @@ def _add_grouping_symbols(op, layer_str, wire_map):
     return layer_str
 
 
-def _add_op(op, layer_str, wire_map, decimals):
+def _add_op(op, layer_str, wire_map, decimals, cache):
     """Updates ``layer_str`` with ``op`` operation."""
     layer_str = _add_grouping_symbols(op, layer_str, wire_map)
 
@@ -44,7 +44,7 @@ def _add_op(op, layer_str, wire_map, decimals):
     for w in control_wires:
         layer_str[wire_map[w]] += "C"
 
-    label = op.label(decimals=decimals).replace("\n", "")
+    label = op.label(decimals=decimals, cache=cache).replace("\n", "")
     for w in op.wires:
         if w not in control_wires:
             layer_str[wire_map[w]] += label
@@ -61,11 +61,14 @@ measurement_label_map = {
 }
 
 
-def _add_measurement(m, layer_str, wire_map, decimals):
+def _add_measurement(m, layer_str, wire_map, decimals, cache):
     """Updates ``layer_str`` with the ``m`` measurement."""
     layer_str = _add_grouping_symbols(m, layer_str, wire_map)
 
-    obs_label = None if m.obs is None else m.obs.label(decimals=decimals).replace("\n", "")
+    if m.obs is None:
+        obs_label = None
+    else:
+        obs_label = m.obs.label(decimals=decimals, cache=cache).replace("\n", "")
     meas_label = measurement_label_map[m.return_type](obs_label)
 
     if len(m.wires) == 0:  # state or probability across all wires
@@ -79,7 +82,13 @@ def _add_measurement(m, layer_str, wire_map, decimals):
 
 # pylint: disable=too-many-arguments
 def tape_text(
-    tape, wire_order=None, show_all_wires=False, decimals=None, max_length=100, cache=None
+    tape,
+    wire_order=None,
+    show_all_wires=False,
+    decimals=None,
+    max_length=100,
+    show_matrices=False,
+    cache=None,
 ):
     """Text based diagram for a Quantum Tape.
 
@@ -93,8 +102,9 @@ def tape_text(
             Default ``None`` will omit parameters from operation labels.
         max_length (Int) : Maximum length of a individual line.  After this length, the diagram will
             begin anew beneath the previous lines.
-        cache (dict): Used to store information between recursive calls. Currently only used for
-            numbering nested tapes.
+        show_matrices=False (bool): show matrix valued parameters below all circuit diagrams
+        cache (dict): Used to store information between recursive calls. Necessary keys are ``'tape_offset'``
+            and ``'matrices'``.
 
     Returns:
         str : String based graphic of the circuit.
@@ -177,6 +187,58 @@ def tape_text(
       1: ─├QFT──RY─╰C─┤ ╰Var[Z@Z] ├Probs
       2: ─╰QFT──RZ────┤           ╰Probs
 
+    Matrix valued parameters are always denoted by ``M`` followed by an integer corresponding to
+    unique matrices.  The list of unique matrices can be printed at the end of the diagram by
+    selecting ``show_matrices=True``:
+
+    .. code-block:: python
+
+        with qml.tape.QuantumTape() as tape:
+            qml.QubitUnitary(np.eye(2), wires=0)
+            qml.QubitUnitary(np.eye(2), wires=1)
+            qml.expval(qml.Hermitian(np.eye(4), wires=(0,1)))
+
+    >>> print(tape_text(tape, show_matrices=True))
+    0: ──U(M0)─┤ ╭<𝓗(M1)>
+    1: ──U(M0)─┤ ╰<𝓗(M1)>
+    M0 =
+    [[1. 0.]
+    [0. 1.]]
+    M1 =
+    [[1. 0. 0. 0.]
+    [0. 1. 0. 0.]
+    [0. 0. 1. 0.]
+    [0. 0. 0. 1.]]
+
+    An existing matrix cache can be passed via the ``cache`` keyword. Note that the dictionary
+    passed to ``cache`` will be modified during execution to contain any new matrices and the
+    tape offset.
+
+    >>> cache = {'matrices': [-np.eye(3)]}
+    >>> print(tape_text(tape, show_matrices=True, cache=cache))
+    0: ──U(M1)─┤ ╭<𝓗(M2)>
+    1: ──U(M1)─┤ ╰<𝓗(M2)>
+    M0 =
+    [[-1. -0. -0.]
+    [-0. -1. -0.]
+    [-0. -0. -1.]]
+    M1 =
+    [[1. 0.]
+    [0. 1.]]
+    M2 =
+    [[1. 0. 0. 0.]
+    [0. 1. 0. 0.]
+    [0. 0. 1. 0.]
+    [0. 0. 0. 1.]]
+    >>> cache
+    {'matrices': [tensor([[-1., -0., -0.],
+        [-0., -1., -0.],
+        [-0., -0., -1.]], requires_grad=True), tensor([[1., 0.],
+        [0., 1.]], requires_grad=True), tensor([[1., 0., 0., 0.],
+        [0., 1., 0., 0.],
+        [0., 0., 1., 0.],
+        [0., 0., 0., 1.]], requires_grad=True)], 'tape_offset': 0}
+
     When the provided tape has nested tapes inside, this function is called recursively.
     To maintain numbering of tapes to arbitrary levels of nesting, the ``cache`` keyword
     uses the ``"tape_offset"`` value to determine numbering. Note that the value is updated
@@ -202,8 +264,9 @@ def tape_text(
         New tape offset:  4
 
     """
-    if cache is None:
-        cache = {"tape_offset": 0}
+    cache = cache or {}
+    cache.setdefault("tape_offset", 0)
+    cache.setdefault("matrices", [])
     tape_cache = []
 
     wire_map = convert_wire_order(
@@ -241,7 +304,7 @@ def tape_text(
                         layer_str[wire_map[w]] += label
                     tape_cache.append(op)
                 else:
-                    layer_str = add(op, layer_str, wire_map, decimals)
+                    layer_str = add(op, layer_str, wire_map, decimals, cache)
 
             max_label_len = max(len(s) for s in layer_str)
             layer_str = [s.ljust(max_label_len, filler) for s in layer_str]
@@ -264,7 +327,21 @@ def tape_text(
     cache["tape_offset"] += len(tape_cache)
     for i, nested_tape in enumerate(tape_cache):
         label = f"\nTape:{i+current_tape_offset}"
-        tape_str = tape_text(nested_tape, wire_order, show_all_wires, decimals, max_length, cache)
+        tape_str = tape_text(
+            nested_tape,
+            wire_order,
+            show_all_wires,
+            decimals,
+            max_length,
+            show_matrices=False,
+            cache=cache,
+        )
         tape_totals = "\n".join([tape_totals, label, tape_str])
 
+    if show_matrices:
+        mat_str = ""
+        for i, mat in enumerate(cache["matrices"]):
+            mat_str += f"\nM{i} = \n{mat}"
+        return tape_totals + mat_str
+
     return tape_totals
diff --git a/pennylane/operation.py b/pennylane/operation.py
index f8685e84b..6be69ab92 100644
--- a/pennylane/operation.py
+++ b/pennylane/operation.py
@@ -799,13 +799,15 @@ class Operator(abc.ABC):
     def name(self, value):
         self._name = value
 
-    def label(self, decimals=None, base_label=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         r"""A customizable string representation of the operator.
 
         Args:
             decimals=None (int): If ``None``, no parameters are included. Else,
                 specifies how to round the parameters.
             base_label=None (str): overwrite the non-parameter component of the label
+            cache=None (dict): dictionary that caries information between label calls
+                in the same drawing
 
         Returns:
             str: label to use in drawings
@@ -825,16 +827,58 @@ class Operator(abc.ABC):
         >>> op.label()
         "RX⁻¹"
 
+        If the operation has a matrix-valued parameter and a cache dictionary is provided,
+        unique matrices will be cached in the ``'matrices'`` key list. The label will contain
+        the index of the matrix in the ``'matrices'`` list.
+
+        >>> op2 = qml.QubitUnitary(np.eye(2), wires=0)
+        >>> cache = {'matrices': []}
+        >>> op2.label(cache=cache)
+        'U(M0)'
+        >>> cache['matrices']
+        [tensor([[1., 0.],
+         [0., 1.]], requires_grad=True)]
+        >>> op3 = qml.QubitUnitary(np.eye(4), wires=(0,1))
+        >>> op3.label(cache=cache)
+        'U(M1)'
+        >>> cache['matrices']
+        [tensor([[1., 0.],
+                [0., 1.]], requires_grad=True),
+        tensor([[1., 0., 0., 0.],
+                [0., 1., 0., 0.],
+                [0., 0., 1., 0.],
+                [0., 0., 0., 1.]], requires_grad=True)]
+
         """
         op_label = base_label or self.__class__.__name__
 
-        if decimals is None or self.num_params == 0:
+        if self.num_params == 0:
             return op_label
 
         params = self.parameters
 
-        # matrix parameters not rendered
         if len(qml.math.shape(params[0])) != 0:
+            # assume that if the first parameter is matrix-valued, there is only a single parameter
+            # this holds true for all current operations and templates
+            if (
+                cache is None
+                or not isinstance(cache.get("matrices", None), list)
+                or len(params) != 1
+            ):
+                return op_label
+
+            for i, mat in enumerate(cache["matrices"]):
+                if qml.math.shape(params[0]) == qml.math.shape(mat) and qml.math.allclose(
+                    params[0], mat
+                ):
+                    return f"{op_label}(M{i})"
+
+            # matrix not in cache
+            mat_num = len(cache["matrices"])
+            cache["matrices"].append(params[0])
+            return f"{op_label}(M{mat_num})"
+
+        if decimals is None:
             return op_label
 
         def _format(x):
@@ -844,9 +888,6 @@ class Operator(abc.ABC):
                 # If the parameter can't be displayed as a float
                 return format(x)
 
-        if self.num_params == 1:
-            return op_label + f"\n({_format(params[0])})"
-
         param_string = ",\n".join(_format(p) for p in params)
         return op_label + f"\n({param_string})"
 
@@ -1342,11 +1383,11 @@ class Operation(Operator):
         """Name of the operator."""
         return self._name + ".inv" if self.inverse else self._name
 
-    def label(self, decimals=None, base_label=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         if self.inverse:
             base_label = base_label or self.__class__.__name__
             base_label += "⁻¹"
-        return super().label(decimals=decimals, base_label=base_label)
+        return super().label(decimals=decimals, base_label=base_label, cache=cache)
 
     def __init__(self, *params, wires=None, do_queue=True, id=None):
 
@@ -1605,7 +1646,7 @@ class Tensor(Observable):
         self._args = args
         self.queue(init=True)
 
-    def label(self, decimals=None, base_label=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         r"""How the operator is represented in diagrams and drawings.
 
         Args:
@@ -1613,6 +1654,8 @@ class Tensor(Observable):
                 how to round the parameters.
             base_label=None (Iterable[str]): overwrite the non-parameter component of the label.
                 Must be same length as ``obs`` attribute.
+            cache=None (dict): dictionary that caries information between label calls
+                in the same drawing
 
         Returns:
             str: label to use in drawings
diff --git a/pennylane/ops/cv.py b/pennylane/ops/cv.py
index ad1a3271a..55eb663b0 100644
--- a/pennylane/ops/cv.py
+++ b/pennylane/ops/cv.py
@@ -115,8 +115,8 @@ class Rotation(CVOperation):
     def adjoint(self, do_queue=True):
         return Rotation(-self.parameters[0], wires=self.wires, do_queue=do_queue)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "R")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "R", cache=cache)
 
 
 class Squeezing(CVOperation):
@@ -173,8 +173,8 @@ class Squeezing(CVOperation):
         new_phi = (phi + np.pi) % (2 * np.pi)
         return Squeezing(r, new_phi, wires=self.wires, do_queue=do_queue)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "S")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "S", cache=cache)
 
 
 class Displacement(CVOperation):
@@ -232,8 +232,8 @@ class Displacement(CVOperation):
         new_phi = (phi + np.pi) % (2 * np.pi)
         return Displacement(a, new_phi, wires=self.wires, do_queue=do_queue)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "D")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "D", cache=cache)
 
 
 class Beamsplitter(CVOperation):
@@ -294,8 +294,8 @@ class Beamsplitter(CVOperation):
         theta, phi = self.parameters
         return Beamsplitter(-theta, phi, wires=self.wires, do_queue=do_queue)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "BS")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "BS", cache=cache)
 
 
 class TwoModeSqueezing(CVOperation):
@@ -364,8 +364,8 @@ class TwoModeSqueezing(CVOperation):
         new_phi = (phi + np.pi) % (2 * np.pi)
         return TwoModeSqueezing(r, new_phi, wires=self.wires, do_queue=do_queue)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "S")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "S", cache=cache)
 
 
 class QuadraticPhase(CVOperation):
@@ -417,8 +417,8 @@ class QuadraticPhase(CVOperation):
         U[2, 1] = p[0]
         return U
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "P")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "P", cache=cache)
 
 
 class ControlledAddition(CVOperation):
@@ -476,8 +476,8 @@ class ControlledAddition(CVOperation):
     def adjoint(self, do_queue=True):
         return ControlledAddition(-self.parameters[0], wires=self.wires, do_queue=do_queue)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "X")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "X", cache=cache)
 
 
 class ControlledPhase(CVOperation):
@@ -535,8 +535,8 @@ class ControlledPhase(CVOperation):
     def adjoint(self, do_queue=True):
         return ControlledPhase(-self.parameters[0], wires=self.wires, do_queue=do_queue)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "Z")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "Z", cache=cache)
 
 
 class Kerr(CVOperation):
@@ -631,8 +631,8 @@ class CubicPhase(CVOperation):
     def adjoint(self, do_queue=True):
         return CubicPhase(-self.parameters[0], wires=self.wires, do_queue=do_queue)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "V")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "V", cache=cache)
 
 
 class InterferometerUnitary(CVOperation):
@@ -699,8 +699,8 @@ class InterferometerUnitary(CVOperation):
             qml_math.T(qml_math.conj(U)), wires=self.wires, do_queue=do_queue
         )
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "U")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "U", cache=cache)
 
 
 # =============================================================================
@@ -822,8 +822,8 @@ class ThermalState(CVOperation):
     def __init__(self, nbar, wires, do_queue=True, id=None):
         super().__init__(nbar, wires=wires, do_queue=do_queue, id=id)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "Thermal")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "Thermal", cache=cache)
 
 
 class GaussianState(CVOperation):
@@ -852,8 +852,8 @@ class GaussianState(CVOperation):
     def __init__(self, V, r, wires, do_queue=True, id=None):
         super().__init__(V, r, wires=wires, do_queue=do_queue, id=id)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "Gaussian")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "Gaussian", cache=cache)
 
 
 class FockState(CVOperation):
@@ -880,13 +880,15 @@ class FockState(CVOperation):
     def __init__(self, n, wires, do_queue=True, id=None):
         super().__init__(n, wires=wires, do_queue=do_queue, id=id)
 
-    def label(self, decimals=None, base_label=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         r"""A customizable string representation of the operator.
 
         Args:
             decimals=None (int): If ``None``, no parameters are included. Else,
                 specifies how to round the parameters.
             base_label=None (str): overwrite the non-parameter component of the label
+            cache=None (dict): dictionary that caries information between label calls
+                in the same drawing
 
         Returns:
             str: label to use in drawings
@@ -970,13 +972,15 @@ class FockStateVector(CVOperation):
     def __init__(self, state, wires, do_queue=True, id=None):
         super().__init__(state, wires=wires, do_queue=do_queue, id=id)
 
-    def label(self, decimals=None, base_label=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         r"""A customizable string representation of the operator.
 
         Args:
             decimals=None (int): If ``None``, no parameters are included. Else,
                 specifies how to round the parameters.
             base_label=None (str): overwrite the non-parameter component of the label
+            cache=None (dict): dictionary that caries information between label calls
+                in the same drawing
 
         Returns:
             str: label to use in drawings
@@ -1100,7 +1104,7 @@ class NumberOperator(CVObservable):
         hbar = 2
         return np.diag([-0.5, 0.5 / hbar, 0.5 / hbar])
 
-    def label(self, decimals=None, base_label=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         return base_label or "n"
 
 
@@ -1157,7 +1161,7 @@ class TensorN(CVObservable):
 
         return super().__new__(cls)
 
-    def label(self, decimals=None, base_label=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         if base_label is not None:
             return base_label
         return "⊗".join("n" for _ in self.wires)
@@ -1268,13 +1272,15 @@ class QuadOperator(CVObservable):
         phi = p[0]
         return np.array([0, math.cos(phi), math.sin(phi)])  # TODO check
 
-    def label(self, decimals=None, base_label=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         r"""A customizable string representation of the operator.
 
         Args:
             decimals=None (int): If ``None``, no parameters are included. Else,
                 specifies how to round the parameters.
             base_label=None (str): overwrite the non-parameter component of the label
+            cache=None (dict): dictionary that caries information between label calls
+                in the same drawing
 
         Returns:
             str: label to use in drawings
@@ -1292,7 +1298,7 @@ class QuadOperator(CVObservable):
         """
 
         if base_label is not None:
-            return super().label(decimals=decimals, base_label=base_label)
+            return super().label(decimals=decimals, base_label=base_label, cache=cache)
 
         if decimals is None:
             p = "φ"
@@ -1401,13 +1407,15 @@ class FockStateProjector(CVObservable):
     def __init__(self, n, wires, do_queue=True, id=None):
         super().__init__(n, wires=wires, do_queue=do_queue, id=id)
 
-    def label(self, decimals=None, base_label=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         r"""A customizable string representation of the operator.
 
         Args:
             decimals=None (int): If ``None``, no parameters are included. Else,
                 specifies how to round the parameters.
             base_label=None (str): overwrite the non-parameter component of the label
+            cache=None (dict): dictionary that caries information between label calls
+                in the same drawing
 
         Returns:
             str: label to use in drawings
@@ -1420,7 +1428,7 @@ class FockStateProjector(CVObservable):
         """
 
         if base_label is not None:
-            return super().label(decimals=decimals, base_label=base_label)
+            return super().label(decimals=decimals, base_label=base_label, cache=cache)
 
         basis_string = "".join(str(int(i)) for i in self.parameters[0])
         return f"|{basis_string}⟩⟨{basis_string}|"
diff --git a/pennylane/ops/identity.py b/pennylane/ops/identity.py
index 81ffea63e..9b5d8b436 100644
--- a/pennylane/ops/identity.py
+++ b/pennylane/ops/identity.py
@@ -41,7 +41,7 @@ class Identity(CVObservable, Operation):
 
     ev_order = 1
 
-    def label(self, decimals=None, base_label=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         return base_label or "I"
 
     @staticmethod
diff --git a/pennylane/ops/qubit/arithmetic_ops.py b/pennylane/ops/qubit/arithmetic_ops.py
index 45e672d9e..b82842415 100644
--- a/pennylane/ops/qubit/arithmetic_ops.py
+++ b/pennylane/ops/qubit/arithmetic_ops.py
@@ -237,8 +237,8 @@ class QubitSum(Operation):
     num_params = 0
     """int: Number of trainable parameters that the operator depends on."""
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "Σ")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "Σ", cache=cache)
 
     @staticmethod
     def compute_matrix():  # pylint: disable=arguments-differ
diff --git a/pennylane/ops/qubit/hamiltonian.py b/pennylane/ops/qubit/hamiltonian.py
index 777e8ef05..986b94c11 100644
--- a/pennylane/ops/qubit/hamiltonian.py
+++ b/pennylane/ops/qubit/hamiltonian.py
@@ -215,9 +215,9 @@ class Hamiltonian(Observable):
         # while H.coeffs is the original tensor
         super().__init__(*coeffs_flat, wires=self._wires, id=id, do_queue=do_queue)
 
-    def label(self, decimals=None, base_label=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         decimals = None if (len(self.parameters) > 3) else decimals
-        return super().label(decimals=decimals, base_label=base_label or "𝓗")
+        return super().label(decimals=decimals, base_label=base_label or "𝓗", cache=cache)
 
     @property
     def coeffs(self):
diff --git a/pennylane/ops/qubit/matrix_ops.py b/pennylane/ops/qubit/matrix_ops.py
index f232fb92d..1b62ca92c 100644
--- a/pennylane/ops/qubit/matrix_ops.py
+++ b/pennylane/ops/qubit/matrix_ops.py
@@ -156,8 +156,8 @@ class QubitUnitary(Operation):
     def _controlled(self, wire):
         ControlledQubitUnitary(*self.parameters, control_wires=wire, wires=self.wires)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "U")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "U", cache=cache)
 
 
 class ControlledQubitUnitary(QubitUnitary):
@@ -439,5 +439,5 @@ class DiagonalQubitUnitary(Operation):
             wires=Wires(control) + self.wires,
         )
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "U")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "U", cache=cache)
diff --git a/pennylane/ops/qubit/non_parametric_ops.py b/pennylane/ops/qubit/non_parametric_ops.py
index f6eb6db7b..357634eb3 100644
--- a/pennylane/ops/qubit/non_parametric_ops.py
+++ b/pennylane/ops/qubit/non_parametric_ops.py
@@ -49,7 +49,7 @@ class Hadamard(Observable, Operation):
     num_params = 0
     """int: Number of trainable parameters that the operator depends on."""
 
-    def label(self, decimals=None, base_label=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         return base_label or "H"
 
     @staticmethod
@@ -181,7 +181,7 @@ class PauliX(Observable, Operation):
 
     basis = "X"
 
-    def label(self, decimals=None, base_label=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         return base_label or "X"
 
     @staticmethod
@@ -318,7 +318,7 @@ class PauliY(Observable, Operation):
 
     basis = "Y"
 
-    def label(self, decimals=None, base_label=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         return base_label or "Y"
 
     @staticmethod
@@ -455,7 +455,7 @@ class PauliZ(Observable, Operation):
 
     basis = "Z"
 
-    def label(self, decimals=None, base_label=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         return base_label or "Z"
 
     @staticmethod
@@ -902,7 +902,7 @@ class CNOT(Operation):
 
     basis = "X"
 
-    def label(self, decimals=None, base_label=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         return base_label or "X"
 
     @staticmethod
@@ -966,7 +966,7 @@ class CZ(Operation):
 
     basis = "Z"
 
-    def label(self, decimals=None, base_label=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         return base_label or "Z"
 
     @staticmethod
@@ -1052,7 +1052,7 @@ class CY(Operation):
 
     basis = "Y"
 
-    def label(self, decimals=None, base_label=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         return base_label or "Y"
 
     @staticmethod
@@ -1478,7 +1478,7 @@ class CSWAP(Operation):
     num_params = 0
     """int: Number of trainable parameters that the operator depends on."""
 
-    def label(self, decimals=None, base_label=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         return base_label or "SWAP"
 
     @staticmethod
@@ -1586,7 +1586,7 @@ class Toffoli(Operation):
 
     basis = "X"
 
-    def label(self, decimals=None, base_label=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         return base_label or "X"
 
     @staticmethod
@@ -1790,6 +1790,9 @@ class MultiControlledX(Operation):
 
         super().__init__(wires=total_wires, do_queue=do_queue)
 
+    def label(self, decimals=None, base_label=None, cache=None):
+        return base_label or "X"
+
     # pylint: disable=unused-argument
     @staticmethod
     def compute_matrix(
@@ -1847,9 +1850,6 @@ class MultiControlledX(Operation):
     def control_wires(self):
         return self.wires[:~0]
 
-    def label(self, decimals=None, base_label=None):
-        return base_label or "X"
-
     def adjoint(self):
         return MultiControlledX(
             control_wires=self.wires[:-1],
@@ -2050,7 +2050,7 @@ class Barrier(Operation):
         """
         return []
 
-    def label(self, decimals=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         return "||"
 
     def adjoint(self):
@@ -2092,7 +2092,7 @@ class WireCut(Operation):
         """
         return []
 
-    def label(self, decimals=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         return "//"
 
     def adjoint(self):
diff --git a/pennylane/ops/qubit/observables.py b/pennylane/ops/qubit/observables.py
index 033619d2b..5f2e2571e 100644
--- a/pennylane/ops/qubit/observables.py
+++ b/pennylane/ops/qubit/observables.py
@@ -63,8 +63,8 @@ class Hermitian(Observable):
     def __init__(self, A, wires, do_queue=True, id=None):
         super().__init__(A, wires=wires, do_queue=do_queue, id=id)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "𝓗")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "𝓗", cache=cache)
 
     @staticmethod
     def compute_matrix(A):  # pylint: disable=arguments-differ
@@ -209,8 +209,8 @@ class SparseHamiltonian(Observable):
             raise TypeError("Observable must be a scipy sparse coo_matrix.")
         super().__init__(H, wires=wires, do_queue=do_queue, id=id)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "𝓗")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "𝓗", cache=cache)
 
     @staticmethod
     def compute_matrix(H):  # pylint: disable=arguments-differ
@@ -330,13 +330,15 @@ class Projector(Observable):
 
         super().__init__(basis_state, wires=wires, do_queue=do_queue, id=id)
 
-    def label(self, decimals=None, base_label=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         r"""A customizable string representation of the operator.
 
         Args:
             decimals=None (int): If ``None``, no parameters are included. Else,
                 specifies how to round the parameters.
             base_label=None (str): overwrite the non-parameter component of the label
+            cache=None (dict): dictionary that caries information between label calls
+                in the same drawing
 
         Returns:
             str: label to use in drawings
diff --git a/pennylane/ops/qubit/parametric_ops.py b/pennylane/ops/qubit/parametric_ops.py
index 37bca18ab..9ddad2dff 100644
--- a/pennylane/ops/qubit/parametric_ops.py
+++ b/pennylane/ops/qubit/parametric_ops.py
@@ -341,8 +341,8 @@ class PhaseShift(Operation):
     def __init__(self, phi, wires, do_queue=True, id=None):
         super().__init__(phi, wires=wires, do_queue=do_queue, id=id)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "Rϕ")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "Rϕ", cache=cache)
 
     @staticmethod
     def compute_matrix(phi):  # pylint: disable=arguments-differ
@@ -483,8 +483,8 @@ class ControlledPhaseShift(Operation):
     def __init__(self, phi, wires, do_queue=True, id=None):
         super().__init__(phi, wires=wires, do_queue=do_queue, id=id)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "Rϕ")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "Rϕ", cache=cache)
 
     @staticmethod
     def compute_matrix(phi):  # pylint: disable=arguments-differ
@@ -952,13 +952,15 @@ class PauliRot(Operation):
                 f"The given Pauli word has length {len(pauli_word)}, length {num_wires} was expected for wires {wires}"
             )
 
-    def label(self, decimals=None, base_label=None):
+    def label(self, decimals=None, base_label=None, cache=None):
         r"""A customizable string representation of the operator.
 
         Args:
             decimals=None (int): If ``None``, no parameters are included. Else,
                 specifies how to round the parameters.
             base_label=None (str): overwrite the non-parameter component of the label
+            cache=None (dict): dictionary that caries information between label calls
+                in the same drawing
 
         Returns:
             str: label to use in drawings
@@ -1225,8 +1227,8 @@ class CRX(Operation):
     def __init__(self, phi, wires, do_queue=True, id=None):
         super().__init__(phi, wires=wires, do_queue=do_queue, id=id)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "RX")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "RX", cache=cache)
 
     @staticmethod
     def compute_matrix(theta):  # pylint: disable=arguments-differ
@@ -1375,8 +1377,8 @@ class CRY(Operation):
     def __init__(self, phi, wires, do_queue=True, id=None):
         super().__init__(phi, wires=wires, do_queue=do_queue, id=id)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "RY")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "RY", cache=cache)
 
     @staticmethod
     def compute_matrix(theta):  # pylint: disable=arguments-differ
@@ -1518,8 +1520,8 @@ class CRZ(Operation):
     def __init__(self, phi, wires, do_queue=True, id=None):
         super().__init__(phi, wires=wires, do_queue=do_queue, id=id)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "RZ")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "RZ", cache=cache)
 
     @staticmethod
     def compute_matrix(theta):  # pylint: disable=arguments-differ
@@ -1679,8 +1681,8 @@ class CRot(Operation):
     def __init__(self, phi, theta, omega, wires, do_queue=True, id=None):
         super().__init__(phi, theta, omega, wires=wires, do_queue=do_queue, id=id)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "Rot")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "Rot", cache=cache)
 
     @staticmethod
     def compute_matrix(phi, theta, omega):  # pylint: disable=arguments-differ
diff --git a/pennylane/ops/qubit/qchem_ops.py b/pennylane/ops/qubit/qchem_ops.py
index 4fcc332c3..882249d15 100644
--- a/pennylane/ops/qubit/qchem_ops.py
+++ b/pennylane/ops/qubit/qchem_ops.py
@@ -166,8 +166,8 @@ class SingleExcitation(Operation):
         (phi,) = self.parameters
         return SingleExcitation(-phi, wires=self.wires)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "G")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "G", cache=cache)
 
 
 class SingleExcitationMinus(Operation):
@@ -305,8 +305,8 @@ class SingleExcitationMinus(Operation):
         (phi,) = self.parameters
         return SingleExcitationMinus(-phi, wires=self.wires)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "G₋")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "G₋", cache=cache)
 
 
 class SingleExcitationPlus(Operation):
@@ -444,8 +444,8 @@ class SingleExcitationPlus(Operation):
         (phi,) = self.parameters
         return SingleExcitationPlus(-phi, wires=self.wires)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "G₊")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "G₊", cache=cache)
 
 
 class DoubleExcitation(Operation):
@@ -644,8 +644,8 @@ class DoubleExcitation(Operation):
         (theta,) = self.parameters
         return DoubleExcitation(-theta, wires=self.wires)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "G²")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "G²", cache=cache)
 
 
 class DoubleExcitationPlus(Operation):
@@ -739,8 +739,8 @@ class DoubleExcitationPlus(Operation):
         (theta,) = self.parameters
         return DoubleExcitationPlus(-theta, wires=self.wires)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "G²₊")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "G²₊", cache=cache)
 
 
 class DoubleExcitationMinus(Operation):
@@ -832,8 +832,8 @@ class DoubleExcitationMinus(Operation):
         (theta,) = self.parameters
         return DoubleExcitationMinus(-theta, wires=self.wires)
 
-    def label(self, decimals=None, base_label=None):
-        return super().label(decimals=decimals, base_label=base_label or "G²₋")
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label=base_label or "G²₋", cache=cache)
 
 
 class OrbitalRotation(Operation):
diff --git a/pennylane/transforms/draw.py b/pennylane/transforms/draw.py
index 72e985946..439cfca02 100644
--- a/pennylane/transforms/draw.py
+++ b/pennylane/transforms/draw.py
@@ -158,6 +158,7 @@ def draw(
     show_all_wires=False,
     decimals=2,
     max_length=100,
+    show_matrices=False,
     expansion_strategy=None,
 ):
     """Create a function that draws the given qnode.
@@ -169,6 +170,7 @@ def draw(
         decimals (int): How many decimal points to include when formatting operation parameters.
             ``None`` will omit parameters from operation labels.
         max_length (int): Maximum string width (columns) when printing the circuit
+        show_matrices=False (bool): show matrix valued parameters below all circuit diagrams
         expansion_strategy (str): The strategy to use when circuit expansions or decompositions
             are required.
 
@@ -203,8 +205,7 @@ def draw(
     .. UsageDetails::
 
 
-    By specifying the ``decimals`` keyword, parameters
-    are displayed to the specified precision. Matrix-valued parameters are never displayed.
+    By specifying the ``decimals`` keyword, parameters are displayed to the specified precision.
 
     >>> print(qml.draw(circuit, decimals=4)(a=2.3, w=[1.2, 3.2, 0.7]))
     0: ──H─╭C─────────────────────────────────────╭C───────────┤ ╭<Z@Z>
@@ -226,6 +227,28 @@ def draw(
     >>> print(qml.draw(circuit2)("x"))
     0: ──RX(x)─┤  <Z>
 
+    When requested with ``show_matrices=True``, matrix valued parameters are printed below the
+    circuit:
+
+    .. code-block:: python
+
+    >>> @qml.qnode(qml.device('default.qubit', wires=2))
+    ... def circuit():
+    ...     qml.QubitUnitary(np.eye(2), wires=0)
+    ...     qml.QubitUnitary(-np.eye(4), wires=(0,1))
+    ...     return qml.expval(qml.Hermitian(np.eye(2), wires=1))
+    >>> print(qml.draw(circuit, show_matrices=True)())
+    0: ──U(M0)─╭U(M1)─┤
+    1: ────────╰U(M1)─┤  <𝓗(M0)>
+    M0 =
+    [[1. 0.]
+    [0. 1.]]
+    M1 =
+    [[-1. -0. -0. -0.]
+    [-0. -1. -0. -0.]
+    [-0. -0. -1. -0.]
+    [-0. -0. -0. -1.]]
+
     The ``max_length`` keyword warps long circuits:
 
     .. code-block:: python
@@ -302,18 +325,24 @@ def draw(
         _wire_order = wire_order or qnode.device.wires
 
         if tapes is not None:
-            cache = {"tape_offset": 0}
+            cache = {"tape_offset": 0, "matrices": []}
             res = [
                 tape_text(
                     t,
                     wire_order=_wire_order,
                     show_all_wires=show_all_wires,
                     decimals=decimals,
+                    show_matrices=False,
                     max_length=max_length,
                     cache=cache,
                 )
                 for t in tapes[0]
             ]
+            if show_matrices:
+                mat_str = "\n"
+                for i, mat in enumerate(cache["matrices"]):
+                    mat_str += f"\nM{i} = \n{mat}"
+                return "\n\n".join(res) + mat_str
             return "\n\n".join(res)
 
         return tape_text(
@@ -321,6 +350,7 @@ def draw(
             wire_order=_wire_order,
             show_all_wires=show_all_wires,
             decimals=decimals,
+            show_matrices=show_matrices,
             max_length=max_length,
         )
 


Test patch:
diff --git a/tests/drawer/test_tape_text.py b/tests/drawer/test_tape_text.py
index a8a2be0ca..efaab56bd 100644
--- a/tests/drawer/test_tape_text.py
+++ b/tests/drawer/test_tape_text.py
@@ -17,6 +17,7 @@ Unit tests for the `pennylane.draw_text` function.
 
 import pytest
 import pennylane as qml
+from pennylane import numpy as np
 
 from pennylane.drawer import tape_text
 from pennylane.drawer.tape_text import _add_grouping_symbols, _add_op, _add_measurement
@@ -54,7 +55,19 @@ class TestHelperFunctions:
     )
     def test_add_measurements(self, op, out):
         """Test private _add_measurement function renders as expected."""
-        assert out == _add_measurement(op, [""] * 4, default_wire_map, None)
+        assert out == _add_measurement(op, [""] * 4, default_wire_map, None, None)
+
+    def test_add_measurements_cache(self):
+        """Test private _add_measurement function with a matrix cache."""
+        cache = {"matrices": []}
+        op = qml.expval(qml.Hermitian(np.eye(2), wires=0))
+        assert _add_measurement(op, ["", ""], {0: 0, 1: 1}, None, cache) == ["<𝓗(M0)>", ""]
+
+        assert qml.math.allclose(cache["matrices"][0], np.eye(2))
+
+        op2 = qml.expval(qml.Hermitian(np.eye(2), wires=1))
+        # new op with same matrix, should have same M0 designation
+        assert _add_measurement(op2, ["", ""], {0: 0, 1: 1}, None, cache) == ["", "<𝓗(M0)>"]
 
     @pytest.mark.parametrize(
         "op, out",
@@ -67,7 +80,7 @@ class TestHelperFunctions:
     )
     def test_add_op(self, op, out):
         """Test adding the first operation to array of strings"""
-        assert out == _add_op(op, ["─"] * 4, default_wire_map, None)
+        assert out == _add_op(op, ["─"] * 4, default_wire_map, None, None)
 
     @pytest.mark.parametrize(
         "op, out",
@@ -79,8 +92,18 @@ class TestHelperFunctions:
     )
     def test_add_second_op(self, op, out):
         """Test adding a second operation to the array of strings"""
-        start = _add_op(qml.PauliX(0), ["─"] * 4, default_wire_map, None)
-        assert out == _add_op(op, start, default_wire_map, None)
+        start = _add_op(qml.PauliX(0), ["─"] * 4, default_wire_map, None, None)
+        assert out == _add_op(op, start, default_wire_map, None, None)
+
+    def test_add_op_cache(self):
+        """Test private _add_op method functions with a matrix cache."""
+        cache = {"matrices": []}
+        op1 = qml.QubitUnitary(np.eye(2), wires=0)
+        assert _add_op(op1, ["", ""], {0: 0, 1: 1}, None, cache) == ["U(M0)", ""]
+
+        assert qml.math.allclose(cache["matrices"][0], np.eye(2))
+        op2 = qml.QubitUnitary(np.eye(2), wires=1)
+        assert _add_op(op2, ["", ""], {0: 0, 1: 1}, None, cache) == ["", "U(M0)"]
 
 
 class TestEmptyTapes:
@@ -220,12 +243,12 @@ single_op_tests_data = [
         qml.DoubleExcitationPlus(1.23, wires=(0, 1, 2, 3)),
         "0: ─╭G²₊(1.23)─┤  \n1: ─├G²₊(1.23)─┤  \n2: ─├G²₊(1.23)─┤  \n3: ─╰G²₊(1.23)─┤  ",
     ),
-    (qml.QubitUnitary(qml.numpy.eye(4), wires=(0, 1)), "0: ─╭U─┤  \n1: ─╰U─┤  "),
+    (qml.QubitUnitary(qml.numpy.eye(4), wires=(0, 1)), "0: ─╭U(M0)─┤  \n1: ─╰U(M0)─┤  "),
     (qml.QubitSum(wires=(0, 1, 2)), "0: ─╭Σ─┤  \n1: ─├Σ─┤  \n2: ─╰Σ─┤  "),
     (qml.AmplitudeDamping(0.98, wires=0), "0: ──AmplitudeDamping(0.98)─┤  "),
     (
         qml.QubitStateVector([0, 1, 0, 0], wires=(0, 1)),
-        "0: ─╭QubitStateVector─┤  \n1: ─╰QubitStateVector─┤  ",
+        "0: ─╭QubitStateVector(M0)─┤  \n1: ─╰QubitStateVector(M0)─┤  ",
     ),
     (qml.Kerr(1.234, wires=0), "0: ──Kerr(1.23)─┤  "),
     (
@@ -295,6 +318,55 @@ class TestLayering:
         assert tape_text(tape, wire_order=[0, 1, 2]) == expected
 
 
+with qml.tape.QuantumTape() as tape_matrices:
+    qml.QubitStateVector([1.0, 0.0], wires=(0, 1))
+    qml.QubitUnitary(np.eye(2), wires=0)
+    qml.expval(qml.Hermitian(np.eye(2), wires=0))
+
+
+class TestShowMatrices:
+    """Test the handling of matrix-valued parameters."""
+
+    def test_default_no_matrix_parameters(self):
+        """Test matrices numbered but not included by default."""
+
+        expected = (
+            "0: ─╭QubitStateVector(M0)──U(M1)─┤  <𝓗(M1)>\n"
+            "1: ─╰QubitStateVector(M0)────────┤         "
+        )
+
+        assert tape_text(tape_matrices) == expected
+
+    def test_show_matrices(self):
+        """Test matrices included when requested."""
+
+        expected = (
+            "0: ─╭QubitStateVector(M0)──U(M1)─┤  <𝓗(M1)>\n"
+            "1: ─╰QubitStateVector(M0)────────┤         \n"
+            "M0 = \n[1.0, 0.0]\n"
+            "M1 = \n[[1. 0.]\n [0. 1.]]"
+        )
+
+        assert tape_text(tape_matrices, show_matrices=True) == expected
+
+    def test_matrix_parameters_provided_cache(self):
+        """Providing an existing matrix cache determines numbering order of matrices.
+        All matrices printed out regardless of use."""
+
+        cache = {"matrices": [np.eye(2), -np.eye(3)]}
+
+        expected = (
+            "0: ─╭QubitStateVector(M2)──U(M0)─┤  <𝓗(M0)>\n"
+            "1: ─╰QubitStateVector(M2)────────┤         \n"
+            "M0 = \n[[1. 0.]\n [0. 1.]]\n"
+            "M1 = \n[[-1. -0. -0.]\n [-0. -1. -0.]\n [-0. -0. -1.]]\n"
+            "M2 = \n[1.0, 0.0]"
+        )
+
+        assert tape_text(tape_matrices, show_matrices=True, cache=cache) == expected
+        assert cache["matrices"][2] == [1.0, 0.0]
+
+
 class TestNestedTapes:
     """Test situations with nested tapes."""
 
diff --git a/tests/ops/qubit/test_matrix_ops.py b/tests/ops/qubit/test_matrix_ops.py
index 6ce360fb2..66e00403b 100644
--- a/tests/ops/qubit/test_matrix_ops.py
+++ b/tests/ops/qubit/test_matrix_ops.py
@@ -444,18 +444,90 @@ class TestControlledQubitUnitary:
 
 
 label_data = [
-    (qml.QubitUnitary(X, wires=0), "U"),
-    (qml.DiagonalQubitUnitary([1, 1], wires=1), "U"),
-    (qml.ControlledQubitUnitary(X, control_wires=0, wires=1), "U"),
+    (X, qml.QubitUnitary(X, wires=0)),
+    (X, qml.ControlledQubitUnitary(X, control_wires=0, wires=1)),
+    ([1, 1], qml.DiagonalQubitUnitary([1, 1], wires=0)),
 ]
 
 
-@pytest.mark.parametrize("op, label", label_data)
-def test_label(op, label):
-    assert op.label() == label
-    assert op.label(decimals=5) == label
-    op.inv()
-    assert op.label() == label + "⁻¹"
+@pytest.mark.parametrize("mat, op", label_data)
+class TestUnitaryLabels:
+    def test_no_cache(self, mat, op):
+        """Test labels work without a provided cache."""
+        assert op.label() == "U"
+
+    def test_matrices_not_in_cache(self, mat, op):
+        """Test provided cache doesn't have a 'matrices' keyword."""
+        assert op.label(cache={}) == "U"
+
+    def test_cache_matrices_not_list(self, mat, op):
+        """Test 'matrices' key pair is not a list."""
+        assert op.label(cache={"matrices": 0}) == "U"
+
+    def test_empty_cache_list(self, mat, op):
+        """Test matrices list is provided, but empty. Operation should have `0` label and matrix
+        should be added to cache."""
+        cache = {"matrices": []}
+        assert op.label(cache=cache) == "U(M0)"
+        assert qml.math.allclose(cache["matrices"][0], mat)
+
+    def test_something_in_cache_list(self, mat, op):
+        """If something exists in the matrix list, but parameter is not in the list, then parameter
+        added to list and label given number of its position."""
+        cache = {"matrices": [Z]}
+        assert op.label(cache=cache) == "U(M1)"
+
+        assert len(cache["matrices"]) == 2
+        assert qml.math.allclose(cache["matrices"][1], mat)
+
+    def test_matrix_already_in_cache_list(self, mat, op):
+        """If the parameter already exists in the matrix cache, then the label uses that index and the
+        matrix cache is unchanged."""
+        cache = {"matrices": [Z, mat, S]}
+        assert op.label(cache=cache) == "U(M1)"
+
+        assert len(cache["matrices"]) == 3
+
+
+class TestInterfaceMatricesLabel:
+    """Test different interface matrices with qubit."""
+
+    def check_interface(self, mat):
+        """Interface independent helper method."""
+
+        op = qml.QubitUnitary(mat, wires=0)
+
+        cache = {"matrices": []}
+        assert op.label(cache=cache) == "U(M0)"
+        assert qml.math.allclose(cache["matrices"][0], mat)
+
+        cache = {"matrices": [0, mat, 0]}
+        assert op.label(cache=cache) == "U(M1)"
+        assert len(cache["matrices"]) == 3
+
+    def test_labelling_torch_tensor(self):
+        """Test matrix cache labelling with torch interface."""
+
+        torch = pytest.importorskip("torch")
+
+        mat = torch.tensor([[1, 0], [0, -1]])
+        self.check_interface(mat)
+
+    def test_labelling_tf_variable(self):
+        """Test matrix cache labelling with tf interface."""
+
+        tf = pytest.importorskip("tensorflow")
+        mat = tf.Variable([[1, 0], [0, -1]])
+
+        self.check_interface(mat)
+
+    def test_labelling_jax_variable(self):
+        """Test matrix cache labelling with jax interface."""
+
+        jnp = pytest.importorskip("jax.numpy")
+        mat = jnp.array([[1, 0], [0, -1]])
+
+        self.check_interface(mat)
 
 
 control_data = [
diff --git a/tests/ops/qubit/test_observables.py b/tests/ops/qubit/test_observables.py
index 30484b523..364e553d6 100644
--- a/tests/ops/qubit/test_observables.py
+++ b/tests/ops/qubit/test_observables.py
@@ -480,6 +480,21 @@ label_data = [
 
 @pytest.mark.parametrize("op, label", label_data)
 def test_label_method(op, label):
+    """Test non-cache label functionality."""
     assert op.label() == label
     assert op.label(decimals=5) == label
     assert op.label(base_label="obs") == "obs"
+
+
+def test_hermitian_labelling_w_cache():
+    """Test hermitian matrix interacts with matrix cache provided to label."""
+
+    op = qml.Hermitian(X, wires=0)
+
+    cache = {"matrices": [Z]}
+    op.label(cache=cache) == "𝓗(M1)"
+    assert qml.math.allclose(cache["matrices"][1], X)
+
+    cache = {"matrices": [Z, Y, X]}
+    op.label(cache=cache) == "𝓗(M2)"
+    assert len(cache["matrices"]) == 3
diff --git a/tests/ops/qubit/test_state_prep.py b/tests/ops/qubit/test_state_prep.py
index 3fa62f62e..dcca51c39 100644
--- a/tests/ops/qubit/test_state_prep.py
+++ b/tests/ops/qubit/test_state_prep.py
@@ -21,18 +21,42 @@ from pennylane import numpy as np
 from pennylane.wires import Wires
 
 
-class TestOperations:
-    @pytest.mark.parametrize(
-        "op",
-        [
-            qml.BasisState(np.array([0, 1]), wires=0),
-            qml.QubitStateVector(np.array([1.0, 0.0]), wires=0),
-            qml.QubitDensityMatrix(np.array([[1.0, 0.0], [0.0, 0.0]]), wires=0),
-        ],
-    )
-    def test_adjoint_error_exception(self, op, tol):
-        with pytest.raises(qml.ops.AdjointError):
-            op.adjoint()
+densitymat0 = np.array([[1.0, 0.0], [0.0, 0.0]])
+
+
+@pytest.mark.parametrize(
+    "op",
+    [
+        qml.BasisState(np.array([0, 1]), wires=0),
+        qml.QubitStateVector(np.array([1.0, 0.0]), wires=0),
+        qml.QubitDensityMatrix(densitymat0, wires=0),
+    ],
+)
+def test_adjoint_error_exception(op):
+    with pytest.raises(qml.ops.AdjointError):
+        op.adjoint()
+
+
+@pytest.mark.parametrize(
+    "op, mat, base",
+    [
+        (qml.BasisState(np.array([0, 1]), wires=0), [0, 1], "BasisState"),
+        (qml.QubitStateVector(np.array([1.0, 0.0]), wires=0), [1.0, 0.0], "QubitStateVector"),
+        (qml.QubitDensityMatrix(densitymat0, wires=0), densitymat0, "QubitDensityMatrix"),
+    ],
+)
+def test_labelling_matrix_cache(op, mat, base):
+    """Test state prep matrix parameters interact with labelling matrix cache"""
+
+    assert op.label() == base
+
+    cache = {"matrices": []}
+    assert op.label(cache=cache) == base + "(M0)"
+    assert qml.math.allclose(cache["matrices"][0], mat)
+
+    cache = {"matrices": [0, mat, 0]}
+    assert op.label(cache=cache) == base + "(M1)"
+    assert len(cache["matrices"]) == 3
 
 
 class TestDecomposition:
diff --git a/tests/transforms/test_draw.py b/tests/transforms/test_draw.py
index c378bd2bb..bf66fd108 100644
--- a/tests/transforms/test_draw.py
+++ b/tests/transforms/test_draw.py
@@ -1,4 +1,4 @@
-# Copyright 2018-2020 Xanadu Quantum Technologies Inc.
+# Copyright 2018-2022 Xanadu Quantum Technologies Inc.
 
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
@@ -148,6 +148,63 @@ class TestDecimals:
         assert draw(circuit)("x", "y", "z") == expected
 
 
+class TestMatrixParameters:
+    def test_matrix_parameters(self):
+        """Test matrix valued parameters remembered and printed out upon request."""
+
+        @qml.qnode(qml.device("default.qubit", wires=2))
+        def matrices_circuit():
+            qml.QubitStateVector([1.0, 0.0], wires=(0, 1))
+            qml.QubitUnitary(np.eye(2), wires=0)
+            return qml.expval(qml.Hermitian(np.eye(2), wires=0))
+
+        expected1 = (
+            "0: ─╭QubitStateVector(M0)──U(M1)─┤  <𝓗(M1)>\n"
+            "1: ─╰QubitStateVector(M0)────────┤         "
+        )
+
+        assert draw(matrices_circuit)() == expected1
+
+        expected2 = (
+            "0: ─╭QubitStateVector(M0)──U(M1)─┤  <𝓗(M1)>\n"
+            "1: ─╰QubitStateVector(M0)────────┤         \n"
+            "M0 = \n[1.0, 0.0]\n"
+            "M1 = \n[[1. 0.]\n [0. 1.]]"
+        )
+
+        assert draw(matrices_circuit, show_matrices=True)() == expected2
+
+    def test_matrix_parameters_batch_transform(self):
+        """Test matrix parameters only printed once after a batch transform."""
+
+        @qml.gradients.param_shift(shifts=[(0.2,)])
+        @qml.qnode(qml.device("default.qubit", wires=2))
+        def matrices_circuit(x):
+            qml.QubitStateVector([1.0, 0.0], wires=(0, 1))
+            qml.QubitUnitary(np.eye(2, requires_grad=False), wires=0)
+            qml.RX(x, wires=1)
+            return qml.expval(qml.Hermitian(np.eye(2, requires_grad=False), wires=1))
+
+        expected1 = (
+            "0: ─╭QubitStateVector(M0)──U(M1)────┤         \n"
+            "1: ─╰QubitStateVector(M0)──RX(1.20)─┤  <𝓗(M1)>\n\n"
+            "0: ─╭QubitStateVector(M0)──U(M1)────┤         \n"
+            "1: ─╰QubitStateVector(M0)──RX(0.80)─┤  <𝓗(M1)>"
+        )
+
+        assert draw(matrices_circuit)(np.array(1.0, requires_grad=True))
+
+        expected2 = (
+            "0: ─╭QubitStateVector(M0)──U(M1)────┤         \n"
+            "1: ─╰QubitStateVector(M0)──RX(1.20)─┤  <𝓗(M1)>\n\n"
+            "0: ─╭QubitStateVector(M0)──U(M1)────┤         \n"
+            "1: ─╰QubitStateVector(M0)──RX(0.80)─┤  <𝓗(M1)>\n\n"
+            "M0 = \n[1.0, 0.0]\n"
+            "M1 = \n[[1. 0.]\n [0. 1.]]"
+        )
+        assert draw(matrices_circuit, show_matrices=True)(np.array(1.0, requires_grad=True))
+
+
 class TestMaxLength:
     """Test the max_length keyword."""