# This code is part of qtealeaves.
#
# 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 module contains tools for the QUBO solvers in :mod:`qubo_solver`.
"""
# pylint: disable=too-many-locals
# pylint: disable=too-many-arguments
# Modules
# ----------------------------------------------------------------------------
from typing import TYPE_CHECKING, Any, TypeAlias
import numpy as np
import numpy.typing as npt
import scipy.sparse as ssp
from qtealeaves.convergence_parameters import TNConvergenceParameters
from qtealeaves.emulator import MPS, TTN
from qtealeaves.tensors import TensorBackend
__all__ = [
"spins_to_binaries",
"construct_exact_observables",
"construct_exact_spinglass_hamiltonian",
"get_exact_spectrum",
"measure_exact_observables",
"generate_perturbation_transverse_fields",
"construct_spinglass_driver_state",
"read_tn_state",
]
# Resolves sphinx problems at import with |, alternative future / annotations
if TYPE_CHECKING:
NDArrayOrListInt: TypeAlias = npt.NDArray[np.integer] | list[int]
FloatOrArray: TypeAlias = float | npt.NDArray[np.floating]
SparseOrDense: TypeAlias = ssp.csr_matrix | npt.NDArray
ArrayOrNone: TypeAlias = npt.NDArray[np.floating] | None
ArrayOrListArray: TypeAlias = npt.NDArray | list[npt.NDArray]
else:
NDArrayOrListInt: TypeAlias = Any
FloatOrArray: TypeAlias = Any
SparseOrDense: TypeAlias = Any
ArrayOrNone: TypeAlias = Any
ArrayOrListArray: TypeAlias = Any
# ----------------------------------------------------------------------------
# ----------------------------------------------------------------------------
def spins_to_binaries(
configuration: NDArrayOrListInt,
) -> npt.NDArray[np.integer]:
r"""
Spin-to-binary converter.
Arguments
---------
configuration : numpy.ndarray[int] | list[int]
The string of +1 and -1 representing the
quantum expectation value of the Z-Pauli
matrix on each site (+1 means spin-up while
-1 mean spin-down).
Returns
-------
bitstring : numpy.ndarray[int]
The converted string of 0s and 1s.
Raises
------
:exc:`ValueError`
If invalid spin states are found.
Details
-------
Implements the conversion from qubit (spin-:math:`1 \\mathbin{/} 2`)
variables to binary variables, e.g., the decision variables of a
combinatorial optimization problem.
The applied linear transformation is
.. math::
x = \\dfrac{1 + \\sigma}{2}
x \\in \\left\{0, 1\\right\\} (binary)
\\sigma \\in \\left\{-1, +1\\right\\} (qubit)
that maps bit :math:`x = 0` into spin-down and bit :math:`x = 1`
into spin-up.
"""
# Read the input spin state
bitstring = np.asarray(configuration, dtype=int)
# Check non-valid values
if not np.all(np.isin(bitstring, [-1, 1])):
raise ValueError("Invalid measurement: spin-1/2 state must be up or down.")
# Output
return (bitstring + 1) // 2
# ==================================================
# Spin glass Hamiltonians: exact ground-state search
# ==================================================
def construct_exact_observables(
number_of_spins: int,
) -> dict[str, list[ssp.csr_matrix]]:
"""
Constructs the set of local operators to describe spin glass
observables to be measured.
Arguments
---------
number_of_spins : int
The number of spins in the system.
Returns
-------
dict[str, list[scipy.sparse.csr_matrix]]
Spin glass observables as a dictionary.
The name of the observables is the key, the matrix
representing the observable is the value.
Values are lists as we have local observables per
physical site, for a total of ``number_of_spins``
sites.
Details
-------
Operators are represented as sparse matrices whenever possible.
"""
# Single-qubit local operators
zlocal = ssp.csr_matrix([[1, 0], [0, -1]])
# Logical qubit-state observable
zproj = [
ssp.kron(
A=ssp.kron(ssp.identity(2**site, format="csr"), B=zlocal, format="csr"),
B=ssp.identity(2 ** (number_of_spins - site - 1), format="csr"),
format="csr",
)
for site in range(number_of_spins)
]
# Other observables to be implemented if needed
# For now, we need to measure only the projections
# onto the computational basis.
# Output
return {"sz": zproj}
def construct_exact_spinglass_hamiltonian(
spinglass_couplings: dict[str, FloatOrArray],
one_body_terms_prefactor: float | None = 1.0,
two_body_terms_prefactor: float | None = 1.0,
) -> ssp.csr_matrix:
"""
Constructs the exact matrix representation of a spin glass Hamiltonian
operator characterized by ``spinglass_couplings``.
Arguments
---------
spinglass_couplings : dict[str, float | numpy.ndarray[float]]
The set of couplings defining the spin glass Hamiltonian.
Specifically:
- **'offset'**: the constant term proportional to the identity.
This energy offset does not influence the spectrum of the spin
glass system, but it is necessary to reconstruct the exact
energy values of the eigenstates encoding spin configurations;
- **'one-body'**: the set of single-body couplings, i.e., the set
of local longitudinal magnetic fields (biases), one for each
spin-1/2 in the spin glass system;
- **'two-body'**: the set of two-body (in general all-to-all)
couplings describing the interactions between pairs of spin-1/2.
one_body_terms_prefactor : float, optional
The multiplicative prefactor of the spin glass
local biases terms.
Default to 1.0.
two_body_terms_prefactor : float, optional
The multiplicative prefactor of the spin-spin
interaction-strength terms.
Default to 1.0.
Returns
-------
ham_matrix : scipy.sparse.csr_matrix
The 2D scipy sparse matrix representing the spin glass
Hamiltonian matrix as a many-body operator.
"""
n_spins = spinglass_couplings["one-body"].size
hilbert_space_dim = 2**n_spins
zlocal = ssp.csr_matrix([[1, 0], [0, -1]])
# Build the Hamiltonian matrix
ham_matrix = ssp.csr_matrix((hilbert_space_dim, hilbert_space_dim))
# Add one-body terms
for site, bias in enumerate(spinglass_couplings["one-body"]):
ham_matrix += (
one_body_terms_prefactor
* bias
* ssp.kron(
A=ssp.kron(
A=ssp.identity(2**site, format="csr"), B=zlocal, format="csr"
),
B=ssp.identity(2 ** (n_spins - site - 1), format="csr"),
format="csr",
)
)
# Add two-spin interactions
jmat = two_body_terms_prefactor * spinglass_couplings["two-body"]
for site1 in range(n_spins - 1):
for site2 in range(site1 + 1, n_spins):
mel = jmat[site1, site2]
if mel != 0:
ham_matrix += mel * ssp.kron(
A=ssp.kron(
A=ssp.kron(
A=ssp.kron(
A=ssp.identity(2**site1, format="csr"),
B=zlocal,
format="csr",
),
B=ssp.identity(n=2 ** (site2 - site1 - 1), format="csr"),
format="csr",
),
B=zlocal,
format="csr",
),
B=ssp.identity(2 ** (n_spins - site2 - 1), format="csr"),
format="csr",
)
# Output
return ham_matrix
[docs]
def get_exact_spectrum(
ham_matrix: SparseOrDense,
number_of_eigenstates: int | None = 1,
full_spectrum: bool | None = False,
) -> list[list[float, npt.NDArray]]:
"""
Computes exact eigenstates and eigenvalues of a many-body **classical**
Hamiltonian.
Arguments
---------
ham_matrix : scipy.sparse.csr_matrix | numpy.ndarray
The matrix representing the Hamiltonian operator.
number_of_eigenstates : int, optional
The number of eigenstates to be computed.
Default to 1, i.e., only the ground state.
full_spectrum : bool, optional
If :obj:`True`, computes the full Hamiltonian spectrum,
ignoring ``number_of_eigenstates``.
Default to :obj:`False`.
Returns
-------
spectrum : list[list[float, numpy.ndarray]]
The target number of eigenstates and associated exact energies
in the Hamiltonian spectrum.
The output is a list of exact [energy, eigenstate] sorted in
ascending order w.r.t. the energy values.
Notes
-----
The diagonalization implemented here is straightforward (a.k.a.,
sorting), as it directly follows from the classical nature of the
mapped spin glass QUBO Hamiltonian, which is diagonal in the
computational basis.
"""
# Get the diagonal elements of the Hamiltonian
energies = ham_matrix.diagonal()
hilbert_space_dim = energies.size
# Compute eigenstates and eigenvalues
# Partial spectrum
if not full_spectrum:
if not 0 < number_of_eigenstates <= hilbert_space_dim:
raise ValueError(
"Number of eigenstates must be between 1 and the dimension of"
f" the Hilbert space. Got {number_of_eigenstates}."
)
eigs_idxs = np.argpartition(energies, number_of_eigenstates)[
:number_of_eigenstates
]
eigs_idxs = eigs_idxs[np.argsort(energies[eigs_idxs])]
# Full spectrum
else:
eigs_idxs = np.arange(hilbert_space_dim)
spectrum = [
[energies[ii], np.bincount([ii], minlength=hilbert_space_dim)]
for ii in eigs_idxs
]
# Output
return spectrum
def measure_exact_observables(
observables: dict[str, list[SparseOrDense]],
exact_spectrum: list[list[float, npt.NDArray]],
) -> dict[str, list[float]]:
r"""
Measures observables on the exact eigenstates of a classical
many-body Hamiltonian.
Arguments
---------
observables: dict[str, list[scipy.sparse.csr_matrix | numpy.ndarray]]
Observables to measure.
See :func:`~construct_exact_observables` for more details.
exact_spectrum : list[list[float, numpy.ndarray]]
The list of [energy, eigenvector] obtained via the exact
computation of a classical many-body Hamiltonian.
See :func:`~get_exact_spectrum` for more details.
Returns
-------
expectation_values : dict[str, list[float]]
Quantum average
.. math:
\\langle \\psi \\vert \\mathcal{O} \\vert \\psi \\rangle
of each observable in ``observables`` on the exact computed
spectrum.
Keys of the resulting dictionary correspond to observable keys.
"""
# Loop over requested observables
expectation_values = {
obs_name: [
[(np.conj(eigv[1]).T @ op @ eigv[1]).item() for op in obs_ops]
for eigv in exact_spectrum
]
for obs_name, obs_ops in observables.items()
}
# Output
return expectation_values
# =========================================
# Tensor network ground-state search: utils
# =========================================
def generate_perturbation_transverse_fields(
couplings: dict[str, FloatOrArray], ratio: float | None = None
) -> ArrayOrNone:
"""
Generates a set of local random inhomogeneous transverse
one-body couplings to add a non-commuting (X-Pauli like)
term at each site to the classical MP4EO Hamiltonian.
Arguments
---------
couplings : dict[str, float | numpy.ndarray[float]]
The set of couplings defining the quantum Hamiltonian.
Expected keys:
- **'offset'**: the constant term proportional to the identity.
- **'one-body'**: the set of single-body couplings, one for each
quantum degrees of freedom in the mapped quantum system.
- **'two-body'**: the set of two-body (in general all-to-all)
couplings describing the interactions between pairs of degrees
of freedom.
The method also works in the presence of other terms
beyond pairwise interaction, which may arise in more
sophisticated quantum-inspired mappings.
ratio : float
The ratio of the magnitude of the classical coupling
strengths to the quantum off-diagonal transverse
couplings to be added.
If :obj:`None`, the method returns :obj:`None`.
Returns
-------
transverse_couplings : numpy.ndarray[float] | None
The set of random local transverse couplings to be added
to the classical spin glass Hamiltonian, or :obj:`None`
is ``ratio`` is :obj:`None`.
Details
-------
These perturbations are randomly drawn, with their magnitude
determined by the average value of the Hamiltonian couplings.
These terms could improve the ground-state search by aiding in
escaping local minima during tensor-network optimization.
"""
if ratio is None:
return None
# Validate arguments
if not isinstance(couplings, dict):
raise TypeError(
"The spin glass Hamiltonian couplings must be represented"
" as a dictionary."
)
assert ratio > 0, ValueError("Ratio must be a positive number.")
# Determine the order of magnitude of the perturbation
abs_values = [
np.abs(coupling)
for coupling in couplings.values()
if isinstance(coupling, np.ndarray) and np.issubdtype(coupling.dtype, np.number)
]
coupling_mean_value = np.mean([np.mean(value) for value in abs_values])
order_of_magnitude = coupling_mean_value / ratio
# Generate local random one-body transverse couplings
transverse_couplings = np.random.uniform(
low=-order_of_magnitude,
high=order_of_magnitude,
size=couplings["one-body"].shape,
)
# Output
return transverse_couplings
def construct_spinglass_driver_state(
nsites: int,
bond_dimension: int,
numerical_noise: float,
tn_ansatz: str,
tn_backend: TensorBackend,
) -> MPS | TTN:
r"""
Constructs the product state
.. math::
{\\left\\vert + \\rangle}^{\\otimes n}
to be used as an initial state for the tensor-network
ground-state search. Here, :math:`n` is the number of
physical sites in the tensor network ansatz.
Arguments
---------
nsites : int
Number of physical sites in the system.
bond_dimension : int
Bond dimension of the created product state.
This allows to match the desidered bond dimension
for, e.g., a ground-state search starting from this
initial state.
numerical_noise : float
Numerical noise for creating a product state
with bond dimension greater than 1.
tn_ansatz : str
Ansatz for the representation of the product
state.
tn_backend : qtealeaves.tensors.TensorBackend
The backend for the tensors in the initial driver
product state.
This should match the backend used for the ground-state
search.
See :mod:`qtealeaves` and :mod:`qredtea` for more details.
Returns
-------
init_psi : qtealeaves.emulator.TTN | qtealeaves.emulator.MPS
The created driver product state.
Details
-------
The product state is exactly constructed by hand
but it can be represented as a tensor network state
with a given topology and a given bond dimension,
determined by ``bond_dimension``, by padding properly
the tensors.
"""
# Create the local state |+>
local_dimensions = [2] * nsites
local_plus_states = [
np.full(local_dim, 1 / np.sqrt(local_dim)) for local_dim in local_dimensions
]
# Create the product state depending on the topology
init_psi = construct_tn_state_from_local(
local_states=local_plus_states,
bond_dimension=bond_dimension,
numerical_noise=numerical_noise,
tn_ansatz=tn_ansatz,
tn_backend=tn_backend,
)
# Output
return init_psi
def construct_tn_state_from_local(
local_states: ArrayOrListArray,
bond_dimension: int,
numerical_noise: float,
tn_ansatz: str,
tn_backend: TensorBackend,
) -> MPS | TTN:
r"""
Constructs a tensor network state from the given product state
``local_prod_state`` depending on the selected topology.
Arguments
---------
local_states : numpy.ndarray | list[numpy.ndarray]
The local states to be represented as a tensor-network state.
You can pass either a rank-2 :mod:`numpy` array where each
row represents the local state of each site (homogeneous local
dimensions) or a :class:`list` of rank-1 :mod:`numpy` arrays
whenever each site has a different Hilbert space local dimension.
In this case, the length of the list gives the number of sites
in the system.
bond_dimension : int
Bond dimension of the created product state.
numerical_noise : float
Numerical noise for creating an initial tensor-network
state with bond dimension greater than 1.
If 0, no noise will be insert in the tensors entries
when padding.
tn_ansatz : str
Ansatz for the representation of the product state.
tn_backend : qtealeaves.tensors.TensorBackend
The backend for the tensors in the initial product state.
This should match the backend used for, e.g., a ground-state
search that starts from this initial state.
Returns
-------
qtealeaves.emulator.TTN | qtealeaves.emulator.MPS
The created tensor-network product state.
Notes
-----
``local_states`` is a matrix with the :math:`j`-th row being a
normalized local state of the :math:`j`-th physical site (qudit).
The number of rows is therefore equal to the number of
sites, and number of columns corresponds to the local
dimension of the quantum degrees of freedom (qudits).
"""
conv_params = TNConvergenceParameters(max_bond_dimension=bond_dimension)
# Select the correct tensor-network class
ansatz_cls = TTN if tn_ansatz == "TTN" else MPS
# Output
return ansatz_cls.product_state_from_local_states(
mat=local_states,
padding=[bond_dimension, numerical_noise],
convergence_parameters=conv_params,
tensor_backend=tn_backend,
)
def read_tn_state(
filepath: str, tn_ansatz: str, tn_backend: TensorBackend
) -> MPS | TTN:
"""
Reads and returns the wave function of an Hamiltonian, e.g.,
after the ground-state search, as a tensor-network state,
depending on the chosen tensor-network topology ``tn_ansatz``.
Arguments
---------
filepath : str
The name of the file where the TN state to be read is saved.
The full path to the file is required.
tn_ansatz : str
Ansatz for the representation of the loaded state.
tn_backend : qtealeaves.tensors.TensorBackend
The linear algebra backend to use at the tensor level.
Returns
-------
qtealeaves.emulator.MPS | qtealeaves.emulator.TTN
The loaded tensor network state.
"""
# Select the correct tensor-network class
ansatz_cls = TTN if tn_ansatz == "TTN" else MPS
# Check filepath compatibility with qtealeaves
required_suffix = f".pkl{tn_ansatz.lower()}"
if not filepath.endswith(required_suffix):
filepath += required_suffix
# Output
return ansatz_cls.read(filename=filepath, tensor_backend=tn_backend)