Repository: PennyLaneAI/pennylane

Issue:
[BUG] `qml.center` does not work for all linear combinations of Paulis
### Expected behavior

`qml.center` always computes the center of an algebra.

### Actual behavior

If the algebra basis is such that the center elements are linear combinations of them, `qml.center` does not find them.

### Additional information

We need to overlap kernels of adjoint representations instead of commuting basis elements alone.

### Source code

```shell
>>> dla = qml.lie_closure([qml.X(0), qml.X(1), (qml.X(0)@qml.X(1) + qml.Y(0) @ qml.Y(1) + qml.Z(0)@ qml.Z(1))])
>>> print(dla)
[X(0),
 X(1),
 X(0) @ X(1) + Y(0) @ Y(1) + Z(0) @ Z(1),
 Z(0) @ Y(1) + -1.0 * (Y(0) @ Z(1)),
 -1.0 * (Y(0) @ Y(1)) + -1.0 * (Z(0) @ Z(1))]
>>> qml.center(dla)
[] # Should contain X(0) + X(1)
```


### Tracebacks

_No response_

### System information

```shell
pl dev
```


### Existing GitHub issues

- [X] I have searched existing GitHub issues to make sure the issue does not already exist.

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/6049?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 `95.00000%` with `1 line` in your changes missing coverage. Please review.
> Project coverage is 99.65%. Comparing base [(`9808d0c`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/9808d0ca19f104033970792bbbfe019e236679fc?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) to head [(`3e9e148`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/3e9e148bb213e1b3e07b545fd04518e8ffce1375?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).
> Report is 2 commits behind head on master.

| [Files](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/6049?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/center.py](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/6049?src=pr&el=tree&filepath=pennylane%2Fpauli%2Fdla%2Fcenter.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3BhdWxpL2RsYS9jZW50ZXIucHk=) | 95.00% | [1 Missing :warning: ](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/6049?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    #6049      +/-   ##
==========================================
- Coverage   99.66%   99.65%   -0.01%     
==========================================
  Files         430      430              
  Lines       41488    41198     -290     
==========================================
- Hits        41348    41056     -292     
- Misses        140      142       +2     
```

</details>

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

dwierichs: @Qottmann 
>Could you add a bit more description of the involved math, since this is not super trivial and not explained in the doc string, would be good to at least have some explanation in the code as comment.

Yep, good point! Will do!

>That being said, is it really necessary to compute three kernels here?

I'm not sure, maybe we can do Gram-Schmidt instead :thinking: I'll give it another attempt!
dwierichs: @Qottmann Found a (probably) easier way for the kernel intersections. Added derivation for that new technique to docstring and put the intersection code into a helper function with dev docstring.
Qottmann: Is this a bug or am I once again not fully getting the definition of the center 🫠 

```
gg = [
    Z(0) @ Z(1),
    Y(0) @ Y(1),
    Z(0) @ Y(1),
    Y(0) @ Z(1)
]
c = qml.center(gg) 
print(c)
# [Z(0) @ Z(1), Y(0) @ Y(1), Z(0) @ Y(1), Y(0) @ Z(1)]


gg = [g_i.pauli_rep for g_i in gg]
for g1 in gg:
    for g2 in gg:
        com = g1.commutator(g2)
        com.simplify()
        assert len(com) == 0 # fails with e.g. third output being -2j * X(1)
```

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 4eb3ec849..08688461f 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -277,6 +277,10 @@
 * Fix `jax.grad` + `jax.jit` not working for `AmplitudeEmbedding`, `StatePrep` and `MottonenStatePreparation`.
   [(#5620)](https://github.com/PennyLaneAI/pennylane/pull/5620) 
 
+* Fixed a bug in `qml.center` that omitted elements from the center if they were
+  linear combinations of input elements.
+  [(#6049)](https://github.com/PennyLaneAI/pennylane/pull/6049)
+
 * Fix a bug where the global phase returned by `one_qubit_decomposition` gained a broadcasting dimension.
   [(#5923)](https://github.com/PennyLaneAI/pennylane/pull/5923)
 
@@ -333,4 +337,4 @@ Anurav Modak,
 Mudit Pandey,
 Erik Schultheis,
 nate stemen,
-David Wierichs,
\ No newline at end of file
+David Wierichs,
diff --git a/pennylane/pauli/dla/center.py b/pennylane/pauli/dla/center.py
index 6fa5fc663..d810b429b 100644
--- a/pennylane/pauli/dla/center.py
+++ b/pennylane/pauli/dla/center.py
@@ -16,9 +16,50 @@ from itertools import combinations
 from typing import Union
 
 import numpy as np
+from scipy.linalg import norm, null_space
 
 from pennylane.operation import Operator
 from pennylane.pauli import PauliSentence, PauliWord
+from pennylane.pauli.dla import structure_constants
+
+
+def _intersect_bases(basis_0, basis_1):
+    r"""Compute the intersection of two vector spaces that are given by a basis each.
+    This is done by constructing a matrix [basis_0 | -basis_1] and computing its null space
+    in form of vectors (u, v)^T, which is equivalent to solving the equation
+    ``basis_0 @ u = basis_1 @ v``.
+    Given a basis for this null space, the vectors ``basis_0 @ u`` (or equivalently
+    ``basis_1 @ v``) form a basis for the intersection of the vector spaces.
+
+    Also see https://math.stackexchange.com/questions/25371/how-to-find-a-basis-for-the-intersection-of-two-vector-spaces-in-mathbbrn
+    """
+    # Compute (orthonormal) basis for the null space of the augmented matrix [basis_0, -basis_1]
+    augmented_basis = null_space(np.hstack([basis_0, -basis_1]))
+    # Compute basis_0 @ u for each vector u from the basis (u, v)^T in the augmented basis
+    intersection_basis = basis_0 @ augmented_basis[: basis_0.shape[1]]
+    # Normalize the output for cleaner results, because the augmented kernel was normalized
+    intersection_basis = intersection_basis / norm(intersection_basis, axis=0)
+    return intersection_basis
+
+
+def _center_pauli_words(g, pauli):
+    """Compute the center of an algebra given in a PauliWord basis."""
+    # Guarantees all operators to be given as a PauliWord
+    # If `pauli=True` we know that they are PauliWord or PauliSentence instances
+    g_pws = [o if isinstance(o, PauliWord) else next(iter(o.pauli_rep.keys())) for o in g]
+    d = len(g_pws)
+    commutators = np.zeros((d, d), dtype=int)
+    for (j, op1), (k, op2) in combinations(enumerate(g_pws), r=2):
+        if not op1.commutes_with(op2):
+            commutators[j, k] = 1  # dummy value to indicate operators dont commute
+            commutators[k, j] = 1
+
+    ids = np.where(np.all(commutators == 0, axis=0))[0]
+    res = [g[idx] for idx in ids]
+
+    if not pauli:
+        res = [op.operation() if isinstance(op, (PauliWord, PauliSentence)) else op for op in res]
+    return res
 
 
 def center(
@@ -33,9 +74,11 @@ def center(
     .. math:: \mathfrak{\xi}(\mathfrak{g}) := \{h \in \mathfrak{g} | [h, h_i]=0 \ \forall h_i \in \mathfrak{g} \}
 
     Args:
-        g (List[Union[Operator, PauliSentence, PauliWord]]): List of operators for which to find the center.
-        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`
+        g (List[Union[Operator, PauliSentence, PauliWord]]): List of operators that spans
+            the algebra for which to find the center.
+        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``.
 
     Returns:
@@ -56,23 +99,94 @@ def center(
     >>> qml.center(g)
     [X(0)]
 
+    .. details::
+        :title: Derivation
+        :href: derivation
+
+        The center :math:`\mathfrak{\xi}(\mathfrak{g})` of an algebra :math:`\mathfrak{g}`
+        can be computed in the following steps. First, compute the
+        :func:`~.pennylane.structure_constants`, or adjoint representation, of the algebra
+        with respect to some basis :math:`\mathbb{B}` of :math:`\mathfrak{g}`.
+        The center of :math:`\mathfrak{g}` is then given by
+
+        .. math::
+
+            \mathfrak{\xi}(\mathfrak{g}) = \operatorname{span}\left\{\bigcap_{x\in\mathbb{B}}
+            \operatorname{ker}(\operatorname{ad}_x)\right\},
+
+        i.e., the intersection of the kernels, or null spaces, of all basis elements in the
+        adjoint representation.
+
+        The kernel can be computed with ``scipy.linalg.null_space``, and vector space
+        intersections are computed recursively from pairwise intersections. The intersection
+        between two vectors spaces :math:`V_1` and :math:`V_2` given by (orthonormal) bases
+        :math:`\mathbb{B}_i` can be computed from the kernel of the matrix that has all
+        basis vectors from :math:`\mathbb{B}_1` and :math:`-\mathbb{B}_2` as columns, i.e.,
+        :math:`\operatorname{ker}(\left[\ \mathbb{B}_1 \ | -\mathbb{B}_2\ \right])`. For an
+        (orthonormal) basis of this kernel, consisting of two stacked column vectors
+        :math:`u^{(i)}_1` and :math:`u^{(i)}_2` for each basis, a basis of the
+        intersection space :math:`V_1 \cap V_2` is given by :math:`\{\mathbb{B}_1 u_1^{(i)}\}_i`
+        (or equivalently by :math:`\{\mathbb{B}_2 u_2^{(i)}\}_i`).
+        Also see `this post <https://math.stackexchange.com/questions/25371/how-to-find-a-basis-for-the-intersection-of-two-vector-spaces-in-mathbbrn>`_
+        for details.
+
+        If the input consists of :class:`~.pennylane.PauliWord` instances only, we can
+        instead compute pairwise commutators and know that the center consists solely of
+        basis elements that commute with all other basis elements. This can be seen in the
+        following way.
+
+        Assume that the center elements identified based on the basis have been removed
+        already and we are left with a basis :math:`\mathbb{B}=\{p_i\}_i` of Pauli
+        words such that :math:`\forall i\ \exists j:\ [p_i, p_j] \neq 0`. Assume that there is
+        another center element :math:`x\neq 0`, which was missed before because it is a linear
+        combination of Pauli words:
+
+        .. math::
+
+            \forall j: \ [x, p_j] = [\sum_i x_i p_i, p_j] = 0.
+
+        As products of Paulis are unique when fixing one of the factors (:math:`p_j` is fixed
+        above), we then know that
+
+        .. math::
+
+            &\forall j: \ 0 = \sum_i x_i [p_i, p_j] = 2 \sum_i x_i \chi_{i,j} p_ip_j\\
+            \Rightarrow &\forall i,j \text{ s.t. } \chi_{i,j}\neq 0: x_i = 0,
+
+        where :math:`\chi_{i,j}` denotes an indicator that is :math:`0` if the commutator
+        :math:`[p_i, p_j]` vanishes and :math:`1` otherwise.
+        However, we know that for each :math:`i` there is at least one :math:`j` such that
+        :math:`\chi_{i,j}\neq 0`. This means that :math:`x_i = 0` is guaranteed for all
+        :math:`i` by at least one :math:`j`. Therefore :math:`x=0`, which is a contradiction
+        to our initial assumption that :math:`x\neq 0`.
     """
+    if len(g) < 2:
+        # A length-zero list has zero center, a length-one list has full center
+        return g
+    if all(isinstance(x, PauliWord) or len(x.pauli_rep) == 1 for x in g):
+        return _center_pauli_words(g, pauli)
+
+    adjoint_repr = structure_constants(g, pauli)
+    # Start kernels intersection with kernel of first DLA element
+    kernel_intersection = null_space(adjoint_repr[0])
+    for ad_x in adjoint_repr[1:]:
+        # Compute the next kernel and intersect it with previous intersection
+        next_kernel = null_space(ad_x)
+        kernel_intersection = _intersect_bases(kernel_intersection, next_kernel)
+
+        # If the intersection is zero-dimensional, exit early
+        if kernel_intersection.shape[1] == 0:
+            return []
+
+    # Construct operators from numerical output and convert to desired format
+    res = [sum(c * x for c, x in zip(c_coeffs, g)) for c_coeffs in kernel_intersection.T]
+
+    have_paulis = all(isinstance(x, (PauliWord, PauliSentence)) for x in res)
+    if pauli or have_paulis:
+        _ = [el.simplify() for el in res]
+        if not pauli:
+            res = [el.operation() for el in res]
+    else:
+        res = [el.simplify() for el in res]
 
-    if not pauli:
-        g = [o.pauli_rep for o in g]
-
-    d = len(g)
-    commutators = np.zeros((d, d), dtype=int)
-    for (j, op1), (k, op2) in combinations(enumerate(g), r=2):
-        res = op1.commutator(op2)
-        res.simplify()
-        if res != PauliSentence({}):
-            commutators[j, k] = 1  # dummy value to indicate operators dont commute
-            commutators[k, j] = 1
-
-    mask = np.all(commutators == 0, axis=0)
-    res = list(np.array(g)[mask])
-
-    if not pauli:
-        res = [op.operation() for op in res]
     return res


Test patch:
diff --git a/tests/pauli/dla/test_center.py b/tests/pauli/dla/test_center.py
index 07d724b80..efa0045c5 100644
--- a/tests/pauli/dla/test_center.py
+++ b/tests/pauli/dla/test_center.py
@@ -13,17 +13,16 @@
 # limitations under the License.
 """Tests for pennylane/dla/center.py functionality"""
 
+import numpy as np
 import pytest
 
 import pennylane as qml
-from pennylane.pauli import PauliSentence, center
+from pennylane.pauli import PauliSentence, PauliWord, center
 
 
 def test_trivial_center():
     """Test that the center of an empty list of generators is an empty list of generators."""
-    ops = []
-    res = center(ops)
-    assert res == []
+    assert center([]) == []
 
 
 DLA_CENTERS = (
@@ -32,14 +31,17 @@ DLA_CENTERS = (
     ([qml.X(0), qml.X(1)], [qml.X(0), qml.X(1)]),  # two non-overlapping wires
     ([qml.X(0), qml.Y(1)], [qml.X(0), qml.Y(1)]),  # two non-overlapping wires with different ops
     ([qml.X(0), qml.Y(0), qml.Z(0), qml.I()], [qml.I()]),  # non-trivial DLA, but trivial center
+    (
+        [qml.X(0) + qml.Y(0), qml.X(0) - qml.Y(0), qml.Z(0)],
+        [],
+    ),  # non-trivial DLA, but trivial center
 )
 
 
 @pytest.mark.parametrize("ops, true_res", DLA_CENTERS)
 def test_center(ops, true_res):
     """Test centers with Identity operators or non-overlapping wires"""
-    res = center(ops)
-    assert res == true_res
+    assert center(ops) == true_res
 
 
 @pytest.mark.parametrize("ops, true_res", DLA_CENTERS)
@@ -49,30 +51,53 @@ def test_center_pauli(ops, true_res):
     res = center(ops, pauli=True)
 
     assert all(isinstance(op, PauliSentence) for op in res)
-    true_res = [op.pauli_rep for op in true_res]
-    assert res == true_res
+    assert res == [op.pauli_rep for op in true_res]
 
 
 @pytest.mark.parametrize("pauli", [False, True])
-def test_center_pauli_word_pauli_True(pauli):
+def test_center_pauli_word(pauli):
     """Test that PauliWord instances can be passed for both pauli=True/False"""
-    ops = [
-        qml.pauli.PauliWord({0: "X"}),
-        qml.pauli.PauliWord({0: "X", 1: "X"}),
-        qml.pauli.PauliWord({1: "Y"}),
-    ]
+    words = [{0: "X"}, {0: "X", 1: "X"}, {1: "Y"}, {0: "X", 1: "Z"}]
+    ops = list(map(PauliWord, words))
     if pauli:
-        assert qml.center(ops, pauli=pauli) == [qml.pauli.PauliWord({0: "X"})]
+        assert qml.center(ops, pauli=pauli) == [PauliWord({0: "X"})]
     else:
         assert qml.center(ops, pauli=pauli) == [qml.X(0)]
 
 
+@pytest.mark.parametrize("pauli", [False, True])
+def test_center_pauli_sentence(pauli):
+    """Test that PauliSentence instances can be passed for both pauli=True/False"""
+    words = [{0: "X"}, {0: "X", 1: "X"}, {1: "Y"}, {0: "X", 1: "Z"}]
+    words = list(map(PauliWord, words))
+    sentences = [
+        {words[0]: 0.5, words[1]: 3.2},
+        {words[0]: -0.2, words[2]: 2.5},
+        {words[2]: 1.2, words[3]: 0.72, words[1]: 0.6},
+        {words[1]: 0.9, words[2]: 1.8},
+    ]
+    sentences = list(map(PauliSentence, sentences))
+    if pauli:
+        cent = qml.center(sentences, pauli=pauli)
+        assert isinstance(cent, list) and len(cent) == 1
+        assert isinstance(cent[0], PauliSentence)
+        assert PauliWord({0: "X"}) in cent[0]
+    else:
+        cent = qml.center(sentences, pauli=pauli)
+        assert isinstance(cent, list) and len(cent) == 1
+        assert isinstance(cent[0], qml.ops.op_math.SProd)
+        assert cent[0].base == qml.X(0)
+
+
+c = 1 / np.sqrt(2)
+
 GENERATOR_CENTERS = (
     ([qml.X(0), qml.X(0) @ qml.X(1), qml.Y(1)], [qml.X(0)]),
     ([qml.X(0) @ qml.X(1), qml.Y(1), qml.X(0)], [qml.X(0)]),
     ([qml.X(0) @ qml.X(1), qml.Y(1), qml.X(1)], []),
     ([qml.X(0) @ qml.X(1), qml.Y(1), qml.Z(0)], []),
     ([p(0) @ p(1) for p in [qml.X, qml.Y, qml.Z]], [p(0) @ p(1) for p in [qml.X, qml.Y, qml.Z]]),
+    ([qml.X(0), qml.X(1), sum(p(0) @ p(1) for p in [qml.Y, qml.Z])], [c * qml.X(0) + c * qml.X(1)]),
 )
 
 
@@ -80,5 +105,4 @@ GENERATOR_CENTERS = (
 def test_center_dla(generators, true_res):
     """Test computing the center for a non-trivial DLA"""
     g = qml.pauli.lie_closure(generators)
-    res = center(g)
-    assert res == true_res
+    assert center(g) == true_res