Repository: PennyLaneAI/pennylane
Pull request #1483 (2021-07-29 07:49:20)

Make Hamiltonian differentiable with respect to coefficients, and inherit from Observable
**Context:**

This is the second PR for the new VQE design. 

**Description of the Change:**

The Hamiltonian class now inherits from `Observable`, and all processing of the coefficients is done interface-agnostic, which makes them differentiable.

``` python
  from pennylane import numpy as np
  
  dev = qml.device("default.qubit", wires=2)
  @qml.qnode(dev)
  def circuit(coeffs, param):
    qml.RX(param, wires=0)
    qml.RY(param, wires=0)
    return qml.expval(
      qml.Hamiltonian(coeffs, [qml.PauliX(0), qml.PauliZ(0)], simplify=True)
    )
    
  coeffs = np.array([-0.05, 0.17])
  param = np.array(1.7)
  grad_fn = qml.grad(circuit)
  ```
  ``` pycon
  >>> grad_fn(coeffs, param)
  (array([-0.12777055,  0.0166009 ]), array(0.0917819))
  ```
In addition, I made `group_observables` differentiable, because if coefficients are tensors this will be important. 

  ``` python
  import jax
  from jax import numpy as jnp
  
  coeffs = jnp.array([1., 2., 3.])
  obs = [PauliX(wires=0), PauliX(wires=1), PauliZ(wires=1)]

  def group(coeffs, select=None):
    _, grouped_coeffs = qml.grouping.group_observables(obs, coeffs)
    # in this example, grouped_coeffs is a list of two jax tensors
    # [DeviceArray([1., 2.], dtype=float32), DeviceArray([3.], dtype=float32)]
    return grouped_coeffs[select]

  jac_fn = jax.jacobian(group)
  ```
  ```pycon
  >>> jac_fn(coeffs, select=0)
  [[1. 0. 0.]
  [0. 1. 0.]]
  
  >>> jac_fn(coeffs, select=1)
  [[0., 0., 1.]]
  ```

**Benefits:**

Hamiltonians can now be part of differentiation pipelines!!

**Possible Drawbacks:**

