Commit 0c5d2bc5 authored by Manon Eugénie Voisin--Leprince's avatar Manon Eugénie Voisin--Leprince
Browse files

correction week4

parent 8853ae1e
#!/usr/bin/env python3
import numpy as np
def solve_GMRES(A, b, x0, max_iter=30, tol=1e-15, callback=None):
# Initialisation
r = b - np.einsum('ij, j->i', A, x0)
h = np.zeros((max_iter+1, max_iter))
q = np.array([0 for i in range(max_iter)], dtype=object)
q[0] = r/np.linalg.norm(r)
n = min(max_iter, len(b))
for k in range(n):
# Arnoldi iteration
Q, h = arnoldi_iteration(A, q, h, k, max_iter)
# least square solution
x, sol, dependent_variable = least_square_solution(
r, max_iter, h, q, x0, callback=callback)
# Calcul of the residual
residual = np.linalg.norm(np.einsum('ij, j->i', h, sol)-np.linalg.norm(dependent_variable))/np.linalg.norm(dependent_variable)
if residual < tol:
return x, 0
return x, 1
def arnoldi_iteration(A, q, h, k, max_iter):
r_new = np.einsum('ij, j -> i', A, q[k])
for _l in range(k+1):
h[_l, k] = np.einsum('i, i', q[_l], r_new)
r_new = r_new - q[_l]*h[_l, k]
h[k+1, k] = np.linalg.norm(r_new)
if (h[k + 1, k] != 0 and k != max_iter - 1):
q[k + 1] = r_new/h[k + 1, k]
return q, h
def least_square_solution(r, max_iter, h, q, x0, callback=None):
dependent_variable = np.zeros(max_iter + 1)
dependent_variable[0] = np.linalg.norm(r)
sol = np.linalg.lstsq(h, dependent_variable, rcond=None)[0]
x = np.dot(np.transpose(q), sol) + x0
if callback is not None:
callback(x)
return x, sol, dependent_variable
# Instructions to launch the program
- In order to launch the program, simply launch:
```python
python3 ./main.py a_11 a_12 a_21 a_22 b_1 b_2
```
With the coefficients $a_ij$ and $b_i$ the terms of the matrix $A$ and vector $b$ in the quadratic form $1/2 x^T A x - x^T b$, as advertised by the help:
```python
python3 ./main.py -h
```
This will compute the minimum of the quadratic function using one optimizer.
- In order to change the method you can add an additional argument:
```python
python3 ./main.py a_11 a_12 a_21 a_22 b_1 b_2 -m BFGS
python3 ./main.py a_11 a_12 a_21 a_22 b_1 b_2 -m GMRES
python3 ./main.py a_11 a_12 a_21 a_22 b_1 b_2 -m my_GMRES
```
- In order to make a plot, use the following option:
```python
python3 ./main.py a_11 a_12 a_21 a_22 b_1 b_2 -p 1
```
Example of launch:
```python
python3 ./main.py -m my_GMRES 8 1 1 3 2 4 -p 1
```
\ No newline at end of file
#!/usr/bin/env python3
from optimizer import optimize
import numpy as np
import argparse
################################################################
# main function: parsing of program arguments
################################################################
if __name__ == '__main__':
parser = argparse.ArgumentParser(
description='Optimization exercise in python')
parser.add_argument(
'coeffs', type=float,
nargs='+',
help=('Specify the coefficients of the quadratic function,'
' in the following order: a_11, a_12, a_21, a_22, b_1, b_2.'
' The coefficients are the terms of the matrix A'
' and the vector b in the quadratic form 1/2 x^T A x - x^T b.'))
parser.add_argument(
'-m', '--method', type=str,
help=('Specify the method to use for the optimization.'
' Can be BFGS/GMRES/my_GMRES'), default='BFGS')
parser.add_argument(
'-p', '--plot', type=bool,
help=('Choose to plot or not the solution.'
' True/False'), default=False)
args = parser.parse_args()
A = np.array(args.coeffs[:4]).reshape(2, 2)
b = np.array(args.coeffs[4:])
# transform to dictionary
kwargs = vars(args)
del kwargs['coeffs']
# definitions of the quadratic form
def func(X, A, b):
"functor calculation of the quadratic form"
X = X.reshape(X.flatten().shape[0]//2, 2)
# computes the quadratic form with einsum
res = 0.5*np.einsum('ai,ij,aj->a', X, A, X)
res -= np.einsum('i,ai->a', b, X)
return res
optimize(A, b, func=lambda x: func(x, A, b),
jac=lambda x: A@x-b,
hess=lambda x: A,
**kwargs
)
#!/usr/bin/env python3
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import minimize as scipySolve
from scipy.sparse.linalg import lgmres as gmresSolve
from GMRES import solve_GMRES
def plotFunction(func, iterations, title):
'Plotting function for plane and iterations'
x = np.linspace(-3, 3, 100)
y = np.linspace(-3, 3, 100)
x, y = np.meshgrid(x, y)
X = np.stack([x.flatten(), y.flatten()], axis=1)
z = func(X).reshape((100, 100))
xy = np.array(iterations)
zs = func(xy)
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(x, y, z, linewidth=0, antialiased=True,
cmap=plt.cm.viridis, alpha=0.2)
ax.contour(x, y, z, 10, colors="k", linestyles="solid")
ax.plot(xy[:, 0], xy[:, 1], zs, 'ro--')
ax.grid(False)
ax.view_init(elev=62., azim=143)
ax.set_title(title)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.show()
################################################################
# generic optimization program
################################################################
def optimize(A, b, method='BFGS', plot=False, func=None, **kwargs):
if func is None:
raise RuntimeError('Need a function to optimize')
tol = 1e-6
x0 = np.array([0, 0])
if method == 'BFGS':
optimizer = scipySolve
# avoids a warning
if 'hess' in kwargs:
del kwargs['hess']
elif method == 'GMRES':
optimizer = gmresSolve
elif method == 'my_GMRES':
optimizer = solve_GMRES
else:
raise RuntimeError('unknown solve method: ' + str(method))
if method == 'BFGS':
iterations = []
iterations.append(x0)
res = optimizer(func, x0, tol=tol,
callback=lambda Xi: iterations.append(Xi), **kwargs)
if method == 'GMRES' or method == 'my_GMRES':
iterations = []
iterations.append(x0)
res, info = optimizer(A, b, x0, tol=tol,
callback=lambda Xi: iterations.append(Xi))
print('Method: ', method, ' Solution :', res)
if plot is True:
plotFunction(func, iterations, method)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment