# This code is part of a Qiskit project.
#
# (C) Copyright IBM 2021, 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.
"""The QuadraticHamiltonian class."""
from __future__ import annotations
from typing import cast
import numpy as np
import scipy.linalg
from qiskit_nature.second_q.operators import FermionicOp, PolynomialTensor
from .hamiltonian import Hamiltonian
def _is_hermitian(mat: np.ndarray, rtol: float = 1e-5, atol: float = 1e-8) -> bool:
return np.allclose(mat, mat.T.conj(), rtol=rtol, atol=atol)
def _is_antisymmetric(mat: np.ndarray, rtol: float = 1e-5, atol: float = 1e-8) -> bool:
return np.allclose(mat, -mat.T, rtol=rtol, atol=atol)
def _infer_num_modes(
hermitian_part: np.ndarray | None, antisymmetric_part: np.ndarray | None, num_modes: int | None
) -> int:
if num_modes is not None:
return num_modes
if hermitian_part is not None:
return hermitian_part.shape[0]
if antisymmetric_part is not None:
return antisymmetric_part.shape[0]
else:
raise ValueError(
"Either Hermitian part, antisymmetric part, or number of modes must be specified."
)
def _validate(
hermitian_part: np.ndarray | None,
antisymmetric_part: np.ndarray | None,
num_modes: int | None,
rtol: float,
atol: float,
) -> None:
if (
hermitian_part is not None
and antisymmetric_part is not None
and hermitian_part.shape != antisymmetric_part.shape
):
raise ValueError(
"Hermitian part and antisymmetric part must have same shape. "
f"Got shapes {hermitian_part.shape} and {antisymmetric_part.shape}."
)
if hermitian_part is not None:
if hermitian_part.shape[0] != num_modes:
raise ValueError(
"Hermitian part must have shape num_modes x num_modes. "
f"Got shape {hermitian_part.shape}, while num_modes={num_modes}."
)
if not _is_hermitian(hermitian_part, rtol=rtol, atol=atol):
raise ValueError("Hermitian part must be Hermitian.")
if antisymmetric_part is not None:
if antisymmetric_part.shape[0] != num_modes:
raise ValueError(
"Antisymmetric part must have shape num_modes x num_modes. "
f"Got shape {antisymmetric_part.shape}, while num_modes={num_modes}."
)
if not _is_antisymmetric(antisymmetric_part, rtol=rtol, atol=atol):
raise ValueError("Antisymmetric part must be antisymmetric.")
[docs]class QuadraticHamiltonian(PolynomialTensor, Hamiltonian):
r"""A Hamiltonian that is quadratic in the fermionic ladder operators.
A quadratic Hamiltonian is an operator of the form
.. math::
\sum_{p, q} M_{pq} a^\dagger_p a_q
+ \frac12 \sum_{p, q}
(\Delta_{pq} a^\dagger_p a^\dagger_q + \text{h.c.})
+ \text{constant},
where :math:`M` is a Hermitian matrix and :math:`\Delta` is an antisymmetric matrix.
Note:
The :class:`~.FermionicOp` class can also be used to represent any quadratic Hamiltonian.
The reason to have a class specifically for quadratic Hamiltonians is that they
support special numerical routines that involve performing linear algebra on the
matrices :math:`M` and :math:`\Delta`. The internal representation format
of :class:`~.FermionicOp` is not suitable for these routines.
"""
def __init__(
self,
hermitian_part: np.ndarray | None = None,
antisymmetric_part: np.ndarray | None = None,
constant: float = 0.0,
*,
num_modes: int | None = None,
validate: bool = True,
rtol: float | None = None,
atol: float | None = None,
) -> None:
r"""
Args:
hermitian_part: The matrix :math:`M` containing the coefficients of the terms
that conserve particle number.
antisymmetric_part: The matrix :math:`\Delta` containing the coefficients of
the terms that do not conserve particle number.
constant: An additive constant term.
num_modes: Number of fermionic modes. This should be consistent with ``hermitian_part``
and ``antisymmetric_part`` if they are specified.
validate: Whether to validate the inputs.
rtol: Relative numerical tolerance for input validation.
The default behavior is to use ``self.rtol``.
atol: Absolute numerical tolerance for input validation.
The default behavior is to use ``self.atol``.
Raises:
ValueError: Either Hermitian part, antisymmetric part, or number of modes must
be specified.
ValueError: Hermitian part and antisymmetric part must have same shape.
ValueError: Hermitian part must have shape num_modes x num_modes.
ValueError: Hermitian part must be Hermitian.
ValueError: Antisymmetric part must have shape num_modes x num_modes.
ValueError: Antisymmetric part must be antisymmetric.
Note:
The ``rtol`` and ``atol`` arguments are only used for input validation and
are discarded afterwards. They do not affect the class attributes
``QuadraticHamiltonian.rtol`` and ``QuadraticHamiltonian.atol``.
"""
if rtol is None:
rtol = self.rtol
if atol is None:
atol = self.atol
self._num_modes = _infer_num_modes(hermitian_part, antisymmetric_part, num_modes)
if validate:
_validate(hermitian_part, antisymmetric_part, self._num_modes, rtol, atol)
data: dict[str, float | np.ndarray] = {"": constant}
if hermitian_part is not None:
data["+-"] = hermitian_part
if antisymmetric_part is not None:
data["++"] = 0.5 * antisymmetric_part
data["--"] = -0.5 * antisymmetric_part.conj()
super().__init__(data, validate=validate)
def __repr__(self) -> str:
return (
f"QuadraticHamiltonian("
f"hermitian_part={repr(self.hermitian_part)}, "
f"antisymmetric_part={repr(self.antisymmetric_part)}, "
f"constant={self.constant}, "
f"num_modes={self._num_modes})"
)
@property
def hermitian_part(self) -> np.ndarray:
"""The matrix of coefficients of terms that conserve particle number."""
if "+-" in self:
return cast(np.ndarray, self["+-"])
return np.zeros((self._num_modes, self._num_modes), dtype=complex)
@property
def antisymmetric_part(self) -> np.ndarray:
"""The matrix of coefficients of terms that do not conserve particle number."""
if "++" in self:
return 2 * cast(np.ndarray, self["++"])
return np.zeros((self._num_modes, self._num_modes), dtype=complex)
@property
def constant(self) -> float:
"""The constant."""
return cast(float, self[""])
@property
def num_modes(self) -> float:
"""The number of modes this operator acts on."""
return self._num_modes
@property
def register_length(self) -> int:
return self._num_modes
[docs] def second_q_op(self) -> FermionicOp:
"""Convert to FermionicOp."""
return FermionicOp.from_polynomial_tensor(self)
[docs] def conserves_particle_number(self) -> bool:
"""Whether the Hamiltonian conserves particle number."""
return np.allclose(self.antisymmetric_part, 0.0)
def _particle_num_conserving_bogoliubov_transform(self) -> tuple[np.ndarray, np.ndarray, float]:
orbital_energies, basis_change = np.linalg.eigh(self.hermitian_part)
return basis_change.T, orbital_energies, self.constant
def _non_particle_num_conserving_bogoliubov_transform(
self,
) -> tuple[np.ndarray, np.ndarray, float]:
matrix, constant = self.majorana_form()
canonical_form, basis_change = _antisymmetric_canonical_form(matrix)
orbital_energies = canonical_form[
range(self._num_modes), range(self._num_modes, 2 * self._num_modes)
]
constant -= 0.5 * np.sum(orbital_energies)
eye = np.eye(self._num_modes, dtype=complex)
majorana_basis = np.block([[eye, eye], [1j * eye, -1j * eye]]) / np.sqrt(2)
diagonalizing_unitary = majorana_basis.T.conj() @ basis_change @ majorana_basis
return diagonalizing_unitary[: self._num_modes], orbital_energies, constant
def _antisymmetric_canonical_form(matrix: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
"""Put an antisymmetric matrix into canonical form.
The input is an antisymmetric matrix A with even dimension.
Its canonical form is
```
A = R^T C R
```
where R is a real orthogonal matrix and C has the form
```
C = [ 0 D ]
[ -D 0 ]
```
where D is a real diagonal matrix with non-negative entries
sorted in ascending order. This form is equivalent to the Schur form
up to a permutation.
Returns:
The matrices C and R.
"""
dim, _ = matrix.shape
canonical_form, basis_change = scipy.linalg.schur(matrix, output="real")
# shift 2x2 blocks so they lie on even indices
permutation = np.arange(dim)
for i in range(1, dim - 1, 2):
if not np.isclose(canonical_form[i + 1, i], 0.0):
_swap_indices(permutation, i - 1, i + 1)
_permute_indices(canonical_form, permutation)
basis_change = basis_change[:, permutation]
# permute matrix into [[0, A], [-A, 0]] form
permutation = _schur_permutation(dim)
_permute_indices(canonical_form, permutation)
basis_change = basis_change[:, permutation]
# permute matrix so that the upper right block is non-negative
n = dim // 2
permutation = np.arange(dim)
for i in range(n):
if canonical_form[i, n + i] < 0.0:
_swap_indices(permutation, i, n + i)
_permute_indices(canonical_form, permutation)
basis_change = basis_change[:, permutation]
# permute matrix so that the energies are sorted
permutation = np.argsort(canonical_form[range(n), range(n, 2 * n)])
permutation = np.concatenate([permutation, permutation + n])
_permute_indices(canonical_form, permutation)
basis_change = basis_change[:, permutation]
return canonical_form, basis_change.T
def _schur_permutation(dim: int) -> np.ndarray:
n = dim // 2
permutation = np.arange(dim)
for i in range(1, n, 2):
_swap_indices(permutation, i, n + i - 1)
if n % 2 != 0:
_swap_indices(permutation, n - 1, n + i)
return permutation
def _swap_indices(array: np.ndarray, i: int, j: int) -> None:
array[i], array[j] = array[j], array[i]
def _permute_indices(matrix: np.ndarray, permutation: np.ndarray) -> None:
n, _ = matrix.shape
for i in range(n):
matrix[i, :] = matrix[i, permutation]
for i in range(n):
matrix[:, i] = matrix[permutation, i]