From 7bdf71b44318176e937206e61ccdb47a78e38ce7 Mon Sep 17 00:00:00 2001 From: Korbinian Kottmann <43949391+Qottmann@users.noreply.github.com> Date: Wed, 5 Mar 2025 11:11:00 +0100 Subject: [PATCH] Upgrade `qml.structure_constants` to handle dense matrices (#6861) Integrate structure constants dense into [sc-81966] Todo - [x] update and unify documentation --------- Co-authored-by: David Wierichs Co-authored-by: Yushao Chen (Jerry) Co-authored-by: ANTH0NY <39093564+AntonNI8@users.noreply.github.com> Co-authored-by: Diego <67476785+DSGuala@users.noreply.github.com> --- doc/releases/changelog-dev.md | 6 + pennylane/labs/dla/__init__.py | 2 - pennylane/labs/dla/cartan_subalgebra.py | 3 +- .../labs/dla/structure_constants_dense.py | 123 ------------ .../dla/test_structure_constants_dense.py | 121 ------------ .../labs/tests/dla/test_variational_kak.py | 7 +- pennylane/pauli/dla/structure_constants.py | 177 +++++++++++++++++- tests/pauli/dla/test_structure_constants.py | 72 ++++++- 8 files changed, 249 insertions(+), 262 deletions(-) delete mode 100644 pennylane/labs/dla/structure_constants_dense.py delete mode 100644 pennylane/labs/tests/dla/test_structure_constants_dense.py diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md index 0c83b8e413e..7a875cf8cd3 100644 --- a/doc/releases/changelog-dev.md +++ b/doc/releases/changelog-dev.md @@ -41,6 +41,9 @@ Also added ``qml.pauli.trace_inner_product`` that can handle batches of dense matrices. [(#6811)](https://github.com/PennyLaneAI/pennylane/pull/6811) +* ``qml.structure_constants`` now accepts and outputs matrix inputs using the ``matrix`` keyword. + [(#6861)](https://github.com/PennyLaneAI/pennylane/pull/6861) +

Improvements 🛠

* Added a class `qml.capture.transforms.MergeAmplitudeEmbedding` that merges `qml.AmplitudeEmbedding` operators @@ -326,6 +329,9 @@ * ``pennylane.labs.dla.lie_closure_dense`` is removed and integrated into ``qml.lie_closure`` using the new ``dense`` keyword. [(#6811)](https://github.com/PennyLaneAI/pennylane/pull/6811) +* ``pennylane.labs.dla.structure_constants_dense`` is removed and integrated into ``qml.structure_constants`` using the new ``matrix`` keyword. + [(#6861)](https://github.com/PennyLaneAI/pennylane/pull/6861) + * ``ResourceOperator.resource_params`` is changed to a property. [(#6973)](https://github.com/PennyLaneAI/pennylane/pull/6973) diff --git a/pennylane/labs/dla/__init__.py b/pennylane/labs/dla/__init__.py index 20fc68798b8..7b27024be7e 100644 --- a/pennylane/labs/dla/__init__.py +++ b/pennylane/labs/dla/__init__.py @@ -20,7 +20,6 @@ .. autosummary:: :toctree: api - ~structure_constants_dense ~cartan_decomp ~recursive_cartan_decomp ~cartan_subalgebra @@ -80,7 +79,6 @@ """ -from .structure_constants_dense import structure_constants_dense from .cartan import ( cartan_decomp, recursive_cartan_decomp, diff --git a/pennylane/labs/dla/cartan_subalgebra.py b/pennylane/labs/dla/cartan_subalgebra.py index e43ebe3205d..36da4a6cf22 100644 --- a/pennylane/labs/dla/cartan_subalgebra.py +++ b/pennylane/labs/dla/cartan_subalgebra.py @@ -173,9 +173,8 @@ def cartan_subalgebra( all remaining operators from ``m``. >>> import numpy as np - >>> from pennylane.labs.dla import structure_constants_dense >>> g = np.vstack([k, m]) # re-order g to separate k and m operators - >>> adj = structure_constants_dense(g) # compute adjoint representation of g + >>> adj = qml.structure_constants(g, matrix=True) # compute adjoint representation of g Finally, we can compute a Cartan subalgebra :math:`\mathfrak{a}`, a maximal Abelian subalgebra of :math:`\mathfrak{m}`. diff --git a/pennylane/labs/dla/structure_constants_dense.py b/pennylane/labs/dla/structure_constants_dense.py deleted file mode 100644 index a76cfa01ed8..00000000000 --- a/pennylane/labs/dla/structure_constants_dense.py +++ /dev/null @@ -1,123 +0,0 @@ -# Copyright 2024 Xanadu Quantum Technologies Inc. - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -"""A function to compute the structure constants that make up -the adjoint representation of a Lie algebra in a dense matrix representation.""" - -import numpy as np - -from pennylane.typing import TensorLike - - -def structure_constants_dense(g: TensorLike, is_orthonormal: bool = True) -> TensorLike: - r""" - Compute the structure constants that make up the adjoint representation of a Lie algebra. - - This function computes the structure constants of a Lie algebra provided by their dense matrix representation, - obtained from, e.g., :func:`~lie_closure` with `matrix=True`. - This is sometimes more efficient than using the sparse Pauli representations of :class:`~PauliWord` and - :class:`~PauliSentence` that are employed in :func:`~structure_constants`, e.g., when there are few generators - that are sums of many Paulis. - - .. seealso:: For details on the mathematical definitions, see :func:`~structure_constants` and the section "Lie algebra basics" in our :doc:`g-sim demo `. - - Args: - g (np.array): The (dynamical) Lie algebra provided as dense matrices, as generated from :func:`~lie_closure` with `matrix=True`. - ``g`` should have shape ``(d, 2**n, 2**n)`` where ``d`` is the dimension of the algebra and ``n`` is the number of qubits. Each matrix ``g[i]`` should be Hermitian. - is_orthonormal (bool): Whether or not the matrices in ``g`` are orthonormal with respect to the Hilbert-Schmidt inner product on - (skew-)Hermitian matrices. If the inputs are orthonormal, it is recommended to set ``is_orthonormal`` to ``True`` to reduce - computational cost. Defaults to ``True``. - - Returns: - TensorLike: The adjoint representation of shape ``(d, d, d)``, corresponding to - the indices ``(gamma, alpha, beta)``. - - **Example** - - Let us generate the DLA of the transverse field Ising model using :func:`~lie_closure`. - - >>> n = 4 - >>> gens = [qml.X(i) @ qml.X(i+1) + qml.Y(i) @ qml.Y(i+1) + qml.Z(i) @ qml.Z(i+1) for i in range(n-1)] - >>> g = qml.lie_closure(gens, matrix=True) - >>> g.shape - (12, 16, 16) - - The DLA is represented by a collection of twelve :math:`2^4 \times 2^4` matrices. - Hence, the dimension of the DLA is :math:`d = 12` and the structure constants have shape ``(12, 12, 12)``. - - >>> from pennylane.labs.dla import structure_constants_dense - >>> adj = structure_constants_dense(g) - >>> adj.shape - (12, 12, 12) - - **Internal representation** - - As mentioned above, the input is assumed to be a batch of Hermitian matrices, even though - algebra elements are usually skew-Hermitian. That is, the input should represent the operators - :math:`G_\alpha` for an algebra basis :math:`\{iG_\alpha\}_\alpha`. - In an orthonormal basis of this form, the structure constants can then be computed simply via - - .. math:: - - f^\gamma_{\alpha, \beta} = \text{tr}[-i G_\gamma[iG_\alpha, iG_\beta]] = i\text{tr}[G_\gamma [G_\alpha, G_\beta]] \in \mathbb{R}. - - **Structure constants in non-orthonormal bases** - - Structure constants are often discussed using an orthonormal basis of the algebra. - This function can deal with non-orthonormal bases as well. For this, the Gram - matrix :math:`g` between the basis elements is taken into account when computing the overlap - of a commutator :math:`[iG_\alpha, iG_\beta]` with all algebra elements :math:`iG_\gamma`. - The resulting formula reads - - .. math:: - - f^\gamma_{\alpha, \beta} &= \sum_\eta g^{-1}_{\gamma\eta} i \text{tr}[G_\eta [G_\alpha, G_\beta]]\\ - g_{\gamma \eta} &= \text{tr}[G_\gamma G_\eta] \quad(\in\mathbb{R}) - - Internally, the commutators are computed by evaluating all operator products and subtracting - suitable pairs of products from each other. These products can be reused to evaluate the - Gram matrix as well. - """ - - if isinstance(g, list): - 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) - - # compute all commutators by computing all products first. - # Axis ordering is (dimg, chi, _chi_) x (dimg, _chi_, chi) -> (dimg, chi, dimg, chi) - prod = np.tensordot(g, g, axes=[[2], [1]]) - # The commutators now are the difference of prod with itself, with dimg axes swapped - all_coms = prod - np.transpose(prod, (2, 1, 0, 3)) - - # 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) - # 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) - # 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 diff --git a/pennylane/labs/tests/dla/test_structure_constants_dense.py b/pennylane/labs/tests/dla/test_structure_constants_dense.py deleted file mode 100644 index ece54de9089..00000000000 --- a/pennylane/labs/tests/dla/test_structure_constants_dense.py +++ /dev/null @@ -1,121 +0,0 @@ -# Copyright 2024 Xanadu Quantum Technologies Inc. - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -"""Tests for pennylane/labs/dla/structure_constants.py functionality""" -# pylint: disable=no-self-use - -import numpy as np -import pytest - -import pennylane as qml -from pennylane import X, Y, Z -from pennylane.labs.dla import ( - batched_pauli_decompose, - check_orthonormal, - orthonormalize, - structure_constants_dense, -) -from pennylane.pauli import PauliSentence, PauliWord, trace_inner_product - -## 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) - -# XXZ-type DLA, i.e. with true PauliSentences -gens2 = [ - PauliSentence( - { - PauliWord({i: "X", i + 1: "X"}): 1.0, - PauliWord({i: "Y", i + 1: "Y"}): 1.0, - } - ) - for i in range(2) -] -gens2 += [PauliSentence({PauliWord({i: "Z"}): 1.0}) for i in range(3)] -XXZ3 = qml.lie_closure(gens2, pauli=True) - -gens3 = [X(i) @ X(i + 1) + Y(i) @ Y(i + 1) + Z(i) @ Z(i + 1) for i in range(2)] -Heisenberg3_sum = qml.lie_closure(gens3) -Heisenberg3_sum = [op.pauli_rep for op in Heisenberg3_sum] - -coeffs = np.random.random((len(XXZ3), len(XXZ3))) -sum_XXZ3 = [qml.sum(*(c * op for c, op in zip(_coeffs, XXZ3))).pauli_rep for _coeffs in coeffs] - - -class TestAdjointRepr: - """Tests for structure_constants""" - - def test_structure_constants_dim(self): - """Test the dimension of the adjoint repr""" - d = len(Ising3) - Ising3_dense = np.array([qml.matrix(op, wire_order=range(3)) for op in Ising3]) - adjoint = structure_constants_dense(Ising3_dense) - assert adjoint.shape == (d, d, d) - assert adjoint.dtype == float - - def test_structure_constants_with_is_orthonormal(self): - """Test that the structure constants with is_orthonormal=True/False match for - orthonormal inputs.""" - - 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", [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_α, 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_dense = orthonormalize(dla_dense) - assert check_orthonormal(dla_dense, trace_inner_product) - dla = batched_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): - - comm_res = 1j * dla[i].commutator(dla[j]) - comm_res.simplify() - res = sum((c + 0j) * dla[gamma] for gamma, c in enumerate(ad_rep[:, i, j])) - res.simplify() - assert set(comm_res) == set(res) # Assert equal keys - if len(comm_res) > 0: - assert np.allclose(*zip(*[(comm_res[key], res[key]) for key in res.keys()])) - - @pytest.mark.parametrize("dla", [Ising3, XXZ3]) - @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""" - if use_orthonormal: - dla = orthonormalize(dla) - - ops = np.array([qml.matrix(op.operation(), wire_order=range(3)) for op in 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) diff --git a/pennylane/labs/tests/dla/test_variational_kak.py b/pennylane/labs/tests/dla/test_variational_kak.py index ff83edf6d3a..081bd50ae75 100644 --- a/pennylane/labs/tests/dla/test_variational_kak.py +++ b/pennylane/labs/tests/dla/test_variational_kak.py @@ -25,7 +25,6 @@ check_cartan_decomp, concurrence_involution, orthonormalize, - structure_constants_dense, validate_kak, variational_kak_adj, ) @@ -54,7 +53,7 @@ def test_kak_Ising(n, dense): adj = qml.structure_constants(g) else: g = np.vstack([k, m]) - adj = structure_constants_dense(g) + adj = qml.structure_constants(g, matrix=True) g, k, mtilde, h, adj = cartan_subalgebra(g, k, m, adj, tol=1e-10, start_idx=0) @@ -92,7 +91,7 @@ def test_kak_Heisenberg(n, dense): adj = qml.structure_constants(g) else: g = np.vstack([k, m]) - adj = structure_constants_dense(g) + adj = qml.structure_constants(g, matrix=True) g, k, mtilde, h, adj = cartan_subalgebra(g, k, m, adj, tol=1e-10, start_idx=0) @@ -132,7 +131,7 @@ def test_kak_Heisenberg_summed(is_orthogonal, dense): adj = qml.structure_constants(g, is_orthogonal=is_orthogonal) else: g = np.vstack([k, m]) - adj = structure_constants_dense(g, is_orthonormal=is_orthogonal) + adj = qml.structure_constants(g, matrix=True, is_orthogonal=is_orthogonal) g, k, mtilde, h, adj = cartan_subalgebra( g, k, m, adj, tol=1e-10, start_idx=0, is_orthogonal=is_orthogonal diff --git a/pennylane/pauli/dla/structure_constants.py b/pennylane/pauli/dla/structure_constants.py index a8f1d101f16..8e10862097c 100644 --- a/pennylane/pauli/dla/structure_constants.py +++ b/pennylane/pauli/dla/structure_constants.py @@ -17,6 +17,7 @@ import numpy as np +import pennylane as qml from pennylane.operation import Operator from pennylane.typing import TensorLike @@ -34,8 +35,9 @@ def _all_commutators(ops): def structure_constants( - g: list[Union[Operator, PauliWord, PauliSentence]], + g: list[Union[Operator, PauliWord, PauliSentence, TensorLike]], pauli: bool = False, + matrix: bool = False, is_orthogonal: bool = True, ) -> TensorLike: r""" @@ -60,6 +62,7 @@ def structure_constants( pauli (bool): Indicates whether it is assumed that :class:`~.PauliSentence` or :class:`~.PauliWord` instances are input. This can help with performance to avoid unnecessary conversions to :class:`~pennylane.operation.Operator` and vice versa. Default is ``False``. + matrix (bool): Whether or not matrix representations are used in the structure constants computation. This can help speed up the computation when using sums of Paulis with many terms. Default is ``False``. is_orthogonal (bool): Whether the set of operators in ``g`` is orthogonal with respect to the trace inner product. Default is ``True``. @@ -116,6 +119,28 @@ def structure_constants( >>> adjoint_rep[:, 0, 1] # commutator of X_0 and Y_0 consists of first and last operator array([-2., 0., 2.]) + We can also use matrix representations for the computation, which is sometimes faster, in particular for sums of many Pauli words. + This only affects how the structure constants are computed internally, it does not change the result. + + >>> adjoint_rep2 = qml.structure_constants(dla, is_orthogonal=False, matrix=True) + >>> qml.math.allclose(adjoint_rep, adjoint_rep2) + True + + We can also input the DLA as a list of matrices. For that we use :func:`~lie_closure` with ``matrix=True``. + + >>> n = 4 + >>> gens = [qml.X(i) @ qml.X(i+1) + qml.Y(i) @ qml.Y(i+1) + qml.Z(i) @ qml.Z(i+1) for i in range(n-1)] + >>> g = qml.lie_closure(gens, matrix=True) + >>> g.shape + (12, 16, 16) + + The DLA is represented by a collection of twelve :math:`2^4 \times 2^4` matrices. + Hence, the dimension of the DLA is :math:`d = 12` and the structure constants have shape ``(12, 12, 12)``. + + >>> adj = qml.structure_constants(g, matrix=True) + >>> adj.shape + (12, 12, 12) + .. details:: :title: Mathematical details @@ -162,9 +187,12 @@ def structure_constants( be skipped. """ - if any((op.pauli_rep is None) for op in g): + if matrix: + return _structure_constants_matrix(g, is_orthogonal) + + if any((getattr(op, "pauli_rep", None) is None) for op in g): raise ValueError( - f"Cannot compute adjoint representation of non-pauli operators. Received {g}." + f"Cannot compute adjoint representation of non-pauli operators. Received {g}. If you want to use matrices, use structure_constants(.., matrix=True)" ) if not pauli: @@ -191,3 +219,146 @@ def structure_constants( rep = np.tensordot(np.linalg.pinv(gram), rep, axes=[[-1], [0]]) return rep + + +def _structure_constants_matrix(g: TensorLike, is_orthogonal: bool = True) -> TensorLike: + r""" + Compute the structure constants that make up the adjoint representation of a Lie algebra. + + This function computes the structure constants of a Lie algebra provided by their matrix representation, + obtained from, e.g., :func:`~lie_closure`. + This is sometimes more efficient than using the sparse Pauli representations of :class:`~PauliWord` and + :class:`~PauliSentence` that are employed in :func:`~structure_constants`, e.g., when there are few generators + that are sums of many Paulis. + + .. seealso:: For details on the mathematical definitions, see :func:`~structure_constants` and the section "Lie algebra basics" in our `g-sim demo `__. + + Args: + g (np.array): The (dynamical) Lie algebra provided as matrix matrices, as generated from :func:`~lie_closure`. + ``g`` should have shape ``(d, 2**n, 2**n)`` where ``d`` is the dimension of the algebra and ``n`` is the number of qubits. Each matrix ``g[i]`` should be Hermitian. + is_orthogonal (bool): Whether or not the matrices in ``g`` are orthogonal with respect to the Hilbert-Schmidt inner product on + (skew-)Hermitian matrices. If the inputs are orthogonal, it is recommended to set ``is_orthogonal`` to ``True`` to reduce + computational cost. Defaults to ``True``. + + Returns: + TensorLike: The adjoint representation of shape ``(d, d, d)``, corresponding to + the indices ``(gamma, alpha, beta)``. + + **Example** + + Let us generate the DLA of the transverse field Ising model using :func:`~lie_closure`. + + >>> n = 4 + >>> gens = [qml.X(i) @ qml.X(i+1) + qml.Y(i) @ qml.Y(i+1) + qml.Z(i) @ qml.Z(i+1) for i in range(n-1)] + >>> g = qml.lie_closure(gens, matrix=True) + >>> g.shape + (12, 16, 16) + + The DLA is represented by a collection of twelve :math:`2^4 \times 2^4` matrices. + Hence, the dimension of the DLA is :math:`d = 12` and the structure constants have shape ``(12, 12, 12)``. + + >>> adj = qml.structure_constants(g, matrix=True) + >>> adj.shape + (12, 12, 12) + + **Internal representation** + + As mentioned above, the input is assumed to be a batch of Hermitian matrices, even though + algebra elements are usually skew-Hermitian. That is, the input should represent the operators + :math:`G_\alpha` for an algebra basis :math:`\{iG_\alpha\}_\alpha`. + In an orthonormal basis of this form, the structure constants can then be computed simply via + + .. math:: + + f^\gamma_{\alpha, \beta} = \text{tr}[-i G_\gamma[iG_\alpha, iG_\beta]] = i\text{tr}[G_\gamma [G_\alpha, G_\beta]] \in \mathbb{R}. + + Possible deviations of an orthogonal basis from normalization is taken into account in a + reduced version of the step for non-orthogonal bases below. + + **Structure constants in non-orthogonal bases** + + Structure constants are often discussed using an orthogonal basis of the algebra. + This function can deal with non-orthogonal bases as well. For this, the Gram + matrix :math:`g` between the basis elements is taken into account when computing the overlap + of a commutator :math:`[iG_\alpha, iG_\beta]` with all algebra elements :math:`iG_\gamma`. + The resulting formula reads + + .. math:: + + f^\gamma_{\alpha, \beta} &= \sum_\eta g^{-1}_{\gamma\eta} i \text{tr}[G_\eta [G_\alpha, G_\beta]]\\ + g_{\gamma \eta} &= \text{tr}[G_\gamma G_\eta] \quad(\in\mathbb{R}) + + Internally, the commutators are computed by evaluating all operator products and subtracting + suitable pairs of products from each other. These products can be reused to evaluate the + Gram matrix as well. + For orthogonal but not normalized bases, a reduced version of this step is used, only + computing (and inverting) the diagonal of the Gram matrix. + """ + + if getattr(g[0], "wires", False): + # operator input + all_wires = qml.wires.Wires.all_wires([_.wires for _ in g]) + n = len(all_wires) + assert all_wires.toset() == set(range(n)) + + g = qml.math.array( + [qml.matrix(op, wire_order=range(n)) for op in g], dtype=complex, like=g[0] + ) + chi = 2**n + assert np.shape(g) == (len(g), chi, chi) + + interface = qml.math.get_interface(g[0]) + + if isinstance(g[0], TensorLike) and isinstance(g, (list, tuple)): + # list of matrices + g = qml.math.stack(g, like=interface) + + chi = qml.math.shape(g[0])[0] + assert qml.math.shape(g) == (len(g), chi, chi) + assert qml.math.allclose( + qml.math.transpose(qml.math.conj(g), (0, 2, 1)), g + ), "Input matrices to structure_constants not Hermitian" + + # compute all commutators by computing all products first. + # Axis ordering is (dimg, chi, _chi_) x (dimg, _chi_, chi) -> (dimg, chi, dimg, chi) + prod = qml.math.tensordot(g, g, axes=[[2], [1]]) + # The commutators now are the difference of prod with itself, with dimg axes swapped + all_coms = prod - qml.math.transpose(prod, (2, 1, 0, 3)) + + # 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) + # Normalize trace inner product by dimension chi + adj = qml.math.real(1j * qml.math.tensordot(g / chi, all_coms, axes=[[1, 2], [3, 1]])) + + if is_orthogonal: + # Orthogonal but not normalized inputs. Need to correct by (diagonal) Gram matrix + + if interface == "tensorflow": + import keras # pylint: disable=import-outside-toplevel + + pre_diag = keras.ops.diagonal( + keras.ops.diagonal(prod, axis1=1, axis2=3), axis1=0, axis2=1 + ) + else: + # offset, axis1, axis2 arguments are called differently in torch, use positional arguments + pre_diag = qml.math.diagonal(qml.math.diagonal(prod, 0, 1, 3), 0, 0, 1) + + gram_diag = qml.math.real(qml.math.sum(pre_diag, axis=0)) + + adj = (chi / gram_diag[:, None, None]) * adj + else: + # Non-orthogonal inputs. Need to correct by (full) Gram matrix + # 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 = qml.math.linalg.pinv( + qml.math.real(qml.math.sum(qml.math.diagonal(prod, axis1=1, axis2=3), axis=-1)) + ) + # Axis ordering for contraction with gamma axis of raw structure constants: + # (dimg, _dimg_), (_dimg_, dimg, dimg) -> (dimg, dimg, dim) + # Here we add the missing normalization factor of the trace inner product (after inversion) + adj = qml.math.tensordot(gram_inv * chi, adj, axes=1) + + return adj diff --git a/tests/pauli/dla/test_structure_constants.py b/tests/pauli/dla/test_structure_constants.py index 5fc27fb670a..4b7f871a651 100644 --- a/tests/pauli/dla/test_structure_constants.py +++ b/tests/pauli/dla/test_structure_constants.py @@ -51,8 +51,9 @@ def test_structure_constants_dim(self): @pytest.mark.parametrize("dla", [Ising3, XXZ3]) @pytest.mark.parametrize("assume_orthogonal", [True, False]) @pytest.mark.parametrize("change_norms", [False, True]) + @pytest.mark.parametrize("matrix", [False, True]) def test_structure_constants_elements_with_orthogonal_basis( - self, dla, assume_orthogonal, change_norms + self, dla, assume_orthogonal, change_norms, matrix ): r"""Test relation :math:`[i G_α, i G_β] = \sum_{γ=0}^{d-1} f^γ_{α,β} iG_γ_` with orthogonal bases. @@ -68,7 +69,9 @@ def test_structure_constants_elements_with_orthogonal_basis( # Sample some random new coefficients between 0.5 and 1.5 coeffs = np.random.random(d) + 0.5 dla = [c * op for c, op in zip(coeffs, dla)] - ad_rep = structure_constants(dla, pauli=True, is_orthogonal=assume_orthogonal) + ad_rep = structure_constants( + dla, pauli=True, matrix=matrix, is_orthogonal=assume_orthogonal + ) for alpha in range(d): for beta in range(d): @@ -80,7 +83,8 @@ def test_structure_constants_elements_with_orthogonal_basis( assert all(np.isclose(comm_res[k], res[k]) for k in res) @pytest.mark.parametrize("ortho_dla", [Ising3, XXZ3]) - def test_structure_constants_elements_with_non_orthogonal(self, ortho_dla): + @pytest.mark.parametrize("matrix", [False, True]) + def test_structure_constants_elements_with_non_orthogonal(self, ortho_dla, matrix): r"""Test relation :math:`[i G_α, i G_β] = \sum_{γ=0}^{d-1} f^γ_{α,β} iG_γ_` with non-orthogonal bases. """ @@ -88,7 +92,7 @@ def test_structure_constants_elements_with_non_orthogonal(self, ortho_dla): coeffs = np.random.random((d, d)) + 0.5 dla = [sum(c * op for c, op in zip(_coeffs, ortho_dla)) for _coeffs in coeffs] - ad_rep = structure_constants(dla, pauli=True, is_orthogonal=False) + ad_rep = structure_constants(dla, pauli=True, matrix=matrix, is_orthogonal=False) for alpha in range(d): for beta in range(d): @@ -110,14 +114,24 @@ def test_structure_constants_elements_with_non_orthogonal(self, ortho_dla): assert np.allclose(transf_ortho_ad_rep, ad_rep) @pytest.mark.parametrize("dla", [Ising3, XXZ3]) - def test_use_operators(self, dla): + @pytest.mark.parametrize("matrix", [False, True]) + def test_use_operators(self, dla, matrix): """Test that operators can be passed and lead to the same result""" - ad_rep_true = structure_constants(dla, pauli=True) + ad_rep_true = structure_constants(dla, pauli=True, matrix=matrix) ops = [op.operation() for op in dla] - ad_rep = structure_constants(ops, pauli=False) + ad_rep = structure_constants(ops, pauli=False, matrix=matrix) assert qml.math.allclose(ad_rep, ad_rep_true) + @pytest.mark.parametrize("dla", [Ising3, XXZ3]) + def test_matrix_input(self, dla): + """Test structure constants work as expected for matrix inputs""" + dla_m = [qml.matrix(op, wire_order=range(3)) for op in dla] + adj = qml.structure_constants(dla, matrix=True, is_orthogonal=False) + adj_m = qml.structure_constants(dla_m, matrix=True, is_orthogonal=False) + + assert np.allclose(adj, adj_m) + def test_raise_error_for_non_paulis(self): """Test that an error is raised when passing operators that do not have a pauli_rep""" generators = [qml.Hadamard(0), qml.X(0)] @@ -125,3 +139,47 @@ def test_raise_error_for_non_paulis(self): ValueError, match="Cannot compute adjoint representation of non-pauli operators" ): qml.pauli.structure_constants(generators) + + +dla0 = qml.lie_closure([qml.X(0) @ qml.X(1), qml.Z(0), qml.Z(1)], matrix=True) +adj0 = qml.structure_constants(dla0, matrix=True) + + +class TestInterfacesStructureConstants: + """Test interfaces jax, torch and tensorflow with structure constants""" + + @pytest.mark.jax + def test_jax_structure_constants(self): + """Test jax interface for structure constants""" + + import jax.numpy as jnp + + dla_jax = jnp.array(dla0) + adj_jax = qml.structure_constants(dla_jax, matrix=True) + + assert qml.math.allclose(adj_jax, adj0) + assert qml.math.get_interface(adj_jax) == "jax" + + @pytest.mark.torch + def test_torch_structure_constants(self): + """Test torch interface for structure constants""" + + import torch + + dla_torch = torch.tensor(dla0) + adj_torch = qml.structure_constants(dla_torch, matrix=True) + + assert qml.math.allclose(adj_torch, adj0) + assert qml.math.get_interface(adj_torch) == "torch" + + @pytest.mark.tf + def test_tf_structure_constants(self): + """Test tf interface for structure constants""" + + import tensorflow as tf + + dla_tf = tf.constant(dla0) + adj_tf = qml.structure_constants(dla_tf, matrix=True) + + assert qml.math.allclose(adj_tf, adj0) + assert qml.math.get_interface(adj_tf) == "tensorflow"