French
Languages
English
Bengali
French
German
Japanese
Korean
Portuguese
Spanish
Tamil

```# This code is part of Qiskit.
#
# (C) Copyright IBM 2018, 2023.
#
# obtain a copy of this license in the LICENSE.txt file in the root directory
#
# 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.

from collections.abc import Iterable
from typing import List, Tuple, Callable, Optional, Union
import functools
import numpy as np

from qiskit.circuit.quantumcircuit import _compare_parameters
from qiskit.circuit import ParameterVector, ParameterExpression
from qiskit.utils import optionals as _optionals
from qiskit.utils.deprecation import deprecate_func
from ..operator_base import OperatorBase
from ..list_ops.list_op import ListOp
from ..list_ops.composed_op import ComposedOp
from ..state_fns.circuit_state_fn import CircuitStateFn
from .circuit_qfis import CircuitQFI
from .qfi import QFI

# Error tolerance variable
ETOL = 1e-8
# Cut-off ratio for small singular values for least square solver
RCOND = 1e-2

r"""Deprecated: Convert an operator expression to the first-order gradient.

Given an ill-posed inverse problem

x = arg min{||Ax-C||^2} (1)

one can use regularization schemes can be used to stabilize the system and find a numerical
solution

x_lambda = arg min{||Ax-C||^2 + lambda*R(x)} (2)

where R(x) represents the penalization term.
"""

@deprecate_func(
since="0.24.0",
additional_msg="For code migration guidelines, visit https://qisk.it/opflow_migration.",
)
def __init__(
self,
qfi_method: Union[str, CircuitQFI] = "lin_comb_full",
regularization: Optional[str] = None,
**kwargs,
):
r"""
Args:
grad_method: The method used to compute the state gradient. Can be either
``'param_shift'`` or ``'lin_comb'`` or ``'fin_diff'``.
qfi_method: The method used to compute the QFI. Can be either
``'lin_comb_full'`` or ``'overlap_block_diag'`` or ``'overlap_diag'``.
regularization: Use the following regularization with a least square method to solve the
underlying system of linear equations
Can be either None or ``'ridge'`` or ``'lasso'`` or ``'perturb_diag'``
``'ridge'`` and ``'lasso'`` use an automatic optimal parameter search
If regularization is None but the metric is ill-conditioned or singular then
a least square solver is used without regularization
kwargs (dict): Optional parameters for a CircuitGradient
"""

self._qfi_method = QFI(qfi_method)
self._regularization = regularization
self._epsilon = kwargs.get("epsilon", 1e-6)

[docs]    def convert(
self,
operator: OperatorBase,
params: Optional[
Union[ParameterVector, ParameterExpression, List[ParameterExpression]]
] = None,
) -> OperatorBase:
r"""
Args:
operator: The operator we are taking the gradient of.
params: The parameters we are taking the gradient with respect to. If not explicitly
passed, they are inferred from the operator and sorted by name.

Returns:
An operator whose evaluation yields the NaturalGradient.

Raises:
TypeError: If ``operator`` does not represent an expectation value or the quantum
state is not ``CircuitStateFn``.
ValueError: If ``params`` contains a parameter not present in ``operator``.
ValueError: If ``operator`` is not parameterized.
"""
if not isinstance(operator, ComposedOp):
if not (isinstance(operator, ListOp) and len(operator.oplist) == 1):
raise TypeError(
"Please provide the operator either as ComposedOp or as ListOp of "
"a CircuitStateFn potentially with a combo function."
)

if not isinstance(operator[-1], CircuitStateFn):
raise TypeError(
"Please make sure that the operator for which you want to compute "
"Quantum Fisher Information represents an expectation value or a "
"loss function and that the quantum state is given as "
"CircuitStateFn."
)
if len(operator.parameters) == 0:
raise ValueError("The operator we are taking the gradient of is not parameterized!")
if params is None:
params = sorted(operator.parameters, key=functools.cmp_to_key(_compare_parameters))
if not isinstance(params, Iterable):
params = [params]
# Instantiate the QFI metric which is used to re-scale the gradient
metric = self._qfi_method.convert(operator[-1], params) * 0.25

def combo_fn(x):

# Define the ListOp which combines the gradient and the QFI according to the combination
# function defined above.

[docs]    @staticmethod
def nat_grad_combo_fn(x: tuple, regularization: Optional[str] = None) -> np.ndarray:
r"""

Args:
x: Iterable consisting of Gradient, Quantum Fisher Information.
regularization: Regularization method.

Returns:

Raises:
ValueError: If the gradient has imaginary components that are non-negligible.

"""
metric = x
raise ValueError(
"The imaginary part of the gradient are non-negligible. The largest absolute "
"Please increase the number of shots."
)

if np.amax(np.abs(np.imag(metric))) > ETOL:
raise ValueError(
"The imaginary part of the metric are non-negligible. The largest "
"absolute imaginary value in the gradient is "
"increase the number of shots."
)
metric = np.real(metric)

if regularization is not None:
# If a regularization method is chosen then use a regularized solver to
)
else:
# Check if numerical instabilities lead to a metric which is not positive semidefinite
w, v = np.linalg.eigh(metric)

if not all(ew >= (-1) * ETOL for ew in w):
raise ValueError(
f"The underlying metric has at least one Eigenvalue < -{ETOL}. "
f"The smallest Eigenvalue is {np.amin(w)} "
"Please use a regularized least-square solver for this problem or "
"increase the number of backend shots.",
)
if not all(ew >= 0 for ew in w):
# If not all eigenvalues are non-negative, set them to a small positive
# value
w = [max(ETOL, ew) for ew in w]
# Recompose the adapted eigenvalues with the eigenvectors to get a new metric
metric = np.real(v @ np.diag(w) @ np.linalg.inv(v))

@property
def qfi_method(self) -> CircuitQFI:
"""Returns ``CircuitQFI``.

Returns: ``CircuitQFI``.

"""
return self._qfi_method.qfi_method

@property
def regularization(self) -> Optional[str]:
"""Returns the regularization option.

Returns: the regularization option.

"""
return self._regularization

@staticmethod
def _reg_term_search(
metric: np.ndarray,
reg_method: Callable[[np.ndarray, np.ndarray, float], float],
lambda1: float = 1e-3,
lambda4: float = 1.0,
tol: float = 1e-8,
) -> Tuple[float, np.ndarray]:
"""
This method implements a search for a regularization parameter lambda by finding for the
corner of the L-curve.
More explicitly, one has to evaluate a suitable lambda by finding a compromise between
the error in the solution and the norm of the regularization.
This function implements a method presented in
`A simple algorithm to find the L-curve corner in the regularization of inverse problems
<https://arxiv.org/pdf/1608.04571.pdf>`

