mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-05 01:08:53 +08:00
Update IRLS directive to allow multiple GN iterations.
Remove modifications to the ProjGN solver. Update IRLS example.
This commit is contained in:
D Fournier
+82
-74
@@ -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
|
||||
|
||||
@@ -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
|
||||
|
||||
|
||||
+1
-71
@@ -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
|
||||
@@ -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.)
|
||||
|
||||
|
||||
Reference in New Issue
Block a user