From fc435feb9e2ce03ff3c62611ce5f4fd05ca14c01 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 19 Nov 2013 15:59:37 -0800 Subject: [PATCH 1/3] 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 From 4472b39d6e66ced02c7985b164d07ce8874f19e0 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 19 Nov 2013 16:00:36 -0800 Subject: [PATCH 2/3] Iterative Solvers in SimPEG.Solver --- SimPEG/tests/HTMLTestRunner.py | 3 +- SimPEG/tests/TestUtils.py | 8 ++-- SimPEG/utils/Solver.py | 78 ++++++++++++++++++++++++++++------ 3 files changed, 72 insertions(+), 17 deletions(-) diff --git a/SimPEG/tests/HTMLTestRunner.py b/SimPEG/tests/HTMLTestRunner.py index af384971..05ae09df 100644 --- a/SimPEG/tests/HTMLTestRunner.py +++ b/SimPEG/tests/HTMLTestRunner.py @@ -241,7 +241,8 @@ function showClassDetail(cid, count) { for (var i = 0; i < count; i++) { tid = id_list[i]; if (toHide) { - document.getElementById('div_'+tid).style.display = 'none' + var divTid = document.getElementById('div_'+tid); + if(divTid !== null){divTid.style.display = 'none';} document.getElementById(tid).className = 'hiddenRow'; } else { diff --git a/SimPEG/tests/TestUtils.py b/SimPEG/tests/TestUtils.py index 1cc2bbae..bfb150c3 100644 --- a/SimPEG/tests/TestUtils.py +++ b/SimPEG/tests/TestUtils.py @@ -276,16 +276,16 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None): -def getQuadratic(A, b): +def getQuadratic(A, b, c=0): """ - Given A and b, this returns a quadratic, Q + Given A, b and c, this returns a quadratic, Q .. math:: - \mathbf{Q( x ) = 0.5 x A x + b x} + \mathbf{Q( x ) = 0.5 x A x + b x} + c """ def Quadratic(x, return_g=True, return_H=True): - f = 0.5 * x.dot( A.dot(x)) + b.dot( x ) + f = 0.5 * x.dot( A.dot(x)) + b.dot( x ) + c out = (f,) if return_g: g = A.dot(x) + b diff --git a/SimPEG/utils/Solver.py b/SimPEG/utils/Solver.py index a1194487..82c4b709 100644 --- a/SimPEG/utils/Solver.py +++ b/SimPEG/utils/Solver.py @@ -1,9 +1,10 @@ import numpy as np -import scipy.sparse as sparse +import scipy.sparse as sp import scipy.sparse.linalg as linalg -from SimPEG.utils import mkvc +from SimPEG.utils import mkvc, sdiag +import warnings -DEFAULTS = {'direct':'scipy', 'forward':'fortran', 'backward':'fortran', 'diagonal':'python'} +DEFAULTS = {'direct':'scipy', 'iter':'scipy', 'forward':'fortran', 'backward':'fortran', 'diagonal':'python'} try: import TriSolve @@ -45,13 +46,41 @@ class Solver(object): def __init__(self, A, doDirect=True, flag=None, options={}): assert type(doDirect) is bool, 'doDirect must be a boolean' assert flag in [None, 'L', 'U', 'D'], "flag must be set to None, 'L', 'U', or 'D'" - + assert type(options) is dict, 'options must be a dictionary object' self.A = A self.dsolve = None self.doDirect = doDirect self.flag = flag self.options = options + if doDirect: return + + # Now deal with iterative stuff only + if 'M' not in options: + warnings.warn("You should provide a preconditioner, M.", UserWarning) + return + M = options['M'] + if type(M) is sp.linalg.LinearOperator: + return + elif type(M) is tuple: + PreconditionerList = ['J','GS'] + 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()) + options['M'] = Jacobi + elif M[0] is 'GS': + LL = sp.tril(M[1]) + UU = sp.triu(M[1]) + DD = sdiag(M[1].diagonal()) + Uinv = Solver(UU, flag='U') + Linv = Solver(LL, flag='L') + def GS(f): + return Uinv.solve(DD*Linv.solve(f)) + options['M'] = sp.linalg.LinearOperator( A.shape, GS, dtype=A.dtype ) + + else: + raise Exception('M must be a LinearOperator or a tuple') + def solve(self, b): """ @@ -118,8 +147,20 @@ class Solver(object): return X - def solveIter(self, b, M=None, iterSolver='CG'): - pass + def solveIter(self, b, backend=None, M=None, iterSolver='CG', tol=1e-6, maxIter=50): + if backend is None: backend = DEFAULTS['iter'] + + algorithms = {'CG':sp.linalg.cg} + assert iterSolver in algorithms, "iterSolver must be 'CG', or implement it yourself and add it here!" + alg = algorithms[iterSolver] + + if len(b.shape) == 1 or b.shape[1] == 1: + x, self.info = alg(self.A, b, M=M, tol=tol, maxiter=maxIter) + else: + x = np.empty_like(b) + for i in range(b.shape[1]): + x[:,i], self.info = alg(self.A, b[:,i], M=M, tol=tol, maxiter=maxIter) + return x def solveBackward(self, b, backend=None): """ @@ -132,9 +173,8 @@ class Solver(object): :return: x """ if backend is None: backend = DEFAULTS['backward'] - if type(self.A) is not sparse.csr.csr_matrix: - from scipy.sparse import csr_matrix - self.A = csr_matrix(self.A) + if type(self.A) is not sp.csr.csr_matrix: + self.A = sp.csr_matrix(self.A) vals = self.A.data rowptr = self.A.indptr colind = self.A.indices @@ -164,7 +204,7 @@ class Solver(object): :return: x """ if backend is None: backend = DEFAULTS['forward'] - if type(self.A) is not sparse.csr.csr_matrix: + if type(self.A) is not sp.csr.csr_matrix: from scipy.sparse import csr_matrix self.A = csr_matrix(self.A) vals = self.A.data @@ -240,13 +280,13 @@ if __name__ == '__main__': print np.linalg.norm(e-x,np.inf) - n = 6000 + n = 600 A_dense = np.random.random((n,n)) L = np.tril(np.dot(A_dense, A_dense)) # Positive definite is better conditioned. e = np.ones(n) b = np.dot(L, e) - A = sparse.csr_matrix(L) + A = sp.csr_matrix(L) pSolve = Solver(A,flag='L',options={'backend':'python'}); fSolve = Solver(A,flag='L',options={'backend':'fortran'}) tic = time() @@ -257,3 +297,17 @@ if __name__ == '__main__': x = fSolve.solve(b) toc = time() - tic print 'Error Forward Fortran = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc + + + + A = -D*D.T + A[0,0] *= 10 # remove the constant null space from the matrix + e = np.ones(M.nC) + b = A.dot(e) + + iSolve = Solver(A, doDirect=False,options={'M':('GS',A)}) + tic = time() + x = iSolve.solve(b) + toc = time() - tic + print x + print 'Error CG = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc, 'Info: ', iSolve.info From 500844ccb5f182f987355ef212b759be49f66747 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 19 Nov 2013 17:41:51 -0800 Subject: [PATCH 3/3] Brought BFGS into Optimize. updated InexactGaussNewton to use BFGS as a preconditioner. --- SimPEG/inverse/Optimize.py | 114 +++++++++++++++++++++++++++++++++++-- SimPEG/utils/Solver.py | 7 ++- notebooks/BFGS.ipynb | 105 +++++++++++----------------------- 3 files changed, 147 insertions(+), 79 deletions(-) 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": {}