Repository: PennyLaneAI/pennylane
Pull request #5161 (2024-02-06 13:49:51)

[DLA 1] Add qml.dla.VSpace class for keeping track of linearly independent vectors
This class enables easy implementation for `qml.dla.lie_closure` realized in the follow-up PR.
The class VSpace is currently not (yet) user-facing and acts as a utility to lie_closure. We could consider making it user-facing though. Open to discussing this! The main reason I am holding back is that the UI is experimental and some things may be subject to change / improvement, so I dont want to lock in just yet.

The main point of this class is to efficiently keep track of linearly independent operators that span a vector space, specifically for PauliSentences.

This is done by storing and processing ``M``, which
is a dictionary-of-keys (DOK) style sparse representation of the set of basis vectors. You can
think of it as the numpy-equivalent of a PauliSentence: each ``PauliWord`` (key of ``PauliSentence``)
represents one row of ``M`` with the coefficient (value of ``PauliSentence``).
For example the set of 3 linearly independent generators ``X(0) + X(1), X(0) + X(2), X(0) + 0.5 * Y(0)``
can be represented as

```
[
    [1, 1, 1],
    [1, 0, 0],
    [0, 1, 0],
    [0, 0, 0.5]
]
```

where each column represents one sentence, and each row represents the coefficient of the respective word in the sentence.
To make sense of this representation one additionally needs to keep track of the mapping between keys and rows. In this case we have

```
pw_to_idx = {
    X(0) : 0,
    X(1) : 1,
    X(2) : 2,
    Y(0) : 3
}
``` 

where we have set the numbering based on appearance in the list of generators. This mapping is in general not unique.

ToDo

- [x] VSpace class init
- [x] VSpace.add() method
- [x] test init
- [x] test VSpace.add()
- [x] add cheaper check for proportionality before performing rank computation
- [x] add test for util accordingly
- [x] increase test coverage (i.e. add examples)
- [x] check edge cases
- [x] changelog (no changelog since this is private functionality)
- [x] Fix equality check
- [x] Overhaul `VSpace.add` logic (inplace vs not)

Based off @dwierichs 's [prototype](https://github.com/PennyLaneAI/pennylane/pull/4852)
Follow-up PR adding lie_closure https://github.com/PennyLaneAI/pennylane/pull/5169
[sc-51430]

https://github.com/PennyLaneAI/pennylane/pull/5161 (DLA1) https://github.com/PennyLaneAI/pennylane/pull/5169 (DLA2) https://github.com/PennyLaneAI/pennylane/pull/5406 (DLA3)

