fix normalization in structure_constants_dense and its tests

pennylane/labs/dla/structure_constants_dense.py

@@ -92,6 +92,7 @@ def structure_constants_dense(g: TensorLike, is_orthonormal: bool = True) -> Ten
         g = np.array(g)
 
     assert g.shape[2] == g.shape[1]
+    chi = g.shape[1]
     # Assert Hermiticity of the input. Otherwise we'll get the sign wrong
     assert np.allclose(g.conj().transpose((0, 2, 1)), g)
 
@@ -103,16 +104,19 @@ def structure_constants_dense(g: TensorLike, is_orthonormal: bool = True) -> Ten
 
     # project commutators on the basis of g, see docstring for details.
     # Axis ordering is (dimg, _chi_, *chi*) x (dimg, *chi*, dimg, _chi_) -> (dimg, dimg, dimg)
-    adj = (1j * np.tensordot(g, all_coms, axes=[[1, 2], [3, 1]])).real
+    # Normalize trace inner product by dimension chi
+    adj = (1j * np.tensordot(g / chi, all_coms, axes=[[1, 2], [3, 1]])).real
 
     if not is_orthonormal:
         # If the basis is not orthonormal, compute the Gram matrix and apply its
         # (pseudo-)inverse to the obtained projections. See the docstring for details.
         # The Gram matrix is just one additional diagonal contraction of the ``prod`` tensor,
         # across the Hilbert space dimensions. (dimg, _chi_, dimg, _chi_) -> (dimg, dimg)
+        # This contraction is missing the normalization factor 1/chi of the trace inner product.
         gram_inv = np.linalg.pinv(np.sum(np.diagonal(prod, axis1=1, axis2=3), axis=-1).real)
         # Axis ordering for contraction with gamma axis of raw structure constants:
         # (dimg, _dimg_), (_dimg_, dimg, dimg) -> (dimg, dimg, dim)
-        adj = np.tensordot(gram_inv, adj, axes=1)
+        # Here we add the missing normalization factor of the trace inner product (after inversion)
+        adj = np.tensordot(gram_inv * chi, adj, axes=1)
 
     return adj

pennylane/labs/tests/dla/test_structure_constants_dense.py

@@ -16,20 +16,24 @@
 
 import numpy as np
 import pytest
-import scipy as sp
 
 import pennylane as qml
 from pennylane import X, Y, Z
 from pennylane.labs.dla import (
     check_orthonormal,
     orthonormalize,
+    pauli_decompose,
     structure_constants_dense,
     trace_inner_product,
 )
 from pennylane.pauli import PauliSentence, PauliWord
 
 ## Construct some example DLAs
 # TFIM
+gens = [PauliSentence({PauliWord({i: "X", i + 1: "X"}): 1.0}) for i in range(1)]
+gens += [PauliSentence({PauliWord({i: "Z"}): 1.0}) for i in range(2)]
+Ising2 = qml.lie_closure(gens, pauli=True)
+
 gens = [PauliSentence({PauliWord({i: "X", i + 1: "X"}): 1.0}) for i in range(2)]
 gens += [PauliSentence({PauliWord({i: "Z"}): 1.0}) for i in range(3)]
 Ising3 = qml.lie_closure(gens, pauli=True)
@@ -68,39 +72,31 @@ def test_structure_constants_dim(self):
 
     def test_structure_constants_with_is_orthonormal(self):
         """Test that the structure constants with is_orthonormal=True/False match for
-        orthogonal inputs."""
+        orthonormal inputs."""
 
-        Ising3_dense = np.array([qml.matrix(op, wire_order=range(3)) for op in Ising3]) / np.sqrt(8)
+        Ising3_dense = np.array([qml.matrix(op, wire_order=range(3)) for op in Ising3])
+        assert check_orthonormal(Ising3_dense, trace_inner_product)
         adjoint_true = structure_constants_dense(Ising3_dense, is_orthonormal=True)
         adjoint_false = structure_constants_dense(Ising3_dense, is_orthonormal=False)
         assert np.allclose(adjoint_true, adjoint_false)
 
-    @pytest.mark.parametrize(
-        "dla, use_orthonormal",
-        [
-            (Ising3, True),
-            (Ising3, False),
-            (XXZ3, True),
-            (XXZ3, False),
-            (Heisenberg3_sum, True),
-            (Heisenberg3_sum, False),
-            (sum_XXZ3, True),
-            (sum_XXZ3, False),
-        ],
-    )
+    @pytest.mark.parametrize("dla", [Ising2, Ising3, XXZ3, Heisenberg3_sum, sum_XXZ3])
+    @pytest.mark.parametrize("use_orthonormal", [True, False])
     def test_structure_constants_elements(self, dla, use_orthonormal):
-        r"""Test relation :math:`[i G_\alpha, i G_\beta] = \sum_{\gamma = 0}^{\mathfrak{d}-1} f^\gamma_{\alpha, \beta} iG_\gamma`."""
+        r"""Test relation :math:`[i G_α, i G_β] = \sum_{γ=0}^{d-1} f^γ_{α,β} iG_γ_`."""
 
         d = len(dla)
         dla_dense = np.array([qml.matrix(op, wire_order=range(3)) for op in dla])
 
         if use_orthonormal:
-            dla = orthonormalize(dla)
-            assert check_orthonormal(dla, trace_inner_product)
             dla_dense = orthonormalize(dla_dense)
             assert check_orthonormal(dla_dense, trace_inner_product)
+            dla = pauli_decompose(dla_dense, pauli=True)
+            assert check_orthonormal(dla, trace_inner_product)
 
+        ad_rep_non_dense = qml.structure_constants(dla, is_orthogonal=False)
         ad_rep = structure_constants_dense(dla_dense, is_orthonormal=use_orthonormal)
+        assert np.allclose(ad_rep, ad_rep_non_dense)
         for i in range(d):
             for j in range(d):
 
@@ -116,15 +112,11 @@ def test_structure_constants_elements(self, dla, use_orthonormal):
     @pytest.mark.parametrize("use_orthonormal", [False, True])
     def test_use_operators(self, dla, use_orthonormal):
         """Test that operators can be passed and lead to the same result"""
-        ops = np.array([qml.matrix(op.operation(), wire_order=range(3)) for op in dla])
-
         if use_orthonormal:
-            gram_inv = sp.linalg.sqrtm(
-                np.linalg.pinv(np.tensordot(ops, ops, axes=[[1, 2], [2, 1]]).real)
-            )
-            ops = np.tensordot(gram_inv, ops, axes=1)
-            dla = [(scale * op).pauli_rep for scale, op in zip(np.diag(gram_inv), dla)]
+            dla = orthonormalize(dla)
+
+        ops = np.array([qml.matrix(op.operation(), wire_order=range(3)) for op in dla])
 
-        ad_rep_true = qml.pauli.dla.structure_constants(dla)
+        ad_rep_true = qml.structure_constants(dla)
         ad_rep = structure_constants_dense(ops, is_orthonormal=use_orthonormal)
         assert qml.math.allclose(ad_rep, ad_rep_true)