Source code for qiskit_nature.second_q.circuit.library.bogoliubov_transform

# This code is part of a Qiskit project.
#
# (C) Copyright IBM 2022, 2023.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

"""Bogoliubov transform."""

from __future__ import annotations

from collections.abc import Iterator

import numpy as np
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.circuit import Gate, Qubit
from qiskit.circuit.library import RZGate, XXPlusYYGate
from qiskit_nature.second_q.mappers import QubitMapper
from qiskit_nature.second_q.mappers import JordanWignerMapper
from qiskit_nature.utils import apply_matrix_to_slices, givens_matrix
from qiskit_nature.utils.linalg import fermionic_gaussian_decomposition_jw


def _rows_are_orthonormal(mat: np.ndarray, rtol: float = 1e-5, atol: float = 1e-8) -> bool:
    m, _ = mat.shape
    return np.allclose(mat @ mat.T.conj(), np.eye(m), rtol=rtol, atol=atol)


def _validate_transformation_matrix(
    mat: np.ndarray, rtol: float = 1e-5, atol: float = 1e-8
) -> None:
    if not len(mat.shape) == 2:
        raise ValueError(
            "transformation_matrix must be a 2-dimensional array. "
            f"Instead, got shape {mat.shape}."
        )
    n, p = mat.shape  # pylint: disable=invalid-name
    if p == n:
        if not _rows_are_orthonormal(mat, rtol=rtol, atol=atol):
            raise ValueError("transformation_matrix must have orthonormal rows.")
    elif p == n * 2:
        left = mat[:, :n]
        right = mat[:, n:]
        comm1 = left @ left.T.conj() + right @ right.T.conj()
        comm2 = left @ right.T + right @ left.T
        if not np.allclose(comm1, np.eye(n), rtol=rtol, atol=atol) or not np.allclose(
            comm2, 0.0, atol=atol
        ):
            raise ValueError(
                "transformation_matrix does not describe a valid transformation "
                "of fermionic ladder operators. A valid matrix should have the block form "
                "[W1 W2] where W1 @ W1.T.conj() + W2 @ W2.T.conj() = I and "
                "W1 @ W2.T + W2 @ W1.T = 0."
            )
    else:
        raise ValueError(
            f"transformation_matrix must be N x N or N x 2N. Instead, got shape {mat.shape}."
        )


