new file: drawFunc.py

	new file:   drawFunc.pyc
	new file:   func.py
	new file:   func.pyc
	new file:   main.py
	new file:   optMethods.py
	new file:   optMethods.pyc

drawFunc.py

@@ -0,0 +1,19 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from func import *
+
+def draw(f_x):
+# allocate grids
+    dx = np.arange(-1.5, 2, .02)
+    dy = np.arange(-3.0, 3.5, .02)
+    X, Y = np.meshgrid(dx, dy)
+
+# eval function f_x
+    Z = [f_x(np.matrix([[x1], [x2]])) \
+        for (x1, x2) in zip(np.hstack(X), np.hstack(Y))]
+    Z = np.array(Z)
+    Z = Z.reshape(np.shape(X))
+
+# take this pencil
+    cs = plt.contour(X, Y, Z, 20, colors='k')
+    plt.clabel(cs, inline=1)

func.py

@@ -0,0 +1,46 @@
+import numpy as np
+
+class quadratic:
+
+# trivial (solving linear system)
+    A = np.zeros([2, 2])
+    b = np.zeros([2, 1])
+
+    def __init__(self):
+        pass
+
+    def __init__(self, A, b):
+        self.A = A
+        self.b = b
+
+    def f_x(self, x):
+        f = 0.5 * np.dot(np.dot(x.T, self.A), x) - np.dot(self.b.T, x)
+        return f[0, 0]
+
+    def g_x(self, x):
+        return np.dot(self.A, x) - self.b
+
+    def G_x(self, x):
+        return np.asmatrix(self.A)
+
+
+class rosenbrock:
+
+    def __init__(self):
+        pass
+
+    def f_x(self, x):
+        f = 100 * (x[1] - (x[0])**2)**2 + (1-x[0])*2
+        return f[0, 0]
+
+    def g_x(self, x):
+        g1 = -400.0 * (x[0]*x[1] - x[0]**3) - 2.0 + 2.0*x[0]
+        g2 = 200.0 * (x[1] - x[0]**2)
+        return np.matrix([[g1[0,0]], [g2[0,0]]])
+
+    def G_x(self, x):
+        G11 = 1200.0*x[0]**2 - 400.0*x[1] + 2
+        G12 = -400.0*x[0]
+        G21 = -400.0*x[0]
+        G22 = 200.0
+        return np.matrix([[G11[0,0], G12[0,0]], [G21[0,0], G22]])

main.py

@@ -0,0 +1,63 @@
+import optMethods as opt
+import drawFunc
+import func
+import numpy as np
+import matplotlib.pyplot as plt
+import pdb
+
+def quadratic():
+
+    A = np.matrix([
+        [3., 2.],
+        [2., 6.]])
+    b = np.matrix([[1.], [-5.]])
+    obj = func.quadratic(A,b)
+
+    plt.figure()
+    drawFunc.draw(obj.f_x)
+
+    x = np.matrix([[0], [2]])
+    track = opt.newton(x, obj, 0, 1) # just one iteration
+    plt.plot(track[0,:].tolist()[0], track[1,:].tolist()[0])
+
+    plt.show()
+
+def rosen():
+    obj = func.rosenbrock()
+
+    plt.figure()
+    drawFunc.draw(obj.f_x)
+    start_point = np.matrix([[-1.0], [2.0]])
+    plt.annotate('Start', xy=(-1, 2), xytext=(-1.2, 2.2),
+            arrowprops=dict(facecolor='black', shrink=0.05))
+    plt.annotate('Optimal', xy=(1, 1), xytext=(0.5, 1.2),
+            arrowprops=dict(facecolor='black', shrink=0.05))
+
+    v = 0.0
+    str1 = "Newton"
+    track = opt.newton(start_point, obj, v)
+    p1, = plt.plot(track[0,:].tolist()[0], track[1,:].tolist()[0],
+            'o-', linewidth=2.0)
+
+    v = 0.1
+    str2 = "Modified_Newton %.2f"%v
+    track = opt.newton(start_point, obj, v)
+    p2, = plt.plot(track[0,:].tolist()[0], track[1,:].tolist()[0],
+            'o-', linewidth=2.0)
+
+    v = 1.5
+    str3 = "Modified_Newton %.2f"%v
+    track = opt.newton(start_point, obj, v, 50)
+    p3, = plt.plot(track[0,:].tolist()[0], track[1,:].tolist()[0],
+            'o-', linewidth=2.0)
+
+    str4 = "Quasi_Newton_Rank_One"
+    track = opt.quasi_newton(start_point, obj, 100)
+    p4, = plt.plot(track[0,:].tolist()[0], track[1,:].tolist()[0],
+            'o-', linewidth=2.0)
+
+    plt.legend([p1, p2, p3, p4], [str1, str2, str3, str4])
+    plt.show()
+
+#quadratic()
+rosen()

optMethods.py

@@ -0,0 +1,70 @@
+import numpy as np
+
+def newton(start_point, obj_fun, modified=0, iterations=10):
+
+# with second order information
+    x = start_point
+    track = x
+
+    k = 0
+    while k < iterations:
+
+        if modified > 0:
+            vI = modified * np.matrix([[1, 0], [0, 1]])
+            G_ = np.linalg.inv(obj_fun.G_x(x) + vI) # inverse?
+        else:
+            G_ = np.linalg.inv(obj_fun.G_x(x)) # inverse?
+        g_ = obj_fun.g_x(x)
+        delta = -1.0 * np.dot(G_, g_)
+        x = x + delta
+
+        track = np.concatenate((track, x), axis=1)
+        k += 1
+
+    return track
+
+def quasi_newton(start_point, obj_fun, iteration=10):
+
+# with first order information
+    x = start_point
+    track = x
+
+    k = 0
+    H = np.matrix([[1, 0], [0, 1]])
+
+    while k < iteration:
+
+        s = -1.0 * np.dot(H, obj_fun.g_x(x))
+
+        alpha_k = _backtracking_line_search(obj_fun, x, s)
+        delta_k = alpha_k * s
+        x_k_1 = x + delta_k
+
+        gamma_k = obj_fun.g_x(x_k_1) - obj_fun.g_x(x)
+        u = delta_k - np.dot(H, gamma_k)
+
+        scale_a = np.dot(u.T, gamma_k)
+        if scale_a == 0:
+            scale_a = 0.000001
+        H = H + np.dot(u, u.T) / scale_a
+
+        track = np.concatenate((track, x), axis=1)
+        x = x_k_1
+        k += 1
+
+    return track
+
+def _backtracking_line_search(obj_fun, x, s):
+    alpha = 1.0
+    p = 0.86 # magic number
+    c = 0.0001
+    diff = 1
+
+    i = 0
+    while (diff > 0 & i < 10000):
+        f = obj_fun.f_x(x + alpha * s)
+        f_x = obj_fun.f_x(x)
+        diff = f - f_x - c * alpha * np.dot(obj_fun.g_x(x).T, s)
+        alpha *= p
+        i += 1
+    return alpha