From 21f78dc5ee236d51697e65db916ccd04f4dfc3bc Mon Sep 17 00:00:00 2001 From: Julien Gacon Date: Wed, 13 Jan 2021 23:06:47 +0100 Subject: [PATCH] Amplitude Estimation algorithms (#5517) * ae codes * add tests * update docstrings in result objects * add reno * fix old aqua dependencies * fix cyclic imports * try fixing windows tests no. 1 * fix test * try fixing windows attempt #2 * stefan's first comments * accumulate -> evaluate in AE * prob to meas 1 -> prob good state * stefans comments pt. 2 * fix lint & indent * add faster ae + bernoulli tests * fix > 1 qubit case * add comments * add FAE specific tests * raise error for rescale on custom q operator * play with the oracles * fix classicalregister > 1 * use correct reflection and rm X gates around oracle * fix construction of rescaled est. problem * fix phase in grover op * adjust result objects * adjust fae result object * fix minus sign and lint * rm trailing whitespaces * fix tests * fix lint * rm utf8 header * move rescale to problem, add evalschedule to MLAE * code review suggestions * state and obj.qubits not optional in estimation problem * add QAE algos to documentation * allow setting post processing and is good state to None to use default * rm quantum instance logic from interface * fix typehints for is good state * add Faster AE to reno * reset shots to user defined shots in FAE * fix doc build (hopefully) * add result classes to docs * make grover op optional in setter * rm reno -- there will be a global one Co-authored-by: Matthew Treinish Co-authored-by: mergify[bot] <37929162+mergify[bot]@users.noreply.github.com> --- qiskit/algorithms/__init__.py | 38 ++ .../amplitude_estimators/__init__.py | 34 + qiskit/algorithms/amplitude_estimators/ae.py | 585 ++++++++++++++++++ .../amplitude_estimators/ae_utils.py | 249 ++++++++ .../amplitude_estimator.py | 129 ++++ .../estimation_problem.py | 256 ++++++++ qiskit/algorithms/amplitude_estimators/fae.py | 304 +++++++++ qiskit/algorithms/amplitude_estimators/iae.py | 563 +++++++++++++++++ .../algorithms/amplitude_estimators/mlae.py | 564 +++++++++++++++++ qiskit/circuit/library/grover_operator.py | 3 + .../algorithms/test_amplitude_estimators.py | 498 +++++++++++++++ 11 files changed, 3223 insertions(+) create mode 100644 qiskit/algorithms/amplitude_estimators/__init__.py create mode 100644 qiskit/algorithms/amplitude_estimators/ae.py create mode 100644 qiskit/algorithms/amplitude_estimators/ae_utils.py create mode 100644 qiskit/algorithms/amplitude_estimators/amplitude_estimator.py create mode 100644 qiskit/algorithms/amplitude_estimators/estimation_problem.py create mode 100644 qiskit/algorithms/amplitude_estimators/fae.py create mode 100644 qiskit/algorithms/amplitude_estimators/iae.py create mode 100644 qiskit/algorithms/amplitude_estimators/mlae.py create mode 100644 test/python/algorithms/test_amplitude_estimators.py diff --git a/qiskit/algorithms/__init__.py b/qiskit/algorithms/__init__.py index d0efcbe5b005..a45ebe40bfed 100644 --- a/qiskit/algorithms/__init__.py +++ b/qiskit/algorithms/__init__.py @@ -48,6 +48,25 @@ Grover GroverResult +Amplitude Estimators +++++++++++++++++++++ + +.. autosummary:: + :toctree: ../stubs/ + :nosignatures: + + AmplitudeEstimator + AmplitudeEstimatorResult + AmplitudeEstimation + AmplitudeEstimationResult + EstimationProblem + FasterAmplitudeEstimation + FasterAmplitudeEstimationResult + IterativeAmplitudeEstimation + IterativeAmplitudeEstimationResult + MaximumLikelihoodAmplitudeEstimation + MaximumLikelihoodAmplitudeEstimationResult + Eigensolvers ++++++++++++ Algorithms to find eigenvalues of an operator. For chemistry these can be used to find excited @@ -109,6 +128,14 @@ from .algorithm_result import AlgorithmResult from .variational_algorithm import VariationalAlgorithm, VariationalResult from .amplitude_amplifiers import Grover, GroverResult +from .amplitude_estimators import ( + AmplitudeEstimator, AmplitudeEstimatorResult, + AmplitudeEstimation, AmplitudeEstimationResult, + FasterAmplitudeEstimation, FasterAmplitudeEstimationResult, + IterativeAmplitudeEstimation, IterativeAmplitudeEstimationResult, + MaximumLikelihoodAmplitudeEstimation, MaximumLikelihoodAmplitudeEstimationResult, + EstimationProblem +) from .eigen_solvers import NumPyEigensolver, Eigensolver, EigensolverResult from .factorizers import Shor, ShorResult from .minimum_eigen_solvers import (VQE, VQEResult, QAOA, @@ -122,6 +149,17 @@ 'VariationalResult', 'Grover', 'GroverResult', + 'AmplitudeEstimator', + 'AmplitudeEstimatorResult', + 'AmplitudeEstimation', + 'AmplitudeEstimationResult', + 'FasterAmplitudeEstimation', + 'FasterAmplitudeEstimationResult', + 'IterativeAmplitudeEstimation', + 'IterativeAmplitudeEstimationResult', + 'MaximumLikelihoodAmplitudeEstimation', + 'MaximumLikelihoodAmplitudeEstimationResult', + 'EstimationProblem', 'NumPyEigensolver', 'Eigensolver', 'EigensolverResult', diff --git a/qiskit/algorithms/amplitude_estimators/__init__.py b/qiskit/algorithms/amplitude_estimators/__init__.py new file mode 100644 index 000000000000..a9274e707d82 --- /dev/null +++ b/qiskit/algorithms/amplitude_estimators/__init__.py @@ -0,0 +1,34 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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 Amplitude Estimators package.""" + +from .amplitude_estimator import AmplitudeEstimator, AmplitudeEstimatorResult +from .ae import AmplitudeEstimation, AmplitudeEstimationResult +from .fae import FasterAmplitudeEstimation, FasterAmplitudeEstimationResult +from .iae import IterativeAmplitudeEstimation, IterativeAmplitudeEstimationResult +from .mlae import MaximumLikelihoodAmplitudeEstimation, MaximumLikelihoodAmplitudeEstimationResult +from .estimation_problem import EstimationProblem + +__all__ = [ + 'AmplitudeEstimator', + 'AmplitudeEstimatorResult', + 'AmplitudeEstimation', + 'AmplitudeEstimationResult', + 'FasterAmplitudeEstimation', + 'FasterAmplitudeEstimationResult', + 'IterativeAmplitudeEstimation', + 'IterativeAmplitudeEstimationResult', + 'MaximumLikelihoodAmplitudeEstimation', + 'MaximumLikelihoodAmplitudeEstimationResult', + 'EstimationProblem', +] diff --git a/qiskit/algorithms/amplitude_estimators/ae.py b/qiskit/algorithms/amplitude_estimators/ae.py new file mode 100644 index 000000000000..641c8c5751dd --- /dev/null +++ b/qiskit/algorithms/amplitude_estimators/ae.py @@ -0,0 +1,585 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2018, 2020. +# +# 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 Quantum Phase Estimation-based Amplitude Estimation algorithm.""" + +from typing import Optional, Union, List, Tuple, Dict +from collections import OrderedDict +import numpy as np +from scipy.stats import chi2, norm +from scipy.optimize import bisect + +from qiskit import QuantumCircuit, ClassicalRegister +from qiskit.providers import BaseBackend, Backend +from qiskit.utils import QuantumInstance +from .amplitude_estimator import AmplitudeEstimator, AmplitudeEstimatorResult +from .ae_utils import pdf_a, derivative_log_pdf_a, bisect_max +from .estimation_problem import EstimationProblem + + +class AmplitudeEstimation(AmplitudeEstimator): + r"""The Quantum Phase Estimation-based Amplitude Estimation algorithm. + + This class implements the original Quantum Amplitude Estimation (QAE) algorithm, introduced by + [1]. This canonical version uses quantum phase estimation along with a set of :math:`m` + additional evaluation qubits to find an estimate :math:`\tilde{a}`, that is restricted to the + grid + + .. math:: + + \tilde{a} \in \{\sin^2(\pi y / 2^m) : y = 0, ..., 2^{m-1}\} + + More evaluation qubits produce a finer sampling grid, therefore the accuracy of the algorithm + increases with :math:`m`. + + Using a maximum likelihood post processing, this grid constraint can be circumvented. + This improved estimator is implemented as well, see [2] Appendix A for more detail. + + References: + [1]: Brassard, G., Hoyer, P., Mosca, M., & Tapp, A. (2000). + Quantum Amplitude Amplification and Estimation. + `arXiv:quant-ph/0005055 `_. + [2]: Grinko, D., Gacon, J., Zoufal, C., & Woerner, S. (2019). + Iterative Quantum Amplitude Estimation. + `arXiv:1912.05559 `_. + """ + + def __init__(self, num_eval_qubits: int, + phase_estimation_circuit: Optional[QuantumCircuit] = None, + iqft: Optional[QuantumCircuit] = None, + quantum_instance: Optional[ + Union[QuantumInstance, BaseBackend, Backend]] = None, + ) -> None: + r""" + Args: + num_eval_qubits: The number of evaluation qubits. + phase_estimation_circuit: The phase estimation circuit used to run the algorithm. + Defaults to the standard phase estimation circuit from the circuit library, + `qiskit.circuit.library.PhaseEstimation` when None. + iqft: The inverse quantum Fourier transform component, defaults to using a standard + implementation from `qiskit.circuit.library.QFT` when None. + quantum_instance: The backend (or `QuantumInstance`) to execute the circuits on. + + Raises: + ValueError: If the number of evaluation qubits is smaller than 1. + """ + if num_eval_qubits < 1: + raise ValueError('The number of evaluation qubits must at least be 1.') + + super().__init__() + + # set quantum instance + self.quantum_instance = quantum_instance + + # get parameters + self._m = num_eval_qubits # pylint: disable=invalid-name + self._M = 2 ** num_eval_qubits # pylint: disable=invalid-name + + self._iqft = iqft + self._pec = phase_estimation_circuit + + @property + def quantum_instance(self) -> Optional[QuantumInstance]: + """Get the quantum instance. + + Returns: + The quantum instance used to run this algorithm. + """ + return self._quantum_instance + + @quantum_instance.setter + def quantum_instance(self, quantum_instance: Union[QuantumInstance, + BaseBackend, Backend]) -> None: + """Set quantum instance. + + Args: + quantum_instance: The quantum instance used to run this algorithm. + """ + if isinstance(quantum_instance, (BaseBackend, Backend)): + quantum_instance = QuantumInstance(quantum_instance) + self._quantum_instance = quantum_instance + + def construct_circuit(self, + estimation_problem: EstimationProblem, + measurement: bool = False) -> QuantumCircuit: + """Construct the Amplitude Estimation quantum circuit. + + Args: + estimation_problem: The estimation problem for which to construct the QAE circuit. + measurement: Boolean flag to indicate if measurements should be included in the circuit. + + Returns: + The QuantumCircuit object for the constructed circuit. + """ + # use custom Phase Estimation circuit if provided + if self._pec is not None: + pec = self._pec + + # otherwise use the circuit library -- note that this does not include the A operator + else: + from qiskit.circuit.library import PhaseEstimation + pec = PhaseEstimation(self._m, estimation_problem.grover_operator, iqft=self._iqft) + + # combine the Phase Estimation circuit with the A operator + circuit = QuantumCircuit(*pec.qregs) + circuit.compose(estimation_problem.state_preparation, + list(range(self._m, circuit.num_qubits)), + inplace=True) + circuit.compose(pec, inplace=True) + + # add measurements if necessary + if measurement: + cr = ClassicalRegister(self._m) + circuit.add_register(cr) + circuit.measure(list(range(self._m)), list(range(self._m))) + + return circuit + + def evaluate_measurements(self, circuit_results: Union[Dict[str, int], np.ndarray], + threshold: float = 1e-6, + ) -> Tuple[Dict[int, float], Dict[float, float]]: + """Evaluate the results from the circuit simulation. + + Given the probabilities from statevector simulation of the QAE circuit, compute the + probabilities that the measurements y/gridpoints a are the best estimate. + + Args: + circuit_results: The circuit result from the QAE circuit. Can be either a counts dict + or a statevector. + threshold: Measurements with probabilities below the threshold are discarded. + + Returns: + Dictionaries containing the a gridpoints with respective probabilities and + y measurements with respective probabilities, in this order. + """ + # compute grid sample and measurement dicts + if isinstance(circuit_results, dict): + samples, measurements = self._evaluate_count_results(circuit_results) + else: + samples, measurements = self._evaluate_statevector_results(circuit_results) + + # cutoff probabilities below the threshold + samples = {a: p for a, p in samples.items() if p > threshold} + measurements = {y: p for y, p in measurements.items() if p > threshold} + + return samples, measurements + + def _evaluate_statevector_results(self, statevector): + # map measured results to estimates + measurements = OrderedDict() # type: OrderedDict + num_qubits = int(np.log2(len(statevector))) + for i, amplitude in enumerate(statevector): + b = bin(i)[2:].zfill(num_qubits)[::-1] + y = int(b[:self._m], 2) # chop off all except the evaluation qubits + measurements[y] = measurements.get(y, 0) + np.abs(amplitude) ** 2 + + samples = OrderedDict() # type: OrderedDict + for y, probability in measurements.items(): + if y >= int(self._M / 2): + y = self._M - y + # due to the finite accuracy of the sine, we round the result to 7 decimals + a = np.round(np.power(np.sin(y * np.pi / 2 ** self._m), 2), decimals=7) + samples[a] = samples.get(a, 0) + probability + + return samples, measurements + + def _evaluate_count_results(self, counts): + # construct probabilities + measurements = OrderedDict() + samples = OrderedDict() + shots = self._quantum_instance._run_config.shots + + for state, count in counts.items(): + y = int(state.replace(' ', '')[:self._m][::-1], 2) + probability = count / shots + measurements[y] = probability + a = np.round(np.power(np.sin(y * np.pi / 2 ** self._m), 2), decimals=7) + samples[a] = samples.get(a, 0.0) + probability + + return samples, measurements + + @staticmethod + def compute_mle(result: 'AmplitudeEstimationResult', apply_post_processing: bool = False + ) -> float: + """Compute the Maximum Likelihood Estimator (MLE). + + Args: + result: An amplitude estimation result object. + apply_post_processing: If True, apply the post processing to the MLE before returning + it. + + Returns: + The MLE for the provided result object. + """ + m = result.num_evaluation_qubits + M = 2 ** m # pylint: disable=invalid-name + qae = result.estimation + + # likelihood function + a_i = np.asarray(list(result.samples.keys())) + p_i = np.asarray(list(result.samples.values())) + + def loglikelihood(a): + return np.sum(result.shots * p_i * np.log(pdf_a(a_i, a, m))) + + # y is pretty much an integer, but to map 1.9999 to 2 we must first + # use round and then int conversion + y = int(np.round(M * np.arcsin(np.sqrt(qae)) / np.pi)) + + # Compute the two intervals in which are candidates for containing + # the maximum of the log-likelihood function: the two bubbles next to + # the QAE estimate + if y == 0: + right_of_qae = np.sin(np.pi * (y + 1) / M) ** 2 + bubbles = [qae, right_of_qae] + + elif y == int(M / 2): # remember, M = 2^m is a power of 2 + left_of_qae = np.sin(np.pi * (y - 1) / M) ** 2 + bubbles = [left_of_qae, qae] + + else: + left_of_qae = np.sin(np.pi * (y - 1) / M) ** 2 + right_of_qae = np.sin(np.pi * (y + 1) / M) ** 2 + bubbles = [left_of_qae, qae, right_of_qae] + + # Find global maximum amongst the two local maxima + a_opt = qae + loglik_opt = loglikelihood(a_opt) + for a, b in zip(bubbles[:-1], bubbles[1:]): + locmax, val = bisect_max(loglikelihood, a, b, retval=True) + if val > loglik_opt: + a_opt = locmax + loglik_opt = val + + if apply_post_processing: + return result.post_processing(a_opt) + + return a_opt + + def estimate(self, estimation_problem: EstimationProblem) -> 'AmplitudeEstimationResult': + """Run the amplitude estimation algorithm on provided estimation problem. + + Args: + estimation_problem: The estimation problem. + + Returns: + An amplitude estimation results object. + + Raises: + ValueError: If `state_preparation` or `objective_qubits` are not set in the + `estimation_problem`. + """ + # check if A factory or state_preparation has been set + if estimation_problem.state_preparation is None: + raise ValueError('The state_preparation property of the estimation problem must be ' + 'set.') + + if estimation_problem.objective_qubits is None: + raise ValueError('The objective_qubits property of the estimation problem must be ' + 'set.') + + result = AmplitudeEstimationResult() + result.num_evaluation_qubits = self._m + result.post_processing = estimation_problem.post_processing + + if self._quantum_instance.is_statevector: + circuit = self.construct_circuit(estimation_problem, measurement=False) + # run circuit on statevector simulator + statevector = self._quantum_instance.execute(circuit).get_statevector() + result.circuit_results = statevector + + # store number of shots: convention is 1 shot for statevector, + # needed so that MLE works! + result.shots = 1 + else: + # run circuit on QASM simulator + circuit = self.construct_circuit(estimation_problem, measurement=True) + counts = self._quantum_instance.execute(circuit).get_counts() + result.circuit_results = counts + + # store shots + result.shots = sum(counts.values()) + + samples, measurements = self.evaluate_measurements(result.circuit_results) + + result.samples = samples + result.samples_processed = {estimation_problem.post_processing(a): p + for a, p in samples.items()} + result.measurements = measurements + + # determine the most likely estimate + result.max_probability = 0 + for amplitude, (mapped, prob) in zip(samples.keys(), result.samples_processed.items()): + if prob > result.max_probability: + result.max_probability = prob + result.estimation = amplitude + result.estimation_processed = mapped + + # store the number of oracle queries + result.num_oracle_queries = result.shots * (self._M - 1) + + # run the MLE post processing + mle = self.compute_mle(result) + result.mle = mle + result.mle_processed = estimation_problem.post_processing(mle) + + result.confidence_interval = self.compute_confidence_interval(result) + result.confidence_interval_processed = tuple(estimation_problem.post_processing(value) + for value in result.confidence_interval) + + return result + + @staticmethod + def compute_confidence_interval(result: 'AmplitudeEstimationResult', + alpha: float = 0.05, + kind: str = 'likelihood_ratio' + ) -> Tuple[float, float]: + """Compute the (1 - alpha) confidence interval. + + Args: + result: An amplitude estimation result for which to compute the confidence interval. + alpha: Confidence level: compute the (1 - alpha) confidence interval. + kind: The method to compute the confidence interval, can be 'fisher', 'observed_fisher' + or 'likelihood_ratio' (default) + + Returns: + The (1 - alpha) confidence interval of the specified kind. + + Raises: + AquaError: If 'mle' is not in self._ret.keys() (i.e. `run` was not called yet). + NotImplementedError: If the confidence interval method `kind` is not implemented. + """ + # if statevector simulator the estimate is exact + if isinstance(result.circuit_results, (list, np.ndarray)): + return (result.mle, result.mle) + + if kind in ['likelihood_ratio', 'lr']: + return _likelihood_ratio_confint(result, alpha) + + if kind in ['fisher', 'fi']: + return _fisher_confint(result, alpha, observed=False) + + if kind in ['observed_fisher', 'observed_information', 'oi']: + return _fisher_confint(result, alpha, observed=True) + + raise NotImplementedError('CI `{}` is not implemented.'.format(kind)) + + +class AmplitudeEstimationResult(AmplitudeEstimatorResult): + """The ``AmplitudeEstimation`` result object.""" + + def __init__(self) -> None: + super().__init__() + self._num_evaluation_qubits = None + self._mle = None + self._mle_processed = None + self._samples = None + self._samples_processed = None + self._y_measurements = None + self._max_probability = None + + @property + def num_evaluation_qubits(self) -> int: + """Returns the number of evaluation qubits.""" + return self._num_evaluation_qubits + + @num_evaluation_qubits.setter + def num_evaluation_qubits(self, num_evaluation_qubits: int) -> None: + """Set the number of evaluation qubits.""" + self._num_evaluation_qubits = num_evaluation_qubits + + @property + def mle_processed(self) -> float: + """Return the post-processed MLE for the amplitude.""" + return self._mle_processed + + @mle_processed.setter + def mle_processed(self, value: float) -> None: + """Set the post-processed MLE for the amplitude.""" + self._mle_processed = value + + @property + def samples_processed(self) -> Dict[float, float]: + """Return the post-processed measurement samples with their measurement probability.""" + return self._samples_processed + + @samples_processed.setter + def samples_processed(self, value: Dict[float, float]) -> None: + """Set the post-processed measurement samples.""" + self._samples_processed = value + + @property + def mle(self) -> float: + r"""Return the MLE for the amplitude, in $[0, 1]$.""" + return self._mle + + @mle.setter + def mle(self, value: float) -> None: + r"""Set the MLE for the amplitude, in $[0, 1]$.""" + self._mle = value + + @property + def samples(self) -> Dict[float, float]: + """Return the measurement samples with their measurement probability.""" + return self._samples + + @samples.setter + def samples(self, value: Dict[float, float]) -> None: + """Set the measurement samples with their measurement probability.""" + self._samples = value + + @property + def measurements(self) -> Dict[int, float]: + """Return the measurements as integers with their measurement probability.""" + return self._y_measurements + + @measurements.setter + def measurements(self, value: Dict[int, float]) -> None: + """Set the measurements as integers with their measurement probability.""" + self._y_measurements = value + + @property + def max_probability(self) -> float: + """Return the maximum sampling probability.""" + return self._max_probability + + @max_probability.setter + def max_probability(self, value: float) -> None: + """Set the maximum sampling probability.""" + self._max_probability = value + + +def _compute_fisher_information(result: AmplitudeEstimationResult, observed: bool = False) -> float: + """Computes the Fisher information for the output of the previous run. + + Args: + result: An amplitude estimation result for which to compute the confidence interval. + observed: If True, the observed Fisher information is returned, otherwise + the expected Fisher information. + + Returns: + The Fisher information. + """ + fisher_information = None + mlv = result.mle # MLE in [0,1] + m = result.num_evaluation_qubits + M = 2 ** m # pylint: disable=invalid-name + + if observed: + a_i = np.asarray(list(result.samples.keys())) + p_i = np.asarray(list(result.samples.values())) + + # Calculate the observed Fisher information + fisher_information = sum(p * derivative_log_pdf_a(a, mlv, m) ** 2 + for p, a in zip(p_i, a_i)) + else: + def integrand(x): + return (derivative_log_pdf_a(x, mlv, m))**2 * pdf_a(x, mlv, m) + + grid = np.sin(np.pi * np.arange(M / 2 + 1) / M) ** 2 + fisher_information = sum(integrand(x) for x in grid) + + return fisher_information + + +def _fisher_confint(result: AmplitudeEstimationResult, alpha: float, observed: bool = False + ) -> List[float]: + """Compute the Fisher information confidence interval for the MLE of the previous run. + + Args: + result: An amplitude estimation result for which to compute the confidence interval. + alpha: Specifies the (1 - alpha) confidence level (0 < alpha < 1). + observed: If True, the observed Fisher information is used to construct the + confidence interval, otherwise the expected Fisher information. + + Returns: + The Fisher information confidence interval. + """ + # approximate the standard deviation of the MLE and construct the confidence interval + std = np.sqrt(result.shots * _compute_fisher_information(result, observed)) + confint = result.mle + norm.ppf(1 - alpha / 2) / std * np.array([-1, 1]) + + # transform the confidence interval from [0, 1] to the target interval + return tuple(result.post_processing(bound) for bound in confint) + + +def _likelihood_ratio_confint(result: AmplitudeEstimationResult, + alpha: float) -> List[float]: + """Compute the likelihood ratio confidence interval for the MLE of the previous run. + + Args: + result: An amplitude estimation result for which to compute the confidence interval. + alpha: Specifies the (1 - alpha) confidence level (0 < alpha < 1). + + Returns: + The likelihood ratio confidence interval. + """ + # Compute the two intervals in which we the look for values above + # the likelihood ratio: the two bubbles next to the QAE estimate + m = result.num_evaluation_qubits + M = 2 ** m # pylint: disable=invalid-name + qae = result.estimation + + y = int(np.round(M * np.arcsin(np.sqrt(qae)) / np.pi)) + if y == 0: + right_of_qae = np.sin(np.pi * (y + 1) / M)**2 + bubbles = [qae, right_of_qae] + + elif y == int(M / 2): # remember, M = 2^m is a power of 2 + left_of_qae = np.sin(np.pi * (y - 1) / M)**2 + bubbles = [left_of_qae, qae] + + else: + left_of_qae = np.sin(np.pi * (y - 1) / M)**2 + right_of_qae = np.sin(np.pi * (y + 1) / M)**2 + bubbles = [left_of_qae, qae, right_of_qae] + + # likelihood function + a_i = np.asarray(list(result.samples.keys())) + p_i = np.asarray(list(result.samples.values())) + + def loglikelihood(a): + return np.sum(result.shots * p_i * np.log(pdf_a(a_i, a, m))) + + # The threshold above which the likelihoods are in the + # confidence interval + loglik_mle = loglikelihood(result.mle) + thres = loglik_mle - chi2.ppf(1 - alpha, df=1) / 2 + + def cut(x): + return loglikelihood(x) - thres + + # Store the boundaries of the confidence interval + # It's valid to start off with the zero-width confidence interval, since the maximum + # of the likelihood function is guaranteed to be over the threshold, and if alpha = 0 + # that's the valid interval + lower = upper = result.mle + + # Check the two intervals/bubbles: check if they surpass the + # threshold and if yes add the part that does to the CI + for a, b in zip(bubbles[:-1], bubbles[1:]): + # Compute local maximum and perform a bisect search between + # the local maximum and the bubble boundaries + locmax, val = bisect_max(loglikelihood, a, b, retval=True) + if val >= thres: + # Bisect pre-condition is that the function has different + # signs at the boundaries of the interval we search in + if cut(a) * cut(locmax) < 0: + left = bisect(cut, a, locmax) + lower = np.minimum(lower, left) + if cut(locmax) * cut(b) < 0: + right = bisect(cut, locmax, b) + upper = np.maximum(upper, right) + + # Put together CI + confint = [lower, upper] + return tuple(result.post_processing(bound) for bound in confint) diff --git a/qiskit/algorithms/amplitude_estimators/ae_utils.py b/qiskit/algorithms/amplitude_estimators/ae_utils.py new file mode 100644 index 000000000000..aeebe50d655a --- /dev/null +++ b/qiskit/algorithms/amplitude_estimators/ae_utils.py @@ -0,0 +1,249 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2018, 2020. +# +# 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. + +"""Utils for the Maximum-Likelihood estimation used in ``AmplitudeEstimation``.""" + +import logging +import numpy as np + +logger = logging.getLogger(__name__) + +# pylint: disable=invalid-name + + +def bisect_max(f, a, b, steps=50, minwidth=1e-12, retval=False): + """Find the maximum of the real-valued function f in the interval [a, b] using bisection. + + Args: + f (callable): the function to find the maximum of + a (float): the lower limit of the interval + b (float): the upper limit of the interval + steps (int): the maximum number of steps in the bisection + minwidth (float): if the current interval is smaller than minwidth stop + the search + retval (bool): return value + + Returns: + float: The maximum of f in [a,b] according to this algorithm. + """ + it = 0 + m = (a + b) / 2 + fm = 0 + while it < steps and b - a > minwidth: + l, r = (a + m) / 2, (m + b) / 2 + fl, fm, fr = f(l), f(m), f(r) + + # fl is the maximum + if fl > fm and fl > fr: + b = m + m = l + # fr is the maximum + elif fr > fm and fr > fl: + a = m + m = r + # fm is the maximum + else: + a = l + b = r + + it += 1 + + if it == steps: + logger.warning("-- Warning, bisect_max didn't converge after %s steps", steps) + + if retval: + return m, fm + + return m + + +def _circ_dist(x, p): + r"""Circumferential distance function. + + For two angles :math:`x` and :math:`p` on the unit circuit this function is defined as + + .. math:: + + d(x, p) = \min_{z \in [-1, 0, 1]} |z + p - x| + + Args: + x (float): first angle + p (float): second angle + + Returns: + float: d(x, p) + """ + t = p - x + # Since x and p \in [0,1] it suffices to check not all integers + # but only -1, 0 and 1 + z = np.array([-1, 0, 1]) + + if hasattr(t, "__len__"): + d = np.empty_like(t) + for idx, ti in enumerate(t): + d[idx] = np.min(np.abs(z + ti)) + return d + + return np.min(np.abs(z + t)) + + +def _derivative_circ_dist(x, p): + """Derivative of circumferential distance function. + + Args: + x (float): first angle + p (float): second angle + + Returns: + float: The derivative. + """ + # pylint: disable=chained-comparison,misplaced-comparison-constant + t = p - x + if t < -0.5 or (0 < t and t < 0.5): + return -1 + if t > 0.5 or (-0.5 < t and t < 0): + return 1 + return 0 + + +def _amplitude_to_angle(a): + r"""Transform from the amplitude :math:`a \in [0, 1]` to the generating angle. + + In QAE, the amplitude can be written from a generating angle :math:`\omega` as + + .. math: + + a = \sin^2(\pi \omega) + + This returns the :math:`\omega` for a given :math:`a`. + + Args: + a (float): A value in :math:`[0,1]`. + + Returns: + float: :math:`\sin^{-1}(\sqrt{a}) / \pi` + """ + return np.arcsin(np.sqrt(a)) / np.pi + + +def _derivative_amplitude_to_angle(a): + """Compute the derivative of ``amplitude_to_angle``.""" + return 1 / (2 * np.pi * np.sqrt((1 - a) * a)) + + +def _alpha(x, p): + """Helper function for `pdf_a`, alpha = pi * d(omega(x), omega(p)). + + Here, omega(x) is `_amplitude_to_angle(x)`. + """ + omega = _amplitude_to_angle + return np.pi * _circ_dist(omega(x), omega(p)) + + +def _derivative_alpha(x, p): + """Compute the derivative of alpha.""" + omega = _amplitude_to_angle + d_omega = _derivative_amplitude_to_angle + return np.pi * _derivative_circ_dist(omega(x), omega(p)) * d_omega(p) + + +def _beta(x, p): + """Helper function for `pdf_a`, beta = pi * d(1 - omega(x), omega(p)).""" + omega = _amplitude_to_angle + return np.pi * _circ_dist(1 - omega(x), omega(p)) + + +def _derivative_beta(x, p): + """Compute the derivative of beta.""" + omega = _amplitude_to_angle + d_omega = _derivative_amplitude_to_angle + return np.pi * _derivative_circ_dist(1 - omega(x), omega(p)) * d_omega(p) + + +def _pdf_a_single_angle(x, p, m, pi_delta): + """Helper function for `pdf_a`.""" + M = 2**m + + d = pi_delta(x, p) + res = np.sin(M * d)**2 / (M * np.sin(d))**2 if d != 0 else 1 + + return res + + +def pdf_a(x, p, m): + """ + Return the PDF of a, i.e. the probability of getting the estimate x + (in [0, 1]) if p (in [0, 1]) is the true value, given that we use m qubits. + + Args: + x (float): the grid point + p (float): the true value + m (float): the number of evaluation qubits + + Returns: + float: PDF(x|p) + """ + # We'll use list comprehension, so the input should be a list + scalar = False + if not hasattr(x, "__len__"): + scalar = True + x = np.asarray([x]) + + # Compute the probabilities: Add up both angles that produce the given + # value, except for the angles 0 and 0.5, which map to the unique a-values, + # 0 and 1, respectively + pr = np.array([_pdf_a_single_angle(xi, p, m, _alpha) + _pdf_a_single_angle(xi, p, m, _beta) + if (xi not in [0, 1]) else _pdf_a_single_angle(xi, p, m, _alpha) + for xi in x + ]).flatten() + + # If is was a scalar return scalar otherwise the array + return pr[0] if scalar else pr + + +def derivative_log_pdf_a(x, p, m): + """ + Return the derivative of the logarithm of the PDF of a. + + Args: + x (float): the grid point + p (float): the true value + m (float): the number of evaluation qubits + + Returns: + float: d/dp log(PDF(x|p)) + """ + M = 2**m + + if x not in [0, 1]: + num_p1 = 0 + for A, dA, B, dB in zip([_alpha, _beta], + [_derivative_alpha, _derivative_beta], + [_beta, _alpha], [_derivative_beta, _derivative_alpha]): + num_p1 += 2 * M * np.sin(M * A(x, p)) * np.cos(M * A(x, p)) \ + * dA(x, p) * np.sin(B(x, p))**2 \ + + 2 * np.sin(M * A(x, p))**2 * np.sin(B(x, p)) * np.cos(B(x, p)) * dB(x, p) + + den_p1 = np.sin(M * _alpha(x, p))**2 * np.sin(_beta(x, p))**2 + \ + np.sin(M * _beta(x, p))**2 * np.sin(_alpha(x, p))**2 + + num_p2 = 0 + for A, dA, B in zip([_alpha, _beta], + [_derivative_alpha, _derivative_beta], + [_beta, _alpha]): + num_p2 += 2 * np.cos(A(x, p)) * dA(x, p) * np.sin(B(x, p)) + + den_p2 = np.sin(_alpha(x, p)) * np.sin(_beta(x, p)) + + return num_p1 / den_p1 - num_p2 / den_p2 + + return 2 * _derivative_alpha(x, p) * (M / np.tan(M * _alpha(x, p)) - 1 / np.tan(_alpha(x, p))) diff --git a/qiskit/algorithms/amplitude_estimators/amplitude_estimator.py b/qiskit/algorithms/amplitude_estimators/amplitude_estimator.py new file mode 100644 index 000000000000..7cdcfa29973c --- /dev/null +++ b/qiskit/algorithms/amplitude_estimators/amplitude_estimator.py @@ -0,0 +1,129 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2018, 2020. +# +# 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 Amplitude Estimation interface.""" + +from abc import abstractmethod +from typing import Union, Optional, Dict, Callable, Tuple +import numpy as np + +from .estimation_problem import EstimationProblem +from ..algorithm_result import AlgorithmResult + + +class AmplitudeEstimator: + """The Amplitude Estimation interface.""" + + @abstractmethod + def estimate(self, estimation_problem: EstimationProblem) -> 'AmplitudeEstimatorResult': + """Run the amplitude estimation algorithm. + + Args: + estimation_problem: An ``EstimationProblem`` containing all problem-relevant information + such as the state preparation and the objective qubits. + """ + raise NotImplementedError + + +class AmplitudeEstimatorResult(AlgorithmResult): + """The results object for amplitude estimation algorithms.""" + + def __init__(self) -> None: + super().__init__() + self._circuit_results = None + self._shots = None + self._estimation = None + self._estimation_processed = None + self._num_oracle_queries = None + self._post_processing = None + self._confidence_interval = None + self._confidence_interval_processed = None + + @property + def circuit_results(self) -> Optional[Union[np.ndarray, Dict[str, int]]]: + """Return the circuit results. Can be a statevector or counts dictionary.""" + return self._circuit_results + + @circuit_results.setter + def circuit_results(self, value: Union[np.ndarray, Dict[str, int]]) -> None: + """Set the circuit results.""" + self._circuit_results = value + + @property + def shots(self) -> int: + """Return the number of shots used. Is 1 for statevector-based simulations.""" + return self._shots + + @shots.setter + def shots(self, value: int) -> None: + """Set the number of shots used.""" + self._shots = value + + @property + def estimation(self) -> float: + r"""Return the estimation for the amplitude in :math:`[0, 1]`.""" + return self._estimation + + @estimation.setter + def estimation(self, value: float) -> None: + r"""Set the estimation for the amplitude in :math:`[0, 1]`.""" + self._estimation = value + + @property + def estimation_processed(self) -> float: + """Return the estimation for the amplitude after the post-processing has been applied.""" + return self._estimation_processed + + @estimation_processed.setter + def estimation_processed(self, value: float) -> None: + """Set the estimation for the amplitude after the post-processing has been applied.""" + self._estimation_processed = value + + @property + def num_oracle_queries(self) -> int: + """Return the number of Grover oracle queries.""" + return self._num_oracle_queries + + @num_oracle_queries.setter + def num_oracle_queries(self, value: int) -> None: + """Set the number of Grover oracle queries.""" + self._num_oracle_queries = value + + @property + def post_processing(self) -> Callable[[float], float]: + """Return a handle to the post processing function.""" + return self._post_processing + + @post_processing.setter + def post_processing(self, post_processing: Callable[[float], float]) -> None: + """Set a handle to the post processing function.""" + self._post_processing = post_processing + + @property + def confidence_interval(self) -> Tuple[float, float]: + """Return the confidence interval for the amplitude (95% interval by default).""" + return self._confidence_interval + + @confidence_interval.setter + def confidence_interval(self, confidence_interval: Tuple[float, float]) -> None: + """Set the confidence interval for the amplitude (95% interval by default).""" + self._confidence_interval = confidence_interval + + @property + def confidence_interval_processed(self) -> Tuple[float, float]: + """Return the post-processed confidence interval (95% interval by default).""" + return self._confidence_interval_processed + + @confidence_interval_processed.setter + def confidence_interval_processed(self, confidence_interval: Tuple[float, float]) -> None: + """Set the post-processed confidence interval (95% interval by default).""" + self._confidence_interval_processed = confidence_interval diff --git a/qiskit/algorithms/amplitude_estimators/estimation_problem.py b/qiskit/algorithms/amplitude_estimators/estimation_problem.py new file mode 100644 index 000000000000..8819f26e355f --- /dev/null +++ b/qiskit/algorithms/amplitude_estimators/estimation_problem.py @@ -0,0 +1,256 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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 Estimation problem class.""" + +import warnings +from typing import Optional, List, Callable +import numpy + +from qiskit.circuit import QuantumCircuit, QuantumRegister +from qiskit.circuit.library import GroverOperator + + +class EstimationProblem: + """The estimation problem is the input to amplitude estimation algorithm. + + This class contains all problem-specific information required to run an amplitude estimation + algorithm. That means, it minimally contains the state preparation and the specification + of the good state. It can further hold some post processing on the estimation of the amplitude + or a custom Grover operator. + """ + + def __init__(self, + state_preparation: QuantumCircuit, + objective_qubits: List[int], + grover_operator: Optional[QuantumCircuit] = None, + post_processing: Optional[Callable[[float], float]] = None, + is_good_state: Optional[Callable[[str], bool]] = None, + ) -> None: + r""" + Args: + state_preparation: A circuit preparing the input state, referred to as + :math:`\mathcal{A}`. + objective_qubits: A list of qubit indices to specify the oracle in the Grover operator, + if the Grover operator is not supplied. A measurement outcome is classified as + 'good' state if all objective qubits are in state :math:`|1\rangle`, otherwise it + is classified as 'bad'. + grover_operator: The Grover operator :math:`\mathcal{Q}` used as unitary in the + phase estimation circuit. + post_processing: A mapping applied to the result of the algorithm + :math:`0 \leq a \leq 1`, usually used to map the estimate to a target interval. + Defaults to the identity. + is_good_state: A function to check whether a string represents a good state. Defaults + to all objective qubits being in state 1. + """ + self._state_preparation = state_preparation + self._objective_qubits = objective_qubits + self._grover_operator = grover_operator + self._post_processing = post_processing + self._is_good_state = is_good_state + + @property + def state_preparation(self) -> Optional[QuantumCircuit]: + r"""Get the :math:`\mathcal{A}` operator encoding the amplitude :math:`a`. + + Returns: + The :math:`\mathcal{A}` operator as `QuantumCircuit`. + """ + return self._state_preparation + + @state_preparation.setter + def state_preparation(self, state_preparation: QuantumCircuit) -> None: + r"""Set the :math:`\mathcal{A}` operator, that encodes the amplitude to be estimated. + + Args: + state_preparation: The new :math:`\mathcal{A}` operator. + """ + self._state_preparation = state_preparation + + @property + def objective_qubits(self) -> List[int]: + """Get the criterion for a measurement outcome to be in a 'good' state. + + Returns: + The criterion as list of qubit indices. + """ + return self._objective_qubits + + @objective_qubits.setter + def objective_qubits(self, objective_qubits: List[int]) -> None: + """Set the criterion for a measurement outcome to be in a 'good' state. + + Args: + objective_qubits: The criterion as callable of list of qubit indices. + """ + self._objective_qubits = objective_qubits + + @property + def post_processing(self) -> Callable[[float], float]: + """Apply post processing to the input value. + + Returns: + A handle to the post processing function. Acts as identity by default. + """ + if self._post_processing is None: + return lambda x: x + + return self._post_processing + + @post_processing.setter + def post_processing(self, post_processing: Optional[Callable[[float], float]]) -> None: + """Set the post processing function. + + Args: + post_processing: A handle to the post processing function. If set to ``None``, the + identity will be used as post processing. + """ + self._post_processing = post_processing + + @property + def is_good_state(self) -> Callable[[str], bool]: + """Checks whether a bitstring represents a good state. + + Returns: + Handle to the ``is_good_state`` callable. + """ + if self._is_good_state is None: + return lambda x: all(bit == '1' for bit in x) + + return self._is_good_state + + @is_good_state.setter + def is_good_state(self, is_good_state: Optional[Callable[[str], bool]]) -> None: + """Set the ``is_good_state`` function. + + Args: + is_good_state: A function to determine whether a bitstring represents a good state. + If set to ``None``, the good state will be defined as all bits being one. + """ + self._is_good_state = is_good_state + + @property + def grover_operator(self) -> Optional[QuantumCircuit]: + r"""Get the :math:`\mathcal{Q}` operator, or Grover operator. + + If the Grover operator is not set, we try to build it from the :math:`\mathcal{A}` operator + and `objective_qubits`. This only works if `objective_qubits` is a list of integers. + + Returns: + The Grover operator, or None if neither the Grover operator nor the + :math:`\mathcal{A}` operator is set. + """ + if self._grover_operator is not None: + return self._grover_operator + + # build the reflection about the bad state: a MCZ with open controls (thus X gates + # around the controls) and X gates around the target to change from a phaseflip on + # |1> to a phaseflip on |0> + num_state_qubits = self.state_preparation.num_qubits \ + - self.state_preparation.num_ancillas + + oracle = QuantumCircuit(num_state_qubits) + oracle.h(self.objective_qubits[-1]) + if len(self.objective_qubits) == 1: + oracle.x(self.objective_qubits[0]) + else: + oracle.mcx(self.objective_qubits[:-1], self.objective_qubits[-1]) + oracle.h(self.objective_qubits[-1]) + + # construct the grover operator + return GroverOperator(oracle, self.state_preparation) + + @grover_operator.setter + def grover_operator(self, grover_operator: Optional[QuantumCircuit]) -> None: + r"""Set the :math:`\mathcal{Q}` operator. + + Args: + grover_operator: The new :math:`\mathcal{Q}` operator. If set to ``None``, + the default construction via ``qiskit.circuit.library.GroverOperator`` is used. + """ + self._grover_operator = grover_operator + + def rescale(self, scaling_factor: float) -> 'EstimationProblem': + """Rescale the good state amplitude in the estimation problem. + + Args: + scaling_factor: The scaling factor in [0, 1]. + + Returns: + A rescaled estimation problem. + """ + if self._grover_operator is not None: + warnings.warn('Rescaling discards the Grover operator.') + + # rescale the amplitude by a factor of 1/4 by adding an auxiliary qubit + rescaled_stateprep = _rescale_amplitudes(self.state_preparation, scaling_factor) + num_qubits = self.state_preparation.num_qubits + objective_qubits = self.objective_qubits + [num_qubits] + + # add the scaling qubit to the good state qualifier + def is_good_state(bitstr): + # pylint:disable=not-callable + return self.is_good_state(bitstr[1:]) and bitstr[0] == '1' + + # rescaled estimation problem + problem = EstimationProblem(rescaled_stateprep, + objective_qubits=objective_qubits, + post_processing=self.post_processing, + is_good_state=is_good_state) + + return problem + + +def _rescale_amplitudes(circuit: QuantumCircuit, scaling_factor: float) -> QuantumCircuit: + r"""Uses an auxiliary qubit to scale the amplitude of :math:`|1\rangle` by ``scaling_factor``. + + Explained in Section 2.1. of [1]. + + For example, for a scaling factor of 0.25 this turns this circuit + + .. code-block:: + + ┌───┐ + state_0: ─────┤ H ├─────────■──── + ┌───┴───┴───┐ ┌───┴───┐ + obj_0: ─┤ RY(0.125) ├─┤ RY(1) ├ + └───────────┘ └───────┘ + + into + + .. code-block:: + + ┌───┐ + state_0: ─────┤ H ├─────────■──── + ┌───┴───┴───┐ ┌───┴───┐ + obj_0: ─┤ RY(0.125) ├─┤ RY(1) ├ + ┌┴───────────┴┐└───────┘ + scaling_0: ┤ RY(0.50536) ├───────── + └─────────────┘ + + References: + + [1]: K. Nakaji. Faster Amplitude Estimation, 2020; + `arXiv:2002.02417 `_ + + Args: + circuit: The circuit whose amplitudes to rescale. + scaling_factor: The rescaling factor. + + Returns: + A copy of the circuit with an additional qubit and RY gate for the rescaling. + """ + qr = QuantumRegister(1, 'scaling') + rescaled = QuantumCircuit(*circuit.qregs, qr) + rescaled.compose(circuit, circuit.qubits, inplace=True) + rescaled.ry(2 * numpy.arcsin(scaling_factor), qr) + return rescaled diff --git a/qiskit/algorithms/amplitude_estimators/fae.py b/qiskit/algorithms/amplitude_estimators/fae.py new file mode 100644 index 000000000000..70665ab0e6f8 --- /dev/null +++ b/qiskit/algorithms/amplitude_estimators/fae.py @@ -0,0 +1,304 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2017, 2020. +# +# 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. + +"""Faster Amplitude Estimation.""" + +from typing import Optional, Union, List, Tuple +import numpy as np + +from qiskit.circuit import QuantumCircuit, ClassicalRegister +from qiskit.providers import BaseBackend, Backend +from qiskit.utils import QuantumInstance +from qiskit.algorithms.exceptions import AlgorithmError + +from .amplitude_estimator import AmplitudeEstimator, AmplitudeEstimatorResult +from .estimation_problem import EstimationProblem + + +class FasterAmplitudeEstimation(AmplitudeEstimator): + """The Faster Amplitude Estimation algorithm. + + The Faster Amplitude Estimation (FAE) [1] algorithm is a variant of Quantum Amplitude + Estimation (QAE), where the Quantum Phase Estimation (QPE) by an iterative Grover search, + similar to [2]. + + Due to the iterative version of the QPE, this algorithm does not require any additional + qubits, as the originally proposed QAE [3] and thus the resulting circuits are less complex. + + References: + + [1]: K. Nakaji. Faster Amplitude Estimation, 2020; + `arXiv:2002.02417 `_ + [2]: D. Grinko et al. Iterative Amplitude Estimation, 2019; + `arXiv:1912.05559 `_ + [3]: G. Brassard et al. Quantum Amplitude Amplification and Estimation, 2000; + `arXiv:quant-ph/0005055 `_ + + """ + + def __init__(self, + delta: float, + maxiter: int, + rescale: bool = True, + quantum_instance: Optional[Union[QuantumInstance, BaseBackend, Backend]] = None + ) -> None: + r""" + Args: + delta: The probability that the true value is outside of the final confidence interval. + maxiter: The number of iterations, the maximal power of Q is `2 ** (maxiter - 1)`. + rescale: Whether to rescale the problem passed to `estimate`. + quantum_instance: The quantum instance or backend to run the circuits. + + .. note:: + + This algorithm overwrites the number of shots set in the ``quantum_instance`` + argument, but will reset them to the initial number after running. + + """ + super().__init__() + self.quantum_instance = quantum_instance + self._shots = (int(1944 * np.log(2 / delta)), int(972 * np.log(2 / delta))) + self._rescale = rescale + self._delta = delta + self._maxiter = maxiter + self._num_oracle_calls = 0 + + @property + def quantum_instance(self) -> Optional[QuantumInstance]: + """Get the quantum instance. + + Returns: + The quantum instance used to run this algorithm. + """ + return self._quantum_instance + + @quantum_instance.setter + def quantum_instance(self, quantum_instance: Union[QuantumInstance, + BaseBackend, Backend]) -> None: + """Set quantum instance. + + Args: + quantum_instance: The quantum instance used to run this algorithm. + """ + if isinstance(quantum_instance, (BaseBackend, Backend)): + quantum_instance = QuantumInstance(quantum_instance) + self._quantum_instance = quantum_instance + + def _cos_estimate(self, estimation_problem, k, shots): + if self._quantum_instance is None: + raise AlgorithmError('Quantum instance must be set.') + + if self._quantum_instance.is_statevector: + circuit = self.construct_circuit(estimation_problem, k, measurement=False) + statevector = self._quantum_instance.execute(circuit).get_statevector() + + # sum over all amplitudes where the objective qubits are 1 + prob = 0 + for i, amplitude in enumerate(statevector): + # get bitstring of objective qubits + full_state = bin(i)[2:].zfill(circuit.num_qubits)[::-1] + state = ''.join([full_state[i] for i in estimation_problem.objective_qubits]) + + # check if it is a good state + if estimation_problem.is_good_state(state[::-1]): + prob = prob + np.abs(amplitude) ** 2 + + cos_estimate = 1 - 2 * prob + else: + circuit = self.construct_circuit(estimation_problem, k, measurement=True) + + self._quantum_instance.run_config.shots = shots + counts = self._quantum_instance.execute(circuit).get_counts() + self._num_oracle_calls += (2 * k + 1) * shots + + good_counts = 0 + for state, count in counts.items(): + if estimation_problem.is_good_state(state): + good_counts += count + + cos_estimate = 1 - 2 * good_counts / shots + + return cos_estimate + + def _chernoff(self, cos, shots): + width = np.sqrt(np.log(2 / self._delta) * 12 / shots) + confint = [np.maximum(-1, cos - width), np.minimum(1, cos + width)] + return confint + + def construct_circuit(self, estimation_problem: EstimationProblem, k: int, + measurement: bool = False + ) -> Union[QuantumCircuit, Tuple[QuantumCircuit, List[int]]]: + r"""Construct the circuit :math:`Q^k X |0\rangle>`. + + The A operator is the unitary specifying the QAE problem and Q the associated Grover + operator. + + Args: + estimation_problem: The estimation problem for which to construct the circuit. + k: The power of the Q operator. + measurement: Boolean flag to indicate if measurements should be included in the + circuits. + + Returns: + The circuit :math:`Q^k X |0\rangle`. + """ + num_qubits = max(estimation_problem.state_preparation.num_qubits, + estimation_problem.grover_operator.num_qubits) + circuit = QuantumCircuit(num_qubits, name='circuit') + + # add classical register if needed + if measurement: + c = ClassicalRegister(len(estimation_problem.objective_qubits)) + circuit.add_register(c) + + # add A operator + circuit.compose(estimation_problem.state_preparation, inplace=True) + + # add Q^k + if k != 0: + circuit.compose(estimation_problem.grover_operator.power(k), inplace=True) + + # add optional measurement + if measurement: + # real hardware can currently not handle operations after measurements, which might + # happen if the circuit gets transpiled, hence we're adding a safeguard-barrier + circuit.barrier() + circuit.measure(estimation_problem.objective_qubits, c[:]) + + return circuit + + def estimate(self, estimation_problem: EstimationProblem) -> 'FasterAmplitudeEstimationResult': + self._num_oracle_calls = 0 + user_defined_shots = self.quantum_instance._run_config.shots + + if self._rescale: + problem = estimation_problem.rescale(0.25) + else: + problem = estimation_problem + + if self._quantum_instance.is_statevector: + cos = self._cos_estimate(problem, k=0, shots=1) + theta = np.arccos(cos) / 2 + theta_ci = [theta, theta] + theta_cis = [theta_ci] + num_steps = num_first_stage_steps = 1 + + else: + theta_ci = [0, np.arcsin(0.25)] + first_stage = True + j_0 = self._maxiter + + theta_cis = [theta_ci] + num_first_stage_steps = 0 + num_steps = 0 + + def cos_estimate(power, shots): + return self._cos_estimate(problem, power, shots) + + for j in range(1, self._maxiter + 1): + num_steps += 1 + if first_stage: + num_first_stage_steps += 1 + c = cos_estimate(2**(j - 1), self._shots[0]) + chernoff_ci = self._chernoff(c, self._shots[0]) + theta_ci = [np.arccos(x) / (2 ** (j + 1) + 2) for x in chernoff_ci[::-1]] + + if 2 ** (j + 1) * theta_ci[1] >= 3 * np.pi / 8 and j < self._maxiter: + j_0 = j + v = 2 ** j * np.sum(theta_ci) + first_stage = False + else: + cos = cos_estimate(2**(j - 1), self._shots[1]) + cos_2 = cos_estimate(2 ** (j - 1) + 2 ** (j_0 - 1), self._shots[1]) + sin = (cos * np.cos(v) - cos_2) / np.sin(v) + rho = np.arctan2(sin, cos) + n = int(((2 ** (j + 1) + 2) * theta_ci[1] - rho + np.pi / 3) / (2 * np.pi)) + + theta_ci = [(2 * np.pi * n + rho + sign * np.pi / 3) / (2 ** (j + 1) + 2) + for sign in [-1, 1]] + theta_cis.append(theta_ci) + + theta = np.mean(theta_ci) + rescaling = 4 if self._rescale else 1 + value = (rescaling * np.sin(theta)) ** 2 + value_ci = [(rescaling * np.sin(x)) ** 2 for x in theta_ci] + + result = FasterAmplitudeEstimationResult() + result.num_oracle_queries = self._num_oracle_calls + result.num_steps = num_steps + result.num_first_state_steps = num_first_stage_steps + if self._quantum_instance.is_statevector: + result.success_probability = 1 + else: + result.success_probability = 1 - (2 * self._maxiter - j_0) * self._delta + + result.estimation = value + result.estimation_processed = problem.post_processing(value) + result.confidence_interval = value_ci + result.confidence_interval_processed = tuple(problem.post_processing(x) + for x in value_ci) + result.theta_intervals = theta_cis + + # reset shots to what the user had defined + self.quantum_instance._run_config.shots = user_defined_shots + return result + + +class FasterAmplitudeEstimationResult(AmplitudeEstimatorResult): + """The result object for the Faster Amplitude Estimation algorithm.""" + + def __init__(self) -> None: + super().__init__() + self._success_probability = None + self._num_steps = None + self._num_first_state_steps = None + self._theta_intervals = None + + @property + def success_probability(self) -> int: + """Return the success probability of the algorithm.""" + return self._success_probability + + @success_probability.setter + def success_probability(self, probability: int) -> None: + """Set the success probability of the algorithm.""" + self._success_probability = probability + + @property + def num_steps(self) -> int: + """Return the total number of steps taken in the algorithm.""" + return self._num_steps + + @num_steps.setter + def num_steps(self, num_steps: int) -> None: + """Set the total number of steps taken in the algorithm.""" + self._num_steps = num_steps + + @property + def num_first_state_steps(self) -> int: + """Return the number of steps taken in the first step of algorithm.""" + return self._num_first_state_steps + + @num_first_state_steps.setter + def num_first_state_steps(self, num_steps: int) -> None: + """Set the number of steps taken in the first step of algorithm.""" + self._num_first_state_steps = num_steps + + @property + def theta_intervals(self) -> List[List[float]]: + """Return the confidence intervals for the angles in each iteration.""" + return self._theta_intervals + + @theta_intervals.setter + def theta_intervals(self, value: List[List[float]]) -> None: + """Set the confidence intervals for the angles in each iteration.""" + self._theta_intervals = value diff --git a/qiskit/algorithms/amplitude_estimators/iae.py b/qiskit/algorithms/amplitude_estimators/iae.py new file mode 100644 index 000000000000..8112440e681a --- /dev/null +++ b/qiskit/algorithms/amplitude_estimators/iae.py @@ -0,0 +1,563 @@ + +# This code is part of Qiskit. +# +# (C) Copyright IBM 2018, 2020. +# +# 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 Iterative Quantum Amplitude Estimation Algorithm.""" + +from typing import Optional, Union, List, Tuple, Dict, cast +import numpy as np +from scipy.stats import beta + +from qiskit import ClassicalRegister, QuantumCircuit +from qiskit.providers import BaseBackend, Backend +from qiskit.utils import QuantumInstance + +from .amplitude_estimator import AmplitudeEstimator, AmplitudeEstimatorResult +from .estimation_problem import EstimationProblem +from ..exceptions import AlgorithmError + + +class IterativeAmplitudeEstimation(AmplitudeEstimator): + r"""The Iterative Amplitude Estimation algorithm. + + This class implements the Iterative Quantum Amplitude Estimation (IQAE) algorithm, proposed + in [1]. The output of the algorithm is an estimate that, + with at least probability :math:`1 - \alpha`, differs by epsilon to the target value, where + both alpha and epsilon can be specified. + + It differs from the original QAE algorithm proposed by Brassard [2] in that it does not rely on + Quantum Phase Estimation, but is only based on Grover's algorithm. IQAE iteratively + applies carefully selected Grover iterations to find an estimate for the target amplitude. + + References: + [1]: Grinko, D., Gacon, J., Zoufal, C., & Woerner, S. (2019). + Iterative Quantum Amplitude Estimation. + `arXiv:1912.05559 `_. + [2]: Brassard, G., Hoyer, P., Mosca, M., & Tapp, A. (2000). + Quantum Amplitude Amplification and Estimation. + `arXiv:quant-ph/0005055 `_. + """ + + def __init__(self, + epsilon_target: float, + alpha: float, + confint_method: str = 'beta', + min_ratio: float = 2, + quantum_instance: Optional[Union[QuantumInstance, BaseBackend, Backend]] = None + ) -> None: + r""" + The output of the algorithm is an estimate for the amplitude `a`, that with at least + probability 1 - alpha has an error of epsilon. The number of A operator calls scales + linearly in 1/epsilon (up to a logarithmic factor). + + Args: + epsilon_target: Target precision for estimation target `a`, has values between 0 and 0.5 + alpha: Confidence level, the target probability is 1 - alpha, has values between 0 and 1 + confint_method: Statistical method used to estimate the confidence intervals in + each iteration, can be 'chernoff' for the Chernoff intervals or 'beta' for the + Clopper-Pearson intervals (default) + min_ratio: Minimal q-ratio (:math:`K_{i+1} / K_i`) for FindNextK + quantum_instance: Quantum Instance or Backend + + Raises: + AlgorithmError: if the method to compute the confidence intervals is not supported + ValueError: If the target epsilon is not in (0, 0.5] + ValueError: If alpha is not in (0, 1) + ValueError: If confint_method is not supported + """ + # validate ranges of input arguments + if not 0 < epsilon_target <= 0.5: + raise ValueError(f'The target epsilon must be in (0, 0.5], but is {epsilon_target}.') + + if not 0 < alpha < 1: + raise ValueError(f'The confidence level alpha must be in (0, 1), but is {alpha}') + + if confint_method not in {'chernoff', 'beta'}: + raise ValueError('The confidence interval method must be chernoff or beta, but ' + f'is {confint_method}.') + + super().__init__() + + # set quantum instance + self.quantum_instance = quantum_instance + + # store parameters + self._epsilon = epsilon_target + self._alpha = alpha + self._min_ratio = min_ratio + self._confint_method = confint_method + + @property + def quantum_instance(self) -> Optional[QuantumInstance]: + """Get the quantum instance. + + Returns: + The quantum instance used to run this algorithm. + """ + return self._quantum_instance + + @quantum_instance.setter + def quantum_instance(self, quantum_instance: Union[QuantumInstance, + BaseBackend, Backend]) -> None: + """Set quantum instance. + + Args: + quantum_instance: The quantum instance used to run this algorithm. + """ + if isinstance(quantum_instance, (BaseBackend, Backend)): + quantum_instance = QuantumInstance(quantum_instance) + self._quantum_instance = quantum_instance + + @property + def epsilon_target(self) -> float: + """Returns the target precision ``epsilon_target`` of the algorithm. + + Returns: + The target precision (which is half the width of the confidence interval). + """ + return self._epsilon + + @epsilon_target.setter + def epsilon_target(self, epsilon: float) -> None: + """Set the target precision of the algorithm. + + Args: + epsilon: Target precision for estimation target `a`. + """ + self._epsilon = epsilon + + def _find_next_k(self, k: int, upper_half_circle: bool, theta_interval: Tuple[float, float], + min_ratio: float = 2.0) -> Tuple[int, bool]: + """Find the largest integer k_next, such that the interval (4 * k_next + 2)*theta_interval + lies completely in [0, pi] or [pi, 2pi], for theta_interval = (theta_lower, theta_upper). + + Args: + k: The current power of the Q operator. + upper_half_circle: Boolean flag of whether theta_interval lies in the + upper half-circle [0, pi] or in the lower one [pi, 2pi]. + theta_interval: The current confidence interval for the angle theta, + i.e. (theta_lower, theta_upper). + min_ratio: Minimal ratio K/K_next allowed in the algorithm. + + Returns: + The next power k, and boolean flag for the extrapolated interval. + + Raises: + AlgorithmError: if min_ratio is smaller or equal to 1 + """ + if min_ratio <= 1: + raise AlgorithmError('min_ratio must be larger than 1 to ensure convergence') + + # initialize variables + theta_l, theta_u = theta_interval + old_scaling = 4 * k + 2 # current scaling factor, called K := (4k + 2) + + # the largest feasible scaling factor K cannot be larger than K_max, + # which is bounded by the length of the current confidence interval + max_scaling = int(1 / (2 * (theta_u - theta_l))) + scaling = max_scaling - (max_scaling - 2) % 4 # bring into the form 4 * k_max + 2 + + # find the largest feasible scaling factor K_next, and thus k_next + while scaling >= min_ratio * old_scaling: + theta_min = scaling * theta_l - int(scaling * theta_l) + theta_max = scaling * theta_u - int(scaling * theta_u) + + if theta_min <= theta_max <= 0.5 and theta_min <= 0.5: + # the extrapolated theta interval is in the upper half-circle + upper_half_circle = True + return int((scaling - 2) / 4), upper_half_circle + + elif theta_max >= 0.5 and theta_max >= theta_min >= 0.5: + # the extrapolated theta interval is in the upper half-circle + upper_half_circle = False + return int((scaling - 2) / 4), upper_half_circle + + scaling -= 4 + + # if we do not find a feasible k, return the old one + return int(k), upper_half_circle + + def construct_circuit(self, estimation_problem: EstimationProblem, + k: int = 0, measurement: bool = False) -> QuantumCircuit: + r"""Construct the circuit :math:`\mathcal{Q}^k \mathcal{A} |0\rangle`. + + The A operator is the unitary specifying the QAE problem and Q the associated Grover + operator. + + Args: + estimation_problem: The estimation problem for which to construct the QAE circuit. + k: The power of the Q operator. + measurement: Boolean flag to indicate if measurements should be included in the + circuits. + + Returns: + The circuit implementing :math:`\mathcal{Q}^k \mathcal{A} |0\rangle`. + """ + num_qubits = max(estimation_problem.state_preparation.num_qubits, + estimation_problem.grover_operator.num_qubits) + circuit = QuantumCircuit(num_qubits, name='circuit') + + # add classical register if needed + if measurement: + c = ClassicalRegister(len(estimation_problem.objective_qubits)) + circuit.add_register(c) + + # add A operator + circuit.compose(estimation_problem.state_preparation, inplace=True) + + # add Q^k + if k != 0: + circuit.compose(estimation_problem.grover_operator.power(k), inplace=True) + + # add optional measurement + if measurement: + # real hardware can currently not handle operations after measurements, which might + # happen if the circuit gets transpiled, hence we're adding a safeguard-barrier + circuit.barrier() + circuit.measure(estimation_problem.objective_qubits, c[:]) + + return circuit + + def _good_state_probability(self, + problem: EstimationProblem, + counts_or_statevector: Union[Dict[str, int], np.ndarray], + num_state_qubits: int, + ) -> Union[Tuple[int, float], float]: + """Get the probability to measure '1' in the last qubit. + + Args: + problem: The estimation problem, used to obtain the number of objective qubits and + the ``is_good_state`` function. + counts_or_statevector: Either a counts-dictionary (with one measured qubit only!) or + the statevector returned from the statevector_simulator. + num_state_qubits: The number of state qubits. + + Returns: + If a dict is given, return (#one-counts, #one-counts/#all-counts), + otherwise Pr(measure '1' in the last qubit). + """ + if isinstance(counts_or_statevector, dict): + one_counts = 0 + for state, counts in counts_or_statevector.items(): + if problem.is_good_state(state): + one_counts += counts + + return int(one_counts), one_counts / sum(counts_or_statevector.values()) + else: + statevector = counts_or_statevector + num_qubits = int(np.log2(len(statevector))) # the total number of qubits + + # sum over all amplitudes where the objective qubit is 1 + prob = 0 + for i, amplitude in enumerate(statevector): + # consider only state qubits and revert bit order + bitstr = bin(i)[2:].zfill(num_qubits)[-num_state_qubits:][::-1] + objectives = [bitstr[index] for index in problem.objective_qubits] + if problem.is_good_state(objectives): + prob = prob + np.abs(amplitude)**2 + + return prob + + def estimate(self, estimation_problem: EstimationProblem + ) -> 'IterativeAmplitudeEstimationResult': + # initialize memory variables + powers = [0] # list of powers k: Q^k, (called 'k' in paper) + ratios = [] # list of multiplication factors (called 'q' in paper) + theta_intervals = [[0, 1 / 4]] # a priori knowledge of theta / 2 / pi + a_intervals = [[0.0, 1.0]] # a priori knowledge of the confidence interval of the estimate + num_oracle_queries = 0 + num_one_shots = [] + + # maximum number of rounds + max_rounds = int(np.log(self._min_ratio * np.pi / 8 + / self._epsilon) / np.log(self._min_ratio)) + 1 + upper_half_circle = True # initially theta is in the upper half-circle + + # for statevector we can directly return the probability to measure 1 + # note, that no iterations here are necessary + if self._quantum_instance.is_statevector: + # simulate circuit + circuit = self.construct_circuit(estimation_problem, k=0, measurement=False) + ret = self._quantum_instance.execute(circuit) + + # get statevector + statevector = ret.get_statevector(circuit) + + # calculate the probability of measuring '1' + num_qubits = circuit.num_qubits - circuit.num_ancillas + prob = self._good_state_probability(estimation_problem, statevector, num_qubits) + prob = cast(float, prob) # tell MyPy it's a float and not Tuple[int, float ] + + a_confidence_interval = [prob, prob] # type: List[float] + a_intervals.append(a_confidence_interval) + + theta_i_interval = [np.arccos(1 - 2 * a_i) / 2 / np.pi # type: ignore + for a_i in a_confidence_interval] + theta_intervals.append(theta_i_interval) + num_oracle_queries = 0 # no Q-oracle call, only a single one to A + + else: + num_iterations = 0 # keep track of the number of iterations + shots = self._quantum_instance._run_config.shots # number of shots per iteration + + # do while loop, keep in mind that we scaled theta mod 2pi such that it lies in [0,1] + while theta_intervals[-1][1] - theta_intervals[-1][0] > self._epsilon / np.pi: + num_iterations += 1 + + # get the next k + k, upper_half_circle = self._find_next_k(powers[-1], upper_half_circle, + theta_intervals[-1], # type: ignore + min_ratio=self._min_ratio) + + # store the variables + powers.append(k) + ratios.append((2 * powers[-1] + 1) / (2 * powers[-2] + 1)) + + # run measurements for Q^k A|0> circuit + circuit = self.construct_circuit(estimation_problem, k, measurement=True) + ret = self._quantum_instance.execute(circuit) + + # get the counts and store them + counts = ret.get_counts(circuit) + + # calculate the probability of measuring '1', 'prob' is a_i in the paper + num_qubits = circuit.num_qubits - circuit.num_ancillas + # type: ignore + one_counts, prob = self._good_state_probability(estimation_problem, counts, + num_qubits) + + num_one_shots.append(one_counts) + + # track number of Q-oracle calls + num_oracle_queries += shots * k + + # if on the previous iterations we have K_{i-1} == K_i, we sum these samples up + j = 1 # number of times we stayed fixed at the same K + round_shots = shots + round_one_counts = one_counts + if num_iterations > 1: + while powers[num_iterations - j] == powers[num_iterations] \ + and num_iterations >= j + 1: + j = j + 1 + round_shots += shots + round_one_counts += num_one_shots[-j] + + # compute a_min_i, a_max_i + if self._confint_method == 'chernoff': + a_i_min, a_i_max = _chernoff_confint(prob, round_shots, max_rounds, + self._alpha) + else: # 'beta' + a_i_min, a_i_max = _clopper_pearson_confint(round_one_counts, round_shots, + self._alpha / max_rounds) + + # compute theta_min_i, theta_max_i + if upper_half_circle: + theta_min_i = np.arccos(1 - 2 * a_i_min) / 2 / np.pi + theta_max_i = np.arccos(1 - 2 * a_i_max) / 2 / np.pi + else: + theta_min_i = 1 - np.arccos(1 - 2 * a_i_max) / 2 / np.pi + theta_max_i = 1 - np.arccos(1 - 2 * a_i_min) / 2 / np.pi + + # compute theta_u, theta_l of this iteration + scaling = 4 * k + 2 # current K_i factor + theta_u = (int(scaling * theta_intervals[-1][1]) + theta_max_i) / scaling + theta_l = (int(scaling * theta_intervals[-1][0]) + theta_min_i) / scaling + theta_intervals.append([theta_l, theta_u]) + + # compute a_u_i, a_l_i + a_u = np.sin(2 * np.pi * theta_u)**2 + a_l = np.sin(2 * np.pi * theta_l)**2 + a_u = cast(float, a_u) + a_l = cast(float, a_l) + a_intervals.append([a_l, a_u]) + + # get the latest confidence interval for the estimate of a + confidence_interval = tuple(a_intervals[-1]) + + # the final estimate is the mean of the confidence interval + estimation = np.mean(confidence_interval) + + result = IterativeAmplitudeEstimationResult() + result.alpha = self._alpha + result.post_processing = estimation_problem.post_processing + result.num_oracle_queries = num_oracle_queries + + result.estimation = estimation + result.epsilon_estimated = (confidence_interval[1] - confidence_interval[0]) / 2 + result.confidence_interval = confidence_interval + + result.estimation_processed = estimation_problem.post_processing(estimation) + confidence_interval = tuple(estimation_problem.post_processing(x) + for x in confidence_interval) + result.confidence_interval_processed = confidence_interval + result.epsilon_estimated_processed = (confidence_interval[1] - confidence_interval[0]) / 2 + result.estimate_intervals = a_intervals + result.theta_intervals = theta_intervals + result.powers = powers + result.ratios = ratios + + return result + + +class IterativeAmplitudeEstimationResult(AmplitudeEstimatorResult): + """The ``IterativeAmplitudeEstimation`` result object.""" + + def __init__(self) -> None: + super().__init__() + self._alpha = None + self._epsilon_target = None + self._epsilon_estimated = None + self._epsilon_estimated_processed = None + self._estimate_intervals = None + self._theta_intervals = None + self._powers = None + self._ratios = None + self._confidence_interval_processed = None + + @property + def alpha(self) -> float: + r"""Return the confidence level :math:`\alpha`.""" + return self._alpha + + @alpha.setter + def alpha(self, value: float) -> None: + r"""Set the confidence level :math:`\alpha`.""" + self._alpha = value + + @property + def epsilon_target(self) -> float: + """Return the target half-width of the confidence interval.""" + return self._epsilon_target + + @epsilon_target.setter + def epsilon_target(self, value: float) -> None: + """Set the target half-width of the confidence interval.""" + self._epsilon_target = value + + @property + def epsilon_estimated(self) -> float: + """Return the estimated half-width of the confidence interval.""" + return self._epsilon_estimated + + @epsilon_estimated.setter + def epsilon_estimated(self, value: float) -> None: + """Set the estimated half-width of the confidence interval.""" + self._epsilon_estimated = value + + @property + def epsilon_estimated_processed(self) -> float: + """Return the post-processed estimated half-width of the confidence interval.""" + return self._epsilon_estimated_processed + + @epsilon_estimated_processed.setter + def epsilon_estimated_processed(self, value: float) -> None: + """Set the post-processed estimated half-width of the confidence interval.""" + self._epsilon_estimated_processed = value + + @property + def estimate_intervals(self) -> List[List[float]]: + """Return the confidence intervals for the estimate in each iteration.""" + return self._estimate_intervals + + @estimate_intervals.setter + def estimate_intervals(self, value: List[List[float]]) -> None: + """Set the confidence intervals for the estimate in each iteration.""" + self._estimate_intervals = value + + @property + def theta_intervals(self) -> List[List[float]]: + """Return the confidence intervals for the angles in each iteration.""" + return self._theta_intervals + + @theta_intervals.setter + def theta_intervals(self, value: List[List[float]]) -> None: + """Set the confidence intervals for the angles in each iteration.""" + self._theta_intervals = value + + @property + def powers(self) -> List[int]: + """Return the powers of the Grover operator in each iteration.""" + return self._powers + + @powers.setter + def powers(self, value: List[int]) -> None: + """Set the powers of the Grover operator in each iteration.""" + self._powers = value + + @property + def ratios(self) -> List[float]: + r"""Return the ratios :math:`K_{i+1}/K_{i}` for each iteration :math:`i`.""" + return self._ratios + + @ratios.setter + def ratios(self, value: List[float]) -> None: + r"""Set the ratios :math:`K_{i+1}/K_{i}` for each iteration :math:`i`.""" + self._ratios = value + + @property + def confidence_interval_processed(self) -> Tuple[float, float]: + """Return the post-processed confidence interval.""" + return self._confidence_interval_processed + + @confidence_interval_processed.setter + def confidence_interval_processed(self, value: Tuple[float, float]) -> None: + """Set the post-processed confidence interval.""" + self._confidence_interval_processed = value + + +def _chernoff_confint(value: float, shots: int, max_rounds: int, alpha: float + ) -> Tuple[float, float]: + """Compute the Chernoff confidence interval for `shots` i.i.d. Bernoulli trials. + + The confidence interval is + + [value - eps, value + eps], where eps = sqrt(3 * log(2 * max_rounds/ alpha) / shots) + + but at most [0, 1]. + + Args: + value: The current estimate. + shots: The number of shots. + max_rounds: The maximum number of rounds, used to compute epsilon_a. + alpha: The confidence level, used to compute epsilon_a. + + Returns: + The Chernoff confidence interval. + """ + eps = np.sqrt(3 * np.log(2 * max_rounds / alpha) / shots) + lower = np.maximum(0, value - eps) + upper = np.minimum(1, value + eps) + return lower, upper + + +def _clopper_pearson_confint(counts: int, shots: int, alpha: float) -> Tuple[float, float]: + """Compute the Clopper-Pearson confidence interval for `shots` i.i.d. Bernoulli trials. + + Args: + counts: The number of positive counts. + shots: The number of shots. + alpha: The confidence level for the confidence interval. + + Returns: + The Clopper-Pearson confidence interval. + """ + lower, upper = 0, 1 + + # if counts == 0, the beta quantile returns nan + if counts != 0: + lower = beta.ppf(alpha / 2, counts, shots - counts + 1) + + # if counts == shots, the beta quantile returns nan + if counts != shots: + upper = beta.ppf(1 - alpha / 2, counts + 1, shots - counts) + + return lower, upper diff --git a/qiskit/algorithms/amplitude_estimators/mlae.py b/qiskit/algorithms/amplitude_estimators/mlae.py new file mode 100644 index 000000000000..20dd258c4634 --- /dev/null +++ b/qiskit/algorithms/amplitude_estimators/mlae.py @@ -0,0 +1,564 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2018, 2020. +# +# 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 Maximum Likelihood Amplitude Estimation algorithm.""" + +from typing import Optional, List, Union, Tuple, Dict, Callable +import logging +import numpy as np +from scipy.optimize import brute +from scipy.stats import norm, chi2 + +from qiskit.providers import BaseBackend +from qiskit.providers import Backend +from qiskit import ClassicalRegister, QuantumRegister, QuantumCircuit +from qiskit.utils import QuantumInstance + +from .amplitude_estimator import AmplitudeEstimator, AmplitudeEstimatorResult +from .estimation_problem import EstimationProblem +from ..exceptions import AlgorithmError + +logger = logging.getLogger(__name__) + +MINIMIZER = Callable[[Callable[[float], float], List[Tuple[float, float]]], float] + + +class MaximumLikelihoodAmplitudeEstimation(AmplitudeEstimator): + """The Maximum Likelihood Amplitude Estimation algorithm. + + This class implements the quantum amplitude estimation (QAE) algorithm without phase + estimation, as introduced in [1]. In comparison to the original QAE algorithm [2], + this implementation relies solely on different powers of the Grover operator and does not + require additional evaluation qubits. + Finally, the estimate is determined via a maximum likelihood estimation, which is why this + class in named ``MaximumLikelihoodAmplitudeEstimation``. + + References: + [1]: Suzuki, Y., Uno, S., Raymond, R., Tanaka, T., Onodera, T., & Yamamoto, N. (2019). + Amplitude Estimation without Phase Estimation. + `arXiv:1904.10246 `_. + [2]: Brassard, G., Hoyer, P., Mosca, M., & Tapp, A. (2000). + Quantum Amplitude Amplification and Estimation. + `arXiv:quant-ph/0005055 `_. + """ + + def __init__(self, evaluation_schedule: Union[List[int], int], + minimizer: Optional[MINIMIZER] = None, + quantum_instance: Optional[ + Union[QuantumInstance, BaseBackend, Backend]] = None) -> None: + r""" + Args: + evaluation_schedule: If a list, the powers applied to the Grover operator. The list + element must be non-negative. If a non-negative integer, an exponential schedule is + used where the highest power is 2 to the integer minus 1: + `[id, Q^2^0, ..., Q^2^(evaluation_schedule-1)]`. + minimizer: A minimizer used to find the minimum of the likelihood function. + Defaults to a brute search where the number of evaluation points is determined + according to ``evaluation_schedule``. The minimizer takes a function as first + argument and a list of (float, float) tuples (as bounds) as second argument and + returns a single float which is the found minimum. + quantum_instance: Quantum Instance or Backend + + Raises: + ValueError: If the number of oracle circuits is smaller than 1. + """ + + super().__init__() + + # set quantum instance + self.quantum_instance = quantum_instance + + # get parameters + if isinstance(evaluation_schedule, int): + if evaluation_schedule < 0: + raise ValueError('The evaluation schedule cannot be < 0.') + + self._evaluation_schedule = [0] + [2**j for j in range(evaluation_schedule)] + else: + if any([value < 0 for value in evaluation_schedule]): + raise ValueError('The elements of the evaluation schedule cannot be < 0.') + + self._evaluation_schedule = evaluation_schedule + + if minimizer is None: + # default number of evaluations is max(10^5, pi/2 * 10^3 * 2^(m)) + nevals = max(10000, int(np.pi / 2 * 1000 * 2 * self._evaluation_schedule[-1])) + + def default_minimizer(objective_fn, bounds): + return brute(objective_fn, bounds, Ns=nevals)[0] + + self._minimizer = default_minimizer + else: + self._minimizer = minimizer + + @property + def quantum_instance(self) -> Optional[QuantumInstance]: + """Get the quantum instance. + + Returns: + The quantum instance used to run this algorithm. + """ + return self._quantum_instance + + @quantum_instance.setter + def quantum_instance(self, quantum_instance: Union[QuantumInstance, + BaseBackend, Backend]) -> None: + """Set quantum instance. + + Args: + quantum_instance: The quantum instance used to run this algorithm. + """ + if isinstance(quantum_instance, (BaseBackend, Backend)): + quantum_instance = QuantumInstance(quantum_instance) + self._quantum_instance = quantum_instance + + def construct_circuits(self, estimation_problem: EstimationProblem, + measurement: bool = False) -> List[QuantumCircuit]: + """Construct the Amplitude Estimation w/o QPE quantum circuits. + + Args: + estimation_problem: The estimation problem for which to construct the QAE circuit. + measurement: Boolean flag to indicate if measurement should be included in the circuits. + + Returns: + A list with the QuantumCircuit objects for the algorithm. + """ + # keep track of the Q-oracle queries + circuits = [] + + num_qubits = max(estimation_problem.state_preparation.num_qubits, + estimation_problem.grover_operator.num_qubits) + q = QuantumRegister(num_qubits, 'q') + qc_0 = QuantumCircuit(q, name='qc_a') # 0 applications of Q, only a single A operator + + # add classical register if needed + if measurement: + c = ClassicalRegister(len(estimation_problem.objective_qubits)) + qc_0.add_register(c) + + qc_0.compose(estimation_problem.state_preparation, inplace=True) + + for k in self._evaluation_schedule: + qc_k = qc_0.copy(name='qc_a_q_%s' % k) + + if k != 0: + qc_k.compose(estimation_problem.grover_operator.power(k), inplace=True) + + if measurement: + # real hardware can currently not handle operations after measurements, + # which might happen if the circuit gets transpiled, hence we're adding + # a safeguard-barrier + qc_k.barrier() + qc_k.measure(estimation_problem.objective_qubits, c[:]) + + circuits += [qc_k] + + return circuits + + @staticmethod + def compute_confidence_interval(result: 'MaximumLikelihoodAmplitudeEstimationResult', + alpha: float, + kind: str = 'fisher', + apply_post_processing: bool = False) -> Tuple[float, float]: + # pylint: disable=wrong-spelling-in-docstring + """Compute the `alpha` confidence interval using the method `kind`. + + The confidence level is (1 - `alpha`) and supported kinds are 'fisher', + 'likelihood_ratio' and 'observed_fisher' with shorthand + notations 'fi', 'lr' and 'oi', respectively. + + Args: + result: A maximum likelihood amplitude estimation result. + alpha: The confidence level. + kind: The method to compute the confidence interval. Defaults to 'fisher', which + computes the theoretical Fisher information. + apply_post_processing: If True, apply post-processing to the confidence interval. + + Returns: + The specified confidence interval. + + Raises: + AlgorithmError: If `run()` hasn't been called yet. + NotImplementedError: If the method `kind` is not supported. + """ + interval = None + + # if statevector simulator the estimate is exact + if all(isinstance(data, (list, np.ndarray)) for data in result.circuit_results): + interval = 2 * [result.estimation] + + elif kind in ['likelihood_ratio', 'lr']: + interval = _likelihood_ratio_confint(result, alpha) + + elif kind in ['fisher', 'fi']: + interval = _fisher_confint(result, alpha, observed=False) + + elif kind in ['observed_fisher', 'observed_information', 'oi']: + interval = _fisher_confint(result, alpha, observed=True) + + if interval is None: + raise NotImplementedError('CI `{}` is not implemented.'.format(kind)) + + if apply_post_processing: + return tuple(result.post_processing(value) for value in interval) + + return interval + + def compute_mle(self, + circuit_results: Union[List[Dict[str, int]], List[np.ndarray]], + estimation_problem: EstimationProblem, + num_state_qubits: Optional[int] = None, + return_counts: bool = False) -> Union[float, Tuple[float, List[float]]]: + """Compute the MLE via a grid-search. + + This is a stable approach if sufficient gridpoints are used. + + Args: + circuit_results: A list of circuit outcomes. Can be counts or statevectors. + estimation_problem: The estimation problem containing the evaluation schedule and the + number of likelihood function evaluations used to find the minimum. + num_state_qubits: The number of state qubits, required for statevector simulations. + return_counts: If True, returns the good counts. + + Returns: + The MLE for the provided result object. + """ + good_counts, all_counts = _get_counts(circuit_results, estimation_problem, num_state_qubits) + + # search range + eps = 1e-15 # to avoid invalid value in log + search_range = [0 + eps, np.pi / 2 - eps] + + def loglikelihood(theta): + # loglik contains the first `it` terms of the full loglikelihood + loglik = 0 + for i, k in enumerate(self._evaluation_schedule): + angle = (2 * k + 1) * theta + loglik += np.log(np.sin(angle) ** 2) * good_counts[i] + loglik += np.log(np.cos(angle) ** 2) * (all_counts[i] - good_counts[i]) + return -loglik + + est_theta = self._minimizer(loglikelihood, [search_range]) + + if return_counts: + return est_theta, good_counts + return est_theta + + def estimate(self, estimation_problem: EstimationProblem + ) -> 'MaximumLikelihoodAmplitudeEstimationResult': + if estimation_problem.state_preparation is None: + raise AlgorithmError('Either the state_preparation variable or the a_factory ' + '(deprecated) must be set to run the algorithm.') + + result = MaximumLikelihoodAmplitudeEstimationResult() + result.evaluation_schedule = self._evaluation_schedule + result.minimizer = self._minimizer + result.evaluation_schedule = self._evaluation_schedule + result.post_processing = estimation_problem.post_processing + + if self._quantum_instance.is_statevector: + # run circuit on statevector simulator + circuits = self.construct_circuits(estimation_problem, measurement=False) + ret = self._quantum_instance.execute(circuits) + + # get statevectors and construct MLE input + statevectors = [np.asarray(ret.get_statevector(circuit)) for circuit in circuits] + result.circuit_results = statevectors + + # to count the number of Q-oracle calls (don't count shots) + result.shots = 1 + + else: + # run circuit on QASM simulator + circuits = self.construct_circuits(estimation_problem, measurement=True) + ret = self._quantum_instance.execute(circuits) + + # get counts and construct MLE input + result.circuit_results = [ret.get_counts(circuit) for circuit in circuits] + + # to count the number of Q-oracle calls + result.shots = self._quantum_instance._run_config.shots + + # run maximum likelihood estimation + num_state_qubits = circuits[0].num_qubits - circuits[0].num_ancillas + theta, good_counts = self.compute_mle(result.circuit_results, estimation_problem, + num_state_qubits, True) + + # store results + result.theta = theta + result.good_counts = good_counts + result.estimation = np.sin(result.theta)**2 + + # not sure why pylint complains, this is a callable and the tests pass + # pylint: disable=not-callable + result.estimation_processed = result.post_processing(result.estimation) + + result.fisher_information = _compute_fisher_information(result) + result.num_oracle_queries = result.shots * sum(k for k in result.evaluation_schedule) + + # compute and store confidence interval + confidence_interval = self.compute_confidence_interval(result, alpha=0.05, kind='fisher') + result.confidence_interval = confidence_interval + result.confidence_interval_processed = tuple(estimation_problem.post_processing(value) + for value in confidence_interval) + + return result + + +class MaximumLikelihoodAmplitudeEstimationResult(AmplitudeEstimatorResult): + """The ``MaximumLikelihoodAmplitudeEstimation`` result object.""" + + def __init__(self) -> None: + super().__init__() + self._theta = None + self._minimizer = None + self._good_counts = None + self._evaluation_schedule = None + self._fisher_information = None + + @property + def theta(self) -> float: + r"""Return the estimate for the angle :math:`\theta`.""" + return self._theta + + @theta.setter + def theta(self, value: float) -> None: + r"""Set the estimate for the angle :math:`\theta`.""" + self._theta = value + + @property + def minimizer(self) -> callable: + """Return the minimizer used for the search of the likelihood function.""" + return self._minimizer + + @minimizer.setter + def minimizer(self, value: callable) -> None: + """Set the number minimizer used for the search of the likelihood function.""" + self._minimizer = value + + @property + def good_counts(self) -> List[float]: + """Return the percentage of good counts per circuit power.""" + return self._good_counts + + @good_counts.setter + def good_counts(self, counts: List[float]) -> None: + """Set the percentage of good counts per circuit power.""" + self._good_counts = counts + + @property + def evaluation_schedule(self) -> List[int]: + """Return the evaluation schedule for the powers of the Grover operator.""" + return self._evaluation_schedule + + @evaluation_schedule.setter + def evaluation_schedule(self, evaluation_schedule: List[int]) -> None: + """Set the evaluation schedule for the powers of the Grover operator.""" + self._evaluation_schedule = evaluation_schedule + + @property + def fisher_information(self) -> float: + """Return the Fisher information for the estimated amplitude.""" + return self._fisher_information + + @fisher_information.setter + def fisher_information(self, value: float) -> None: + """Set the Fisher information for the estimated amplitude.""" + self._fisher_information = value + + +def _safe_min(array, default=0): + if len(array) == 0: + return default + return np.min(array) + + +def _safe_max(array, default=(np.pi / 2)): + if len(array) == 0: + return default + return np.max(array) + + +def _compute_fisher_information(result: 'MaximumLikelihoodAmplitudeEstimationResult', + num_sum_terms: Optional[int] = None, + observed: bool = False) -> float: + """Compute the Fisher information. + + Args: + result: A maximum likelihood amplitude estimation result. + num_sum_terms: The number of sum terms to be included in the calculation of the + Fisher information. By default all values are included. + observed: If True, compute the observed Fisher information, otherwise the theoretical + one. + + Returns: + The computed Fisher information, or np.inf if statevector simulation was used. + + Raises: + KeyError: Call run() first! + """ + a = result.estimation + + # Corresponding angle to the value a (only use real part of 'a') + theta_a = np.arcsin(np.sqrt(np.real(a))) + + # Get the number of hits (shots_k) and one-hits (h_k) + one_hits = result.good_counts + all_hits = [result.shots] * len(one_hits) + + # Include all sum terms or just up to a certain term? + evaluation_schedule = result.evaluation_schedule + if num_sum_terms is not None: + evaluation_schedule = evaluation_schedule[:num_sum_terms] + # not necessary since zip goes as far as shortest list: + # all_hits = all_hits[:num_sum_terms] + # one_hits = one_hits[:num_sum_terms] + + # Compute the Fisher information + fisher_information = None + if observed: + # Note, that the observed Fisher information is very unreliable in this algorithm! + d_loglik = 0 + for shots_k, h_k, m_k in zip(all_hits, one_hits, evaluation_schedule): + tan = np.tan((2 * m_k + 1) * theta_a) + d_loglik += (2 * m_k + 1) * (h_k / tan + (shots_k - h_k) * tan) + + d_loglik /= np.sqrt(a * (1 - a)) + fisher_information = d_loglik ** 2 / len(all_hits) + + else: + fisher_information = sum(shots_k * (2 * m_k + 1)**2 + for shots_k, m_k in zip(all_hits, evaluation_schedule)) + fisher_information /= a * (1 - a) + + return fisher_information + + +def _fisher_confint(result: MaximumLikelihoodAmplitudeEstimationResult, + alpha: float = 0.05, + observed: bool = False) -> Tuple[float, float]: + """Compute the `alpha` confidence interval based on the Fisher information. + + Args: + result: A maximum likelihood amplitude estimation results object. + alpha: The level of the confidence interval (must be <= 0.5), default to 0.05. + observed: If True, use observed Fisher information. + + Returns: + float: The alpha confidence interval based on the Fisher information + Raises: + AssertionError: Call run() first! + """ + # Get the (observed) Fisher information + fisher_information = None + try: + fisher_information = result.fisher_information + except KeyError as ex: + raise AssertionError("Call run() first!") from ex + + if observed: + fisher_information = _compute_fisher_information(result, observed=True) + + normal_quantile = norm.ppf(1 - alpha / 2) + confint = np.real(result.estimation) + \ + normal_quantile / np.sqrt(fisher_information) * np.array([-1, 1]) + mapped_confint = tuple(result.post_processing(bound) for bound in confint) + return mapped_confint + + +def _likelihood_ratio_confint(result: MaximumLikelihoodAmplitudeEstimationResult, + alpha: float = 0.05, + nevals: Optional[int] = None) -> List[float]: + """Compute the likelihood-ratio confidence interval. + + Args: + result: A maximum likelihood amplitude estimation results object. + alpha: The level of the confidence interval (< 0.5), defaults to 0.05. + nevals: The number of evaluations to find the intersection with the loglikelihood + function. Defaults to an adaptive value based on the maximal power of Q. + + Returns: + The alpha-likelihood-ratio confidence interval. + """ + if nevals is None: + nevals = max(10000, int(np.pi / 2 * 1000 * 2 * result.evaluation_schedule[-1])) + + def loglikelihood(theta, one_counts, all_counts): + loglik = 0 + for i, k in enumerate(result.evaluation_schedule): + loglik += np.log(np.sin((2 * k + 1) * theta) ** 2) * one_counts[i] + loglik += np.log(np.cos((2 * k + 1) * theta) ** 2) * (all_counts[i] - one_counts[i]) + return loglik + + one_counts = result.good_counts + all_counts = [result.shots] * len(one_counts) + + eps = 1e-15 # to avoid invalid value in log + thetas = np.linspace(0 + eps, np.pi / 2 - eps, nevals) + values = np.zeros(len(thetas)) + for i, theta in enumerate(thetas): + values[i] = loglikelihood(theta, one_counts, all_counts) + + loglik_mle = loglikelihood(result.theta, one_counts, all_counts) + chi2_quantile = chi2.ppf(1 - alpha, df=1) + thres = loglik_mle - chi2_quantile / 2 + + # the (outer) LR confidence interval + above_thres = thetas[values >= thres] + + # it might happen that the `above_thres` array is empty, + # to still provide a valid result use safe_min/max which + # then yield [0, pi/2] + confint = [_safe_min(above_thres, default=0), + _safe_max(above_thres, default=np.pi / 2)] + mapped_confint = tuple(result.post_processing(np.sin(bound) ** 2) for bound in confint) + + return mapped_confint + + +def _get_counts(circuit_results: List[Union[np.ndarray, List[float], Dict[str, int]]], + estimation_problem: EstimationProblem, + num_state_qubits: int, + ) -> Tuple[List[float], List[int]]: + """Get the good and total counts. + + Returns: + A pair of two lists, ([1-counts per experiment], [shots per experiment]). + + Raises: + AlgorithmError: If self.run() has not been called yet. + """ + one_hits = [] # h_k: how often 1 has been measured, for a power Q^(m_k) + all_hits = [] # shots_k: how often has been measured at a power Q^(m_k) + if all(isinstance(data, (list, np.ndarray)) for data in circuit_results): + probabilities = [] + num_qubits = int(np.log2(len(circuit_results[0]))) # the total number of qubits + for statevector in circuit_results: + p_k = 0.0 + for i, amplitude in enumerate(statevector): + probability = np.abs(amplitude) ** 2 + # consider only state qubits and revert bit order + bitstr = bin(i)[2:].zfill(num_qubits)[-num_state_qubits:][::-1] + objectives = [bitstr[index] for index in estimation_problem.objective_qubits] + if estimation_problem.is_good_state(objectives): + p_k += probability + probabilities += [p_k] + + one_hits = probabilities + all_hits = np.ones_like(one_hits) + else: + for counts in circuit_results: + all_hits.append(sum(counts.values())) + one_hits.append(sum(count for bitstr, count in counts.items() + if estimation_problem.is_good_state(bitstr))) + + return one_hits, all_hits diff --git a/qiskit/circuit/library/grover_operator.py b/qiskit/circuit/library/grover_operator.py index 3b79afafbba5..49e173dd3dc1 100644 --- a/qiskit/circuit/library/grover_operator.py +++ b/qiskit/circuit/library/grover_operator.py @@ -257,6 +257,9 @@ def _build(self): self.compose(self.state_preparation, list(range(self.state_preparation.num_qubits)), inplace=True) + # minus sign + self.global_phase = numpy.pi + # TODO use the oracle compiler or the bit string oracle def _zero_reflection(num_state_qubits: int, qubits: List[int], mcx_mode: Optional[str] = None diff --git a/test/python/algorithms/test_amplitude_estimators.py b/test/python/algorithms/test_amplitude_estimators.py new file mode 100644 index 000000000000..e185643e99d3 --- /dev/null +++ b/test/python/algorithms/test_amplitude_estimators.py @@ -0,0 +1,498 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2018, 2020. +# +# 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. + +"""Test the quantum amplitude estimation algorithm.""" + +import unittest +from test.python.algorithms import QiskitAlgorithmsTestCase +import numpy as np +from ddt import ddt, idata, data, unpack +from qiskit import QuantumRegister, QuantumCircuit, BasicAer +from qiskit.circuit.library import QFT, GroverOperator +from qiskit.utils import QuantumInstance +from qiskit.algorithms import ( + AmplitudeEstimation, MaximumLikelihoodAmplitudeEstimation, IterativeAmplitudeEstimation, + FasterAmplitudeEstimation, EstimationProblem +) +from qiskit.quantum_info import Operator, Statevector + + +class BernoulliStateIn(QuantumCircuit): + """A circuit preparing sqrt(1 - p)|0> + sqrt(p)|1>.""" + + def __init__(self, probability): + super().__init__(1) + angle = 2 * np.arcsin(np.sqrt(probability)) + self.ry(angle, 0) + + +class BernoulliGrover(QuantumCircuit): + """The Grover operator corresponding to the Bernoulli A operator.""" + + def __init__(self, probability): + super().__init__(1, global_phase=np.pi) + self.angle = 2 * np.arcsin(np.sqrt(probability)) + self.ry(2 * self.angle, 0) + + def power(self, power, matrix_power=False): + if matrix_power: + return super().power(power, True) + + powered = QuantumCircuit(1) + powered.ry(power * 2 * self.angle, 0) + return powered + + +class SineIntegral(QuantumCircuit): + r"""Construct the A operator to approximate the integral + + \int_0^1 \sin^2(x) d x + + with a specified number of qubits. + """ + + def __init__(self, num_qubits): + qr_state = QuantumRegister(num_qubits, 'state') + qr_objective = QuantumRegister(1, 'obj') + super().__init__(qr_state, qr_objective) + + # prepare 1/sqrt{2^n} sum_x |x>_n + self.h(qr_state) + + # apply the sine/cosine term + self.ry(2 * 1 / 2 / 2 ** num_qubits, qr_objective[0]) + for i, qubit in enumerate(qr_state): + self.cry(2 * 2 ** i / 2 ** num_qubits, qubit, qr_objective[0]) + + +@ddt +class TestBernoulli(QiskitAlgorithmsTestCase): + """Tests based on the Bernoulli A operator. + + This class tests + * the estimation result + * the constructed circuits + """ + + def setUp(self): + super().setUp() + + self._statevector = QuantumInstance(backend=BasicAer.get_backend('statevector_simulator'), + seed_simulator=2, seed_transpiler=2) + self._unitary = QuantumInstance(backend=BasicAer.get_backend('unitary_simulator'), shots=1, + seed_simulator=42, seed_transpiler=91) + + def qasm(shots=100): + return QuantumInstance(backend=BasicAer.get_backend('qasm_simulator'), shots=shots, + seed_simulator=2, seed_transpiler=2) + self._qasm = qasm + + @idata([ + [0.2, AmplitudeEstimation(2), {'estimation': 0.5, 'mle': 0.2}], + [0.49, AmplitudeEstimation(3), {'estimation': 0.5, 'mle': 0.49}], + [0.2, MaximumLikelihoodAmplitudeEstimation([0, 1, 2]), {'estimation': 0.2}], + [0.49, MaximumLikelihoodAmplitudeEstimation(3), {'estimation': 0.49}], + [0.2, IterativeAmplitudeEstimation(0.1, 0.1), {'estimation': 0.2}], + [0.49, IterativeAmplitudeEstimation(0.001, 0.01), {'estimation': 0.49}], + [0.2, FasterAmplitudeEstimation(0.1, 3, rescale=False), {'estimation': 0.2}], + [0.12, FasterAmplitudeEstimation(0.1, 2, rescale=False), {'estimation': 0.12}] + ]) + @unpack + def test_statevector(self, prob, qae, expect): + """ statevector test """ + qae.quantum_instance = self._statevector + problem = EstimationProblem(BernoulliStateIn(prob), [0], BernoulliGrover(prob)) + + result = qae.estimate(problem) + self.assertGreaterEqual(self._statevector.time_taken, 0.) + self._statevector.reset_execution_results() + for key, value in expect.items(): + self.assertAlmostEqual(value, getattr(result, key), places=3, + msg="estimate `{}` failed".format(key)) + + @idata([ + [0.2, 100, AmplitudeEstimation(4), {'estimation': 0.14644, 'mle': 0.193888}], + [0.0, 1000, AmplitudeEstimation(2), {'estimation': 0.0, 'mle': 0.0}], + [0.2, 100, MaximumLikelihoodAmplitudeEstimation([0, 1, 2, 4, 8]), {'estimation': 0.199606}], + [0.8, 10, IterativeAmplitudeEstimation(0.1, 0.05), {'estimation': 0.811711}], + [0.2, 1000, FasterAmplitudeEstimation(0.1, 3, rescale=False), {'estimation': 0.198640}], + [0.12, 100, FasterAmplitudeEstimation(0.01, 3, rescale=False), {'estimation': 0.119037}], + ]) + @ unpack + def test_qasm(self, prob, shots, qae, expect): + """ qasm test """ + qae.quantum_instance = self._qasm(shots) + problem = EstimationProblem(BernoulliStateIn(prob), [0], BernoulliGrover(prob)) + + result = qae.estimate(problem) + for key, value in expect.items(): + self.assertAlmostEqual(value, getattr(result, key), places=3, + msg="estimate `{}` failed".format(key)) + + @ data(True, False) + def test_qae_circuit(self, efficient_circuit): + """Test circuits resulting from canonical amplitude estimation. + + Build the circuit manually and from the algorithm and compare the resulting unitaries. + """ + prob = 0.5 + + problem = EstimationProblem(BernoulliStateIn(prob), objective_qubits=[0]) + for m in [2, 5]: + qae = AmplitudeEstimation(m) + angle = 2 * np.arcsin(np.sqrt(prob)) + + # manually set up the inefficient AE circuit + qr_eval = QuantumRegister(m, 'a') + qr_objective = QuantumRegister(1, 'q') + circuit = QuantumCircuit(qr_eval, qr_objective) + + # initial Hadamard gates + for i in range(m): + circuit.h(qr_eval[i]) + + # A operator + circuit.ry(angle, qr_objective) + + if efficient_circuit: + qae.grover_operator = BernoulliGrover(prob) + for power in range(m): + circuit.cry(2 * 2 ** power * angle, qr_eval[power], qr_objective[0]) + else: + oracle = QuantumCircuit(1) + oracle.z(0) + + state_preparation = QuantumCircuit(1) + state_preparation.ry(angle, 0) + grover_op = GroverOperator(oracle, state_preparation) + grover_op.global_phase = np.pi + for power in range(m): + circuit.compose(grover_op.power(2 ** power).control(), + qubits=[qr_eval[power], qr_objective[0]], + inplace=True) + + # fourier transform + iqft = QFT(m, do_swaps=False).inverse() + circuit.append(iqft.to_instruction(), qr_eval) + + actual_circuit = qae.construct_circuit(problem, measurement=False) + + self.assertEqual(Operator(circuit), Operator(actual_circuit)) + + @ data(True, False) + def test_iqae_circuits(self, efficient_circuit): + """Test circuits resulting from iterative amplitude estimation. + + Build the circuit manually and from the algorithm and compare the resulting unitaries. + """ + prob = 0.5 + problem = EstimationProblem(BernoulliStateIn(prob), objective_qubits=[0]) + + for k in [2, 5]: + qae = IterativeAmplitudeEstimation(0.01, 0.05) + angle = 2 * np.arcsin(np.sqrt(prob)) + + # manually set up the inefficient AE circuit + q_objective = QuantumRegister(1, 'q') + circuit = QuantumCircuit(q_objective) + + # A operator + circuit.ry(angle, q_objective) + + if efficient_circuit: + qae.grover_operator = BernoulliGrover(prob) + circuit.ry(2 * k * angle, q_objective[0]) + + else: + oracle = QuantumCircuit(1) + oracle.z(0) + state_preparation = QuantumCircuit(1) + state_preparation.ry(angle, 0) + grover_op = GroverOperator(oracle, state_preparation) + grover_op.global_phase = np.pi + for _ in range(k): + circuit.compose(grover_op, inplace=True) + + actual_circuit = qae.construct_circuit(problem, k, measurement=False) + self.assertEqual(Operator(circuit), Operator(actual_circuit)) + + @ data(True, False) + def test_mlae_circuits(self, efficient_circuit): + """ Test the circuits constructed for MLAE """ + prob = 0.5 + problem = EstimationProblem(BernoulliStateIn(prob), objective_qubits=[0]) + + for k in [2, 5]: + qae = MaximumLikelihoodAmplitudeEstimation(k) + angle = 2 * np.arcsin(np.sqrt(prob)) + + # compute all the circuits used for MLAE + circuits = [] + + # 0th power + q_objective = QuantumRegister(1, 'q') + circuit = QuantumCircuit(q_objective) + circuit.ry(angle, q_objective) + circuits += [circuit] + + # powers of 2 + for power in range(k): + q_objective = QuantumRegister(1, 'q') + circuit = QuantumCircuit(q_objective) + + # A operator + circuit.ry(angle, q_objective) + + # Q^(2^j) operator + if efficient_circuit: + qae.grover_operator = BernoulliGrover(prob) + circuit.ry(2 * 2 ** power * angle, q_objective[0]) + + else: + oracle = QuantumCircuit(1) + oracle.z(0) + state_preparation = QuantumCircuit(1) + state_preparation.ry(angle, 0) + grover_op = GroverOperator(oracle, state_preparation) + grover_op.global_phase = np.pi + for _ in range(2**power): + circuit.compose(grover_op, inplace=True) + + actual_circuits = qae.construct_circuits(problem, measurement=False) + + for actual, expected in zip(actual_circuits, circuits): + self.assertEqual(Operator(actual), Operator(expected)) + + +@ ddt +class TestSineIntegral(QiskitAlgorithmsTestCase): + """Tests based on the A operator to integrate sin^2(x). + + This class tests + * the estimation result + * the confidence intervals + """ + + def setUp(self): + super().setUp() + + self._statevector = QuantumInstance(backend=BasicAer.get_backend('statevector_simulator'), + seed_simulator=123, + seed_transpiler=41) + + def qasm(shots=100): + return QuantumInstance(backend=BasicAer.get_backend('qasm_simulator'), shots=shots, + seed_simulator=7192, seed_transpiler=90000) + + self._qasm = qasm + + @ idata([ + [2, AmplitudeEstimation(2), {'estimation': 0.5, 'mle': 0.270290}], + [4, MaximumLikelihoodAmplitudeEstimation(4), {'estimation': 0.272675}], + [3, IterativeAmplitudeEstimation(0.1, 0.1), {'estimation': 0.272082}], + [3, FasterAmplitudeEstimation(0.01, 1), {'estimation': 0.272082}], + ]) + @ unpack + def test_statevector(self, n, qae, expect): + """ Statevector end-to-end test """ + # construct factories for A and Q + # qae.state_preparation = SineIntegral(n) + qae.quantum_instance = self._statevector + estimation_problem = EstimationProblem(SineIntegral(n), objective_qubits=[n]) + + # result = qae.run(self._statevector) + result = qae.estimate(estimation_problem) + self.assertGreaterEqual(self._statevector.time_taken, 0.) + self._statevector.reset_execution_results() + for key, value in expect.items(): + self.assertAlmostEqual(value, getattr(result, key), places=3, + msg="estimate `{}` failed".format(key)) + + @ idata([ + [4, 10, AmplitudeEstimation(2), {'estimation': 0.5, 'mle': 0.333333}], + [3, 10, MaximumLikelihoodAmplitudeEstimation(2), {'estimation': 0.256878}], + [3, 1000, IterativeAmplitudeEstimation(0.01, 0.01), {'estimation': 0.271790}], + [3, 1000, FasterAmplitudeEstimation(0.1, 4), {'estimation': 0.274168}], + ]) + @ unpack + def test_qasm(self, n, shots, qae, expect): + """QASM simulator end-to-end test.""" + # construct factories for A and Q + qae.quantum_instance = self._qasm(shots) + estimation_problem = EstimationProblem(SineIntegral(n), objective_qubits=[n]) + + result = qae.estimate(estimation_problem) + for key, value in expect.items(): + self.assertAlmostEqual(value, getattr(result, key), places=3, + msg="estimate `{}` failed".format(key)) + + @ idata([ + [AmplitudeEstimation(3), 'mle', + {'likelihood_ratio': (0.2494734, 0.3003771), + 'fisher': (0.2486176, 0.2999286), + 'observed_fisher': (0.2484562, 0.3000900)} + ], + [MaximumLikelihoodAmplitudeEstimation(3), 'estimation', + {'likelihood_ratio': (0.2598794, 0.2798536), + 'fisher': (0.2584889, 0.2797018), + 'observed_fisher': (0.2659279, 0.2722627)}], + ]) + @ unpack + def test_confidence_intervals(self, qae, key, expect): + """End-to-end test for all confidence intervals.""" + n = 3 + qae.quantum_instance = self._statevector + estimation_problem = EstimationProblem(SineIntegral(n), objective_qubits=[n]) + + # statevector simulator + result = qae.estimate(estimation_problem) + self.assertGreater(self._statevector.time_taken, 0.) + self._statevector.reset_execution_results() + methods = ['lr', 'fi', 'oi'] # short for likelihood_ratio, fisher, observed_fisher + alphas = [0.1, 0.00001, 0.9] # alpha shouldn't matter in statevector + for alpha, method in zip(alphas, methods): + confint = qae.compute_confidence_interval(result, alpha, method) + # confidence interval based on statevector should be empty, as we are sure of the result + self.assertAlmostEqual(confint[1] - confint[0], 0.0) + self.assertAlmostEqual(confint[0], getattr(result, key)) + + # qasm simulator + shots = 100 + alpha = 0.01 + qae.quantum_instance = self._qasm(shots) + result = qae.estimate(estimation_problem) + for method, expected_confint in expect.items(): + confint = qae.compute_confidence_interval(result, alpha, method) + np.testing.assert_array_almost_equal(confint, expected_confint) + self.assertTrue(confint[0] <= getattr(result, key) <= confint[1]) + + def test_iqae_confidence_intervals(self): + """End-to-end test for the IQAE confidence interval.""" + n = 3 + qae = IterativeAmplitudeEstimation(0.1, 0.01, quantum_instance=self._statevector) + expected_confint = (0.1984050, 0.3511015) + estimation_problem = EstimationProblem(SineIntegral(n), objective_qubits=[n]) + + # statevector simulator + result = qae.estimate(estimation_problem) + self.assertGreaterEqual(self._statevector.time_taken, 0.) + self._statevector.reset_execution_results() + confint = result.confidence_interval + # confidence interval based on statevector should be empty, as we are sure of the result + self.assertAlmostEqual(confint[1] - confint[0], 0.0) + self.assertAlmostEqual(confint[0], result.estimation) + + # qasm simulator + shots = 100 + qae.quantum_instance = self._qasm(shots) + result = qae.estimate(estimation_problem) + confint = result.confidence_interval + np.testing.assert_array_almost_equal(confint, expected_confint) + self.assertTrue(confint[0] <= result.estimation <= confint[1]) + + +@ ddt +class TestFasterAmplitudeEstimation(QiskitAlgorithmsTestCase): + """Specific tests for Faster AE.""" + + def test_rescaling(self): + """Test the rescaling.""" + amplitude = 0.8 + scaling = 0.25 + circuit = QuantumCircuit(1) + circuit.ry(2 * np.arcsin(amplitude), 0) + problem = EstimationProblem(circuit, objective_qubits=[0]) + + rescaled = problem.rescale(scaling) + rescaled_amplitude = Statevector.from_instruction(rescaled.state_preparation).data[3] + + self.assertAlmostEqual(scaling * amplitude, rescaled_amplitude) + + def test_run_without_rescaling(self): + """Run Faster AE without rescaling if the amplitude is in [0, 1/4].""" + # construct estimation problem + prob = 0.11 + a_op = QuantumCircuit(1) + a_op.ry(2 * np.arcsin(np.sqrt(prob)), 0) + problem = EstimationProblem(a_op, objective_qubits=[0]) + + # construct algo without rescaling + backend = BasicAer.get_backend('statevector_simulator') + fae = FasterAmplitudeEstimation(0.1, 1, rescale=False, quantum_instance=backend) + + # run the algo + result = fae.estimate(problem) + + # assert the result is correct + self.assertAlmostEqual(result.estimation, prob) + + # assert no rescaling was used + theta = np.mean(result.theta_intervals[-1]) + value_without_scaling = np.sin(theta) ** 2 + self.assertAlmostEqual(result.estimation, value_without_scaling) + + def test_rescaling_with_custom_grover_raises(self): + """Test that the rescaling option fails if a custom Grover operator is used.""" + prob = 0.8 + a_op = BernoulliStateIn(prob) + q_op = BernoulliGrover(prob) + problem = EstimationProblem(a_op, objective_qubits=[0], grover_operator=q_op) + + # construct algo without rescaling + backend = BasicAer.get_backend('statevector_simulator') + fae = FasterAmplitudeEstimation(0.1, 1, quantum_instance=backend) + + # run the algo + with self.assertWarns(Warning): + _ = fae.estimate(problem) + + @ data( + ('statevector_simulator', 0.2), + ('qasm_simulator', 0.199440) + ) + @ unpack + def test_good_state(self, backend_str, expect): + """Test with a good state function.""" + def is_good_state(bitstr): + return bitstr[1] == '1' + + # construct the estimation problem where the second qubit is ignored + a_op = QuantumCircuit(2) + a_op.ry(2 * np.arcsin(np.sqrt(0.2)), 0) + + # oracle only affects first qubit + oracle = QuantumCircuit(2) + oracle.z(0) + + # reflect only on first qubit + q_op = GroverOperator(oracle, a_op, reflection_qubits=[0]) + + # but we measure both qubits (hence both are objective qubits) + problem = EstimationProblem(a_op, objective_qubits=[0, 1], + grover_operator=q_op, + is_good_state=is_good_state) + + # construct algo + backend = QuantumInstance(BasicAer.get_backend(backend_str), + seed_simulator=2, seed_transpiler=2) + # cannot use rescaling with a custom grover operator + fae = FasterAmplitudeEstimation(0.01, 5, rescale=False, quantum_instance=backend) + + # run the algo + result = fae.estimate(problem) + + # assert the result is correct + self.assertAlmostEqual(result.estimation, expect, places=5) + + +if __name__ == '__main__': + unittest.main()