HW1 uploaded

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

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

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