diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 8a9abaed..ff38d67b 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -588,7 +588,7 @@ class ProjectedGradient(Minimize, Remember): def reduceHess(v): # Z is tall and skinny return Z.T*(self.H*(Z*v)) - operator = sp.linalg.LinearOperator( (shape[1], shape[1]), reduceHess, dtype=float ) + operator = sp.linalg.LinearOperator( (shape[1], shape[1]), reduceHess, dtype=self.xc.dtype ) p, info = sp.linalg.cg(operator, -Z.T*self.g, tol=self.tolCG, maxiter=self.maxIterCG) p = Z*p # bring up to full size # aSet_after = self.activeSet(self.xc+p) @@ -622,21 +622,125 @@ class ProjectedGradient(Minimize, Remember): if self.debug: print 'doEndIteration.ProjGrad, f_decrease_max: ', self.f_decrease_max if self.debug: print 'doEndIteration.ProjGrad, stopDoingSD: ', self.stopDoingSD + + +class BFGS(Minimize, Remember): + name = 'BFGS' + nbfgs = 10 + + @property + def bfgsH0(self): + """ + Approximate Hessian used in preconditioning the problem. + + Must be a SimPEG.Solver + """ + _bfgsH0 = getattr(self,'_bfgsH0',None) + if _bfgsH0 is None: + return Solver(sp.identity(self.xc.size).tocsc(), flag='D') + return _bfgsH0 + @bfgsH0.setter + def bfgsH0(self, value): + assert type(value) is Solver, 'bfgsH0 must be a SimPEG.Solver' + self._bfgsH0 = value + + def _startup_BFGS(self,x0): + self._bfgscnt = -1 + self._bfgsY = np.zeros((x0.size, self.nbfgs)) + self._bfgsS = np.zeros((x0.size, self.nbfgs)) + if not np.any([p is IterationPrinters.comment for p in self.printers]): + self.printers.append(IterationPrinters.comment) + + def bfgs(self, d): + n = self._bfgscnt + nn = ktop = min(self._bfgsS.shape[1],n) + return self.bfgsrec(ktop,n,nn,self._bfgsS,self._bfgsY,d) + + def bfgsrec(self,k,n,nn,S,Y,d): + """BFGS recursion""" + if k < 0: + d = self.bfgsH0.solve(d) + else: + khat = np.mod(n-nn+k,nn) + gamma = np.vdot(S[:,khat],d)/np.vdot(Y[:,khat],S[:,khat]) + d = d - gamma*Y[:,khat] + d = self.bfgsrec(k-1,n,nn,S,Y,d) + d = d + (gamma - np.vdot(Y[:,khat],d)/np.vdot(Y[:,khat],S[:,khat]))*S[:,khat] + return d + + def findSearchDirection(self): + return self.bfgs(-self.g) + + def _doEndIteration_BFGS(self, xt): + if self._iter is 0: + self.g_last = self.g + return + + yy = self.g - self.g_last; + ss = self.xc - xt; + self.g_last = self.g + + if yy.dot(ss) > 0: + self._bfgscnt += 1 + ktop = np.mod(self._bfgscnt,self.nbfgs) + self._bfgsY[:,ktop] = yy + self._bfgsS[:,ktop] = ss + self.comment = '' + else: + self.comment = 'Skip BFGS' + + class GaussNewton(Minimize, Remember): name = 'Gauss Newton' def findSearchDirection(self): return Solver(self.H).solve(-self.g) -class InexactGaussNewton(Minimize, Remember): +class InexactGaussNewton(BFGS, Minimize, Remember): + """ + Minimizes using CG as the inexact solver of + + .. math:: + + \mathbf{H p = -g} + + By default BFGS is used as the preconditioner. + + Use *nbfgs* to set the memory limitation of BFGS. + + To set the initial H0 to be used in BFGS, set *bfgsH0* to be a SimPEG.Solver + + """ + + def __init__(self, **kwargs): + Minimize.__init__(self, **kwargs) + name = 'Inexact Gauss Newton' maxIterCG = 10 - tolCG = 1e-5 + tolCG = 1e-3 + + @property + def approxHinv(self): + """ + The approximate Hessian inverse is used to precondition CG. + + Default uses BFGS, with an initial H0 of *bfgsH0*. + + Must be a scipy.sparse.linalg.LinearOperator + """ + _approxHinv = getattr(self,'_approxHinv',None) + if _approxHinv is None: + M = sp.linalg.LinearOperator( (self.xc.size, self.xc.size), self.bfgs, dtype=self.xc.dtype ) + return M + return _approxHinv + @approxHinv.setter + def approxHinv(self, value): + self._approxHinv = value def findSearchDirection(self): - # TODO: use BFGS as a preconditioner or gauss sidel of the WtW or solve WtW directly - p, info = sp.linalg.cg(self.H, -self.g, tol=self.tolCG, maxiter=self.maxIterCG) + Hinv = Solver(self.H, doDirect=False, options={'iterSolver': 'CG', 'M': self.approxHinv, 'tol': self.tolCG, 'maxIter': self.maxIterCG}) + p = Hinv.solve(-self.g) return p diff --git a/SimPEG/utils/Solver.py b/SimPEG/utils/Solver.py index 82c4b709..16db872f 100644 --- a/SimPEG/utils/Solver.py +++ b/SimPEG/utils/Solver.py @@ -62,8 +62,11 @@ class Solver(object): M = options['M'] if type(M) is sp.linalg.LinearOperator: return - elif type(M) is tuple: - PreconditionerList = ['J','GS'] + PreconditionerList = ['J','GS'] + if type(M) is str: + assert M in PreconditionerList, "M must be in the known preconditioner list. ['J','GS']" + M = (M,A) # use A as the base for the preconditioner. + if type(M) is tuple: assert type(M[0]) is str and M[0] in PreconditionerList, "M as a tuple must be (str, Matrix) where str is in ['J','GS']: e.g. ('J', WtW) where J stands for Jacobi, and WtW is a sparse matrix." if M[0] is 'J': Jacobi = sdiag(1.0/M[1].diagonal()) diff --git a/notebooks/BFGS.ipynb b/notebooks/BFGS.ipynb index 23de6e06..c2f68c99 100644 --- a/notebooks/BFGS.ipynb +++ b/notebooks/BFGS.ipynb @@ -23,7 +23,7 @@ "language": "python", "metadata": {}, "outputs": [], - "prompt_number": 1 + "prompt_number": 2 }, { "cell_type": "code", @@ -32,77 +32,8 @@ "FUN = SimPEG.tests.Rosenbrock\n", "FUN = SimPEG.tests.getQuadratic(sp.csr_matrix(([100,1],([0,1],[0,1])),shape=(2,2)),np.array([-5,-5]),100)\n", "\n", - "\n", - "class BFGS(Minimize, Remember):\n", - " name = 'BFGS'\n", - " nbfgs = 10\n", - " \n", - " @property\n", - " def H0(self):\n", - " \"\"\"\n", - " Approximate Hessian used in preconditioning the problem.\n", - " \n", - " Must be a SimPEG.Solver\n", - " \"\"\"\n", - " _H0 = getattr(self,'_H0',None)\n", - " if _H0 is None:\n", - " return Solver(sp.identity(self.xc.size).tocsc())\n", - " return _H0\n", - " @H0.setter\n", - " def H0(self, value):\n", - " self._H0 = value\n", - " \n", - " def _startup_BFGS(self,x0):\n", - " self._bfgscnt = -1\n", - " self._bfgsY = np.zeros((x0.size, self.nbfgs))\n", - " self._bfgsS = np.zeros((x0.size, self.nbfgs))\n", - " if not np.any([p is IterationPrinters.comment for p in self.printers]):\n", - " self.printers.append(IterationPrinters.comment)\n", - " \n", - " def bfgs(self,n,d):\n", - " nn = min(self._bfgsS.shape[1],n)\n", - " ktop = nn\n", - " d = self.bfgsrec(ktop,n,nn,self._bfgsS,self._bfgsY,d)\n", - " return d\n", - "\n", - " def bfgsrec(self,k,n,nn,S,Y,d):\n", - " \"\"\"BFGS recursion\"\"\"\n", - " if k < 0:\n", - " d = self.H0.solve(d)\n", - " else:\n", - " khat = mod(n-nn+k,nn)\n", - " gamma = np.vdot(S[:,khat],d)/np.vdot(Y[:,khat],S[:,khat])\n", - " d = d - gamma*Y[:,khat]\n", - " d = self.bfgsrec(k-1,n,nn,S,Y,d)\n", - " d = d + (gamma - np.vdot(Y[:,khat],d)/np.vdot(Y[:,khat],S[:,khat]))*S[:,khat]\n", - " \n", - " return d\n", - " \n", - " def findSearchDirection(self):\n", - " return self.bfgs(self._bfgscnt,-self.g)\n", - " \n", - " def _doEndIteration_BFGS(self, xt):\n", - " if self._iter is 0: \n", - " self.g_last = self.g\n", - " return\n", - " \n", - " yy = self.g - self.g_last;\n", - " ss = self.xc - xt;\n", - " self.g_last = self.g\n", - " \n", - " if yy.dot(ss) > 0:\n", - " self._bfgscnt += 1\n", - " ktop = np.mod(self._bfgscnt,self.nbfgs)\n", - " self._bfgsY[:,ktop] = yy\n", - " self._bfgsS[:,ktop] = ss\n", - " self.comment = ''\n", - " else:\n", - " self.comment = 'Skip BFGS'\n", - " \n", - "\n", - "\n", "x0 = np.array([1,0])\n", - "opt = BFGS()\n", + "opt = inverse.BFGS()\n", "xopt = opt.minimize(FUN,x0)\n", "print xopt\n", "opt = inverse.GaussNewton()\n", @@ -186,9 +117,39 @@ "------------------------- DONE! -------------------------\n", "[ 0.07777107 1.6849632 ]\n" ] + }, + { + "output_type": "stream", + "stream": "stderr", + "text": [ + "/Users/rowan/git/simpeg/SimPEG/inverse/Optimize.py:664: RuntimeWarning: divide by zero encountered in remainder\n", + " khat = np.mod(n-nn+k,nn)\n" + ] } ], - "prompt_number": 14 + "prompt_number": 3 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "A = sp.identity(2)\n", + "S = Solver(A)\n", + "\n", + "assert type(S) is Solver" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 6 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [], + "language": "python", + "metadata": {}, + "outputs": [] } ], "metadata": {}