HW4 Done

HW4/BFGSDescent.py

@@ -0,0 +1,95 @@
+# Optimization for Engineers - Dr.Johannes Hild
+# global BFGS descent
+
+# Purpose: Find xmin to satisfy norm(gradf(xmin))<=eps
+# Iteration: x_k = x_k + t_k * d_k
+# d_k is the BFGS direction. If a descent direction check fails, d_k is set to steepest descent and the inverse BFGS matrix is reset.
+# t_k results from Wolfe-Powell
+
+# Input Definition:
+# f: objective class with methods .objective() and .gradient()
+# x0: column vector in R ** n(domain point)
+# eps: tolerance for termination. Default value: 1.0e-3
+# verbose: bool, if set to true, verbose information is displayed
+
+# Output Definition:
+# xmin: column vector in R ** n(domain point)
+
+# Required files:
+# t = WolfePowellSearch(f, x, d) from WolfePowellSearch.py
+
+# Test cases:
+# myObjective = nonlinearObjective()
+# x0 = np.array([[-0.01], [0.01]])
+# eps = 1.0e-6
+# xmin = BFGSDescent(myObjective, x0, eps, 1)
+# should return
+# xmin close to [[0.26],[-0.21]] with the inverse BFGS matrix being close to [[0.0078, 0.0005], [0.0005, 0.0080]]
+
+# myObjective = nonlinearObjective()
+# x0 = np.array([[0.6], [-0.6]])
+# eps = 1.0e-3
+# xmin = BFGSDescent(myObjective, x0, eps, 1)
+# should return
+# xmin close to [[-0.26],[0.21]] with the inverse BFGS matrix being close to [[0.0150, 0.0012], [0.0012, 0.0156]]
+
+# myObjective = bananaValleyObjective()
+# x0 = np.array([[0], [1]])
+# eps = 1.0e-6
+# xmin = BFGSDescent(myObjective, x0, eps, 1)
+# should return
+# xmin close to [[1],[1]] in less than 100 iterations with the inverse BFGS matrix being almost singular and close to [[0.4996, 0.9993], [0.9993, 2.0040]]
+
+
+import numpy as np
+import WolfePowellSearch as WP
+
+
+def matrnr():
+    # set your matriculation number here
+    matrnr = 23062789
+    return matrnr
+
+
+def BFGSDescent(f, x0: np.array, eps=1.0e-3, verbose=0):
+    if eps <= 0:
+        raise TypeError('range of eps is wrong!')
+
+    if verbose:
+        print('Start BFGSDescent...')
+
+    countIter = 0
+
+    x = x0
+    B = np.eye(x.shape[0])
+    def update(xk):
+        gradx = f.gradient(xk)
+        gradnormx = np.linalg.norm(gradx)
+        return gradx, gradnormx
+
+    gradx, gradnormx = update(x)
+    while gradnormx > eps:
+        dk = -np.dot(B, gradx)
+        descent = gradx.T @ dk
+        if descent > 0:
+            dk = -gradx
+            B = np.eye(x.shape[0])
+        t = WP.WolfePowellSearch(f, x, dk)
+        ex = x.copy()
+        ex_grad = gradx.copy()
+        x = x + t*dk
+        gradx, gradnormx = update(x)
+        dg = gradx - ex_grad
+        dx = x - ex
+
+        rk = dx - (B @ dg)
+        B += ((rk @ dx.T) + (dx @ rk.T)) / (dg.T @ dx)
+        B -= (rk.T @ dg) * (dx @ dx.T) / ((dg.T @ dx) @ (dg.T @ dx))
+        countIter += 1
+
+    if verbose:
+        gradx = f.gradient(x)
+        print('BFGSDescent terminated after ', countIter, ' steps with norm of gradient =', np.linalg.norm(gradx), 'and the inverse BFGS matrix is')
+        print(B)
+
+    return x

HW4/CheckHW00.py

