diff --git a/Homework1/README.md b/Homework1/README.md index 0b3c0d4..03bd31b 100644 --- a/Homework1/README.md +++ b/Homework1/README.md @@ -1,145 +1,154 @@ # SP4E - Homework 1 ---- ## General Info This file provides a brief documentation and information related to the first Homework of the course "Scientific Programming for Engineers", fall 2019. This homework is done by **O. Ashtari** and **A. Sieber**. Last update: 16.10.2019 ---- ## Project Description This project is intended to solve a minimization problem on n-dimensional quadratic functions defined as S(X)= XᵀAX - Xᵀb where Xᵀ =[x1, x2, ..., xn], T represents the transpose operation, A is an n by n square matrix, and b a column vector of size n. ---- ## Requirements The project is created using Python 3.7. It moreover requires the following libraries to work properly: * numpy * argparse * scipy * matplotlib * sys * math ---- -## How to Use Program 1: optimizer.py +## Program 1: optimizer.py +Based on the problem description in the homework sheet, the objective function here is: + +S(X)= XᵀAX - Xᵀb + +Therefore, the user should note that definition of matrix A here differs from that in Program 2 in a factor 1/2. The goal here is using `scipy.optimize.minimize` routine to optimize the aforementioned function for a specific matrix A and vector b: + +A=[[4, 0], [1, 3]] , bᵀ=[0, 1] -This script solves minimization problems using the scipy.optimize.minimize routine. The minimization solver can be specified as an input argumeent by entering the following command line: +### How to Use +The minimization solver can be specified as an input argumeent by entering the following command line: ``` $ python3 optimizer.py method ``` Where the variable method is one of the following: * Nelder-Mead * Powell * CG (default) * BFGS * L-BFGS-B * TNC * SLSQP The initial guess of the minimization process has to be specified directly in the optimizer.py file by changing the value of the variable X0. The matrix A and vector B can also be modified directly in the file. ### conjugate_gradient.py ### post_processing.py Post-processing file that takes as inputs the value of A, B as well as the intermediate solutions of the iterative minimization process and the method used. The file generates a 3D plot displaying the function to be minimized as well as the intermediate solutions. ---- -## How to Use Program 2: conjugate_gradient.py +## Program 2: conjugate_gradient.py Based on the problem description in the homework sheet, the objective function here is: S(X)= ½(XᵀAX) - Xᵀb -Therefore, the user should note the factor 1/2 here, which means definition of matrix A here is different from that in Program 1. +Therefore, the user should note the factor 1/2 here, which means definition of matrix A here is different from that in Program 1. The goal here is implementing conjugate gradients algorithm such that takes A and b of any (compatible) dimension and finds its extremeum (i.e. gradient(S)=0). The problem is thus solving linear system ½(A+Aᵀ)X=b. +### How to Use To run the program, the objective function should be defined by getting matrix A and vector b. The code should be provided with a starting point x₀ as well. These three are mandatory arguments, followed after `-A`, `-b`, and `-x0` respectively which should be entered in rowmajor format separated by space. For example the command below runs the program for: -A=[[4, 1], [1, 3]] , bᵀ=[1,2] , x₀ᵀ=[2,1] +A=[[4, 1], [1, 3]] , bᵀ=[1, 2] , x₀ᵀ=[2, 1] ``` $ python3 conjugate_gradient.py -A 4 1 1 3 -b 1 2 -x0 2 1 ``` There are also some optional arguments for controlling the program, visualization, and verification: * `-Nmax`, followed by an integer number determines maximum number of iterations. By default `Nmax=20`. * `-criteria`, followed by a float number determines the convergence criteria to stop iterations. By default, `criteria=1e-6`. * `-detail` is a flag which tells the program show the result after each step. * `-plot` is a flag which tells the program show, graphically, evolution of the to-be-found optimum point on the contour of the objective function. This flag only works for 2D problems, i.e. S=S(x1, x2). Using this option, the program saves the 3D view as well as its 2D projection as .pdf files. * `-scipy` is a flag which asks the program to solve the same problem using `scipy`'s built-in function for conjugate gradient minimization and show the results for comparison. * `-h` is a flag which displays the help and list of arguments and flags. ### Example: The previous command line using optional flags and arguments will lead to: ``` $ python conjugate_gradient.py -A 4 1 1 3 -b 1 2 -x0 2 1 -Nmax 5 -criteria 1e-12 -detail -plot -scipy ======================Summary====================== Objective function: S = 0.5 * transpose(x) * A * x - transpose(x) * b A = [[4. 1.] [1. 3.]] b = [1. 2.] --------------------------------------------------- Step: 0 Residual: 8.54400374531753 Solution(x): [2. 1.] Objective: 7.5 --------------------------------------------------- Step: 1 Residual: 0.8001937042442401 Solution(x): [0.23564955 0.33836858] Objective: -0.5498489425981874 ==================Final Result===================== Number of steps: 2 Residual: 5.551115123125783e-17 Optimum x: [0.09090909 0.63636364] Optimum objective: -0.6818181818181817 =================================================== ==================Scipy Results===================== fun: -0.6818181818181795 jac: array([-9.68575478e-08, -4.47034836e-08]) message: 'Optimization terminated successfully.' nfev: 20 nit: 2 njev: 5 status: 0 success: True x: array([0.09090906, 0.63636363]) ==================================================== ``` By running the program, four separate blocks are displayed: * First, under "Summary", the objective function expression is shown to emphasize again on the notation used. Then `A` and `b` are shown to let the user check the inputs which define the problem. * Then, if the flag `-detail` is entered, status of the sought after point is shown at all steps which include iteration number, residual (L2-norm of the gradient vector), position of the point, and the corresponding value of the original function S(X). * After that, under "Final Result", number of iterations, final residual, optimum xᵀ, and the extremum value of the objective function are displayed. * Finally, if the flag `-scipy` is entered, the result from solving the same problem using conjugate gradient (CG) algorithm contained in `scipy.optimize.minimize` is also shown to let the user compare results from the written code with that given by the built-in solver to verify correctness of the present code ---- diff --git a/Homework1/conjugate_gradient.py b/Homework1/conjugate_gradient.py index 2ea53a7..8b94b77 100644 --- a/Homework1/conjugate_gradient.py +++ b/Homework1/conjugate_gradient.py @@ -1,122 +1,130 @@ +""" +Created in October 2019 +Authors: Omid Ashtari and Armand Sieber +Description: Scripts intended to minimize a quadratic function using conjugate gradient algorithm. + The function to be minimized is defined as S(X)= 0.5(X'AX) - X'B, with X = [x, y] and X' being its transpose + A is an nxn matrix and B nx1 vector +""" + import numpy as np import argparse import math import post_processing from scipy.optimize import minimize # Arguments parser my_parser = argparse.ArgumentParser() my_parser.add_argument('-A', action='store', type=float, nargs='+', required=True, help='The matrix A in rowmajor format: numbers separated by space') my_parser.add_argument('-b', action='store', type=float, nargs='+', required=True, help='The vector b in rowmajor format: numbers separated by space') my_parser.add_argument('-x0', action='store', type=float, nargs='+', required=True, help='The initial guess x0 in rowmajor format: numbers separated by space') my_parser.add_argument('-Nmax', action='store', type=int, default=20, help='Maximum number of iterations. Default: 20') my_parser.add_argument('-criteria', action='store', type=float, default=1e-6, help='Convergence criteria. Default: 1e-6') my_parser.add_argument('-detail', action='store_true', default=False, help='Showing details of every step during the convergence.') my_parser.add_argument('-plot', action='store_true', default=False, help='Showing and saving the evolution governed by CG.') my_parser.add_argument('-scipy', action='store_true', default=False, help='Using scipy function for minimization using CG method for comparison.') args = my_parser.parse_args() # Extracting matrix A, vector b, and the initial guess x0 from the input arguments A = np.array(args.A) b = np.array(args.b) x0 = np.array(args.x0) # Checking the compatibility of size of A with b (n: dimension of the problem) n=np.size(b) if np.size(A) != n*n: print('\n=====ERROR=====\n' 'Matrix A and vector b are not compatible in size.\n' '===============\n') exit() # Checking the compatibility of size of x0 with A and b if np.size(x0) != n: print('\n=====ERROR=====\n' 'Initial guess does not match to matrix A and vector b in size.\n' '===============\n') exit() # Converting matrix A from rowmajor to square form A=A.reshape((n,n)) # The to-be-optimized function def objective(X, *args): #objective function to be optimized AA = args[0] BB = args[1] S = 0.5*np.einsum('j,j',X,np.einsum('ij,j',AA,X))-np.einsum('j,j',X,BB) return S # Converting A to the generalized form: (A+transpose(A))/2 A_CG=0.5*(A+np.transpose(A)) # Converting optional settings from the arguments N_max=args.Nmax criteria=args.criteria # Iteration procedure step=0 x=x0 r=b-np.einsum('ij,j',A_CG,x0) norm_residual=math.sqrt(np.einsum('j,j',r,r)) p=r X_pp=x0; # The vertical stack for post-processing print('\n\n\n======================Summary======================') print('Objective function:\n'+ 'S = 0.5 * transpose(x) * A * x - transpose(x) * b\n\n'+ 'A =\n'+str(A)+'\n\n'+ 'b =\n'+str(b)) while norm_residual>criteria and step