codecov-commenter: ## [Codecov](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5161?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
Attention: `5 lines` in your changes are missing coverage. Please review.
> Comparison is base [(`0428922`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/0428922017c61c6bc0405582653eae2a9ecbd8cb?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) 99.68% compared to head [(`b3121ca`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5161?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) 99.66%.

| [Files](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5161?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Patch % | Lines |
|---|---|---|
| [pennylane/dla/lie\_closure.py](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5161?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2RsYS9saWVfY2xvc3VyZS5weQ==) | 92.85% | [5 Missing :warning: ](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5161?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) |

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##           master    #5161      +/-   ##
==========================================
- Coverage   99.68%   99.66%   -0.02%     
==========================================
  Files         394      396       +2     
  Lines       36298    36098     -200     
==========================================
- Hits        36183    35977     -206     
- Misses        115      121       +6     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5161?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).

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.
Qottmann: > One thing I was wondering, do we want to preserve the original PauliSentences in the basis, or can we build an equivalent orthonormal basis? Some operations may be simplified in the latter case.

I was wondering about this. If we allow for this we could also be more efficient adding multiple operators at once by looking at the combined subspace and taking the maximal linear independent set of those as the new basis after addition.

On the technical side this requires some book-keeping on the pw_to_idx mapping and updating the vspace.basis accordingly.

On the UI/UX side I could see it be confusing if your operators suddenly change, in particular when you are looking at the Lie algebra generated by some specific Hamiltonian that you are looking at. That is why I am a bit reluctant to change the generators/basis.
Qottmann: Non-urgent, but if you find the time to re-review this that would be great.
I have been extensively using this and the follow-up PR with the actual lie closure (https://github.com/PennyLaneAI/pennylane/pull/5169)
Qottmann: > @Qottmann, PauliVSpace was previously showing up in the docs but now it isn't. Is that intentional?

I am still seeing it listed in `pauli` under _Classes_:
![image](https://github.com/PennyLaneAI/pennylane/assets/43949391/e5597b4f-bc81-4af3-acff-9d2a7d16943a)


Patch:
diff --git a/pennylane/pauli/__init__.py b/pennylane/pauli/__init__.py
index ff08ec228..26ab2ea89 100644
--- a/pennylane/pauli/__init__.py
+++ b/pennylane/pauli/__init__.py
@@ -49,3 +49,5 @@ from .grouping import (
     optimize_measurements,
     graph_colouring,
 )
+
+from .dla import PauliVSpace
diff --git a/pennylane/pauli/dla/__init__.py b/pennylane/pauli/dla/__init__.py
new file mode 100644
index 000000000..204de97fa
--- /dev/null
+++ b/pennylane/pauli/dla/__init__.py
@@ -0,0 +1,18 @@
+# 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.
+"""
+This subpackage defines functions and classes for dynamical Lie algebra functionality
+"""
+
+from .lie_closure import PauliVSpace
diff --git a/pennylane/pauli/dla/lie_closure.py b/pennylane/pauli/dla/lie_closure.py
new file mode 100644
index 000000000..6215d3867
--- /dev/null
+++ b/pennylane/pauli/dla/lie_closure.py
@@ -0,0 +1,289 @@
+# Copyright 2024 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.
+"""A function to compute the Lie closure of a set of operators"""
+# pylint: disable=too-many-arguments
+from functools import reduce
+
+from copy import copy
+import numpy as np
+
+import pennylane as qml
+
+from pennylane.operation import Operator
+from ..pauli_arithmetic import PauliSentence
+
+
+class PauliVSpace:
+    """
+    Class representing the linearly independent basis of a vector space.
+
+    The main purpose of this class is to store and process ``M``, which
+    is a dictionary-of-keys (DOK) style sparse representation of the set of basis vectors. You can
+    think of it as the numpy-equivalent of a PauliSentence: each :class:`~pennylane.pauli.PauliWord` (key of :class:`~pennylane.pauli.PauliSentence`)
+    represents one row of ``M`` with the coefficient (value of :class:`~pennylane.pauli.PauliSentence`).
+    For example the set of 3 linearly independent generators ``X(0) + X(1), X(0) + X(2), X(0) + 0.5 * Y(0)``
+    can be represented as
+
+    .. code-block:: python3
+
+        [
+            [1, 1, 1],
+            [1, 0, 0],
+            [0, 1, 0],
+            [0, 0, 0.5]
+        ]
+
+    where each column represents one sentence, and each row represents the coefficient of the respective word in the sentence.
+    To make sense of this representation one additionally needs to keep track of the mapping between keys and rows. In this case we have
+
+    .. code-block:: python3
+
+        pw_to_idx = {
+            X(0) : 0,
+            X(1) : 1,
+            X(2) : 2,
+            Y(0) : 3
+        }
+
+    where we have set the numbering based on appearance in the list of generators. This mapping is in general not unique.
+
+    Args:
+        generators (Iterable[Union[PauliWord, PauliSentence, Operator]]): Operators that span the vector space.
+        dtype (type): ``dtype`` of the underlying DOK sparse matrix ``M``. Default is ``float``.
+
+    **Example**
+
+    Take a linearly dependent set of operators and span the PauliVSpace.
+
+    .. code-block:: python3
+
+        ops = [
+            X(0) @ X(1) + Y(0) @ Y(1),
+            X(0) @ X(1),
+            Y(0) @ Y(1)
+        ]
+
+        vspace = PauliVSpace(ops)
+
+    It automatically detects that the third operator is linearly dependent on the former two, so it does not add the third operator to the basis.
+
+    >>> vspace.basis
+    [1.0 * X(0) @ X(1)
+     + 1.0 * Y(0) @ Y(1),
+     1.0 * X(0) @ X(1)]
+
+    We can also retrospectively add operators.
+
+    >>> vspace.add(qml.X(0))
+    [1.0 * X(0) @ X(1)
+     + 1.0 * Y(0) @ Y(1),
+     1.0 * X(0) @ X(1),
+     1.0 * X(0)]
+
+    Again, checks of linear independence are always performed. So in the following example no operator is added.
+
+    >>> vspace.add(Y(0) @ Y(1))
+    [1.0 * X(0) @ X(1)
+     + 1.0 * Y(0) @ Y(1),
+     1.0 * X(0) @ X(1),
+     1.0 * X(0)]
+    """
+
+    def __init__(self, generators, dtype=float):
+
+        self.dtype = dtype
+
+        if any(not isinstance(g, PauliSentence) for g in generators):
+            generators = [
+                qml.pauli.pauli_sentence(g) if not isinstance(g, PauliSentence) else g
+                for g in generators
+            ]
+
+        # Get all Pauli words that are present in at least one Pauli sentence
+        all_pws = list(reduce(set.__or__, [set(ps.keys()) for ps in generators]))
+        num_pw = len(all_pws)
+        # Create a dictionary mapping from PauliWord to row index
+        self._pw_to_idx = {pw: i for i, pw in enumerate(all_pws)}
+
+        # Initialize PauliVSpace properties trivially
+        self._basis = []
+        rank = 0
+
+        self._M = np.zeros((num_pw, rank), dtype=self.dtype)
+        self._rank = rank
+        self._num_pw = num_pw
+
+        # Add all generators that are linearly independent
+        self.add(generators)
+
+    @property
+    def basis(self):
+        """List of basis operators of PauliVSpace"""
+        return self._basis
+
+    def __len__(self):
+        return len(self.basis)
+
+    def add(self, other, tol=1e-15):
+        r"""Adding Pauli sentences if they are linearly independent.
+
+        Args:
+            other (List[:class:`~.PauliWord`, :class:`~.PauliSentence`, :class:`~.Operator`]): List of candidate operators to add to the ``PauliVSpace``, if they are linearly independent.
+            tol (float): Numerical tolerance for linear independence check. Defaults to ``1e-15``.
+
+        Returns:
+            List: New basis vectors after adding the linearly independent ones from ``other``.
+
+        **Example**
+
+        We can generate a ``PauliVSpace`` and add a linearly independent operator to its basis.
+
+        >>> ops = [X(0), X(1)]
+        >>> vspace = qml.pauli.PauliVSpace(ops)
+        >>> vspace.add(Y(0))
+        >>> vspace
+        [1.0 * X(0), 1.0 * X(1), 1.0 * Y(0)]
+
+        We can add a list of operators at once. Only those that are linearly dependent with the current ``PauliVSpace`` are added.
+
+        >>> vspace.add([Z(0), X(0)])
+        [1.0 * X(0), 1.0 * X(1), 1.0 * Y(0), 1.0 * Z(0)]
+
+        """
+        if isinstance(other, (qml.pauli.PauliWord, qml.pauli.PauliSentence, Operator)):
+            other = [other]
+
+        other = [
+            qml.pauli.pauli_sentence(op) if not isinstance(op, qml.pauli.PauliSentence) else op
+            for op in other
+        ]
+
+        for ps in other:
+            # TODO: Potential speed-up by computing the maximal linear independent set for all current basis vectors + other, essentially algorithm1 in https://arxiv.org/abs/1012.5256
+            self._M, self._pw_to_idx, self._rank, self._num_pw, is_independent = (
+                self._check_independence(
+                    self._M, ps, self._pw_to_idx, self._rank, self._num_pw, tol
+                )
+            )
+            if is_independent:
+                self._basis.append(ps)
+        return self._basis
+
+    def is_independent(self, pauli_sentence, tol=1e-15):
+        r"""Check if the ``pauli_sentence`` is linearly independent of the basis of ``PauliVSpace``.
+
+        Args:
+            pauli_sentence (`~.PauliSentence`): Candidate Pauli sentence to check against the ``PauliVSpace`` basis for linear independence.
+            tol (float): Numerical tolerance for linear independence check. Defaults to ``1e-15``.
+
+        Returns:
+            bool: whether ``pauli_sentence`` was linearly independent
+
+        **Example**
+
+        >>> ops = [X(0), X(1)]
+        >>> vspace = PauliVSpace([op.pauli_rep for op in ops])
+        >>> vspace.is_independent(X(0).pauli_rep)
+        False
+        >>> vspace.is_independent(Y(0).pauli_rep)
+        True
+
+        """
+        _, _, _, _, is_independent = self._check_independence(
+            self._M, pauli_sentence, self._pw_to_idx, self._rank, self._num_pw, tol
+        )
+        return is_independent
+
+    @staticmethod
+    def _check_independence(M, pauli_sentence, pw_to_idx, rank, num_pw, tol=1e-15):
+        r"""
+        Checks if :class:`~PauliSentence` ``pauli_sentence`` is linearly independent and provides the updated class attributes in case the vector is added.
+
+        This is done in the following way: ``M`` (see description in class) is extended by ``pauli_sentence``.
+        If the added operator has a PauliWord (key) that is new to ``pw_to_idx``, then we have to add a new row
+        and already know that it has to be linearly independent.
+        If it contains the same PauliWords, we need to compute the new rank and compare it with the old rank.
+        If the rank is the same, the operator is linearly dependent and not added. Else, the rank is incrased by 1
+        and the extended M becomes our new M.
+
+        Args:
+            M (ndarray): coefficient matrix for current LIS
+            pauli_sentence (`~.PauliSentence`): Pauli sentence for which to add a column if independent
+            pw_to_idx (dict): map from :class:`~pennylane.pauli.PauliWord` to row index in ``M``
+            rank (int): current rank of ``M``, equal to its number of columns
+            num_pw (int): current number of :class:`~pennylane.pauli.PauliWord`\ s, equal to the number of rows in ``M``
+            tol (float): Numerical tolerance for linear independence check.
+
+        Returns:
+            ndarray: updated coefficient matrix for the LIS
+            dict: updated map from :class:`~pennylane.pauli.PauliWord` to row index in ``M``. Includes new :class:`~pennylane.pauli.PauliWord` keys
+                from the input ``pauli_sentence`` if it was linearly independent
+            int: updated rank/number of columns of ``M``
+            int: updated number of :class:`~pennylane.pauli.PauliWord`\ s/number of rows of ``M``
+            bool: whether ``pauli_sentence`` was linearly independent and whether its column was added to ``M``
+        """
+        new_pws = [pw for pw in pauli_sentence.keys() if pw not in pw_to_idx]
+        new_num_pw = num_pw + len(new_pws)
+
+        if new_num_pw < rank + 1:
+            # Can't span rank+1 independent vectors in fewer than rank+1 dimensions
+            # The input PauliSentence must have been linearly dependent
+            return M, pw_to_idx, rank, num_pw, False
+
+        M = np.pad(M, ((0, new_num_pw - num_pw), (0, 1)))
+
+        # If there are new PauliWords (i.e. new basis vectors), the candidate vector must be linearly independent
+        if new_num_pw > num_pw:
+            new_pw_to_idx = copy(pw_to_idx)
+            for i, pw in enumerate(new_pws, start=num_pw):
+                new_pw_to_idx[pw] = i
+            # Add new PauliSentence entries to matrix
+            for pw, value in pauli_sentence.items():
+                M[new_pw_to_idx[pw], rank] = value
+
+            return M, new_pw_to_idx, rank + 1, new_num_pw, True
+
+        # Add new PauliSentence entries to matrix
+        for pw, value in pauli_sentence.items():
+            M[pw_to_idx[pw], rank] = value
+
+        # Check if new vector is linearly dependent on the current basis
+        v = M[:, -1].copy()  # remove copy to normalize M
+        v /= np.linalg.norm(v)
+        A = M[:, :-1]
+        v = v - A @ qml.math.linalg.inv(qml.math.conj(A.T) @ A) @ qml.math.conj(A).T @ v
+
+        if np.linalg.norm(v) > tol:
+            return M, pw_to_idx, rank + 1, new_num_pw, True
+
+        return M[:num_pw, :rank], pw_to_idx, rank, num_pw, False
+
+    def __repr__(self):
+        return str(self.basis)
+
+    def __eq__(self, other):
+        """
+        Two PauliVSpaces are equivalent when they span the same dimensional space.
+        This is checked here by having matching PauliWord keys in the sparse DOK representation and having the same rank.
+        """
+        if not self._num_pw == other._num_pw:
+            return False
+        if not set(self._pw_to_idx.keys()) == set(other._pw_to_idx.keys()):
+            return False
+
+        rank1 = np.linalg.matrix_rank(self._M)
+        rank2 = np.linalg.matrix_rank(other._M)
+        rank3 = np.linalg.matrix_rank(np.concatenate([self._M, other._M], axis=1))
+
+        return rank1 == rank2 and rank2 == rank3


Test patch:
diff --git a/tests/pauli/dla/test_lie_closure.py b/tests/pauli/dla/test_lie_closure.py
new file mode 100644
index 000000000..528c2e0c7
--- /dev/null
+++ b/tests/pauli/dla/test_lie_closure.py
@@ -0,0 +1,285 @@
+# Copyright 2024 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 pennylane/dla/lie_closure.py functionality"""
+# pylint: disable=too-few-public-methods, protected-access
+from copy import copy
+import pytest
+import numpy as np
+
+import pennylane as qml
+
+from pennylane.pauli import PauliWord, PauliSentence, PauliVSpace
+
+ops1 = [
+    PauliSentence({PauliWord({0: "X", 1: "X"}): 1.0, PauliWord({0: "Y", 1: "Y"}): 1.0}),
+    PauliSentence(
+        {
+            PauliWord({0: "X", 1: "X"}): 1.0,
+        }
+    ),
+    PauliSentence(
+        {
+            PauliWord({0: "Y", 1: "Y"}): 2.0,
+        }
+    ),
+]
+
+ops2 = [
+    PauliSentence({PauliWord({0: "X", 1: "X"}): 1.0, PauliWord({0: "Y", 1: "Y"}): 1.0}),
+    PauliSentence(
+        {
+            PauliWord({0: "X", 1: "X"}): 1.0,
+        }
+    ),
+]
+
+ops2plusY10 = ops2 + [PauliSentence({PauliWord({10: "Y"}): 1.0})]
+
+
+class TestPauliVSpace:
+    """Unit and integration tests for PauliVSpace class"""
+
+    def test_init(self):
+        """Unit tests for initialization"""
+        vspace = PauliVSpace(ops1)
+
+        assert all(isinstance(op, PauliSentence) for op in vspace.basis)
+        assert np.allclose(vspace._M, [[1.0, 1.0], [1.0, 0.0]]) or np.allclose(
+            vspace._M, [[1.0, 0.0], [1.0, 1.0]]
+        )  # the ordering is random as it is taken from a dictionary that has no natural ordering
+        assert vspace.basis == ops1[:-1]
+        assert vspace._rank == 2
+        assert vspace._num_pw == 2
+        assert len(vspace._pw_to_idx) == 2
+
+    @pytest.mark.parametrize("dtype", [float, complex])
+    def test_dtype(self, dtype):
+        vspace = PauliVSpace(ops1, dtype=dtype)
+        assert vspace._M.dtype == dtype
+
+    def test_init_with_ops(self):
+        """Test that initialization with PennyLane operators, PauliWord and PauliSentence works"""
+        ops = [qml.X(0), PauliWord({1: "X"})]
+        vspace = qml.pauli.PauliVSpace(ops)
+        vspace.add(qml.Y(0))
+
+        true_res = qml.pauli.PauliVSpace([qml.X(0), qml.X(1), qml.Y(0)])
+        assert vspace == true_res
+
+    ADD_LINEAR_INDEPENDENT = (
+        (ops2, PauliWord({10: "Y"}), ops2plusY10),
+        (ops2, PauliSentence({PauliWord({10: "Y"}): 1.0}), ops2plusY10),
+        (ops2, qml.PauliY(10), ops2plusY10),
+    )
+
+    def test_repr(self):
+        """Test the string representation of the PauliVSpace instance"""
+        assert (
+            repr(PauliVSpace(ops1)) == "[1.0 * X(0) @ X(1)\n+ 1.0 * Y(0) @ Y(1), 1.0 * X(0) @ X(1)]"
+        )
+
+    @pytest.mark.parametrize("li_tol", [1e-15, 1e-14])
+    @pytest.mark.parametrize("ops, op, true_new_basis", ADD_LINEAR_INDEPENDENT)
+    def test_add_linear_independent(self, ops, op, true_new_basis, li_tol):
+        """Test that adding new (linearly independent) operators works as expected"""
+        vspace = PauliVSpace(ops)
+        new_basis = vspace.add(op, li_tol)
+        assert new_basis == true_new_basis
+
+    ADD_LINEAR_DEPENDENT = (
+        (ops2, PauliWord({0: "Y", 1: "Y"}), ops2),
+        (ops2, PauliSentence({PauliWord({0: "Y", 1: "Y"}): 1.0}), ops2),
+        (ops2, qml.PauliY(0) @ qml.PauliY(1), ops2),
+        (ops2, 0.5 * ops2[0], ops2),
+        (ops2, 0.5 * ops2[1], ops2),
+    )
+
+    @pytest.mark.parametrize("ops, op, true_new_basis", ADD_LINEAR_DEPENDENT)
+    def test_add_lin_dependent(self, ops, op, true_new_basis):
+        """Test that adding linearly dependent operators works as expected"""
+        vspace = PauliVSpace(ops)
+        new_basis = vspace.add(op)
+        assert new_basis == true_new_basis
+
+    @pytest.mark.parametrize("ops, op, true_new_basis", ADD_LINEAR_INDEPENDENT)
+    def test_len(self, ops, op, true_new_basis):
+        """Test the length of the PauliVSpace instance with inplace addition to the basis"""
+        vspace = PauliVSpace(ops)
+        len_before_adding = len(vspace)
+        len_basis_before_adding = len(vspace.basis)
+
+        _ = vspace.add(op)
+        len_after_adding = len(vspace)
+        assert len_after_adding != len_before_adding
+        assert len_after_adding == len(true_new_basis)
+        assert len_before_adding == len_basis_before_adding
+
+    def test_eq_True(self):
+        """Test that equivalent vspaces are correctly determined"""
+        gens1 = [
+            PauliSentence({PauliWord({0: "X"}): 1.0}),
+            PauliSentence({PauliWord({0: "Y"}): 1.0}),
+            PauliSentence({PauliWord({0: "Z"}): 1.0}),
+        ]
+
+        gens2 = [
+            PauliSentence({PauliWord({0: "Z"}): 1.0, PauliWord({0: "Y"}): 1.0}),
+            PauliSentence({PauliWord({0: "X"}): 1.0, PauliWord({0: "Z"}): 1.0}),
+            PauliSentence({PauliWord({0: "Y"}): 1.0}),
+        ]
+
+        vspace1 = PauliVSpace(gens1)
+        vspace2 = PauliVSpace(gens2)
+        assert vspace1 == vspace2
+
+    def test_eq_False0(self):
+        """Test that non-equivalent vspaces are correctly determined"""
+        # Different _num_pw
+        gens1 = [
+            PauliSentence({PauliWord({0: "X"}): 1.0}),
+            PauliSentence({PauliWord({0: "Y"}): 1.0}),
+            PauliSentence({PauliWord({0: "Z"}): 1.0}),
+        ]
+
+        gens2 = [
+            PauliSentence({PauliWord({0: "Z"}): 1.0}),
+            PauliSentence({PauliWord({0: "X"}): 1.0, PauliWord({0: "Z"}): 1.0}),
+        ]
+
+        vspace1 = PauliVSpace(gens1)
+        vspace2 = PauliVSpace(gens2)
+        assert vspace1 != vspace2
+
+    def test_eq_False1(self):
+        """Test that non-equivalent vspaces are correctly determined"""
+        # Same num_pw but acting on different wires
+        gens1 = [
+            PauliSentence({PauliWord({0: "X"}): 1.0}),
+            PauliSentence({PauliWord({0: "Y"}): 1.0}),
+            PauliSentence({PauliWord({0: "Z"}): 1.0}),
+        ]
+
+        gens2 = [
+            PauliSentence({PauliWord({1: "X"}): 1.0}),
+            PauliSentence({PauliWord({0: "Y"}): 1.0}),
+            PauliSentence({PauliWord({0: "Z"}): 1.0}),
+        ]
+
+        vspace1 = PauliVSpace(gens1)
+        vspace2 = PauliVSpace(gens2)
+        assert vspace1 != vspace2
+
+    def test_eq_False2(self):
+        """Test that non-equivalent vspaces are correctly determined"""
+        # Same num_pw but acting on different wires
+        gens1 = [
+            PauliSentence({PauliWord({0: "X"}): 1.0}),
+            PauliSentence({PauliWord({0: "Y"}): 1.0}),
+            PauliSentence({PauliWord({0: "Z"}): 1.0}),
+        ]
+
+        gens2 = [
+            PauliSentence({PauliWord({1: "Z"}): 1.0, PauliWord({0: "Y"}): 1.0}),
+            PauliSentence({PauliWord({0: "X"}): 1.0, PauliWord({1: "Z"}): 1.0}),
+            PauliSentence({PauliWord({0: "Y"}): 1.0}),
+        ]
+
+        vspace1 = PauliVSpace(gens1)
+        vspace2 = PauliVSpace(gens2)
+        assert vspace1 != vspace2
+
+    def test_eq_False3(self):
+        """Test that non-equivalent vspaces are correctly determined"""
+        # Same num_pw, even same pws, but not spanning the same space
+        # vector equivalent of ((1,1,0), (0, 0, 1)) and ((1,0,0), (0,1,0), (0,0,1))
+        # I.e. different rank
+        gens1 = [
+            PauliSentence({PauliWord({0: "X"}): 1.0}),
+            PauliSentence({PauliWord({0: "Y"}): 1.0}),
+            PauliSentence({PauliWord({0: "Z"}): 1.0}),
+        ]
+
+        gens2 = [
+            PauliSentence({PauliWord({1: "X"}): 1.0, PauliWord({0: "Y"}): 1.0}),
+            PauliSentence({PauliWord({0: "Z"}): 1.0}),
+        ]
+
+        vspace1 = PauliVSpace(gens1)
+        vspace2 = PauliVSpace(gens2)
+        assert vspace1 != vspace2
+
+    def test_eq_False4(self):
+        """Test case where both vspaces have the same rank, same PauliWords but span different spaces"""
+        v1 = PauliVSpace(
+            [
+                PauliSentence({PauliWord({0: "X"}): 1.0, PauliWord({0: "Y"}): 1.0}),
+                PauliSentence({PauliWord({0: "Z"}): 1.0, PauliWord({0: "Y"}): 1.0}),
+            ]
+        )
+        v2 = PauliVSpace(
+            [
+                PauliSentence({PauliWord({0: "X"}): 1.0, PauliWord({0: "Z"}): 1.0}),
+                PauliSentence({PauliWord({0: "Y"}): 1.0}),
+            ]
+        )
+        assert v1 != v2
+
+    IS_INDEPENDENT_TEST = (
+        (
+            [
+                PauliSentence({PauliWord({0: "X"}): 1.0}),
+                PauliSentence({PauliWord({0: "Z"}): 1.0}),
+            ],
+            PauliSentence({PauliWord({0: "Y"}): 1.0}),
+            True,
+        ),
+        (
+            [
+                PauliSentence({PauliWord({0: "X"}): 1.0}),
+                PauliSentence({PauliWord({0: "Z"}): 1.0}),
+            ],
+            PauliSentence({PauliWord({0: "Z"}): 1.0}),
+            False,
+        ),
+        (
+            [
+                PauliSentence({PauliWord({0: "X"}): 1.0, PauliWord({1: "Y"}): 1.0}),
+                PauliSentence({PauliWord({0: "Z", 1: "Z"}): 1.0}),
+            ],
+            PauliSentence({PauliWord({0: "Z"}): 1.0}),
+            True,
+        ),
+        (
+            [
+                PauliSentence({PauliWord({0: "X"}): 1.0, PauliWord({1: "Y"}): 1.0}),
+                PauliSentence({PauliWord({0: "Z", 1: "Z"}): 1.0}),
+            ],
+            PauliSentence({PauliWord({0: "Z", 1: "Z"}): 1.0}),
+            False,
+        ),
+    )
+
+    @pytest.mark.parametrize("tol", [1e-15, 1e-8])
+    @pytest.mark.parametrize("ops, op, is_independent_true", IS_INDEPENDENT_TEST)
+    def test_is_independent(self, ops, op, is_independent_true, tol):
+        """Test the `is_independent` method returns correct results and leaves class attributes intact"""
+        v1 = PauliVSpace(ops)
+        vcopy = copy(v1)
+
+        is_independent = v1.is_independent(op, tol=tol)
+        assert is_independent == is_independent_true
+        assert qml.math.allclose(v1._M, vcopy._M)
+        assert v1._pw_to_idx == vcopy._pw_to_idx
+        assert v1._rank == vcopy._rank
+        assert v1._num_pw == vcopy._num_pw