* To see the changes, the `Hamiltonian` class is left in `vqe.py`, but will be moved to `ops` in a follow-up PR.
* Hamiltonian arithmetic is _partially_ supported, for example `H1 + H2` works with tensor coefficients, but `tf.Variable(2.) * qml.PauliX(0)` does not work (since overwriting ``__mul__`` does not overwrite tensorflow's multiplication, as far as I understand)
* One open question marked by a comment is how to avoid all sorts of operations to be queued when the Hamiltonian does things like `simplify`.

*Edit: Locally all tests are passing, but for the environment with Python 3.7 on CI there seems to be no `torch.kron` function. Does that mean I cannot dispatch to kron if it wasn't present in previous torch versions?* 

codecov[bot]: # [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1483?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#1483](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1483?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (9b77ae3) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/ef494a4c7e82d0f21097a1738366636afb468303?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (ef494a4) will **decrease** coverage by `8.54%`.
> The diff coverage is `100.00%`.

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

```diff
@@            Coverage Diff             @@
##           master    #1483      +/-   ##
==========================================
- Coverage   98.22%   89.67%   -8.55%     
==========================================
  Files         184      179       -5     
  Lines       13203    12663     -540     
==========================================
- Hits        12968    11355    -1613     
- Misses        235     1308    +1073     
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1483?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Coverage Δ | |
|---|---|---|
| [pennylane/qnode.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1483/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3Fub2RlLnB5) | `81.25% <ø> (-16.84%)` | :arrow_down: |
| [pennylane/tape/tape.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1483/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RhcGUvdGFwZS5weQ==) | `99.15% <ø> (+0.54%)` | :arrow_up: |
| [pennylane/grouping/group\_observables.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1483/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2dyb3VwaW5nL2dyb3VwX29ic2VydmFibGVzLnB5) | `100.00% <100.00%> (ø)` | |
| [pennylane/math/single\_dispatch.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1483/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL21hdGgvc2luZ2xlX2Rpc3BhdGNoLnB5) | `52.35% <100.00%> (-47.12%)` | :arrow_down: |
| [pennylane/operation.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1483/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wZXJhdGlvbi5weQ==) | `99.77% <100.00%> (+3.98%)` | :arrow_up: |
| [pennylane/transforms/hamiltonian\_expand.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1483/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RyYW5zZm9ybXMvaGFtaWx0b25pYW5fZXhwYW5kLnB5) | `100.00% <100.00%> (ø)` | |
| [pennylane/vqe/vqe.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1483/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3ZxZS92cWUucHk=) | `99.49% <100.00%> (+7.25%)` | :arrow_up: |
| [pennylane/beta/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1483/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2JldGEvX19pbml0X18ucHk=) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| [pennylane/devices/tf\_ops.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1483/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2RldmljZXMvdGZfb3BzLnB5) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| [pennylane/devices/jax\_ops.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1483/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2RldmljZXMvamF4X29wcy5weQ==) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| ... and [58 more](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1483/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/1483?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/1483?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Last update [ef494a4...9b77ae3](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1483?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).

mariaschuld: @glassnotes I implemented some ideas from Josh, maybe could you try another look? 

Patch:
diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index b6e2f1029..c5bf87e98 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -2,8 +2,61 @@
 
 <h3>New features since last release</h3>
 
+* Hamiltonians are now trainable with respect to their coefficients.
+  [(#1483)](https://github.com/PennyLaneAI/pennylane/pull/1483)
+
+  ``` python
+  from pennylane import numpy as np
+  
+  dev = qml.device("default.qubit", wires=2)
+  @qml.qnode(dev)
+  def circuit(coeffs, param):
+      qml.RX(param, wires=0)
+      qml.RY(param, wires=0)
+      return qml.expval(
+          qml.Hamiltonian(coeffs, [qml.PauliX(0), qml.PauliZ(0)], simplify=True)
+      )
+    
+  coeffs = np.array([-0.05, 0.17])
+  param = np.array(1.7)
+  grad_fn = qml.grad(circuit)
+  ```
+  ``` pycon
+  >>> grad_fn(coeffs, param)
+  (array([-0.12777055,  0.0166009 ]), array(0.0917819))
+  ```
+
 <h3>Improvements</h3>
 
+* The `group_observables` transform is now differentiable.
+  [(#1483)](https://github.com/PennyLaneAI/pennylane/pull/1483)
+ 
+  For example:
+
+  ``` python
+  import jax
+  from jax import numpy as jnp
+  
+  coeffs = jnp.array([1., 2., 3.])
+  obs = [PauliX(wires=0), PauliX(wires=1), PauliZ(wires=1)]
+
+  def group(coeffs, select=None):
+    _, grouped_coeffs = qml.grouping.group_observables(obs, coeffs)
+    # in this example, grouped_coeffs is a list of two jax tensors
+    # [DeviceArray([1., 2.], dtype=float32), DeviceArray([3.], dtype=float32)]
+    return grouped_coeffs[select]
+
+  jac_fn = jax.jacobian(group)
+  ```
+  ```pycon
+  >>> jac_fn(coeffs, select=0)
+  [[1. 0. 0.]
+  [0. 1. 0.]]
+  
+  >>> jac_fn(coeffs, select=1)
+  [[0., 0., 1.]]
+  ```
+
 * The tape does not verify any more that all Observables have owners in the annotated queue.
   [(#1505)](https://github.com/PennyLaneAI/pennylane/pull/1505)
 
@@ -153,7 +206,7 @@ Maria Schuld.
   Total shots:  200
   Total shots:  300
   ```
-
+  
 * VQE problems can now intuitively been set up by passing the Hamiltonian 
   as an observable. [(#1474)](https://github.com/PennyLaneAI/pennylane/pull/1474)
 
@@ -491,6 +544,14 @@ Maria Schuld.
 
 <h3>Improvements</h3>
 
+* The tape does not verify any more that all Observables have owners in the annotated queue.
+  [(#1505)](https://github.com/PennyLaneAI/pennylane/pull/1505)
+
+  This allows manipulation of Observables inside a tape context. An example is 
+  `expval(Tensor(qml.PauliX(0), qml.Identity(1)).prune())` which makes the expval an owner 
+  of the pruned tensor and its constituent observables, but leaves the original tensor in 
+  the queue without an owner.
+
 * The `step` and `step_and_cost` methods of `QNGOptimizer` now accept a custom `grad_fn`
   keyword argument to use for gradient computations.
   [(#1487)](https://github.com/PennyLaneAI/pennylane/pull/1487)
diff --git a/doc/requirements.txt b/doc/requirements.txt
index a3fe65cc5..45d12c4b6 100644
--- a/doc/requirements.txt
+++ b/doc/requirements.txt
@@ -18,5 +18,5 @@ sphinxcontrib-bibtex==0.4.2
 tensorflow==2.5.0
 tensornetwork==0.3
 toml
-torch==1.7.0+cpu
-torchvision==0.8.1+cpu
+torch==1.9.0+cpu
+torchvision==0.10.0+cpu
diff --git a/pennylane/grouping/group_observables.py b/pennylane/grouping/group_observables.py
index 5491721e3..ecbad3f5b 100644
--- a/pennylane/grouping/group_observables.py
+++ b/pennylane/grouping/group_observables.py
@@ -18,6 +18,7 @@ This module contains the high-level Pauli-word-partitioning functionality used i
 from copy import copy
 
 import numpy as np
+import pennylane as qml
 
 from pennylane.grouping.graph_colouring import largest_first, recursive_largest_first
 from pennylane.grouping.utils import (
@@ -191,7 +192,7 @@ def group_observables(observables, coefficients=None, grouping_type="qwc", metho
     Args:
         observables (list[Observable]): a list of Pauli word ``Observable`` instances (Pauli
             operation instances and :class:`~.Tensor` instances thereof)
-        coefficients (list[float]): A list of float coefficients. If not specified,
+        coefficients (tensor_like): A tensor or list of coefficients. If not specified,
             output ``partitioned_coeffs`` is not returned.
         grouping_type (str): The type of binary relation between Pauli words.
             Can be ``'qwc'``, ``'commuting'``, or ``'anticommuting'``.
@@ -203,9 +204,9 @@ def group_observables(observables, coefficients=None, grouping_type="qwc", metho
 
            * list[list[Observable]]: A list of the obtained groupings. Each grouping
              is itself a list of Pauli word ``Observable`` instances.
-           * list[list[float]]: A list of coefficient groupings. Each coefficient
-             grouping is itself a list of the grouping's corresponding coefficients. This is only
-             output if coefficients are specified.
+           * list[tensor_like]: A list of coefficient groupings. Each coefficient
+             grouping is itself a tensor or list of the grouping's corresponding coefficients. This is only
+             returned if coefficients are specified.
 
     Raises:
         IndexError: if the input list of coefficients is not of the same length as the input list
@@ -224,7 +225,7 @@ def group_observables(observables, coefficients=None, grouping_type="qwc", metho
     """
 
     if coefficients is not None:
-        if len(coefficients) != len(observables):
+        if qml.math.shape(coefficients)[0] != len(observables):
             raise IndexError(
                 "The coefficients list must be the same length as the observables list."
             )
@@ -237,18 +238,32 @@ def group_observables(observables, coefficients=None, grouping_type="qwc", metho
     if coefficients is None:
         return partitioned_paulis
 
-    partitioned_coeffs = [[0] * len(g) for g in partitioned_paulis]
+    partitioned_coeffs = [
+        qml.math.cast_like([0] * len(g), coefficients) for g in partitioned_paulis
+    ]
 
     observables = copy(observables)
-    coefficients = copy(coefficients)
+    # we cannot delete elements from the coefficients tensor, so we
+    # use a proxy list memorising the indices for this logic
+    coeff_indices = list(range(qml.math.shape(coefficients)[0]))
     for i, partition in enumerate(partitioned_paulis):
-        for j, pauli_word in enumerate(partition):
+        indices = []
+        for pauli_word in partition:
+            # find index of this pauli word in remaining original observables,
             for observable in observables:
                 if are_identical_pauli_words(pauli_word, observable):
                     ind = observables.index(observable)
-                    partitioned_coeffs[i][j] = coefficients[ind]
+                    indices.append(coeff_indices[ind])
                     observables.pop(ind)
-                    coefficients.pop(ind)
+                    coeff_indices.pop(ind)
                     break
 
+        # add a tensor of coefficients to the grouped coefficients
+        partitioned_coeffs[i] = qml.math.take(coefficients, indices, axis=0)
+
+    # make sure the output is of the same format as the input
+    # for these two frequent cases
+    if isinstance(coefficients, list):
+        partitioned_coeffs = [list(p) for p in partitioned_coeffs]
+
     return partitioned_paulis, partitioned_coeffs
diff --git a/pennylane/math/single_dispatch.py b/pennylane/math/single_dispatch.py
index 18464c4e5..8963b9f06 100644
--- a/pennylane/math/single_dispatch.py
+++ b/pennylane/math/single_dispatch.py
@@ -133,6 +133,7 @@ ar.autoray._SUBMODULE_ALIASES["tensorflow", "arccos"] = "tensorflow.math"
 ar.autoray._SUBMODULE_ALIASES["tensorflow", "arctan"] = "tensorflow.math"
 ar.autoray._SUBMODULE_ALIASES["tensorflow", "arctan2"] = "tensorflow.math"
 ar.autoray._SUBMODULE_ALIASES["tensorflow", "diag"] = "tensorflow.linalg"
+ar.autoray._SUBMODULE_ALIASES["tensorflow", "kron"] = "tensorflow.experimental.numpy"
 
 
 ar.autoray._FUNC_ALIASES["tensorflow", "arcsin"] = "asin"
diff --git a/pennylane/operation.py b/pennylane/operation.py
index 5fdc762f3..ea842caa0 100644
--- a/pennylane/operation.py
+++ b/pennylane/operation.py
@@ -1116,10 +1116,10 @@ class Observable(Operator):
         >>> ob1.compare(ob2)
         False
         """
-        if isinstance(other, (Tensor, Observable)):
-            return other._obs_data() == self._obs_data()
         if isinstance(other, qml.Hamiltonian):
             return other.compare(self)
+        if isinstance(other, (Tensor, Observable)):
+            return other._obs_data() == self._obs_data()
 
         raise ValueError(
             "Can only compare an Observable/Tensor, and a Hamiltonian/Observable/Tensor."
@@ -1127,12 +1127,10 @@ class Observable(Operator):
 
     def __add__(self, other):
         r"""The addition operation between Observables/Tensors/qml.Hamiltonian objects."""
-        if isinstance(other, (Observable, Tensor)):
-            return qml.Hamiltonian([1, 1], [self, other], simplify=True)
-
         if isinstance(other, qml.Hamiltonian):
             return other + self
-
+        if isinstance(other, (Observable, Tensor)):
+            return qml.Hamiltonian([1, 1], [self, other], simplify=True)
         raise ValueError(f"Cannot add Observable and {type(other)}")
 
     def __mul__(self, a):
diff --git a/pennylane/qnode.py b/pennylane/qnode.py
index 29e6b3425..f6542056b 100644
--- a/pennylane/qnode.py
+++ b/pennylane/qnode.py
@@ -617,8 +617,8 @@ class QNode:
                 tapes, fn = qml.transforms.hamiltonian_expand(self.qtape, group=False)
             except ValueError as e:
                 raise ValueError(
-                    "At the moment, only single expectations of Hamiltonian observables can be returned"
-                    "on the {} device.".format(self.device.name)
+                    "Only a single expectation of a Hamiltonian observable can be returned"
+                    "when using the {} device.".format(self.device.name)
                 ) from e
             results = [tape.execute(device=self.device) for tape in tapes]
             res = fn(results)
diff --git a/pennylane/transforms/hamiltonian_expand.py b/pennylane/transforms/hamiltonian_expand.py
index 927375fef..590c8dd53 100644
--- a/pennylane/transforms/hamiltonian_expand.py
+++ b/pennylane/transforms/hamiltonian_expand.py
@@ -97,10 +97,15 @@ def hamiltonian_expand(tape, group=True):
 
     hamiltonian = tape.measurements[0].obs
 
-    if not isinstance(hamiltonian, qml.Hamiltonian) or len(tape.measurements) > 1:
+    if (
+        not isinstance(hamiltonian, qml.Hamiltonian)
+        or len(tape.measurements) > 1
+        or tape.measurements[0].return_type != qml.measure.Expectation
+    ):
         raise ValueError(
             "Passed tape must end in `qml.expval(H)`, where H is of type `qml.Hamiltonian`"
         )
+
     if group:
         hamiltonian.simplify()
         return qml.transforms.measurement_grouping(tape, hamiltonian.ops, hamiltonian.coeffs)
@@ -108,10 +113,13 @@ def hamiltonian_expand(tape, group=True):
     # create tapes that measure the Pauli-words in the Hamiltonian
     tapes = []
     for ob in hamiltonian.ops:
-        new_tape = tape.copy()
-        new_tape._measurements = [
-            qml.measure.MeasurementProcess(return_type=qml.operation.Expectation, obs=ob)
-        ]
+        # we need to create a new tape here, because
+        # updating metadata of a copied tape is error-prone
+        # when the observables were changed
+        with tape.__class__() as new_tape:
+            for op in tape.operations:
+                qml.apply(op)
+            qml.expval(ob)
         tapes.append(new_tape)
 
     # create processing function that performs linear recombination
diff --git a/pennylane/vqe/vqe.py b/pennylane/vqe/vqe.py
index fef96c5e2..069e68155 100644
--- a/pennylane/vqe/vqe.py
+++ b/pennylane/vqe/vqe.py
@@ -17,8 +17,9 @@ computations using PennyLane.
 """
 # pylint: disable=too-many-arguments, too-few-public-methods
 from collections.abc import Sequence
-import itertools
 import warnings
+import itertools
+from copy import copy
 
 import pennylane as qml
 from pennylane import numpy as np
@@ -29,23 +30,18 @@ from pennylane.wires import Wires
 OBS_MAP = {"PauliX": "X", "PauliY": "Y", "PauliZ": "Z", "Hadamard": "H", "Identity": "I"}
 
 
-class Hamiltonian:
-    r"""Lightweight class for representing Hamiltonians for Variational Quantum
-    Eigensolver problems.
-
-    Hamiltonians can be expressed as linear combinations of observables, e.g.,
-    :math:`\sum_{k=0}^{N-1} c_k O_k`.
+class Hamiltonian(qml.operation.Observable):
+    r"""Operator representing a Hamiltonian.
 
-    This class keeps track of the terms (coefficients and observables) separately.
+    The Hamiltonian is represented as a linear combination of other operators, e.g.,
+    :math:`\sum_{k=0}^{N-1} c_k O_k`, where the :math:`c_k` are trainable parameters.
 
     Args:
-        coeffs (Iterable[float]): coefficients of the Hamiltonian expression
-        observables (Iterable[Observable]): observables in the Hamiltonian expression
+        coeffs (tensor_like): coefficients of the Hamiltonian expression
+        observables (Iterable[Observable]): observables in the Hamiltonian expression, of same length as coeffs
         simplify (bool): Specifies whether the Hamiltonian is simplified upon initialization
                          (like-terms are combined). The default value is `False`.
 
-    .. seealso:: :class:`~.ExpvalCost`, :func:`~.molecular_hamiltonian`
-
     **Example:**
 
     A Hamiltonian can be created by simply passing the list of coefficients
@@ -58,6 +54,15 @@ class Hamiltonian:
       (-0.543) [Z0 H2]
     + (0.2) [X0 Z1]
 
+    The coefficients can be a trainable tensor, for example:
+
+    >>> coeffs = tf.Variable([0.2, -0.543], dtype=tf.double)
+    >>> obs = [qml.PauliX(0) @ qml.PauliZ(1), qml.PauliZ(0) @ qml.Hadamard(2)]
+    >>> H = qml.Hamiltonian(coeffs, obs)
+    >>> print(H)
+      (-0.543) [Z0 H2]
+    + (0.2) [X0 Z1]
+
     The user can also provide custom observables:
 
     >>> obs_matrix = np.array([[0.5, 1.0j, 0.0, -3j],
@@ -73,54 +78,62 @@ class Hamiltonian:
     :doc:`/introduction/chemistry` module can be used to generate a molecular
     Hamiltonian.
 
-    .. Warning::
+    In many cases, Hamiltonians can be constructed using Pythonic arithmetic operations.
+    For example:
 
-        Hamiltonians can be constructed using Pythonic arithmetic operations. For example:
+    >>> qml.Hamiltonian([1.], [qml.PauliX(0)]) + 2 * qml.PauliZ(0) @ qml.PauliZ(1)
 
-        >>> qml.PauliX(0) + 2 * qml.PauliZ(0) @ qml.PauliZ(1)
+    is equivalent to the following Hamiltonian:
 
-        is equivalent to the following Hamiltonian:
+    >>> qml.Hamiltonian([1, 2], [qml.PauliX(0), qml.PauliZ(0) @ qml.PauliZ(1)])
 
-        >>> qml.Hamiltonian([1, 2], [qml.PauliX(0), qml.PauliZ(0) @ qml.PauliZ(1)])
+    While scalar multiplication requires native python floats or integer types,
+    addition, subtraction, and tensor multiplication of Hamiltonians with Hamiltonians or
+    other observables is possible with tensor-valued coefficients, i.e.,
 
-        When Hamiltonians are defined using arithmetic operations **inside of QNodes**, constituent observables
-        may be queued as operations/an error may be thrown. Thus, Hamiltonians must be defined either outside of QNodes,
-        or inside of QNodes using the conventional method.
-
-        Note that this issue also arises when calling the ``simplify()`` method.
+    >>> H1 = qml.Hamiltonian(torch.tensor([1.]), [qml.PauliX(0)])
+    >>> H2 = qml.Hamiltonian(torch.tensor([2., 3.]), [qml.PauliY(0), qml.PauliX(1)])
+    >>> H3 = qml.Hamiltonian(torch.tensor([1., 2., 3.]), [qml.PauliX(0), qml.PauliY(0), qml.PauliX(1)])
+    >>> H3.compare(H1 + H2)
+    True
     """
-    # Todo: this is a temporary solution to make the circuit drawer work
-    num_params = 0
 
-    def __init__(self, coeffs, observables, simplify=False):
+    num_wires = qml.operation.AnyWires
+    num_params = 1
+    par_domain = "A"
+    grad_method = "A"  # supports analytic gradients
+
+    def __init__(self, coeffs, observables, simplify=False, id=None, do_queue=True):
 
-        if len(coeffs) != len(observables):
+        if qml.math.shape(coeffs)[0] != len(observables):
             raise ValueError(
                 "Could not create valid Hamiltonian; "
                 "number of coefficients and operators does not match."
             )
 
-        if any(np.imag(coeffs) != 0):
-            raise ValueError(
-                "Could not create valid Hamiltonian; " "coefficients are not real-valued."
-            )
-
         for obs in observables:
             if not isinstance(obs, Observable):
                 raise ValueError(
                     "Could not create circuits. Some or all observables are not valid."
                 )
 
-        self._coeffs = list(coeffs)
+        self._coeffs = coeffs
         self._ops = list(observables)
+        self._wires = qml.wires.Wires.all_wires([op.wires for op in self.ops], sort=True)
 
-        self.data = []
         self.return_type = None
 
         if simplify:
             self.simplify()
 
-        self.queue()
+        coeffs_flat = [self._coeffs[i] for i in range(qml.math.shape(self._coeffs)[0])]
+        # overwrite this attribute, now that we have the correct info
+        self.num_params = qml.math.shape(self._coeffs)[0]
+
+        # create the operator using each coefficient as a separate parameter;
+        # this causes H.data to be a list of tensor scalars,
+        # while H.coeffs is the original tensor
+        super().__init__(*coeffs_flat, wires=self._wires, id=id, do_queue=do_queue)
 
     @property
     def coeffs(self):
@@ -156,7 +169,7 @@ class Hamiltonian:
         Returns:
             (Wires): Combined wires present in all terms, sorted.
         """
-        return qml.wires.Wires.all_wires([op.wires for op in self.ops], sort=True)
+        return self._wires
 
     @property
     def name(self):
@@ -174,35 +187,39 @@ class Hamiltonian:
           (-1) [X0]
         + (1) [Y2]
         """
-        coeffs = []
+        data = []
         ops = []
 
-        for c, op in zip(self.coeffs, self.ops):
+        for i in range(len(self.ops)):  # pylint: disable=consider-using-enumerate
+            op = self.ops[i]
+            c = self.coeffs[i]
             op = op if isinstance(op, Tensor) else Tensor(op)
 
             ind = None
-            for i, other in enumerate(ops):
-                if op.compare(other):
-                    ind = i
+            for j, o in enumerate(ops):
+                if op.compare(o):
+                    ind = j
                     break
 
             if ind is not None:
-                coeffs[ind] += c
-                if np.allclose([coeffs[ind]], [0]):
-                    del coeffs[ind]
+                data[ind] += c
+                if np.isclose(qml.math.toarray(data[ind]), np.array(0.0)):
+                    del data[ind]
                     del ops[ind]
             else:
                 ops.append(op.prune())
-                coeffs.append(c)
+                data.append(c)
 
-        self._coeffs = coeffs
+        self._coeffs = qml.math.stack(data) if data else []
+        self.data = data
         self._ops = ops
 
     def __str__(self):
         # Lambda function that formats the wires
         wires_print = lambda ob: ",".join(map(str, ob.wires.tolist()))
 
-        paired_coeff_obs = list(zip(self.coeffs, self.ops))
+        list_of_coeffs = self.data  # list of scalar tensors
+        paired_coeff_obs = list(zip(list_of_coeffs, self.ops))
         paired_coeff_obs.sort(key=lambda pair: (len(pair[1].wires), pair[0]))
 
         terms_ls = []
@@ -223,7 +240,7 @@ class Hamiltonian:
 
     def __repr__(self):
         # Constructor-call-like representation
-        return f"<Hamiltonian: terms={len(self.coeffs)}, wires={self.wires.tolist()}>"
+        return f"<Hamiltonian: terms={qml.math.shape(self.coeffs)[0]}, wires={self.wires.tolist()}>"
 
     def _obs_data(self):
         r"""Extracts the data from a Hamiltonian and serializes it in an order-independent fashion.
@@ -247,12 +264,13 @@ class Hamiltonian:
         """
         data = set()
 
-        for co, op in zip(*self.terms):
+        coeffs_arr = qml.math.toarray(self.coeffs)
+        for co, op in zip(coeffs_arr, self.ops):
             obs = op.non_identity_obs if isinstance(op, Tensor) else [op]
             tensor = []
             for ob in obs:
                 parameters = tuple(
-                    param.tostring() for param in ob.parameters
+                    str(param) for param in ob.parameters
                 )  # Converts params into immutable type
                 tensor.append((ob.name, ob.wires, parameters))
             data.add((co, frozenset(tensor)))
@@ -314,8 +332,8 @@ class Hamiltonian:
 
     def __matmul__(self, H):
         r"""The tensor product operation between a Hamiltonian and a Hamiltonian/Tensor/Observable."""
-        coeffs1 = self.coeffs.copy()
-        terms1 = self.ops.copy()
+        coeffs1 = copy(self.coeffs)
+        ops1 = self.ops.copy()
 
         if isinstance(H, Hamiltonian):
             shared_wires = Wires.shared_wires([self.wires, H.wires])
@@ -326,34 +344,35 @@ class Hamiltonian:
                 )
 
             coeffs2 = H.coeffs
-            terms2 = H.ops
+            ops2 = H.ops
 
-            coeffs = [c[0] * c[1] for c in itertools.product(coeffs1, coeffs2)]
-            term_list = itertools.product(terms1, terms2)
-            terms = [qml.operation.Tensor(t[0], t[1]) for t in term_list]
+            coeffs = qml.math.kron(coeffs1, coeffs2)
+            ops_list = itertools.product(ops1, ops2)
+            terms = [qml.operation.Tensor(t[0], t[1]) for t in ops_list]
 
             return qml.Hamiltonian(coeffs, terms, simplify=True)
 
         if isinstance(H, (Tensor, Observable)):
-            coeffs = coeffs1
-            terms = [term @ H for term in terms1]
+            terms = [op @ H for op in ops1]
 
-            return qml.Hamiltonian(coeffs, terms, simplify=True)
+            return qml.Hamiltonian(coeffs1, terms, simplify=True)
 
         raise ValueError(f"Cannot tensor product Hamiltonian and {type(H)}")
 
     def __add__(self, H):
         r"""The addition operation between a Hamiltonian and a Hamiltonian/Tensor/Observable."""
-        coeffs = self.coeffs.copy()
         ops = self.ops.copy()
+        self_coeffs = copy(self.coeffs)
 
         if isinstance(H, Hamiltonian):
-            coeffs.extend(H.coeffs.copy())
+            coeffs = qml.math.concatenate([self_coeffs, copy(H.coeffs)], axis=0)
             ops.extend(H.ops.copy())
             return qml.Hamiltonian(coeffs, ops, simplify=True)
 
         if isinstance(H, (Tensor, Observable)):
-            coeffs.append(1)
+            coeffs = qml.math.concatenate(
+                [self_coeffs, qml.math.cast_like([1.0], self_coeffs)], axis=0
+            )
             ops.append(H)
             return qml.Hamiltonian(coeffs, ops, simplify=True)
 
@@ -362,7 +381,8 @@ class Hamiltonian:
     def __mul__(self, a):
         r"""The scalar multiplication operation between a scalar and a Hamiltonian."""
         if isinstance(a, (int, float)):
-            coeffs = [a * c for c in self.coeffs.copy()]
+            self_coeffs = copy(self.coeffs)
+            coeffs = qml.math.multiply(qml.math.cast_like([a], self_coeffs), self_coeffs)
             return qml.Hamiltonian(coeffs, self.ops.copy())
 
         raise ValueError(f"Cannot multiply Hamiltonian by {type(a)}")
@@ -378,13 +398,15 @@ class Hamiltonian:
     def __iadd__(self, H):
         r"""The inplace addition operation between a Hamiltonian and a Hamiltonian/Tensor/Observable."""
         if isinstance(H, Hamiltonian):
-            self._coeffs.extend(H.coeffs.copy())
+            self._coeffs = qml.math.concatenate([self._coeffs, H.coeffs], axis=0)
             self._ops.extend(H.ops.copy())
             self.simplify()
             return self
 
         if isinstance(H, (Tensor, Observable)):
-            self._coeffs.append(1)
+            self._coeffs = qml.math.concatenate(
+                [self._coeffs, qml.math.cast_like([1.0], self._coeffs)], axis=0
+            )
             self._ops.append(H)
             self.simplify()
             return self
@@ -394,7 +416,7 @@ class Hamiltonian:
     def __imul__(self, a):
         r"""The inplace scalar multiplication operation between a scalar and a Hamiltonian."""
         if isinstance(a, (int, float)):
-            self._coeffs = [a * c for c in self.coeffs]
+            self._coeffs = qml.math.multiply(qml.math.cast_like([a], self._coeffs), self._coeffs)
             return self
 
         raise ValueError(f"Cannot multiply Hamiltonian by {type(a)}")
@@ -547,7 +569,7 @@ class ExpvalCost:
         self._multiple_devices = isinstance(device, Sequence)
         """Bool: Records if multiple devices are input"""
 
-        if all(c == 0 for c in coeffs) or not coeffs:
+        if np.isclose(qml.math.toarray(qml.math.count_nonzero(coeffs)), 0):
             self.cost_fn = lambda *args, **kwargs: np.array(0)
             return
 


Test patch:
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index bb8c152a5..e0ea98730 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -8,7 +8,7 @@ on:
 
 env:
   TF_VERSION: 2.4
-  TORCH_VERSION: 1.7
+  TORCH_VERSION: 1.9.0+cpu
   COVERAGE_FLAGS: "--cov=pennylane --cov-report=term-missing --cov-report=xml --no-flaky-report -p no:warnings --tb=native"
 
 
diff --git a/tests/grouping/test_group_observables.py b/tests/grouping/test_group_observables.py
index 72156c8dd..ffb399fa5 100644
--- a/tests/grouping/test_group_observables.py
+++ b/tests/grouping/test_group_observables.py
@@ -16,9 +16,11 @@ Unit tests for ``PauliGroupingStrategy`` and ``group_observables`` in ``grouping
 """
 import pytest
 import numpy as np
+import pennylane as qml
 from pennylane import Identity, PauliX, PauliY, PauliZ
 from pennylane.grouping.utils import are_identical_pauli_words
 from pennylane.grouping.group_observables import PauliGroupingStrategy, group_observables
+from pennylane import numpy as pnp
 
 
 class TestPauliGroupingStrategy:
@@ -301,3 +303,93 @@ class TestGroupObservables:
         _ = grouping_instance.binary_repr(n_qubits, wire_map)
 
         assert grouping_instance._wire_map == wire_map
+
+    def test_return_list_coefficients(self):
+        """Tests that if the coefficients are given as a list, the groups
+        are likewise lists."""
+        obs = [qml.PauliX(0), qml.PauliX(1)]
+        coeffs = [1.0, 2.0]
+        _, grouped_coeffs = group_observables(obs, coeffs)
+        assert isinstance(grouped_coeffs[0], list)
+
+
+class TestDifferentiable:
+    """Tests that grouping observables is differentiable with respect to the coefficients."""
+
+    def test_differentiation_autograd(self, tol):
+        """Test that grouping is differentiable with autograd tensors as coefficient"""
+        coeffs = pnp.array([1.0, 2.0, 3.0], requires_grad=True)
+        obs = [PauliX(wires=0), PauliX(wires=1), PauliZ(wires=1)]
+
+        def group(coeffs, select=None):
+            _, grouped_coeffs = qml.grouping.group_observables(obs, coeffs)
+            return grouped_coeffs[select]
+
+        jac_fn = qml.jacobian(group)
+        assert pnp.allclose(
+            jac_fn(coeffs, select=0), pnp.array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]), atol=tol
+        )
+        assert pnp.allclose(jac_fn(coeffs, select=1), pnp.array([[0.0, 0.0, 1.0]]), atol=tol)
+
+    def test_differentiation_jax(self, tol):
+        """Test that grouping is differentiable with jax tensors as coefficient"""
+        jax = pytest.importorskip("jax")
+        jnp = pytest.importorskip("jax.numpy")
+        coeffs = jnp.array([1.0, 2.0, 3.0])
+        obs = [PauliX(wires=0), PauliX(wires=1), PauliZ(wires=1)]
+
+        def group(coeffs, select=None):
+            _, grouped_coeffs = qml.grouping.group_observables(obs, coeffs)
+            return grouped_coeffs[select]
+
+        jac_fn = jax.jacobian(group)
+        assert np.allclose(
+            jac_fn(coeffs, select=0), pnp.array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]), atol=tol
+        )
+        assert np.allclose(jac_fn(coeffs, select=1), pnp.array([[0.0, 0.0, 1.0]]), atol=tol)
+
+    def test_differentiation_torch(self, tol):
+        """Test that grouping is differentiable with torch tensors as coefficient"""
+        torch = pytest.importorskip("torch")
+        obs = [PauliX(wires=0), PauliX(wires=1), PauliZ(wires=1)]
+
+        def group(coeffs, select_group=None, select_index=None):
+            # we return a scalar, since torch is best at computing gradients
+            _, grouped_coeffs = qml.grouping.group_observables(obs, coeffs)
+            return grouped_coeffs[select_group][select_index]
+
+        coeffs = torch.tensor([1.0, 2.0, 3.0], requires_grad=True)
+        res = group(coeffs, select_group=0, select_index=0)
+        res.backward()
+        assert np.allclose(coeffs.grad, [1.0, 0.0, 0.0], atol=tol)
+
+        coeffs = torch.tensor([1.0, 2.0, 3.0], requires_grad=True)
+        res = group(coeffs, select_group=0, select_index=1)
+        res.backward()
+        assert np.allclose(coeffs.grad, [0.0, 1.0, 0.0], atol=tol)
+
+        coeffs = torch.tensor([1.0, 2.0, 3.0], requires_grad=True)
+        res = group(coeffs, select_group=1, select_index=0)
+        res.backward()
+        assert np.allclose(coeffs.grad, [0.0, 0.0, 1.0], atol=tol)
+
+    def test_differentiation_tf(self, tol):
+        """Test that grouping is differentiable with tf tensors as coefficient"""
+        tf = pytest.importorskip("tf")
+        obs = [PauliX(wires=0), PauliX(wires=1), PauliZ(wires=1)]
+
+        def group(coeffs, select=None):
+            _, grouped_coeffs = qml.grouping.group_observables(obs, coeffs)
+            return grouped_coeffs[select]
+
+        coeffs = tf.Variable([1.0, 2.0, 3.0], dtype=tf.double)
+
+        with tf.GradientTape() as tape:
+            res = group(coeffs, select=0)
+        grad = tape.jacobian(res, [coeffs])
+        assert np.allclose(grad, pnp.array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]), atol=tol)
+
+        with tf.GradientTape() as tape:
+            res = group(coeffs, select=1)
+        grad = tape.jacobian(res, [coeffs])
+        assert np.allclose(grad, pnp.array([[0.0, 0.0, 1.0]]), atol=tol)
diff --git a/tests/grouping/test_optimize_measurements.py b/tests/grouping/test_optimize_measurements.py
index a8ab748f1..14220b89e 100644
--- a/tests/grouping/test_optimize_measurements.py
+++ b/tests/grouping/test_optimize_measurements.py
@@ -15,6 +15,7 @@
 Unit tests for ``optimize_measurements`` function in ``grouping/optimize_measurements.py``.
 """
 import pytest
+import numpy as np
 from pennylane import Identity, PauliX, PauliY, PauliZ
 from pennylane.grouping.utils import are_identical_pauli_words
 from pennylane.grouping.optimize_measurements import optimize_measurements
@@ -108,7 +109,8 @@ class TestOptimizeMeasurements:
         assert len(grouped_coeffs) == len(grouped_coeffs)
 
         assert all(
-            grouped_coeffs[i] == grouped_coeffs_sol[i] for i in range(len(grouped_coeffs_sol))
+            np.allclose(grouped_coeffs[i], grouped_coeffs_sol[i])
+            for i in range(len(grouped_coeffs_sol))
         )
 
     @pytest.mark.parametrize(
@@ -140,7 +142,8 @@ class TestOptimizeMeasurements:
         assert len(grouped_coeffs) == len(grouped_coeffs)
 
         assert all(
-            grouped_coeffs[i] == grouped_coeffs_sol[i] for i in range(len(grouped_coeffs_sol))
+            np.allclose(grouped_coeffs[i], grouped_coeffs_sol[i])
+            for i in range(len(grouped_coeffs_sol))
         )
 
     def test_optimize_measurements_not_implemented_catch(self):
diff --git a/tests/ops/test_hamiltonian.py b/tests/ops/test_hamiltonian.py
new file mode 100644
index 000000000..2d429b9dc
--- /dev/null
+++ b/tests/ops/test_hamiltonian.py
@@ -0,0 +1,642 @@
+# Copyright 2018-2021 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.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+Tests for the Hamiltonian class.
+"""
+import numpy as np
+import pytest
+
+import pennylane as qml
+from pennylane import numpy as pnp
+
+# Make test data in different interfaces, if installed
+COEFFS_PARAM_INTERFACE = [
+    ([-0.05, 0.17], 1.7, "autograd"),
+    (np.array([-0.05, 0.17]), np.array(1.7), "autograd"),
+    (pnp.array([-0.05, 0.17], requires_grad=True), pnp.array(1.7, requires_grad=True), "autograd"),
+]
+
+try:
+    from jax import numpy as jnp
+
+    COEFFS_PARAM_INTERFACE.append((jnp.array([-0.05, 0.17]), jnp.array(1.7), "jax"))
+except ImportError:
+    pass
+
+try:
+    import tf
+
+    COEFFS_PARAM_INTERFACE.append(
+        (tf.Variable([-0.05, 0.17], dtype=tf.double), tf.Variable(1.7, dtype=tf.double), "tf")
+    )
+except ImportError:
+    pass
+
+try:
+    import torch
+
+    COEFFS_PARAM_INTERFACE.append((torch.tensor([-0.05, 0.17]), torch.tensor([1.7]), "torch"))
+except ImportError:
+    pass
+
+
+def circuit1(param):
+    """First Pauli subcircuit"""
+    qml.RX(param, wires=0)
+    qml.RY(param, wires=0)
+    return qml.expval(qml.PauliX(0))
+
+
+def circuit2(param):
+    """Second Pauli subcircuit"""
+    qml.RX(param, wires=0)
+    qml.RY(param, wires=0)
+    return qml.expval(qml.PauliZ(0))
+
+
+dev = qml.device("default.qubit", wires=2)
+
+
+class TestHamiltonianCoefficients:
+    """Test the creation of a Hamiltonian"""
+
+    @pytest.mark.parametrize("coeffs", [el[0] for el in COEFFS_PARAM_INTERFACE])
+    def test_creation_different_coeff_types(self, coeffs):
+        """Check that Hamiltonian's coefficients and data attributes are set correctly."""
+        H = qml.Hamiltonian(coeffs, [qml.PauliX(0), qml.PauliZ(0)])
+        assert np.allclose(coeffs, H.coeffs)
+        assert np.allclose([coeffs[i] for i in range(qml.math.shape(coeffs)[0])], H.data)
+
+    @pytest.mark.parametrize("coeffs", [el[0] for el in COEFFS_PARAM_INTERFACE])
+    def test_simplify(self, coeffs):
+        """Test that simplify works with different coefficient types."""
+        H1 = qml.Hamiltonian(coeffs, [qml.PauliX(0), qml.PauliZ(1)])
+        H2 = qml.Hamiltonian(coeffs, [qml.PauliX(0), qml.Identity(0) @ qml.PauliZ(1)])
+        H2.simplify()
+        assert H1.compare(H2)
+        assert H1.data == H2.data
+
+
+class TestHamiltonianArithmeticTF:
+    """Tests creation of Hamiltonians using arithmetic
+    operations with TensorFlow tensor coefficients."""
+
+    def test_hamiltonian_equal(self):
+        """Tests equality"""
+        tf = pytest.importorskip("tensorflow")
+
+        coeffs = tf.Variable([0.5, -1.6])
+        obs = [qml.PauliX(0), qml.PauliY(1)]
+        H1 = qml.Hamiltonian(coeffs, obs)
+
+        coeffs2 = tf.Variable([-1.6, 0.5])
+        obs2 = [qml.PauliY(1), qml.PauliX(0)]
+        H2 = qml.Hamiltonian(coeffs2, obs2)
+
+        assert H1.compare(H2)
+
+    def test_hamiltonian_add(self):
+        """Tests that Hamiltonians are added correctly"""
+        tf = pytest.importorskip("tensorflow")
+
+        coeffs = tf.Variable([0.5, -1.6])
+        obs = [qml.PauliX(0), qml.PauliY(1)]
+        H1 = qml.Hamiltonian(coeffs, obs)
+
+        coeffs2 = tf.Variable([0.5, -0.4])
+        H2 = qml.Hamiltonian(coeffs2, obs)
+
+        coeffs_expected = tf.Variable([1.0, -2.0])
+        H = qml.Hamiltonian(coeffs_expected, obs)
+
+        assert H.compare(H1 + H2)
+
+        H1 += H2
+        assert H.compare(H1)
+
+    def test_hamiltonian_sub(self):
+        """Tests that Hamiltonians are subtracted correctly"""
+        tf = pytest.importorskip("tensorflow")
+
+        coeffs = tf.Variable([1.0, -2.0])
+        obs = [qml.PauliX(0), qml.PauliY(1)]
+        H1 = qml.Hamiltonian(coeffs, obs)
+
+        coeffs2 = tf.Variable([0.5, -0.4])
+        H2 = qml.Hamiltonian(coeffs2, obs)
+
+        coeffs_expected = tf.Variable([0.5, -1.6])
+        H = qml.Hamiltonian(coeffs_expected, obs)
+
+        assert H.compare(H1 - H2)
+
+        H1 -= H2
+        assert H.compare(H1)
+
+    def test_hamiltonian_matmul(self):
+        """Tests that Hamiltonians are tensored correctly"""
+        tf = pytest.importorskip("tensorflow")
+
+        coeffs = tf.Variable([1.0, 2.0])
+        obs = [qml.PauliX(0), qml.PauliY(1)]
+        H1 = qml.Hamiltonian(coeffs, obs)
+
+        coeffs2 = tf.Variable([-1.0, -2.0])
+        obs2 = [qml.PauliX(2), qml.PauliY(3)]
+        H2 = qml.Hamiltonian(coeffs2, obs2)
+
+        coeffs_expected = tf.Variable([-4.0, -2.0, -2.0, -1.0])
+        obs_expected = [
+            qml.PauliY(1) @ qml.PauliY(3),
+            qml.PauliX(0) @ qml.PauliY(3),
+            qml.PauliX(2) @ qml.PauliY(1),
+            qml.PauliX(0) @ qml.PauliX(2),
+        ]
+        H = qml.Hamiltonian(coeffs_expected, obs_expected)
+
+        assert H.compare(H1 @ H2)
+
+
+class TestHamiltonianArithmeticTorch:
+    """Tests creation of Hamiltonians using arithmetic
+    operations with torch tensor coefficients."""
+
+    def test_hamiltonian_equal(self):
+        """Tests equality"""
+        torch = pytest.importorskip("torch")
+
+        coeffs = torch.tensor([0.5, -1.6])
+        obs = [qml.PauliX(0), qml.PauliY(1)]
+        H1 = qml.Hamiltonian(coeffs, obs)
+
+        coeffs2 = torch.tensor([-1.6, 0.5])
+        obs2 = [qml.PauliY(1), qml.PauliX(0)]
+        H2 = qml.Hamiltonian(coeffs2, obs2)
+
+        assert H1.compare(H2)
+
+    def test_hamiltonian_add(self):
+        """Tests that Hamiltonians are added correctly"""
+        torch = pytest.importorskip("torch")
+
+        coeffs = torch.tensor([0.5, -1.6])
+        obs = [qml.PauliX(0), qml.PauliY(1)]
+        H1 = qml.Hamiltonian(coeffs, obs)
+
+        coeffs2 = torch.tensor([0.5, -0.4])
+        H2 = qml.Hamiltonian(coeffs2, obs)
+
+        coeffs_expected = torch.tensor([1.0, -2.0])
+        H = qml.Hamiltonian(coeffs_expected, obs)
+
+        assert H.compare(H1 + H2)
+
+        H1 += H2
+        assert H.compare(H1)
+
+    def test_hamiltonian_sub(self):
+        """Tests that Hamiltonians are subtracted correctly"""
+        torch = pytest.importorskip("torch")
+
+        coeffs = torch.tensor([1.0, -2.0])
+        obs = [qml.PauliX(0), qml.PauliY(1)]
+        H1 = qml.Hamiltonian(coeffs, obs)
+
+        coeffs2 = torch.tensor([0.5, -0.4])
+        H2 = qml.Hamiltonian(coeffs2, obs)
+
+        coeffs_expected = torch.tensor([0.5, -1.6])
+        H = qml.Hamiltonian(coeffs_expected, obs)
+
+        assert H.compare(H1 - H2)
+
+        H1 -= H2
+        assert H.compare(H1)
+
+    def test_hamiltonian_matmul(self):
+        """Tests that Hamiltonians are tensored correctly"""
+        torch = pytest.importorskip("torch")
+
+        coeffs = torch.tensor([1.0, 2.0])
+        obs = [qml.PauliX(0), qml.PauliY(1)]
+        H1 = qml.Hamiltonian(coeffs, obs)
+
+        coeffs2 = torch.tensor([-1.0, -2.0])
+        obs2 = [qml.PauliX(2), qml.PauliY(3)]
+        H2 = qml.Hamiltonian(coeffs2, obs2)
+
+        coeffs_expected = torch.tensor([-4.0, -2.0, -2.0, -1.0])
+        obs_expected = [
+            qml.PauliY(1) @ qml.PauliY(3),
+            qml.PauliX(0) @ qml.PauliY(3),
+            qml.PauliX(2) @ qml.PauliY(1),
+            qml.PauliX(0) @ qml.PauliX(2),
+        ]
+        H = qml.Hamiltonian(coeffs_expected, obs_expected)
+
+        assert H.compare(H1 @ H2)
+
+
+class TestHamiltonianArithmeticAutograd:
+    """Tests creation of Hamiltonians using arithmetic
+    operations with autograd tensor coefficients."""
+
+    def test_hamiltonian_equal(self):
+        """Tests equality"""
+        coeffs = pnp.array([0.5, -1.6])
+        obs = [qml.PauliX(0), qml.PauliY(1)]
+        H1 = qml.Hamiltonian(coeffs, obs)
+
+        coeffs2 = pnp.array([-1.6, 0.5])
+        obs2 = [qml.PauliY(1), qml.PauliX(0)]
+        H2 = qml.Hamiltonian(coeffs2, obs2)
+
+        assert H1.compare(H2)
+
+    def test_hamiltonian_add(self):
+        """Tests that Hamiltonians are added correctly"""
+        coeffs = pnp.array([0.5, -1.6])
+        obs = [qml.PauliX(0), qml.PauliY(1)]
+        H1 = qml.Hamiltonian(coeffs, obs)
+
+        coeffs2 = pnp.array([0.5, -0.4])
+        H2 = qml.Hamiltonian(coeffs2, obs)
+
+        coeffs_expected = pnp.array([1.0, -2.0])
+        H = qml.Hamiltonian(coeffs_expected, obs)
+
+        assert H.compare(H1 + H2)
+
+        H1 += H2
+        assert H.compare(H1)
+
+    def test_hamiltonian_sub(self):
+        """Tests that Hamiltonians are subtracted correctly"""
+        coeffs = pnp.array([1.0, -2.0])
+        obs = [qml.PauliX(0), qml.PauliY(1)]
+        H1 = qml.Hamiltonian(coeffs, obs)
+
+        coeffs2 = pnp.array([0.5, -0.4])
+        H2 = qml.Hamiltonian(coeffs2, obs)
+
+        coeffs_expected = pnp.array([0.5, -1.6])
+        H = qml.Hamiltonian(coeffs_expected, obs)
+
+        assert H.compare(H1 - H2)
+
+        H1 -= H2
+        assert H.compare(H1)
+
+    def test_hamiltonian_matmul(self):
+        """Tests that Hamiltonians are tensored correctly"""
+        coeffs = pnp.array([1.0, 2.0])
+        obs = [qml.PauliX(0), qml.PauliY(1)]
+        H1 = qml.Hamiltonian(coeffs, obs)
+
+        coeffs2 = pnp.array([-1.0, -2.0])
+        obs2 = [qml.PauliX(2), qml.PauliY(3)]
+        H2 = qml.Hamiltonian(coeffs2, obs2)
+
+        coeffs_expected = pnp.array([-4.0, -2.0, -2.0, -1.0])
+        obs_expected = [
+            qml.PauliY(1) @ qml.PauliY(3),
+            qml.PauliX(0) @ qml.PauliY(3),
+            qml.PauliX(2) @ qml.PauliY(1),
+            qml.PauliX(0) @ qml.PauliX(2),
+        ]
+        H = qml.Hamiltonian(coeffs_expected, obs_expected)
+
+        assert H.compare(H1 @ H2)
+
+
+class TestHamiltonianArithmeticJax:
+    """Tests creation of Hamiltonians using arithmetic
+    operations with jax tensor coefficients."""
+
+    def test_hamiltonian_equal(self):
+        """Tests equality"""
+        jax = pytest.importorskip("jax")
+        from jax import numpy as jnp
+
+        coeffs = jnp.array([0.5, -1.6])
+        obs = [qml.PauliX(0), qml.PauliY(1)]
+        H1 = qml.Hamiltonian(coeffs, obs)
+
+        coeffs2 = jnp.array([-1.6, 0.5])
+        obs2 = [qml.PauliY(1), qml.PauliX(0)]
+        H2 = qml.Hamiltonian(coeffs2, obs2)
+
+        assert H1.compare(H2)
+
+    def test_hamiltonian_add(self):
+        """Tests that Hamiltonians are added correctly"""
+        jax = pytest.importorskip("jax")
+        from jax import numpy as jnp
+
+        coeffs = jnp.array([0.5, -1.6])
+        obs = [qml.PauliX(0), qml.PauliY(1)]
+        H1 = qml.Hamiltonian(coeffs, obs)
+
+        coeffs2 = jnp.array([0.5, -0.4])
+        H2 = qml.Hamiltonian(coeffs2, obs)
+
+        coeffs_expected = jnp.array([1.0, -2.0])
+        H = qml.Hamiltonian(coeffs_expected, obs)
+
+        assert H.compare(H1 + H2)
+
+        H1 += H2
+        assert H.compare(H1)
+
+    def test_hamiltonian_sub(self):
+        """Tests that Hamiltonians are subtracted correctly"""
+        jax = pytest.importorskip("jax")
+        from jax import numpy as jnp
+
+        coeffs = jnp.array([1.0, -2.0])
+        obs = [qml.PauliX(0), qml.PauliY(1)]
+        H1 = qml.Hamiltonian(coeffs, obs)
+
+        coeffs2 = jnp.array([0.5, -0.4])
+        H2 = qml.Hamiltonian(coeffs2, obs)
+
+        coeffs_expected = jnp.array([0.5, -1.6])
+        H = qml.Hamiltonian(coeffs_expected, obs)
+
+        assert H.compare(H1 - H2)
+
+        H1 -= H2
+        assert H.compare(H1)
+
+    def test_hamiltonian_matmul(self):
+        """Tests that Hamiltonians are tensored correctly"""
+        jax = pytest.importorskip("jax")
+        from jax import numpy as jnp
+
+        coeffs = jnp.array([1.0, 2.0])
+        obs = [qml.PauliX(0), qml.PauliY(1)]
+        H1 = qml.Hamiltonian(coeffs, obs)
+
+        coeffs2 = jnp.array([-1.0, -2.0])
+        obs2 = [qml.PauliX(2), qml.PauliY(3)]
+        H2 = qml.Hamiltonian(coeffs2, obs2)
+
+        coeffs_expected = jnp.array([-4.0, -2.0, -2.0, -1.0])
+        obs_expected = [
+            qml.PauliY(1) @ qml.PauliY(3),
+            qml.PauliX(0) @ qml.PauliY(3),
+            qml.PauliX(2) @ qml.PauliY(1),
+            qml.PauliX(0) @ qml.PauliX(2),
+        ]
+        H = qml.Hamiltonian(coeffs_expected, obs_expected)
+
+        assert H.compare(H1 @ H2)
+
+
+class TestHamiltonianEvaluation:
+    """Test the usage of a Hamiltonian as an observable"""
+
+    @pytest.mark.parametrize("coeffs, param, interface", COEFFS_PARAM_INTERFACE)
+    def test_vqe_forward_different_coeff_types(self, coeffs, param, interface):
+        """Check that manually splitting a Hamiltonian expectation has the same
+        result as passing the Hamiltonian as an observable"""
+        dev = qml.device("default.qubit", wires=2)
+        H = qml.Hamiltonian(coeffs, [qml.PauliX(0), qml.PauliZ(0)])
+
+        @qml.qnode(dev, interface=interface)
+        def circuit():
+            qml.RX(param, wires=0)
+            qml.RY(param, wires=0)
+            return qml.expval(H)
+
+        @qml.qnode(dev, interface=interface)
+        def circuit1():
+            qml.RX(param, wires=0)
+            qml.RY(param, wires=0)
+            return qml.expval(qml.PauliX(0))
+
+        @qml.qnode(dev, interface=interface)
+        def circuit2():
+            qml.RX(param, wires=0)
+            qml.RY(param, wires=0)
+            return qml.expval(qml.PauliZ(0))
+
+        res = circuit()
+        res_expected = coeffs[0] * circuit1() + coeffs[1] * circuit2()
+        assert np.isclose(res, res_expected)
+
+    def test_simplify_reduces_tape_parameters(self):
+        """Test that simplifying a Hamiltonian reduces the number of parameters on a tape"""
+        dev = qml.device("default.qubit", wires=2)
+
+        @qml.qnode(dev)
+        def circuit():
+            qml.RY(0.1, wires=0)
+            return qml.expval(
+                qml.Hamiltonian([1.0, 2.0], [qml.PauliX(1), qml.PauliX(1)], simplify=True)
+            )
+
+        circuit()
+        pars = circuit.qtape.get_parameters(trainable_only=False)
+        # simplify worked and added 1. and 2.
+        assert pars == [0.1, 3.0]
+
+
+class TestHamiltonianDifferentiation:
+    """Test that the Hamiltonian coefficients are differentiable"""
+
+    @pytest.mark.parametrize("simplify", [True, False])
+    def test_vqe_differentiation_paramshift(self, simplify):
+        """Test the parameter-shift method by comparing the differentiation of linearly combined subcircuits
+        with the differentiation of a Hamiltonian expectation"""
+        coeffs = np.array([-0.05, 0.17])
+        param = np.array(1.7)
+
+        # differentiating a circuit with measurement expval(H)
+        @qml.qnode(dev, diff_method="parameter-shift")
+        def circuit(coeffs, param):
+            qml.RX(param, wires=0)
+            qml.RY(param, wires=0)
+            return qml.expval(
+                qml.Hamiltonian(coeffs, [qml.PauliX(0), qml.PauliZ(0)], simplify=simplify)
+            )
+
+        grad_fn = qml.grad(circuit)
+        grad = grad_fn(coeffs, param)
+
+        # differentiating a cost that combines circuits with
+        # measurements expval(Pauli)
+        half1 = qml.QNode(circuit1, dev, diff_method="parameter-shift")
+        half2 = qml.QNode(circuit2, dev, diff_method="parameter-shift")
+
+        def combine(coeffs, param):
+            return coeffs[0] * half1(param) + coeffs[1] * half2(param)
+
+        grad_fn_expected = qml.grad(combine)
+        grad_expected = grad_fn_expected(coeffs, param)
+
+        assert np.allclose(grad[0], grad_expected[0])
+        assert np.allclose(grad[1], grad_expected[1])
+
+    @pytest.mark.parametrize("simplify", [True, False])
+    def test_vqe_differentiation_autograd(self, simplify):
+        """Test the autograd interface by comparing the differentiation of linearly combined subcircuits
+        with the differentiation of a Hamiltonian expectation"""
+        coeffs = pnp.array([-0.05, 0.17], requires_grad=True)
+        param = pnp.array(1.7, requires_grad=True)
+
+        # differentiating a circuit with measurement expval(H)
+        @qml.qnode(dev, interface="autograd")
+        def circuit(coeffs, param):
+            qml.RX(param, wires=0)
+            qml.RY(param, wires=0)
+            return qml.expval(
+                qml.Hamiltonian(coeffs, [qml.PauliX(0), qml.PauliZ(0)], simplify=simplify)
+            )
+
+        grad_fn = qml.grad(circuit)
+        grad = grad_fn(coeffs, param)
+
+        # differentiating a cost that combines circuits with
+        # measurements expval(Pauli)
+        half1 = qml.QNode(circuit1, dev, interface="autograd")
+        half2 = qml.QNode(circuit2, dev, interface="autograd")
+
+        def combine(coeffs, param):
+            return coeffs[0] * half1(param) + coeffs[1] * half2(param)
+
+        grad_fn_expected = qml.grad(combine)
+        grad_expected = grad_fn_expected(coeffs, param)
+
+        assert np.allclose(grad[0], grad_expected[0])
+        assert np.allclose(grad[1], grad_expected[1])
+
+    @pytest.mark.parametrize("simplify", [True, False])
+    def test_vqe_differentiation_jax(self, simplify):
+        """Test the jax interface by comparing the differentiation of linearly combined subcircuits
+        with the differentiation of a Hamiltonian expectation"""
+
+        jax = pytest.importorskip("jax")
+        jnp = pytest.importorskip("jax.numpy")
+        coeffs = jnp.array([-0.05, 0.17])
+        param = jnp.array(1.7)
+
+        # differentiating a circuit with measurement expval(H)
+        @qml.qnode(dev, interface="jax", diff_method="backprop")
+        def circuit(coeffs, param):
+            qml.RX(param, wires=0)
+            qml.RY(param, wires=0)
+            return qml.expval(
+                qml.Hamiltonian(coeffs, [qml.PauliX(0), qml.PauliZ(0)], simplify=simplify)
+            )
+
+        grad_fn = jax.grad(circuit)
+        grad = grad_fn(coeffs, param)
+
+        # differentiating a cost that combines circuits with
+        # measurements expval(Pauli)
+        half1 = qml.QNode(circuit1, dev, interface="jax", diff_method="backprop")
+        half2 = qml.QNode(circuit2, dev, interface="jax", diff_method="backprop")
+
+        def combine(coeffs, param):
+            return coeffs[0] * half1(param) + coeffs[1] * half2(param)
+
+        grad_fn_expected = jax.grad(combine)
+        grad_expected = grad_fn_expected(coeffs, param)
+
+        assert np.allclose(grad[0], grad_expected[0])
+        assert np.allclose(grad[1], grad_expected[1])
+
+    @pytest.mark.parametrize("simplify", [True, False])
+    def test_vqe_differentiation_torch(self, simplify):
+        """Test the torch interface by comparing the differentiation of linearly combined subcircuits
+        with the differentiation of a Hamiltonian expectation"""
+
+        torch = pytest.importorskip("torch")
+        coeffs = torch.tensor([-0.05, 0.17], requires_grad=True)
+        param = torch.tensor(1.7, requires_grad=True)
+
+        # differentiating a circuit with measurement expval(H)
+        @qml.qnode(dev, interface="torch")
+        def circuit(coeffs, param):
+            qml.RX(param, wires=0)
+            qml.RY(param, wires=0)
+            return qml.expval(
+                qml.Hamiltonian(coeffs, [qml.PauliX(0), qml.PauliZ(0)], simplify=simplify)
+            )
+
+        res = circuit(coeffs, param)
+        res.backward()
+        grad = (coeffs.grad, param.grad)
+
+        # differentiating a cost that combines circuits with
+        # measurements expval(Pauli)
+
+        # we need to create new tensors here
+        coeffs2 = torch.tensor([-0.05, 0.17], requires_grad=True)
+        param2 = torch.tensor(1.7, requires_grad=True)
+
+        half1 = qml.QNode(circuit1, dev, interface="torch")
+        half2 = qml.QNode(circuit2, dev, interface="torch")
+
+        def combine(coeffs, param):
+            return coeffs[0] * half1(param) + coeffs[1] * half2(param)
+
+        res_expected = combine(coeffs2, param2)
+        res_expected.backward()
+        grad_expected = (coeffs2.grad, param2.grad)
+
+        assert np.allclose(grad[0], grad_expected[0])
+        assert np.allclose(grad[1], grad_expected[1])
+
+    @pytest.mark.parametrize("simplify", [True, False])
+    def test_vqe_differentiation_tf(self, simplify):
+        """Test the tf interface by comparing the differentiation of linearly combined subcircuits
+        with the differentiation of a Hamiltonian expectation"""
+
+        tf = pytest.importorskip("tf")
+        coeffs = tf.Variable([-0.05, 0.17], dtype=tf.double)
+        param = tf.Variable(1.7, dtype=tf.double)
+
+        # differentiating a circuit with measurement expval(H)
+        @qml.qnode(dev, interface="tf", diff_method="backprop")
+        def circuit(coeffs, param):
+            qml.RX(param, wires=0)
+            qml.RY(param, wires=0)
+            return qml.expval(
+                qml.Hamiltonian(coeffs, [qml.PauliX(0), qml.PauliZ(0)], simplify=simplify)
+            )
+
+        with tf.GradientTape() as tape:
+            res = circuit(coeffs, param)
+        grad = tape.gradient(res, [coeffs, param])
+
+        # differentiating a cost that combines circuits with
+        # measurements expval(Pauli)
+
+        # we need to create new tensors here
+        coeffs2 = tf.Variable([-0.05, 0.17], dtype=tf.double)
+        param2 = tf.Variable(1.7, dtype=tf.double)
+        half1 = qml.QNode(circuit1, dev, interface="tf", diff_method="backprop")
+        half2 = qml.QNode(circuit2, dev, interface="tf", diff_method="backprop")
+
+        def combine(coeffs, param):
+            return coeffs[0] * half1(param) + coeffs[1] * half2(param)
+
+        with tf.GradientTape() as tape2:
+            res_expected = combine(coeffs2, param2)
+        grad_expected = tape2.gradient(res_expected, [coeffs2, param2])
+
+        assert np.allclose(grad[0], grad_expected[0])
+        assert np.allclose(grad[1], grad_expected[1])
diff --git a/tests/test_qaoa.py b/tests/test_qaoa.py
index 5fb1c04ed..8b904b23d 100644
--- a/tests/test_qaoa.py
+++ b/tests/test_qaoa.py
@@ -54,7 +54,7 @@ for k, v in complete_edge_weight_data.items():
 
 def decompose_hamiltonian(hamiltonian):
 
-    coeffs = hamiltonian.coeffs
+    coeffs = list(qml.math.toarray(hamiltonian.coeffs))
     ops = [i.name for i in hamiltonian.ops]
     wires = [i.wires for i in hamiltonian.ops]
 
@@ -623,7 +623,7 @@ MWC = list(zip(DIGRAPHS, MWC_CONSTRAINED, COST_HAMILTONIANS, MIXER_HAMILTONIANS,
 
 def decompose_hamiltonian(hamiltonian):
 
-    coeffs = hamiltonian.coeffs
+    coeffs = list(qml.math.toarray(hamiltonian.coeffs))
     ops = [i.name for i in hamiltonian.ops]
     wires = [i.wires for i in hamiltonian.ops]
 
@@ -1114,7 +1114,7 @@ class TestCycles:
         ]
         expected_coeffs = [np.log(0.5), np.log(1), np.log(1.5), np.log(2), np.log(2.5), np.log(3)]
 
-        assert expected_coeffs == h.coeffs
+        assert np.allclose(expected_coeffs, h.coeffs)
         assert all([op.wires == exp.wires for op, exp in zip(h.ops, expected_ops)])
         assert all([type(op) is type(exp) for op, exp in zip(h.ops, expected_ops)])
 
@@ -1160,7 +1160,7 @@ class TestCycles:
             np.log(7),
         ]
 
-        assert expected_coeffs == h.coeffs
+        assert np.allclose(expected_coeffs, h.coeffs)
         assert all([op.wires == exp.wires for op, exp in zip(h.ops, expected_ops)])
         assert all([type(op) is type(exp) for op, exp in zip(h.ops, expected_ops)])
 
@@ -1252,7 +1252,7 @@ class TestCycles:
 
         expected_coeffs = [2, 2, -2, -2]
 
-        assert expected_coeffs == h.coeffs
+        assert np.allclose(expected_coeffs, h.coeffs)
         for i, expected_op in enumerate(expected_ops):
             assert str(h.ops[i]) == str(expected_op)
         assert all([op.wires == exp.wires for op, exp in zip(h.ops, expected_ops)])
@@ -1274,7 +1274,7 @@ class TestCycles:
         ]
         expected_coeffs = [4, 2, -2, -2, -2, -2, 2]
 
-        assert expected_coeffs == h.coeffs
+        assert np.allclose(expected_coeffs, h.coeffs)
         for i, expected_op in enumerate(expected_ops):
             assert str(h.ops[i]) == str(expected_op)
         assert all([op.wires == exp.wires for op, exp in zip(h.ops, expected_ops)])
@@ -1290,7 +1290,7 @@ class TestCycles:
         expected_ops = [qml.PauliZ(wires=[0])]
         expected_coeffs = [0]
 
-        assert expected_coeffs == h.coeffs
+        assert np.allclose(expected_coeffs, h.coeffs)
         for i, expected_op in enumerate(expected_ops):
             assert str(h.ops[i]) == str(expected_op)
         assert all([op.wires == exp.wires for op, exp in zip(h.ops, expected_ops)])
@@ -1313,7 +1313,7 @@ class TestCycles:
         ]
         expected_coeffs = [4, -2, -2, 2, 2, -2, -2]
 
-        assert expected_coeffs == h.coeffs
+        assert np.allclose(expected_coeffs, h.coeffs)
         for i, expected_op in enumerate(expected_ops):
             assert str(h.ops[i]) == str(expected_op)
         assert all([op.wires == exp.wires for op, exp in zip(h.ops, expected_ops)])
diff --git a/tests/test_vqe.py b/tests/test_vqe.py
index 6fe2d6f99..8754804ab 100644
--- a/tests/test_vqe.py
+++ b/tests/test_vqe.py
@@ -649,7 +649,8 @@ class TestHamiltonian:
         """Tests that the Hamiltonian object is created with
         the correct attributes"""
         H = qml.vqe.Hamiltonian(coeffs, ops)
-        assert H.terms == (list(coeffs), list(ops))
+        assert np.allclose(H.terms[0], coeffs)
+        assert H.terms[1] == list(ops)
 
     @pytest.mark.parametrize("coeffs, ops", invalid_hamiltonians)
     def test_hamiltonian_invalid_init_exception(self, coeffs, ops):
@@ -658,15 +659,6 @@ class TestHamiltonian:
         with pytest.raises(ValueError, match="number of coefficients and operators does not match"):
             H = qml.vqe.Hamiltonian(coeffs, ops)
 
-    @pytest.mark.parametrize("coeffs", [[0.2, -1j], [0.5j, 2 - 1j]])
-    def test_hamiltonian_complex(self, coeffs):
-        """Tests that an exception is raised when
-        a complex Hamiltonian is given"""
-        obs = [qml.PauliX(0), qml.PauliZ(1)]
-
-        with pytest.raises(ValueError, match="coefficients are not real-valued"):
-            H = qml.vqe.Hamiltonian(coeffs, obs)
-
     @pytest.mark.parametrize(
         "obs", [[qml.PauliX(0), qml.CNOT(wires=[0, 1])], [qml.PauliZ, qml.PauliZ(0)]]
     )
@@ -843,13 +835,15 @@ class TestHamiltonian:
         with qml.tape.QuantumTape() as tape:
             qml.Hadamard(wires=1)
             qml.PauliX(wires=0)
-            H.queue()
+            qml.expval(H)
 
         assert np.all([q1.compare(q2) for q1, q2 in zip(tape.queue, queue)])
 
         # Inside of tape
 
         queue = [
+            qml.Hadamard(wires=1),
+            qml.PauliX(wires=0),
             qml.PauliX(1),
             qml.PauliZ(0),
             qml.PauliZ(2),
@@ -858,16 +852,16 @@ class TestHamiltonian:
             qml.Hamiltonian(
                 [1, 3, 1], [qml.PauliX(1), qml.PauliZ(0) @ qml.PauliZ(2), qml.PauliZ(1)]
             ),
-            qml.Hadamard(wires=1),
-            qml.PauliX(wires=0),
         ]
 
         with qml.tape.QuantumTape() as tape:
-            H = qml.Hamiltonian(
-                [1, 3, 1], [qml.PauliX(1), qml.PauliZ(0) @ qml.PauliZ(2), qml.PauliZ(1)]
-            )
             qml.Hadamard(wires=1)
             qml.PauliX(wires=0)
+            qml.expval(
+                qml.Hamiltonian(
+                    [1, 3, 1], [qml.PauliX(1), qml.PauliZ(0) @ qml.PauliZ(2), qml.PauliZ(1)]
+                )
+            )
 
         assert np.all([q1.compare(q2) for q1, q2 in zip(tape.queue, queue)])
 
@@ -1323,14 +1317,14 @@ class TestNewVQE:
             return qml.expval(H1)
 
         res = qml.draw(circuit1)()
-        expected = " 0: ──H──╭Hamiltonian──┤  \n" + " 2: ─────╰Hamiltonian──┤  \n"
+        expected = " 0: ──H──╭┤ ⟨Hamiltonian(1, 1, 1)⟩ \n" + " 2: ─────╰┤ ⟨Hamiltonian(1, 1, 1)⟩ \n"
         assert res == expected
 
     def test_error_multiple_expvals(self):
         """Tests that error is thrown if more than one expval is evaluated."""
         observables = [qml.PauliZ(0), qml.PauliY(0), qml.PauliZ(1)]
         coeffs = [1.0] * len(observables)
-        dev = qml.device("default.qubit", wires=3)
+        dev = qml.device("default.qubit", wires=4)
         H = qml.Hamiltonian(coeffs, observables)
         w = qml.init.strong_ent_layers_uniform(2, 4, seed=1967)
 
@@ -1339,31 +1333,34 @@ class TestNewVQE:
             qml.templates.StronglyEntanglingLayers(w, wires=range(4))
             return qml.expval(H), qml.expval(qml.PauliX(3))
 
-        with pytest.raises(ValueError, match="At the moment"):
+        with pytest.raises(
+            ValueError, match="Only a single expectation of a Hamiltonian observable"
+        ):
             circuit()
 
     def test_error_non_expval_measurement(self):
         """Tests that error is thrown if sample() or var() is used."""
         observables = [qml.PauliZ(0), qml.PauliY(0), qml.PauliZ(1)]
         coeffs = [1.0] * len(observables)
-        dev = qml.device("default.qubit", wires=3)
+        dev = qml.device("default.qubit", wires=4)
         H = qml.Hamiltonian(coeffs, observables)
-        w = qml.init.strong_ent_layers_uniform(2, 4, seed=1967)
 
         @qml.qnode(dev)
         def circuit():
-            qml.templates.StronglyEntanglingLayers(w, wires=range(4))
             return qml.sample(H)
 
-        with pytest.raises(qml.QuantumFunctionError, match="Hamiltonian is not an observable"):
+        with pytest.raises(
+            ValueError, match="Only a single expectation of a Hamiltonian observable "
+        ):
             circuit()
 
         @qml.qnode(dev)
         def circuit():
-            qml.templates.StronglyEntanglingLayers(w, wires=range(4))
             return qml.var(H)
 
-        with pytest.raises(qml.QuantumFunctionError, match="Hamiltonian is not an observable"):
+        with pytest.raises(
+            ValueError, match="Only a single expectation of a Hamiltonian observable"
+        ):
             circuit()
 
     @pytest.mark.parametrize("diff_method", ["parameter-shift", "best"])
diff --git a/tests/transforms/test_hamiltonian_expand.py b/tests/transforms/test_hamiltonian_expand.py
index a2533a8d2..e1cde7c5d 100644
--- a/tests/transforms/test_hamiltonian_expand.py
+++ b/tests/transforms/test_hamiltonian_expand.py
@@ -34,7 +34,6 @@ with pennylane.tape.QuantumTape() as tape2:
     qml.Hadamard(1)
     qml.PauliZ(1)
     qml.PauliX(2)
-
     H2 = qml.Hamiltonian(
         [1, 3, -2, 1, 1],
         [
@@ -130,7 +129,18 @@ class TestHamiltonianExpval:
         H = qml.Hamiltonian(
             [-0.2, 0.5, 1], [qml.PauliX(1), qml.PauliZ(1) @ qml.PauliY(2), qml.PauliZ(0)]
         )
-        var = np.array([0.1, 0.67, 0.3, 0.4, -0.5, 0.7])
+
+        var = [
+            np.array(0.1),
+            np.array(0.67),
+            np.array(0.3),
+            np.array(0.4),
+            np.array(-0.5),
+            np.array(0.7),
+            np.array(0.4),
+            np.array(-0.5),
+            np.array(0.7),
+        ]
         output = 0.42294409781940356
         output2 = [
             9.68883500e-02,
@@ -139,6 +149,9 @@ class TestHamiltonianExpval:
             -1.94289029e-09,
             3.50307411e-01,
             -3.41123470e-01,
+            0.0,  # these three are the Hamiltonian parameters
+            0.0,
+            0.0,
         ]
 
         with qml.tape.JacobianTape() as tape:
@@ -161,7 +174,9 @@ class TestHamiltonianExpval:
             return fn(res)
 
         assert np.isclose(cost(var), output)
-        assert np.allclose(qml.grad(cost)(var), output2)
+        grad = qml.grad(cost)(var)
+        for g, o in zip(grad, output2):
+            assert np.allclose(g, o)
 
     def test_hamiltonian_dif_tensorflow(self):
         """Tests that the hamiltonian_expand tape transform is differentiable with the Tensorflow interface"""