Source code for pycif.plugins.minimizers.scipy_conjugate.minimize

import numpy as np
import scipy
from scipy.sparse.linalg import LinearOperator, lobpcg, cg
from logging import info




def minimize(self, finit, gradinit, chi0, **kwargs):
    """Solve the Gauss-Newton linear system via sparse conjugate gradient.

    Builds a :class:`~scipy.sparse.linalg.LinearOperator` whose
    matrix-vector product is evaluated by calling the TL/forward simulator,
    then dispatches to ``cg`` or ``lobpcg`` according to the ``method``
    parameter.

    Args:
        self (Plugin): scipy conjugate-gradient plugin instance.
        finit (float): initial cost value (unused; kept for interface
            consistency with other minimizers).
        gradinit (np.ndarray): initial gradient :math:`\\nabla J_0`,
            shape ``(n,)``; used as the right-hand side of the linear system.
        chi0 (np.ndarray): initial iterate (starting guess for the solver),
            shape ``(n,)``.
        **kwargs: forwarded to the simulator.

    Returns:
        np.ndarray: solution :math:`\chi_\mathrm{opt}`, shape ``(n,)``.
    """

    # Initialize iteration number and initial norm of the gradient
    self.iter_id = 0
    self.grad_norm_0 = np.sqrt(np.dot(gradinit, gradinit))

    # Define A.v
    def Amatvec(x):
        fb, fo, gb, go = self.simulator.simul(
            x.flatten(), run_id=self.iter_id,
            linear=self.force_linearize,
            split_obs_vs_control=True,
            **kwargs
        )
        self.iter_id += 1
        self.current_function = fo
        self.current_gradient = go[:]
        return go - gradinit

    # Defines the A matrix
    grad_dim = gradinit.shape[0]
    A = LinearOperator(
        (grad_dim, grad_dim),
        matvec=Amatvec
    )

    # Defines the preconditioner matrix

    # Running scipy minimizer
    try:
        results = cg(
            A, gradinit,
            maxiter=self.maxiter,
            x0=chi0[:, np.newaxis]
            # verbosityLevel=10
        )

        # Final verbose and output
        info(results)

        # Outputs
        xopt = results[0]
        # gradopt = results.jac
        # fopt = results.fun

    except StopIteration:
        info("The convergence criterion on epsg has been reached")
        xopt = self.xopt
        gradopt = self.current_gradient
        fopt = self.current_function

    r1 = np.sqrt(np.dot(xopt, xopt))
    # r2 = np.sqrt(np.dot(gradopt, gradopt))

    info("norm of x = " + str(r1))
    # info("f         = " + str(fopt))
    # info("norm of g = " + str(r2))

    return xopt