diff --git a/README.md b/README.md index b4d7adb0..6452df24 100644 --- a/README.md +++ b/README.md @@ -1,33 +1,33 @@ # SP4E Homeworks Students: * Lars B. * Bertil T. In this repository, we will keep all homeworks for SP4E ## Homework 1 -See folder named *hw1-conjugate-gradient*. - +See folder named *hw1-conjugate-gradient*. + #### Requirements * Numpy * Scipy * Matplotlib * Argparse The code as been tested with Python 3. #### Usage Run the *main_minimize.py* with the following arguments: * -A < matrix elements in row major order > * -b < vector elements > * -x0 < initial guess elements > * -method < CG-scipy or CG-ours > * -plot < True or False > #### Example The following command will run the program with the coefficients given in Exercise 1: -python main_minimize.py -A 8 0 2 6 -b 0 1 -x0 4 4 -method CG-scipy -plot True +python3 main_minimize.py -A 8 0 2 6 -b 0 1 -x0 4 4 -method CG-scipy -plot True -NB: Note the scaling of S(x) with 1/2 between Ex 1 and 2. +NB: Note that A must be s.p.d. in order for the conjugate gradient method to correctly solve Ax=b. Thus, the matrix A in Exercise 1 should not be used to compare the scipy version (CG-scipy) againts our version (CG-ours). diff --git a/hw1-conjugate-gradient/conjugate_gradient.py b/hw1-conjugate-gradient/conjugate_gradient.py index dc1547b7..b2642986 100644 --- a/hw1-conjugate-gradient/conjugate_gradient.py +++ b/hw1-conjugate-gradient/conjugate_gradient.py @@ -1,54 +1,59 @@ import numpy as np def conjgrad(A, b, x, tolerance, max_iterations): # Dimension of problem dim = len(b) # Convert input to column vectors b.reshape(dim, -1) x.reshape(dim, -1) # Save steps taken my minimizer in this array (can be used for plotting) steps = x.reshape(-1, dim) + # Precompute tolerance squared + tolerance_sq = tolerance*tolerance + # Initialization - #k = 0 + k = 0 r = b - np.einsum('ij,j', A, x) # b - A @ x - p = r if (np.sqrt(np.sum(r**2)) < tolerance): return (x, steps) + p = r + r_sq = np.einsum('j,j',r,r) + # CG Algorithm for i in range(0, max_iterations): - alpha = np.einsum('j,j',r,r) / np.einsum('j,j',p,np.einsum('ij,j', A, p)) - x = x + alpha*p - steps = np.vstack((steps, x.reshape(-1, dim))) # keep track of steps for plotting - r_next = r - alpha*np.einsum('ij,j', A, p) - beta = np.einsum('j,j',r_next,r_next)/np.einsum('j,j',r,r) - r = r_next - p = r + beta*p - #k = k + 1 - #print(k) - if (np.sqrt(np.sum(r**2)) < tolerance): + Ap = np.einsum('ij,j', A, p) + alpha = r_sq / np.einsum('j,j', p, Ap) + x = x + alpha * p + steps = np.vstack((steps, x.reshape(-1, dim))) # keep track of steps for plotting + r = r - alpha * Ap + r_sq_next = np.einsum('j,j', r, r) + if (r_sq_next < tolerance_sq): break - + p = r + (r_sq_next / r_sq) * p + r_sq = r_sq_next + k = k + 1 + print(k) return (x, steps) # want to return x-array that minimizes S(x) ########## TESTING AND VALIDATING CONJGRAD FUNCTION ########### #max_iterations = 1000 #tolerance = 0.001 ## Input parameters for 3D: #A = np.array([[2.0, 3.0, 5.0], [2.0, 3.0 , 6.0], [3.0, 6.0, 9.0]]) #b = np.array([2.0, 3.0, 5.0]) #x = np.array([10.0, 20.5, 0.0]).T # initial guess ## Input parameters for 2D: #A = np.array([[4.0, 0], [1.0, 3.0]]) #b = np.array([0.0, 1.0]) #x = np.array([4, 4]) # initial guess #result, steps = conjgrad(A, b, x, tolerance, max_iterations) #print("Minimum found at position ", result, " after ", len(steps)-1, " iterations.") diff --git a/hw1-conjugate-gradient/main_minimize.py b/hw1-conjugate-gradient/main_minimize.py index a628c335..ed03624b 100644 --- a/hw1-conjugate-gradient/main_minimize.py +++ b/hw1-conjugate-gradient/main_minimize.py @@ -1,66 +1,73 @@ ### IMPORT LIBRARIES import numpy as np import scipy as sp import argparse import optimizer from conjugate_gradient import conjgrad from quadratic_function import quadratic_function from plotting import plotting ### DEFINE PARAMETERS / COEFFICIENTS -max_iterations = 1000 -tolerance = 0.001 +max_iterations = 1000000 +tolerance = 1e-3 ### FUNCTION MAP FOR OPTIMIZER METHOD function_map = {'CG-ours' : conjgrad, 'CG-scipy' : optimizer.find_minimum} ### ARGPARSE ### Sets matrix coefficients + optimization Method parser = argparse.ArgumentParser(description='Process Minimization Schemes') parser.add_argument('-A', action='append', nargs='+', help='Input matrix (row major format) : -A matrix') parser.add_argument('-b', action='append', nargs='+', help='Input vector: -b vector') parser.add_argument('-x0', action='append', nargs='+', help='Initial guess: -x0 vector') parser.add_argument('-method', choices=function_map.keys(), help='Choose method of minimization (CG-ours or CG-scipy)') parser.add_argument('-plot', choices=["True", "False"], help='Plot result: True or False') args = parser.parse_args() A = np.array([float(k) for k in args.A[0]]) b = np.array([float(k) for k in args.b[0]]) x = np.array([float(k) for k in args.x0[0]]) x0 = np.copy(x) # saved for later use in defining plotting bounds dim = len(b) A = np.reshape(A, (dim, dim)) b.reshape(dim, -1) x.reshape(dim, -1) ### PRINT INITIAL STATE print("Matrix A: \n", A) print("Vector b: \n", b) print("Initial guess x0: \n", x) print("Tolerance: ", tolerance) print("Max Iterations: ", max_iterations) ### RUN THE OPTIMIZATION func = function_map[args.method] if args.method == 'CG-ours': print("You chose method: ", args.method) result, steps = func(A, b, x, tolerance, max_iterations) elif args.method == 'CG-scipy': print("You chose method: ", args.method) optimizer.STEPS = x optimizer.DIM = dim result = func(quadratic_function, x, tolerance, max_iterations, (A,b)) steps = optimizer.STEPS + +### EXACT SOLUTION +exact = np.linalg.solve(A, b) + ### PRINT RESULTS -print("The (approximated) minimum value is found at the coordinate ", result) -print("The (approximated) minimum value is ", quadratic_function(result, A, b)) -print("These are the steps taken by the minimizer: ") +print("The (approximated) solution of the LSE is: ", result) +print("The exact solution (by solving the LSE) is: ", exact) +print("The norm of the residual is: ", np.sqrt(np.dot((result-exact),(result-exact)))) +print("The (approximated) min value of the quadratic function is ", quadratic_function(result, A, b)) +print("Number of steps performed in the iterative algorithm: ", steps.shape[0]-1) +print("These are the steps taken: ") print(steps) ### PLOT RESULTS if (args.plot == "True"): plotting(dim, steps, x0, A, b) diff --git a/hw1-conjugate-gradient/optimizer.py b/hw1-conjugate-gradient/optimizer.py index 2afcdef2..68e32a8f 100644 --- a/hw1-conjugate-gradient/optimizer.py +++ b/hw1-conjugate-gradient/optimizer.py @@ -1,27 +1,27 @@ # Import libraries import numpy as np import scipy.optimize ### Global variables (only necessary for plotting of steps taken by optimizer) # Dimension of problem, default value = 2 will be overwritten by main_minimize DIM = 2 # Initial guess to minimizer, default value = (4,4) will be overwritten by main_minimize STEPS = np.array([4, 4]) # This array will contain the steps of the minimizer on its way to the minimum STEPS = STEPS.reshape(-1, DIM) ### This function fills the STEPS array, and is used as input to scipy.optimize.minimize def my_callback(xk): global STEPS STEPS = np.vstack((STEPS, np.array(xk))) #### This function returns the minimum value of the function func -# x: initial guess of the minimum +# x: initial guess of the minimum # tolerance: termination criteria for algortihm # max_iterations: maximum number of iterations as termination criteria for algortihm # extra_args: tuple containing matrix A and b to be used if not default A and b used def find_minimum( func, x, tolerance, max_iterations, extra_args=() ): result = scipy.optimize.minimize(func, x0 = x, args=extra_args, method='CG', tol=tolerance, options={ "maxiter" : max_iterations}, callback=my_callback) return result.x