@@ -0,0 +1,54 @@
+# Optimization for Engineers - Dr.Johannes Hild
+# Mock Homework to check setup
+# Do not change this file
+
+print('Welcome to Optimization for Engineers.\n')
+print('If this script fails, then your setup is not working correctly.')
+print('First we check if the math package numpy is installed.')
+
+import numpy as np
+
+X = np.power(2, 3)
+Y = 2**3
+if X == Y:
+    print('numpy seems to work.\n')
+
+print('Next we check if the function definitions in BenchmarkObjective.py are available.')
+
+import benchmarkObjective as BM
+
+p = np.array([[3], [2], [16]])
+x = np.array([[1.05], [-0.54], [-0.03]])
+myObjective = BM.benchmarkObjective(p)
+f = myObjective.objective(x)
+fg = myObjective.gradient(x)
+fh = myObjective.hessian(x)
+xdata = myObjective.getXData()
+fdata = myObjective.getFData()
+
+print('The benchmark objective without noise returns')
+print(f)
+print('at x and the gradient is')
+print(fg)
+print('and the hessian is')
+print(fh,'\n')
+
+print('The benchmark objective measure points are')
+print(xdata)
+print('with measure results')
+print(fdata)
+
+import LLTSolver as LLT
+L = np.array([[3, 0], [1, 2]], dtype=float)
+
+print('Now we use LLTSolver. The matrix A = ')
+print(L @ L.T)
+print('can obviously be decomposed into A = L @ L.T with triangle matrix L = ')
+print(L)
+
+print('We can use this to quickly solve A @ x = [[14], [9]] by calling LLTSolver. Solution: x = ') 
+x = LLT.LLTSolver(L, np.array([[14], [9]]))
+print(x)
+
+
+print('\nEverything seems to be fine, please return your files in StudOn')

HW4/CheckHW01.py

@@ -0,0 +1,82 @@
+# Optimization for Engineers - Dr.Johannes Hild
+# Programming Homework Check Script
+# Do not change this file
+
+print('Welcome to Optimization for Engineers.\n')
+print('If this script fails, then your programming homework is not working correctly.')
+
+import numpy as np
+import incompleteCholesky as IC
+
+print('Checking IncompleteCholesky...')
+alpha = 0
+delta = -1
+verbose = 1
+A = np.array([[5, 4, 3, 2, 1], [4, 5, 2, 1, 0], [3, 2, 5, 0, 0], [2, 1, 0, 5, 0], [1, 0, 0, 0, 5]], dtype=float)
+print('If you get a warning about negative delta, this is okay!')
+L = IC.incompleteCholesky(A, alpha, delta, verbose)
+if np.linalg.norm(A - L @ L.T) > 1.0e-3:
+    raise Exception('Your full Cholesky decomposition is not working correctly.')
+else:
+    print('Check 01 okay')
+
+B = np.array([[4, 1, 0], [1, 4, 0], [0, 0, 4]], dtype=float)
+LBe = np.array([[2, 0, 0], [0.5, 1.9364, 0], [0, 0, 2]], dtype=float)
+LB = IC.incompleteCholesky(B)
+if np.linalg.norm(LB - LBe) > 1.0e-3:
+    raise Exception('Your incomplete Cholesky decomposition is not working correctly.')
+else:
+    print('Check 02 okay')
+
+alpha = 4
+LCe = np.array([[2, 0, 0], [0.5, 2, 0], [0, 0, 2]], dtype=float)
+LC = IC.incompleteCholesky(B, alpha)
+if np.linalg.norm(LC - LCe) > 1.0e-3:
+    raise Exception('Your incomplete Cholesky decomposition is not working correctly for alpha = 4.')
+else:
+    print('Check 03 okay')
+
+alpha = 1.0e-3
+delta = 1
+LDe = np.array([[2, 0, 0], [0, 2, 0], [0, 0, 2]], dtype=float)
+LD = IC.incompleteCholesky(B, alpha, delta)
+if np.linalg.norm(LD - LDe) > 1.0e-3:
+    raise Exception('Your incomplete Cholesky decomposition is not working correctly for delta = 1.')
+else:
+    print('Check 04 okay')
+
+print('Checking PrecCGSolver...')
+
+import PrecCGSolver as PCG
+
+A = np.array([[4, 1, 0], [1, 7, 0], [ 0, 0, 3]], dtype=float)
+b = np.array([[5], [8], [3]], dtype=float)
+delta = 1.0e-6
+x = PCG.PrecCGSolver(A, b, delta, 1)
+xe = np.array([[1], [1], [1]])
+if np.linalg.norm(x - xe) > 1.0e-3:
+    raise Exception('Your PrecCGSolver is not working correctly.')
+else:
+    print('Check 05 okay')
+
+A = np.array([[484, 374, 286, 176, 88], [374, 458, 195, 84, 3], [286, 195, 462, -7, -6], [176, 84, -7, 453, -10], [88, 3, -6, -10, 443]], dtype=float)
+b = np.array([[1320], [773], [1192], [132], [1405]], dtype=float)
+x = PCG.PrecCGSolver(A, b, delta, 1)
+xe = np.array([[1], [0], [2], [0], [3]])
+
+if np.linalg.norm(x - xe) > 1.0e-3:
+    raise Exception('Your PrecCGSolver is not working correctly for other dimensions.')
+else:
+    print('Check 06 okay')
+
+A = np.array([[1, 2, 3, 4], [2, 4, -100, -100], [3, -100, 7, 0], [4, -100, 0, 3]], dtype=float)
+b = np.array([[0], [5], [8], [3]], dtype=float)
+x = PCG.PrecCGSolver(A, b, delta, 1)
+print('Check 07 is okay, if you got a message that PrecCGSolver took more than 4 iterations')
+
+if IC.matrnr() == 0:
+    raise Exception('Please set your matriculation number in incompleteCholesky.py!')
+elif PCG.matrnr() == 0:
+    raise Exception('Please set your matriculation number in PrecCGSolver.py!')
+else:
+    print('Everything seems to be fine, please return your files in StudOn')

