-
Notifications
You must be signed in to change notification settings - Fork 595
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
1 changed file
with
291 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,291 @@ | ||
# 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 | ||
import numpy as np | ||
from scipy.linalg import null_space | ||
|
||
import pennylane as qml | ||
from pennylane.operation import Operator | ||
from pennylane.pauli import PauliSentence, PauliVSpace | ||
from pennylane.pauli.dla.center import _intersect_bases | ||
from pennylane.typing import TensorLike | ||
|
||
|
||
def _gram_schmidt(X): | ||
"""Orthogonalize basis of column vectors in X""" | ||
_, R = np.linalg.qr(X.T, mode="complete") | ||
return R.T | ||
|
||
|
||
def _is_independent(v, A, tol=1e-14): | ||
v /= np.linalg.norm(v) | ||
v = v - A @ qml.math.linalg.solve(qml.math.conj(A.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, _ = np.linalg.qr(h.T) | ||
|
||
# Step 2: Project each vector in m onto the orthogonal complement of span(h) | ||
projections = m - np.dot(np.dot(m, Q), Q.T) | ||
|
||
# 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] | ||
|
||
return basis.T # Transpose to get row vectors | ||
|
||
|
||
def cartan_subalgebra(g, k, m, ad, which=0, tol=1e-14, verbose=0): | ||
r""" | ||
Compute Cartan subalgebra of g given m, the odd parity subspace. I.e. the maximal Abelian subalgebra in m. | ||
.. seealso:: :func:`~cartan_decomposition`, :func:`~structure_constants` | ||
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{m}` 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}` | ||
Returns: | ||
np_newg, np_k, np_mtilde, np_h, new_adj, basis_change | ||
np.ndarray: Adjoint vectors for new ordered :math:`\mathfrak{g}` of dimension :math:`|\mathfrak{g}| \times |\mathfrak{g}|`, ordered as :math:`\mathfrak{g} = \mathfrak{k} \oplus \tilde{\mathfrak{m}} \oplus \mathfrak{h}` | ||
np.ndarray: Adjoint vectors for :math:`\mathfrak{k}` of dimension :math:`|\mathfrak{k}| \times |\mathfrak{g}|`. | ||
np.ndarray: Adjoint vectors for :math:`\tilde{\mathfrak{m}}` of dimension :math:`|\tilde{\mathfrak{m}}| \times |\mathfrak{g}|`. | ||
np.ndarray: Adjoint vectors for :math:`\mathfrak{h}` of dimension :math:`|\mathfrak{h}| \times |\mathfrak{g}|`. | ||
np.ndarray: The new :math:`|\mathfrak{g}| \times |\mathfrak{g}| \times |\mathfrak{g}|` dimensional adjoint representation of :math:`\mathfrak{g}_\text{new}`, according to the ordering of new vectors. | ||
**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_decomposition, cartan_subalgebra, even_odd_involution | ||
>>> g = list(qml.pauli.pauli_group(2)) # u(4) | ||
>>> g = g[1:] # remove identity -> su(4) | ||
>>> k, m = cartan_decomposition(g, even_odd_involution) | ||
>>> g = k + m # re-order g to separate k and m | ||
>>> adj = qml.structure_constants(g) | ||
>>> np_newg, np_k, np_mtilde, np_h, new_adj = cartan_subalgebra(g, k, m, adj) | ||
All results are 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``. | ||
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 = qml.simplify(op) | ||
>>> op | ||
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 | ||
[Z(0) @ Z(1), (-1+0j) * (Y(0) @ X(1)), (-1+0j) * (X(0) @ Y(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_decomposition` and :func:`~structure_constants`. | ||
As an example, we take the Lie algebra of the Heisenberg model with generators :math:`\{X_i X_j, Y_i Y_j, Z_i, Z_j\}`. | ||
>>> from pennylane.labs.dla import lie_closure_dense, cartan_decomposition | ||
>>> 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, can perform the Cartan decomposition. We take the :math:`~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}] \subeq \mathfrak{k}`, :math:`[\mathfrak{k}, \mathfrak{m}] \subeq \mathfrak{m}` and :math:`[\mathfrak{m}, \mathfrak{m}] \subeq \mathfrak{k}`, | ||
which we can check using the helper function :math:`~check_cartan_decomp`. | ||
>>> from pennylane.labs.dla import even_odd_involution, check_cartan_decomp | ||
>>> k, m = cartan_decomposition(g, even_odd_involution) | ||
>>> assert check_cartan_decomp(k, m) # check commutation relations to be valid Cartan decomposition | ||
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 :math:`\mathfrak{k}` first, and then | ||
all remaining operators from :math:`\mathfrak{m}`. | ||
>>> 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}`. | ||
>>> np_newg, np_k, np_mtilde, np_h, new_adj, basis_change = cartan_subalgebra(g, k, m, adj, which=3) | ||
The final results are all returned as what we call adjoint representation vectors. | ||
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 | ||
[0.35355339059327384 * Y(1) @ Y(2), | ||
0.3535533905932738 * Y(0) @ Y(2), | ||
-0.3535533905932738 * Y(0) @ Y(1)] | ||
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``. | ||
In particular, we can compute commutators between | ||
""" | ||
|
||
np_m = op_to_adjvec(m, g) | ||
np_h = op_to_adjvec([m[which]], g) | ||
|
||
iteration = 1 | ||
while True: | ||
if verbose: | ||
print(f"iteration: {iteration}") | ||
kernel_intersection = np_m.T | ||
for h_i in np_h: | ||
|
||
# obtain adjoint rep of candidate h_i | ||
adjoint_of_h_i = np.einsum("gab,a->gb", ad, h_i) | ||
# 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, new_kernel, rcond=tol) | ||
|
||
if kernel_intersection.shape[1] == len(np_h): | ||
# No new vector was added from all the kernels | ||
break | ||
|
||
kernel_intersection = _gram_schmidt(kernel_intersection) # orthogonalize | ||
for vec in kernel_intersection.T: | ||
if _is_independent(vec, np.array(np_h).T, tol): | ||
np_h = np.vstack([np_h, vec]) | ||
break | ||
|
||
iteration += 1 | ||
|
||
np_h = _gram_schmidt(np_h.T).T # orthogonalize Abelian subalgebra | ||
np_k = op_to_adjvec(k, g) # adjoint vectors of k space for re-ordering | ||
|
||
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 | ||
# implementation to be tested: | ||
basis_change = np.linalg.lstsq(np.vstack([np_k, np_m]).T, np_newg.T, rcond=None)[ | ||
0 | ||
] # basis change old g @ X = new g => adj_new = contact(adj, X, X, X) | ||
new_adj = np.einsum("kmn,ki,mj,nl->ijl", ad, basis_change, basis_change, basis_change) | ||
|
||
return np_newg, np_k, np_mtilde, np_h, new_adj | ||
|
||
|
||
def adjvec_to_op(adj_vecs, basis): | ||
"""Transform adjoint vectors representing operators back into a basis of choice | ||
Args: | ||
adj_vecs (np.ndarray): collection of vectors with shape ``(batch, dim_basis)`` | ||
basis (List[Union[PauliSentence, Operator, np.ndarray]]): collection of basis operators | ||
""" | ||
res = [] | ||
# currently agnostic to dense or op input, but could be vectorized in the dense case (TODO?) | ||
if all(isinstance(op, PauliSentence) for op in basis): | ||
for vec in adj_vecs: | ||
op_j = sum(c * op for c, op in zip(vec, basis)) | ||
op_j.simplify() | ||
res.append(op_j) | ||
return res | ||
|
||
if all(isinstance(op, Operator) for op in basis): | ||
for vec in adj_vecs: | ||
op_j = sum(c * op for c, op in zip(vec, basis)) | ||
op_j = qml.simplify(op_j) | ||
res.append(op_j) | ||
return res | ||
|
||
if all(isinstance(op, np.ndarray) for op in basis): | ||
res = np.einsum("id,dmn->imn", adj_vecs, basis) # i, dimg | dimg m n | ||
return res | ||
|
||
return NotImplementedError | ||
|
||
|
||
def op_to_adjvec(ops, basis): | ||
"""Project a batch of ops onto a given basis""" | ||
if isinstance(basis, PauliVSpace): | ||
basis = basis.basis | ||
|
||
if all(isinstance(op, Operator) for op in ops) and all( | ||
isinstance(op, Operator) for op in basis | ||
): | ||
ops = [op.pauli_rep for op in ops] | ||
basis = [op.pauli_rep for op in basis] | ||
|
||
# PauliSentence branch | ||
if all(isinstance(op, PauliSentence) for op in ops) and all( | ||
isinstance(op, PauliSentence) for op in basis | ||
): | ||
res = [] | ||
for op in ops: | ||
rep = np.zeros((len(basis),), dtype=complex) | ||
for i, basis_i in enumerate(basis): | ||
# v = ∑ (v · e_j / ||e_j||^2) * e_j | ||
value = (basis_i @ op).trace() | ||
value = value / (basis_i @ basis_i).trace() | ||
|
||
rep[i] = value | ||
|
||
res.append(rep) | ||
|
||
return np.array(res) | ||
|
||
# dense branch | ||
if all(isinstance(op, TensorLike) for op in ops) and all( | ||
isinstance(op, TensorLike) for op in basis | ||
): | ||
basis = np.array(basis) | ||
# if len(ops.shape) == 2: | ||
# ops = np.array([ops]) | ||
ops = np.array(ops) | ||
|
||
res = np.einsum("bij,cji->bc", ops, basis) | ||
|
||
return res | ||
|
||
return NotImplemented |