# 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 ---- ## Program 1: optimizer.py Based on the first 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 the factor 1/2. The goal here is using the `scipy.optimize.minimize` routine to optimize the aforementioned function for a specific matrix A=[[4, 0], [1, 3]] and a specific vector bᵀ=[0, 1]: ### How to Use Basically, the program works with the simple following command line for a default method (conjugate gradients (CG)) and a default initial guess (x₀ᵀ=[3, 2]): ``` $ python optimizer.py ``` However, both method and initial guess can be entered by the user. To that end, `-method` should be followed by the name of the method of choice, and `-x0` should be followed by first and second elements of the initial guess, separated by space. An example is shown below: ``` $ python optimizer.py -method BFGS -x0 0 4 ``` Method can be each of the following items: * Nelder-Mead * Powell * CG (default) * BFGS * L-BFGS-B * TNC * SLSQP The program, finally, draws the evolution of the minimization process (red dots) on top of the contour plot of the objective function. The program saves the 3D view as well as its 2D projection as .pdf files. **N.B.**: some systems may require the user to type in `python3` instead of `python` in the command line. ### Example: The previous command line using optional flags and arguments will display the following information on the terminal: ``` $ python optimizer.py -method BFGS -x0 0 4 ======================Summary====================== Objective function: S = transpose(x) * A * x - transpose(x) * b A = [[4 0] [1 3]] b = [0 1] method: BFGS --------------------------------------------------- fun: -0.08510638297872285 hess_inv: array([[ 0.12765957, -0.0212766 ], [-0.0212766 , 0.17021277]]) jac: array([-1.02445483e-08, -1.86264515e-08]) message: 'Optimization terminated successfully.' nfev: 24 nit: 4 njev: 6 status: 0 success: True x: array([-0.0212766 , 0.17021276]) ``` Under "Summary", the objective function expression is shown to emphasize again on the notation used. Then `A`, `b`, and the chosen optimization method are displayed. Finally, details of the iterative optimization process are shown. Those are the default outputs of the `scipy.optimize.minimize` routine. ---- ## Program 2: conjugate_gradient.py Based on the second problem description in the homework sheet, the objective function here is: S(X)= ½(XᵀAX) - Xᵀb The user should note the factor 1/2. As such, the definition of matrix A here is different from the one of program 1. The goal of this program is to implement a conjugate gradients algorithm that takes A and b of any (compatible) dimension and finds its extremum (i.e. gradient(S)=0). The problem thus consists in solving the linear system ½(A+Aᵀ)X=b. ### How to Use To run the program, the objective function must first be defined through the matrix A and the vector b. The code must also be provided with a starting point x₀ (initial guess). These are mandatory arguments. They are entered in the program through the command line following their respective flags `-A`, `-b`, and `-x0`, in a rowmajor format separated by spaces. For instance, the command line below runs the program for: A=[[4, 1], [1, 3]], bᵀ=[1, 2], x₀ᵀ=[2, 1] ``` $ python conjugate_gradient.py -A 4 1 1 3 -b 1 2 -x0 2 1 ``` Additionally, optional arguments can be provided by the user to perform tasks such as controlling the program progression, visualization, and verification. Those additional functionalities are accessed by using the following flags: * `-Nmax`, followed by an integer number to determine the maximum number of iterations. By default `Nmax=20`. * `-criteria`, followed by a float number to determine 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 whether or not to use post-processing capabilities through the `post_processing.py` routine. 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. **N.B.**: some systems may require the user to type in `python3` instead of `python` in the command line. ### 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 defining the function to minimize. * Then, if the flag `-detail` is entered, the status of the sought solution is displayed for each iterative step. This status provides information concerning the iteration number, the residual (L2-norm of the gradient vector), the position of the current optimum, and the corresponding value of the original function S(X). * After that, under "Final Result", the total number of iterations, the final residual, the optimum xᵀ, and the extremum value of the objective function are displayed. * Finally, if the flag `-scipy` is entered, the same initial problem is computed using the conjugate gradient (CG) algorithm provided by `scipy.optimize.minimize` and the result is displayed. This is meant to provide a mean of comparison between our implementation of the conjugate gradient algorithm and the one provided by `scipy`. ---- ## Program 3: post_processing.py The function `plot_results` in this file gets the matrix A and the vector b as inputs (those are used to construct the surface S(x,y)) as well as the intermediate solutions of the iterative minimization process(to mark location of the to-be-found optimum point at each step). This function also gets the method used as a string to display it in the title of each plot. The function generates a 3D plot displaying the to-be-minimized function S as well as the intermediate solutions. This plot is saved as `3D_representation.pdf`. This function also saves the 2D projection of the discussed 3D plot on the x-y plane as `2D_projection.pdf`. There is no specific instruction for this program since it is used inside programs 1 and 2 and is not designed to be used independently.