From d892b10d109b9cb4da48ecf8c272dc7265b8d926 Mon Sep 17 00:00:00 2001
From: Korbinian Kottmann <43949391+Qottmann@users.noreply.github.com>
Date: Wed, 2 Oct 2024 14:20:36 +0200
Subject: [PATCH] [BUGFIX] Fix ``qml.lie_closure`` by using SVD based linear
 independence check (#6232)

fixes https://github.com/PennyLaneAI/pennylane/issues/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 <pietropaolo.frisoni@xanadu.ai>
---
 doc/releases/changelog-dev.md       |  4 ++++
 pennylane/pauli/dla/lie_closure.py  | 17 ++++++++++-------
 tests/pauli/dla/test_lie_closure.py |  6 ++++--
 3 files changed, 18 insertions(+), 9 deletions(-)

diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 2c5dec1ec70..3a92a78bec1 100644
--- a/doc/releases/changelog-dev.md
+++ b/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,
diff --git a/pennylane/pauli/dla/lie_closure.py b/pennylane/pauli/dla/lie_closure.py
index 7ed9a939314..cf7176df586 100644
--- a/pennylane/pauli/dla/lie_closure.py
+++ b/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
 
diff --git a/tests/pauli/dla/test_lie_closure.py b/tests/pauli/dla/test_lie_closure.py
index 50babf8e000..413f6d3a376 100644
--- a/tests/pauli/dla/test_lie_closure.py
+++ b/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