diff --git a/qiskit/algorithms/__init__.py b/qiskit/algorithms/__init__.py index 2ef602089cf1..b9dc8b8a4ed9 100644 --- a/qiskit/algorithms/__init__.py +++ b/qiskit/algorithms/__init__.py @@ -179,19 +179,6 @@ time_evolvers.trotterization -Factorizers ------------ - -Algorithms to find factors of a number. - -.. autosummary:: - :toctree: ../stubs/ - :nosignatures: - - Shor - ShorResult - - Gradients ---------- @@ -203,18 +190,6 @@ gradients -Linear Solvers --------------- - -Algorithms to solve linear systems of equations. - -.. autosummary:: - :toctree: ../stubs/ - :nosignatures: - - linear_solvers - - Minimum Eigensolvers --------------------- @@ -343,8 +318,6 @@ EstimationProblem, ) from .eigen_solvers import NumPyEigensolver, Eigensolver, EigensolverResult, VQD, VQDResult -from .factorizers import Shor, ShorResult -from .linear_solvers import HHL, LinearSolver, NumPyLinearSolver, LinearSolverResult from .minimum_eigen_solvers import ( VQE, VQEResult, @@ -411,17 +384,11 @@ "EvolutionProblem", "TimeEvolutionResult", "TimeEvolutionProblem", - "LinearSolverResult", "Eigensolver", "EigensolverResult", - "Shor", - "ShorResult", "VQE", "VQEResult", "QAOA", - "LinearSolver", - "HHL", - "NumPyLinearSolver", "NumPyMinimumEigensolver", "MinimumEigensolver", "MinimumEigensolverResult", diff --git a/qiskit/algorithms/factorizers/__init__.py b/qiskit/algorithms/factorizers/__init__.py deleted file mode 100644 index 939614a346eb..000000000000 --- a/qiskit/algorithms/factorizers/__init__.py +++ /dev/null @@ -1,20 +0,0 @@ -# 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. - -"""Factorizers Package""" - -from .shor import Shor, ShorResult - -__all__ = [ - "Shor", - "ShorResult", -] diff --git a/qiskit/algorithms/factorizers/shor.py b/qiskit/algorithms/factorizers/shor.py deleted file mode 100644 index 67d2722eb881..000000000000 --- a/qiskit/algorithms/factorizers/shor.py +++ /dev/null @@ -1,531 +0,0 @@ -# This code is part of Qiskit. -# -# (C) Copyright IBM 2019, 2022. -# -# 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. - -"""Shor's factoring algorithm.""" - -import array -import fractions -import logging -import math -import sys -from typing import Optional, Union, List, Tuple - -import numpy as np - -from qiskit import ClassicalRegister, QuantumCircuit, QuantumRegister -from qiskit.circuit import Gate, Instruction, ParameterVector -from qiskit.circuit.library import QFT -from qiskit.providers import Backend -from qiskit.quantum_info import partial_trace -from qiskit.utils import summarize_circuits -from qiskit.utils.arithmetic import is_power -from qiskit.utils.quantum_instance import QuantumInstance -from qiskit.utils.validation import validate_min -from qiskit.utils.deprecation import deprecate_function -from ..algorithm_result import AlgorithmResult -from ..exceptions import AlgorithmError - -logger = logging.getLogger(__name__) - - -# pylint: disable=invalid-name - - -class Shor: - """The deprecated Shor's factoring algorithm. - - The Shor class is deprecated as of Qiskit Terra 0.22.0 - and will be removed no sooner than 3 months after the release date. - It is replaced by the tutorial at - `Shor `_ - - Shor's Factoring algorithm is one of the most well-known quantum algorithms and finds the - prime factors for input integer :math:`N` in polynomial time. - - Adapted from https://github.com/ttlion/ShorAlgQiskit - - See also https://arxiv.org/abs/quant-ph/0205095 - """ - - @deprecate_function( - "The Shor class is deprecated as of Qiskit Terra 0.22.0 and will be removed " - "no sooner than 3 months after the release date. It is replaced by the tutorial " - "at https://qiskit.org/textbook/ch-algorithms/shor.html", - since="0.22.0", - ) - def __init__(self, quantum_instance: Optional[Union[QuantumInstance, Backend]] = None) -> None: - """ - Args: - quantum_instance: Quantum Instance or Backend - - """ - self._quantum_instance = None - if quantum_instance: - self.quantum_instance = quantum_instance - - @property - def quantum_instance(self) -> Optional[QuantumInstance]: - """Returns quantum instance.""" - return self._quantum_instance - - @quantum_instance.setter - def quantum_instance(self, quantum_instance: Union[QuantumInstance, Backend]) -> None: - """Sets quantum instance.""" - if isinstance(quantum_instance, Backend): - quantum_instance = QuantumInstance(quantum_instance) - self._quantum_instance = quantum_instance - - @staticmethod - def _get_angles(a: int, n: int) -> np.ndarray: - """Calculates the array of angles to be used in the addition in Fourier Space.""" - bits_little_endian = (bin(int(a))[2:].zfill(n))[::-1] - - angles = np.zeros(n) - for i in range(n): - for j in range(i + 1): - k = i - j - if bits_little_endian[j] == "1": - angles[i] += pow(2, -k) - - return angles * np.pi - - @staticmethod - def _phi_add_gate(angles: Union[np.ndarray, ParameterVector]) -> Gate: - """Gate that performs addition by a in Fourier Space.""" - circuit = QuantumCircuit(len(angles), name="phi_add_a") - for i, angle in enumerate(angles): - circuit.p(angle, i) - return circuit.to_gate() - - def _double_controlled_phi_add_mod_N( - self, - angles: Union[np.ndarray, ParameterVector], - c_phi_add_N: Gate, - iphi_add_N: Gate, - qft: Gate, - iqft: Gate, - ) -> QuantumCircuit: - """Creates a circuit which implements double-controlled modular addition by a.""" - ctrl_qreg = QuantumRegister(2, "ctrl") - b_qreg = QuantumRegister(len(angles), "b") - flag_qreg = QuantumRegister(1, "flag") - - circuit = QuantumCircuit(ctrl_qreg, b_qreg, flag_qreg, name="ccphi_add_a_mod_N") - - cc_phi_add_a = self._phi_add_gate(angles).control(2) - cc_iphi_add_a = cc_phi_add_a.inverse() - - circuit.append(cc_phi_add_a, [*ctrl_qreg, *b_qreg]) - - circuit.append(iphi_add_N, b_qreg) - - circuit.append(iqft, b_qreg) - circuit.cx(b_qreg[-1], flag_qreg[0]) - circuit.append(qft, b_qreg) - - circuit.append(c_phi_add_N, [*flag_qreg, *b_qreg]) - - circuit.append(cc_iphi_add_a, [*ctrl_qreg, *b_qreg]) - - circuit.append(iqft, b_qreg) - circuit.x(b_qreg[-1]) - circuit.cx(b_qreg[-1], flag_qreg[0]) - circuit.x(b_qreg[-1]) - circuit.append(qft, b_qreg) - - circuit.append(cc_phi_add_a, [*ctrl_qreg, *b_qreg]) - - return circuit - - def _controlled_multiple_mod_N( - self, n: int, N: int, a: int, c_phi_add_N: Gate, iphi_add_N: Gate, qft: Gate, iqft: Gate - ) -> Instruction: - """Implements modular multiplication by a as an instruction.""" - ctrl_qreg = QuantumRegister(1, "ctrl") - x_qreg = QuantumRegister(n, "x") - b_qreg = QuantumRegister(n + 1, "b") - flag_qreg = QuantumRegister(1, "flag") - - circuit = QuantumCircuit(ctrl_qreg, x_qreg, b_qreg, flag_qreg, name="cmult_a_mod_N") - - angle_params = ParameterVector("angles", length=n + 1) - modulo_adder = self._double_controlled_phi_add_mod_N( - angle_params, c_phi_add_N, iphi_add_N, qft, iqft - ) - - def append_adder(adder: QuantumCircuit, constant: int, idx: int): - partial_constant = (pow(2, idx, N) * constant) % N - angles = self._get_angles(partial_constant, n + 1) - bound = adder.assign_parameters({angle_params: angles}) - circuit.append(bound, [*ctrl_qreg, x_qreg[idx], *b_qreg, *flag_qreg]) - - circuit.append(qft, b_qreg) - - # perform controlled addition by a on the aux register in Fourier space - for i in range(n): - append_adder(modulo_adder, a, i) - - circuit.append(iqft, b_qreg) - - # perform controlled subtraction by a in Fourier space on both the aux and down register - for i in range(n): - circuit.cswap(ctrl_qreg, x_qreg[i], b_qreg[i]) - - circuit.append(qft, b_qreg) - - a_inv = pow(a, -1, mod=N) if sys.version_info >= (3, 8) else self.modinv(a, N) - - modulo_adder_inv = modulo_adder.inverse() - for i in reversed(range(n)): - append_adder(modulo_adder_inv, a_inv, i) - - circuit.append(iqft, b_qreg) - - return circuit.to_instruction() - - def _power_mod_N(self, n: int, N: int, a: int) -> Instruction: - """Implements modular exponentiation a^x as an instruction.""" - up_qreg = QuantumRegister(2 * n, name="up") - down_qreg = QuantumRegister(n, name="down") - aux_qreg = QuantumRegister(n + 2, name="aux") - - circuit = QuantumCircuit(up_qreg, down_qreg, aux_qreg, name=f"{a}^x mod {N}") - - qft = QFT(n + 1, do_swaps=False).to_gate() - iqft = qft.inverse() - - # Create gates to perform addition/subtraction by N in Fourier Space - phi_add_N = self._phi_add_gate(self._get_angles(N, n + 1)) - iphi_add_N = phi_add_N.inverse() - c_phi_add_N = phi_add_N.control(1) - - # Apply the multiplication gates as showed in - # the report in order to create the exponentiation - for i in range(2 * n): - partial_a = pow(a, pow(2, i), N) - modulo_multiplier = self._controlled_multiple_mod_N( - n, N, partial_a, c_phi_add_N, iphi_add_N, qft, iqft - ) - circuit.append(modulo_multiplier, [up_qreg[i], *down_qreg, *aux_qreg]) - - return circuit.to_instruction() - - @staticmethod - def _validate_input(N: int, a: int): - """Check parameters of the algorithm. - - Args: - N: The odd integer to be factored, has a min. value of 3. - a: Any integer that satisfies 1 < a < N and gcd(a, N) = 1. - - Raises: - ValueError: Invalid input - - """ - validate_min("N", N, 3) - validate_min("a", a, 2) - - if N < 1 or N % 2 == 0: - raise ValueError("The input needs to be an odd integer greater than 1.") - - if a >= N or math.gcd(a, N) != 1: - raise ValueError("The integer a needs to satisfy a < N and gcd(a, N) = 1.") - - def construct_circuit(self, N: int, a: int = 2, measurement: bool = False) -> QuantumCircuit: - """Construct quantum part of the algorithm. - - Args: - N: The odd integer to be factored, has a min. value of 3. - a: Any integer that satisfies 1 < a < N and gcd(a, N) = 1. - measurement: Boolean flag to indicate if measurement should be included in the circuit. - - Returns: - Quantum circuit. - - """ - self._validate_input(N, a) - - # Get n value used in Shor's algorithm, to know how many qubits are used - n = N.bit_length() - - # quantum register where the sequential QFT is performed - up_qreg = QuantumRegister(2 * n, name="up") - # quantum register where the multiplications are made - down_qreg = QuantumRegister(n, name="down") - # auxiliary quantum register used in addition and multiplication - aux_qreg = QuantumRegister(n + 2, name="aux") - - # Create Quantum Circuit - circuit = QuantumCircuit(up_qreg, down_qreg, aux_qreg, name=f"Shor(N={N}, a={a})") - - # Create maximal superposition in top register - circuit.h(up_qreg) - - # Initialize down register to 1 - circuit.x(down_qreg[0]) - - # Apply modulo exponentiation - modulo_power = self._power_mod_N(n, N, a) - circuit.append(modulo_power, circuit.qubits) - - # Apply inverse QFT - iqft = QFT(len(up_qreg)).inverse().to_gate() - circuit.append(iqft, up_qreg) - - if measurement: - up_cqreg = ClassicalRegister(2 * n, name="m") - circuit.add_register(up_cqreg) - circuit.measure(up_qreg, up_cqreg) - - logger.info(summarize_circuits(circuit)) - - return circuit - - @staticmethod - def modinv(a: int, m: int) -> int: - """Returns the modular multiplicative inverse of a with respect to the modulus m.""" - - def egcd(a: int, b: int) -> Tuple[int, int, int]: - if a == 0: - return b, 0, 1 - else: - g, y, x = egcd(b % a, a) - return g, x - (b // a) * y, y - - g, x, _ = egcd(a, m) - if g != 1: - raise ValueError( - "The greatest common divisor of {} and {} is {}, so the " - "modular inverse does not exist.".format(a, m, g) - ) - return x % m - - def _get_factors(self, N: int, a: int, measurement: str) -> Optional[List[int]]: - """Apply the continued fractions to find r and the gcd to find the desired factors.""" - x_final = int(measurement, 2) - logger.info("In decimal, x_final value for this result is: %s.", x_final) - - if x_final <= 0: - fail_reason = "x_final value is <= 0, there are no continued fractions." - else: - fail_reason = None - logger.debug("Running continued fractions for this case.") - - # Calculate T and x/T - T_upper = len(measurement) - T = pow(2, T_upper) - x_over_T = x_final / T - - # Cycle in which each iteration corresponds to putting one more term in the - # calculation of the Continued Fraction (CF) of x/T - - # Initialize the first values according to CF rule - i = 0 - b = array.array("i") - t = array.array("f") - - b.append(math.floor(x_over_T)) - t.append(x_over_T - b[i]) - - exponential = 0.0 - while i < N and fail_reason is None: - # From the 2nd iteration onwards, calculate the new terms of the CF based - # on the previous terms as the rule suggests - if i > 0: - b.append(math.floor(1 / t[i - 1])) - t.append((1 / t[i - 1]) - b[i]) # type: ignore - - # Calculate the denominator of the CF using the known terms - denominator = self._calculate_continued_fraction(b) - - # Increment i for next iteration - i += 1 - - if denominator % 2 == 1: - logger.debug("Odd denominator, will try next iteration of continued fractions.") - continue - - # Denominator is even, try to get factors of N - # Get the exponential a^(r/2) - - if denominator < 1000: - exponential = pow(a, denominator / 2) - - # Check if the value is too big or not - if exponential > 1000000000: - fail_reason = "denominator of continued fraction is too big." - else: - # The value is not too big, - # get the right values and do the proper gcd() - putting_plus = int(exponential + 1) - putting_minus = int(exponential - 1) - one_factor = math.gcd(putting_plus, N) - other_factor = math.gcd(putting_minus, N) - - # Check if the factors found are trivial factors or are the desired factors - if any(factor in {1, N} for factor in (one_factor, other_factor)): - logger.debug("Found just trivial factors, not good enough.") - # Check if the number has already been found, - # (use i - 1 because i was already incremented) - if t[i - 1] == 0: - fail_reason = "the continued fractions found exactly x_final/(2^(2n))." - else: - # Successfully factorized N - return sorted((one_factor, other_factor)) - - # Search for factors failed, write the reason for failure to the debug logs - logger.debug( - "Cannot find factors from measurement %s because %s", - measurement, - fail_reason or "it took too many attempts.", - ) - return None - - @staticmethod - def _calculate_continued_fraction(b: array.array) -> int: - """Calculate the continued fraction of x/T from the current terms of expansion b.""" - - x_over_T = 0 - - for i in reversed(range(len(b) - 1)): - x_over_T = 1 / (b[i + 1] + x_over_T) - - x_over_T += b[0] - - # Get the denominator from the value obtained - frac = fractions.Fraction(x_over_T).limit_denominator() - - logger.debug("Approximation number %s of continued fractions:", len(b)) - logger.debug("Numerator:%s \t\t Denominator: %s.", frac.numerator, frac.denominator) - return frac.denominator - - def factor( - self, - N: int, - a: int = 2, - ) -> "ShorResult": - """Execute the algorithm. - - The input integer :math:`N` to be factored is expected to be odd and greater than 2. - Even though this implementation is general, its capability will be limited by the - capacity of the simulator/hardware. Another input integer :math:`a` can also be supplied, - which needs to be a co-prime smaller than :math:`N` . - - Args: - N: The odd integer to be factored, has a min. value of 3. - a: Any integer that satisfies 1 < a < N and gcd(a, N) = 1. - - Returns: - ShorResult: results of the algorithm. - - Raises: - ValueError: Invalid input - AlgorithmError: If a quantum instance or backend has not been provided - - """ - self._validate_input(N, a) - - if self.quantum_instance is None: - raise AlgorithmError( - "A QuantumInstance or Backend must be supplied to run the quantum algorithm." - ) - - result = ShorResult() - - # check if the input integer N is a power - tf, b, p = is_power(N, return_decomposition=True) - if tf: - logger.info("The input integer is a power: %s=%s^%s.", N, b, p) - result.factors.append(b) - - if not result.factors: - logger.debug("Running with N=%s and a=%s.", N, a) - - if self._quantum_instance.is_statevector: - # Get n value used in Shor's algorithm, to know how many qubits are used - n = N.bit_length() - - circuit = self.construct_circuit(N=N, a=a, measurement=False) - logger.warning( - "The statevector_simulator might lead to " - "subsequent computation using too much memory." - ) - result = self._quantum_instance.execute(circuit) - complete_state_vec = result.get_statevector(circuit) - # TODO: this uses too much memory - up_qreg_density_mat = partial_trace(complete_state_vec, range(2 * n, 4 * n + 2)) - up_qreg_density_mat_diag = np.diag(up_qreg_density_mat) - - counts = {} - for i, v in enumerate(up_qreg_density_mat_diag): - if not v == 0: - counts[bin(int(i))[2:].zfill(2 * n)] = v**2 - else: - circuit = self.construct_circuit(N=N, a=a, measurement=True) - counts = self._quantum_instance.execute(circuit).get_counts(circuit) - - result.total_counts = len(counts) - - # For each simulation result, print proper info to user - # and try to calculate the factors of N - for measurement in list(counts.keys()): - # Get the x_final value from the final state qubits - logger.info("------> Analyzing result %s.", measurement) - factors = self._get_factors(N, a, measurement) - - if factors: - logger.info("Found factors %s from measurement %s.", factors, measurement) - result.successful_counts = result.successful_counts + 1 - if factors not in result.factors: - result.factors.append(factors) - - return result - - -class ShorResult(AlgorithmResult): - """The deprecated Shor Result.""" - - def __init__(self) -> None: - super().__init__() - self._factors = [] - self._total_counts = 0 - self._successful_counts = 0 - - @property - def factors(self) -> List[List[int]]: - """returns factors""" - return self._factors - - @factors.setter - def factors(self, value: List[List[int]]) -> None: - """set factors""" - self._factors = value - - @property - def total_counts(self) -> int: - """returns total counts""" - return self._total_counts - - @total_counts.setter - def total_counts(self, value: int) -> None: - """set total counts""" - self._total_counts = value - - @property - def successful_counts(self) -> int: - """returns successful counts""" - return self._successful_counts - - @successful_counts.setter - def successful_counts(self, value: int) -> None: - """set successful counts""" - self._successful_counts = value diff --git a/qiskit/algorithms/linear_solvers/__init__.py b/qiskit/algorithms/linear_solvers/__init__.py deleted file mode 100644 index 746fb8e40496..000000000000 --- a/qiskit/algorithms/linear_solvers/__init__.py +++ /dev/null @@ -1,92 +0,0 @@ -# This code is part of Qiskit. -# -# (C) Copyright IBM 2020, 2021. -# -# 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 deprecated Linear solvers (:mod:`qiskit.algorithms.linear_solvers`) -======================================================================= -It contains classical and quantum algorithms to solve systems of linear equations such as -:class:`~qiskit.algorithms.HHL`. -Although the quantum algorithm accepts a general Hermitian matrix as input, Qiskit's default -Hamiltonian evolution is exponential in such cases and therefore the quantum linear solver will -not achieve an exponential speedup. -Furthermore, the quantum algorithm can find a solution exponentially faster in the size of the -system than their classical counterparts (i.e. logarithmic complexity instead of polynomial), -meaning that reading the full solution vector would kill such speedup (since this would take -linear time in the size of the system). -Therefore, to achieve an exponential speedup we can only compute functions from the solution -vector (the so called observables) to learn information about the solution. -Known efficient implementations of Hamiltonian evolutions or observables are contained in the -following subfolders: - -`Matrices`_ - A placeholder for efficient implementations of the Hamiltonian evolution of particular types of - matrices. - -`Observables`_ - A placeholder for efficient implementations of functions that can be computed from the solution - vector to a system of linear equations. - -.. currentmodule:: qiskit.algorithms.linear_solvers - -Linear Solvers -============== - -.. autosummary:: - :toctree: ../stubs/ - :nosignatures: - - LinearSolver - LinearSolverResult - HHL - NumPyLinearSolver - -Matrices -======== - -.. autosummary:: - :toctree: ../stubs/ - :nosignatures: - - LinearSystemMatrix - NumPyMatrix - TridiagonalToeplitz - -Observables -=========== - -.. autosummary:: - :toctree: ../stubs/ - - LinearSystemObservable - AbsoluteAverage - MatrixFunctional - -""" - -from .hhl import HHL -from .numpy_linear_solver import NumPyLinearSolver -from .linear_solver import LinearSolver, LinearSolverResult -from .matrices import LinearSystemMatrix, NumPyMatrix, TridiagonalToeplitz -from .observables import LinearSystemObservable, AbsoluteAverage, MatrixFunctional - -__all__ = [ - "HHL", - "NumPyLinearSolver", - "LinearSolver", - "LinearSolverResult", - "LinearSystemMatrix", - "NumPyMatrix", - "TridiagonalToeplitz", - "LinearSystemObservable", - "AbsoluteAverage", - "MatrixFunctional", -] diff --git a/qiskit/algorithms/linear_solvers/hhl.py b/qiskit/algorithms/linear_solvers/hhl.py deleted file mode 100644 index 97cdf09cb86d..000000000000 --- a/qiskit/algorithms/linear_solvers/hhl.py +++ /dev/null @@ -1,548 +0,0 @@ -# This code is part of Qiskit. -# -# (C) Copyright IBM 2020, 2022. -# -# 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 HHL algorithm.""" - -from typing import Optional, Union, List, Callable, Tuple -import numpy as np - -from qiskit.circuit import QuantumCircuit, QuantumRegister, AncillaRegister -from qiskit.circuit.library import PhaseEstimation -from qiskit.circuit.library.arithmetic.piecewise_chebyshev import PiecewiseChebyshev -from qiskit.circuit.library.arithmetic.exact_reciprocal import ExactReciprocal -from qiskit.opflow import ( - Z, - I, - StateFn, - TensoredOp, - ExpectationBase, - CircuitSampler, - ListOp, - ExpectationFactory, -) -from qiskit.providers import Backend -from qiskit.quantum_info.operators.base_operator import BaseOperator -from qiskit.utils import QuantumInstance -from qiskit.utils.deprecation import deprecate_function - -from .linear_solver import LinearSolver, LinearSolverResult -from .matrices.numpy_matrix import NumPyMatrix -from .observables.linear_system_observable import LinearSystemObservable - - -class HHL(LinearSolver): - r"""The deprecated systems of linear equations arise naturally in many real-life applications - in a wide range - of areas, such as in the solution of Partial Differential Equations, the calibration of - financial models, fluid simulation or numerical field calculation. The problem can be defined - as, given a matrix :math:`A\in\mathbb{C}^{N\times N}` and a vector - :math:`\vec{b}\in\mathbb{C}^{N}`, find :math:`\vec{x}\in\mathbb{C}^{N}` satisfying - :math:`A\vec{x}=\vec{b}`. - - A system of linear equations is called :math:`s`-sparse if :math:`A` has at most :math:`s` - non-zero entries per row or column. Solving an :math:`s`-sparse system of size :math:`N` with - a classical computer requires :math:`\mathcal{ O }(Ns\kappa\log(1/\epsilon))` running time - using the conjugate gradient method. Here :math:`\kappa` denotes the condition number of the - system and :math:`\epsilon` the accuracy of the approximation. - - The deprecated HHL is a quantum algorithm to estimate a function of the solution with running time - complexity of :math:`\mathcal{ O }(\log(N)s^{2}\kappa^{2}/\epsilon)` when - :math:`A` is a Hermitian matrix under the assumptions of efficient oracles for loading the - data, Hamiltonian simulation and computing a function of the solution. This is an exponential - speed up in the size of the system, however one crucial remark to keep in mind is that the - classical algorithm returns the full solution, while the HHL can only approximate functions of - the solution vector. - - The HHL class is deprecated as of Qiskit Terra 0.22.0 - and will be removed no sooner than 3 months after the release date. - It is replaced by the tutorial at - `HHL `_ - - Examples:: - - import warnings - import numpy as np - from qiskit import QuantumCircuit - from qiskit.algorithms.linear_solvers.hhl import HHL - from qiskit.algorithms.linear_solvers.matrices import TridiagonalToeplitz - from qiskit.algorithms.linear_solvers.observables import MatrixFunctional - - with warnings.catch_warnings(): - warnings.simplefilter('ignore') - matrix = TridiagonalToeplitz(2, 1, 1 / 3, trotter_steps=2) - right_hand_side = [1.0, -2.1, 3.2, -4.3] - observable = MatrixFunctional(1, 1 / 2) - rhs = right_hand_side / np.linalg.norm(right_hand_side) - - # Initial state circuit - num_qubits = matrix.num_state_qubits - qc = QuantumCircuit(num_qubits) - qc.isometry(rhs, list(range(num_qubits)), None) - - with warnings.catch_warnings(): - warnings.simplefilter('ignore') - hhl = HHL() - solution = hhl.solve(matrix, qc, observable) - approx_result = solution.observable - - References: - - [1]: Harrow, A. W., Hassidim, A., Lloyd, S. (2009). - Quantum algorithm for linear systems of equations. - `Phys. Rev. Lett. 103, 15 (2009), 1–15. `_ - - [2]: Carrera Vazquez, A., Hiptmair, R., & Woerner, S. (2022). - Enhancing the Quantum Linear Systems Algorithm Using Richardson Extrapolation. - `ACM Transactions on Quantum Computing 3, 1, Article 2 `_ - - """ - - @deprecate_function( - "The HHL class is deprecated as of Qiskit Terra 0.22.0 and will be removed " - "no sooner than 3 months after the release date. It is replaced by the tutorial at " - "https://qiskit.org/textbook/ch-applications/hhl_tutorial.html", - since="0.22.0", - ) - def __init__( - self, - epsilon: float = 1e-2, - expectation: Optional[ExpectationBase] = None, - quantum_instance: Optional[Union[Backend, QuantumInstance]] = None, - ) -> None: - r""" - Args: - epsilon: Error tolerance of the approximation to the solution, i.e. if :math:`x` is the - exact solution and :math:`\tilde{x}` the one calculated by the algorithm, then - :math:`||x - \tilde{x}|| \le epsilon`. - expectation: The expectation converter applied to the expectation values before - evaluation. If None then PauliExpectation is used. - quantum_instance: Quantum Instance or Backend. If None, a Statevector calculation is - done. - """ - super().__init__() - - self._epsilon = epsilon - # Tolerance for the different parts of the algorithm as per [1] - self._epsilon_r = epsilon / 3 # conditioned rotation - self._epsilon_s = epsilon / 3 # state preparation - self._epsilon_a = epsilon / 6 # hamiltonian simulation - - self._scaling = None # scaling of the solution - - self._sampler = None - self.quantum_instance = quantum_instance - - self._expectation = expectation - - # For now the default reciprocal implementation is exact - self._exact_reciprocal = True - # Set the default scaling to 1 - self.scaling = 1 - - @property - def quantum_instance(self) -> Optional[QuantumInstance]: - """Get the quantum instance. - - Returns: - The quantum instance used to run this algorithm. - """ - return None if self._sampler is None else self._sampler.quantum_instance - - @quantum_instance.setter - def quantum_instance(self, quantum_instance: Optional[Union[QuantumInstance, Backend]]) -> None: - """Set quantum instance. - - Args: - quantum_instance: The quantum instance used to run this algorithm. - If None, a Statevector calculation is done. - """ - if quantum_instance is not None: - self._sampler = CircuitSampler(quantum_instance) - else: - self._sampler = None - - @property - def scaling(self) -> float: - """The scaling of the solution vector.""" - return self._scaling - - @scaling.setter - def scaling(self, scaling: float) -> None: - """Set the new scaling of the solution vector.""" - self._scaling = scaling - - @property - def expectation(self) -> ExpectationBase: - """The expectation value algorithm used to construct the expectation measurement from - the observable.""" - return self._expectation - - @expectation.setter - def expectation(self, expectation: ExpectationBase) -> None: - """Set the expectation value algorithm.""" - self._expectation = expectation - - def _get_delta(self, n_l: int, lambda_min: float, lambda_max: float) -> float: - """Calculates the scaling factor to represent exactly lambda_min on nl binary digits. - - Args: - n_l: The number of qubits to represent the eigenvalues. - lambda_min: the smallest eigenvalue. - lambda_max: the largest eigenvalue. - - Returns: - The value of the scaling factor. - """ - formatstr = "#0" + str(n_l + 2) + "b" - lambda_min_tilde = np.abs(lambda_min * (2**n_l - 1) / lambda_max) - # floating point precision can cause problems - if np.abs(lambda_min_tilde - 1) < 1e-7: - lambda_min_tilde = 1 - binstr = format(int(lambda_min_tilde), formatstr)[2::] - lamb_min_rep = 0 - for i, char in enumerate(binstr): - lamb_min_rep += int(char) / (2 ** (i + 1)) - return lamb_min_rep - - def _calculate_norm(self, qc: QuantumCircuit) -> float: - """Calculates the value of the euclidean norm of the solution. - - Args: - qc: The quantum circuit preparing the solution x to the system. - - Returns: - The value of the euclidean norm of the solution. - """ - # Calculate the number of qubits - nb = qc.qregs[0].size - nl = qc.qregs[1].size - na = qc.num_ancillas - - # Create the Operators Zero and One - zero_op = (I + Z) / 2 - one_op = (I - Z) / 2 - - # Norm observable - observable = one_op ^ TensoredOp((nl + na) * [zero_op]) ^ (I ^ nb) - norm_2 = (~StateFn(observable) @ StateFn(qc)).eval() - - return np.real(np.sqrt(norm_2) / self.scaling) - - def _calculate_observable( - self, - solution: QuantumCircuit, - observable: Optional[Union[LinearSystemObservable, BaseOperator]] = None, - observable_circuit: Optional[QuantumCircuit] = None, - post_processing: Optional[ - Callable[[Union[float, List[float]]], Union[float, List[float]]] - ] = None, - ) -> Tuple[Union[float, List[float]], Union[float, List[float]]]: - """Calculates the value of the observable(s) given. - - Args: - solution: The quantum circuit preparing the solution x to the system. - observable: Information to be extracted from the solution. - observable_circuit: Circuit to be applied to the solution to extract information. - post_processing: Function to compute the value of the observable. - - Returns: - The value of the observable(s) and the circuit results before post-processing as a - tuple. - """ - # Get the number of qubits - nb = solution.qregs[0].size - nl = solution.qregs[1].size - na = solution.num_ancillas - - # if the observable is given construct post_processing and observable_circuit - if observable is not None: - observable_circuit = observable.observable_circuit(nb) - post_processing = observable.post_processing - - if isinstance(observable, LinearSystemObservable): - observable = observable.observable(nb) - - # in the other case use the identity as observable - else: - observable = I ^ nb - - # Create the Operators Zero and One - zero_op = (I + Z) / 2 - one_op = (I - Z) / 2 - - is_list = True - if not isinstance(observable_circuit, list): - is_list = False - observable_circuit = [observable_circuit] - observable = [observable] - - expectations = [] - for circ, obs in zip(observable_circuit, observable): - circuit = QuantumCircuit(solution.num_qubits) - circuit.append(solution, circuit.qubits) - circuit.append(circ, range(nb)) - - ob = one_op ^ TensoredOp((nl + na) * [zero_op]) ^ obs - expectations.append(~StateFn(ob) @ StateFn(circuit)) - - if is_list: - # execute all in a list op to send circuits in batches - expectations = ListOp(expectations) - else: - expectations = expectations[0] - - # check if an expectation converter is given - if self._expectation is not None: - expectations = self._expectation.convert(expectations) - # if otherwise a backend was specified, try to set the best expectation value - elif self._sampler is not None: - if is_list: - op = expectations.oplist[0] - else: - op = expectations - self._expectation = ExpectationFactory.build(op, self._sampler.quantum_instance) - - if self._sampler is not None: - expectations = self._sampler.convert(expectations) - - # evaluate - expectation_results = expectations.eval() - - # apply post_processing - result = post_processing(expectation_results, nb, self.scaling) - - return result, expectation_results - - def construct_circuit( - self, - matrix: Union[List, np.ndarray, QuantumCircuit], - vector: Union[List, np.ndarray, QuantumCircuit], - neg_vals: Optional[bool] = True, - ) -> QuantumCircuit: - """Construct the HHL circuit. - - Args: - matrix: The matrix specifying the system, i.e. A in Ax=b. - vector: The vector specifying the right hand side of the equation in Ax=b. - neg_vals: States whether the matrix has negative eigenvalues. If False the - computation becomes cheaper. - - Returns: - The HHL circuit. - - Raises: - ValueError: If the input is not in the correct format. - ValueError: If the type of the input matrix is not supported. - """ - # State preparation circuit - default is qiskit - if isinstance(vector, QuantumCircuit): - nb = vector.num_qubits - vector_circuit = vector - elif isinstance(vector, (list, np.ndarray)): - if isinstance(vector, list): - vector = np.array(vector) - nb = int(np.log2(len(vector))) - vector_circuit = QuantumCircuit(nb) - vector_circuit.isometry(vector / np.linalg.norm(vector), list(range(nb)), None) - - # If state preparation is probabilistic the number of qubit flags should increase - nf = 1 - - # Hamiltonian simulation circuit - default is Trotterization - if isinstance(matrix, QuantumCircuit): - matrix_circuit = matrix - elif isinstance(matrix, (list, np.ndarray)): - if isinstance(matrix, list): - matrix = np.array(matrix) - - if matrix.shape[0] != matrix.shape[1]: - raise ValueError("Input matrix must be square!") - if np.log2(matrix.shape[0]) % 1 != 0: - raise ValueError("Input matrix dimension must be 2^n!") - if not np.allclose(matrix, matrix.conj().T): - raise ValueError("Input matrix must be hermitian!") - if matrix.shape[0] != 2**vector_circuit.num_qubits: - raise ValueError( - "Input vector dimension does not match input " - "matrix dimension! Vector dimension: " - + str(vector_circuit.num_qubits) - + ". Matrix dimension: " - + str(matrix.shape[0]) - ) - matrix_circuit = NumPyMatrix(matrix, evolution_time=2 * np.pi) - else: - raise ValueError(f"Invalid type for matrix: {type(matrix)}.") - - # Set the tolerance for the matrix approximation - if hasattr(matrix_circuit, "tolerance"): - matrix_circuit.tolerance = self._epsilon_a - - # check if the matrix can calculate the condition number and store the upper bound - if ( - hasattr(matrix_circuit, "condition_bounds") - and matrix_circuit.condition_bounds() is not None - ): - kappa = matrix_circuit.condition_bounds()[1] - else: - kappa = 1 - # Update the number of qubits required to represent the eigenvalues - # The +neg_vals is to register negative eigenvalues because - # e^{-2 \pi i \lambda} = e^{2 \pi i (1 - \lambda)} - nl = max(nb + 1, int(np.ceil(np.log2(kappa + 1)))) + neg_vals - - # check if the matrix can calculate bounds for the eigenvalues - if hasattr(matrix_circuit, "eigs_bounds") and matrix_circuit.eigs_bounds() is not None: - lambda_min, lambda_max = matrix_circuit.eigs_bounds() - # Constant so that the minimum eigenvalue is represented exactly, since it contributes - # the most to the solution of the system. -1 to take into account the sign qubit - delta = self._get_delta(nl - neg_vals, lambda_min, lambda_max) - # Update evolution time - matrix_circuit.evolution_time = 2 * np.pi * delta / lambda_min / (2**neg_vals) - # Update the scaling of the solution - self.scaling = lambda_min - else: - delta = 1 / (2**nl) - print("The solution will be calculated up to a scaling factor.") - - if self._exact_reciprocal: - reciprocal_circuit = ExactReciprocal(nl, delta, neg_vals=neg_vals) - # Update number of ancilla qubits - na = matrix_circuit.num_ancillas - else: - # Calculate breakpoints for the reciprocal approximation - num_values = 2**nl - constant = delta - a = int(round(num_values ** (2 / 3))) - - # Calculate the degree of the polynomial and the number of intervals - r = 2 * constant / a + np.sqrt(np.abs(1 - (2 * constant / a) ** 2)) - degree = min( - nb, - int( - np.log( - 1 - + ( - 16.23 - * np.sqrt(np.log(r) ** 2 + (np.pi / 2) ** 2) - * kappa - * (2 * kappa - self._epsilon_r) - ) - / self._epsilon_r - ) - ), - ) - num_intervals = int(np.ceil(np.log((num_values - 1) / a) / np.log(5))) - - # Calculate breakpoints and polynomials - breakpoints = [] - for i in range(0, num_intervals): - # Add the breakpoint to the list - breakpoints.append(a * (5**i)) - - # Define the right breakpoint of the interval - if i == num_intervals - 1: - breakpoints.append(num_values - 1) - - reciprocal_circuit = PiecewiseChebyshev( - lambda x: np.arcsin(constant / x), degree, breakpoints, nl - ) - na = max(matrix_circuit.num_ancillas, reciprocal_circuit.num_ancillas) - - # Initialise the quantum registers - qb = QuantumRegister(nb) # right hand side and solution - ql = QuantumRegister(nl) # eigenvalue evaluation qubits - if na > 0: - qa = AncillaRegister(na) # ancilla qubits - qf = QuantumRegister(nf) # flag qubits - - if na > 0: - qc = QuantumCircuit(qb, ql, qa, qf) - else: - qc = QuantumCircuit(qb, ql, qf) - - # State preparation - qc.append(vector_circuit, qb[:]) - # QPE - phase_estimation = PhaseEstimation(nl, matrix_circuit) - if na > 0: - qc.append(phase_estimation, ql[:] + qb[:] + qa[: matrix_circuit.num_ancillas]) - else: - qc.append(phase_estimation, ql[:] + qb[:]) - # Conditioned rotation - if self._exact_reciprocal: - qc.append(reciprocal_circuit, ql[::-1] + [qf[0]]) - else: - qc.append( - reciprocal_circuit.to_instruction(), - ql[:] + [qf[0]] + qa[: reciprocal_circuit.num_ancillas], - ) - # QPE inverse - if na > 0: - qc.append(phase_estimation.inverse(), ql[:] + qb[:] + qa[: matrix_circuit.num_ancillas]) - else: - qc.append(phase_estimation.inverse(), ql[:] + qb[:]) - return qc - - def solve( - self, - matrix: Union[List, np.ndarray, QuantumCircuit], - vector: Union[List, np.ndarray, QuantumCircuit], - observable: Optional[ - Union[ - LinearSystemObservable, - BaseOperator, - List[LinearSystemObservable], - List[BaseOperator], - ] - ] = None, - observable_circuit: Optional[Union[QuantumCircuit, List[QuantumCircuit]]] = None, - post_processing: Optional[ - Callable[[Union[float, List[float]]], Union[float, List[float]]] - ] = None, - ) -> LinearSolverResult: - """Tries to solve the given linear system of equations. - - Args: - matrix: The matrix specifying the system, i.e. A in Ax=b. - vector: The vector specifying the right hand side of the equation in Ax=b. - observable: Optional information to be extracted from the solution. - Default is the probability of success of the algorithm. - observable_circuit: Optional circuit to be applied to the solution to extract - information. Default is `None`. - post_processing: Optional function to compute the value of the observable. - Default is the raw value of measuring the observable. - - Raises: - ValueError: If an invalid combination of observable, observable_circuit and - post_processing is passed. - - Returns: - The result object containing information about the solution vector of the linear - system. - """ - # verify input - if observable is not None: - if observable_circuit is not None or post_processing is not None: - raise ValueError( - "If observable is passed, observable_circuit and post_processing cannot be set." - ) - - solution = LinearSolverResult() - solution.state = self.construct_circuit(matrix, vector) - solution.euclidean_norm = self._calculate_norm(solution.state) - - if observable is not None or observable_circuit is not None: - solution.observable, solution.circuit_results = self._calculate_observable( - solution.state, observable, observable_circuit, post_processing - ) - - return solution diff --git a/qiskit/algorithms/linear_solvers/linear_solver.py b/qiskit/algorithms/linear_solvers/linear_solver.py deleted file mode 100644 index bdf1bfa451ab..000000000000 --- a/qiskit/algorithms/linear_solvers/linear_solver.py +++ /dev/null @@ -1,145 +0,0 @@ -# This code is part of Qiskit. -# -# (C) Copyright IBM 2020, 2022. -# -# 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. - -"""An abstract class for linear systems solvers.""" - -from abc import ABC, abstractmethod -from typing import Union, Optional, List, Callable -import numpy as np - -from qiskit import QuantumCircuit -from qiskit.result import Result -from qiskit.quantum_info.operators.base_operator import BaseOperator -from qiskit.utils.deprecation import deprecate_function - -from .observables.linear_system_observable import LinearSystemObservable -from ..algorithm_result import AlgorithmResult - - -class LinearSolverResult(AlgorithmResult): - """The deprecated base class for linear systems results. - - The linear systems algorithms return an object of the type ``LinearSystemsResult`` - with the information about the solution obtained. - """ - - @deprecate_function( - "The LinearSolverResult class is deprecated as of Qiskit Terra 0.22.0 " - "and will be removed no sooner than 3 months after the release date.", - since="0.22.0", - ) - def __init__(self) -> None: - super().__init__() - - # Set the default to None, if the algorithm knows how to calculate it can override it. - self._state = None - self._observable = None - self._euclidean_norm = None - self._circuit_results = None - - @property - def observable(self) -> Union[float, List[float]]: - """return the (list of) calculated observable(s)""" - return self._observable - - @observable.setter - def observable(self, observable: Union[float, List[float]]) -> None: - """Set the value(s) of the observable(s). - - Args: - observable: The new value(s) of the observable(s). - """ - self._observable = observable - - @property - def state(self) -> Union[QuantumCircuit, np.ndarray]: - """return either the circuit that prepares the solution or the solution as a vector""" - return self._state - - @state.setter - def state(self, state: Union[QuantumCircuit, np.ndarray]) -> None: - """Set the solution state as either the circuit that prepares it or as a vector. - - Args: - state: The new solution state. - """ - self._state = state - - @property - def euclidean_norm(self) -> float: - """return the euclidean norm if the algorithm knows how to calculate it""" - return self._euclidean_norm - - @euclidean_norm.setter - def euclidean_norm(self, norm: float) -> None: - """Set the euclidean norm of the solution. - - Args: - norm: The new euclidean norm of the solution. - """ - self._euclidean_norm = norm - - @property - def circuit_results(self) -> Union[List[float], List[Result]]: - """return the results from the circuits""" - return self._circuit_results - - @circuit_results.setter - def circuit_results(self, results: Union[List[float], List[Result]]): - self._circuit_results = results - - -class LinearSolver(ABC): - """The deprecated abstract class for linear system solvers in Qiskit.""" - - @deprecate_function( - "The LinearSolver class is deprecated as of Qiskit Terra 0.22.0 " - "and will be removed no sooner than 3 months after the release date.", - since="0.22.0", - ) - def __init__(self) -> None: - pass - - @abstractmethod - def solve( - self, - matrix: Union[np.ndarray, QuantumCircuit], - vector: Union[np.ndarray, QuantumCircuit], - observable: Optional[ - Union[ - LinearSystemObservable, - BaseOperator, - List[LinearSystemObservable], - List[BaseOperator], - ] - ] = None, - observable_circuit: Optional[Union[QuantumCircuit, List[QuantumCircuit]]] = None, - post_processing: Optional[ - Callable[[Union[float, List[float]]], Union[float, List[float]]] - ] = None, - ) -> LinearSolverResult: - """Solve the system and compute the observable(s) - - Args: - matrix: The matrix specifying the system, i.e. A in Ax=b. - vector: The vector specifying the right hand side of the equation in Ax=b. - observable: Optional information to be extracted from the solution. - Default is the probability of success of the algorithm. - observable_circuit: Optional circuit to be applied to the solution to extract - information. Default is ``None``. - post_processing: Optional function to compute the value of the observable. - Default is the raw value of measuring the observable. - - Returns: - The result of the linear system. - """ - raise NotImplementedError diff --git a/qiskit/algorithms/linear_solvers/matrices/__init__.py b/qiskit/algorithms/linear_solvers/matrices/__init__.py deleted file mode 100644 index 97b34d9c624e..000000000000 --- a/qiskit/algorithms/linear_solvers/matrices/__init__.py +++ /dev/null @@ -1,19 +0,0 @@ -# This code is part of Qiskit. -# -# (C) Copyright IBM 2020, 2021. -# -# 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. - -"""System matrices for Qiskit's linear solvers.""" - -from .linear_system_matrix import LinearSystemMatrix -from .numpy_matrix import NumPyMatrix -from .tridiagonal_toeplitz import TridiagonalToeplitz - -__all__ = ["LinearSystemMatrix", "NumPyMatrix", "TridiagonalToeplitz"] diff --git a/qiskit/algorithms/linear_solvers/matrices/linear_system_matrix.py b/qiskit/algorithms/linear_solvers/matrices/linear_system_matrix.py deleted file mode 100644 index 5a66ed9e1d53..000000000000 --- a/qiskit/algorithms/linear_solvers/matrices/linear_system_matrix.py +++ /dev/null @@ -1,140 +0,0 @@ -# This code is part of Qiskit. -# -# (C) Copyright IBM 2020, 2022. -# -# 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. - -"""An abstract class for matrices input to the linear systems solvers in Qiskit.""" - -from abc import ABC, abstractmethod -from typing import Tuple - -from qiskit import QuantumCircuit -from qiskit.circuit.library import BlueprintCircuit -from qiskit.utils.deprecation import deprecate_function - - -class LinearSystemMatrix(BlueprintCircuit, ABC): - """The deprecated base class for linear system matrices.""" - - @deprecate_function( - "The LinearSystemMatrix class is deprecated as of Qiskit Terra 0.22.0 " - "and will be removed no sooner than 3 months after the release date.", - since="0.22.0", - ) - def __init__( - self, - num_state_qubits: int, - tolerance: float, - evolution_time: float, - name: str = "ls_matrix", - ) -> None: - """ - Args: - num_state_qubits: the number of qubits where the unitary acts. - tolerance: the accuracy desired for the approximation - evolution_time: the time of the Hamiltonian simulation - name: The name of the object. - """ - super().__init__(name=name) - - # define internal parameters - self._num_state_qubits = None - self._tolerance = None - self._evolution_time = None # makes sure the eigenvalues are contained in [0,1) - - # store parameters - self.num_state_qubits = num_state_qubits - self.tolerance = tolerance - self.evolution_time = evolution_time - - @property - def num_state_qubits(self) -> int: - r"""The number of state qubits representing the state :math:`|x\rangle`. - - Returns: - The number of state qubits. - """ - return self._num_state_qubits - - @num_state_qubits.setter - def num_state_qubits(self, num_state_qubits: int) -> None: - """Set the number of state qubits. - - Note that this may change the underlying quantum register, if the number of state qubits - changes. - - Args: - num_state_qubits: The new number of qubits. - """ - if num_state_qubits != self._num_state_qubits: - self._invalidate() - self._num_state_qubits = num_state_qubits - self._reset_registers(num_state_qubits) - - @property - def tolerance(self) -> float: - """Return the error tolerance""" - return self._tolerance - - @tolerance.setter - def tolerance(self, tolerance: float) -> None: - """Set the error tolerance - Args: - tolerance: The new error tolerance. - """ - self._tolerance = tolerance - - @property - def evolution_time(self) -> float: - """Return the time of the evolution.""" - return self._evolution_time - - @evolution_time.setter - def evolution_time(self, evolution_time: float) -> None: - """Set the time of the evolution. - - Args: - evolution_time: The new time of the evolution. - """ - self._evolution_time = evolution_time - - @abstractmethod - def eigs_bounds(self) -> Tuple[float, float]: - """Return lower and upper bounds on the eigenvalues of the matrix.""" - raise NotImplementedError - - @abstractmethod - def condition_bounds(self) -> Tuple[float, float]: - """Return lower and upper bounds on the condition number of the matrix.""" - raise NotImplementedError - - @abstractmethod - def _reset_registers(self, num_state_qubits: int) -> None: - """Reset the registers according to the new number of state qubits. - - Args: - num_state_qubits: The new number of qubits. - """ - raise NotImplementedError - - @abstractmethod - def power(self, power: int, matrix_power: bool = False) -> QuantumCircuit: - """Build powers of the circuit. - - Args: - power: The power to raise this circuit to. - matrix_power: If True, the circuit is converted to a matrix and then the - matrix power is computed. If False, and ``power`` is a positive integer, - the implementation defaults to ``repeat``. - - Returns: - The quantum circuit implementing powers of the unitary. - """ - raise NotImplementedError diff --git a/qiskit/algorithms/linear_solvers/matrices/numpy_matrix.py b/qiskit/algorithms/linear_solvers/matrices/numpy_matrix.py deleted file mode 100644 index bd441096a578..000000000000 --- a/qiskit/algorithms/linear_solvers/matrices/numpy_matrix.py +++ /dev/null @@ -1,221 +0,0 @@ -# This code is part of Qiskit. -# -# (C) Copyright IBM 2020, 2022. -# -# 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. - -"""Hamiltonian simulation of matrices given as numpy arrays.""" - -from typing import Tuple -import numpy as np -import scipy as sp - -from qiskit import QuantumCircuit, QuantumRegister -from qiskit.utils.deprecation import deprecate_function - -from .linear_system_matrix import LinearSystemMatrix - - -class NumPyMatrix(LinearSystemMatrix): - """The deprecated class of matrices given as a numpy array. - - Examples:: - - import warnings - import numpy as np - from qiskit import QuantumCircuit - from qiskit.algorithms.linear_solvers.matrices.numpy_matrix import NumPyMatrix - - with warnings.catch_warnings(): - warnings.simplefilter('ignore') - matrix = NumPyMatrix(np.array([[1 / 2, 1 / 6, 0, 0], [1 / 6, 1 / 2, 1 / 6, 0], - [0, 1 / 6, 1 / 2, 1 / 6], [0, 0, 1 / 6, 1 / 2]])) - power = 2 - - num_qubits = matrix.num_state_qubits - # Controlled power (as used within QPE) - pow_circ = matrix.power(power).control() - circ_qubits = pow_circ.num_qubits - qc = QuantumCircuit(circ_qubits) - qc.append(matrix.power(power).control(), list(range(circ_qubits))) - """ - - @deprecate_function( - "The NumPyMatrix class is deprecated as of Qiskit Terra 0.22.0 " - "and will be removed no sooner than 3 months after the release date.", - since="0.22.0", - ) - def __init__( - self, - matrix: np.ndarray, - tolerance: float = 1e-2, - evolution_time: float = 1.0, - name: str = "np_matrix", - ) -> None: - """ - Args: - matrix: The matrix defining the linear system problem. - tolerance: The accuracy desired for the approximation. - evolution_time: The time of the Hamiltonian simulation. - name: The name of the object. - """ - - # define internal parameters - self._num_state_qubits = None - self._tolerance = None - self._evolution_time = None # makes sure the eigenvalues are contained in [0,1) - self._matrix = None - - super().__init__( - num_state_qubits=int(np.log2(matrix.shape[0])), - tolerance=tolerance, - evolution_time=evolution_time, - name=name, - ) - - # store parameters - self.num_state_qubits = int(np.log2(matrix.shape[0])) - self.tolerance = tolerance - self.evolution_time = evolution_time - self.matrix = matrix - - @property - def num_state_qubits(self) -> int: - r"""The number of state qubits representing the state :math:`|x\rangle`. - - Returns: - The number of state qubits. - """ - return self._num_state_qubits - - @num_state_qubits.setter - def num_state_qubits(self, num_state_qubits: int) -> None: - """Set the number of state qubits. - - Note that this may change the underlying quantum register, if the number of state qubits - changes. - - Args: - num_state_qubits: The new number of qubits. - """ - if num_state_qubits != self._num_state_qubits: - self._invalidate() - self._num_state_qubits = num_state_qubits - self._reset_registers(num_state_qubits) - - @property - def tolerance(self) -> float: - """Return the error tolerance""" - return self._tolerance - - @tolerance.setter - def tolerance(self, tolerance: float) -> None: - """Set the error tolerance - Args: - tolerance: The new error tolerance. - """ - self._tolerance = tolerance - - @property - def evolution_time(self) -> float: - """Return the time of the evolution.""" - return self._evolution_time - - @evolution_time.setter - def evolution_time(self, evolution_time: float) -> None: - """Set the time of the evolution. - - Args: - evolution_time: The new time of the evolution. - """ - self._evolution_time = evolution_time - - @property - def matrix(self) -> np.ndarray: - """Return the matrix.""" - return self._matrix - - @matrix.setter - def matrix(self, matrix: np.ndarray) -> None: - """Set the matrix. - - Args: - matrix: The new matrix. - """ - self._matrix = matrix - - def eigs_bounds(self) -> Tuple[float, float]: - """Return lower and upper bounds on the eigenvalues of the matrix.""" - matrix_array = self.matrix - lambda_max = max(np.abs(np.linalg.eigvals(matrix_array))) - lambda_min = min(np.abs(np.linalg.eigvals(matrix_array))) - return lambda_min, lambda_max - - def condition_bounds(self) -> Tuple[float, float]: - """Return lower and upper bounds on the condition number of the matrix.""" - matrix_array = self.matrix - kappa = np.linalg.cond(matrix_array) - return kappa, kappa - - def _check_configuration(self, raise_on_failure: bool = True) -> bool: - """Check if the current configuration is valid.""" - valid = True - - if self.matrix.shape[0] != self.matrix.shape[1]: - if raise_on_failure: - raise AttributeError("Input matrix must be square!") - return False - if np.log2(self.matrix.shape[0]) % 1 != 0: - if raise_on_failure: - raise AttributeError("Input matrix dimension must be 2^n!") - return False - if not np.allclose(self.matrix, self.matrix.conj().T): - if raise_on_failure: - raise AttributeError("Input matrix must be hermitian!") - return False - - return valid - - def _reset_registers(self, num_state_qubits: int) -> None: - """Reset the quantum registers. - - Args: - num_state_qubits: The number of qubits to represent the matrix. - """ - qr_state = QuantumRegister(num_state_qubits, "state") - self.qregs = [qr_state] - - def _build(self) -> None: - """If not already built, build the circuit.""" - if self._is_built: - return - - super()._build() - - self.compose(self.power(1), inplace=True) - - def inverse(self): - return NumPyMatrix(self.matrix, evolution_time=-1 * self.evolution_time) - - def power(self, power: int, matrix_power: bool = False) -> QuantumCircuit: - """Build powers of the circuit. - - Args: - power: The power to raise this circuit to. - matrix_power: If True, the circuit is converted to a matrix and then the - matrix power is computed. If False, and ``power`` is a positive integer, - the implementation defaults to ``repeat``. - - Returns: - The quantum circuit implementing powers of the unitary. - """ - qc = QuantumCircuit(self.num_state_qubits) - evolved = sp.linalg.expm(1j * self.matrix * self.evolution_time) - qc.unitary(evolved, qc.qubits) - return qc.power(power) diff --git a/qiskit/algorithms/linear_solvers/matrices/tridiagonal_toeplitz.py b/qiskit/algorithms/linear_solvers/matrices/tridiagonal_toeplitz.py deleted file mode 100644 index d29e93dd3a03..000000000000 --- a/qiskit/algorithms/linear_solvers/matrices/tridiagonal_toeplitz.py +++ /dev/null @@ -1,486 +0,0 @@ -# This code is part of Qiskit. -# -# (C) Copyright IBM 2020, 2022. -# -# 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. - -"""Hamiltonian simulation of tridiagonal Toeplitz symmetric matrices.""" - -from typing import Tuple -import numpy as np -from scipy.sparse import diags - -from qiskit.circuit import QuantumCircuit, QuantumRegister, AncillaRegister -from qiskit.circuit.library import UGate, MCMTVChain -from qiskit.utils.deprecation import deprecate_function - -from .linear_system_matrix import LinearSystemMatrix - - -class TridiagonalToeplitz(LinearSystemMatrix): - r"""The deprecated class of tridiagonal Toeplitz symmetric matrices. - - Given the main entry, :math:`a`, and the off diagonal entry, :math:`b`, the :math:`4\times 4` - dimensional tridiagonal Toeplitz symmetric matrix is - - .. math:: - - \begin{pmatrix} - a & b & 0 & 0 \\ - b & a & b & 0 \\ - 0 & b & a & b \\ - 0 & 0 & b & a - \end{pmatrix}. - - Examples:: - - import warnings - import numpy as np - from qiskit import QuantumCircuit - from qiskit.algorithms.linear_solvers.matrices import TridiagonalToeplitz - - with warnings.catch_warnings(): - warnings.simplefilter('ignore') - matrix = TridiagonalToeplitz(2, 1, -1 / 3) - power = 3 - - # Controlled power (as within QPE) - num_qubits = matrix.num_state_qubits - pow_circ = matrix.power(power).control() - circ_qubits = pow_circ.num_qubits - qc = QuantumCircuit(circ_qubits) - qc.append(matrix.power(power).control(), list(range(circ_qubits))) - """ - - @deprecate_function( - "The TridiagonalToeplitz class is deprecated as of Qiskit Terra 0.22.0 " - "and will be removed no sooner than 3 months after the release date.", - since="0.22.0", - ) - def __init__( - self, - num_state_qubits: int, - main_diag: float, - off_diag: float, - tolerance: float = 1e-2, - evolution_time: float = 1.0, - trotter_steps: int = 1, - name: str = "tridi", - ) -> None: - """ - Args: - num_state_qubits: the number of qubits where the unitary acts. - main_diag: the main diagonal entry - off_diag: the off diagonal entry - tolerance: the accuracy desired for the approximation - evolution_time: the time of the Hamiltonian simulation - trotter_steps: the number of Trotter steps - name: The name of the object. - """ - # define internal parameters - self._num_state_qubits = None - self._main_diag = None - self._off_diag = None - self._tolerance = None - self._evolution_time = None # makes sure the eigenvalues are contained in [0,1) - self._trotter_steps = None - - # store parameters - self.main_diag = main_diag - self.off_diag = off_diag - super().__init__( - num_state_qubits=num_state_qubits, - tolerance=tolerance, - evolution_time=evolution_time, - name=name, - ) - self.trotter_steps = trotter_steps - - @property - def num_state_qubits(self) -> int: - r"""The number of state qubits representing the state :math:`|x\rangle`. - - Returns: - The number of state qubits. - """ - return self._num_state_qubits - - @num_state_qubits.setter - def num_state_qubits(self, num_state_qubits: int) -> None: - """Set the number of state qubits. - - Note that this may change the underlying quantum register, if the number of state qubits - changes. - - Args: - num_state_qubits: The new number of qubits. - """ - if num_state_qubits != self._num_state_qubits: - self._invalidate() - self._num_state_qubits = num_state_qubits - self._reset_registers(num_state_qubits) - - @property - def main_diag(self) -> float: - """Return the entry in the main diagonal.""" - return self._main_diag - - @main_diag.setter - def main_diag(self, main_diag: float) -> None: - """Set the entry in the main diagonal. - Args: - main_diag: The new entry in the main diagonal. - """ - self._main_diag = main_diag - - @property - def off_diag(self) -> float: - """Return the entry in the off diagonals.""" - return self._off_diag - - @off_diag.setter - def off_diag(self, off_diag: float) -> None: - """Set the entry in the off diagonals. - Args: - off_diag: The new entry in the main diagonal. - """ - self._off_diag = off_diag - - @property - def tolerance(self) -> float: - """Return the error tolerance""" - return self._tolerance - - @tolerance.setter - def tolerance(self, tolerance: float) -> None: - """Set the error tolerance. - Args: - tolerance: The new error tolerance. - """ - self._tolerance = tolerance - - @property - def evolution_time(self) -> float: - """Return the time of the evolution.""" - return self._evolution_time - - @evolution_time.setter - def evolution_time(self, evolution_time: float) -> None: - """Set the time of the evolution and update the number of Trotter steps because the error - tolerance is a function of the evolution time and the number of trotter steps. - - Args: - evolution_time: The new time of the evolution. - """ - self._evolution_time = evolution_time - # Update the number of trotter steps. Max 7 for now, upper bounds too loose. - self.trotter_steps = int( - np.ceil(np.sqrt(((evolution_time * np.abs(self.off_diag)) ** 3) / 2 / self.tolerance)) - ) - - @property - def trotter_steps(self) -> int: - """Return the number of trotter steps.""" - return self._trotter_steps - - @trotter_steps.setter - def trotter_steps(self, trotter_steps: int) -> None: - """Set the number of trotter steps. - Args: - trotter_steps: The new number of trotter steps. - """ - self._trotter_steps = trotter_steps - - @property - def matrix(self) -> np.ndarray: - """Returns the tridiagonal Toeplitz matrix built according to the main and off diagonal - entries.""" - matrix = diags( - [self.off_diag, self.main_diag, self.off_diag], - [-1, 0, 1], - shape=(2**self.num_state_qubits, 2**self.num_state_qubits), - ).toarray() - return matrix - - def eigs_bounds(self) -> Tuple[float, float]: - """Return lower and upper bounds on the absolute eigenvalues of the matrix.""" - n_b = 2**self.num_state_qubits - - # Calculate minimum and maximum of absolute value of eigenvalues - # according to the formula for Toeplitz 3-diagonal matrices - - # For maximum it's enough to check border points of segment [1, n_b] - candidate_eig_ids = [1, n_b] - - # Trying to add candidates near the minimum value of absolute eigenvalues - # function abs(main_diag - 2 * off_diag * cos(i * pi / (nb + 1)) - if abs(self.main_diag) < 2 * abs(self.off_diag): - optimal_index = int(np.arccos(self.main_diag / 2 / self.off_diag) / np.pi * (n_b + 1)) - - def add_candidate_index_if_valid(index_to_add: int) -> None: - if 1 <= index_to_add <= n_b: - candidate_eig_ids.append(index_to_add) - - add_candidate_index_if_valid(optimal_index - 1) - add_candidate_index_if_valid(optimal_index) - add_candidate_index_if_valid(optimal_index + 1) - - candidate_abs_eigs = np.abs( - [ - self.main_diag - 2 * self.off_diag * np.cos(eig_id * np.pi / (n_b + 1)) - for eig_id in candidate_eig_ids - ] - ) - - lambda_min = np.min(candidate_abs_eigs) - lambda_max = np.max(candidate_abs_eigs) - return lambda_min, lambda_max - - def condition_bounds(self) -> Tuple[float, float]: - """Return lower and upper bounds on the condition number of the matrix.""" - matrix_array = self.matrix - kappa = np.linalg.cond(matrix_array) - return kappa, kappa - - def _check_configuration(self, raise_on_failure: bool = True) -> bool: - """Check if the current configuration is valid.""" - valid = True - - if self.trotter_steps < 1: - valid = False - if raise_on_failure: - raise AttributeError("The number of trotter steps should be a positive integer.") - return False - - return valid - - def _reset_registers(self, num_state_qubits: int) -> None: - """Reset the quantum registers. - - Args: - num_state_qubits: The number of qubits to represent the matrix. - """ - qr_state = QuantumRegister(num_state_qubits, "state") - self.qregs = [qr_state] - self._ancillas = [] - self._qubits = qr_state[:] - - if num_state_qubits > 1: - qr_ancilla = AncillaRegister(max(1, num_state_qubits - 1)) - self.add_register(qr_ancilla) - - def _build(self) -> None: - """If not already built, build the circuit.""" - if self._is_built: - return - - super()._build() - - self.compose(self.power(1), inplace=True) - - def _main_diag_circ(self, theta: float = 1) -> QuantumCircuit: - """Circuit implementing the matrix consisting of entries in the main diagonal. - - Args: - theta: Scale factor for the main diagonal entries (e.g. evolution_time/trotter_steps). - - Returns: - The quantum circuit implementing the matrix consisting of entries in the main diagonal. - """ - theta *= self.main_diag - qc = QuantumCircuit(self.num_state_qubits, name="main_diag") - qc.x(0) - qc.p(theta, 0) - qc.x(0) - qc.p(theta, 0) - - # pylint: disable=unused-argument - def control(num_ctrl_qubits=1, label=None, ctrl_state=None): - qc_control = QuantumCircuit(self.num_state_qubits + 1, name="main_diag") - qc_control.p(theta, 0) - return qc_control - - qc.control = control - return qc - - def _off_diag_circ(self, theta: float = 1) -> QuantumCircuit: - """Circuit implementing the matrix consisting of entries in the off diagonals. - - Args: - theta: Scale factor for the off diagonal entries (e.g. evolution_time/trotter_steps). - - Returns: - The quantum circuit implementing the matrix consisting of entries in the off diagonals. - """ - theta *= self.off_diag - - qr = QuantumRegister(self.num_state_qubits) - if self.num_state_qubits > 1: - qr_ancilla = AncillaRegister(max(1, self.num_state_qubits - 2)) - qc = QuantumCircuit(qr, qr_ancilla, name="off_diags") - else: - qc = QuantumCircuit(qr, name="off_diags") - qr_ancilla = None - - qc.u(-2 * theta, 3 * np.pi / 2, np.pi / 2, qr[0]) - - for i in range(0, self.num_state_qubits - 1): - q_controls = [] - qc.cx(qr[i], qr[i + 1]) - q_controls.append(qr[i + 1]) - - # Now we want controlled by 0 - qc.x(qr[i]) - for j in range(i, 0, -1): - qc.cx(qr[i], qr[j - 1]) - q_controls.append(qr[j - 1]) - qc.x(qr[i]) - - # Multicontrolled rotation - if len(q_controls) > 1: - ugate = UGate(-2 * theta, 3 * np.pi / 2, np.pi / 2) - qc.append( - MCMTVChain(ugate, len(q_controls), 1), - q_controls[:] + [qr[i]] + qr_ancilla[: len(q_controls) - 1], - ) - else: - qc.cu(-2 * theta, 3 * np.pi / 2, np.pi / 2, 0, q_controls[0], qr[i]) - - # Uncompute - qc.x(qr[i]) - for j in range(0, i): - qc.cx(qr[i], qr[j]) - qc.x(qr[i]) - qc.cx(qr[i], qr[i + 1]) - - # pylint: disable=unused-argument - def control(num_ctrl_qubits=1, label=None, ctrl_state=None): - qr_state = QuantumRegister(self.num_state_qubits + 1) - if self.num_state_qubits > 1: - qr_ancilla = AncillaRegister(max(1, self.num_state_qubits - 1)) - qc_control = QuantumCircuit(qr_state, qr_ancilla, name="off_diags") - else: - qc_control = QuantumCircuit(qr_state, name="off_diags") - qr_ancilla = None - # Control will be qr[0] - q_control = qr_state[0] - qr = qr_state[1:] - qc_control.cu(-2 * theta, 3 * np.pi / 2, np.pi / 2, 0, q_control, qr[0]) - - for i in range(0, self.num_state_qubits - 1): - q_controls = [] - q_controls.append(q_control) - qc_control.cx(qr[i], qr[i + 1]) - q_controls.append(qr[i + 1]) - - # Now we want controlled by 0 - qc_control.x(qr[i]) - for j in range(i, 0, -1): - qc_control.cx(qr[i], qr[j - 1]) - q_controls.append(qr[j - 1]) - qc_control.x(qr[i]) - - # Multicontrolled x rotation - if len(q_controls) > 1: - ugate = UGate(-2 * theta, 3 * np.pi / 2, np.pi / 2) - qc_control.append( - MCMTVChain(ugate, len(q_controls), 1).to_gate(), - q_controls[:] + [qr[i]] + qr_ancilla[: len(q_controls) - 1], - ) - else: - qc_control.cu(-2 * theta, 3 * np.pi / 2, np.pi / 2, 0, q_controls[0], qr[i]) - - # Uncompute - qc_control.x(qr[i]) - for j in range(0, i): - qc_control.cx(qr[i], qr[j]) - qc_control.x(qr[i]) - qc_control.cx(qr[i], qr[i + 1]) - return qc_control - - qc.control = control - return qc - - def inverse(self): - return TridiagonalToeplitz( - self.num_state_qubits, - self.main_diag, - self.off_diag, - evolution_time=-1 * self.evolution_time, - ) - - def power(self, power: int, matrix_power: bool = False) -> QuantumCircuit: - """Build powers of the circuit. - - Args: - power: The power to raise this circuit to. - matrix_power: If True, the circuit is converted to a matrix and then the - matrix power is computed. If False, and ``power`` is a positive integer, - the implementation defaults to ``repeat``. - - Returns: - The quantum circuit implementing powers of the unitary. - """ - qc_raw = QuantumCircuit(self.num_state_qubits) - - # pylint: disable=unused-argument - def control(num_ctrl_qubits=1, label=None, ctrl_state=None): - qr_state = QuantumRegister(self.num_state_qubits + 1, "state") - if self.num_state_qubits > 1: - qr_ancilla = AncillaRegister(max(1, self.num_state_qubits - 1)) - qc = QuantumCircuit(qr_state, qr_ancilla, name="exp(iHk)") - else: - qc = QuantumCircuit(qr_state, name="exp(iHk)") - qr_ancilla = None - # Control will be qr[0] - q_control = qr_state[0] - qr = qr_state[1:] - # A1 commutes, so one application with evolution_time*2^{j} to the last qubit is enough - qc.append( - self._main_diag_circ(self.evolution_time * power).control().to_gate(), - [q_control] + qr[:], - ) - - # Update trotter steps to compensate the error - trotter_steps_new = int(np.ceil(np.sqrt(power) * self.trotter_steps)) - - # exp(iA2t/2m) - qc.u( - self.off_diag * self.evolution_time * power / trotter_steps_new, - 3 * np.pi / 2, - np.pi / 2, - qr[0], - ) - # for _ in range(power): - for _ in range(0, trotter_steps_new): - if qr_ancilla: - qc.append( - self._off_diag_circ(self.evolution_time * power / trotter_steps_new) - .control() - .to_gate(), - [q_control] + qr[:] + qr_ancilla[:], - ) - else: - qc.append( - self._off_diag_circ(self.evolution_time * power / trotter_steps_new) - .control() - .to_gate(), - [q_control] + qr[:], - ) - # exp(-iA2t/2m) - qc.u( - -self.off_diag * self.evolution_time * power / trotter_steps_new, - 3 * np.pi / 2, - np.pi / 2, - qr[0], - ) - return qc - - qc_raw.control = control - return qc_raw diff --git a/qiskit/algorithms/linear_solvers/numpy_linear_solver.py b/qiskit/algorithms/linear_solvers/numpy_linear_solver.py deleted file mode 100644 index b0f6f80c200b..000000000000 --- a/qiskit/algorithms/linear_solvers/numpy_linear_solver.py +++ /dev/null @@ -1,109 +0,0 @@ -# This code is part of Qiskit. -# -# (C) Copyright IBM 2020, 2022. -# -# 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 Numpy LinearSolver algorithm (classical).""" - -from typing import List, Union, Optional, Callable -import numpy as np - -from qiskit import QuantumCircuit -from qiskit.quantum_info import Operator, Statevector -from qiskit.quantum_info.operators.base_operator import BaseOperator -from qiskit.utils.deprecation import deprecate_function - -from .linear_solver import LinearSolverResult, LinearSolver -from .observables.linear_system_observable import LinearSystemObservable - - -class NumPyLinearSolver(LinearSolver): - """The deprecated Numpy Linear Solver algorithm (classical). - - This linear system solver computes the exact value of the given observable(s) or the full - solution vector if no observable is specified. - - Examples:: - - import warnings - import numpy as np - from qiskit.algorithms import NumPyLinearSolver - from qiskit.algorithms.linear_solvers.matrices import TridiagonalToeplitz - from qiskit.algorithms.linear_solvers.observables import MatrixFunctional - - with warnings.catch_warnings(): - warnings.simplefilter('ignore') - matrix = TridiagonalToeplitz(2, 1, 1 / 3, trotter_steps=2) - right_hand_side = [1.0, -2.1, 3.2, -4.3] - observable = MatrixFunctional(1, 1 / 2) - rhs = right_hand_side / np.linalg.norm(right_hand_side) - - with warnings.catch_warnings(): - warnings.simplefilter('ignore') - np_solver = NumPyLinearSolver() - solution = np_solver.solve(matrix, rhs, observable) - result = solution.observable - """ - - @deprecate_function( - "The NumPyLinearSolver class is deprecated as of Qiskit Terra 0.22.0 " - "and will be removed no sooner than 3 months after the release date. ", - since="0.22.0", - ) - def __init__(self) -> None: - super().__init__() - - def solve( - self, - matrix: Union[np.ndarray, QuantumCircuit], - vector: Union[np.ndarray, QuantumCircuit], - observable: Optional[ - Union[LinearSystemObservable, BaseOperator, List[BaseOperator]] - ] = None, - observable_circuit: Optional[Union[QuantumCircuit, List[QuantumCircuit]]] = None, - post_processing: Optional[ - Callable[[Union[float, List[float]]], Union[float, List[float]]] - ] = None, - ) -> LinearSolverResult: - """Solve classically the linear system and compute the observable(s) - - Args: - matrix: The matrix specifying the system, i.e. A in Ax=b. - vector: The vector specifying the right hand side of the equation in Ax=b. - observable: Optional information to be extracted from the solution. - Default is the probability of success of the algorithm. - observable_circuit: Optional circuit to be applied to the solution to extract - information. Default is ``None``. - post_processing: Optional function to compute the value of the observable. - Default is the raw value of measuring the observable. - - Returns: - The result of the linear system. - """ - # Check if either matrix or vector are QuantumCircuits and get the array from them - if isinstance(vector, QuantumCircuit): - vector = Statevector(vector).data - if isinstance(matrix, QuantumCircuit): - if hasattr(matrix, "matrix"): - matrix = matrix.matrix - else: - matrix = Operator(matrix).data - - solution_vector = np.linalg.solve(matrix, vector) - solution = LinearSolverResult() - solution.state = solution_vector - if observable is not None: - if isinstance(observable, list): - solution.observable = [] - for obs in observable: - solution.observable.append(obs.evaluate_classically(solution_vector)) - else: - solution.observable = observable.evaluate_classically(solution_vector) - solution.euclidean_norm = np.linalg.norm(solution_vector) - return solution diff --git a/qiskit/algorithms/linear_solvers/observables/__init__.py b/qiskit/algorithms/linear_solvers/observables/__init__.py deleted file mode 100644 index f4ff2b7c9aea..000000000000 --- a/qiskit/algorithms/linear_solvers/observables/__init__.py +++ /dev/null @@ -1,19 +0,0 @@ -# This code is part of Qiskit. -# -# (C) Copyright IBM 2020, 2021. -# -# 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. - -"""Observables for Qiskit's linear solvers.""" - -from .linear_system_observable import LinearSystemObservable -from .absolute_average import AbsoluteAverage -from .matrix_functional import MatrixFunctional - -__all__ = ["LinearSystemObservable", "AbsoluteAverage", "MatrixFunctional"] diff --git a/qiskit/algorithms/linear_solvers/observables/absolute_average.py b/qiskit/algorithms/linear_solvers/observables/absolute_average.py deleted file mode 100644 index 2a7bf25fd11f..000000000000 --- a/qiskit/algorithms/linear_solvers/observables/absolute_average.py +++ /dev/null @@ -1,133 +0,0 @@ -# This code is part of Qiskit. -# -# (C) Copyright IBM 2020, 2022. -# -# 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 absolute value of the average of a linear system of equations solution.""" - -from typing import Union, List -import numpy as np - -from qiskit import QuantumCircuit -from qiskit.opflow import I, Z, TensoredOp -from qiskit.quantum_info import Statevector -from qiskit.utils.deprecation import deprecate_function - -from .linear_system_observable import LinearSystemObservable - - -class AbsoluteAverage(LinearSystemObservable): - r"""The deprecated observable for the absolute average of a linear system of equations solution. - - For a vector :math:`x=(x_1,...,x_N)`, the absolute average is defined as - :math:`\abs{\frac{1}{N}\sum_{i=1}^{N}x_i}`. - - Examples:: - - import warnings - import numpy as np - from qiskit import QuantumCircuit - from qiskit.algorithms.linear_solvers.observables.absolute_average import \ - AbsoluteAverage - from qiskit.opflow import StateFn - - with warnings.catch_warnings(): - warnings.simplefilter('ignore') - observable = AbsoluteAverage() - vector = [1.0, -2.1, 3.2, -4.3] - - init_state = vector / np.linalg.norm(vector) - num_qubits = int(np.log2(len(vector))) - - qc = QuantumCircuit(num_qubits) - qc.isometry(init_state, list(range(num_qubits)), None) - qc.append(observable.observable_circuit(num_qubits), list(range(num_qubits))) - - # Observable operator - observable_op = observable.observable(num_qubits) - state_vec = (~StateFn(observable_op) @ StateFn(qc)).eval() - - # Obtain result - result = observable.post_processing(state_vec, num_qubits) - - # Obtain analytical evaluation - exact = observable.evaluate_classically(init_state) - """ - - @deprecate_function( - "The AbsoluteAverage class is deprecated as of Qiskit Terra 0.22.0 " - "and will be removed no sooner than 3 months after the release date.", - since="0.22.0", - ) - def __init__(self) -> None: - super().__init__() - - def observable(self, num_qubits: int) -> Union[TensoredOp, List[TensoredOp]]: - """The observable operator. - - Args: - num_qubits: The number of qubits on which the observable will be applied. - - Returns: - The observable as a sum of Pauli strings. - """ - zero_op = (I + Z) / 2 - return TensoredOp(num_qubits * [zero_op]) - - def observable_circuit(self, num_qubits: int) -> Union[QuantumCircuit, List[QuantumCircuit]]: - """The circuit implementing the absolute average observable. - - Args: - num_qubits: The number of qubits on which the observable will be applied. - - Returns: - The observable as a QuantumCircuit. - """ - qc = QuantumCircuit(num_qubits) - qc.h(qc.qubits) - return qc - - def post_processing( - self, solution: Union[float, List[float]], num_qubits: int, scaling: float = 1 - ) -> float: - """Evaluates the absolute average on the solution to the linear system. - - Args: - solution: The probability calculated from the circuit and the observable. - num_qubits: The number of qubits where the observable was applied. - scaling: Scaling of the solution. - - Returns: - The value of the absolute average. - - Raises: - ValueError: If the input is not in the correct format. - """ - if isinstance(solution, list): - if len(solution) == 1: - solution = solution[0] - else: - raise ValueError("Solution probability must be given as a single value.") - - return np.real(np.sqrt(solution / (2**num_qubits)) / scaling) - - def evaluate_classically(self, solution: Union[np.array, QuantumCircuit]) -> float: - """Evaluates the given observable on the solution to the linear system. - - Args: - solution: The solution to the system as a numpy array or the circuit that prepares it. - - Returns: - The value of the observable. - """ - # Check if it is QuantumCircuits and get the array from them - if isinstance(solution, QuantumCircuit): - solution = Statevector(solution).data - return np.abs(np.mean(solution)) diff --git a/qiskit/algorithms/linear_solvers/observables/linear_system_observable.py b/qiskit/algorithms/linear_solvers/observables/linear_system_observable.py deleted file mode 100644 index 5cb6557b7e4c..000000000000 --- a/qiskit/algorithms/linear_solvers/observables/linear_system_observable.py +++ /dev/null @@ -1,86 +0,0 @@ -# This code is part of Qiskit. -# -# (C) Copyright IBM 2020, 2022. -# -# 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. - -"""An abstract class for linear systems solvers in Qiskit's aqua module.""" - -from abc import ABC, abstractmethod -from typing import Union, List -import numpy as np - -from qiskit import QuantumCircuit -from qiskit.opflow import TensoredOp -from qiskit.utils.deprecation import deprecate_function - - -class LinearSystemObservable(ABC): - """The deprecated abstract class for linear system observables in Qiskit.""" - - @deprecate_function( - "The LinearSystemObservable class is deprecated as of Qiskit Terra 0.22.0 " - "and will be removed no sooner than 3 months after the release date.", - since="0.22.0", - ) - def __init__(self) -> None: - pass - - @abstractmethod - def observable(self, num_qubits: int) -> Union[TensoredOp, List[TensoredOp]]: - """The observable operator. - - Args: - num_qubits: The number of qubits on which the observable will be applied. - - Returns: - The observable as a sum of Pauli strings. - """ - raise NotImplementedError - - @abstractmethod - def observable_circuit(self, num_qubits: int) -> Union[QuantumCircuit, List[QuantumCircuit]]: - """The circuit implementing the observable. - - Args: - num_qubits: The number of qubits on which the observable will be applied. - - Returns: - The observable as a QuantumCircuit. - """ - raise NotImplementedError - - @abstractmethod - def post_processing( - self, solution: Union[float, List[float]], num_qubits: int, scaling: float = 1 - ) -> float: - """Evaluates the given observable on the solution to the linear system. - - Args: - solution: The probability calculated from the circuit and the observable. - num_qubits: The number of qubits where the observable was applied. - scaling: Scaling of the solution. - - Returns: - The value of the observable. - """ - raise NotImplementedError - - @abstractmethod - def evaluate_classically(self, solution: Union[np.array, QuantumCircuit]) -> float: - """Calculates the analytical value of the given observable from the solution vector to the - linear system. - - Args: - solution: The solution to the system as a numpy array or the circuit that prepares it. - - Returns: - The value of the observable. - """ - raise NotImplementedError diff --git a/qiskit/algorithms/linear_solvers/observables/matrix_functional.py b/qiskit/algorithms/linear_solvers/observables/matrix_functional.py deleted file mode 100644 index 03caf9ecc7ec..000000000000 --- a/qiskit/algorithms/linear_solvers/observables/matrix_functional.py +++ /dev/null @@ -1,190 +0,0 @@ -# This code is part of Qiskit. -# -# (C) Copyright IBM 2020, 2022. -# -# 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 matrix functional of the vector solution to the linear systems.""" - -from typing import Union, List -import warnings -import numpy as np -from scipy.sparse import diags - -from qiskit import QuantumCircuit -from qiskit.quantum_info import Statevector -from qiskit.opflow import I, Z, TensoredOp -from qiskit.utils.deprecation import deprecate_function - -from .linear_system_observable import LinearSystemObservable - - -class MatrixFunctional(LinearSystemObservable): - """The deprecated class for the matrix functional of the vector solution to the linear systems. - - Examples:: - - import warnings - import numpy as np - from qiskit import QuantumCircuit - from qiskit.algorithms.linear_solvers.observables.matrix_functional import \ - MatrixFunctional - from qiskit.transpiler.passes import RemoveResetInZeroState - from qiskit.opflow import StateFn - - tpass = RemoveResetInZeroState() - - vector = [1.0, -2.1, 3.2, -4.3] - with warnings.catch_warnings(): - warnings.simplefilter('ignore') - observable = MatrixFunctional(1, -1 / 3) - - init_state = vector / np.linalg.norm(vector) - num_qubits = int(np.log2(len(vector))) - - # Get observable circuits - obs_circuits = observable.observable_circuit(num_qubits) - qcs = [] - for obs_circ in obs_circuits: - qc = QuantumCircuit(num_qubits) - qc.isometry(init_state, list(range(num_qubits)), None) - qc.append(obs_circ, list(range(num_qubits))) - qcs.append(tpass(qc.decompose())) - - # Get observables - observable_ops = observable.observable(num_qubits) - state_vecs = [] - # First is the norm - state_vecs.append((~StateFn(observable_ops[0]) @ StateFn(qcs[0])).eval()) - for i in range(1, len(observable_ops), 2): - state_vecs += [(~StateFn(observable_ops[i]) @ StateFn(qcs[i])).eval(), - (~StateFn(observable_ops[i + 1]) @ StateFn(qcs[i + 1])).eval()] - - # Obtain result - result = observable.post_processing(state_vecs, num_qubits) - - # Obtain analytical evaluation - exact = observable.evaluate_classically(init_state) - """ - - @deprecate_function( - "The MatrixFunctional class is deprecated as of Qiskit Terra 0.22.0 " - "and will be removed no sooner than 3 months after the release date.", - since="0.22.0", - ) - def __init__(self, main_diag: float, off_diag: int) -> None: - """ - Args: - main_diag: The main diagonal of the tridiagonal Toeplitz symmetric matrix to compute - the functional. - off_diag: The off diagonal of the tridiagonal Toeplitz symmetric matrix to compute - the functional. - """ - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - super().__init__() - self._main_diag = main_diag - self._off_diag = off_diag - - def observable(self, num_qubits: int) -> Union[TensoredOp, List[TensoredOp]]: - """The observable operators. - - Args: - num_qubits: The number of qubits on which the observable will be applied. - - Returns: - The observable as a list of sums of Pauli strings. - """ - zero_op = (I + Z) / 2 - one_op = (I - Z) / 2 - observables = [] - # First we measure the norm of x - observables.append(I ^ num_qubits) - for i in range(num_qubits): - j = num_qubits - i - 1 - - # TODO this if can be removed once the bug in Opflow is fixed where - # TensoredOp([X, TensoredOp([])]).eval() ends up in infinite recursion - if i > 0: - observables += [ - (I ^ j) ^ zero_op ^ TensoredOp(i * [one_op]), - (I ^ j) ^ one_op ^ TensoredOp(i * [one_op]), - ] - else: - observables += [(I ^ j) ^ zero_op, (I ^ j) ^ one_op] - - return observables - - def observable_circuit(self, num_qubits: int) -> Union[QuantumCircuit, List[QuantumCircuit]]: - """The circuits to implement the matrix functional observable. - - Args: - num_qubits: The number of qubits on which the observable will be applied. - - Returns: - The observable as a list of QuantumCircuits. - """ - qcs = [] - # Again, the first value in the list will correspond to the norm of x - qcs.append(QuantumCircuit(num_qubits)) - for i in range(0, num_qubits): - qc = QuantumCircuit(num_qubits) - for j in range(0, i): - qc.cx(i, j) - qc.h(i) - qcs += [qc, qc] - - return qcs - - def post_processing( - self, solution: Union[float, List[float]], num_qubits: int, scaling: float = 1 - ) -> float: - """Evaluates the matrix functional on the solution to the linear system. - - Args: - solution: The list of probabilities calculated from the circuit and the observable. - num_qubits: The number of qubits where the observable was applied. - scaling: Scaling of the solution. - - Returns: - The value of the absolute average. - - Raises: - ValueError: If the input is not in the correct format. - """ - if not isinstance(solution, list): - raise ValueError("Solution probabilities must be given in list form.") - - # Calculate the value from the off-diagonal elements - off_val = 0 - for i in range(1, len(solution), 2): - off_val += (solution[i] - solution[i + 1]) / (scaling**2) - main_val = solution[0] / (scaling**2) - return np.real(self._main_diag * main_val + self._off_diag * off_val) - - def evaluate_classically(self, solution: Union[np.array, QuantumCircuit]) -> float: - """Evaluates the given observable on the solution to the linear system. - - Args: - solution: The solution to the system as a numpy array or the circuit that prepares it. - - Returns: - The value of the observable. - """ - # Check if it is QuantumCircuits and get the array from them - if isinstance(solution, QuantumCircuit): - solution = Statevector(solution).data - - matrix = diags( - [self._off_diag, self._main_diag, self._off_diag], - [-1, 0, 1], - shape=(len(solution), len(solution)), - ).toarray() - - return np.dot(solution.transpose(), np.dot(matrix, solution)) diff --git a/releasenotes/notes/remove-deprecated-factorizers-linear-solvers-4631870129749624.yaml b/releasenotes/notes/remove-deprecated-factorizers-linear-solvers-4631870129749624.yaml new file mode 100644 index 000000000000..330a99ef9f22 --- /dev/null +++ b/releasenotes/notes/remove-deprecated-factorizers-linear-solvers-4631870129749624.yaml @@ -0,0 +1,10 @@ +--- +upgrade: + - | + The deprecated modules ``factorizers`` and ``linear_solvers``, containing + ``HHL`` and ``Shor`` have been removed from + :mod:`qiskit.algorithms`. These functionalities + were originally deprecated as part of the 0.22.0 release (released on + October 13, 2022). You can access the code through the Qiskit Textbook instead: + `Linear Solvers (HHL) `_ , + `Factorizers (Shor) `_ diff --git a/test/python/algorithms/test_backendv1.py b/test/python/algorithms/test_backendv1.py index 67514cbdc327..ed4bf878495a 100644 --- a/test/python/algorithms/test_backendv1.py +++ b/test/python/algorithms/test_backendv1.py @@ -13,17 +13,15 @@ """Test Providers that support BackendV1 interface""" import unittest -import warnings from test.python.algorithms import QiskitAlgorithmsTestCase from qiskit import QuantumCircuit from qiskit.providers.fake_provider import FakeProvider from qiskit.utils import QuantumInstance, algorithm_globals -from qiskit.algorithms import Shor, VQE, Grover, AmplificationProblem +from qiskit.algorithms import VQE, Grover, AmplificationProblem from qiskit.opflow import X, Z, I from qiskit.algorithms.optimizers import SPSA from qiskit.circuit.library import TwoLocal, EfficientSU2 from qiskit.utils.mitigation import CompleteMeasFitter -from qiskit.test import slow_test class TestBackendV1(QiskitAlgorithmsTestCase): @@ -35,26 +33,6 @@ def setUp(self): self._qasm = self._provider.get_backend("fake_qasm_simulator") self.seed = 50 - @slow_test - def test_shor_factoring(self): - """shor factoring test""" - n_v = 15 - factors = [3, 5] - qasm_simulator = QuantumInstance( - self._qasm, shots=1000, seed_simulator=self.seed, seed_transpiler=self.seed - ) - with warnings.catch_warnings(record=True) as caught_warnings: - warnings.filterwarnings( - "always", - category=DeprecationWarning, - ) - shor = Shor(quantum_instance=qasm_simulator) - self.assertTrue("Shor class is deprecated" in str(caught_warnings[0].message)) - - result = shor.factor(N=n_v) - self.assertListEqual(result.factors[0], factors) - self.assertTrue(result.total_counts >= result.successful_counts) - def test_vqe_qasm(self): """Test the VQE on QASM simulator.""" h2_op = ( diff --git a/test/python/algorithms/test_backendv2.py b/test/python/algorithms/test_backendv2.py index 3945c82f8f01..71712a4a5c7e 100644 --- a/test/python/algorithms/test_backendv2.py +++ b/test/python/algorithms/test_backendv2.py @@ -13,17 +13,15 @@ """Test Providers that support BackendV2 interface""" import unittest -import warnings from test.python.algorithms import QiskitAlgorithmsTestCase from qiskit import QuantumCircuit from qiskit.providers.fake_provider import FakeProvider from qiskit.providers.fake_provider.fake_backend_v2 import FakeBackendSimple from qiskit.utils import QuantumInstance -from qiskit.algorithms import Shor, VQE, Grover, AmplificationProblem +from qiskit.algorithms import VQE, Grover, AmplificationProblem from qiskit.opflow import X, Z, I from qiskit.algorithms.optimizers import SPSA from qiskit.circuit.library import TwoLocal -from qiskit.test import slow_test class TestBackendV2(QiskitAlgorithmsTestCase): @@ -35,26 +33,6 @@ def setUp(self): self._qasm = FakeBackendSimple() self.seed = 50 - @slow_test - def test_shor_factoring(self): - """shor factoring test""" - n_v = 15 - factors = [3, 5] - qasm_simulator = QuantumInstance( - self._qasm, shots=1000, seed_simulator=self.seed, seed_transpiler=self.seed - ) - with warnings.catch_warnings(record=True) as caught_warnings: - warnings.filterwarnings( - "always", - category=DeprecationWarning, - ) - shor = Shor(quantum_instance=qasm_simulator) - self.assertTrue("Shor class is deprecated" in str(caught_warnings[0].message)) - - result = shor.factor(N=n_v) - self.assertListEqual(result.factors[0], factors) - self.assertTrue(result.total_counts >= result.successful_counts) - def test_vqe_qasm(self): """Test the VQE on QASM simulator.""" h2_op = ( diff --git a/test/python/algorithms/test_linear_solvers.py b/test/python/algorithms/test_linear_solvers.py deleted file mode 100644 index c2cb747e2799..000000000000 --- a/test/python/algorithms/test_linear_solvers.py +++ /dev/null @@ -1,358 +0,0 @@ -# This code is part of Qiskit. -# -# (C) Copyright IBM 2020, 2022. -# -# 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 linear system solver algorithm.""" - -import unittest -import warnings -from test.python.algorithms import QiskitAlgorithmsTestCase -from scipy.linalg import expm -import numpy as np -from ddt import ddt, idata, unpack -from qiskit import BasicAer, QuantumCircuit -from qiskit.algorithms.linear_solvers.hhl import HHL -from qiskit.algorithms.linear_solvers.matrices.tridiagonal_toeplitz import TridiagonalToeplitz -from qiskit.algorithms.linear_solvers.matrices.numpy_matrix import NumPyMatrix -from qiskit.algorithms.linear_solvers.observables.absolute_average import AbsoluteAverage -from qiskit.algorithms.linear_solvers.observables.matrix_functional import MatrixFunctional -from qiskit.circuit.library.arithmetic.exact_reciprocal import ExactReciprocal -from qiskit.quantum_info import Operator, partial_trace -from qiskit.opflow import I, Z, StateFn -from qiskit.utils import QuantumInstance -from qiskit import quantum_info -from qiskit.test import slow_test - - -def _factory_tridiagonal_toeplitz( - num_state_qubits: int, main_diag: float, off_diag: float, trotter_steps: int = 1 -): - with warnings.catch_warnings(record=True): - warnings.simplefilter("ignore") - return TridiagonalToeplitz( - num_state_qubits, main_diag, off_diag, trotter_steps=trotter_steps - ) - - -def _factory_numpy_matrix(matrix: np.ndarray): - with warnings.catch_warnings(record=True): - warnings.simplefilter("ignore") - return NumPyMatrix(matrix) - - -@ddt -class TestMatrices(QiskitAlgorithmsTestCase): - """Tests based on the matrices classes. - - This class tests - * the constructed circuits - """ - - @idata( - [ - [_factory_tridiagonal_toeplitz(2, 1, -1 / 3)], - [_factory_tridiagonal_toeplitz(3, 2, 1), 1.1, 3], - [ - _factory_numpy_matrix( - np.array( - [ - [1 / 2, 1 / 6, 0, 0], - [1 / 6, 1 / 2, 1 / 6, 0], - [0, 1 / 6, 1 / 2, 1 / 6], - [0, 0, 1 / 6, 1 / 2], - ] - ) - ) - ], - ] - ) - @unpack - def test_matrices(self, matrix, time=1.0, power=1): - """Test the different matrix classes.""" - matrix.evolution_time = time - - num_qubits = matrix.num_state_qubits - pow_circ = matrix.power(power).control() - circ_qubits = pow_circ.num_qubits - qc = QuantumCircuit(circ_qubits) - qc.append(matrix.power(power).control(), list(range(circ_qubits))) - # extract the parts of the circuit matrix corresponding to TridiagonalToeplitz - zero_op = (I + Z) / 2 - one_op = (I - Z) / 2 - proj = Operator((zero_op ^ pow_circ.num_ancillas) ^ (I ^ num_qubits) ^ one_op).data - circ_matrix = Operator(qc).data - approx_exp = partial_trace( - np.dot(proj, circ_matrix), [0] + list(range(num_qubits + 1, circ_qubits)) - ).data - - exact_exp = expm(1j * matrix.evolution_time * power * matrix.matrix) - np.testing.assert_array_almost_equal(approx_exp, exact_exp, decimal=2) - - @idata( - [ - [_factory_tridiagonal_toeplitz(2, 1.5, 2.5)], - [_factory_tridiagonal_toeplitz(4, -1, 1.6)], - ] - ) - @unpack - def test_eigs_bounds(self, matrix): - """Test the capability of TridiagonalToeplitz matrix class - to find accurate absolute eigenvalues bounds.""" - - matrix_lambda_min, matrix_lambda_max = matrix.eigs_bounds() - - numpy_matrix = matrix.matrix - eigenvalues, _ = np.linalg.eig(numpy_matrix) - abs_eigenvalues = np.abs(eigenvalues) - exact_lambda_min = np.min(abs_eigenvalues) - exact_lambda_max = np.max(abs_eigenvalues) - - np.testing.assert_almost_equal(matrix_lambda_min, exact_lambda_min, decimal=6) - np.testing.assert_almost_equal(matrix_lambda_max, exact_lambda_max, decimal=6) - - -def _factory_absolute_average(): - with warnings.catch_warnings(record=True): - warnings.simplefilter("ignore") - return AbsoluteAverage() - - -def _factory_matrix_functional(main_diag: float, off_diag: int): - with warnings.catch_warnings(record=True): - warnings.simplefilter("ignore") - return MatrixFunctional(main_diag, off_diag) - - -@ddt -class TestObservables(QiskitAlgorithmsTestCase): - """Tests based on the observables classes. - - This class tests - * the constructed circuits - """ - - @idata( - [ - [_factory_absolute_average(), [1.0, -2.1, 3.2, -4.3]], - [_factory_absolute_average(), [-9 / 4, -0.3, 8 / 7, 10, -5, 11.1, 13 / 11, -27 / 12]], - ] - ) - @unpack - def test_absolute_average(self, observable, vector): - """Test the absolute average observable.""" - init_state = vector / np.linalg.norm(vector) - num_qubits = int(np.log2(len(vector))) - - qc = QuantumCircuit(num_qubits) - qc.isometry(init_state, list(range(num_qubits)), None) - qc.append(observable.observable_circuit(num_qubits), list(range(num_qubits))) - - # Observable operator - observable_op = observable.observable(num_qubits) - state_vec = (~StateFn(observable_op) @ StateFn(qc)).eval() - - # Obtain result - result = observable.post_processing(state_vec, num_qubits) - - # Obtain analytical evaluation - exact = observable.evaluate_classically(init_state) - - np.testing.assert_almost_equal(result, exact, decimal=2) - - @idata( - [ - [_factory_matrix_functional(1, -1 / 3), [1.0, -2.1, 3.2, -4.3]], - [ - _factory_matrix_functional(2 / 3, 11 / 7), - [-9 / 4, -0.3, 8 / 7, 10, -5, 11.1, 13 / 11, -27 / 12], - ], - ] - ) - @unpack - def test_matrix_functional(self, observable, vector): - """Test the matrix functional class.""" - from qiskit.transpiler.passes import RemoveResetInZeroState - - tpass = RemoveResetInZeroState() - - init_state = vector / np.linalg.norm(vector) - num_qubits = int(np.log2(len(vector))) - - # Get observable circuits - obs_circuits = observable.observable_circuit(num_qubits) - qcs = [] - for obs_circ in obs_circuits: - qc = QuantumCircuit(num_qubits) - qc.isometry(init_state, list(range(num_qubits)), None) - qc.append(obs_circ, list(range(num_qubits))) - qcs.append(tpass(qc.decompose())) - - # Get observables - observable_ops = observable.observable(num_qubits) - state_vecs = [] - # First is the norm - state_vecs.append((~StateFn(observable_ops[0]) @ StateFn(qcs[0])).eval()) - for i in range(1, len(observable_ops), 2): - state_vecs += [ - (~StateFn(observable_ops[i]) @ StateFn(qcs[i])).eval(), - (~StateFn(observable_ops[i + 1]) @ StateFn(qcs[i + 1])).eval(), - ] - - # Obtain result - result = observable.post_processing(state_vecs, num_qubits) - - # Obtain analytical evaluation - exact = observable.evaluate_classically(init_state) - - np.testing.assert_almost_equal(result, exact, decimal=2) - - -@ddt -class TestReciprocal(QiskitAlgorithmsTestCase): - """Tests based on the reciprocal classes. - - This class tests - * the constructed circuits - """ - - @idata([[2, 0.1, False], [3, 1 / 9, True]]) - @unpack - def test_exact_reciprocal(self, num_qubits, scaling, neg_vals): - """Test the ExactReciprocal class.""" - reciprocal = ExactReciprocal(num_qubits + neg_vals, scaling, neg_vals) - - qc = QuantumCircuit(num_qubits + 1 + neg_vals) - qc.h(list(range(num_qubits))) - # If negative eigenvalues, set the sign qubit to 1 - if neg_vals: - qc.x(num_qubits) - qc.append(reciprocal, list(range(num_qubits + 1 + neg_vals))) - - # Create the operator 0 - state_vec = quantum_info.Statevector.from_instruction(qc).data[-(2**num_qubits) :] - - # Remove the factor from the hadamards - state_vec *= np.sqrt(2) ** num_qubits - - # Analytic value - exact = [] - for i in range(0, 2**num_qubits): - if i == 0: - exact.append(0) - else: - if neg_vals: - exact.append(-scaling / (1 - i / (2**num_qubits))) - else: - exact.append(scaling * (2**num_qubits) / i) - - np.testing.assert_array_almost_equal(state_vec, exact, decimal=2) - - -@ddt -class TestLinearSolver(QiskitAlgorithmsTestCase): - """Tests based on the linear solvers classes. - - This class tests - * the constructed circuits - """ - - @slow_test - @idata( - [ - [ - _factory_tridiagonal_toeplitz(2, 1, 1 / 3, trotter_steps=2), - [1.0, -2.1, 3.2, -4.3], - _factory_matrix_functional(1, 1 / 2), - ], - [ - np.array( - [[0, 0, 1.585, 0], [0, 0, -0.585, 1], [1.585, -0.585, 0, 0], [0, 1, 0, 0]] - ), - [1.0, 0, 0, 0], - _factory_matrix_functional(1, 1 / 2), - ], - [ - [ - [1 / 2, 1 / 6, 0, 0], - [1 / 6, 1 / 2, 1 / 6, 0], - [0, 1 / 6, 1 / 2, 1 / 6], - [0, 0, 1 / 6, 1 / 2], - ], - [1.0, -2.1, 3.2, -4.3], - _factory_matrix_functional(1, 1 / 2), - ], - [ - np.array([[82, 34], [34, 58]]), - np.array([[1], [0]]), - _factory_absolute_average(), - 3, - ], - [ - _factory_tridiagonal_toeplitz(3, 1, -1 / 2, trotter_steps=2), - [-9 / 4, -0.3, 8 / 7, 10, -5, 11.1, 13 / 11, -27 / 12], - _factory_absolute_average(), - ], - ] - ) - @unpack - def test_hhl(self, matrix, right_hand_side, observable, decimal=1): - """Test the HHL class.""" - if isinstance(matrix, QuantumCircuit): - num_qubits = matrix.num_state_qubits - elif isinstance(matrix, (np.ndarray)): - num_qubits = int(np.log2(matrix.shape[0])) - elif isinstance(matrix, list): - num_qubits = int(np.log2(len(matrix))) - - rhs = right_hand_side / np.linalg.norm(right_hand_side) - - # Initial state circuit - qc = QuantumCircuit(num_qubits) - qc.isometry(rhs, list(range(num_qubits)), None) - - with warnings.catch_warnings(record=True) as caught_warnings: - warnings.simplefilter("always") - hhl = HHL() - self.assertTrue("HHL class is deprecated" in str(caught_warnings[0].message)) - solution = hhl.solve(matrix, qc, observable) - approx_result = solution.observable - - # Calculate analytical value - if isinstance(matrix, QuantumCircuit): - exact_x = np.dot(np.linalg.inv(matrix.matrix), rhs) - elif isinstance(matrix, (list, np.ndarray)): - if isinstance(matrix, list): - matrix = np.array(matrix) - exact_x = np.dot(np.linalg.inv(matrix), rhs) - exact_result = observable.evaluate_classically(exact_x) - - np.testing.assert_almost_equal(approx_result, exact_result, decimal=decimal) - - def test_hhl_qi(self): - """Test the HHL quantum instance getter and setter.""" - with warnings.catch_warnings(record=True) as caught_warnings: - warnings.simplefilter("always") - hhl = HHL() - self.assertTrue("HHL class is deprecated" in str(caught_warnings[0].message)) - self.assertIsNone(hhl.quantum_instance) # Defaults to None - - # First set a valid quantum instance and check via getter - qinst = QuantumInstance(backend=BasicAer.get_backend("qasm_simulator")) - hhl.quantum_instance = qinst - self.assertEqual(hhl.quantum_instance, qinst) - - # Now set quantum instance back to None and check via getter - hhl.quantum_instance = None - self.assertIsNone(hhl.quantum_instance) - - -if __name__ == "__main__": - unittest.main() diff --git a/test/python/algorithms/test_shor.py b/test/python/algorithms/test_shor.py deleted file mode 100644 index 742f7eb1dfdd..000000000000 --- a/test/python/algorithms/test_shor.py +++ /dev/null @@ -1,199 +0,0 @@ -# This code is part of Qiskit. -# -# (C) Copyright IBM 2018, 2022. -# -# 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 Shor""" - -import unittest -import warnings -import math -from test.python.algorithms import QiskitAlgorithmsTestCase -from ddt import ddt, data, idata, unpack - -from qiskit import ClassicalRegister -from qiskit.utils import QuantumInstance, optionals -from qiskit.algorithms import Shor -from qiskit.test import slow_test - - -@ddt -class TestShor(QiskitAlgorithmsTestCase): - """test Shor's algorithm""" - - @unittest.skipUnless(optionals.HAS_AER, "qiskit-aer is required to run this test") - def setUp(self): - super().setUp() - from qiskit_aer import Aer - - backend = Aer.get_backend("aer_simulator") - with warnings.catch_warnings(record=True) as caught_warnings: - warnings.filterwarnings( - "always", - category=DeprecationWarning, - ) - self.instance = Shor(quantum_instance=QuantumInstance(backend, shots=1000)) - self.assertTrue("Shor class is deprecated" in str(caught_warnings[0].message)) - - @slow_test - @idata( - [ - [15, "aer_simulator", [3, 5]], - ] - ) - @unpack - def test_shor_factoring(self, n_v, backend, factors): - """shor factoring test for n = log(N) = 4""" - self._test_shor_factoring(backend, factors, n_v) - - @slow_test - @idata( - [ - [21, "aer_simulator", [3, 7]], - ] - ) - @unpack - def test_shor_factoring_5_bit_number(self, n_v, backend, factors): - """shor factoring test for n = log(N) = 5""" - self._test_shor_factoring(backend, factors, n_v) - - def _test_shor_factoring(self, backend, factors, n_v): - """shor factoring test""" - from qiskit_aer import Aer - - with warnings.catch_warnings(record=True) as caught_warnings: - warnings.simplefilter("always") - shor = Shor(quantum_instance=QuantumInstance(Aer.get_backend(backend), shots=1000)) - self.assertTrue("Shor class is deprecated" in str(caught_warnings[0].message)) - result = shor.factor(N=n_v) - self.assertListEqual(result.factors[0], factors) - self.assertTrue(result.total_counts >= result.successful_counts) - - @slow_test - @data(5, 7) - def test_shor_no_factors(self, n_v): - """shor no factors test""" - shor = self.instance - result = shor.factor(N=n_v) - self.assertTrue(result.factors == []) - self.assertTrue(result.successful_counts == 0) - - @idata( - [ - [3, 5], - [5, 3], - ] - ) - @unpack - def test_shor_input_being_power(self, base, power): - """shor input being power test""" - n_v = int(math.pow(base, power)) - shor = self.instance - result = shor.factor(N=n_v) - self.assertTrue(result.factors == [base]) - self.assertTrue(result.total_counts >= result.successful_counts) - - @idata( - [[N, 2] for N in [-1, 0, 1, 2, 4, 16]] + [[15, a] for a in [-1, 0, 1, 3, 5, 15, 16]], - ) - @unpack - def test_shor_bad_input(self, n_v, a_v): - """shor factor bad input test""" - shor = self.instance - with self.assertRaises(ValueError): - _ = shor.factor(N=n_v, a=a_v) - - @slow_test - @idata( - [ - [15, 4, 2], - [15, 7, 4], - ] - ) - @unpack - def test_shor_quantum_result(self, n_v, a_v, order): - """shor quantum result test (for order being power of 2)""" - self._test_quantum_result(a_v, n_v, order) - - @slow_test - @idata( - [ - [17, 8, 8], - [21, 13, 2], - ] - ) - @unpack - def test_shor_quantum_result_for_5_bit_number(self, n_v, a_v, order): - """shor quantum result test (for order being power of 2 and n = log(N) = 5)""" - self._test_quantum_result(a_v, n_v, order) - - def _test_quantum_result(self, a_v, n_v, order): - shor = self.instance - circuit = shor.construct_circuit(N=n_v, a=a_v, measurement=True) - - result = shor.quantum_instance.execute(circuit) - measurements = [int(key, base=2) for key in result.get_counts(circuit).keys()] - - # calculate values that could be measured - values = [i << (2 * n_v.bit_length() - order.bit_length() + 1) for i in range(order)] - - for measurement in measurements: - self.assertTrue(measurement in values) - - @slow_test - @idata( - [ - [15, 4, [1, 4]], - [15, 7, [1, 4, 7, 13]], - ] - ) - @unpack - def test_shor_exponentiation_result(self, n_v, a_v, values): - """shor exponentiation result test (for n = log(N) = 4)""" - self._test_exponentiation_result(a_v, n_v, values) - - @slow_test - @idata( - [ - [5, 21, [1, 4, 5, 16, 17, 20]], - [4, 25, [1, 4, 6, 9, 11, 14, 16, 19, 21, 24]], - ] - ) - @unpack - def test_shor_exponentiation_result_for_5_bit_number(self, a_v, n_v, values): - """shor exponentiation result test (for n = log(N) = 5)""" - self._test_exponentiation_result(a_v, n_v, values) - - def _test_exponentiation_result(self, a_v, n_v, values): - shor = self.instance - - circuit = shor.construct_circuit(N=n_v, a=a_v, measurement=False) - # modify circuit to measure output (down) register - down_qreg = circuit.qregs[1] - down_creg = ClassicalRegister(len(down_qreg), name="m") - circuit.add_register(down_creg) - circuit.measure(down_qreg, down_creg) - - result = shor.quantum_instance.execute(circuit) - measurements = [int(key, base=2) for key in result.get_counts(circuit).keys()] - - for measurement in measurements: - self.assertTrue(measurement in values) - - @idata([[2, 15, 8], [4, 15, 4]]) - @unpack - def test_shor_modinv(self, a_v, m_v, expected): - """shor modular inverse test""" - modinv = Shor.modinv(a_v, m_v) - self.assertTrue(modinv == expected) - - -if __name__ == "__main__": - unittest.main()