wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-11 08:17:48 +08:00
Code Issues Packages Projects Releases Wiki Activity
Files
fc435feb9e2ce03ff3c62611ce5f4fd05ca14c01
simpeg/notebooks/BFGS.ipynb
T
Rowan Cockett fc435feb9e BFGS implementation in the notebook.
2013-11-19 15:59:37 -08:00

197 lines
7.9 KiB
Plaintext
Raw Blame History

{
"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": {}
}
]
}
Reference in New Issue View Git Blame Copy Permalink