Args:
metric: See (1) and (2).
reg_method: Given the metric, gradient and lambda the regularization method must return
``x_lambda`` - see (2).
lambda1: Left starting point for L-curve corner search.
lambda4: Right starting point for L-curve corner search.
tol: Termination threshold.

Returns:
Regularization coefficient which is the solution to the regularization inverse problem.
"""

def _get_curvature(x_lambda: List) -> float:
"""Calculate Menger curvature

Menger, K. (1930).  Untersuchungen  ̈uber Allgemeine Metrik. Math. Ann.,103(1), 466–501

Args:
``x_lambda: [[x_lambdaj], [x_lambdak], [x_lambdal]]``
``lambdaj < lambdak < lambdal``

Returns:
Menger Curvature

"""
eps = []
eta = []
for x in x_lambda:
try:
eps.append(np.log(np.linalg.norm(np.matmul(metric, x) - gradient) ** 2))
except ValueError:
eps.append(
np.log(np.linalg.norm(np.matmul(metric, np.transpose(x)) - gradient) ** 2)
)
eta.append(np.log(max(np.linalg.norm(x) ** 2, ETOL)))
p_temp = 1
c_k = 0
for i in range(3):
p_temp *= (eps[np.mod(i + 1, 3)] - eps[i]) ** 2 + (
eta[np.mod(i + 1, 3)] - eta[i]
) ** 2
c_k += eps[i] * eta[np.mod(i + 1, 3)] - eps[np.mod(i + 1, 3)] * eta[i]
c_k = 2 * c_k / max(1e-4, np.sqrt(p_temp))
return c_k

def get_lambda2_lambda3(lambda1, lambda4):
gold_sec = (1 + np.sqrt(5)) / 2.0
lambda2 = 10 ** ((np.log10(lambda4) + np.log10(lambda1) * gold_sec) / (1 + gold_sec))
lambda3 = 10 ** (np.log10(lambda1) + np.log10(lambda4) - np.log10(lambda2))
return lambda2, lambda3

lambda2, lambda3 = get_lambda2_lambda3(lambda1, lambda4)
lambda_ = [lambda1, lambda2, lambda3, lambda4]
x_lambda = []
for lam in lambda_:
counter = 0
while (lambda_ - lambda_) / lambda_ >= tol:
counter += 1
c_2 = _get_curvature(x_lambda[:-1])
c_3 = _get_curvature(x_lambda[1:])
while c_3 < 0:
lambda_ = lambda_
x_lambda = x_lambda
lambda_ = lambda_
x_lambda = x_lambda
lambda2, _ = get_lambda2_lambda3(lambda_, lambda_)
lambda_ = lambda2
c_3 = _get_curvature(x_lambda[1:])

if c_2 > c_3:
lambda_mc = lambda_
x_mc = x_lambda
lambda_ = lambda_
x_lambda = x_lambda
lambda_ = lambda_
x_lambda = x_lambda
lambda2, _ = get_lambda2_lambda3(lambda_, lambda_)
lambda_ = lambda2
else:
lambda_mc = lambda_
x_mc = x_lambda
lambda_ = lambda_
x_lambda = x_lambda
lambda_ = lambda_
x_lambda = x_lambda
_, lambda3 = get_lambda2_lambda3(lambda_, lambda_)
lambda_ = lambda3
return lambda_mc, x_mc

@staticmethod
@_optionals.HAS_SKLEARN.require_in_call
def _ridge(
metric: np.ndarray,
lambda_: float = 1.0,
lambda1: float = 1e-4,
lambda4: float = 1e-1,
tol_search: float = 1e-8,
fit_intercept: bool = True,
normalize: bool = False,
copy_a: bool = True,
max_iter: int = 1000,
tol: float = 0.0001,
solver: str = "auto",
random_state: Optional[int] = None,
) -> Tuple[float, np.ndarray]:
"""
Ridge Regression with automatic search for a good regularization term lambda
x_lambda = arg min{||Ax-C||^2 + lambda*||x||_2^2} (3)
`Scikit Learn Ridge Regression
<https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Ridge.html>`

