Source code for pysit.optimization.gauss_newton



import time
import copy

import numpy as np
from scipy.sparse.linalg import LinearOperator
from pyamg.krylov import cg, gmres

from pysit.optimization.optimization import OptimizationBase

__all__=['GaussNewton']

__docformat__ = "restructuredtext en"

[docs]class GaussNewton(OptimizationBase): def __init__(self, objective, krylov_maxiter=50, *args, **kwargs): OptimizationBase.__init__(self, objective, *args, **kwargs) self.krylov_maxiter = krylov_maxiter def _select_step(self, shots, current_objective_value, gradient, iteration, objective_arguments, **kwargs): """Compute the Gauss-Newton update for a set of shots. Gives the step s as a function of the gradient vector. Implemented as in p178 of Nocedal and Wright. Parameters ---------- shots : list of pysit.Shot List of Shots for which to compute. grad : ndarray Gradient vector. iteration : int Current time index. """ m0 = self.base_model rhs = -1*gradient.asarray() def matvec(x): m1 = m0.perturbation(data=x) return self.objective_function.apply_hessian(shots, self.base_model, x, **objective_arguments).data A_shape = (len(rhs), len(rhs)) A = LinearOperator(shape=A_shape, matvec=matvec, dtype=gradient.dtype) resid = [] # d, info = cg(A, rhs, maxiter=self.krylov_maxiter, residuals=resid) d, info = gmres(A, rhs, maxiter=self.krylov_maxiter, residuals=resid) d.shape = rhs.shape direction = m0.perturbation(data=d) if info < 0: print("CG Failure") if info == 0: print("CG Converge") if info > 0: print("CG ran {0} iterations".format(info)) alpha0_kwargs = {'reset' : False} if iteration == 0: alpha0_kwargs = {'reset' : True} alpha = self.select_alpha(shots, gradient, direction, objective_arguments, current_objective_value=current_objective_value, alpha0_kwargs=alpha0_kwargs, **kwargs) self._print(' alpha {0}'.format(alpha)) self.store_history('alpha', iteration, alpha) step = alpha * direction return step def _compute_alpha0(self, phi0, grad0, reset=False, upscale_factor=None, **kwargs): return 1.0