Source code for qiskit_algorithms.eigensolvers.numpy_eigensolver

# 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
# 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 NumPy eigensolver algorithm."""

from __future__ import annotations

from import Iterable
from typing import Callable, Union, Tuple, Dict, List, Optional, cast
import logging
import numpy as np
from scipy import sparse as scisparse

from qiskit.quantum_info import SparsePauliOp, Statevector
from qiskit.quantum_info.operators.base_operator import BaseOperator

from qiskit_algorithms.utils.validation import validate_min
from .eigensolver import Eigensolver, EigensolverResult
from ..exceptions import AlgorithmError
from ..list_or_dict import ListOrDict

logger = logging.getLogger(__name__)

FilterType = Callable[
    [Union[List, np.ndarray], float, Optional[ListOrDict[Tuple[float, Dict[str, float]]]]], bool

[docs]class NumPyEigensolver(Eigensolver): r""" The NumPy eigensolver algorithm. The NumPy Eigensolver computes up to the first :math:`k` eigenvalues of a complex-valued square matrix of dimension :math:`n \times n`, with :math:`k \leq n`. Note: Operators are automatically converted to SciPy's ``spmatrix`` as needed and this conversion can be costly in terms of memory and performance as the operator size, mostly in terms of number of qubits it represents, gets larger. """ def __init__( self, k: int = 1, filter_criterion: FilterType | None = None, ) -> None: """ Args: k: Number of eigenvalues are to be computed, with a minimum value of 1. filter_criterion: Callable that allows to filter eigenvalues/eigenstates. Only feasible eigenstates are returned in the results. The callable has the signature ``filter(eigenstate, eigenvalue, aux_values)`` and must return a boolean to indicate whether to keep this value in the final returned result or not. If the number of elements that satisfies the criterion is smaller than ``k``, then the returned list will have fewer elements and can even be empty. """ validate_min("k", k, 1) super().__init__() self._in_k = k self._k = k # pylint: disable=invalid-name self._filter_criterion = filter_criterion @property def k(self) -> int: """Return k (number of eigenvalues requested).""" return self._in_k @k.setter def k(self, k: int) -> None: """Set k (number of eigenvalues requested).""" validate_min("k", k, 1) self._in_k = k self._k = k @property def filter_criterion( self, ) -> FilterType | None: """Return the filter criterion if set.""" return self._filter_criterion @filter_criterion.setter def filter_criterion(self, filter_criterion: FilterType | None) -> None: """Set the filter criterion.""" self._filter_criterion = filter_criterion
[docs] @classmethod def supports_aux_operators(cls) -> bool: return True
def _check_set_k(self, operator: BaseOperator) -> None: if operator is not None: if self._in_k > 2**operator.num_qubits: self._k = 2**operator.num_qubits logger.debug( "WARNING: Asked for %s eigenvalues but max possible is %s.", self._in_k, self._k ) else: self._k = self._in_k def _solve(self, operator: BaseOperator) -> tuple[np.ndarray, np.ndarray]: try: op_matrix = operator.to_matrix(sparse=True) except TypeError: logger.debug( "WARNING: operator of type `%s` does not support sparse matrices. " "Trying dense computation", type(operator), ) try: op_matrix = operator.to_matrix() except AttributeError as ex: raise AlgorithmError(f"Unsupported operator type `{type(operator)}`.") from ex if isinstance(op_matrix, scisparse.csr_matrix): # If matrix is diagonal, the elements on the diagonal are the eigenvalues. Solve by sorting. if scisparse.csr_matrix(op_matrix.diagonal()).nnz == op_matrix.nnz: diag = op_matrix.diagonal() indices = np.argsort(diag)[: self._k] eigval = diag[indices] eigvec = np.zeros((op_matrix.shape[0], self._k)) for i, idx in enumerate(indices): eigvec[idx, i] = 1.0 else: if self._k >= 2**operator.num_qubits - 1: logger.debug( "SciPy doesn't support to get all eigenvalues, using NumPy instead." ) eigval, eigvec = self._solve_dense(operator.to_matrix()) else: eigval, eigvec = self._solve_sparse(op_matrix, self._k) else: # Sparse SciPy matrix not supported, use dense NumPy computation. eigval, eigvec = self._solve_dense(operator.to_matrix()) indices = np.argsort(eigval)[: self._k] eigval = eigval[indices] eigvec = eigvec[:, indices] return eigval, eigvec.T @staticmethod def _solve_sparse(op_matrix: scisparse.csr_matrix, k: int) -> tuple[np.ndarray, np.ndarray]: if (op_matrix != op_matrix.H).nnz == 0: # Operator is Hermitian return scisparse.linalg.eigsh(op_matrix, k=k, which="SA") else: return scisparse.linalg.eigs(op_matrix, k=k, which="SR") @staticmethod def _solve_dense(op_matrix: np.ndarray) -> tuple[np.ndarray, np.ndarray]: if op_matrix.all() == op_matrix.conj().T.all(): # Operator is Hermitian return cast(Tuple[np.ndarray, np.ndarray], np.linalg.eigh(op_matrix)) else: return cast(Tuple[np.ndarray, np.ndarray], np.linalg.eig(op_matrix)) @staticmethod def _eval_aux_operators( aux_operators: ListOrDict[BaseOperator], wavefn: np.ndarray, threshold: float = 1e-12, ) -> ListOrDict[tuple[float, dict[str, float]]]: values: ListOrDict[tuple[float, dict[str, float]]] # As a list, aux_operators can contain None operators for which None values are returned. # As a dict, the None operators in aux_operators have been dropped in compute_eigenvalues. key_op_iterator: Iterable[tuple[str | int, BaseOperator]] if isinstance(aux_operators, list): values = [None] * len(aux_operators) key_op_iterator = enumerate(aux_operators) else: values = {} key_op_iterator = aux_operators.items() for key, operator in key_op_iterator: if operator is None: continue if operator.num_qubits is None or operator.num_qubits < 1: "The number of qubits of the %s operator must be greater than zero.", key ) continue op_matrix = None try: op_matrix = operator.to_matrix(sparse=True) except TypeError: logger.debug( "WARNING: operator of type `%s` does not support sparse matrices. " "Trying dense computation", type(operator), ) try: op_matrix = operator.to_matrix() except AttributeError as ex: raise AlgorithmError(f"Unsupported operator type {type(operator)}.") from ex if isinstance(op_matrix, scisparse.csr_matrix): value = elif isinstance(op_matrix, np.ndarray): value = Statevector(wavefn).expectation_value(operator) else: value = 0.0 value = value if np.abs(value) > threshold else 0.0 # The value gets wrapped into a tuple: (mean, metadata). # The metadata includes variance (and, for other eigensolvers, shots). # Since this is an exact computation, there are no shots # and the variance is known to be zero. values[key] = (value, {"variance": 0.0}) # type: ignore[index] return values
[docs] def compute_eigenvalues( self, operator: BaseOperator, aux_operators: ListOrDict[BaseOperator] | None = None, ) -> NumPyEigensolverResult: super().compute_eigenvalues(operator, aux_operators) if operator.num_qubits is None or operator.num_qubits < 1: raise AlgorithmError("The number of qubits of the operator must be greater than zero.") self._check_set_k(operator) zero_op = SparsePauliOp(["I" * operator.num_qubits], coeffs=[0.0]) if isinstance(aux_operators, list) and len(aux_operators) > 0: # For some reason Chemistry passes aux_ops with 0 qubits and paulis sometimes. aux_operators = [zero_op if op == 0 else op for op in aux_operators] elif isinstance(aux_operators, dict) and len(aux_operators) > 0: aux_operators = { key: zero_op if op == 0 else op # Convert zero values to zero operators for key, op in aux_operators.items() if op is not None # Discard None values } else: aux_operators = None k_orig = self._k if self._filter_criterion: # need to consider all elements if a filter is set self._k = 2**operator.num_qubits eigvals, eigvecs = self._solve(operator) # compute energies before filtering, as this also evaluates the aux operators if aux_operators is not None: aux_op_vals = [ self._eval_aux_operators(aux_operators, eigvecs[i]) for i in range(self._k) ] else: aux_op_vals = None # if a filter is set, loop over the given values and only keep if self._filter_criterion: filt_eigvals = [] filt_eigvecs = [] filt_aux_op_vals = [] count = 0 for i, (eigval, eigvec) in enumerate(zip(eigvals, eigvecs)): if aux_op_vals is not None: aux_op_val = aux_op_vals[i] else: aux_op_val = None if self._filter_criterion(eigvec, eigval, aux_op_val): count += 1 filt_eigvecs.append(eigvec) filt_eigvals.append(eigval) if aux_op_vals is not None: filt_aux_op_vals.append(aux_op_val) if count == k_orig: break eigvals = np.array(filt_eigvals) eigvecs = np.array(filt_eigvecs) aux_op_vals = filt_aux_op_vals self._k = k_orig result = NumPyEigensolverResult() result.eigenvalues = eigvals result.eigenstates = [Statevector(vec) for vec in eigvecs] result.aux_operators_evaluated = aux_op_vals logger.debug("NumpyEigensolverResult:\n%s", result) return result
[docs]class NumPyEigensolverResult(EigensolverResult): """NumPy eigensolver result.""" def __init__(self) -> None: super().__init__() self._eigenstates: list[Statevector] | None = None @property def eigenstates(self) -> list[Statevector] | None: """Return eigenstates.""" return self._eigenstates @eigenstates.setter def eigenstates(self, value: list[Statevector]) -> None: """Set eigenstates.""" self._eigenstates = value