From 022e1f7660dd0b20e44b00ebb919a0a2cf77b036 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Fri, 27 May 2016 13:11:31 -0700 Subject: [PATCH] Update IRLS directive to allow multiple GN iterations. Remove modifications to the ProjGN solver. Update IRLS example. --- SimPEG/Directives.py | 156 ++++++++++++++++-------------- SimPEG/Examples/Inversion_IRLS.py | 14 +-- SimPEG/Optimization.py | 72 +------------- SimPEG/Regularization.py | 38 ++++++-- 4 files changed, 121 insertions(+), 159 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index f5c25249..b84d78d1 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -144,34 +144,6 @@ class BetaSchedule(InversionDirective): if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter self.invProb.beta /= self.coolingFactor -#class BetaSchedule_PGN_CG(InversionDirective): -# """BetaSchedule""" -# -# coolingFactor = 5. -# coolingRate = 1 -# GN_step_last = None -# GN_step_c = None -# -# def endIter(self): -# -# """ Compute the change in GN step, and proceed with cooling if below tol""" -# if self.opt.iter == 1: -# self.GN_step_last = np.linalg.norm(self.opt.xc - self.opt.x_last) -# d_GN_step = 1. -# -# else: -# self.GN_step_c = np.linalg.norm(self.opt.xc - self.opt.x_last) -# d_GN_step = self.GN_step_c / self.GN_step_last -# -# # Re-initiate last GN step -# self.GN_step_last = self.GN_step_c -# -# print "GN_step_last: ", self.GN_step_last -# print "d_GN_step: ", d_GN_step -# -# if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0: -# if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter -# self.invProb.beta /= self.coolingFactor class TargetMisfit(InversionDirective): @@ -286,10 +258,16 @@ class Update_IRLS(InversionDirective): gamma = None phi_m_last = None phi_d_last = None - - + f_old = None + f_min_change = 1e-1 + coolingRate = 3 + maxIRLSiter = 10 + + def initialize(self): - + + self.IRLSiter = 0 + # Scale the regularization for changes in norm if getattr(self, 'phi_m_last', None) is not None: @@ -310,48 +288,75 @@ class Update_IRLS(InversionDirective): self.reg._Wx = None self.reg._Wy = None self.reg._Wz = None - + if getattr(self, 'phi_d_last', None) is None: self.phi_d_last = self.invProb.phi_d + if getattr(self, 'f_last', None) is None: + self.f_old = self.invProb.evalFunction(self.reg.curModel, return_g=False, return_H=False) + print self.f_old def endIter(self): - # Cool the threshold parameter if required - if getattr(self, 'factor', None) is not None: - eps = self.reg.eps / self.factor + + + # Only update after GN iterations + if self.opt.iter % self.coolingRate == 0: + + self.IRLSiter += 1 - if getattr(self, 'eps_min', None) is not None: - self.reg.eps = np.max([self.eps_min,eps]) - else: - self.reg.eps = eps - - # Get phi_m at the end of current iteration - self.phi_m_last = self.invProb.phi_m_last - - # Reset the regularization matrices so that it is - # recalculated for current model - self.reg._Wsmall = None - self.reg._Wx = None - self.reg._Wy = None - self.reg._Wz = None - - # Update the model used for the IRLS weights - self.reg.curModel = self.invProb.curModel - - # Temporarely set gamma to 1. to get raw phi_m - self.reg.gamma = 1. - - # Compute new model objective function value - phim_new = self.reg.eval(self.invProb.curModel) - - # Update gamma to scale the regularization between IRLS iterations - self.reg.gamma = self.phi_m_last / phim_new - - # Reset the regularization matrices again for new gamma - self.reg._Wsmall = None - self.reg._Wx = None - self.reg._Wy = None - self.reg._Wz = None + self.f_change = np.abs(self.f_old - self.opt.f_last) / self.f_old + + print "Function decrease" + str(self.f_change) + + # Check for maximum number of IRLS cycles + if self.IRLSiter == self.maxIRLSiter: + self.opt.stopNextIteration = True + return + + # Check if the function has changed enough + if self.f_change < self.f_min_change: + self.opt.stopNextIteration = True + return + else: + self.f_old = self.opt.f_last + + # Cool the threshold parameter if required + if getattr(self, 'factor', None) is not None: + eps = self.reg.eps / self.factor + + if getattr(self, 'eps_min', None) is not None: + self.reg.eps = np.max([self.eps_min,eps]) + else: + self.reg.eps = eps + + # Get phi_m at the end of current iteration + self.phi_m_last = self.invProb.phi_m_last + + # Reset the regularization matrices so that it is + # recalculated for current model + self.reg._Wsmall = None + self.reg._Wx = None + self.reg._Wy = None + self.reg._Wz = None + + # Update the model used for the IRLS weights + self.reg.curModel = self.invProb.curModel + + # Temporarely set gamma to 1. to get raw phi_m + self.reg.gamma = 1. + + # Compute new model objective function value + phim_new = self.reg.eval(self.invProb.curModel) + + # Update gamma to scale the regularization between IRLS iterations + self.reg.gamma = self.phi_m_last / phim_new + + # Reset the regularization matrices again for new gamma + self.reg._Wsmall = None + self.reg._Wx = None + self.reg._Wy = None + self.reg._Wz = None + # Compute the change in class Update_lin_PreCond(InversionDirective): """ Create a Jacobi preconditioner for the linear problem @@ -411,11 +416,14 @@ class Scale_Beta(InversionDirective): update is done only if the misfit is outside some threshold bounds. """ tol = 0.05 - + coolingRate=5 + def endIter(self): - - # Check if misfit is within the tolerance, otherwise adjust beta - val = self.invProb.phi_d / (self.survey.nD*0.5) - - if np.abs(1.-val) > self.tol: - self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d + + # Only update after GN iterations + if self.opt.iter % self.coolingRate == 0: + # Check if misfit is within the tolerance, otherwise adjust beta + val = self.invProb.phi_d / (self.survey.nD*0.5) + + if np.abs(1.-val) > self.tol: + self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d diff --git a/SimPEG/Examples/Inversion_IRLS.py b/SimPEG/Examples/Inversion_IRLS.py index 6551bf21..d1da048f 100644 --- a/SimPEG/Examples/Inversion_IRLS.py +++ b/SimPEG/Examples/Inversion_IRLS.py @@ -20,7 +20,7 @@ def run(N=200, plotIt=True): m0 = np.ones(mesh.nC) * 1e-4 mref = np.zeros(mesh.nC) - nk = 10 + nk = 15 jk = np.linspace(1.,nk,nk) p = -2. q = 1. @@ -59,7 +59,7 @@ def run(N=200, plotIt=True): dmis = DataMisfit.l2_DataMisfit(survey) dmis.Wd = 1./wd - opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, maxIterGN=1, tolCG = 1e-4) + opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, tolCG = 1e-4) invProb = InvProblem.BaseInvProblem(dmis, reg, opt) invProb.curModel = m0 @@ -83,18 +83,18 @@ def run(N=200, plotIt=True): reg.cell_weights = wr reg.mref = np.zeros(mesh.nC) - reg.eps_p = 5e-2 + reg.eps_p = 1e-3 reg.eps_q = 1e-2 reg.norms = [0., 0., 2., 2.] - opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3) + opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3) invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.) beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1) - #betaest = Directives.BetaEstimate_ByEig() + update_beta = Directives.Scale_Beta(tol = 0.05, coolingRate=5) target = Directives.TargetMisfit() - IRLS =Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid ) + IRLS = Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid, coolingRate=5 ) - inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS]) + inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS,update_beta]) m0 = mrec diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py index 95f53320..05214dd5 100644 --- a/SimPEG/Optimization.py +++ b/SimPEG/Optimization.py @@ -893,10 +893,6 @@ class ProjectedGNCG(BFGS, Minimize, Remember): lower = -np.inf upper = np.inf - # Variables to control inner GN iterations - rdm_tol = 1e-2 # Tolerance for largest change in step (Default 1%) - maxIterGN = 3 # Maximum number of GN inner iterations - def _startup(self, x0): # ensure bound vectors are the same size as the model if type(self.lower) is not np.ndarray: @@ -942,72 +938,6 @@ class ProjectedGNCG(BFGS, Minimize, Remember): def approxHinv(self, value): self._approxHinv = value - @Utils.timeIt - def minimize(self, evalFunction, x0): - """minimize(evalFunction, x0) - - Minimizes the function (evalFunction) starting at the location x0. - - :param def evalFunction: function handle that evaluates: f, g, H = F(x) - :param numpy.ndarray x0: starting location - :rtype: numpy.ndarray - :return: x, the last iterate of the optimization algorithm - - The GN newton steps are repeated until it the maxIterGN is reached or the - relative step change falls below some tolrerance. - - """ - self.evalFunction = evalFunction - self.startup(x0) - self.printInit() - - - while True: - self.doStartIteration() - self.f, self.g, self.H = evalFunction(self.xc, return_g=True, return_H=True) - self.printIter() - if self.stoppingCriteria(): break - - # Inner GN iterations, stop on maximum number of iterations - # or on tolerance for step length change - GN_count = 0 - dm0 = None # Initial GN step length - dmc = None # Current GN step length - rdm = 1. # Relative change in step length - - while rdm > self.rdm_tol and GN_count < self.maxIterGN: - - GN_count += 1 - self.searchDirection = self.findSearchDirection() - p = self.scaleSearchDirection(self.searchDirection) - xt, passLS = self.modifySearchDirection(p) - if not passLS: - xt, caught = self.modifySearchDirectionBreak(p) - if not caught: return self.xc - - if GN_count == 1: - dm0 = np.linalg.norm(self.xc - xt) - dmc = dm0 - - else: - dmc = np.linalg.norm(self.xc - xt) - - rdm = dmc / dm0 - self.xc = xt # Update current model - - # Form system for next iteration - self.f, self.g, self.H = evalFunction(self.xc, return_g=True, return_H=True) - - print "GN iter: %i,\t dm: %8.5e,\t rdm: %8.5e"% (GN_count, dmc, rdm) - - self.doEndIteration(xt) - if self.stopNextIteration: break - - self.printDone() - self.finish() - - return self.xc - @Utils.timeIt def findSearchDirection(self): @@ -1078,4 +1008,4 @@ class ProjectedGNCG(BFGS, Minimize, Remember): indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0)) delx[indx] = 0. - return delx + return delx \ No newline at end of file diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 422c8f9c..c3b1c9e5 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -890,14 +890,37 @@ class Tikhonov(Simple): class Sparse(Simple): - + """ + The regularization is: + + .. math:: + + R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top R^\\top R W(m-m_\\text{ref})} + + where the IRLS weight + + .. math:: + + R = \eta TO FINISH LATER!!! + + So the derivative is straight forward: + + .. math:: + + R(m) = \mathbf{W^\\top R^\\top R W (m-m_\\text{ref})} + + The IRLS weights are recomputed after each beta solves. + It is strongly recommended to do a few Gauss-Newton iterations + before updating. + """ + # set default values - eps_p = 1e-1 - eps_q = 1e-1 - curModel = None # use a model to compute the weights - gamma = 1. - norms = [0., 2., 2., 2.] - cell_weights = 1. + eps_p = 1e-1 # Threshold value for the model norm + eps_q = 1e-1 # Threshold value for the model gradient norm + curModel = None # Requires model to compute the weights + gamma = 1. # Model norm scaling to smooth out convergence + norms = [0., 2., 2., 2.] # Values for norm on (m, dmdx, dmdy, dmdz) + cell_weights = 1. # Consider overwriting with sensitivity weights def __init__(self, mesh, mapping=None, indActive=None, **kwargs): Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs) @@ -971,6 +994,7 @@ class Sparse(Simple): def R(self, f_m , eps, exponent): + # Eta scaling is important for mix-norms...do not mess with it eta = (eps**(1.-exponent/2.))**0.5 r = eta / (f_m**2.+ eps**2.)**((1.-exponent/2.)/2.)