From: Bjørn Rustad <rustadbjornen@gmail.com>
Date: Wed, 28 Sep 2011 12:57:31 +0000 (+0200)
Subject: Improvements
X-Git-Url: http://git.rustad.me/?a=commitdiff_plain;h=bf4d23318c741cef778fda2de1cdcfabf951d9db;p=nummat1

Improvements
---

diff --git a/function.py b/function.py
index 9d352e0..6a4dafe 100644
--- a/function.py
+++ b/function.py
@@ -1,17 +1,21 @@
 import numpy as np
 
+# Calculate function value in point x
 def function_g(x, H, b, c):
     return -b.T.dot(x) + 0.5 * x.T.dot(H.dot(x)) + (1/12) * x.T.dot(_C(x,c).dot(x))
 
+# Create matrix C(x)
 def _C(x, c):
     A = np.zeros( (len(c), len(c)) )
     for i in xrange(0, len(c)):
         A[i][i] = c[i]*(x[i])**2
     return A
 
+# Calculate gradient of g in point x
 def gradient_g(x, H, b, c):
     return -b.T + 0.5 * H.dot(x) + 0.5 * H.T.dot(x) + (1/3) * x.T.dot(_C(x,c))
 
+# Calculate hessian of g in point x
 def hessian_g(x, H, b, c):
     return 0.5 * H.T + 0.5 * H + _C(x,c)
 
diff --git a/generate.py b/generate.py
index 7013e28..4af4c0b 100644
--- a/generate.py
+++ b/generate.py
@@ -6,7 +6,7 @@ def spdmatrix(n, limit):
     A = np.zeros( (n,n) )
     for y in xrange(0, n):
         A[y][y] = math.sqrt(random.uniform(0, limit))
-        for x in xrange(y+1, n):
+        for x in xrange(0, y):
             A[y][x] = random.uniform(-limit, limit)
     return A.T.dot(A)
 
diff --git a/steepest.py b/steepest.py
index 7e10fd2..857cb15 100644
--- a/steepest.py
+++ b/steepest.py
@@ -1,39 +1,72 @@
 from enthought.mayavi import mlab
-from function import function_g, gradient_g
+from function import function_g, gradient_g, _C
 import generate
 import surf
 import numpy as np
 import matplotlib.pyplot as plt
+import time
 from numpy.linalg import norm
 
 def steepest_iter(x, alpha, grad_f):
     return x - alpha * grad_f(x)
 
-H = generate.spdmatrix(2, 3)
-b = generate.vector(2, -1, 1)
-c = generate.vector(2, -3, 3)
+def optimal_step(x, H, b, c, grad_f):
+    pk = -grad_f(x)
+    Cpk = _C(pk, c)
+    Cxk = _C(x, c)
 
-alpha = 0.1
-xiter = np.array([0, 0])
+    coeff = [0] * 4
+
+    coeff[0] = (1.0/3.0) * pk.T.dot(Cpk.dot(pk))
+    coeff[1] = x.T.dot(Cpk.dot(pk))
+    coeff[2] = pk.T.dot(H.dot(pk)) + pk.T.dot(Cxk.dot(pk))
+    coeff[3] = (1.0/3.0) * pk.T.dot(Cxk.dot(x)) + pk.T.dot(H.dot(x)) - b.T.dot(pk)
+
+    roots = np.roots(coeff)
+
+    for a in roots:
+        if abs(a.imag) < 0.0000001:
+            return a.real
+
+    return 0.1
+
+#H = generate.spdmatrix(2, 3)
+#b = generate.vector(2, 0.01, 1)
+#c = generate.vector(2, 0.01, 3)
+
+H = np.array([ (2.06810668,-1.51172618),  (-1.51172618, 3.95629853)])
+b = np.array([-0.09229406, 0.38657143])
+c = np.array([ 0.12659127, 2.42818636])
+
+xiter = np.array([1.5, 1.5])
 x0_norm = norm(gradient_g(xiter, H, b, c))
 
-tolerance = 0.01
+tolerance = 0.001
 
 relative_residual = norm(gradient_g(xiter, H, b, c)) / x0_norm
 relative_residuals = []
 function_values = []
 
-while relative_residual > tolerance:
+maxiter = 0
+
+while relative_residual > tolerance and maxiter < 100:
+    maxiter += 1
+
+    alpha = optimal_step(xiter, H, b, c, lambda x: gradient_g(x, H, b, c))
+
     relative_residual = norm(gradient_g(xiter, H, b, c)) / x0_norm
     relative_residuals.append(relative_residual)
     function_values.append(function_g(xiter, H, b, c))
 
     xiter = steepest_iter(xiter, alpha, lambda x: gradient_g(x, H, b, c))
+
     z = function_g(xiter, H, b, c)
-    p = mlab.points3d([xiter[0]], [xiter[1]], [z], scale_factor=0.05)
+    p = mlab.points3d([xiter[0]], [xiter[1]], [z], scale_factor=0.05,
+            color=(maxiter % 2, maxiter % 2, maxiter % 2))
+
     print xiter
 
-X, Y, Z = surf.surf(H, b, c, xiter[0]-2, xiter[0]+2, xiter[1]-2, xiter[1]+2, -10, 3)
+X, Y, Z = surf.surf(H, b, c, xiter[0]-1, xiter[0]+1, xiter[1]-1, xiter[1]+1, -10, 3)
 
 plt.figure(1)
 ax = plt.subplot(211)