diff --git a/HW1/CheckHW00.py b/HW1/CheckHW00.py new file mode 100644 index 0000000..69f8334 --- /dev/null +++ b/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') diff --git a/HW1/LLTSolver.py b/HW1/LLTSolver.py new file mode 100644 index 0000000..a4957bc --- /dev/null +++ b/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 diff --git a/HW1/PrecCGSolver.py b/HW1/PrecCGSolver.py index 503b2c8..080f58e 100644 --- a/HW1/PrecCGSolver.py +++ b/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)) diff --git a/HW1/benchmarkObjective.py b/HW1/benchmarkObjective.py new file mode 100644 index 0000000..b703f74 --- /dev/null +++ b/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') diff --git a/HW1/incompleteCholesky.py b/HW1/incompleteCholesky.py index 0795298..0ba9330 100644 --- a/HW1/incompleteCholesky.py +++ b/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