new

HW1/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')

HW1/LLTSolver.py

@@ -0,0 +1,65 @@
+# Optimization for Engineers - Dr.Johannes Hild
+# LLT Solver
+
+# Purpose: LLTSolver solves  (L @ L.T)*y=r for y using forward and backward substitution
+
+# Input Definition:
+# L: real valued lower triangle matrix nxn with nonzero diagonal elements
+# r: column vector in R ** n
+# verbose: bool, if set to true, verbose information is displayed
+
+# Output Definition:
+# y: column vector in R ** n (solution in domain space)
+
+# Required files:
+# < none >
+
+# Test cases:
+# L = np.array([[2, 0, 0], [0.5, np.sqrt(15 / 4), 0], [0, 0, 2]], dtype=float)
+# r = np.array([[5], [5], [4]], dtype=float)
+# y = LLTSolver(L,r)
+# should return y = [[1], [1], [1]]
+
+# L = np.array([[22, 0, 0, 0, 0], [17, 13, 0, 0, 0], [13, -2, 17, 0, 0], [8, -4, -7, 18, 0], [4, -5, -4, -5, 19]], dtype=float)
+# r = np.array([[1320], [773], [1192], [132], [1405]], dtype=float)
+# y = LLTSolver(L,r)
+# should return y = [[1],[0],[2],[0],[3]]
+
+import numpy as np
+
+
+def matrnr():
+    # set your matriculation number here
+    matrnr = 0
+    return matrnr
+
+
+def LLTSolver(L: np.array, r: np.array, verbose=0):
+
+    if verbose:
+        print('Start LLTSolver...')
+
+    n = np.size(r)
+    s = r.copy()
+    for i in range(n):
+        for j in range(i):
+            s[i, 0] = s[i, 0] - L[i, j] * s[j, 0]
+
+        if L[i, i] == 0:
+            raise Exception('Zero diagonal element detected...')
+
+        s[i, 0] = s[i, 0] / L[i, i]
+
+    y = s
+
+    for i in range(n-1, -1, -1):
+        for j in range(n-1, i, -1):
+            y[i, 0] = y[i, 0] - L[j, i] * y[j, 0]
+
+        y[i, 0] = y[i, 0] / L[i, i]
+
+    if verbose:
+        residual = (L@L.T)@y-r
+        print('LLTSolver terminated with residual:')
+        print(residual)
+    return y

HW1/PrecCGSolver.py

@@ -44,7 +44,7 @@
 
 def matrnr():
     # set your matriculation number here
-    matrnr = 0
+    matrnr = 23062789
     return matrnr
 
 
@@ -56,14 +56,27 @@ def PrecCGSolver(A: np.array, b: np.array, delta=1.0e-6, verbose=0):
     countIter = 0
 
     L = IC.incompleteCholesky(A)
-    x = LLT.LLTSolver(L, b)
-
-    while MISSING STATEMENT:
-        MISSING CODE
+    x_j = LLT.LLTSolver(L, b)
+
+    r = (A @ x_j) - b
+    d = -LLT.LLTSolver(L, r)
+
+    while np.linalg.norm(r) > delta:
+        # MISSING CODE
+        d_temp = A @ d
+        rho = d.T @ d_temp
+        t = (r.T @ LLT.LLTSolver(L, r)) / rho
+        x_j = x_j + t @ d
+        r_old = r
+        r = r_old + t @ d_temp
+        beta = (r.T @ LLT.LLTSolver(L, r)) / (r_old.T @ LLT.LLTSolver(L, r_old))
+        d = -LLT.LLTSolver(L, r) + beta @ d
         countIter = countIter + 1
         if verbose:
             print('STEP ', countIter, ': norm of residual is ', np.linalg.norm(r))
-
+
+    x = x_j
+
     if verbose:
         print('precCGSolver terminated after ', countIter, ' steps with norm of residual being ', np.linalg.norm(r))
 

HW1/benchmarkObjective.py

