diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md index 741b53c7250..d7bf78fd266 100644 --- a/doc/releases/changelog-dev.md +++ b/doc/releases/changelog-dev.md @@ -17,41 +17,41 @@ * Developers of plugin devices now have the option of providing a TOML-formatted configuration file to declare the capabilities of the device. See [Device Capabilities](https://docs.pennylane.ai/en/latest/development/plugins.html#device-capabilities) for details. - * An internal module `pennylane.devices.capabilities` is added that defines a new `DeviceCapabilites` - data class, as well as functions that load and parse the TOML-formatted configuration files. - [(#6407)](https://github.com/PennyLaneAI/pennylane/pull/6407) +* An internal module `pennylane.devices.capabilities` is added that defines a new `DeviceCapabilites` + data class, as well as functions that load and parse the TOML-formatted configuration files. + [(#6407)](https://github.com/PennyLaneAI/pennylane/pull/6407) - ```pycon - >>> from pennylane.devices.capabilities import DeviceCapabilities - >>> capabilities = DeviceCapabilities.from_toml_file("my_device.toml") - >>> isinstance(capabilities, DeviceCapabilities) - True - ``` + ```pycon + >>> from pennylane.devices.capabilities import DeviceCapabilities + >>> capabilities = DeviceCapabilities.from_toml_file("my_device.toml") + >>> isinstance(capabilities, DeviceCapabilities) + True + ``` - * Devices that extends `qml.devices.Device` now has an optional class attribute `capabilities` - that is an instance of the `DeviceCapabilities` data class, constructed from the configuration - file if it exists. Otherwise, it is set to `None`. - [(#6433)](https://github.com/PennyLaneAI/pennylane/pull/6433) +* Devices that extends `qml.devices.Device` now has an optional class attribute `capabilities` + that is an instance of the `DeviceCapabilities` data class, constructed from the configuration + file if it exists. Otherwise, it is set to `None`. + [(#6433)](https://github.com/PennyLaneAI/pennylane/pull/6433) - ```python - from pennylane.devices import Device + ```python + from pennylane.devices import Device - class MyDevice(Device): + class MyDevice(Device): - config_filepath = "path/to/config.toml" + config_filepath = "path/to/config.toml" - ... - ``` - ```pycon - >>> isinstance(MyDevice.capabilities, DeviceCapabilities) - True - ``` + ... + ``` + ```pycon + >>> isinstance(MyDevice.capabilities, DeviceCapabilities) + True + ``` - * Default implementations of `Device.setup_execution_config` and `Device.preprocess_transforms` - are added to the device API for devices that provides a TOML configuration file and thus have - a `capabilities` property. - [(#6632)](https://github.com/PennyLaneAI/pennylane/pull/6632) - [(#6653)](https://github.com/PennyLaneAI/pennylane/pull/6653) +* Default implementations of `Device.setup_execution_config` and `Device.preprocess_transforms` + are added to the device API for devices that provides a TOML configuration file and thus have + a `capabilities` property. + [(#6632)](https://github.com/PennyLaneAI/pennylane/pull/6632) + [(#6653)](https://github.com/PennyLaneAI/pennylane/pull/6653) * Support is added for `if`/`else` statements and `for` and `while` loops in circuits executed with `qml.capture.enabled`, via Autograph. Autograph conversion is now used by default in `make_plxpr`, but can be skipped with the keyword arg `autograph=False`. @@ -60,10 +60,20 @@ [(#6426)](https://github.com/PennyLaneAI/pennylane/pull/6426) [(#6645)](https://github.com/PennyLaneAI/pennylane/pull/6645) - * New `qml.GQSP` template has been added to perform Generalized Quantum Signal Processing (GQSP). +* New `qml.GQSP` template has been added to perform Generalized Quantum Signal Processing (GQSP). The functionality `qml.poly_to_angles` has been also extended to support GQSP. [(#6565)](https://github.com/PennyLaneAI/pennylane/pull/6565) +* Added `unary_mapping()` function to map `BoseWord` and `BoseSentence` to qubit operators, using unary mapping. + [(#6576)](https://github.com/PennyLaneAI/pennylane/pull/6576) + +* Added `binary_mapping()` function to map `BoseWord` and `BoseSentence` to qubit operators, using standard-binary mapping. + [(#6564)](https://github.com/PennyLaneAI/pennylane/pull/6564) + +* New functionality to calculate angles for QSP and QSVT has been added. This includes the function `qml.poly_to_angles` + to obtain angles directly and the function `qml.transform_angles` to convert angles from one subroutine to another. + [(#6483)](https://github.com/PennyLaneAI/pennylane/pull/6483) + * Added a function `qml.trotterize` to generalize the Suzuki-Trotter product to arbitrary quantum functions. [(#6627)](https://github.com/PennyLaneAI/pennylane/pull/6627) @@ -96,32 +106,31 @@ b: ──RY(-0.17)─╰X─╰X──RY(-0.17)─┤ State ``` -

New `labs` module `dla` for handling dynamical Lie algebras (DLAs)

+

New `pennylane.labs.dla` module for handling (dynamical) Lie algebras (DLAs)

* Added a dense implementation of computing the Lie closure in a new function `lie_closure_dense` in `pennylane.labs.dla`. [(#6371)](https://github.com/PennyLaneAI/pennylane/pull/6371) [(#6695)](https://github.com/PennyLaneAI/pennylane/pull/6695) -* New functionality to calculate angles for QSP and QSVT has been added. This includes the function `qml.poly_to_angles` - to obtain angles directly and the function `qml.transform_angles` to convert angles from one subroutine to another. - [(#6483)](https://github.com/PennyLaneAI/pennylane/pull/6483) - * Added a dense implementation of computing the structure constants in a new function `structure_constants_dense` in `pennylane.labs.dla`. [(#6376)](https://github.com/PennyLaneAI/pennylane/pull/6376) -* Added utility functions for handling dense matrices in the Lie theory context. +* Added utility functions for handling dense matrices and advanced functionality in the Lie theory context. [(#6563)](https://github.com/PennyLaneAI/pennylane/pull/6563) [(#6392)](https://github.com/PennyLaneAI/pennylane/pull/6392) + [(#6396)](https://github.com/PennyLaneAI/pennylane/pull/6396) * Added a ``cartan_decomp`` function along with two standard involutions ``even_odd_involution`` and ``concurrence_involution``. [(#6392)](https://github.com/PennyLaneAI/pennylane/pull/6392) -* Added `unary_mapping()` function to map `BoseWord` and `BoseSentence` to qubit operators, using unary mapping - [(#6576)](https://github.com/PennyLaneAI/pennylane/pull/6576); -added `binary_mapping()` function to map `BoseWord` and `BoseSentence` to qubit operators, using standard-binary mapping. - [(#6564)](https://github.com/PennyLaneAI/pennylane/pull/6564) +* Added a `recursive_cartan_decomp` function and all canonical Cartan involutions. + [(#6396)](https://github.com/PennyLaneAI/pennylane/pull/6396) + +* Added a `cartan_subalgebra` function to compute the (horizontal) Cartan subalgebra of a Cartan decomposition. + [(#6403)](https://github.com/PennyLaneAI/pennylane/pull/6403) + [(#6396)](https://github.com/PennyLaneAI/pennylane/pull/6396)

New API for Qubit Mixed

diff --git a/pennylane/labs/dla/__init__.py b/pennylane/labs/dla/__init__.py index eb4c4ccbd76..de12d7cbadc 100644 --- a/pennylane/labs/dla/__init__.py +++ b/pennylane/labs/dla/__init__.py @@ -23,6 +23,8 @@ ~lie_closure_dense ~structure_constants_dense ~cartan_decomp + ~recursive_cartan_decomp + ~cartan_subalgebra Utility functions @@ -34,20 +36,28 @@ :toctree: api ~adjvec_to_op - ~change_basis_ad_rep - ~check_orthonormal - ~check_commutation - ~check_all_commuting - ~check_cartan_decomp ~op_to_adjvec + ~trace_inner_product ~orthonormalize ~pauli_coefficients ~batched_pauli_decompose - ~trace_inner_product + ~check_orthonormal + ~check_commutation + ~check_all_commuting + ~check_cartan_decomp + ~change_basis_ad_rep Involutions ~~~~~~~~~~~ +A map :math:`\theta: \mathfrak{g} \rightarrow \mathfrak{g}` from the Lie algebra :math:`\mathfrak{g}` to itself is called an involution +when it fulfills :math:`\theta(\theta(g)) = g \ \forall g \in \mathfrak{g}` and is compatible with commutators, +:math:`[\theta(g), \theta(g')]=\theta([g, g']).` Involutions are used to construct a :func:`~cartan_decomp`. There are seven canonical +Cartan involutions of real simple Lie algebras (``AI, AII, AIII, BDI, CI, CII, DIII``), +see `Wikipedia `__. +In addition, there is a canonical Cartan involution for real semisimple algebras that consist of +two isomorphic simple components (``ClassB``), see `here `__. + .. currentmodule:: pennylane.labs.dla .. autosummary:: @@ -55,23 +65,51 @@ ~even_odd_involution ~concurrence_involution + ~khaneja_glaser_involution + ~AI + ~AII + ~AIII + ~BDI + ~CI + ~CII + ~DIII + ~ClassB """ from .lie_closure_dense import lie_closure_dense from .structure_constants_dense import structure_constants_dense -from .cartan import cartan_decomp, even_odd_involution, concurrence_involution +from .cartan import ( + cartan_decomp, + recursive_cartan_decomp, +) from .dense_util import ( + adjvec_to_op, change_basis_ad_rep, + check_all_commuting, + check_cartan_decomp, + check_commutation, + check_orthonormal, pauli_coefficients, batched_pauli_decompose, - orthonormalize, - check_orthonormal, trace_inner_product, - adjvec_to_op, op_to_adjvec, - check_commutation, - check_all_commuting, - check_cartan_decomp, + orthonormalize, ) + +from .involutions import ( + khaneja_glaser_involution, + even_odd_involution, + concurrence_involution, + AI, + AII, + AIII, + BDI, + CI, + CII, + DIII, + ClassB, +) + +from .cartan_subalgebra import cartan_subalgebra diff --git a/pennylane/labs/dla/cartan.py b/pennylane/labs/dla/cartan.py index 1bc1af3ee60..097006b480f 100644 --- a/pennylane/labs/dla/cartan.py +++ b/pennylane/labs/dla/cartan.py @@ -12,16 +12,19 @@ # See the License for the specific language governing permissions and # limitations under the License. """Functionality for Cartan decomposition""" -from functools import singledispatch +# pylint: disable= missing-function-docstring +from functools import partial from typing import List, Tuple, Union import numpy as np -import pennylane as qml -from pennylane import Y +from pennylane import QubitUnitary from pennylane.operation import Operator from pennylane.pauli import PauliSentence +from .dense_util import check_cartan_decomp +from .involutions import int_log2 + def cartan_decomp( g: List[Union[PauliSentence, Operator]], involution: callable @@ -47,8 +50,8 @@ def cartan_decomp( Returns: Tuple(List[Union[PauliSentence, Operator]], List[Union[PauliSentence, Operator]]): Tuple ``(k, m)`` containing the even - parity subspace :math:`\Theta(\mathfrak{k}) = \mathfrak{k}` and the odd - parity subspace :math:`\Theta(\mathfrak{m}) = -\mathfrak{m}`. + parity subspace :math:`\Theta(\mathfrak{k}) = \mathfrak{k}` and the odd + parity subspace :math:`\Theta(\mathfrak{m}) = -\mathfrak{m}`. .. seealso:: :func:`~even_odd_involution`, :func:`~concurrence_involution`, :func:`~check_cartan_decomp` @@ -99,155 +102,246 @@ def cartan_decomp( k = [] for op in g: - if involution(op): # odd parity + if involution(op): # odd parity theta(k) = k k.append(op) - else: # even parity + else: # even parity theta(m) = -m m.append(op) return k, m -# dispatch to different input types -def even_odd_involution(op: Union[PauliSentence, np.ndarray, Operator]) -> bool: - r"""The Even-Odd involution - - This is defined in `quant-ph/0701193 `__. - For Pauli words and sentences, it comes down to counting non-trivial Paulis in Pauli words. - - Args: - op ( Union[PauliSentence, np.ndarray, Operator]): Input operator - - Returns: - bool: Boolean output ``True`` or ``False`` for odd (:math:`\mathfrak{k}`) and even parity subspace (:math:`\mathfrak{m}`), respectively - - .. seealso:: :func:`~cartan_decomp` - - **Example** - - >>> from pennylane import X, Y, Z - >>> from pennylane.labs.dla import even_odd_involution - >>> ops = [X(0), X(0) @ Y(1), X(0) @ Y(1) @ Z(2)] - >>> [even_odd_involution(op) for op in ops] - [True, False, True] - - Operators with an odd number of non-identity Paulis yield ``1``, whereas even ones yield ``0``. - - The function also works with dense matrix representations. - - >>> ops_m = [qml.matrix(op, wire_order=range(3)) for op in ops] - >>> [even_odd_involution(op_m) for op_m in ops_m] - [True, False, True] - - """ - return _even_odd_involution(op) - - -@singledispatch -def _even_odd_involution(op): # pylint:disable=unused-argument - return NotImplementedError(f"Involution not defined for operator {op} of type {type(op)}") - - -@_even_odd_involution.register(PauliSentence) -def _even_odd_involution_ps(op: PauliSentence): - # Generalization to sums of Paulis: check each term and assert they all have the same parity - parity = [] - for pw in op.keys(): - parity.append(len(pw) % 2) - - # only makes sense if parity is the same for all terms, e.g. Heisenberg model - assert all( - parity[0] == p for p in parity - ), f"The Even-Odd involution is not well-defined for operator {op} as individual terms have different parity" - return bool(parity[0]) - - -@_even_odd_involution.register(np.ndarray) -def _even_odd_involution_matrix(op: np.ndarray): - """see Table CI in https://arxiv.org/abs/2406.04418""" - n = int(np.round(np.log2(op.shape[-1]))) - YYY = qml.prod(*[Y(i) for i in range(n)]) - YYY = qml.matrix(YYY, range(n)) - - transformed = YYY @ op.conj() @ YYY - if np.allclose(transformed, op): - return False - if np.allclose(transformed, -op): - return True - raise ValueError(f"The Even-Odd involution is not well-defined for operator {op}.") - - -@_even_odd_involution.register(Operator) -def _even_odd_involution_op(op: Operator): - """use pauli representation""" - return _even_odd_involution_ps(op.pauli_rep) - +IDENTITY = object() + + +def pauli_y_eigenbasis(wire, num_wires): + V = np.array([[1, 1], [1j, -1j]]) / np.sqrt(2) + return QubitUnitary(V, wire).matrix(wire_order=range(num_wires)) + + +_basis_change_constructors = { + ("AI", "BDI"): IDENTITY, + ("AI", "DIII"): IDENTITY, + ("AII", "CI"): IDENTITY, + ("AII", "CII"): IDENTITY, + ("AIII", "ClassB"): IDENTITY, + ("BDI", "ClassB"): IDENTITY, + ("CI", "AI"): pauli_y_eigenbasis, + ("CI", "AII"): pauli_y_eigenbasis, + ("CI", "AIII"): IDENTITY, + ("CII", "ClassB"): IDENTITY, + ("DIII", "AI"): pauli_y_eigenbasis, + ("DIII", "AII"): pauli_y_eigenbasis, + ("DIII", "AIII"): pauli_y_eigenbasis, + ("ClassB", "AI"): IDENTITY, + ("ClassB", "AII"): IDENTITY, + ("ClassB", "AIII"): IDENTITY, + ("ClassB", "BDI"): IDENTITY, + ("ClassB", "DIII"): IDENTITY, + ("ClassB", "CI"): IDENTITY, + ("ClassB", "CII"): IDENTITY, +} + + +def _check_classb_sequence(before, after): + if before == "AIII" and after.startswith("A"): + return + if before == "BDI" and after in ("BDI", "DIII"): + return + if before == "CII" and after.startswith("C"): + return + raise ValueError( + f"The 3-sequence ({before}, ClassB, {after}) of involutions is not a valid sequence." + ) -# dispatch to different input types -def concurrence_involution(op: Union[PauliSentence, np.ndarray, Operator]) -> bool: - r"""The Concurrence Canonical Decomposition :math:`\Theta(g) = -g^T` as a Cartan involution function - This is defined in `quant-ph/0701193 `__. - For Pauli words and sentences, it comes down to counting Pauli-Y operators. +def _check_chain(chain, num_wires): + """Validate a chain of involutions for a recursive Cartan decomposition.""" + # Take the function name or its `func` attribute if it exists (e.g., for `partial` of an involution) + names = [getattr(phi, "func", phi).__name__ for phi in chain] + + # Assume some standard behaviour regarding the wires on which we need to perform basis changes + basis_changes = [] + wire = 0 + for i, name in enumerate(names[:-1]): + invol_pair = (name, names[i + 1]) + if invol_pair not in _basis_change_constructors: + raise ValueError( + f"The specified chain contains the pair {'-->'.join(invol_pair)}, " + "which is not a valid pair." + ) + # Run specific check for sequence of three involutions where ClassB is the middle one + if name == "ClassB" and i > 0: + _check_classb_sequence(names[i - 1], names[i + 1]) + bc_constructor = _basis_change_constructors[invol_pair] + if bc_constructor is IDENTITY: + bc = bc_constructor + else: + bc = bc_constructor(wire, num_wires) + # Next assumption: The wire is only incremented if a basis change is applied. + wire += 1 + basis_changes.append(bc) + + basis_changes.append(IDENTITY) # Do not perform any basis change after last involution. + return names, basis_changes + + +def _apply_basis_change(change_op, targets): + """Helper function for recursive Cartan decompositions that applies a basis change matrix + ``change_op`` to a batch of matrices ``targets`` (with leading batch dimension).""" + # Compute x V^\dagger for all x in ``targets``. ``moveaxis`` brings the batch axis to the front + out = np.moveaxis(np.tensordot(change_op, targets, axes=[[1], [1]]), 1, 0) + out = np.tensordot(out, change_op.conj().T, axes=[[2], [0]]) + return out + + +def recursive_cartan_decomp(g, chain, validate=True, verbose=True): + r"""Apply a recursive Cartan decomposition specified by a chain of decomposition types. + The decompositions will use canonical involutions and hardcoded basis transformations + between them in order to obtain a valid recursion. + + This function tries to make sure that only sensible involution sequences are applied, + and to raise an error otherwise. However, the involutions still need to be configured + properly, regarding the wires their conjugation operators act on. Args: - op ( Union[PauliSentence, np.ndarray, Operator]): Input operator + g (tensor_like): Basis of the algebra to be decomposed. + chain (Iterable[Callable]): Sequence of involutions. Each callable should be + one of + :func:`~.pennylane.labs.dla.AI`, + :func:`~.pennylane.labs.dla.AII`, + :func:`~.pennylane.labs.dla.AIII`, + :func:`~.pennylane.labs.dla.BDI`, + :func:`~.pennylane.labs.dla.CI + :func:`~.pennylane.labs.dla.CII`, + :func:`~.pennylane.labs.dla.DIII`, or + :func:`~.pennylane.labs.dla.ClassB`, + or a partial evolution thereof. + validate (bool): Whether or not to verify that the involutions return a subalgebra. + verbose (bool): Whether or not to print status updates during the computation. Returns: - bool: Boolean output ``True`` or ``False`` for odd (:math:`\mathfrak{k}`) and even parity subspace (:math:`\mathfrak{m}`), respectively - - .. seealso:: :func:`~cartan_decomp` - - **Example** - - >>> from pennylane import X, Y, Z - >>> from pennylane.labs.dla import concurrence_involution - >>> ops = [X(0), X(0) @ Y(1), X(0) @ Y(1) @ Z(2), Y(0) @ Y(2)] - >>> [concurrence_involution(op) for op in ops] - [False, True, True, False] - - Operators with an odd number of ``Y`` operators yield ``1``, whereas even ones yield ``0``. - - The function also works with dense matrix representations. - - >>> ops_m = [qml.matrix(op, wire_order=range(3)) for op in ops] - >>> [even_odd_involution(op_m) for op_m in ops_m] - [False, True, True, False] + dict: The decompositions at each level. The keys are (zero-based) integers for the + different levels of the recursion, the values are tuples ``(k, m)`` with subalgebra + ``k`` and horizontal space ``m``. For each level, ``k`` and ``m`` combine into + ``k`` from the previous recursion level. + + **Examples** + + Let's set up the special unitary algebra on 2 qubits. Note that we are using the Hermitian + matrices that correspond to the skew-Hermitian algebra elements via multiplication + by :math:`i`. Also note that :func:`~.pauli.pauli_group` returns the identity as the first + element, which is not part of the special unitary algebra of traceless matrices. + + >>> g = [qml.matrix(op, wire_order=range(2)) for op in qml.pauli.pauli_group(2)] # u(4) + >>> g = g[1:] # Remove identity: u(4) -> su(4) + + Now we can apply Cartan decompositions of type AI and DIII in sequence: + + >>> from pennylane.labs.dla import recursive_cartan_decomp, AI, DIII + >>> chain = [AI, DIII] + >>> decompositions = recursive_cartan_decomp(g, chain) + Iteration 0: 15 -----AI----> 6, 9 + Iteration 1: 6 ----DIII---> 4, 2 + + The function prints the progress of the decompositions by default, which can be deactivated by + setting ``verbose=False``. Here we see how the initial :math:`\mathfrak{g}=\mathfrak{su(4)}` + was decomposed by AI into the six-dimensional :math:`\mathfrak{k}_1=\mathfrak{so(4)}` and a + horizontal space of dimension nine. Then, :math:`\mathfrak{k}_1` was further decomposed + by the DIII decomposition into the four-dimensional :math:`\mathfrak{k}_2=\mathfrak{u}(2)` + and a two-dimensional horizontal space. + + In a more elaborate example, let's apply a chain of decompositions AII, CI, AI, BDI, and DIII + to the four-qubit unitary algebra. While we took care of the global phase term of :math:`u(4)` + explicitly above, we leave it in the algebra here, and see that it does not cause problems. + We discuss the ``wire`` keyword argument below. + + >>> from pennylane.labs.dla import AII, CI, BDI, ClassB + >>> from functools import partial + >>> chain = [ + ... AII, + ... CI, + ... AI, + ... partial(BDI, wire=1), + ... partial(ClassB, wire=1), + ... partial(DIII, wire=2), + ... ] + >>> g = [qml.matrix(op, wire_order=range(4)) for op in qml.pauli.pauli_group(4)] # u(16) + >>> decompositions = recursive_cartan_decomp(g, chain) + Iteration 0: 256 ----AII----> 136, 120 + Iteration 1: 136 -----CI----> 64, 72 + Iteration 2: 64 -----AI----> 28, 36 + Iteration 3: 28 ----BDI----> 12, 16 + Iteration 4: 12 ---ClassB--> 6, 6 + Iteration 5: 6 ----DIII---> 4, 2 + + The obtained chain of algebras is + + .. math:: + + \mathfrak{u}(16) + \rightarrow \mathfrak{sp}(8) + \rightarrow \mathfrak{u}(8) + \rightarrow \mathfrak{so}(8) + \rightarrow \mathfrak{so}(4) \oplus \mathfrak{so}(4) + \rightarrow \mathfrak{so}(4) + \rightarrow \mathfrak{u}(2). + + What about the wire keyword argument to the used involutions? + A good rule of thumb is that it should start at ``0`` and increment by one every second + involution. For the involution :func:`~.pennylane.labs.dla.CI` it should additionally be + increased by one. As ``0`` is the default for ``wire``, it usually does not have to be + provided explicitly for the first two involutions, unless ``CI`` is among them. + + .. note:: + + A typical effect of setting the wire wrongly is that the decomposition does not + split the subalgebra from the previous step but keeps it intact and returns a + zero-dimensional horizontal space. For example: + + >>> g = [qml.matrix(op, wire_order=range(2)) for op in qml.pauli.pauli_group(2)] # u(4) + >>> chain = [AI, DIII, AII] + >>> decompositions = recursive_cartan_decomp(g, chain) + Iteration 0: 16 -----AI----> 6, 10 + Iteration 1: 6 ----DIII---> 4, 2 + Iteration 2: 4 ----AII----> 4, 0 + + We see that the ``AII`` decomposition did not further decompose :math:`\mathfrak{u}(2)`. + It works if we provide the correct ``wire`` argument: + + >>> chain = [AI, DIII, partial(AII, wire=1)] + >>> decompositions = recursive_cartan_decomp(g, chain) + Iteration 0: 16 -----AI----> 6, 10 + Iteration 1: 6 ----DIII---> 4, 2 + Iteration 2: 4 ----AII----> 3, 1 + + We obtain :math:`\mathfrak{sp}(1)` as expected from the decomposition of type AII. """ - return _concurrence_involution(op) - - -@singledispatch -def _concurrence_involution(op): - return NotImplementedError(f"Involution not defined for operator {op} of type {type(op)}") - - -@_concurrence_involution.register(PauliSentence) -def _concurrence_involution_pauli(op: PauliSentence): - # Generalization to sums of Paulis: check each term and assert they all have the same parity - parity = [] - for pw in op.keys(): - result = sum(1 if el == "Y" else 0 for el in pw.values()) - parity.append(result % 2) - - # only makes sense if parity is the same for all terms, e.g. Heisenberg model - assert all( - parity[0] == p for p in parity - ), f"The concurrence canonical decomposition is not well-defined for operator {op} as individual terms have different parity" - return bool(parity[0]) - -@_concurrence_involution.register(Operator) -def _concurrence_involution_operator(op: Operator): - return _concurrence_involution_matrix(op.matrix()) - - -@_concurrence_involution.register(np.ndarray) -def _concurrence_involution_matrix(op: np.ndarray): - if np.allclose(op, -op.T): - return True - if np.allclose(op, op.T): - return False - raise ValueError( - f"The concurrence canonical decomposition is not well-defined for operator {op}" - ) + # Prerun the validation by obtaining the required basis changes and raising an error if + # an invalid pair is found. + num_wires = int_log2(np.shape(g)[-1]) + names, basis_changes = _check_chain(chain, num_wires) + + decompositions = {} + for i, (phi, name, bc) in enumerate(zip(chain, names, basis_changes)): + try: + k, m = cartan_decomp(g, phi) + except ValueError as e: + if "please specify p and q for the involution" in str(e): + phi = partial(phi, p=2 ** (num_wires - 1), q=2 ** (num_wires - 1)) + k, m = cartan_decomp(g, phi) + else: + raise ValueError from e + + if validate: + check_cartan_decomp(k, m, verbose=verbose) + if verbose: + print(f"Iteration {i}: {len(g):>4} -{name:-^10}> {len(k):>4},{len(m):>4}") + decompositions[i] = (k, m) + if not bc is IDENTITY: + k = _apply_basis_change(bc, k) + m = _apply_basis_change(bc, m) + g = k + + return decompositions diff --git a/pennylane/labs/dla/cartan_subalgebra.py b/pennylane/labs/dla/cartan_subalgebra.py new file mode 100644 index 00000000000..cc97dceb0f7 --- /dev/null +++ b/pennylane/labs/dla/cartan_subalgebra.py @@ -0,0 +1,271 @@ +# 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. +"""Functionality to compute the Cartan subalgebra""" +# pylint: disable=too-many-arguments, too-many-positional-arguments +import copy + +import numpy as np +from scipy.linalg import null_space + +import pennylane as qml +from pennylane.pauli.dla.center import _intersect_bases + +from .dense_util import adjvec_to_op, change_basis_ad_rep, op_to_adjvec + + +def _gram_schmidt(X): + """Orthogonalize basis of column vectors in X""" + Q, _ = np.linalg.qr(X.T, mode="reduced") + return Q.T + + +def _is_independent(v, A, tol=1e-14): + """Check whether ``v`` is independent of the columns of A.""" + v /= np.linalg.norm(v) + v = v - A @ qml.math.linalg.solve(A.conj().T @ A, A.conj().T) @ v + return np.linalg.norm(v) > tol + + +def _orthogonal_complement_basis(h, m, tol): + """find mtilde = m - h""" + # Step 1: Find the span of h + h = np.array(h) + m = np.array(m) + + # Compute the orthonormal basis of h using QR decomposition + + Q = _gram_schmidt(h) + + # Step 2: Project each vector in m onto the orthogonal complement of span(h) + projections = m - np.dot(np.dot(m, Q.T), Q) + assert np.allclose( + np.tensordot(h, projections, axes=[[1], [1]]), 0.0 + ), f"{np.tensordot(h, projections, axes=[[1], [1]])}" + + # Step 3: Find a basis for the non-zero projections + # We'll use SVD to find the basis + U, S, _ = np.linalg.svd(projections.T) + + # Choose columns of U corresponding to non-zero singular values + rank = np.sum(S > tol) + basis = U[:, :rank] + assert np.allclose( + np.tensordot(h, basis, axes=[[1], [0]]), 0.0 + ), f"{np.tensordot(h, basis, axes=[[1], [0]])}" + + return basis.T # Transpose to get row vectors + + +def cartan_subalgebra( + g, k, m, ad, start_idx=0, tol=1e-10, verbose=0, return_adjvec=False, is_orthogonal=True +): + r""" + Compute a Cartan subalgebra (CSA) :math:`\mathfrak{h} \subseteq \mathfrak{m}`. + + A non-unique CSA is a maximal Abelian subalgebra in the horizontal subspace :math:`\mathfrak{m}` of a Cartan decomposition. + Note that this is sometimes called a horizontal CSA, and is different from the definition of a CSA on `Wikipedia `__. + + .. seealso:: :func:`~cartan_decomp`, :func:`~structure_constants`, `The KAK theorem (demo) `__ + + Args: + g (List[Union[PauliSentence, np.ndarray]]): Lie algebra :math:`\mathfrak{g}`, which is assumed to be ordered as :math:`\mathfrak{g} = \mathfrak{k} \oplus \mathfrak{m}` + k (List[Union[PauliSentence, np.ndarray]]): Vertical space :math:`\mathfrak{k}` from Cartan decomposition :math:`\mathfrak{g} = \mathfrak{k} \oplus \mathfrak{m}` + m (List[Union[PauliSentence, np.ndarray]]): Horizontal space :math:`\mathfrak{m}` from Cartan decomposition :math:`\mathfrak{g} = \mathfrak{k} \oplus \mathfrak{m}` + ad (Array): The :math:`|\mathfrak{g}| \times |\mathfrak{g}| \times |\mathfrak{g}|` dimensional adjoint representation of :math:`\mathfrak{g}` (see :func:`~structure_constants`) + start_idx (bool): Indicates from which element in ``m`` the CSA computation starts. + tol (float): Numerical tolerance for linear independence check + verbose (bool): Whether or not to output progress during computation + return_adjvec (bool): The output format. If ``False``, returns operators in their original + input format (matrices or :class:`~PauliSentence`). If ``True``, returns the spaces as adjoint representation vectors. + is_orthogonal (bool): Whether the basis elements are all orthogonal, both within + and between ``g``, ``k`` and ``m``. + + Returns: + Tuple(np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray): A tuple of adjoint vector representations + ``(newg, k, mtilde, h, new_adj)``, corresponding to + :math:`\mathfrak{g}`, :math:`\mathfrak{k}`, :math:`\tilde{\mathfrak{m}}`, :math:`\mathfrak{h}` and the new adjoint representation. + The dimensions are ``(|g|, |g|)``, ``(|k|, |g|)``, ``(|mtilde|, |g|)``, ``(|h|, |g|)`` and ``(|g|, |g|, |g|)``, respectively. + + **Example** + + A quick example computing a Cartan subalgebra of :math:`\mathfrak{su}(4)` using the Cartan involution :func:`~even_odd_involution`. + + >>> import pennylane as qml + >>> from pennylane.labs.dla import cartan_decomp, cartan_subalgebra, even_odd_involution + >>> g = list(qml.pauli.pauli_group(2)) # u(4) + >>> g = g[1:] # remove identity -> su(4) + >>> g = [op.pauli_rep for op in g] # optional; turn to PauliSentence for convenience + >>> k, m = cartan_decomp(g, even_odd_involution) + >>> g = k + m # re-order g to separate k and m + >>> adj = qml.structure_constants(g) + >>> newg, k, mtilde, h, new_adj = cartan_subalgebra(g, k, m, adj) + >>> newg == k + mtilde + h + True + >>> h + [-1.0 * Z(0) @ Z(1), -1.0 * Y(0) @ Y(1), 1.0 * X(0) @ X(1)] + + We can confirm that these all commute with each other, as the CSA is Abelian (= all operators commute). + + >>> from pennylane.labs.dla import check_all_commuting + >>> check_all_commuting(h) + True + + We can opt-in to return what we call adjoint vectors of dimension :math:`|\mathfrak{g}|`, where each component corresponds to an entry in (the ordered) ``g``. + The adjoint vectors for the Cartan subalgebra are in ``np_h``. + + >>> np_newg, np_k, np_mtilde, np_h, new_adj = cartan_subalgebra(g, k, m, adj, return_adjvec=True) + + We can reconstruct an operator by computing :math:`\hat{O}_v = \sum_i v_i g_i` for an adjoint vector :math:`v` and :math:`g_i \in \mathfrak{g}`. + + >>> v = np_h[0] + >>> op = sum(v_i * g_i for v_i, g_i in zip(v, g)) + >>> op.simplify() + >>> op + -1.0 * Z(0) @ Z(1) + + For convenience, we provide a helper function :func:`~adjvec_to_op` for the collections of adjoint vectors in the returns. + + >>> from pennylane.labs.dla import adjvec_to_op + >>> h = adjvec_to_op(np_h, g) + >>> h + [-1.0 * Z(0) @ Z(1), -1.0 * Y(0) @ Y(1), 1.0 * X(0) @ X(1)] + + .. details:: + :title: Usage Details + + Let us walk through an example of computing the Cartan subalgebra. The basis for computing + the Cartan subalgebra is having the Lie algebra :math:`\mathfrak{g}`, a Cartan decomposition + :math:`\mathfrak{g} = \mathfrak{k} \oplus \mathfrak{m}` and its adjoint representation. + + We start by computing these ingredients using :func:`~cartan_decomp` and :func:`~structure_constants`. + As an example, we take the Lie algebra of the Heisenberg model with generators :math:`\{X_i X_{i+1}, Y_i Y_{i+1}, Z_i Z_{i+1}\}`. + + >>> from pennylane.labs.dla import lie_closure_dense, cartan_decomp + >>> from pennylane import X, Y, Z + >>> n = 3 + >>> gens = [X(i) @ X(i+1) for i in range(n-1)] + >>> gens += [Y(i) @ Y(i+1) for i in range(n-1)] + >>> gens += [Z(i) @ Z(i+1) for i in range(n-1)] + >>> g = lie_closure_dense(gens) + + Taking the Heisenberg Lie algebra, we can perform the Cartan decomposition. We take the :func:`~even_odd_involution` as a valid Cartan involution. + The resulting vertical and horizontal subspaces :math:`\mathfrak{k}` and :math:`\mathfrak{m}` need to fulfill the commutation relations + :math:`[\mathfrak{k}, \mathfrak{k}] \subseteq \mathfrak{k}`, :math:`[\mathfrak{k}, \mathfrak{m}] \subseteq \mathfrak{m}` and :math:`[\mathfrak{m}, \mathfrak{m}] \subseteq \mathfrak{k}`, + which we can check using the helper function :func:`~check_cartan_decomp`. + + >>> from pennylane.labs.dla import even_odd_involution, check_cartan_decomp + >>> k, m = cartan_decomp(g, even_odd_involution) + >>> check_cartan_decomp(k, m) # check commutation relations of Cartan decomposition + True + + Our life is easier when we use a canonical ordering of the operators. This is why we re-define ``g`` with the new ordering in terms of operators in ``k`` first, and then + 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 + + Finally, we can compute a Cartan subalgebra :math:`\mathfrak{h}`, a maximal Abelian subalgebra of :math:`\mathfrak{m}`. + + >>> newg, k, mtilde, h, new_adj = cartan_subalgebra(g, k, m, adj, start_idx=3) + + The new DLA ``newg`` is just the concatenation of ``k``, ``mtilde``, ``h``. Each component is returned in the original input format. + Here we obtain collections of :math:`8\times 8` matrices (``numpy`` arrays), as this is what we started from. + + >>> newg.shape, k.shape, mtilde.shape, h.shape, new_adj.shape + ((15, 8, 8), (6, 8, 8), (6, 8, 8), (3, 8, 8), (15, 15, 15)) + + We can also let the function return what we call adjoint representation vectors. + + >>> kwargs = {"start_idx": 3, "return_adjvec": True} + >>> np_newg, np_k, np_mtilde, np_h, new_adj = cartan_subalgebra(g, k, m, adj, **kwargs) + >>> np_newg.shape, np_k.shape, np_mtilde.shape, np_h.shape, new_adj.shape + ((15, 15), (6, 15), (6, 15), (3, 15), (15, 15, 15)) + + These are dense vector representations of dimension :math:`|\mathfrak{g}|`, in which each entry corresponds to the respective operator in :math:`\mathfrak{g}`. + Given an adjoint representation vector :math:`v`, we can reconstruct the respective operator by simply computing :math:`\hat{O}_v = \sum_i v_i g_i`, where + :math:`g_i \in \mathfrak{g}` (hence the need for a canonical ordering). + + We provide a convenience function :func:`~adjvec_to_op` that works with both ``g`` represented as dense matrices or PL operators. + Because we used dense matrices in this example, we transform the operators back to PennyLane operators using :func:`~pauli_decompose`. + + >>> from pennylane.labs.dla import adjvec_to_op + >>> h = adjvec_to_op(np_h, g) + >>> h_op = [qml.pauli_decompose(op).pauli_rep for op in h] + >>> h_op + [-1.0 * Y(1) @ Y(2), -1.0 * Z(1) @ Z(2), 1.0 * X(1) @ X(2)] + + In that case we chose a Cartan subalgebra from which we can readily see that it is commuting, but we also provide a small helper function to check that. + + >>> from pennylane.labs.dla import check_all_commuting + >>> assert check_all_commuting(h_op) + + Last but not least, the adjoint representation ``new_adj`` is updated to represent the new basis and its ordering of ``g``. + """ + + g_copy = copy.deepcopy(g) + np_m = op_to_adjvec(m, g, is_orthogonal=is_orthogonal) + np_h = op_to_adjvec([m[start_idx]], g, is_orthogonal=is_orthogonal) + + iteration = 1 + while True: + if verbose: + print(f"iteration {iteration}: Found {len(np_h)} independent Abelian operators.") + kernel_intersection = np_m + for h_i in np_h: + + # obtain adjoint rep of candidate h_i + adjoint_of_h_i = np.tensordot(ad, h_i, axes=[[1], [0]]) + # compute kernel of adjoint + new_kernel = null_space(adjoint_of_h_i, rcond=tol) + + # intersect kernel to stay in m + kernel_intersection = _intersect_bases(kernel_intersection.T, new_kernel, rcond=tol).T + + kernel_intersection = _gram_schmidt(kernel_intersection) # orthogonalize + for vec in kernel_intersection: + if _is_independent(vec, np.array(np_h).T, tol): + np_h = np.vstack([np_h, vec]) + break + else: + # No new vector was added from all the kernels + break + + iteration += 1 + + np_h = _gram_schmidt(np_h) # orthogonalize Abelian subalgebra + np_k = op_to_adjvec( + k, g, is_orthogonal=is_orthogonal + ) # adjoint vectors of k space for re-ordering + np_oldg = np.vstack([np_k, np_m]) + np_k = _gram_schmidt(np_k) + + np_mtilde = _orthogonal_complement_basis(np_h, np_m, tol=tol) # the "rest" of m without h + np_newg = np.vstack([np_k, np_mtilde, np_h]) + + # Instead of recomputing the adjoint representation, take the basis transformation + # oldg -> newg and transform the adjoint representation accordingly + basis_change = np.tensordot(np_newg, np.linalg.pinv(np_oldg), axes=[[1], [0]]) + new_adj = change_basis_ad_rep(ad, basis_change) + + if return_adjvec: + return np_newg, np_k, np_mtilde, np_h, new_adj + + newg, k, mtilde, h = [ + adjvec_to_op(adjvec, g_copy, is_orthogonal=is_orthogonal) + for adjvec in [np_newg, np_k, np_mtilde, np_h] + ] + + return newg, k, mtilde, h, new_adj diff --git a/pennylane/labs/dla/dense_util.py b/pennylane/labs/dla/dense_util.py index c73b93fa2b1..58836046709 100644 --- a/pennylane/labs/dla/dense_util.py +++ b/pennylane/labs/dla/dense_util.py @@ -12,7 +12,7 @@ # See the License for the specific language governing permissions and # limitations under the License. """Utility tools for dense Lie algebra representations""" -# pylint: disable=possibly-used-before-assignment, too-many-return-statements +# pylint: disable=too-many-return-statements, missing-function-docstring, possibly-used-before-assignment from functools import reduce from itertools import combinations, combinations_with_replacement from typing import Iterable, List, Optional, Union @@ -163,7 +163,8 @@ def _idx_to_pw(idx, n): def batched_pauli_decompose(H: TensorLike, tol: Optional[float] = None, pauli: bool = False): - r"""Decomposes a Hermitian matrix into a linear combination of Pauli operators. + r"""Decomposes a Hermitian matrix or a batch of matrices into a linear combination + of Pauli operators. Args: H (tensor_like[complex]): a Hermitian matrix of dimension ``(2**n, 2**n)`` or a collection @@ -235,7 +236,7 @@ def batched_pauli_decompose(H: TensorLike, tol: Optional[float] = None, pauli: b def check_commutation(ops1, ops2, vspace): - r"""Helper function to check :math:`[\text{ops1}, \text{ops2}] \subseteq \text{vspace}` + r"""Helper function to check :math:`[\text{ops1}, \text{ops2}] \subseteq \text{vspace}`. .. warning:: This function is expensive to compute @@ -277,7 +278,7 @@ def check_commutation(ops1, ops2, vspace): def check_all_commuting(ops: List[Union[PauliSentence, np.ndarray, Operator]]): - r"""Helper function to check if all operators in a set of operators commute + r"""Helper function to check if all operators in ``ops`` commute. .. warning:: This function is expensive to compute @@ -326,7 +327,7 @@ def check_all_commuting(ops: List[Union[PauliSentence, np.ndarray, Operator]]): def check_cartan_decomp(k: List[PauliSentence], m: List[PauliSentence], verbose=True): - r"""Helper function to check the validity of a Cartan decomposition :math:`\mathfrak{g} = \mathfrak{k} \oplus \mathfrak{m}` + r"""Helper function to check the validity of a Cartan decomposition :math:`\mathfrak{g} = \mathfrak{k} \oplus \mathfrak{m}.` Check whether of not the following properties are fulfilled. @@ -586,7 +587,7 @@ def trace_inner_product( def change_basis_ad_rep(adj: np.ndarray, basis_change: np.ndarray): - r"""Apply the basis change between bases of operators to the adjoint representation. + r"""Apply a ``basis_change`` between bases of operators to the adjoint representation ``adj``. Assume the adjoint repesentation is given in terms of a basis :math:`\{b_j\}`, :math:`\text{ad_\mu}_{\alpha \beta} \propto \text{tr}\left(b_\mu \cdot [b_\alpha, b_\beta] \right)`. diff --git a/pennylane/labs/dla/involutions.py b/pennylane/labs/dla/involutions.py new file mode 100644 index 00000000000..7bc9f5ad752 --- /dev/null +++ b/pennylane/labs/dla/involutions.py @@ -0,0 +1,744 @@ +# 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. +"""Cartan involutions""" +from functools import singledispatch + +# pylint: disable=missing-function-docstring +from typing import Optional, Union + +import numpy as np + +import pennylane as qml +from pennylane import Y, Z, matrix +from pennylane.operation import Operator +from pennylane.pauli import PauliSentence + + +def khaneja_glaser_involution(op: Union[np.ndarray, PauliSentence, Operator], wire: int = None): + r"""Khaneja-Glaser involution, which is a type-:func:`~AIII` Cartan involution with p=q. + + Args: + op (Union[np.ndarray, PauliSentence, Operator]): Input operator + wire (int): Qubit wire on which to perform Khaneja-Glaser decomposition + + Returns: + bool: Whether ``op`` should go to the even or odd subspace of the decomposition + + .. seealso:: :func:`~cartan_decomp` + + See `Khaneja and Glaser (2000) `__ for reference. + + **Example** + + Let us first perform a single Khaneja-Glaser decomposition of + :math:`\mathfrak{g} = \mathfrak{su}(8)`, i.e. the Lie algebra of all Pauli words on 3 qubits. + + >>> g = list(qml.pauli.pauli_group(3)) # u(8) + >>> g = g[1:] # remove identity + >>> g = [op.pauli_rep for op in g] + + + We perform the first iteration on the first qubit. We use :func:`~cartan_decomp`. + + >>> from functools import partial + >>> from pennylane.labs.dla import cartan_decomp, khaneja_glaser_involution + >>> k0, m0 = cartan_decomp(g, partial(khaneja_glaser_involution, wire=0)) + >>> print(f"First iteration, AIII: {len(k0)}, {len(m0)}") + First iteration, AIII: 31, 32 + >>> assert qml.labs.dla.check_cartan_decomp(k0, m0) # check Cartan commutation relations + + We see that we split the 63-dimensional algebra :math:`\mathfrak{su}(8)` into a 31-dimensional + and a 32-dimensional subspace, and :func:`~check_cartan_decomp` verified that they satisfy + the commutation relations of a Cartan decomposition. + + To expand on this example, we continue process recursively on the :math:`\mathfrak{k}_0` + subalgebra. Before we can apply the Khaneja-Glaser (type AIII) decomposition a second time, + we need to a) remove the :math:`\mathfrak{u}(1)` center from + :math:`\mathfrak{k}_0=\mathfrak{su}(4)\oplus\mathfrak{su}(4)\oplus\mathfrak{u}(1)` to make it + semisimple and b) perform a different Cartan decomposition, which is sometimes termed class B. + + >>> center_k0 = qml.center(k0, pauli=True) # Compute center of k0 + >>> k0_semi = [op for op in k0 if op not in center_k0] # Remove center from k0 + >>> print(f"Removed operators {center_k0}, new length: {len(k0_semi)}") + Removed operators [1.0 * Z(0)], new length: 30 + + >>> from pennylane.labs.dla import ClassB + >>> k1, m1 = cartan_decomp(k0_semi, partial(ClassB, wire=0)) + >>> assert qml.labs.dla.check_cartan_decomp(k1, m1) + >>> print(f"First iteration, class B: {len(k1)}, {len(m1)}") + First iteration, class B: 15, 15 + + Now we arrived at the subalgebra :math:`\mathfrak{k}_1=\mathfrak{su}(4)` and can perform the + next iteration of the recursive decomposition. + + >>> k2, m2 = cartan_decomp(k1, partial(khaneja_glaser_involution, wire=1)) + >>> assert qml.labs.dla.check_cartan_decomp(k2, m2) + >>> print(f"Second iteration, AIII: {len(k2)}, {len(m2)}") + Second iteration, AIII: 7, 8 + + We can follow up on this by removing the center from :math:`\mathfrak{k}_1` again, and + applying a final class B decomposition: + + >>> center_k2 = qml.center(k2, pauli=True) + >>> k2_semi = [op for op in k2 if op not in center_k2] + >>> print(f"Removed operators {center_k2}, new length: {len(k2_semi)}") + Removed operators [1.0 * Z(1)], new length: 6 + + >>> k3, m3 = cartan_decomp(k2_semi, partial(ClassB, wire=1)) + >>> assert qml.labs.dla.check_cartan_decomp(k3, m3) + >>> print(f"Second iteration, class B: {len(k3)}, {len(m3)}") + Second iteration, class B: 3, 3 + + This concludes our example of the Khaneja-Glaser recursive decomposition, as we arrived at + easy to implement single-qubit operations. + """ + if wire is None: + raise ValueError( + "Please specify the wire for the Khaneja-Glaser involution via " + "functools.partial(khaneja_glaser_involution, wire=wire)" + ) + if isinstance(op, np.ndarray): + p = q = op.shape[0] // 2 + elif isinstance(op, (PauliSentence, Operator)): + p = q = 2 ** (len(op.wires) - 1) + else: + raise ValueError(f"Can't infer p and q from operator of type {type(op)}.") + return AIII(op, p, q, wire) + + +# Canonical involutions +# see https://arxiv.org/abs/2406.04418 appendix C + +# matrices + + +def int_log2(x): + """Compute the integer closest to log_2(x).""" + return int(np.round(np.log2(x))) + + +def is_qubit_case(p, q): + """Return whether p and q are the same and a power of 2.""" + return p == q and 2 ** int_log2(p) == p + + +def J(n, wire=None): + """This is the standard choice for the symplectic transformation operator. + For an :math:`N`-qubit system (:math:`n=2^N`), it equals :math:`Y_0`.""" + N = int_log2(n) + if 2**N == n: + if wire is None: + wire = 0 + return Y(wire).matrix(wire_order=range(N + 1)) + if wire is not None: + raise ValueError("The wire argument is only supported for n=2**N for some integer N.") + zeros = np.zeros((n, n)) + eye = np.eye(n) + return np.block([[zeros, -1j * eye], [1j * eye, zeros]]) + + +def Ipq(p, q, wire=None): + """This is the canonical transformation operator for AIII and BDI Cartan + decompositions. For an :math:`N`-qubit system (:math:`n=2^N`) and + :math:`p=q=n/2`, it equals :math:`Z_0`.""" + # If p = q and they are a power of two, use Pauli representation + if is_qubit_case(p, q): + if wire is None: + wire = 0 + return Z(wire).matrix(wire_order=range(int_log2(p) + 1)) + if wire is not None: + raise ValueError("The wire argument is only supported for p=q=2**N for some integer N.") + return np.block([[np.eye(p), np.zeros((p, q))], [np.zeros((q, p)), -np.eye(q)]]) + + +def Kpq(p, q, wire=None): + """This is the canonical transformation operator for CII Cartan + decompositions. For an :math:`N`-qubit system (:math:`n=2^N`) and + :math:`p=q=n/2`, it equals :math:`Z_1`.""" + # If p = q and they are a power of two, use Pauli representation + if is_qubit_case(p, q): + if wire is None: + wire = 1 + return Z(wire).matrix(wire_order=range(int_log2(p) + 1)) + if wire is not None: + raise ValueError("The wire argument is only supported for p=q=2**N for some integer N.") + zeros = np.zeros((p + q, p + q)) + d = np.diag(np.concatenate([np.ones(p), -np.ones(q)])) + KKm = np.block([[d, zeros], [zeros, d]]) + return KKm + + +def AI(op: Union[np.ndarray, PauliSentence, Operator]) -> bool: + r"""Canonical Cartan decomposition of type AI, given by :math:`\theta: x \mapsto x^\ast`. + + .. note:: Note that we work with Hermitian + operators internally, so that the input will be multiplied by :math:`i` before + evaluating the involution. + + Args: + op (Union[np.ndarray, PauliSentence, Operator]): Operator on which the involution is + evaluated and for which the parity under the involution is returned. + + Returns: + bool: Whether or not the input operator (times :math:`i`) is in the eigenspace of the + involution :math:`\theta` with eigenvalue :math:`+1`. + """ + return _AI(op) + + +@singledispatch +def _AI(op): # pylint:disable=unused-argument + r"""Default implementation of the canonical form of the AI involution + :math:`\theta: x \mapsto x^\ast`. + """ + raise NotImplementedError(f"Involution not implemented for operator {op} of type {type(op)}") + + +@_AI.register(np.ndarray) +def _AI_matrix(op: np.ndarray) -> bool: + r"""Matrix implementation of the canonical form of the AI involution + :math:`\theta: x \mapsto x^\ast`. + """ + op = op * 1j + return np.allclose(op, op.conj()) + + +@_AI.register(PauliSentence) +def _AI_ps(op: PauliSentence) -> bool: + r"""PauliSentence implementation of the canonical form of the AI involution + :math:`\theta: x \mapsto x^\ast`. + """ + parity = [] + for pw in op: + result = sum(el == "Y" for el in pw.values()) + parity.append(bool(result % 2)) + + # only makes sense if parity is the same for all terms + assert all(parity) or not any(parity) + return parity[0] + + +@_AI.register(Operator) +def _AI_op(op: Operator) -> bool: + r"""Operator implementation of the canonical form of the AI involution + :math:`\theta: x \mapsto x^\ast`. + """ + return _AI_ps(op.pauli_rep) + + +def AII(op: Union[np.ndarray, PauliSentence, Operator], wire: Optional[int] = None) -> bool: + r"""Canonical Cartan decomposition of type AII, given by :math:`\theta: x \mapsto Y_0 x^\ast Y_0`. + + .. note:: Note that we work with Hermitian + operators internally, so that the input will be multiplied by :math:`i` before + evaluating the involution. + + Args: + op (Union[np.ndarray, PauliSentence, Operator]): Operator on which the involution is + evaluated and for which the parity under the involution is returned. + wire (int): The wire on which the Pauli-:math:`Y` operator acts to implement the + involution. Will default to ``0`` if ``None``. + + Returns: + bool: Whether or not the input operator (times :math:`i`) is in the eigenspace of the + involution :math:`\theta` with eigenvalue :math:`+1`. + """ + return _AII(op, wire) + + +@singledispatch +def _AII(op, wire=None): # pylint:disable=unused-argument + r"""Default implementation of the canonical form of the AII involution + :math:`\theta: x \mapsto Y_0 x^\ast Y_0`. + """ + raise NotImplementedError(f"Involution not implemented for operator {op} of type {type(op)}") + + +@_AII.register(np.ndarray) +def _AII_matrix(op: np.ndarray, wire: Optional[int] = None) -> bool: + r"""Matrix implementation of the canonical form of the AII involution + :math:`\theta: x \mapsto Y_0 x^\ast Y_0`. + """ + op = op * 1j + + y = J(op.shape[-1] // 2, wire=wire) + return np.allclose(op, y @ op.conj() @ y) + + +@_AII.register(PauliSentence) +def _AII_ps(op: PauliSentence, wire: Optional[int] = None) -> bool: + r"""PauliSentence implementation of the canonical form of the AII involution + :math:`\theta: x \mapsto Y_0 x^\ast Y_0`. + """ + if wire is None: + wire = 0 + parity = [] + for pw in op: + result = sum(el == "Y" for el in pw.values()) + (pw.get(wire, "I") in "XZ") + parity.append(bool(result % 2)) + + # only makes sense if parity is the same for all terms + assert all(parity) or not any(parity) + return parity[0] + + +@_AII.register(Operator) +def _AII_op(op: Operator, wire: Optional[int] = None) -> bool: + r"""Operator implementation of the canonical form of the AII involution + :math:`\theta: x \mapsto Y_0 x^\ast Y_0`. + """ + return _AII_ps(op.pauli_rep, wire) + + +def AIII( + op: Union[np.ndarray, PauliSentence, Operator], + p: int = None, + q: int = None, + wire: Optional[int] = None, +) -> bool: + r"""Canonical Cartan decomposition of type AIII, given by :math:`\theta: x \mapsto I_{p,q} x I_{p,q}`. + + The matrix :math:`I_{p,q}` is given by + + .. math:: + I_{p,q}=\text{diag}(\underset{p \text{times}}{\underbrace{1, \dots 1}}, + \underset{q \text{times}}{\underbrace{-1, \dots -1}}). + + For :math:`p=q=2^N` for some integer :math:`N`, we have :math:`I_{p,q}=Z_0`. + + .. note:: Note that we work with Hermitian operators internally, so that the input will be + multiplied by :math:`i` before evaluating the involution. + + Args: + op (Union[np.ndarray, PauliSentence, Operator]): Operator on which the involution is + evaluated and for which the parity under the involution is returned. + p (int): Dimension of first subspace. + q (int): Dimension of second subspace. + wire (int): The wire on which the Pauli-:math:`Z` operator acts to implement the + involution. Will default to ``0`` if ``None``. + + Returns: + bool: Whether or not the input operator (times :math:`i`) is in the eigenspace of the + involution :math:`\theta` with eigenvalue :math:`+1`. + """ + if p is None or q is None: + raise ValueError( + "please specify p and q for the involution via functools.partial(AIII, p=p, q=q)" + ) + return _AIII(op, p, q, wire) + + +@singledispatch +def _AIII(op, p=None, q=None, wire=None): # pylint:disable=unused-argument + r"""Default implementation of the canonical form of the AIII involution + :math:`\theta: x \mapsto I_{p,q} x I_{p,q}`. + """ + raise NotImplementedError(f"Involution not implemented for operator {op} of type {type(op)}") + + +@_AIII.register(np.ndarray) +def _AIII_matrix(op: np.ndarray, p: int = None, q: int = None, wire: Optional[int] = None) -> bool: + r"""Matrix implementation of the canonical form of the AIII involution + :math:`\theta: x \mapsto I_{p,q} x I_{p,q}`. + """ + op = op * 1j + + z = Ipq(p, q, wire=wire) + return np.allclose(op, z @ op @ z) + + +@_AIII.register(PauliSentence) +def _AIII_ps(op: PauliSentence, p: int = None, q: int = None, wire: Optional[int] = None) -> bool: + r"""PauliSentence implementation of the canonical form of the AIII involution + :math:`\theta: x \mapsto I_{p,q} x I_{p,q}`. + """ + if is_qubit_case(p, q): + + if wire is None: + wire = 0 + parity = [pw.get(wire, "I") in "IZ" for pw in op] + + # only makes sense if parity is the same for all terms + assert all(parity) or not any(parity) + return parity[0] + + # If it is not a qubit case, use the matrix representation + return _AIII_matrix(matrix(op, wire_order=sorted(op.wires)), p, q, wire) + + +@_AIII.register(Operator) +def _AIII_op(op: Operator, p: int = None, q: int = None, wire: Optional[int] = None) -> bool: + r"""Operator implementation of the canonical form of the AIII involution + :math:`\theta: x \mapsto I_{p,q} x I_{p,q}`. + """ + return _AIII_ps(op.pauli_rep, p, q, wire) + + +def BDI( + op: Union[np.ndarray, PauliSentence, Operator], + p: int = None, + q: int = None, + wire: Optional[int] = None, +) -> bool: + r"""Canonical Cartan decomposition of type BDI, given by :math:`\theta: x \mapsto I_{p,q} x I_{p,q}`. + + The matrix :math:`I_{p,q}` is given by + + .. math:: + I_{p,q}=\text{diag}(\underset{p \text{times}}{\underbrace{1, \dots 1}}, + \underset{q \text{times}}{\underbrace{-1, \dots -1}}). + + For :math:`p=q=2^N` for some integer :math:`N`, we have :math:`I_{p,q}=Z_0`. + + .. note:: Note that we work with Hermitian operators internally, so that the input will be + multiplied by :math:`i` before evaluating the involution. + + + Args: + op (Union[np.ndarray, PauliSentence, Operator]): Operator on which the involution is + evaluated and for which the parity under the involution is returned. + p (int): Dimension of first subspace. + q (int): Dimension of second subspace. + wire (int): The wire on which the Pauli-:math:`Z` operator acts to implement the + involution. Will default to ``0`` if ``None``. + + Returns: + bool: Whether or not the input operator (times :math:`i`) is in the eigenspace of the + involution :math:`\theta` with eigenvalue :math:`+1`. + """ + return AIII(op, p, q, wire) + + +def CI(op: Union[np.ndarray, PauliSentence, Operator]) -> bool: + r"""Canonical Cartan decomposition of type CI, given by :math:`\theta: x \mapsto x^\ast`. + + .. note:: Note that we work with Hermitian + operators internally, so that the input will be multiplied by :math:`i` before + evaluating the involution. + + Args: + op (Union[np.ndarray, PauliSentence, Operator]): Operator on which the involution is + evaluated and for which the parity under the involution is returned. + + Returns: + bool: Whether or not the input operator (times :math:`i`) is in the eigenspace of the + involution :math:`\theta` with eigenvalue :math:`+1`. + """ + return AI(op) + + +def CII( + op: Union[np.ndarray, PauliSentence, Operator], + p: int = None, + q: int = None, + wire: Optional[int] = None, +) -> bool: + r"""Canonical Cartan decomposition of type CII, given by :math:`\theta: x \mapsto K_{p,q} x K_{p,q}`. + + The matrix :math:`K_{p,q}` is given by + + .. math:: + + K_{p,q}=\text{diag}( + \underset{p \text{times}}{\underbrace{1, \dots 1}}, + \underset{q \text{times}}{\underbrace{-1, \dots -1}}, + \underset{p \text{times}}{\underbrace{1, \dots 1}}, + \underset{q \text{times}}{\underbrace{-1, \dots -1}}, + ). + + For :math:`p=q=2^N` for some integer :math:`N`, we have :math:`K_{p,q}=Z_1`. + + .. note:: Note that we work with Hermitian operators internally, so that the input will be + multiplied by :math:`i` before evaluating the involution. + + Args: + op (Union[np.ndarray, PauliSentence, Operator]): Operator on which the involution is + evaluated and for which the parity under the involution is returned. + p (int): Dimension of first subspace. + q (int): Dimension of second subspace. + wire (int): The wire on which the Pauli-:math:`Z` operator acts to implement the + involution. Will default to ``1`` if ``None``. + + Returns: + bool: Whether or not the input operator (times :math:`i`) is in the eigenspace of the + involution :math:`\theta` with eigenvalue :math:`+1`. + """ + if p is None or q is None: + raise ValueError( + "please specify p and q for the involution via functools.partial(CII, p=p, q=q)" + ) + return _CII(op, p, q, wire) + + +@singledispatch +def _CII(op, p=None, q=None, wire=None): # pylint:disable=unused-argument + r"""Default implementation of the canonical form of the CII involution + :math:`\theta: x \mapsto K_{p,q} x K_{p,q}`. + """ + raise NotImplementedError(f"Involution not implemented for operator {op} of type {type(op)}") + + +@_CII.register(np.ndarray) +def _CII_matrix(op: np.ndarray, p: int = None, q: int = None, wire: Optional[int] = None) -> bool: + r"""Matrix implementation of the canonical form of the CII involution + :math:`\theta: x \mapsto K_{p,q} x K_{p,q}`. + """ + op = op * 1j + + z = Kpq(p, q, wire=wire) + return np.allclose(op, z @ op @ z) + + +@_CII.register(PauliSentence) +def _CII_ps(op: PauliSentence, p: int = None, q: int = None, wire: Optional[int] = None) -> bool: + r"""PauliSentence implementation of the canonical form of the CII involution + :math:`\theta: x \mapsto K_{p,q} x K_{p,q}`. + """ + if is_qubit_case(p, q): + + if wire is None: + wire = 1 + parity = [pw.get(wire, "I") in "IZ" for pw in op] + + # only makes sense if parity is the same for all terms + assert all(parity) or not any(parity) + return parity[0] + + # If it is not a qubit case, use the matrix representation + return _CII(matrix(op, wire_order=sorted(op.wires)), p, q, wire) + + +@_CII.register(Operator) +def _CII_op(op: Operator, p: int = None, q: int = None, wire: Optional[int] = None) -> bool: + r"""Operator implementation of the canonical form of the CII involution + :math:`\theta: x \mapsto K_{p,q} x K_{p,q}`. + """ + return _CII_ps(op.pauli_rep, p, q, wire) + + +def DIII(op: Union[np.ndarray, PauliSentence, Operator], wire: Optional[int] = None) -> bool: + r"""Canonical Cartan decomposition of type DIII, given by :math:`\theta: x \mapsto Y_0 x Y_0`. + + Args: + op (Union[np.ndarray, PauliSentence, Operator]): Operator on which the involution is + evaluated and for which the parity under the involution is returned. + wire (int): The wire on which the Pauli-:math:`Y` operator acts to implement the + involution. Will default to ``0`` if ``None``. + + Returns: + bool: Whether or not the input operator (times :math:`i`) is in the eigenspace of the + involution :math:`\theta` with eigenvalue :math:`+1`. + """ + return _DIII(op, wire) + + +@singledispatch +def _DIII(op, wire=None): # pylint:disable=unused-argument + r"""Default implementation of the canonical form of the DIII involution + :math:`\theta: x \mapsto Y_0 x Y_0`. + """ + raise NotImplementedError(f"Involution not implemented for operator {op} of type {type(op)}") + + +@_DIII.register(np.ndarray) +def _DIII_matrix(op: np.ndarray, wire: Optional[int] = None) -> bool: + r"""Matrix implementation of the canonical form of the DIII involution + :math:`\theta: x \mapsto Y_0 x Y_0`. + """ + y = J(op.shape[-1] // 2, wire=wire) + return np.allclose(op, y @ op @ y) + + +@_DIII.register(PauliSentence) +def _DIII_ps(op: PauliSentence, wire: Optional[int] = None) -> bool: + r"""PauliSentence implementation of the canonical form of the DIII involution + :math:`\theta: x \mapsto Y_0 x Y_0`. + """ + if wire is None: + wire = 0 + parity = [pw.get(wire, "I") in "IY" for pw in op] + + # only makes sense if parity is the same for all terms + assert all(parity) or not any(parity) + return parity[0] + + +@_DIII.register(Operator) +def _DIII_op(op: Operator, wire: Optional[int] = None) -> bool: + r"""Operator implementation of the canonical form of the DIII involution + :math:`\theta: x \mapsto Y_0 x Y_0`. + """ + return _DIII_ps(op.pauli_rep, wire) + + +def ClassB(op: Union[np.ndarray, PauliSentence, Operator], wire: Optional[int] = None) -> bool: + r"""Canonical Cartan decomposition of class B, given by :math:`\theta: x \mapsto Y_0 x Y_0`. + + Args: + op (Union[np.ndarray, PauliSentence, Operator]): Operator on which the involution is + evaluated and for which the parity under the involution is returned. + wire (int): The wire on which the Pauli-:math:`Y` operator acts to implement the + involution. Will default to ``0`` if ``None``. + + Returns: + bool: Whether or not the input operator (times :math:`i`) is in the eigenspace of the + involution :math:`\theta` with eigenvalue :math:`+1`. + """ + return DIII(op, wire) + + +# dispatch to different input types +def even_odd_involution(op: Union[PauliSentence, np.ndarray, Operator]) -> bool: + r"""The Even-Odd involution. + + This is defined in `quant-ph/0701193 `__. + For Pauli words and sentences, it comes down to counting non-trivial Paulis in Pauli words. + + Args: + op ( Union[PauliSentence, np.ndarray, Operator]): Input operator + + Returns: + bool: Boolean output ``True`` or ``False`` for odd (:math:`\mathfrak{k}`) and even parity subspace (:math:`\mathfrak{m}`), respectively + + .. seealso:: :func:`~cartan_decomp` + + **Example** + + >>> from pennylane import X, Y, Z + >>> from pennylane.labs.dla import even_odd_involution + >>> ops = [X(0), X(0) @ Y(1), X(0) @ Y(1) @ Z(2)] + >>> [even_odd_involution(op) for op in ops] + [True, False, True] + + Operators with an odd number of non-identity Paulis yield ``1``, whereas even ones yield ``0``. + + The function also works with dense matrix representations. + + >>> ops_m = [qml.matrix(op, wire_order=range(3)) for op in ops] + >>> [even_odd_involution(op_m) for op_m in ops_m] + [True, False, True] + + """ + return _even_odd_involution(op) + + +@singledispatch +def _even_odd_involution(op): # pylint:disable=unused-argument, missing-function-docstring + return NotImplementedError(f"Involution not defined for operator {op} of type {type(op)}") + + +@_even_odd_involution.register(PauliSentence) +def _even_odd_involution_ps(op: PauliSentence): + # Generalization to sums of Paulis: check each term and assert they all have the same parity + parity = [] + for pw in op.keys(): + parity.append(len(pw) % 2) + + # only makes sense if parity is the same for all terms, e.g. Heisenberg model + assert all( + parity[0] == p for p in parity + ), f"The Even-Odd involution is not well-defined for operator {op} as individual terms have different parity" + return bool(parity[0]) + + +@_even_odd_involution.register(np.ndarray) +def _even_odd_involution_matrix(op: np.ndarray): + """see Table CI in https://arxiv.org/abs/2406.04418""" + n = int(np.round(np.log2(op.shape[-1]))) + YYY = qml.prod(*[Y(i) for i in range(n)]) + YYY = qml.matrix(YYY, range(n)) + + transformed = YYY @ op.conj() @ YYY + if np.allclose(transformed, op): + return False + if np.allclose(transformed, -op): + return True + raise ValueError(f"The Even-Odd involution is not well-defined for operator {op}.") + + +@_even_odd_involution.register(Operator) +def _even_odd_involution_op(op: Operator): + """use pauli representation""" + return _even_odd_involution_ps(op.pauli_rep) + + +# dispatch to different input types +def concurrence_involution(op: Union[PauliSentence, np.ndarray, Operator]) -> bool: + r"""The Concurrence Canonical Decomposition :math:`\Theta(g) = -g^T` as a Cartan involution function. + + This is defined in `quant-ph/0701193 `__. + For Pauli words and sentences, it comes down to counting Pauli-Y operators. + + Args: + op ( Union[PauliSentence, np.ndarray, Operator]): Input operator + + Returns: + bool: Boolean output ``True`` or ``False`` for odd (:math:`\mathfrak{k}`) and even parity subspace (:math:`\mathfrak{m}`), respectively + + .. seealso:: :func:`~cartan_decomp` + + **Example** + + >>> from pennylane import X, Y, Z + >>> from pennylane.labs.dla import concurrence_involution + >>> ops = [X(0), X(0) @ Y(1), X(0) @ Y(1) @ Z(2), Y(0) @ Y(2)] + >>> [concurrence_involution(op) for op in ops] + [False, True, True, False] + + Operators with an odd number of ``Y`` operators yield ``1``, whereas even ones yield ``0``. + + The function also works with dense matrix representations. + + >>> ops_m = [qml.matrix(op, wire_order=range(3)) for op in ops] + >>> [even_odd_involution(op_m) for op_m in ops_m] + [False, True, True, False] + + """ + return _concurrence_involution(op) + + +@singledispatch +def _concurrence_involution(op): + return NotImplementedError(f"Involution not defined for operator {op} of type {type(op)}") + + +@_concurrence_involution.register(PauliSentence) +def _concurrence_involution_pauli(op: PauliSentence): + # Generalization to sums of Paulis: check each term and assert they all have the same parity + parity = [] + for pw in op.keys(): + result = sum(1 if el == "Y" else 0 for el in pw.values()) + parity.append(result % 2) + + # only makes sense if parity is the same for all terms, e.g. Heisenberg model + assert all( + parity[0] == p for p in parity + ), f"The concurrence canonical decomposition is not well-defined for operator {op} as individual terms have different parity" + return bool(parity[0]) + + +@_concurrence_involution.register(Operator) +def _concurrence_involution_operator(op: Operator): + return _concurrence_involution_matrix(op.matrix()) + + +@_concurrence_involution.register(np.ndarray) +def _concurrence_involution_matrix(op: np.ndarray): + if np.allclose(op, -op.T): + return True + if np.allclose(op, op.T): + return False + raise ValueError( + f"The concurrence canonical decomposition is not well-defined for operator {op}" + ) diff --git a/pennylane/labs/tests/dla/test_cartan_subalgebra.py b/pennylane/labs/tests/dla/test_cartan_subalgebra.py new file mode 100644 index 00000000000..22689758f6c --- /dev/null +++ b/pennylane/labs/tests/dla/test_cartan_subalgebra.py @@ -0,0 +1,52 @@ +# 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/dla/lie_closure_dense.py functionality""" +import numpy as np + +# pylint: disable=no-self-use,too-few-public-methods,missing-class-docstring +import pytest + +import pennylane as qml +from pennylane import X, Z +from pennylane.labs.dla import ( + cartan_decomp, + cartan_subalgebra, + check_cartan_decomp, + even_odd_involution, +) + + +@pytest.mark.parametrize("n, len_g, len_h, len_mtilde", [(2, 6, 2, 2), (3, 15, 2, 6)]) +def test_Ising(n, len_g, len_h, len_mtilde): + """Test Cartan subalgebra of 2 qubit Ising model""" + gens = [X(w) @ X(w + 1) for w in range(n - 1)] + [Z(w) for w in range(n)] + gens = [op.pauli_rep for op in gens] + g = qml.lie_closure(gens, pauli=True) + + k, m = cartan_decomp(g, even_odd_involution) + assert check_cartan_decomp(k, m) + + g = k + m + assert len(g) == len_g + + adj = qml.structure_constants(g) + + newg, k, mtilde, h, new_adj = cartan_subalgebra(g, k, m, adj, start_idx=0) + assert len(h) == len_h + assert len(mtilde) == len_mtilde + assert len(h) + len(mtilde) == len(m) + + new_adj_re = qml.structure_constants(newg) + + assert np.allclose(new_adj_re, new_adj) diff --git a/pennylane/labs/tests/dla/test_involutions.py b/pennylane/labs/tests/dla/test_involutions.py new file mode 100644 index 00000000000..c828ecfd089 --- /dev/null +++ b/pennylane/labs/tests/dla/test_involutions.py @@ -0,0 +1,229 @@ +# 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/lie_closure_dense.py functionality""" +from functools import partial + +import numpy as np + +# pylint: disable=too-few-public-methods, protected-access, no-self-use +import pytest + +import pennylane as qml +from pennylane import X, Y, Z +from pennylane.labs.dla import AI, AII, AIII, BDI, CI, CII, DIII, ClassB, khaneja_glaser_involution +from pennylane.labs.dla.involutions import Ipq, J, Kpq + + +class TestMatrixConstructors: + """Tests for the matrix constructing methods used in Cartan involutions.""" + + @pytest.mark.parametrize( + "num_wires, wire", + [(N, w) for N in range(5) for w in range(1, N + 1)] + [(N, None) for N in range(5)], + ) + def test_J_qubit_cases(self, num_wires, wire): + """Test J matrix constructor for qubit cases.""" + out = J(2**num_wires, wire) + op = qml.Y(0 if wire is None else wire) + expected = op.matrix(wire_order=range(num_wires + 1)) + assert np.allclose(out, expected) + + @pytest.mark.parametrize("n", [3, 5, 6, 7]) + def test_J_non_qubit_cases(self, n): + """Test J matrix constructor for non-qubit cases.""" + out = J(n) + expected = np.zeros((2 * n, 2 * n), dtype=complex) + for i in range(n): + expected[i, n + i] = -1j + expected[n + i, i] = 1j + assert np.allclose(out, expected) + + @pytest.mark.parametrize( + "num_wires, wire", + [(N, w) for N in range(5) for w in range(1, N + 1)] + [(N, None) for N in range(5)], + ) + def test_I_qubit_cases(self, num_wires, wire): + """Test I_pq matrix constructor for qubit cases.""" + p = 2**num_wires + out = Ipq(p, p, wire) + op = qml.Z(0 if wire is None else wire) + expected = op.matrix(wire_order=range(num_wires + 1)) + assert np.allclose(out, expected) + + @pytest.mark.parametrize("p, q", [(1, 3), (2, 5), (3, 2), (6, 6)]) + def test_Ipq_non_qubit_cases(self, p, q): + """Test I_pq matrix constructor for non-qubit cases.""" + out = Ipq(p, q) + expected = np.diag(np.concatenate([np.ones(p), -np.ones(q)], axis=0)) + assert np.allclose(out, expected) + + @pytest.mark.parametrize( + "num_wires, wire", + [(N, w) for N in range(5) for w in range(1, N + 1)] + [(N, None) for N in range(1, 5)], + ) + def test_Kpq_qubit_cases(self, num_wires, wire): + """Test K_pq matrix constructor for qubit cases.""" + p = 2**num_wires + out = Kpq(p, p, wire) + op = qml.Z(1 if wire is None else wire) + expected = op.matrix(wire_order=range(num_wires + 1)) + assert np.allclose(out, expected) + + @pytest.mark.parametrize("p, q", [(1, 3), (2, 5), (3, 2), (6, 6)]) + def test_Kpq_non_qubit_cases(self, p, q): + """Test K_pq matrix constructor for non-qubit cases.""" + out = Kpq(p, q) + expected = np.diag( + np.concatenate([np.ones(p), -np.ones(q), np.ones(p), -np.ones(q)], axis=0) + ) + assert np.allclose(out, expected) + + +class TestInvolutionExceptions: + """Test exceptions being raised by involutions.""" + + def test_non_qubit_wire_given(self): + """Test that an error is raised if the wire is specified but it is a non-qubit case.""" + with pytest.raises(ValueError, match="wire argument is only supported"): + J(3, wire=0) + + with pytest.raises(ValueError, match="wire argument is only supported"): + Ipq(3, 5, wire=0) + + with pytest.raises(ValueError, match="wire argument is only supported"): + Kpq(3, 5, wire=0) + + def test_p_or_q_missing(self): + """Test that an error is raised if p or q (or both) are not given for the involutions + AIII, BDI, CII.""" + for p, q in [(None, 2), (5, None), (None, None)]: + for invol in [AIII, BDI, CII]: + with pytest.raises(ValueError, match="please specify p and q"): + invol(X(0) @ Y(1), p=p, q=q, wire=0) + + def test_Khaneja_Glaser_exceptions(self): + """Test that the Khaneja-Glaser involution raises custom exceptions related + to wire and infering p and q.""" + op = qml.X(0) @ qml.Y(1) + with pytest.raises(ValueError, match="Please specify the wire for the Khaneja"): + khaneja_glaser_involution(op) + + [op] = op.pauli_rep + with pytest.raises(ValueError, match="Can't infer p and q from operator of type -1 + (X(0) @ Y(1) @ Z(2) - Z(0) @ X(1) @ Y(2), False), # (True, False) -> -1 + (X(1) @ Y(0) @ Z(2) - Y(0) @ Y(1) @ Y(2), True), # (True, True) -> 1 + (Y(0) @ Z(2) @ Y(1), False), # (False, True) -> -1 + (X(1) @ Z(2) @ X(0) + X(0) @ Y(1) @ Y(2), True), # (False, False) -> 1 +] + +AIII_cases = [ # I/Z on wire 0? + (X(0) @ Y(1) @ Z(2), False), + (X(0) @ Y(1) @ Z(2) - Y(0) @ X(1) @ Y(2), False), + (Z(0) @ X(1) @ Z(2) - Z(0) @ Y(1) @ Y(2), True), + (Y(0) @ Y(1) @ Z(2), False), + (Y(1) @ Z(2), True), + (Z(0) + 0.2 * Y(1) @ Y(2), True), +] + +BDI_cases = AIII_cases # BDI = AIII + +CI_cases = AI_cases # CI = AI + +CII_cases = [ # I/Z on wire 1? + (X(0) @ Z(1) @ Z(2), True), + (X(0) @ Z(2) - Z(0) @ Z(1) @ Y(2), True), + (Z(0) @ X(1) @ Z(2) - Y(0) @ Y(1) @ Y(2), False), + (Y(0) @ X(1) @ Z(2), False), + (Y(1) @ Z(2), False), + (Z(0) + 0.2 * Y(0) @ Y(2), True), +] + +DIII_cases = [ # I/Y on wire 0? + (X(0) @ Y(1) @ Z(2), False), + (X(0) @ Y(1) @ Z(2) - Z(0) @ X(1) @ Y(2), False), + (X(1) @ Y(0) @ Z(2) - Y(0) @ Y(1) @ Y(2), True), + (Y(0) @ Y(1) @ Z(2), True), + (X(1) @ Z(2), True), + (Z(0) + 0.2 * X(2) @ X(0), False), +] + +ClassB_cases = DIII_cases # ClassB = DIII + + +class TestInvolutions: + """Test the involutions themselves.""" + + def run_test_case(self, op, expected, invol): + """Run a generic test case for a given operator and involution""" + inputs = [op, op.pauli_rep, qml.matrix(op, wire_order=[0, 1, 2])] + outputs = [invol(_input) for _input in inputs] + if expected: + assert all(outputs) + else: + assert not any(outputs) + + @pytest.mark.parametrize("op, expected", AI_cases) + def test_AI(self, op, expected): + """Test singledispatch for AI involution""" + self.run_test_case(op, expected, AI) + + @pytest.mark.parametrize("op, expected", AII_cases) + def test_AII(self, op, expected): + """Test singledispatch for AII involution""" + self.run_test_case(op, expected, AII) + + @pytest.mark.parametrize("op, expected", AIII_cases) + def test_AIII(self, op, expected): + """Test singledispatch for AIII involution""" + self.run_test_case(op, expected, partial(AIII, p=4, q=4)) + # Khaneja-Glaser is just AIII with automatically inferred p and q. + self.run_test_case(op, expected, partial(khaneja_glaser_involution, wire=0)) + + @pytest.mark.parametrize("op, expected", BDI_cases) + def test_BDI(self, op, expected): + """Test singledispatch for BDI involution""" + self.run_test_case(op, expected, partial(BDI, p=4, q=4)) + + @pytest.mark.parametrize("op, expected", CI_cases) + def test_CI(self, op, expected): + """Test singledispatch for CI involution""" + self.run_test_case(op, expected, CI) + + @pytest.mark.parametrize("op, expected", CII_cases) + def test_CII(self, op, expected): + """Test singledispatch for CII involution""" + self.run_test_case(op, expected, partial(CII, p=4, q=4)) + + @pytest.mark.parametrize("op, expected", DIII_cases) + def test_DIII(self, op, expected): + """Test singledispatch for DIII involution""" + self.run_test_case(op, expected, DIII) + + @pytest.mark.parametrize("op, expected", ClassB_cases) + def test_ClassB(self, op, expected): + """Test singledispatch for ClassB involution""" + self.run_test_case(op, expected, ClassB) diff --git a/pennylane/pauli/dla/center.py b/pennylane/pauli/dla/center.py index d810b429b3e..938cbe145b5 100644 --- a/pennylane/pauli/dla/center.py +++ b/pennylane/pauli/dla/center.py @@ -23,7 +23,7 @@ from pennylane.pauli.dla import structure_constants -def _intersect_bases(basis_0, basis_1): +def _intersect_bases(basis_0, basis_1, rcond=None): 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 @@ -34,7 +34,7 @@ def _intersect_bases(basis_0, basis_1): 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])) + augmented_basis = null_space(np.hstack([basis_0, -basis_1]), rcond=rcond) # 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