wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-06-30 16:50:15 +08:00
Code Issues Packages Projects Releases Wiki Activity

added derivative check

Browse Source
This commit is contained in:
Lars Ruthotto
2013-06-19 17:40:30 -07:00
parent 70a143031b
commit 3ef293b19e
1 changed files with 54 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files
code/GaussNewton.py
+54 -9
View File
@@ -1,8 +1,9 @@
import numpy as np
import matplotlib.pyplot as plt
from pylab import norm
def GaussNewton(fctn, x0,maxIter=20, maxIterLS=10, LSreduction=1e-4, tolJ=1e-3, tolX=1e-3,
tolG=1e-3, eps=1e-16, xStop=np.empty):
tolG=1e-3, eps=1e-16, xStop=[]):
"""
GaussNewton Optimization
@@ -21,7 +22,7 @@ def GaussNewton(fctn, x0,maxIter=20, maxIterLS=10, LSreduction=1e-4, tolJ=1e-3,
print "%s" % '-'*57
# evaluate stopping criteria
if xStop==np.empty:
if xStop==[]:
xStop=x0
Jstop = fctn(xStop)
print "%3d\t%1.2e" % (-1, Jstop[0])
@@ -69,11 +70,11 @@ def GaussNewton(fctn, x0,maxIter=20, maxIterLS=10, LSreduction=1e-4, tolJ=1e-3,
iter = iter +1
print "%s STOP! %s" % ('-'*25,'-'*25)
print "%d : |Jc-Jold|=%1.4e <= tolJ*(1+|Jstop|)=%1.4e" % (STOP[0],abs(Jc[0]-Jold[0]),tolJ*(1+abs(Jstop[0])))
print "%d : |xc-xOld|=%1.4e <= tolX*(1+|x0|) =%1.4e" % (STOP[1],norm(xc-xOld),tolX*(1+norm(x0)))
print "%d : |dJ| =%1.4e <= tolG*(1+|Jstop|)=%1.4e" % (STOP[2],norm(dJ),tolG*(1+abs(Jstop[0])))
print "%d : |dJ| =%1.4e <= 1e3*eps =%1.4e" % (STOP[3],norm(dJ),1e3*eps)
print "%d : iter =%d\t\t <= maxIter =%d" % (STOP[4],iter,maxIter)
print "%d : |Jc-Jold| = %1.4e <= tolJ*(1+|Jstop|) = %1.4e" % (STOP[0],abs(Jc[0]-Jold[0]),tolJ*(1+abs(Jstop[0])))
print "%d : |xc-xOld| = %1.4e <= tolX*(1+|x0|) = %1.4e" % (STOP[1],norm(xc-xOld),tolX*(1+norm(x0)))
print "%d : |dJ| = %1.4e <= tolG*(1+|Jstop|) = %1.4e" % (STOP[2],norm(dJ),tolG*(1+abs(Jstop[0])))
print "%d : |dJ| = %1.4e <= 1e3*eps = %1.4e" % (STOP[3],norm(dJ),1e3*eps)
print "%d : iter = %3d\t <= maxIter\t = %3d" % (STOP[4],iter,maxIter)
print "%s DONE! %s\n" % ('='*25,'='*25)
return xc
@@ -82,15 +83,59 @@ def Rosenbrock(x):
"""
Rosenbrock function for testing GaussNewton scheme
"""
J = 100*(x[1]-x[0]**2)**2+(1-x[0])**2
dJ = np.array([-400*(x[1]*x[0]-x[0]**3)-2*(1-x[0]),200*(x[1]-x[0]**2)])
J = 100*(x[1]-x[0]**2)**2+(1-x[0])**2
dJ = np.array([2*(200*x[0]**3-200*x[0]*x[1]+x[0]-1),200*(x[1]-x[0]**2)])
H = np.array([[-400*x[1]+1200*x[0]**2+2, -400*x[0]],[ -400*x[0], 200]],dtype=float);
return J,dJ,H
def checkDerivative(fctn,x0):
"""
Basic derivative check
Compares error decay of 0th and 1st order Taylor approximation at point
x0 for a randomized search direction.
Input:
------
fctn - function handle
x0 - point at which to check derivative
"""
print "%s checkDerivative %s" % ('='*20,'='*20)
print "iter\th\t\t|J0-Jt|\t\t|J0+h*dJ'*dx-Jt|"
Jc,dJ,H = fctn(x0)
dx = np.random.randn(len(x0),1)
t = np.logspace(-1,-10,10)
E0 = np.zeros(t.shape)
E1 = np.zeros(t.shape)
for i in range(0,10):
Jt = fctn(x0+t[i]*dx)
E0[i] = norm(Jt[0]-Jc[0]) # 0th order Taylor
E1[i] = norm(Jt[0]-Jc[0]-t[i]*np.dot(dJ.T,dx)) # 1st order Taylor
print "%d\t%1.2e\t%1.3e\t%1.3e" % (i,t[i],E0[i],E1[i])
print "%s DONE! %s\n" % ('='*25,'='*25)
plt.figure()
plt.clf()
plt.loglog(t,E0,'b')
plt.loglog(t,E1,'g--')
plt.title('checkDerivative')
plt.xlabel('h')
plt.ylabel('error of Taylor approximation')
plt.legend(['0th order', '1st order'],loc='upper left')
plt.show()
return
if __name__ == '__main__':
x0 = np.array([[2.6],[3.7]])
fctn = lambda x:Rosenbrock(x)
checkDerivative(fctn,x0)
xOpt = GaussNewton(fctn,x0,maxIter=20)
print "xOpt=[%f,%f]" % (xOpt[0],xOpt[1])