Source code for pennylane.math.decomposition

# Copyright 2025 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.

"""Utility functions for decompositions are available from ``qml.math.decomposition``."""

import functools

import numpy as np

from pennylane import math

try:
    import jax
    import jax.numpy as jnp
except ModuleNotFoundError:  # pragma: no cover
    ...


[docs] def zyz_rotation_angles(U, return_global_phase=False): r"""Compute the rotation angles :math:`\phi`, :math:`\theta`, and :math:`\omega` and the phase :math:`\alpha` of a 2x2 unitary matrix as a product of Z and Y rotations in the form :math:`e^{i\alpha} RZ(\omega) RY(\theta) RZ(\phi)` Args: U (array): 2x2 unitary matrix return_global_phase (bool): if True, returns the global phase as well. Returns: tuple: The rotation angles :math:`\phi`, :math:`\theta`, and :math:`\omega` and the global phase :math:`\alpha` if ``return_global_phase=True``. """ U, alpha = math.convert_to_su2(U, return_global_phase=True) # assume U = [[a, b], [c, d]], then here we take U[0, 1] as b abs_b = math.clip(math.abs(U[..., 0, 1]), 0, 1) theta = 2 * math.arcsin(abs_b) EPS = math.finfo(U.dtype).eps half_phi_plus_omega = math.angle(U[..., 1, 1] + EPS) half_omega_minus_phi = math.angle(U[..., 1, 0] + EPS) phi = half_phi_plus_omega - half_omega_minus_phi omega = half_phi_plus_omega + half_omega_minus_phi # Normalize the angles phi = math.squeeze(phi % (4 * np.pi)) theta = math.squeeze(theta % (4 * np.pi)) omega = math.squeeze(omega % (4 * np.pi)) return (phi, theta, omega, alpha) if return_global_phase else (phi, theta, omega)
[docs] def xyx_rotation_angles(U, return_global_phase=False): r"""Compute the rotation angles :math:`\lambda`, :math:`\theta`, and :math:`\phi` and the phase :math:`\alpha` of a 2x2 unitary matrix as a product of X and Y rotations in the form :math:`e^{i\alpha} RX(\phi) RY(\theta) RX(\lambda)`. Args: U (array): 2x2 unitary matrix return_global_phase (bool): if True, returns the global phase as well. Returns: tuple: The rotation angles :math:`\lambda`, :math:`\theta`, and :math:`\phi` and the global phase :math:`\alpha` if ``return_global_phase=True``. """ U, alpha = math.convert_to_su2(U, return_global_phase=True) EPS = math.finfo(U.dtype).eps half_lam_plus_phi = math.arctan2(-math.imag(U[..., 0, 1]), math.real(U[..., 0, 0]) + EPS) half_lam_minus_phi = math.arctan2(math.imag(U[..., 0, 0]), -math.real(U[..., 0, 1]) + EPS) lam = half_lam_plus_phi + half_lam_minus_phi phi = half_lam_plus_phi - half_lam_minus_phi theta = math.where( math.isclose(math.sin(half_lam_plus_phi), math.zeros_like(half_lam_plus_phi)), 2 * math.arccos(math.clip(math.real(U[..., 1, 1]) / math.cos(half_lam_plus_phi), -1, 1)), 2 * math.arccos(math.clip(-math.imag(U[..., 0, 1]) / math.sin(half_lam_plus_phi), -1, 1)), ) phi = math.squeeze(phi % (4 * np.pi)) theta = math.squeeze(theta % (4 * np.pi)) lam = math.squeeze(lam % (4 * np.pi)) return (lam, theta, phi, alpha) if return_global_phase else (lam, theta, phi)
[docs] def xzx_rotation_angles(U, return_global_phase=False): r"""Compute the rotation angles :math:`\lambda`, :math:`\theta`, and :math:`\phi` and the phase :math:`\alpha` of a 2x2 unitary matrix as a product of X and Z rotations in the form :math:`e^{i\alpha} RX(\phi) RZ(\theta) RX(\lambda)`. Args: U (array): 2x2 unitary matrix return_global_phase (bool): if True, returns the global phase as well. Returns: tuple: The rotation angles :math:`\lambda`, :math:`\theta`, and :math:`\phi` and the global phase :math:`\alpha` if ``return_global_phase=True``. """ U, global_phase = math.convert_to_su2(U, return_global_phase=True) EPS = math.finfo(U.dtype).eps # Compute \phi, \theta and \lambda after analytically solving for them from # U = RX(\phi) RZ(\theta) RX(\lambda) sum_diagonal_real = math.real(U[..., 0, 0] + U[..., 1, 1]) sum_off_diagonal_imag = math.imag(U[..., 0, 1] + U[..., 1, 0]) half_phi_plus_lambdas = math.arctan2(-sum_off_diagonal_imag, sum_diagonal_real + EPS) diff_diagonal_imag = math.imag(U[..., 0, 0] - U[..., 1, 1]) diff_off_diagonal_real = math.real(U[..., 0, 1] - U[..., 1, 0]) half_phi_minus_lambdas = math.arctan2(diff_off_diagonal_real, -diff_diagonal_imag + EPS) lam = half_phi_plus_lambdas - half_phi_minus_lambdas phi = half_phi_plus_lambdas + half_phi_minus_lambdas # Compute \theta theta = math.where( math.isclose(math.sin(half_phi_plus_lambdas), math.zeros_like(half_phi_plus_lambdas)), 2 * math.arccos( math.clip(sum_diagonal_real / (2 * math.cos(half_phi_plus_lambdas) + EPS), -1, 1) ), 2 * math.arccos( math.clip(-sum_off_diagonal_imag / (2 * math.sin(half_phi_plus_lambdas) + EPS), -1, 1) ), ) phi = math.squeeze(phi % (4 * np.pi)) theta = math.squeeze(theta % (4 * np.pi)) lam = math.squeeze(lam % (4 * np.pi)) return (lam, theta, phi, global_phase) if return_global_phase else (lam, theta, phi)
[docs] def zxz_rotation_angles(U, return_global_phase=False): r"""Compute the rotation angles :math:`\lambda`, :math:`\theta`, and :math:`\phi` and the phase :math:`\alpha` of a 2x2 unitary matrix as a product of Z and X rotations in the form :math:`e^{i\alpha} RZ(\phi) RX(\theta) RZ(\lambda)`. Args: U (array): 2x2 unitary matrix return_global_phase (bool): if True, returns the global phase as well. Returns: tuple: The rotation angles :math:`\lambda`, :math:`\theta`, and :math:`\phi` and the global phase :math:`\alpha` if ``return_global_phase=True``. """ U, global_phase = math.convert_to_su2(U, return_global_phase=True) EPS = math.finfo(U.dtype).eps abs_a = math.clip(math.abs(U[..., 0, 0]), 0, 1) abs_b = math.clip(math.abs(U[..., 0, 1]), 0, 1) theta = math.where(abs_a < abs_b, 2 * math.arccos(abs_a), 2 * math.arcsin(abs_b)) half_phi_plus_lam = math.angle(U[..., 1, 1] + EPS) half_phi_minus_lam = math.angle(1j * U[..., 1, 0] + EPS) phi = half_phi_plus_lam + half_phi_minus_lam lam = half_phi_plus_lam - half_phi_minus_lam # Normalize the angles phi = math.squeeze(phi % (4 * np.pi)) theta = math.squeeze(theta % (4 * np.pi)) lam = math.squeeze(lam % (4 * np.pi)) return (lam, theta, phi, global_phase) if return_global_phase else (lam, theta, phi)
[docs] def su2su2_to_tensor_products(U): r"""Given a matrix :math:`U = A \otimes B` in SU(2) x SU(2), extract A and B This process has been described in detail in the Appendix of Coffey & Deiotte https://link.springer.com/article/10.1007/s11128-009-0156-3 """ # First, write A = [[a1, a2], [-a2*, a1*]], which we can do for any SU(2) element. # Then, A \otimes B = [[a1 B, a2 B], [-a2*B, a1*B]] = [[C1, C2], [C3, C4]] # where the Ci are 2x2 matrices. C1 = U[0:2, 0:2] C2 = U[0:2, 2:4] C3 = U[2:4, 0:2] C4 = U[2:4, 2:4] # From the definition of A \otimes B, C1 C4^\dag = a1^2 I, so we can extract a1 C14 = math.dot(C1, math.conj(math.T(C4))) a1 = math.sqrt(math.cast_like(C14[0, 0], 1j)) # Similarly, -C2 C3^\dag = a2^2 I, so we can extract a2 C23 = math.dot(C2, math.conj(math.T(C3))) a2 = math.sqrt(-math.cast_like(C23[0, 0], 1j)) # This gets us a1, a2 up to a sign. To resolve the sign, ensure that # C1 C2^dag = a1 a2* I C12 = math.dot(C1, math.conj(math.T(C2))) a2 = math.cond(math.allclose(a1 * math.conj(a2), C12[0, 0]), lambda: a2, lambda: -1 * a2, ()) # Construct A A = math.stack([math.stack([a1, a2]), math.stack([-math.conj(a2), math.conj(a1)])]) # Next, extract B. Can do from any of the C, just need to be careful in # case one of the elements of A is 0. # We use B1 unless division by 0 would cause all elements to be inf. B = math.cond( math.allclose(a1, 0.0, atol=1e-6), lambda: C2 / math.cast_like(a2, 1j), lambda: C1 / math.cast_like(a1, 1j), (), ) return math.convert_like(A, U), math.convert_like(B, U)
[docs] def decomp_int_to_powers_of_two(k: int, n: int) -> list[int]: r"""Decompose an integer :math:`k<=2^{n-1}` into additions and subtractions of the smallest-possible number of powers of two. Args: k (int): Integer to be decomposed n (int): Number of bits to consider Returns: list[int]: A list with length ``n``, with entry :math:`c_i` at position :math:`i`. This function is documented in ``pennylane/ops/qubit/pcphase_decomposition.md``. As an example, consider the number :math:`k=121_{10}=01111001_2`, which can be (trivially) decomposed into a sum of five powers of two by reading off the bits: :math:`k = 2^6 + 2^5 + 2^4 + 2^3 + 2^0 = 64 + 32 + 16 + 8 + 1`. The decomposition here, however, allows for minus signs and achieves the decomposition :math:`k = 2^7 - 2^3 + 2^0 = 128 - 8 + 1`, which only requires three powers of two. """ R = [] assert k <= 2 ** (n - 1) s = 0 powers = 2 ** np.arange(n) for p in powers: # p = 2**(n-1-i) if s & p == k & p: # Equal bit, move on factor = 0 else: # Differing bit, consider pairs of bits if p >= 2 ** (n - 2): # 2**(n-1-i) >= 2**(n-2) is the same condition as i < 2 factor = 1 else: # Table entry from documentation in_middle_rows = (s & (p + 2 * p)).bit_count() == 1 # two bits of s are 01 or 10 in_last_cols = bool(k & (2 * p)) # latter bit of k is 1 if in_middle_rows != in_last_cols: # xor between in_middle_rows and in_last_cols factor = -1 else: factor = 1 s += factor * p R.insert(0, factor) return R
def _set_unitary_matrix(unitary_matrix, index, value, like=None): """Set the values in the ``unitary_matrix`` at the specified index. Args: unitary_matrix (tensor_like): unitary being modified index (Tuple[Int | Ellipsis | List[Int]]): index for slicing the unitary value (tensor_like): new values for the specified index like (str): interface for the unitary matrix Returns: tensor_like: modified unitary Examples: A = np.eye(5) A = _set_unitary_matrix(A, (0, 0), 5) A = _set_unitary_matrix(A, (1, Ellipsis), np.array([1, 2, 3, 4, 5])) A = _set_unitary_matrix(A, (1, [1, 2]), np.array([3, 4])) """ if like is None: like = math.get_interface(unitary_matrix) if like == "jax": return unitary_matrix.at[index[0], index[1]].set( value, indices_are_sorted=True, unique_indices=True ) unitary_matrix[index[0], index[1]] = value return unitary_matrix def _givens_matrix_core(a, b, left=True, tol=1e-8): r"""Build a :math:`2 \times 2` Givens rotation matrix :math:`G`. When the matrix :math:`G` is applied to a vector :math:`[a,\ b]^T` the following would happen: .. math:: G \times \begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} 0 \\ r \end{bmatrix} \quad \quad \quad \begin{bmatrix} a \\ b \end{bmatrix} \times G = \begin{bmatrix} r \\ 0 \end{bmatrix}, where :math:`r` is a complex number. Args: a (float or complex): first element of the vector for which the Givens matrix is being computed b (float or complex): second element of the vector for which the Givens matrix is being computed left (bool): determines if the Givens matrix is being applied from the left side or right side. tol (float): determines tolerance limits for :math:`|a|` and :math:`|b|` under which they are considered as zero. Returns: tensor_like: Givens rotation matrix Raises: TypeError: if a and b have different interfaces """ abs_a, abs_b = math.abs(a), math.abs(b) interface_a, interface_b = math.get_interface(a), math.get_interface(b) if interface_a != interface_b: raise TypeError( f"The interfaces of 'a' and 'b' do not match. Got {interface_a} and {interface_b}." ) interface = interface_a aprod = math.nan_to_num(abs_b * abs_a) hypot = math.hypot(abs_a, abs_b) + 1e-15 # avoid division by zero cosine = math.where(abs_b < tol, 0.0, abs_b / hypot) cosine = math.where(abs_a < tol, 1.0, cosine) sine = math.where(abs_b < tol, 1.0, abs_a / hypot) sine = math.where(abs_a < tol, 0.0, sine) phase = math.where(abs_b < tol, 1.0, (1.0 * b * math.conj(a)) / (aprod + 1e-15)) phase = math.where(abs_a < tol, 1.0, phase) if interface == "jax": return jax.lax.cond( left, lambda phase, cosine, sine: math.array( [[phase * cosine, -sine], [phase * sine, cosine]], like=interface ), lambda phase, cosine, sine: math.array( [[phase * sine, cosine], [-phase * cosine, sine]], like=interface ), phase, cosine, sine, ) if left: return math.array([[phase * cosine, -sine], [phase * sine, cosine]], like=interface) return math.array([[phase * sine, cosine], [-phase * cosine, sine]], like=interface) @functools.lru_cache def _givens_matrix_jax(): @jax.jit def givens_matrix_jax(a, b, left=True, tol=1e-8): return _givens_matrix_core(a, b, left=left, tol=tol) return givens_matrix_jax def _givens_matrix(a, b, left=True, tol=1e-8): interface = math.get_interface(a) if interface == "jax" and isinstance(jnp.array(0), jax.core.Tracer): return _givens_matrix_jax()(a, b, left=left, tol=tol) return _givens_matrix_core(a, b, left=left, tol=tol) def _right_givens_core(indices, unitary, N, j): interface = math.get_interface(unitary) grot_mat = _givens_matrix(*unitary[N - j - 1, indices].T, left=True) unitary = _set_unitary_matrix( unitary, (Ellipsis, indices), unitary[:, indices] @ grot_mat.T, like=interface ) return unitary, math.conj(grot_mat) @functools.lru_cache def _right_givens_jax(): @jax.jit def _right_givens_jax(indices, unitary, N, j): return _right_givens_core(indices, unitary, N, j) return _right_givens_jax def _right_givens(indices, unitary, N, j): interface = math.get_interface(unitary) if interface == "jax" and isinstance(jnp.array(0), jax.core.Tracer): return _right_givens_jax()(indices, unitary, N, j) return _right_givens_core(indices, unitary, N, j) def _left_givens_core(indices, unitary, j): interface = math.get_interface(unitary) grot_mat = _givens_matrix(*unitary[indices, j - 1], left=False) unitary = _set_unitary_matrix( unitary, (indices, Ellipsis), grot_mat @ unitary[indices, :], like=interface ) return unitary, grot_mat @functools.lru_cache def _left_givens_jax(): @jax.jit def _left_givens_jax(indices, unitary, j): return _left_givens_core(indices, unitary, j) return _left_givens_jax def _left_givens(indices, unitary, j): interface = math.get_interface(unitary) if interface == "jax" and isinstance(jnp.array(0), jax.core.Tracer): return _left_givens_jax()(indices, unitary, j) return _left_givens_core(indices, unitary, j) # pylint: disable=too-many-branches
[docs] def givens_decomposition(unitary): r"""Decompose a unitary into a sequence of Givens rotation gates with phase shifts and a diagonal phase matrix. This decomposition is based on the construction scheme given in `Optica, 3, 1460 (2016) <https://opg.optica.org/optica/fulltext.cfm?uri=optica-3-12-1460&id=355743>`_\ , which allows one to write any unitary matrix :math:`U` as: .. math:: U = D \left(\prod_{(m, n) \in G} T_{m, n}(\theta, \phi)\right), where :math:`D` is a diagonal phase matrix, :math:`T(\theta, \phi)` is the Givens rotation gates with phase shifts and :math:`G` defines the specific ordered sequence of the Givens rotation gates acting on wires :math:`(m, n)`. The unitary for the :math:`T(\theta, \phi)` gates appearing in the decomposition is of the following form: .. math:: T(\theta, \phi) = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & e^{i \phi} \cos(\theta) & -\sin(\theta) & 0 \\ 0 & e^{i \phi} \sin(\theta) & \cos(\theta) & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}, where :math:`\theta \in [0, \pi/2]` is the angle of rotation in the :math:`\{|01\rangle, |10\rangle \}` subspace and :math:`\phi \in [0, 2 \pi]` represents the phase shift at the first wire. **Example** .. code-block:: python unitary = np.array([[ 0.73678+0.27511j, -0.5095 +0.10704j, -0.06847+0.32515j], [-0.21271+0.34938j, -0.38853+0.36497j, 0.61467-0.41317j], [ 0.41356-0.20765j, -0.00651-0.66689j, 0.32839-0.48293j]]) phase_mat, ordered_rotations = givens_decomposition(unitary) >>> phase_mat tensor([-0.20604358+0.9785369j , -0.82993272+0.55786114j, 0.56230612-0.82692833j], requires_grad=True) >>> ordered_rotations [(tensor([[-0.65087861-0.63937521j, -0.40933651-0.j ], [-0.29201359-0.28685265j, 0.91238348-0.j ]], requires_grad=True), (0, 1)), (tensor([[ 0.47970366-0.33308926j, -0.8117487 -0.j ], [ 0.66677093-0.46298215j, 0.5840069 -0.j ]], requires_grad=True), (1, 2)), (tensor([[ 0.36147547+0.73779454j, -0.57008306-0.j ], [ 0.2508207 +0.51194108j, 0.82158706-0.j ]], requires_grad=True), (0, 1))] Args: unitary (tensor): unitary matrix on which decomposition will be performed Returns: (tensor_like, list[(tensor_like, tuple)]): diagonal elements of the phase matrix :math:`D` and Givens rotation matrix :math:`T` with their indices Raises: ValueError: if the provided matrix is not square. .. details:: :title: Theory and Pseudocode **Givens Rotation** Applying the Givens rotation :math:`T(\theta, \phi)` performs the following transformation of the basis states: .. math:: &|00\rangle \mapsto |00\rangle\\ &|01\rangle \mapsto e^{i \phi} \cos(\theta) |01\rangle - \sin(\theta) |10\rangle\\ &|10\rangle \mapsto e^{i \phi} \sin(\theta) |01\rangle + \cos(\theta) |10\rangle\\ &|11\rangle \mapsto |11\rangle. **Pseudocode** The algorithm that implements the decomposition is the following: .. code-block:: python def givens_decomposition(U): for i in range(1, N): if i % 2: for j in range(0, i): # Find T^-1(i-j, i-j+1) matrix that nulls element (N-j, i-j) of U # Update U = U @ T^-1(i-j, i-j+1) else: for j in range(1, i): # Find T(N+j-i-1, N+j-i) matrix that nulls element (N+j-i, j) of U # Update U = T(N+j-i-1, N+j-i) @ U """ interface = math.get_deep_interface(unitary) unitary = math.copy(unitary) if interface == "jax" else math.toarray(unitary).copy() M, N = math.shape(unitary) if M != N: raise ValueError(f"The unitary matrix should be of shape NxN, got {unitary.shape}") left_givens, right_givens = [], [] for i in range(1, N): if i % 2: for j in range(0, i): indices = [i - j - 1, i - j] unitary, grot_mat_conj = _right_givens(indices, unitary, N, j) right_givens.append((grot_mat_conj, indices)) else: for j in range(1, i + 1): indices = [N + j - i - 2, N + j - i - 1] unitary, grot_mat = _left_givens(indices, unitary, j) left_givens.append((grot_mat, indices)) nleft_givens = [] for grot_mat, (i, j) in reversed(left_givens): sphase_mat = math.diag(math.diag(unitary)[math.array([i, j])]) decomp_mat = math.conj(grot_mat).T @ sphase_mat givens_mat = _givens_matrix(*decomp_mat[1, :].T) nphase_mat = decomp_mat @ givens_mat.T # check for T_{m,n}^{-1} x D = D x T. if not math.is_abstract(decomp_mat) and not math.allclose( nphase_mat @ math.conj(givens_mat), decomp_mat ): # pragma: no cover raise ValueError("Failed to shift phase transposition.") for diag_idx, diag_val in zip([(i, i), (j, j)], math.diag(nphase_mat)): unitary = _set_unitary_matrix(unitary, diag_idx, diag_val, like=interface) nleft_givens.append((math.conj(givens_mat), (i, j))) phases, ordered_rotations = math.diag(unitary), [] for grot_mat, (i, j) in list(reversed(nleft_givens)) + list(reversed(right_givens)): if not math.is_abstract(grot_mat) and not math.all( math.isreal(grot_mat[0, 1]) and math.isreal(grot_mat[1, 1]) ): # pragma: no cover raise ValueError(f"Incorrect Givens Rotation encountered, {grot_mat}") ordered_rotations.append((grot_mat, (i, j))) return phases, ordered_rotations