Args:
metric: See (1) and (2).
lambda_ : regularization parameter used if auto_search = False
lambda1: left starting point for L-curve corner search
lambda4: right starting point for L-curve corner search
tol_search: termination threshold for regularization parameter search
fit_intercept: if True calculate intercept
normalize: ignored if fit_intercept=False, if True normalize A for regression
copy_a: if True A is copied, else overwritten
max_iter: max. number of iterations if solver is CG
tol: precision of the regression solution
solver: solver {‘auto’, ‘svd’, ‘cholesky’, ‘lsqr’, ‘sparse_cg’, ‘sag’, ‘saga’}
random_state: seed for the pseudo random number generator used when data is shuffled

Returns:
regularization coefficient, solution to the regularization inverse problem

Raises:
MissingOptionalLibraryError: scikit-learn not installed

"""
from sklearn.linear_model import Ridge
from sklearn.preprocessing import StandardScaler

reg = Ridge(
alpha=lambda_,
fit_intercept=fit_intercept,
copy_X=copy_a,
max_iter=max_iter,
tol=tol,
solver=solver,
random_state=random_state,
)

def reg_method(a, c, alpha):
reg.set_params(alpha=alpha)
if normalize:
reg.fit(StandardScaler().fit_transform(a), c)
else:
reg.fit(a, c)
return reg.coef_

metric, gradient, reg_method, lambda1=lambda1, lambda4=lambda4, tol=tol_search
)
return lambda_mc, np.transpose(x_mc)

@staticmethod
@_optionals.HAS_SKLEARN.require_in_call
def _lasso(
metric: np.ndarray,
lambda_: float = 1.0,
lambda1: float = 1e-4,
lambda4: float = 1e-1,
tol_search: float = 1e-8,
fit_intercept: bool = True,
normalize: bool = False,
precompute: Union[bool, Iterable] = False,
copy_a: bool = True,
max_iter: int = 1000,
tol: float = 0.0001,
warm_start: bool = False,
positive: bool = False,
random_state: Optional[int] = None,
selection: str = "random",
) -> Tuple[float, np.ndarray]:
"""
Lasso Regression with automatic search for a good regularization term lambda
x_lambda = arg min{||Ax-C||^2/(2*n_samples) + lambda*||x||_1} (4)
`Scikit Learn Lasso Regression
<https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Lasso.html>`

