Merge branch 'BFGS' of https://bitbucket.org/rcockett/simpeg into richards

This commit is contained in:
Rowan Cockett
2013-11-19 17:46:48 -08:00
5 changed files with 348 additions and 26 deletions
+115 -9
View File
@@ -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)
@@ -586,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)
@@ -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
@@ -620,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
+2 -1
View File
@@ -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 {
+4 -4
View File
@@ -276,16 +276,16 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
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
+69 -12
View File
@@ -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,44 @@ 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
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())
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 +150,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 +176,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 +207,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 +283,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 +300,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
+158
View File
@@ -0,0 +1,158 @@
{
"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": 2
},
{
"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",
"x0 = np.array([1,0])\n",
"opt = inverse.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"
]
},
{
"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": 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": {}
}
]
}