From fc435feb9e2ce03ff3c62611ce5f4fd05ca14c01 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 19 Nov 2013 15:59:37 -0800 Subject: [PATCH] BFGS implementation in the notebook. --- SimPEG/inverse/Optimize.py | 10 +- notebooks/BFGS.ipynb | 197 +++++++++++++++++++++++++++++++++++++ 2 files changed, 203 insertions(+), 4 deletions(-) create mode 100644 notebooks/BFGS.ipynb diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 37c8b296..8a9abaed 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -76,7 +76,7 @@ class IterationPrinters(object): itType = {"title": "itType", "value": lambda M: M._itType, "width": 8, "format": "%s"} aSet = {"title": "aSet", "value": lambda M: np.sum(M.activeSet(M.xc)), "width": 8, "format": "%d"} bSet = {"title": "bSet", "value": lambda M: np.sum(M.bindingSet(M.xc)), "width": 8, "format": "%d"} - comment = {"title": "Comment", "value": lambda M: M.projComment, "width": 7, "format": "%s"} + comment = {"title": "Comment", "value": lambda M: M.comment, "width": 12, "format": "%s"} beta = {"title": "beta", "value": lambda M: M.parent._beta, "width": 10, "format": "%1.2e"} phi_d = {"title": "phi_d", "value": lambda M: M.parent.phi_d, "width": 10, "format": "%1.2e"} @@ -106,6 +106,8 @@ class Minimize(object): debug = False debugLS = False + comment = '' + def __init__(self, **kwargs): self._id = int(np.random.rand()*1e6) # create a unique identifier to this program to be used in pubsub self.stoppers = [StoppingCriteria.tolerance_f, StoppingCriteria.moving_x, StoppingCriteria.tolerance_g, StoppingCriteria.norm_g, StoppingCriteria.iteration] @@ -525,7 +527,7 @@ class ProjectedGradient(Minimize, Remember): self.stopDoingPG = False self._itType = 'SD' - self.projComment = '' + self.comment = '' self.aSet_prev = self.activeSet(x0) @@ -600,7 +602,7 @@ class ProjectedGradient(Minimize, Remember): self.exploreCG = np.all(aSet == bSet) # explore conjugate gradient f_current_decrease = self.f_last - self.f - self.projComment = '' + self.comment = '' if self._iter < 1: # Note that this is reset on every CG iteration. self.f_decrease_max = -np.inf @@ -608,7 +610,7 @@ class ProjectedGradient(Minimize, Remember): self.f_decrease_max = max(self.f_decrease_max, f_current_decrease) self.stopDoingPG = f_current_decrease < 0.25 * self.f_decrease_max if self.stopDoingPG: - self.projComment = 'Stop SD' + self.comment = 'Stop SD' self.explorePG = False self.exploreCG = True # implement 3.8, MoreToraldo91 diff --git a/notebooks/BFGS.ipynb b/notebooks/BFGS.ipynb new file mode 100644 index 00000000..23de6e06 --- /dev/null +++ b/notebooks/BFGS.ipynb @@ -0,0 +1,197 @@ +{ + "metadata": { + "name": "" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "code", + "collapsed": false, + "input": [ + "import SimPEG\n", + "from SimPEG import Solver\n", + "from SimPEG.mesh import TensorMesh\n", + "from SimPEG.regularization import Regularization\n", + "import SimPEG.inverse as inverse\n", + "from SimPEG.inverse import Minimize, Remember, IterationPrinters\n", + "import numpy as np\n", + "import scipy.sparse as sp" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 1 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "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", + "xopt = opt.minimize(FUN,x0)\n", + "print xopt\n", + "opt = inverse.GaussNewton()\n", + "xopt = opt.minimize(FUN,x0)\n", + "print xopt\n", + "opt = inverse.SteepestDescent()\n", + "xopt = opt.minimize(FUN,x0)\n", + "print xopt" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "===================== BFGS =====================\n", + " # f |proj(x-g)-x| LS Comment \n", + "-----------------------------------------------\n", + " 0 1.45e+02 9.51e+01 0 \n", + " 1 1.14e+02 5.37e+01 6 \n", + " 2 1.04e+02 3.04e+01 6 \n", + " 3 8.83e+01 1.37e+01 0 \n", + " 4 8.76e+01 5.97e+00 0 Skip BFGS \n", + " 5 8.74e+01 2.61e+00 0 Skip BFGS \n", + " 6 8.74e+01 1.14e+00 0 Skip BFGS \n", + " 7 8.74e+01 5.01e-01 0 Skip BFGS \n", + " 8 8.74e+01 2.19e-01 0 Skip BFGS \n", + " 9 8.74e+01 9.60e-02 0 Skip BFGS \n", + "------------------------- STOP! -------------------------\n", + "1 : |fc-fOld| = 1.9437e-04 <= tolF*(1+|f0|) = 1.4600e+01\n", + "1 : |xc-x_last| = 1.2663e-03 <= tolX*(1+|x0|) = 2.0000e-01\n", + "1 : |proj(x-g)-x| = 9.5952e-02 <= tolG = 1.0000e-01\n", + "0 : |proj(x-g)-x| = 9.5952e-02 <= 1e3*eps = 1.0000e-02\n", + "0 : maxIter = 20 <= iter = 9\n", + "------------------------- DONE! -------------------------\n", + "[ 0.05095952 4.99977449]\n", + "=========== Gauss Newton ===========\n", + " # f |proj(x-g)-x| LS \n", + "-----------------------------------\n", + " 0 1.45e+02 9.51e+01 0 \n", + " 1 8.74e+01 4.44e-15 0 \n", + "------------------------- STOP! -------------------------\n", + "0 : |fc-fOld| = 5.7625e+01 <= tolF*(1+|f0|) = 1.4600e+01\n", + "0 : |xc-x_last| = 5.0894e+00 <= tolX*(1+|x0|) = 2.0000e-01\n", + "1 : |proj(x-g)-x| = 4.4409e-15 <= tolG = 1.0000e-01\n", + "1 : |proj(x-g)-x| = 4.4409e-15 <= 1e3*eps = 1.0000e-02\n", + "0 : maxIter = 20 <= iter = 1\n", + "------------------------- DONE! -------------------------\n", + "[ 0.05 5. ]\n", + "========= Steepest Descent =========\n", + " # f |proj(x-g)-x| LS \n", + "-----------------------------------\n", + " 0 1.45e+02 9.51e+01 0 \n", + " 1 1.14e+02 5.37e+01 6 \n", + " 2 1.04e+02 3.04e+01 6 \n", + " 3 1.00e+02 1.76e+01 6 \n", + " 4 9.88e+01 1.06e+01 6 \n", + " 5 9.82e+01 7.07e+00 6 \n", + " 6 9.80e+01 1.22e+01 5 \n", + " 7 9.73e+01 7.77e+00 6 \n", + " 8 9.68e+01 5.64e+00 6 \n", + " 9 9.65e+01 8.72e+00 5 \n", + " 10 9.60e+01 5.97e+00 6 \n", + " 11 9.58e+01 9.98e+00 5 \n", + " 12 9.53e+01 6.48e+00 6 \n", + " 13 9.53e+01 1.16e+01 5 \n", + " 14 9.46e+01 7.20e+00 6 \n", + " 15 9.43e+01 5.07e+00 6 \n", + " 16 9.41e+01 8.17e+00 5 \n", + " 17 9.37e+01 5.43e+00 6 \n", + " 18 9.36e+01 9.42e+00 5 \n", + " 19 9.32e+01 5.98e+00 6 \n", + " 20 9.29e+01 4.32e+00 6 \n", + "------------------------- STOP! -------------------------\n", + "1 : |fc-fOld| = 2.5913e-01 <= tolF*(1+|f0|) = 1.4600e+01\n", + "1 : |xc-x_last| = 9.3379e-02 <= tolX*(1+|x0|) = 2.0000e-01\n", + "0 : |proj(x-g)-x| = 4.3246e+00 <= tolG = 1.0000e-01\n", + "0 : |proj(x-g)-x| = 4.3246e+00 <= 1e3*eps = 1.0000e-02\n", + "1 : maxIter = 20 <= iter = 20\n", + "------------------------- DONE! -------------------------\n", + "[ 0.07777107 1.6849632 ]\n" + ] + } + ], + "prompt_number": 14 + } + ], + "metadata": {} + } + ] +} \ No newline at end of file