HW4/CheckHW02.py

@@ -0,0 +1,86 @@
+# Optimization for Engineers - Dr.Johannes Hild
+# Programming Homework Check Script
+# Do not change this file
+
+print('Welcome to Optimization for Engineers.\n')
+print('If this script fails, then your programming homework is not working correctly.')
+
+import numpy as np
+import simpleValleyObjective as SO
+import nonlinearObjective as NO
+import WolfePowellSearch as WP
+
+p = np.array([[0], [1]])
+myObjective = SO.simpleValleyObjective(p)
+x = np.array([[-1.01], [1]])
+d = np.array([[1], [1]])
+sigma = 1.0e-3
+rho = 1.0e-2
+t = WP.WolfePowellSearch(myObjective, x, d, sigma, rho, 1)
+te = 1
+if t != te:
+    raise Exception('Your Wolfe-Powell search is not recognizing t = 1 as valid starting point.')
+else:
+    print('Check 01 okay')
+
+x = np.array([[-1.2], [1]])
+d = np.array([[0.1], [1]])
+t = WP.WolfePowellSearch(myObjective, x, d, sigma, rho, 1)
+te = 16
+if t != te:
+    raise Exception('Your Wolfe-Powell search is not fronttracking correctly.')
+else:
+    print('Check 02 okay')
+
+x = np.array([[-0.2], [1]])
+d = np.array([[1], [1]])
+t = WP.WolfePowellSearch(myObjective, x, d, sigma, rho, 1)
+te = 0.25
+if t != te:
+    raise Exception('Your Wolfe-Powell search is not refining correctly.')
+else:
+    print('Check 03 okay')
+
+myObjective = NO.nonlinearObjective()
+x = np.array([[0.53], [-0.29]])
+d = np.array([[-3.88], [1.43]])
+t = WP.WolfePowellSearch(myObjective, x, d, sigma, rho, 1)
+te = 0.0938
+if np.abs(t-te) > 1.0e-3:
+    raise Exception('Your Wolfe-Powell search is not working for general objective class.')
+else:
+    print('Check 04 okay')
+
+import globalNewtonDescent as GN
+
+x0 = np.array([[-0.01], [0.01]])
+eps = 1.0e-6
+xmin = GN.globalNewtonDescent(myObjective, x0, eps, 1)
+xe = np.array([[0.26], [-0.21]])
+if np.linalg.norm(xmin-xe) > 1.0e-2:
+    raise Exception('Your global Newton Descent is not working correctly.')
+else:
+    print('Check 05 okay')
+
+x0 = np.array([[-0.6], [0.6]])
+xmin = GN.globalNewtonDescent(myObjective, x0, eps, 1)
+xe = np.array([[-0.26], [0.21]])
+if np.linalg.norm(xmin-xe) > 1.0e-2:
+    raise Exception('Your global Newton Descent walks a wrong path, maybe you switch to steepest descent too often?')
+else:
+    print('Check 06 okay')
+
+x0 = np.array([[0.6], [-0.6]])
+xmin = GN.globalNewtonDescent(myObjective, x0, eps, 1)
+xe = np.array([[-0.26], [0.21]])
+if np.linalg.norm(xmin - xe) > 1.0e-2:
+    raise Exception('Your global Newton Descent walks a wrong path, maybe you make mistakes in choosing the descent directions?')
+else:
+    print('Check 07 okay')
+
+if WP.matrnr() == 0:
+    raise Exception('Please set your matriculation number in WolfePowellSearch.py!')
+elif GN.matrnr() == 0:
+    raise Exception('Please set your matriculation number in globalNewtonDescent.py!')
+else:
+    print('Everything seems to be fine, please return your files in StudOn')