# This code is part of a Qiskit project.
#
# (C) Copyright IBM 2022, 2025.
#
# 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 Variational Quantum Deflation Algorithm for computing higher energy states.
See https://arxiv.org/abs/1805.08138.
"""
from __future__ import annotations
import logging
import warnings
from collections.abc import Callable, Sequence, Iterable
from time import time
from typing import Any, cast
import numpy as np
from qiskit.circuit import QuantumCircuit
from qiskit.primitives import BaseEstimatorV2
from qiskit.quantum_info import SparsePauliOp
from qiskit.quantum_info.operators.base_operator import BaseOperator
from qiskit_algorithms.state_fidelities import BaseStateFidelity
from .eigensolver import Eigensolver, EigensolverResult
from ..custom_types import Transpiler
from ..exceptions import AlgorithmError
from ..list_or_dict import ListOrDict
from ..observables_evaluator import estimate_observables
from ..optimizers import Optimizer, Minimizer, OptimizerResult
from ..utils import validate_bounds, validate_initial_point
# private function as we expect this to be updated in the next release
from ..utils.set_batching import _set_default_batchsize
from ..variational_algorithm import VariationalAlgorithm
logger = logging.getLogger(__name__)
[docs]
class VQD(VariationalAlgorithm, Eigensolver):
r"""The Variational Quantum Deflation algorithm. Implementation using primitives.
`VQD <https://arxiv.org/abs/1805.08138>`__ is a quantum algorithm that uses a
variational technique to find
the k lowest eigenvalues of the Hamiltonian :math:`H` of a given system.
The algorithm computes excited state energies of generalised hamiltonians
by optimizing over a modified cost function where each successive eigenvalue
is calculated iteratively by introducing an overlap term with all
the previously computed eigenstates that must be minimised, thus ensuring
higher energy eigenstates are found.
An instance of VQD requires defining three algorithmic subcomponents:
an integer k denoting the number of eigenstates to calculate, a trial
state (a.k.a. ansatz) which is a :class:`QuantumCircuit`,
and one instance (or list of) classical :mod:`~qiskit_algorithms.optimizers`.
The optimizer varies the circuit parameters
The trial state :math:`|\psi(\vec\theta)\rangle` is varied by the optimizer,
which modifies the set of ansatz parameters :math:`\vec\theta`
such that the expectation value of the operator on the corresponding
state approaches a minimum. The algorithm does this by iteratively refining
each excited state to be orthogonal to all the previous excited states.
An optional array of parameter values, via the *initial_point*, may be provided
as the starting point for the search of the minimum eigenvalue. This feature is
particularly useful when there are reasons to believe that the solution point
is close to a particular point.
The length of the *initial_point* list value must match the number of the parameters
expected by the ansatz. If the *initial_point* is left at the default
of ``None``, then VQD will look to the ansatz for a preferred value, based on its
given initial state. If the ansatz returns ``None``,
then a random point will be generated within the parameter bounds set, as per above.
It is also possible to give a list of initial points, one for every kth eigenvalue.
If the ansatz provides ``None`` as the lower bound, then VQD
will default it to :math:`-2\pi`; similarly, if the ansatz returns ``None``
as the upper bound, the default value will be :math:`2\pi`.
The following attributes can be set via the initializer but can also be read and
updated once the VQD object has been constructed.
Attributes:
estimator (BaseEstimatorV2): The primitive instance used to perform the expectation
estimation as indicated in the VQD paper.
fidelity (BaseStateFidelity): The fidelity class instance used to compute the
overlap estimation as indicated in the VQD paper.
optimizer(Optimizer | Sequence[Optimizer]): A classical optimizer or a list of optimizers,
one for every k-th eigenvalue. Can either be a Qiskit optimizer or a callable
that takes an array as input and returns a Qiskit or SciPy optimization result.
k (int): the number of eigenvalues to return. Returns the lowest k eigenvalues.
betas (list[float]): Beta parameters in the VQD paper.
Should have length k - 1, with k the number of excited states.
These hyper-parameters balance the contribution of each overlap term to the cost
function and have a default value computed as the mean square sum of the
coefficients of the observable.
callback (Callable[[int, np.ndarray, float, dict[str, Any]], None] | None):
A callback that can access the intermediate data
during the optimization. Four parameter values are passed to the callback as
follows during each evaluation by the optimizer: the evaluation count,
the optimizer parameters for the ansatz, the estimated value, the estimation
metadata, and the current step.
convergence_threshold: A threshold under which the algorithm is considered to have
converged. It corresponds to the maximal average fidelity an eigenstate is allowed
to have with the previous eigenstates. If set to None, no check is performed.
"""
def __init__(
self,
estimator: BaseEstimatorV2,
fidelity: BaseStateFidelity,
ansatz: QuantumCircuit,
optimizer: Optimizer | Minimizer | Sequence[Optimizer | Minimizer],
*,
k: int = 2,
betas: np.ndarray | None = None,
initial_point: np.ndarray | list[np.ndarray] | None = None,
callback: Callable[[int, np.ndarray, float, dict[str, Any], int], None] | None = None,
convergence_threshold: float | None = None,
transpiler: Transpiler | None = None,
transpiler_options: dict[str, Any] | None = None,
) -> None:
"""
Args:
estimator: The estimator primitive.
fidelity: The fidelity class using primitives.
ansatz: A parameterized circuit used as ansatz for the wave function.
optimizer: A classical optimizer or a list of optimizers, one for every k-th eigenvalue.
Can either be a Qiskit optimizer or a callable
that takes an array as input and returns a Qiskit or SciPy optimization result.
k: The number of eigenvalues to return. Returns the lowest k eigenvalues.
betas: Beta parameters in the VQD paper.
Should have length k - 1, with k the number of excited states.
These hyperparameters balance the contribution of each overlap term to the cost
function and have a default value computed as the mean square sum of the
coefficients of the observable.
initial_point: An optional initial point (i.e. initial parameter values)
or a list of initial points (one for every k-th eigenvalue)
for the optimizer.
If ``None`` then VQD will look to the ansatz for a preferred
point and if not will simply compute a random one.
callback: A callback that can access the intermediate data
during the optimization. Five parameter values are passed to the callback as
follows during each evaluation by the optimizer: the evaluation count,
the optimizer parameters for the ansatz, the estimated value,
the estimation metadata, and the current step.
convergence_threshold: A threshold under which the algorithm is considered to have
converged. It corresponds to the maximal average fidelity an eigenstate is allowed
to have with the previous eigenstates. If set to None, no check is performed.
transpiler: An optional object with a `run` method allowing to transpile the circuits
that are run when using this algorithm. If set to `None`, these won't be
transpiled.
transpiler_options: A dictionary of options to be passed to the transpiler's `run`
method as keyword arguments.
"""
super().__init__()
self.estimator = estimator
self.fidelity = fidelity
self._ansatz = ansatz
self.optimizer = optimizer
self.k = k
self.betas = betas
# this has to go via getters and setters due to the VariationalAlgorithm interface
self.initial_point = initial_point
self.callback = callback
self.convergence_threshold = convergence_threshold
self._transpiler = transpiler
self._transpiler_options = transpiler_options if transpiler_options is not None else {}
if self._transpiler is not None:
self.ansatz = ansatz
self._eval_count = 0
@property
def initial_point(self) -> np.ndarray | list[np.ndarray] | None:
"""Returns initial point."""
return self._initial_point
@initial_point.setter
def initial_point(self, initial_point: np.ndarray | list[np.ndarray] | None):
"""Sets initial point"""
self._initial_point = initial_point
@property
def ansatz(self) -> QuantumCircuit:
"""
A parameterized circuit used as ansatz for the wave function. If a transpiler has been
provided, the ansatz will be automatically transpiled upon being set.
"""
return self._ansatz
@ansatz.setter
def ansatz(self, value: QuantumCircuit | None) -> None:
if self._transpiler is not None:
self._ansatz = self._transpiler.run(value, **self._transpiler_options)
else:
self._ansatz = value
def _check_operator_ansatz(self, operator: BaseOperator):
"""Check that the number of qubits of operator and ansatz match."""
if operator is not None and self.ansatz is not None:
if operator.num_qubits != self.ansatz.num_qubits:
# try to set the number of qubits on the ansatz, if possible
try:
self.ansatz.num_qubits = operator.num_qubits
warnings.warn(
"Previously, it was possible to pass to VQD a BlueprintCircuit as an "
"ansatz without its number of qubits being set, the algorithm taking care "
"of setting it. Since BlueprintCircuits are now deprecated, and those "
"being the only ones that can have their number of qubits set after their "
"initialization, this behavior is now also deprecated, and won't be "
"supported once the oldest supported Qiskit version is 3.0. As such, users "
"that made use of this feature would now need to ensure that the ansatz "
"they pass to these algorithms have their number of qubits set and matching "
"with that of the operator they wish to run the algorithm on.",
category=DeprecationWarning,
)
except AttributeError as exc:
raise AlgorithmError(
"The number of qubits of the ansatz does not match the "
"operator, and the ansatz does not allow setting the "
"number of qubits using `num_qubits`."
) from exc
[docs]
@classmethod
def supports_aux_operators(cls) -> bool:
return True
[docs]
def compute_eigenvalues(
self,
operator: BaseOperator,
aux_operators: ListOrDict[BaseOperator] | None = None,
) -> VQDResult:
super().compute_eigenvalues(operator, aux_operators)
if self.ansatz.layout is not None:
operator = operator.apply_layout(self.ansatz.layout)
# this sets the size of the ansatz, so it must be called before the initial point
# validation
self._check_operator_ansatz(operator)
bounds = validate_bounds(self.ansatz)
# We need to handle the array entries being zero or Optional i.e. having value None
if aux_operators:
if self.ansatz.layout is not None:
# len(self.ansatz.layout.final_index_layout()) is the original number of qubits in the
# ansatz, before transpilation
zero_op = SparsePauliOp.from_list(
[("I" * len(self.ansatz.layout.final_index_layout()), 0)]
)
else:
zero_op = SparsePauliOp.from_list([("I" * self.ansatz.num_qubits, 0)])
# Convert the None and zero values when aux_operators is a list.
# Drop None and convert zero values when aux_operators is a dict.
key_op_iterator: Iterable[tuple[str | int, BaseOperator]]
if isinstance(aux_operators, list):
aux_operators = [op if op is not None else zero_op for op in aux_operators]
key_op_iterator = enumerate(aux_operators)
converted: ListOrDict[BaseOperator] = [zero_op] * len(aux_operators)
else:
key_op_iterator = aux_operators.items()
converted = {}
for key, op in key_op_iterator:
if op is not None:
converted[key] = zero_op if op == 0 else op # type: ignore[index]
if self.ansatz.layout is not None:
converted[key] = converted[key].apply_layout(self.ansatz.layout)
aux_operators = converted
else:
aux_operators = None
betas = self.betas
if self.betas is None:
try:
upper_bound = sum(np.abs(operator.coeffs))
except Exception as exc:
raise NotImplementedError(
r"Beta autoevaluation is not supported for operators"
f"of type {type(operator)}."
) from exc
betas = np.asarray([upper_bound * 10] * self.k)
logger.info("beta autoevaluated to %s", betas[0])
result = self._build_vqd_result()
if aux_operators is not None:
aux_values = []
# We keep a list of the bound circuits with optimal parameters, to avoid re-binding
# the same parameters to the ansatz if we do multiple steps
prev_states = []
# These two variables are defined inside if statements and static analysis, e.g. lint can
# see this as a potential error of them not being defined before use. Following the logic
# they do end up being defined before use so the setting of these here, these values would
# not be used in practice.
initial_point = np.asarray([])
initial_points = np.asarray([])
num_initial_points = 0
if self.initial_point is not None:
initial_points = np.reshape(self.initial_point, (-1, self.ansatz.num_parameters))
num_initial_points = len(initial_points)
# 0 just means the initial point is ``None`` and ``validate_initial_point``
# will select a random point
if num_initial_points <= 1:
initial_point = validate_initial_point(
self.initial_point, self.ansatz # type: ignore[arg-type]
)
current_optimal_point: dict[str, Any] = {"optimal_value": float("inf")}
for step in range(1, self.k + 1):
current_optimal_point["optimal_value"] = float("inf")
if num_initial_points > 1:
initial_point = validate_initial_point(initial_points[step - 1], self.ansatz)
if step > 1:
prev_states.append(self.ansatz.assign_parameters(current_optimal_point["x"]))
self._eval_count = 0
energy_evaluation = self._get_evaluate_energy(
step,
operator,
betas,
prev_states=prev_states,
current_optimal_point=current_optimal_point,
)
start_time = time()
# TODO: add gradient support after FidelityGradients are implemented
if isinstance(self.optimizer, Sequence):
optimizer = self.optimizer[step - 1]
else:
optimizer = self.optimizer # fall back to single optimizer if not list
if callable(optimizer):
opt_result = optimizer( # pylint: disable=not-callable
fun=energy_evaluation, # type: ignore[arg-type]
x0=initial_point,
jac=None,
bounds=bounds,
)
else:
# we always want to submit as many estimations per job as possible for minimal
# overhead on the hardware
was_updated = _set_default_batchsize(optimizer)
opt_result = optimizer.minimize(
fun=energy_evaluation, x0=initial_point, bounds=bounds # type: ignore[arg-type]
)
# reset to original value
if was_updated:
optimizer.set_max_evals_grouped(None)
eval_time = time() - start_time
self._update_vqd_result(
result, opt_result, eval_time, self.ansatz.copy(), current_optimal_point
)
if aux_operators is not None:
aux_value = estimate_observables(
self.estimator, self.ansatz, aux_operators, current_optimal_point["x"]
)
aux_values.append(aux_value)
if step == 1:
logger.info(
"Ground state optimization complete in %s seconds.\n"
"Found opt_params %s in %s evals",
eval_time,
result.optimal_points,
self._eval_count,
)
else:
average_fidelity = current_optimal_point["total_fidelity"][0] / (step - 1)
if (
self.convergence_threshold is not None
and average_fidelity > self.convergence_threshold
):
last_digit = step % 10
if last_digit == 1 and step % 100 != 11:
suffix = "st"
elif last_digit == 2:
suffix = "nd"
elif last_digit == 3:
suffix = "rd"
else:
suffix = "th"
raise AlgorithmError(
f"Convergence threshold is set to {self.convergence_threshold} but an "
f"average (weighted by the betas) fidelity of {average_fidelity:.5f} with "
f"the previous eigenstates has been observed during the evaluation of the "
f"{step}{suffix} lowest eigenvalue."
)
logger.info(
(
"%s excited state optimization complete in %s s.\n"
"Found opt_params %s in %s evals"
),
str(step - 1),
eval_time,
result.optimal_points,
self._eval_count,
)
# To match the signature of EigensolverResult
result.eigenvalues = np.array(result.eigenvalues)
if aux_operators is not None:
result.aux_operators_evaluated = aux_values
return result
def _get_evaluate_energy( # pylint: disable=too-many-positional-arguments
self,
step: int,
operator: BaseOperator,
betas: np.ndarray,
current_optimal_point: dict["str", Any],
prev_states: list[QuantumCircuit] | None = None,
) -> Callable[[np.ndarray], float | np.ndarray]:
"""Returns a function handle to evaluate the ansatz's energy for any given parameters.
This is the objective function to be passed to the optimizer that is used for evaluation.
Args:
step: level of energy being calculated. 1 for ground, 2 for first excited state...
operator: The operator whose energy to evaluate.
betas: Beta parameters in the VQD paper.
prev_states: List of optimal circuits from previous rounds of optimization.
current_optimal_point: A dict to keep track of the current optimal point, which is used
to check the algorithm's convergence.
Returns:
A callable that computes and returns the energy of the hamiltonian
of each parameter.
Raises:
AlgorithmError: If the circuit is not parameterized (i.e. has 0 free parameters).
AlgorithmError: If operator was not provided.
RuntimeError: If the previous states array is of the wrong size.
"""
num_parameters = self.ansatz.num_parameters
if num_parameters == 0:
raise AlgorithmError("The ansatz must be parameterized, but has no free parameters.")
if step > 1 and (len(prev_states) + 1) != step:
raise RuntimeError(
f"Passed previous states of the wrong size."
f"Passed array has length {str(len(prev_states))}"
)
def evaluate_energy(parameters: np.ndarray) -> float | np.ndarray:
# handle broadcasting: ensure parameters is of shape [array, array, ...]
if len(parameters.shape) == 1:
parameters = np.reshape(parameters, (-1, num_parameters))
batch_size = len(parameters)
estimator_job = self.estimator.run([(self.ansatz, operator, parameters)])
total_cost = np.zeros(batch_size)
if step > 1:
# compute overlap cost
batched_prev_states = [state for state in prev_states for _ in range(batch_size)]
fidelity_job = self.fidelity.run(
batch_size * [self.ansatz] * (step - 1),
batched_prev_states,
np.tile(parameters, (step - 1, 1)), # type: ignore[arg-type]
)
costs = fidelity_job.result().fidelities
costs = np.reshape(costs, (step - 1, -1))
for state, cost in enumerate(costs):
total_cost += np.real(betas[state] * cost)
try:
estimator_result = estimator_job.result()[0]
except Exception as exc:
raise AlgorithmError("The primitive job to evaluate the energy failed!") from exc
values = estimator_result.data.evs + total_cost
if self.callback is not None:
for params, value in zip(parameters, values):
self._eval_count += 1
self.callback(self._eval_count, params, value, estimator_result.metadata, step)
else:
self._eval_count += len(values)
for param, value in zip(parameters, values):
if value < current_optimal_point["optimal_value"]:
current_optimal_point["optimal_value"] = value
current_optimal_point["x"] = param
if step > 1:
current_optimal_point["total_fidelity"] = total_cost
current_optimal_point["eigenvalue"] = (value - total_cost)[0]
else:
current_optimal_point["eigenvalue"] = value
return values if len(values) > 1 else values[0]
return evaluate_energy
@staticmethod
def _build_vqd_result() -> VQDResult:
result = VQDResult()
result.optimal_points = np.array([])
result.optimal_parameters = []
result.optimal_values = np.array([])
result.cost_function_evals = np.array([], dtype=int)
result.optimizer_times = np.array([])
result.eigenvalues = [] # type: ignore[assignment]
result.optimizer_results = []
result.optimal_circuits = []
return result
@staticmethod
def _update_vqd_result(
result: VQDResult, opt_result: OptimizerResult, eval_time, ansatz, optimal_point
) -> VQDResult:
result.optimal_points = (
np.concatenate([result.optimal_points, [optimal_point["x"]]])
if len(result.optimal_points) > 0
else np.array([optimal_point["x"]])
)
result.optimal_parameters.append(
dict(zip(ansatz.parameters, cast(np.ndarray, optimal_point["x"])))
)
result.optimal_values = np.concatenate(
[result.optimal_values, [optimal_point["optimal_value"]]]
)
result.cost_function_evals = np.concatenate([result.cost_function_evals, [opt_result.nfev]])
result.optimizer_times = np.concatenate([result.optimizer_times, [eval_time]])
result.eigenvalues.append(optimal_point["eigenvalue"] + 0j) # type: ignore[attr-defined]
result.optimizer_results.append(opt_result)
result.optimal_circuits.append(ansatz)
return result
[docs]
class VQDResult(EigensolverResult):
"""VQD Result."""
def __init__(self) -> None:
super().__init__()
self._cost_function_evals: np.ndarray | None = None
self._optimizer_times: np.ndarray | None = None
self._optimal_values: np.ndarray | None = None
self._optimal_points: np.ndarray | None = None
self._optimal_parameters: list[dict] | None = None
self._optimizer_results: list[OptimizerResult] | None = None
self._optimal_circuits: list[QuantumCircuit] | None = None
@property
def cost_function_evals(self) -> np.ndarray | None:
"""Returns number of cost optimizer evaluations"""
return self._cost_function_evals
@cost_function_evals.setter
def cost_function_evals(self, value: np.ndarray) -> None:
"""Sets number of cost function evaluations"""
self._cost_function_evals = value
@property
def optimizer_times(self) -> np.ndarray | None:
"""Returns time taken for optimization for each step"""
return self._optimizer_times
@optimizer_times.setter
def optimizer_times(self, value: np.ndarray) -> None:
"""Sets time taken for optimization for each step"""
self._optimizer_times = value
@property
def optimal_values(self) -> np.ndarray | None:
"""Returns optimal value for each step"""
return self._optimal_values
@optimal_values.setter
def optimal_values(self, value: np.ndarray) -> None:
"""Sets optimal values"""
self._optimal_values = value
@property
def optimal_points(self) -> np.ndarray | None:
"""Returns optimal point for each step"""
return self._optimal_points
@optimal_points.setter
def optimal_points(self, value: np.ndarray) -> None:
"""Sets optimal points"""
self._optimal_points = value
@property
def optimal_parameters(self) -> list[dict] | None:
"""Returns the optimal parameters for each step"""
return self._optimal_parameters
@optimal_parameters.setter
def optimal_parameters(self, value: list[dict]) -> None:
"""Sets optimal parameters"""
self._optimal_parameters = value
@property
def optimizer_results(self) -> list[OptimizerResult] | None:
"""Returns the optimizer results for each step"""
return self._optimizer_results
@optimizer_results.setter
def optimizer_results(self, value: list[OptimizerResult]) -> None:
"""Sets optimizer results"""
self._optimizer_results = value
@property
def optimal_circuits(self) -> list[QuantumCircuit] | None:
"""The optimal circuits. Along with the optimal parameters,
these can be used to retrieve the different eigenstates."""
return self._optimal_circuits
@optimal_circuits.setter
def optimal_circuits(self, optimal_circuits: list[QuantumCircuit]) -> None:
self._optimal_circuits = optimal_circuits