@@ -0,0 +1,101 @@
+# Optimization for Engineers - Dr.Johannes Hild
+# Benchmark objective
+# Do not change this file
+
+# Nonlinear function R**3 -> R
+# test function depending on parameters p in R**3 and mapping
+# x -> p[0] * (x[1] + 1) * x[0]**2 + exp(p[1] * x[2] + 1) * x[1]**2 + p[2] * sqrt(x[0] + 1) * x[2]**2
+
+# Class parameters:
+# p: vector in R**3 (parameter space)
+
+# Input Definition:
+# x: vector in R**3 (domain space)
+
+# Output Definition:
+# objective(): real number, evaluation at x for parameters p
+# gradient(): vector in R**3, evaluation of gradient wrt x
+# hessian(): matrix in R**3x3, evaluation of hessian wrt x
+# setParameters(): sets p
+# parameterGradient(): vector in R**3, evaluation of gradient wrt p
+
+# Required files:
+# < none >
+
+# Test cases:
+# myObjective = benchmarkObjective([3,2,16]).objective(np.array([[0],[0],[-1 / 2]], dtype=float))
+# should return
+# myObjective = 4
+
+# myGradient = benchmarkObjective([3,2,16]).gradient(np.array([[0],[0],[-1 / 2]], dtype=float))
+# should return
+# myGradient = [[2],[0],[-16]]
+
+# myHessian = benchmarkObjective([3,2,16]).hessian(np.array([[0],[0],[-1 / 2]], dtype=float))
+# should return
+# myHessian = [[5, 0, -8],[0, 2, 0],[-8, 0, 32]]
+
+import numpy as np
+
+
+def matrnr():
+    # set your matriculation number here
+    matrnr = 0
+    return matrnr
+
+
+class benchmarkObjective:
+
+    def __init__(self, p: np.array):
+        self.p = p
+
+    def objective(self, x: np.array):
+        checkargs(x)
+        f = self.p[0, 0] * (x[1, 0] + 1) * x[0, 0] ** 2 + np.exp(self.p[1, 0] * x[2, 0] + 1) * x[1, 0] ** 2 \
+            + self.p[2, 0] * np.sqrt(x[0, 0] + 1) * x[2, 0] ** 2
+        return f
+
+    def gradient(self, x: np.array):
+        checkargs(x)
+        f_dx0 = 2 * self.p[0, 0] * (x[1, 0] + 1) * x[0, 0] + self.p[2, 0] / 2 * x[2, 0] ** 2 / np.sqrt(x[0, 0] + 1)
+        f_dx1 = self.p[0, 0] * x[0, 0] ** 2 + 2 * np.exp(self.p[1, 0] * x[2, 0] + 1) * x[1, 0]
+        f_dx2 = self.p[1, 0] * np.exp(self.p[1, 0] * x[2, 0] + 1) * x[1, 0] ** 2 + 2 * self.p[2, 0] * np.sqrt(x[0, 0] + 1) * x[2, 0]
+        g = np.array([[f_dx0], [f_dx1], [f_dx2]])
+        return g
+
+    def hessian(self, x: np.array):
+        checkargs(x)
+        f_dx00 = 2 * self.p[0, 0] * (x[1, 0] + 1) - self.p[2, 0] / 4 * x[2, 0] ** 2 / np.sqrt((x[0, 0] + 1) ** 3)
+        f_dx01 = 2 * self.p[0, 0] * x[0, 0]
+        f_dx02 = self.p[2, 0] * x[2, 0] / np.sqrt(x[0, 0] + 1)
+        f_dx11 = 2 * np.exp(self.p[1, 0] * x[2, 0] + 1)
+        f_dx12 = 2 * self.p[1, 0] * np.exp(self.p[1, 0] * x[2, 0] + 1) * x[1, 0]
+        f_dx22 = self.p[1, 0] ** 2 * np.exp(self.p[1, 0] * x[2, 0] + 1) * x[1, 0] ** 2 + 2 * self.p[2, 0] * np.sqrt(x[0, 0] + 1)
+        h = np.array([[f_dx00, f_dx01, f_dx02], [f_dx01, f_dx11, f_dx12], [f_dx02, f_dx12, f_dx22]])
+        return h
+
+    def setParameters(self, p: np.array):
+        self.p = p
+
+    def parameterGradient(self, x: np.array):
+        R_dp1 = (x[1, 0] + 1) * x[0, 0] ** 2
+        R_dp2 = x[2, 0] * np.exp(self.p[1, 0] * x[2, 0] + 1) * x[1, 0] ** 2
+        R_dp3 = np.sqrt(x[0, 0] + 1) * x[2, 0] ** 2
+
+        myGradP = np.array([[R_dp1], [R_dp2], [R_dp3]], dtype=float)
+
+        return myGradP
+
+    @staticmethod
+    def getXData():
+        xdata = np.array([[0.0, 8.0, 0.0, 8.0, 0.0, 8.0, 0.0, 8.0, 4.0, 0.0, 8.0, 4.0, 4.0, 4.0, 4.0], [-4.0, -4.0, 4.0, 4.0, -4.0, -4.0, 4.0, 4.0, 0.0, 0.0, 0.0, -4.0, 4.0, 0.0, 0.0], [-1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, 1.0]])
+        return xdata
+
+    @staticmethod
+    def getFData():
+        fdata = np.array([[22.0, -522.0, 22.0, 1014.0, 337.0, -207.0, 337.0, 1329.0, 48.0, 0.0, 192.0, -101.0, 283.0, 84.0, 84.0]])
+        return fdata
+
+def checkargs(x: np.array):
+    if x[0] <= -1:
+        raise ValueError('x[0] is not allowed to be smaller than -1')

HW1/incompleteCholesky.py

@@ -52,7 +52,7 @@
 
 def matrnr():
     # set your matriculation number here
-    matrnr = 0
+    matrnr = 23062789
     return matrnr
 
 
@@ -74,8 +74,29 @@ def incompleteCholesky(A: np.array, alpha=1.0e-3, delta=1.0e-6, verbose=0):
 
     if verbose:
         print('Start incompleteCholesky...')
-
-    MISSING CODE
+
+    # MISSING CODE    
+
+    for k in range(n):
+        L[k,k] = np.sqrt(max(L[k,k],alpha))
+
+        for i in range(k+1,n):
+            if (np.abs(L[i,k]) > delta):
+                L[i,k] = (L[i,k] / L[k,k])
+            else:
+                L[i,k] = 0
+
+        for j in range(k+1,n):
+            for i in range(j,n):
+                if (np.abs(L[i,j]) > delta):
+                    # print("sebze")
+                    L[i,j] = L[i,j] - L[i,k]*L[j,k]
+
+    for i in range(n):
+        for j in range(i+1,n):
+            L[i,j] = 0
+
+    # L = np.copy(A)
 
     if verbose:
         residualmatrix = A - L @ L.T