diff --git a/HW1/CheckHW01.py b/HW1/CheckHW01.py new file mode 100644 index 0000000..6b0e9ce --- /dev/null +++ b/HW1/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') diff --git a/HW1/PrecCGSolver.py b/HW1/PrecCGSolver.py new file mode 100644 index 0000000..503b2c8 --- /dev/null +++ b/HW1/PrecCGSolver.py @@ -0,0 +1,70 @@ +# Optimization for Engineers - Dr.Johannes Hild +# Preconditioned Conjugate Gradient Solver + +# Purpose: CGSolver finds y such that norm(A * y - b) <= delta using incompleteCholesky as preconditioner + +# Input Definition: +# A: real valued matrix nxn +# b: column vector in R ** n +# delta: positive value, tolerance for termination. Default value: 1.0e-6. +# verbose: bool, if set to true, verbose information is displayed + +# Output Definition: +# x: column vector in R ^ n(solution in domain space) + +# Required files: +# L = incompletecholesky(A, 1.0e-3, delta) from IncompleteCholesky.py +# y = LLTSolver(L, r) from LLTSolver.py + + +# Test cases: +# 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 = PrecCGSolver(A, b, delta, 1) +# should return x = [[1], [1], [1]] + +# 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) +# delta = 1.0e-6 +# x = PrecCGSolver(A, b, delta, 1) +# should return approx x = [[1], [0], [2], [0], [3]] + +# 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) +# delta = 1.0e-6 +# x = PrecCGSolver(A, b, delta, 1) +# should require MORE than n = 4 steps(because the matrix is not SPD)! + + +import numpy as np +import incompleteCholesky as IC +import LLTSolver as LLT + + +def matrnr(): + # set your matriculation number here + matrnr = 0 + return matrnr + + +def PrecCGSolver(A: np.array, b: np.array, delta=1.0e-6, verbose=0): + + if verbose: + print('Start PrecCGSolver...') + + countIter = 0 + + L = IC.incompleteCholesky(A) + x = LLT.LLTSolver(L, b) + + while MISSING STATEMENT: + MISSING CODE + countIter = countIter + 1 + if verbose: + print('STEP ', countIter, ': norm of residual is ', np.linalg.norm(r)) + + if verbose: + print('precCGSolver terminated after ', countIter, ' steps with norm of residual being ', np.linalg.norm(r)) + + return x \ No newline at end of file diff --git a/HW1/incompleteCholesky.py b/HW1/incompleteCholesky.py new file mode 100644 index 0000000..0795298 --- /dev/null +++ b/HW1/incompleteCholesky.py @@ -0,0 +1,86 @@ +# Optimization for Engineers - Dr.Johannes Hild +# Incomplete Cholesky decomposition + +# Purpose: incompleteCholesky finds lower triangle matrix L such that A - L * L ^ T is small, but +# eigenvalues are positive and sparsity is preserved + +# Input Definition: +# A: real valued symmetric matrix nxn +# alpha: nonnegative scalar, lower bound for eigenvalues of L * L ^ T.Default value: 1.0e-3. +# delta: scalar, if positive it is tolerance for recognizing nonsparse entry. +# If negative, do complete cholesky.Default value: 1.0e-6. +# verbose: bool, if set to true, verbose information is displayed + +# Output Definition: +# A: real valued lower triangle matrix nxn + +# Required files: +# < none > + +# Test cases: +# alpha = 0 +# delta = -1 +# verbose = true +# L = incompleteCholesky(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]]), alpha, delta, verbose) +# executes complete Cholesky decomposition with norm of residual approx. 8.89e-16 +# the warning is okay. + +# alpha = 1.0e-3 +# delta = 1.0e-6 +# verbose = true +# L = incompleteCholesky(np.array([4, 1, 0],[1, 4, 0],[0, 0, 4]]), alpha, delta, verbose) +# should return approximately +# L = [[2 0 0],[0.5 1.94 0], [ 0 0 2]] + +# alpha = 4 +# delta = 1.0e-6 +# verbose = true +# L = incompleteCholesky(np.array([4, 1, 0], [1, 4, 0], [0, 0, 4]]), alpha, delta, verbose) +# should return approximately +# L = [[2 0 0],[0.5 2 0], [ 0 0 2]] + +# alpha = 1.0e-3 +# delta = 1 +# verbose = true +# L = incompleteCholesky(np.array([[4, 1, 0], [1, 4, 0], [0, 0, 4]]), alpha, delta, verbose) +# should return approximately +# L = [[2 0 0],[0 2 0], [ 0 0 2]] + +import numpy as np + + +def matrnr(): + # set your matriculation number here + matrnr = 0 + return matrnr + + +def incompleteCholesky(A: np.array, alpha=1.0e-3, delta=1.0e-6, verbose=0): + L = np.copy(A) + dim = np.shape(L) + n = dim[0] + if n != dim[1]: + raise ValueError('A has wrong dimension.') + + if np.max(np.abs(A-A.T) > 1.0e-6): + raise ValueError('A is not symmetric.') + + if alpha < 0: + raise ValueError('range of alpha is wrong!') + + if delta < 0: + print('Warning: negative delta detected, sparsity is not preserved.') + + if verbose: + print('Start incompleteCholesky...') + + MISSING CODE + + if verbose: + residualmatrix = A - L @ L.T + residual = np.max(np.abs(residualmatrix)) + print('IncompleteCholesky terminated with norm of residual:') + print(residual) + + return L