11 KiB
11 KiB
In [1]:
# DAOUST Alexandre 242015
In [2]:
import numpy as np
In [3]:
def fct(x, nargout=2):
x1, x2 = x
f = 2*x1**3 + 3*x2**3 + 2*x1**2 + 4*x2**2 + x1*x2 - 3*x1 - 2*x2
g = np.array([6*x1**2 + 4*x1 + x2 - 3,
6*x2**2 + 8*x2 + x1 - 2])
if nargout==3:
H = np.array([[12*x1 + 4, 1 ],
[ 1, 12*x2 + 8]])
return f, g, H
return f, g
In [4]:
# Test de f en 2,-1
x = np.array([2.0, -1.0])
f, g, H = fct(x, nargout=3)
print(f"f = {f:.4f}")
print(f"g = |{g[0]:.4f}|\n |{g[1]:.4f}|")
print(f"H = |{H[0][0]:.4f} {H[0][1]:.4f}|\n |{H[1][0]:.4f} {H[1][1]:.4f}|")
print()
# Test de f en 0,0
x = np.array([0.0, 0.0])
f, g, H = fct(x, nargout=3)
print(f"f = {f:.4f}")
print(f"g = |{g[0]:.4f}|\n |{g[1]:.4f}|")
print(f"H = |{H[0][0]:.4f} {H[0][1]:.4f}|\n |{H[1][0]:.4f} {H[1][1]:.4f}|")
# On a bien le bon résultat
In [5]:
def CoordinateDescent(x0, maxiter=10):
# Formules calculées à la main selon x1 et x2 respectivement
# Elles correspondent à f(xi) = Somme(coeff*xi^n) + C (C constante reprennant le(s) xj!=xi actuel(s) et
# les autres constantes) qu'on dérive et qu'on optimise.
x = x0.copy()
n = len(x)
for i in range(maxiter):
# Mise à jour de x1
delta = 4**2-4*6*(x[1]-3)
x[0] = (-4 + np.sqrt(delta))/12
#Mise à jour de x2
delta = 8**2 - 4*6*(x[0]-2)
x[1] = (-8 + np.sqrt(delta))/12
return x
In [6]:
# Test de CoordinateDescent
# On trouve bien un minimum car la gradient est nul (à une tolérance près) et la matrice hessienne est définie positive
# (on augmente f si on se déplace du point)
x = np.array([2.0, -1.0])
x_opt = CoordinateDescent(x)
print(f"x* = |{x_opt[0]:.4f}|\n |{x_opt[1]:.4f}|")
f, g, H = fct(x_opt, nargout=3)
print(f"f = {f:.4f}")
print(f"g = |{g[0]:.4f}|\n |{g[1]:.4f}|")
print(f"H = |{H[0][0]:.4f} {H[0][1]:.4f}|\n |{H[1][0]:.4f} {H[1][1]:.4f}|")
eigval, eigvec = np.linalg.eig(H)
print(f"Valeurs propres:\n {eigval[0]:.4f}\n {eigval[1]:.4f}")
In [7]:
def w1(fct, x, d, alpha, beta1):
f0, g0 = fct(x)
f1, g1 = fct(x + alpha*d)
return f1 <= f0 + alpha*beta1*(g0@d)
def w2(fct, x, d, alpha, beta2):
f0, g0 = fct(x)
f1, g1 = fct(x + alpha*d)
return g1@d >= beta2*(g0@d)
def wolfebissection(fct, x, d, alpha=1, beta1=0.0001, beta2=0.9):
alphal = 0
alphar = np.inf
wf1 = False
wf2 = False
# k pour ne pas boucler à l'infini
k = 0
while not (wf1 and wf2) and k < 1e4:
wf1 = w1(fct, x, d, alpha, beta1)
wf2 = w2(fct, x, d, alpha, beta2)
if not wf1:
alphar = alpha
alpha = (alphal + alphar)/2
elif not wf2:
alphal = alpha
if alphar < np.inf:
alpha = (alphal + alphar)/2
else:
alpha *= 2
k += 1
return alpha
In [8]:
# Test de wolfebissection
x = np.array([0.0, 0.0])
d = np.array([3.0, 2.0])
alpha = wolfebissection(fct, x, d)
print(f"alpha = {alpha:.4f}")
In [9]:
def MethGradient(x0, maxiter=10):
x = x0.copy()
for k in range(maxiter):
f, g = fct(x)
d = -g
alpha = wolfebissection(fct, x, d)
x = x + alpha*d
if np.linalg.norm(x) < 1e-9:
break
return x
In [10]:
# Test de MethGradient
# La méthode du Gradient se raproche de la solution (x1 est correct, x2 est légèrement différent) mais n'y converge pas à cause
# du calcul du pas alpha qui ne donne pas le pas exact pour arriver au minimum (autrement dit, on n'arrive pas à calculer
# le pas alpha permettant de converger vers le minimum). Au plus on est proche, au plus le calcul d'un alpha admissible prend
# du temps.
x = np.array([0.0, 0.0])
x_opt = MethGradient(x)
print(f"x* = |{x_opt[0]:.4f}|\n |{x_opt[1]:.4f}|")
f, g, H = fct(x_opt, nargout=3)
print(f"f = {f:.4f}")
print(f"g = |{g[0]:.4f}|\n |{g[1]:.4f}|")
print(f"H = |{H[0][0]:.4f} {H[0][1]:.4f}|\n |{H[1][0]:.4f} {H[1][1]:.4f}|")
eigval, eigvec = np.linalg.eig(H)
print(f"Valeurs propres:\n {eigval[0]:.4f}\n {eigval[1]:.4f}")
In [11]:
def MethNewton(x0, maxiter=10):
x = x0.copy()
for k in range(maxiter):
f, g, H = fct(x, nargout=3)
z = np.linalg.solve(H,g)
x = x - z
if np.linalg.norm(x) < 1e-9:
break
return x
In [12]:
# Test de MethNewton
# La méthode de Newton est très sensible au point de départ, cela peut conduire (comme c'est le cas ici) à trouver un point
# qui n'est pas un minimum mais un maximum.
x = np.array([-1.0, -1.0])
x_opt = MethNewton(x)
print(f"x* = |{x_opt[0]:.4f}|\n |{x_opt[1]:.4f}|")
f, g, H = fct(x_opt, nargout=3)
print(f"f = {f:.4f}")
print(f"g = |{g[0]:.4f}|\n |{g[1]:.4f}|")
print(f"H = |{H[0][0]:.4f} {H[0][1]:.4f}|\n |{H[1][0]:.4f} {H[1][1]:.4f}|")
eigval, eigvec = np.linalg.eig(H)
print(f"Valeurs propres:\n {eigval[0]:.4f}\n {eigval[1]:.4f}")
In [13]:
# Validation du point calculé par la CoordinateDescent (et par la MethGradient aussi)
# On trouve bien un minimum car la gradient est nul (à une tolérance près) et la matrice hessienne est définie positive
# (on augmente f si on se déplace du point)
x = np.array([0.0, 0.0])
x_opt = MethNewton(x)
print(f"x* = |{x_opt[0]:.4f}|\n |{x_opt[1]:.4f}|")
f, g, H = fct(x_opt, nargout=3)
print(f"f = {f:.4f}")
print(f"g = |{g[0]:.4f}|\n |{g[1]:.4f}|")
print(f"H = |{H[0][0]:.4f} {H[0][1]:.4f}|\n |{H[1][0]:.4f} {H[1][1]:.4f}|")
eigval, eigvec = np.linalg.eig(H)
print(f"Valeurs propres:\n {eigval[0]:.4f}\n {eigval[1]:.4f}")
In [ ]: