# Source code for pennylane.qchem.givens_decomposition

# Copyright 2018-2023 Xanadu Quantum Technologies Inc.

# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
"""
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,
>>> ordered_rotations
[(tensor([[-0.65087861-0.63937521j, -0.40933651-0.j        ],
(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


Using PennyLane

Release news

Development

API

Internals