Source code for pylops.optimization.leastsquares

__all__ = [
    "normal_equations_inversion",
    "regularized_inversion",
    "preconditioned_inversion",
]

from typing import TYPE_CHECKING, List, Optional, Sequence, Tuple

from pylops.optimization.cls_leastsquares import (
    NormalEquationsInversion,
    PreconditionedInversion,
    RegularizedInversion,
)
from pylops.utils.typing import NDArray, SamplingLike

if TYPE_CHECKING:
    from pylops.linearoperator import LinearOperator


[docs]def normal_equations_inversion( Op: "LinearOperator", y: NDArray, Regs: List["LinearOperator"], x0: Optional[NDArray] = None, Weight: Optional["LinearOperator"] = None, dataregs: Optional[List[NDArray]] = None, epsI: float = 0.0, epsRs: Optional[SamplingLike] = None, NRegs: Optional[Sequence["LinearOperator"]] = None, epsNRs: Optional[SamplingLike] = None, engine: str = "scipy", show: bool = False, **kwargs_solver, ) -> Tuple[NDArray, int]: r"""Inversion of normal equations. Solve the regularized normal equations for a system of equations given the operator ``Op``, a data weighting operator ``Weight`` and optionally a list of regularization terms ``Regs`` and/or ``NRegs``. Parameters ---------- Op : :obj:`pylops.LinearOperator` Operator to invert of size :math:`[N \times M]` y : :obj:`numpy.ndarray` Data of size :math:`[N \times 1]` Regs : :obj:`list` Regularization operators (``None`` to avoid adding regularization) x0 : :obj:`numpy.ndarray`, optional Initial guess of size :math:`[M \times 1]` Weight : :obj:`pylops.LinearOperator`, optional Weight operator dataregs : :obj:`list`, optional Regularization data (must have the same number of elements as ``Regs``) epsI : :obj:`float`, optional Tikhonov damping epsRs : :obj:`list`, optional Regularization dampings (must have the same number of elements as ``Regs``) NRegs : :obj:`list` Normal regularization operators (``None`` to avoid adding regularization). Such operators must apply the chain of the forward and the adjoint in one go. This can be convenient in cases where a faster implementation is available compared to applying the forward followed by the adjoint. epsNRs : :obj:`list`, optional Regularization dampings for normal operators (must have the same number of elements as ``NRegs``) engine : :obj:`str`, optional Solver to use (``scipy`` or ``pylops``) show : :obj:`bool`, optional Display normal equations solver log **kwargs_solver Arbitrary keyword arguments for chosen solver (:py:func:`scipy.sparse.linalg.cg` and :py:func:`pylops.optimization.solver.cg` are used for engine ``scipy`` and ``pylops``, respectively) .. note:: When user does not supply ``atol``, it is set to "legacy". Returns ------- xinv : :obj:`numpy.ndarray` Inverted model. istop : :obj:`int` Convergence information (only when using :py:func:`scipy.sparse.linalg.cg`): ``0``: successful exit ``>0``: convergence to tolerance not achieved, number of iterations ``<0``: illegal input or breakdown See Also -------- RegularizedInversion: Regularized inversion PreconditionedInversion: Preconditioned inversion Notes ----- See :class:`pylops.optimization.cls_leastsquares.NormalEquationsInversion` """ nesolve = NormalEquationsInversion(Op) xinv, istop = nesolve.solve( y, Regs, x0=x0, Weight=Weight, dataregs=dataregs, epsI=epsI, epsRs=epsRs, NRegs=NRegs, epsNRs=epsNRs, engine=engine, show=show, **kwargs_solver, ) return xinv, istop
[docs]def regularized_inversion( Op, y: NDArray, Regs: List["LinearOperator"], x0: Optional[NDArray] = None, Weight: Optional["LinearOperator"] = None, dataregs: Optional[List[NDArray]] = None, epsRs: Optional[SamplingLike] = None, engine: str = "scipy", show: bool = False, **kwargs_solver, ) -> Tuple[NDArray, int, int, float, float]: r"""Regularized inversion. Solve a system of regularized equations given the operator ``Op``, a data weighting operator ``Weight``, and a list of regularization terms ``Regs``. Parameters ---------- Op : :obj:`pylops.LinearOperator` Operator to invert of size :math:`[N \times M]` y : :obj:`numpy.ndarray` Data of size :math:`[N \times 1]` Regs : :obj:`list` Regularization operators (``None`` to avoid adding regularization) x0 : :obj:`numpy.ndarray`, optional Initial guess of size :math:`[M \times 1]` Weight : :obj:`pylops.LinearOperator`, optional Weight operator dataregs : :obj:`list`, optional Regularization data (if ``None`` a zero data will be used for every regularization operator in ``Regs``) epsRs : :obj:`list`, optional Regularization dampings engine : :obj:`str`, optional Solver to use (``scipy`` or ``pylops``) show : :obj:`bool`, optional Display normal equations solver log **kwargs_solver Arbitrary keyword arguments for chosen solver (:py:func:`scipy.sparse.linalg.lsqr` and :py:func:`pylops.optimization.solver.cgls` are used for engine ``scipy`` and ``pylops``, respectively) Returns ------- xinv : :obj:`numpy.ndarray` Inverted model. istop : :obj:`int` Gives the reason for termination ``1`` means :math:`\mathbf{x}` is an approximate solution to :math:`\mathbf{y} = \mathbf{Op}\,\mathbf{x}` ``2`` means :math:`\mathbf{x}` approximately solves the least-squares problem itn : :obj:`int` Iteration number upon termination r1norm : :obj:`float` :math:`||\mathbf{r}||_2^2`, where :math:`\mathbf{r} = \mathbf{y} - \mathbf{Op}\,\mathbf{x}` r2norm : :obj:`float` :math:`\sqrt{\mathbf{r}^T\mathbf{r} + \epsilon^2 \mathbf{x}^T\mathbf{x}}`. Equal to ``r1norm`` if :math:`\epsilon=0` See Also -------- RegularizedOperator: Regularized operator NormalEquationsInversion: Normal equations inversion PreconditionedInversion: Preconditioned inversion Notes ----- See :class:`pylops.optimization.cls_leastsquares.RegularizedInversion` """ rsolve = RegularizedInversion(Op) xinv, istop, itn, r1norm, r2norm = rsolve.solve( y, Regs, x0=x0, Weight=Weight, dataregs=dataregs, epsRs=epsRs, engine=engine, show=show, **kwargs_solver, ) return xinv, istop, itn, r1norm, r2norm
[docs]def preconditioned_inversion( Op: "LinearOperator", y: NDArray, P: "LinearOperator", x0: Optional[NDArray] = None, engine: str = "scipy", show: bool = False, **kwargs_solver, ) -> Tuple[NDArray, int, int, float, float]: r"""Preconditioned inversion. Solve a system of preconditioned equations given the operator ``Op`` and a preconditioner ``P``. Parameters ---------- Op : :obj:`pylops.LinearOperator` Operator to invert of size :math:`[N \times M]` y : :obj:`numpy.ndarray` Data of size :math:`[N \times 1]` P : :obj:`pylops.LinearOperator` Preconditioner x0 : :obj:`numpy.ndarray` Initial guess of size :math:`[M \times 1]` engine : :obj:`str`, optional Solver to use (``scipy`` or ``pylops``) show : :obj:`bool`, optional Display normal equations solver log **kwargs_solver Arbitrary keyword arguments for chosen solver (:py:func:`scipy.sparse.linalg.lsqr` and :py:func:`pylops.optimization.solver.cgls` are used as default for numpy and cupy `data`, respectively) Returns ------- xinv : :obj:`numpy.ndarray` Inverted model. istop : :obj:`int` Gives the reason for termination ``1`` means :math:`\mathbf{x}` is an approximate solution to :math:`\mathbf{y} = \mathbf{Op}\,\mathbf{x}` ``2`` means :math:`\mathbf{x}` approximately solves the least-squares problem itn : :obj:`int` Iteration number upon termination r1norm : :obj:`float` :math:`||\mathbf{r}||_2^2`, where :math:`\mathbf{r} = \mathbf{y} - \mathbf{Op}\,\mathbf{x}` r2norm : :obj:`float` :math:`\sqrt{\mathbf{r}^T\mathbf{r} + \epsilon^2 \mathbf{x}^T\mathbf{x}}`. Equal to ``r1norm`` if :math:`\epsilon=0` See Also -------- RegularizedInversion: Regularized inversion NormalEquationsInversion: Normal equations inversion Notes ----- See :class:`pylops.optimization.cls_leastsquares.PreconditionedInversion` """ psolve = PreconditionedInversion(Op) xinv, istop, itn, r1norm, r2norm = psolve.solve( y, P, x0=x0, engine=engine, show=show, **kwargs_solver ) return xinv, istop, itn, r1norm, r2norm