Page MenuHomec4science
Files F64150747

conjugate_gradient.py
No OneTemporary
Actions

File Metadata

Created
Fri, May 24, 22:28

conjugate_gradient.py
View Options

"""
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<N_max:
if args.detail:
print('---------------------------------------------------\n'
'Step: ' + str(step) + '\n'
'Residual: ' + str(norm_residual) + '\n'
'Solution(x): ' + str(x) + '\n'
'Objective: ' + str(objective(x,A,b)))
step = step+1
alpha=np.einsum('j,j',r,r)/np.einsum('j,j',p,np.einsum('ij,j',A_CG,p))
x=x+alpha*p
X_pp = np.vstack((X_pp,x))
r_backup=r
p_backup=p
r=r-alpha*np.einsum('ij,j',A_CG,p)
beta=np.einsum('j,j',r,r)/np.einsum('j,j',r_backup,r_backup)
p=r+beta*p
norm_residual=math.sqrt(np.einsum('j,j',r,r))
# Showing summary of calculations (last step)
print('\n==================Final Result=====================\n'
'Number of steps: ' + str(step) + '\n'
'Residual: ' + str(norm_residual) + '\n'
'Optimum x: ' + str(x) + '\n'
'Optimum objective: ' + str(objective(x,A,b)) + '\n'
'===================================================\n')
if args.scipy:
# Objective function to be minimized
def objective(X, *args):
A_sci = args[0]
B_sci = args[1].T
S_sci = np.dot(X.T, np.dot(A_sci, X)) - np.dot(X.T, B_sci)
return S_sci
# Minimization process using scipy
solution = minimize(objective, x0, args = (0.5*A, b), method = 'CG', options = {'maxiter':100})
print('==================Scipy Results=====================')
print(solution)
print('====================================================')
# Running post-processing module
if args.plot:
if n==2:
post_processing.plot_results(X_pp, 0.5*A, b, 'CG')
else:
print('Unable to plot!')

Event Timeline

c4science · Help