mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-28 11:26:12 +08:00
Propose change to the Projected_GNCG solver. Add inner GN iterations. Nice improvement to the convergence of IRLS
This commit is contained in:
D Fournier
+40
-9
@@ -144,6 +144,35 @@ 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):
|
||||
|
||||
@property
|
||||
@@ -265,13 +294,16 @@ class Update_IRLS(InversionDirective):
|
||||
if getattr(self, 'phi_m_last', None) is not None:
|
||||
|
||||
self.reg.curModel = self.invProb.curModel
|
||||
self.reg._Wsmall = None
|
||||
self.reg._Wx = None
|
||||
self.reg._Wy = None
|
||||
self.reg._Wz = None
|
||||
self.reg.gamma = 1.
|
||||
phim_new = self.reg.eval(self.invProb.curModel)
|
||||
self.gamma = self.phi_m_last / phim_new
|
||||
|
||||
self.reg.curModel = self.invProb.curModel
|
||||
self.reg.gamma = self.gamma
|
||||
print "Initial gamma ", np.linalg.norm(self.reg.gamma)
|
||||
# Reset the regularization matrices so that it is
|
||||
# recalculated with new gamma
|
||||
self.reg._Wsmall = None
|
||||
@@ -295,12 +327,6 @@ class Update_IRLS(InversionDirective):
|
||||
# Get phi_m at the end of current iteration
|
||||
self.phi_m_last = self.invProb.phi_m_last
|
||||
|
||||
# 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.
|
||||
|
||||
# Reset the regularization matrices so that it is
|
||||
# recalculated for current model
|
||||
self.reg._Wsmall = None
|
||||
@@ -308,13 +334,18 @@ class Update_IRLS(InversionDirective):
|
||||
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
|
||||
print "New gamma ", np.linalg.norm(self.reg.gamma)
|
||||
|
||||
|
||||
# Reset the regularization matrices again for new gamma
|
||||
self.reg._Wsmall = None
|
||||
self.reg._Wx = None
|
||||
|
||||
@@ -18,6 +18,8 @@ def run(N=200, plotIt=True):
|
||||
mesh = Mesh.TensorMesh([N])
|
||||
|
||||
m0 = np.ones(mesh.nC) * 1e-4
|
||||
mref = np.zeros(mesh.nC)
|
||||
|
||||
nk = 10
|
||||
jk = np.linspace(1.,nk,nk)
|
||||
p = -2.
|
||||
@@ -51,12 +53,13 @@ def run(N=200, plotIt=True):
|
||||
wr = ( wr/np.max(wr) )
|
||||
|
||||
reg = Regularization.Simple(mesh)
|
||||
reg.wght = wr
|
||||
reg.mref = mref
|
||||
reg.cell_weights = wr
|
||||
|
||||
dmis = DataMisfit.l2_DataMisfit(survey)
|
||||
dmis.Wd = 1./wd
|
||||
|
||||
opt = Optimization.ProjectedGNCG(maxIter=30,lower=-2.,upper=2., maxIterCG= 20, tolCG = 1e-4)
|
||||
opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, maxIterGN=1, tolCG = 1e-4)
|
||||
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
|
||||
invProb.curModel = m0
|
||||
|
||||
@@ -76,22 +79,15 @@ def run(N=200, plotIt=True):
|
||||
phid = invProb.phi_d
|
||||
|
||||
reg = Regularization.Sparse(mesh)
|
||||
reg.mref = mref
|
||||
reg.cell_weights = wr
|
||||
|
||||
#==============================================================================
|
||||
# fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
|
||||
# dmdx = reg.mesh.cellDiffxStencil * mrec
|
||||
# plt.plot(np.sort(dmdx))
|
||||
#==============================================================================
|
||||
|
||||
#reg.recModel = mrec
|
||||
# reg.cell_weight = np.ones(mesh.nC)
|
||||
reg.mref = np.zeros(mesh.nC)
|
||||
reg.eps_p = 5e-2
|
||||
reg.eps_q = 1e-2
|
||||
reg.norms = [0., 0., 2., 2.]
|
||||
reg.cell_weight = wr
|
||||
|
||||
opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 20, tolCG = 1e-3)
|
||||
opt = Optimization.ProjectedGNCG(maxIter=10 ,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()
|
||||
|
||||
@@ -893,6 +893,10 @@ 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:
|
||||
@@ -938,6 +942,72 @@ 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):
|
||||
|
||||
|
||||
Reference in New Issue
Block a user