[docs]class BogoliubovTransform(QuantumCircuit): r"""A circuit that performs a Bogoliubov transform. A Bogoliubov transform effects a unitary basis change that maps the fermionic ladder operators to a new set of ladder operators that also satisfy the fermionic anticommutation relations. That is, it effects a unitary :math:`U` such that .. math:: U a^\dagger_j U^\dagger = b^\dagger_j, \quad j = 1, \ldots, N where the :math:`\{a_j\}` are the original fermionic creation operators and the :math:`\{b_j\}` are the new fermionic creation operators. The new creation operators are linear combinations of the original ladder operators, and the coefficients of the linear combinations are specified by a matrix :math:`W` which determines the unitary :math:`U`. The matrix :math:`W` is either :math:`N \times N` or :math:`N \times 2N`. If :math:`W` is :math:`N \times N`, then the linear combinations involve only the original creation operators: .. math:: \begin{pmatrix} b^\dagger_1 \\ \vdots \\ b^\dagger_N \\ \end{pmatrix} = W \begin{pmatrix} a^\dagger_1 \\ \vdots \\ a^\dagger_N \\ \end{pmatrix}. If :math:`W` is :math:`N \times 2N`, then the linear combinations involve both the original creation and annihilation operators: .. math:: \begin{pmatrix} b^\dagger_1 \\ \vdots \\ b^\dagger_N \\ \end{pmatrix} = W \begin{pmatrix} a^\dagger_1 \\ \vdots \\ a^\dagger_N \\ a_1 \\ \vdots \\ a_N \end{pmatrix}. The matrix :math:`W` is commonly obtained by calling the :meth:`~.QuadraticHamiltonian.diagonalizing_bogoliubov_transform` method of the :class:`~.QuadraticHamiltonian` class. Currently, only the Jordan-Wigner Transformation is supported. References: - `arXiv:1711.05395`_ - `arXiv:1603.08788`_ .. _arXiv:1711.05395: https://arxiv.org/abs/1711.05395 .. _arXiv:1603.08788: https://arxiv.org/abs/1603.08788 """ def __init__( self, transformation_matrix: np.ndarray, qubit_mapper: QubitMapper | None = None, *, validate: bool = True, rtol: float = 1e-5, atol: float = 1e-8, **circuit_kwargs, ) -> None: # pylint: disable=unused-argument r""" Args: transformation_matrix: The matrix :math:`W` that specifies the coefficients of the new creation operators in terms of the original creation operators. Should be either :math:`N \times N` or :math:`N \times 2N`. qubit_mapper: The ``QubitMapper``. The default behavior is to create one using the call ``JordanWignerMapper()``. validate: Whether to validate the inputs. rtol: Relative numerical tolerance for input validation. atol: Absolute numerical tolerance for input validation. circuit_kwargs: Keyword arguments to pass to the ``QuantumCircuit`` initializer. Raises: ValueError: ``transformation_matrix`` must be a 2-dimensional array. ValueError: ``transformation_matrix`` must have orthonormal rows. ValueError: ``transformation_matrix`` does not describe a valid transformation of fermionic ladder operators. If the transformation matrix is :math:`N \times N`, then it should be unitary. If the transformation matrix is :math:`N \times 2N`, then it should have the block form :math:`(W_1 \quad W_2)` where :math:`W_1 W_1^\dagger + W_2 W_2^\dagger = I` and :math:`W_1 W_2^T + W_2 W_1^T = 0`. NotImplementedError: Currently, only the Jordan-Wigner Transform is supported. Please use the :class:`qiskit_nature.second_q.mappers.JordanWignerMapper`. """ if validate: _validate_transformation_matrix(transformation_matrix, rtol=rtol, atol=atol) if qubit_mapper is None: qubit_mapper = JordanWignerMapper() n, _ = transformation_matrix.shape register = QuantumRegister(n) super().__init__(register, **circuit_kwargs) if isinstance(qubit_mapper, JordanWignerMapper): operations = _bogoliubov_transform_jw(register, transformation_matrix) for gate, qubits in operations: self.append(gate, qubits) else: raise NotImplementedError( "Currently, only the Jordan-Wigner Transform is supported. " "Please use the qiskit_nature.second_q.mappers.JordanWignerMapper." )
def _bogoliubov_transform_jw( register: QuantumRegister, transformation_matrix: np.ndarray ) -> Iterator[tuple[Gate, tuple[Qubit, ...]]]: n, p = transformation_matrix.shape # pylint: disable=invalid-name if p == n: yield from _bogoliubov_transform_num_conserving_jw(register, transformation_matrix) else: yield from _bogoliubov_transform_general_jw(register, transformation_matrix) def _bogoliubov_transform_num_conserving_jw( # pylint: disable=invalid-name register: QuantumRegister, transformation_matrix: np.ndarray ) -> Iterator[tuple[Gate, tuple[Qubit, ...]]]: n, _ = transformation_matrix.shape current_matrix = transformation_matrix left_rotations = [] right_rotations = [] # compute left and right Givens rotations for i in range(n - 1): if i % 2 == 0: # rotate columns by right multiplication for j in range(i + 1): target_index = i - j row = n - j - 1 if not np.isclose(current_matrix[row, target_index], 0.0): # zero out element at target index in given row givens_mat = givens_matrix( current_matrix[row, target_index + 1], current_matrix[row, target_index], ) right_rotations.append((givens_mat, (target_index + 1, target_index))) current_matrix = apply_matrix_to_slices( current_matrix, givens_mat, [(Ellipsis, target_index + 1), (Ellipsis, target_index)], ) else: # rotate rows by left multiplication for j in range(i + 1): target_index = n - i + j - 1 col = j if not np.isclose(current_matrix[target_index, col], 0.0): # zero out element at target index in given column givens_mat = givens_matrix( current_matrix[target_index - 1, col], current_matrix[target_index, col], ) left_rotations.append((givens_mat, (target_index - 1, target_index))) current_matrix = apply_matrix_to_slices( current_matrix, givens_mat, [target_index - 1, target_index] ) # convert left rotations to right rotations for givens_mat, (i, j) in reversed(left_rotations): givens_mat = givens_mat.T.conj().astype(complex, copy=False) givens_mat[:, 0] *= current_matrix[i, i] givens_mat[:, 1] *= current_matrix[j, j] new_givens_mat = givens_matrix(givens_mat[1, 1], givens_mat[1, 0]) right_rotations.append((new_givens_mat.T, (i, j))) phase_matrix = givens_mat @ new_givens_mat current_matrix[i, i] = phase_matrix[0, 0] current_matrix[j, j] = phase_matrix[1, 1] # yield operations for i in range(n): phi = np.angle(current_matrix[i, i]) yield RZGate(phi), (register[i],) for givens_mat, (i, j) in reversed(right_rotations): theta = np.arccos(np.real(givens_mat[0, 0])) phi = np.angle(givens_mat[0, 1]) yield XXPlusYYGate(2 * theta, phi - np.pi / 2), (register[j], register[i]) def _bogoliubov_transform_general_jw( # pylint: disable=invalid-name register: QuantumRegister, transformation_matrix: np.ndarray ) -> Iterator[tuple[Gate, tuple[Qubit, ...]]]: decomposition, left_unitary = fermionic_gaussian_decomposition_jw( register, transformation_matrix ) yield from _bogoliubov_transform_num_conserving_jw(register, left_unitary.T) yield from reversed(decomposition)