Source code for pennylane.qchem.givens_decomposition
# Copyright 2018-2023 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.
"""
This module contains the functions needed for performing basis transformations defined by a set of fermionic ladder operators.
"""
import numpy as np
import pennylane as qml
def _givens_matrix(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:
np.ndarray (or tensor): Givens rotation matrix
"""
abs_a, abs_b = np.abs(a), np.abs(b)
if abs_a < tol:
cosine, sine, phase = 1.0, 0.0, 1.0
elif abs_b < tol:
cosine, sine, phase = 0.0, 1.0, 1.0
else:
hypot = np.hypot(abs_a, abs_b)
cosine = abs_b / hypot
sine = abs_a / hypot
phase = 1.0 * b / abs_b * a.conjugate() / abs_a
if left:
return np.array([[phase * cosine, -sine], [phase * sine, cosine]])
return np.array([[phase * sine, cosine], [-phase * cosine, sine]])
[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:
(np.ndarray, list[(np.ndarray, 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
"""
unitary, (M, N) = qml.math.toarray(unitary).copy(), unitary.shape
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]
grot_mat = _givens_matrix(*unitary[N - j - 1, indices].T, left=True)
unitary[:, indices] = unitary[:, indices] @ grot_mat.T
right_givens.append((grot_mat.conj(), indices))
else:
for j in range(1, i + 1):
indices = [N + j - i - 2, N + j - i - 1]
grot_mat = _givens_matrix(*unitary[indices, j - 1], left=False)
unitary[indices] = grot_mat @ unitary[indices]
left_givens.append((grot_mat, indices))
nleft_givens = []
for grot_mat, (i, j) in reversed(left_givens):
sphase_mat = np.diag(np.diag(unitary)[[i, j]])
decomp_mat = grot_mat.conj().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 np.allclose(nphase_mat @ givens_mat.conj(), decomp_mat): # pragma: no cover
raise ValueError("Failed to shift phase transposition.")
unitary[i, i], unitary[j, j] = np.diag(nphase_mat)
nleft_givens.append((givens_mat.conj(), (i, j)))
phases, ordered_rotations = np.diag(unitary), []
for grot_mat, (i, j) in list(reversed(nleft_givens)) + list(reversed(right_givens)):
if not np.all(np.isreal(grot_mat[0, 1]) and np.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
_modules/pennylane/qchem/givens_decomposition
Download Python script
Download Notebook
View on GitHub