R9482/Homework1ff175498c61emaster
README.md
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 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 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]):
$ python3 optimizer.py
However, both method and initial guess can be entered by the user. For this aim, -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. Method can be each of the following items:
- Nelder-Mead
- Powell
- CG (default)
- BFGS
- L-BFGS-B
- TNC
- SLSQP
The program, finally, draws evolution of the to-be-detected optimum point on the contour of the objective function. The program saves the 3D view as well as its 2D projection as .pdf files.
Example:
The previous command line using optional flags and arguments will lead to:
$ 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 optimization method which have defined the problem are displayed. After that, details of the iterative optization algorithm are shown.
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. 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]
$ 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
Program 3: post_processing.py
The function plot_results in this file gets the matrix A and the vector b as inputs (to construct the surfacae 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 and connect them). This function also gets the method used as a string to display as 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.