[BUGFIX] Fix qml.lie_closure by using SVD based linear independen…

…ce check (#6232)

fixes #6065

Additionally introduces normalization of vectors to avoid exploding
coefficients
Could alternatively also introduce orthonormalization with
`orthonormal_basis = U[:, :rank]` in the SVD step, but I dont think this
is necessary.

---------

Co-authored-by: Pietropaolo Frisoni <[email protected]>

doc/releases/changelog-dev.md

@@ -212,6 +212,9 @@
 * Fixes a bug where a simple circuit with no parameters or only builtin/numpy arrays as parameters returns autograd tensors.
   [(#6225)](https://github.com/PennyLaneAI/pennylane/pull/6225)
 
+* `qml.pauli.PauliVSpace` now uses a more stable SVD-based linear independence check to avoid running into `LinAlgError: Singular matrix`. This stabilizes the usage of `qml.lie_closure`. It also introduces normalization of the basis vector's internal representation `_M` to avoid exploding coefficients.
+  [(#6232)](https://github.com/PennyLaneAI/pennylane/pull/6232)
+
 * Fixes a bug where `csc_dot_product` is used during measurement for `Sum`/`Hamiltonian` that contains observables that does not define a sparse matrix.
   [(#6278)](https://github.com/PennyLaneAI/pennylane/pull/6278)
   [(#6310)](https://github.com/PennyLaneAI/pennylane/pull/6310)
@@ -229,6 +232,7 @@ Lillian M. A. Frederiksen,
 Pietropaolo Frisoni,
 Emiliano Godinez,
 Austin Huang,
+Korbinian Kottmann,
 Christina Lee,
 William Maxwell,
 Lee J. O'Riordan,

pennylane/pauli/dla/lie_closure.py

@@ -21,6 +21,7 @@
 from typing import Union
 
 import numpy as np
+import scipy
 
 import pennylane as qml
 from pennylane.operation import Operator
@@ -139,7 +140,7 @@ def lie_closure(
 
         for ps1, ps2 in itertools.combinations(vspace.basis, 2):
             com = ps1.commutator(ps2)
-            com.simplify()
+            com.simplify(tol=vspace.tol)
 
             if len(com) == 0:  # skip because operators commute
                 continue
@@ -415,20 +416,22 @@ def _check_independence(M, pauli_sentence, pw_to_idx, rank, num_pw, tol):
             for pw, value in pauli_sentence.items():
                 M[new_pw_to_idx[pw], rank] = value
 
+            M[:, rank] /= np.linalg.norm(M[:, rank])
+
             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
 
+        M[:, rank] /= np.linalg.norm(M[:, rank])
+
         # 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.solve(qml.math.conj(A.T) @ A, A.conj().T) @ v
+        s = scipy.linalg.svdvals(M)
+        new_rank = np.count_nonzero(s > tol)
 
-        if np.linalg.norm(v) > tol:
-            return M, pw_to_idx, rank + 1, new_num_pw, True
+        if rank + 1 == new_rank:
+            return M, pw_to_idx, new_rank, new_num_pw, True
 
         return M[:num_pw, :rank], pw_to_idx, rank, num_pw, False
 

tests/pauli/dla/test_lie_closure.py

@@ -56,8 +56,10 @@ def test_init(self):
         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]]
+        assert np.allclose(
+            vspace._M, [[1 / np.sqrt(2), 1.0], [1 / np.sqrt(2), 0.0]]
+        ) or np.allclose(
+            vspace._M, [[1 / np.sqrt(2), 0.0], [1 / np.sqrt(2), 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