Repository: PennyLaneAI/pennylane
Pull request #5627 (2024-05-02 15:25:32)

[Bugfix] Fix numerical stability of `lie_closure`
Noticed that `lie_closure` was not numerically stable after a [user report](https://discuss.pennylane.ai/t/calculation-of-dynamical-lie-algebra/4309/3)

Fixing that by

- [x] using the slightly slower but more stable `qml.math.linalg.solve` over `qml.math.linalg.inv`
- [x] exposing `tol` as kwarg to the user
- [x] setting the default of `tol` higher (two orders of magnitude higher as the machine precision of `_M`'s dtype)

github-actions[bot]: Hello. You may have forgotten to update the changelog!
Please edit [`doc/releases/changelog-dev.md`](/PennyLaneAI/pennylane/blob/master/doc/releases/changelog-dev.md) with:
* A one-to-two sentence description of the change. You may include a small working example for new features.
* A link back to this PR.
* Your name (or GitHub username) in the contributors section.
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5627?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
Attention: Patch coverage is `92.30769%` with `1 lines` in your changes are missing coverage. Please review.
> :exclamation: No coverage uploaded for pull request base (`v0.36.0-rc0@5808548`). [Click here to learn what that means](https://docs.codecov.io/docs/error-reference?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#section-missing-base-commit).

| [Files](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5627?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Patch % | Lines |
|---|---|---|
| [pennylane/pauli/dla/lie\_closure.py](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5627?src=pr&el=tree&filepath=pennylane%2Fpauli%2Fdla%2Flie_closure.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3BhdWxpL2RsYS9saWVfY2xvc3VyZS5weQ==) | 92.30% | [1 Missing :warning: ](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5627?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               @@
##             v0.36.0-rc0    #5627   +/-   ##
==============================================
  Coverage               ?   99.54%           
==============================================
  Files                  ?      412           
  Lines                  ?    38409           
  Branches               ?        0           
==============================================
  Hits                   ?    38234           
  Misses                 ?      175           
  Partials               ?        0           
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5627?dropdown=coverage&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).


Patch:
diff --git a/doc/releases/changelog-0.36.0.md b/doc/releases/changelog-0.36.0.md
index de7721e00..5037bf239 100644
--- a/doc/releases/changelog-0.36.0.md
+++ b/doc/releases/changelog-0.36.0.md
@@ -152,6 +152,7 @@
 * A new `qml.lie_closure` function to compute the Lie closure of a list of operators.
   [(#5161)](https://github.com/PennyLaneAI/pennylane/pull/5161)
   [(#5169)](https://github.com/PennyLaneAI/pennylane/pull/5169)
+  [(#5627)](https://github.com/PennyLaneAI/pennylane/pull/5627)
 
   The Lie closure, pronounced "Lee closure", is a way to compute the so-called dynamical Lie algebra (DLA) of a set of operators.
   For a list of operators `ops = [op1, op2, op3, ..]`, one computes all nested commutators between `ops` until no new operators are generated from commutation.
diff --git a/pennylane/pauli/dla/lie_closure.py b/pennylane/pauli/dla/lie_closure.py
index 5fa03fd78..e57e99b10 100644
--- a/pennylane/pauli/dla/lie_closure.py
+++ b/pennylane/pauli/dla/lie_closure.py
@@ -31,6 +31,7 @@ def lie_closure(
     max_iterations: int = 10000,
     verbose: bool = False,
     pauli: bool = False,
+    tol: float = None,
 ) -> Iterable[Union[PauliWord, PauliSentence, Operator]]:
     r"""Compute the dynamical Lie algebra from a set of generators.
 
@@ -47,6 +48,7 @@ def lie_closure(
         pauli (bool): Indicates whether it is assumed that :class:`~.PauliSentence` or :class:`~.PauliWord` instances are input and returned.
             This can help with performance to avoid unnecessary conversions to :class:`~pennylane.operation.Operator`
             and vice versa. Default is ``False``.
+        tol (float): Numerical tolerance for the linear independence check used in :class:`~.PauliVSpace`.
 
     Returns:
         Union[list[:class:`~.PauliSentence`], list[:class:`~.Operator`]]: a basis of either :class:`~.PauliSentence` or :class:`~.Operator` instances that is closed under
@@ -123,7 +125,7 @@ def lie_closure(
             for op in generators
         ]
 
-    vspace = PauliVSpace(generators)
+    vspace = PauliVSpace(generators, tol=tol)
 
     epoch = 0
     old_length = 0  # dummy value
@@ -141,7 +143,7 @@ def lie_closure(
             # remove common factor 2 with Pauli commutators
             for pw, val in com.items():
                 com[pw] = val.imag / 2
-            vspace.add(com)
+            vspace.add(com, tol=tol)
 
         # Updated number of linearly independent PauliSentences from previous and current step
         old_length = new_length
@@ -159,7 +161,7 @@ def lie_closure(
 
 
 class PauliVSpace:
-    """
+    r"""
     Class representing the linearly independent basis of a vector space.
 
     The main purpose of this class is to store and process ``M``, which
@@ -195,6 +197,8 @@ class PauliVSpace:
     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``.
+        tol (float): Numerical tolerance for the linear independence check. If the norm of the projection of the candidate vector
+            onto :math:`M^\perp` is greater than ``tol``, then it is deemed to be linearly independent.
 
     **Example**
 
@@ -234,7 +238,7 @@ class PauliVSpace:
      1.0 * X(0)]
     """
 
-    def __init__(self, generators, dtype=float):
+    def __init__(self, generators, dtype=float, tol=None):
 
         self.dtype = dtype
 
@@ -258,8 +262,10 @@ class PauliVSpace:
         self._rank = rank
         self._num_pw = num_pw
 
+        self.tol = np.finfo(self._M.dtype).eps * 100 if tol is None else tol
+
         # Add all generators that are linearly independent
-        self.add(generators)
+        self.add(generators, tol=tol)
 
     @property
     def basis(self):
@@ -269,7 +275,7 @@ class PauliVSpace:
     def __len__(self):
         return len(self.basis)
 
-    def add(self, other, tol=1e-15):
+    def add(self, other, tol=None):
         r"""Adding Pauli sentences if they are linearly independent.
 
         Args:
@@ -295,6 +301,9 @@ class PauliVSpace:
         [1.0 * X(0), 1.0 * X(1), 1.0 * Y(0), 1.0 * Z(0)]
 
         """
+        if tol is None:
+            tol = self.tol
+
         if isinstance(other, (qml.pauli.PauliWord, qml.pauli.PauliSentence, Operator)):
             other = [other]
 
@@ -318,7 +327,7 @@ class PauliVSpace:
                 self._basis.append(ps)
         return self._basis
 
-    def is_independent(self, pauli_sentence, tol=1e-15):
+    def is_independent(self, pauli_sentence, tol=None):
         r"""Check if the ``pauli_sentence`` is linearly independent of the basis of ``PauliVSpace``.
 
         Args:
@@ -338,13 +347,16 @@ class PauliVSpace:
         True
 
         """
+        if tol is None:
+            tol = self.tol
+
         _, _, _, _, 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):
+    def _check_independence(M, pauli_sentence, pw_to_idx, rank, num_pw, tol):
         r"""
         Checks if :class:`~PauliSentence` ``pauli_sentence`` is linearly independent and provides the updated class attributes in case the vector is added.
 
@@ -400,7 +412,7 @@ class PauliVSpace:
         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
+        v = v - A @ qml.math.linalg.solve(qml.math.conj(A.T) @ A, A.conj().T) @ v
 
         if np.linalg.norm(v) > tol:
             return M, pw_to_idx, rank + 1, new_num_pw, True


Test patch:
diff --git a/tests/pauli/dla/test_lie_closure.py b/tests/pauli/dla/test_lie_closure.py
index 04aa341b7..dc9dd4faa 100644
--- a/tests/pauli/dla/test_lie_closure.py
+++ b/tests/pauli/dla/test_lie_closure.py
@@ -19,7 +19,7 @@ import numpy as np
 
 import pennylane as qml
 
-from pennylane import X, Y, Z
+from pennylane import X, Y, Z, I
 
 from pennylane.pauli import PauliWord, PauliSentence, PauliVSpace, lie_closure
 
@@ -64,12 +64,17 @@ class TestPauliVSpace:
         assert vspace._rank == 2
         assert vspace._num_pw == 2
         assert len(vspace._pw_to_idx) == 2
+        assert vspace.tol == np.finfo(vspace._M.dtype).eps * 100
 
     @pytest.mark.parametrize("dtype", [float, complex])
     def test_dtype(self, dtype):
         vspace = PauliVSpace(ops1, dtype=dtype)
         assert vspace._M.dtype == dtype
 
+    def test_set_tol(self):
+        vspace = PauliVSpace(ops1, tol=1e-13)
+        assert vspace.tol == 1e-13
+
     def test_init_with_ops(self):
         """Test that initialization with PennyLane operators, PauliWord and PauliSentence works"""
         ops = [qml.X(0), PauliWord({1: "X"})]
@@ -91,7 +96,7 @@ class TestPauliVSpace:
             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("li_tol", [None, 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"""
@@ -272,7 +277,7 @@ class TestPauliVSpace:
         ),
     )
 
-    @pytest.mark.parametrize("tol", [1e-15, 1e-8])
+    @pytest.mark.parametrize("tol", [None, 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"""
@@ -487,3 +492,17 @@ class TestLieClosure:
 
         res = qml.pauli.lie_closure(generators)
         assert len(res) == 4 * ((2 ** (n - 2)) ** 2 - 1)
+
+    def test_universal_gate_set(self):
+        """Test universal gate set"""
+        n = 3
+
+        generators = [Z(i) for i in range(n)]
+        generators += [Y(i) for i in range(n)]
+        generators += [
+            (I(i) - Z(i)) @ (I(i + 1) - X(i + 1)) for i in range(n - 1)
+        ]  # generator of CNOT gate
+
+        vspace = qml.lie_closure(generators)
+
+        assert len(vspace) == 4**3