Args:
metric: Matrix of size mxn.
lambda_ : regularization parameter used if auto_search = False
lambda1: left starting point for L-curve corner search
lambda4: right starting point for L-curve corner search
tol_search: termination threshold for regularization parameter search
fit_intercept: if True calculate intercept
normalize: ignored if fit_intercept=False, if True normalize A for regression
precompute: If True compute and use Gram matrix to speed up calculations.
Gram matrix can also be given explicitly
copy_a: if True A is copied, else overwritten
max_iter: max. number of iterations if solver is CG
tol: precision of the regression solution
warm_start: if True reuse solution from previous fit as initialization
positive: if True force positive coefficients
random_state: seed for the pseudo random number generator used when data is shuffled
selection: {'cyclic', 'random'}

Returns:
regularization coefficient, solution to the regularization inverse problem

Raises:
MissingOptionalLibraryError: scikit-learn not installed

"""
from sklearn.linear_model import Lasso
from sklearn.preprocessing import StandardScaler

reg = Lasso(
alpha=lambda_,
fit_intercept=fit_intercept,
precompute=precompute,
copy_X=copy_a,
max_iter=max_iter,
tol=tol,
warm_start=warm_start,
positive=positive,
random_state=random_state,
selection=selection,
)

def reg_method(a, c, alpha):
reg.set_params(alpha=alpha)
if normalize:
reg.fit(StandardScaler().fit_transform(a), c)
else:
reg.fit(a, c)
return reg.coef_

metric, gradient, reg_method, lambda1=lambda1, lambda4=lambda4, tol=tol_search
)

return lambda_mc, x_mc

@staticmethod
def _regularized_sle_solver(
metric: np.ndarray,
regularization: str = "perturb_diag",
lambda1: float = 1e-3,
lambda4: float = 1.0,
alpha: float = 0.0,
tol_norm_x: Tuple[float, float] = (1e-8, 5.0),
tol_cond_a: float = 1000.0,
) -> np.ndarray:
"""
Solve a linear system of equations with a regularization method and automatic lambda fitting

Args:
metric: Matrix of size mxn.
regularization: Regularization scheme to be used: 'ridge', 'lasso',
'perturb_diag_elements' or 'perturb_diag'
lambda1: left starting point for L-curve corner search (for 'ridge' and 'lasso')
lambda4: right starting point for L-curve corner search (for 'ridge' and 'lasso')
alpha: perturbation coefficient for 'perturb_diag_elements' and 'perturb_diag'
tol_norm_x: tolerance for the norm of x
tol_cond_a: tolerance for the condition number of A

Returns:
solution to the regularized system of linear equations

"""
if regularization == "ridge":
elif regularization == "lasso":
elif regularization == "perturb_diag_elements":
alpha = 1e-7
while np.linalg.cond(metric + alpha * np.diag(metric)) > tol_cond_a:
alpha *= 10
# include perturbation in A to avoid singularity
x, _, _, _ = np.linalg.lstsq(metric + alpha * np.diag(metric), gradient, rcond=None)
elif regularization == "perturb_diag":
alpha = 1e-7
while np.linalg.cond(metric + alpha * np.eye(len(gradient))) > tol_cond_a:
alpha *= 10
# include perturbation in A to avoid singularity
x, _, _, _ = np.linalg.lstsq(
)
else:
# include perturbation in A to avoid singularity
x, _, _, _ = np.linalg.lstsq(metric, gradient, rcond=None)

if np.linalg.norm(x) > tol_norm_x or np.linalg.norm(x) < tol_norm_x:
if regularization == "ridge":
lambda1 = lambda1 / 10.0
elif regularization == "lasso":
lambda1 = lambda1 / 10.0
elif regularization == "perturb_diag_elements":
while np.linalg.cond(metric + alpha * np.diag(metric)) > tol_cond_a:
if alpha == 0:
alpha = 1e-7
else:
alpha *= 10
# include perturbation in A to avoid singularity
x, _, _, _ = np.linalg.lstsq(metric + alpha * np.diag(metric), gradient, rcond=None)
else:
if alpha == 0:
alpha = 1e-7
else:
alpha *= 10
while np.linalg.cond(metric + alpha * np.eye(len(gradient))) > tol_cond_a:
# include perturbation in A to avoid singularity
x, _, _, _ = np.linalg.lstsq(