Optimización II

Existen criterios análiticos para encontrar máximos y mínimos, cuando se tienen una gran cantiad de variables, estos métodos son nada prácticos (¿se recuerda el número de párametros de la red de reconocimiento de dígitos?)
Se quiere encontrar el mínimo de una función , comenzando con un y moviendose en una dirección de busqueda controlando la longitud longitud de paso α (tasa de aprendizaje). Para ello se genera una sucesión
Contenido
Gradiente descendiente (steepest descent, gradient descent) Problema P9.1 Problema P9.2 Problema P9.3 Método de Newton Problema P9.4 Problema P9.5 Problema P9.6 Gradiente conjugado Problema P9.7 Problema P9.8 App Referencias

Gradiente descendiente (steepest descent, gradient descent)

Se elige la dirección de mayor descenso esto ocurre en la direción contraria al gradiente pues, es máxima cuando la dirección es paralela al gradiente y mínima cuando es "antiparalela", así el agoritmo es
En general la elección de la tasa de aprendizaje no es clara, pero si la función F es cuadrática se puede establecer una cota para la tasa de aprendizaje para producir trayectorias estables (los elementos de la sucesión no se mantieen oscilando) este valor esta acotado por el valor propio más grande de la matriz Hessiana de F
Otro critrerio para elegir la tasa de aprendizaje es optimizar respecto a la línea de busqueda, es decir
Observamos la sigueinte operación
Así el gradiente actual (se evalua en cda punto) es ortogonal a las direcciones de busqueda anteriores. Ahora esto es posible si la función es derivable.

Problema P9.1

P9_1.png
Solución
syms x1 x2
F(x1,x2) = 5 * x1^2 - 6 * x1 * x2 + 5 * x2^2 + 4*x1 + 4*x2;
figure
fcontour(F, [-7 7 -7 7])
alpha = 0.01;
x0 = [-1 -2.5]';
X = [x0];
GF(x1,x2) = gradient(F)
GF(x1, x2) = 
x1 = x0 - alpha * GF(x0(1),x0(2))
x1 = 
X = [X x1];
x2 = x1 - alpha * GF(x1(1),x1(2))
x2 = 
X = [X x2];
x3 = x2 - alpha * GF(x2(1),x2(2))
x3 = 
X = [X x3];
x4 = x3 - alpha * GF(x3(1),x3(2))
x4 = 
X = [X x4];
x5 = x4 - alpha * GF(x4(1),x4(2))
x5 = 
X = [X x5];
x6 = x5 - alpha * GF(x5(1),x5(2))
x6 = 
X = [X x6];
hold on
plot(X(1,:),X(2,:))
X = grades(F,x0,alpha,50);
figure
fcontour(F, [-7 7 -7 7])
hold on
plot(X(1,:),X(2,:))
vpa(X(1,end))
ans = 
vpa(X(2,end))
ans = 
X = [X;F(X(1,:),X(2,:))];
figure
fmesh(F,[-5 5 -5 5])
hold on
plot3(X(1,[1 end]),X(2,[1,end]),X(3,[1 end]))
HF = hessian(F)
HF(x1, x2) = 
vp = eig(HF)
vp(x1, x2) = 
digits(3)
vpa(2 / max(vp)) % alpha no debe superar
ans(x1, x2) = 
0.125
en particular tenemos solución análitica mínima (bvalore sporpiso positivos).
E = GF==[0 0]'
E(x1, x2) = 
% Compute analytic solution of a symbolic equation
E
E(x1, x2) = 
solution3 = solve(E,[sym('x1'),sym('x2')],'IgnoreProperties',true);
% Display symbolic solution returned by solve
displaySymSolution(solution3);
solution3 = struct with fields:
x1: [1×1 sym] x2: [1×1 sym]
where

Problema P9.2

P9_2.png
Solución
x0 = [0 -2]';
p0 =-GF(x0(1),x0(2))
p0 = 
g0 = GF(x0(1),x0(2))
g0 = 
A0 = HF(x0(1),x0(2));
a0 = - (g0'*p0)/(p0'*A0*p0)
a0 = 
x1 = x0 + a0*p0
x1 = 
g1 = GF(x1(1),x1(2))
g1 = 
Esto es porque la condición inicial esta apuntando en la dirección del eigenvector de

Problema P9.3

P9_31.png
P9_32.png
Solución
syms w b
x = [w b]';
G = [2 1;-1 1];
t = [0.5 0]';
F(w,b)= (t-G*x)'*(t-G*x)
F(w, b) = 
x0 = [1 1]';
alpha = 0.05;
GF=gradient(F)
GF(w, b) = 
g0 = GF(x0(1),x0(2))
g0 = 
x1 = x0 - alpha * g0
x1 = 
syms x1 x2
X = grades(F(x1,x2),x0,alpha,50);
figure
fcontour(F, [-3 3 -3 3])
hold on
plot(X(1,:),X(2,:))
HF = hessian(F)
HF(w, b) = 
vpa(2/max(eig(HF))) % es un valor de tasa válido
ans(w, b) = 
0.189

Método de Newton

Suponiendo diferenciabilidad de segundo orden, e invertibilidad de la matriz Hessiana.
derivando respecto a la variación
Así el algoritmo de Newton
Puede converger más rápido que gradiente descendiente, pero puede converger a puntos sillas, pues converge a puntos estacionarios, es decir no distingue el tipo de punto salvo que es estacionario. Es posible que exista divergencia, mientras que el gradiente descendiente garantiza convergencia si la tasa de aprendizaje no es muy grande (optimizando en cada etapa podría decidirse como elegir la tasa de aprendizaje). Newton computacinamente es mas exigente (Hesianas e inversas)

Problema P9.4

P9_4.png
Solución
syms x1 x2
x0 = [1 -2]';
F(x1,x2) = exp(x1^2 - x1 + 2 * x2^2 + 4);
figure
fmesh(F,[-1 1 -1 1])
GF = gradient(F)
GF(x1, x2) = 
HF = hessian(F)
HF(x1, x2) = 
g0 = GF(x0(1),x0(2))
g0 = 
A0 = HF(x0(1),x0(2))
A0 = 
x11 = x0 - inv(A0)*g0; % puden proponerse variaciones para calcular la inversa
vpa(x11)
ans = 
La función es una composición con una función monotona (la exponencial) así el mínimo del exponente sera el mínimo del exponente. La iteración avanza poco hacia el mínimo porque partimos de una aporximación de Tayoor y la función no tienen una buena aproximación de segundo orden cerca del incio.
e(x1,x2) = x1^2 - x1 + 2 * x2^2 + 4;
[x1a x2a]=solve(gradient(e)==[0 0]')
x1a = 
x2a = 
0

Problema P9.5

P9_5.png
Solución
syms x1 x2
A = [1 -1;-1 1];
x = [x1 x2]';
x0 = [1 0]';
F(x1,x2) = (1/2)*x' * A * x
F(x1, x2) = 
GF = hessian(F)
GF(x1, x2) = 
La matriz hessiana no es ivertible (no tienen rango completo) así Newton no fuciona
X = grades(F,x0,0.1,50);
figure
fcontour(F, [-3 3 -3 3])
hold on
plot(X(1,:),X(2,:))
vpa(X(:,end))
ans = 
gradient(F)
ans(x1, x2) = 
a = solve(gradient(F)==[0 0]')
Warning: Unable to solve symbolically. Returning a numeric solution using vpasolve.
a = struct with fields:
x1: [1×1 sym] x2: [1×1 sym]
a.x2
ans = 
0
a.x1
ans = 
0
figure
fmesh(F, [-3 3 -3 3])

Problema P9.6

P9_6.png
Solución
syms x1 x2
F(x1,x2) = x1^3 + x1*x2 -x1^2*x2^2;
x0 = [1 1]';
GF =gradient(F);
HF = hessian(F);
g0 = GF(x0(1),x0(2));
A0 = HF(x0(1),x0(2));
x1 = x0 -inv(A0)*g0 ;
vpa(x1)
ans = 
P2 = F(x0(1),x0(2))+g0'*(x-x0) +(1/2)*(x-x0)'*A0*(x-x0)
P2 = 
[xs1,xs2] = solve(gradient(P2)==[0 0]')
Warning: Unable to solve symbolically. Returning a numeric solution using vpasolve.
xs1 = 
0.588
xs2 = 
1.12
Newton converge a un punto estacionario. Pero este punto es un punto silla
vpa(eig(A0))
ans = 

Gradiente conjugado

Se busca una variante de Newton que pueda encontrar mínimos de forma exacta para formas cuadraticas. Pero que tenga menos costo computacional (no calcular segundas derivadas).
Conjugados. Un conjunto de vectores es mutuamente ortogonal respecto a una matriz definida positiva A si (A mapea a al espaci ortogonal de )
Se sabe que si se busca el mínimo en n (dimensión de A) direcciones conjugadas se encontrará el mínimo en n pasos a lo más (siempre y cuando este exista en una forma cuadrática). La idea es construir las direcciones conjugadas pero sin la necesidad de calcular segundas derivadas.
(gradiente de una forma cuadratica)
(Hessiano de una forma cuadrática)
(una relación entre gradiente y Hessiano)
(El incremeto de busqueda)
(condición de conjugado, ya no es necesaria A)
(dirección de busqueda inicial)
Se necesita que que cada dirección de busqueda sea ortogoal a los cambio del gradiente . En un proceso similar a Gram-Schmidt se tiene que
Donde se tienen diferentes propuestas para las coeficientes
(Hestenes y Stiefel)
(Fletcher y Reeves)
(Polak y Ribiére)
Dado , F, ,
  1. , con alguna elección de . Para la primer iteración
  2. o constante (de aquí pude haber variaciones)

Problema P9.7

P9_7.png
Solución
syms w b
x = [w b]';
G = [2 1;-1 1];
A = [10 2;2 4];
t = [0.5 0]';
F(w,b)= (t-G*x)'*(t-G*x)
F(w, b) = 
x0 = [1 1]';
GF=gradient(F)
GF(w, b) = 
% primer iteración
g0 = GF(x0(1),x0(2));
p0 = -g0;
alpha0 = - g0'*p0 / (p0'*A*p0);
x1 = x0 + alpha0 * p0;
digits(3)
vpa(x1)
ans = 
% iteración 2
g1 = GF(x1(1),x1(2));
beta1 = (g1-g0)'* g1 / (g0' * g0); % Polak
p1 = -g1 + beta1 * p0;
alpha1 = - g1'*p1 / (p1'*A*p1);
x2 = x1 + alpha1 * p1;
digits(3)
vpa(x2)
ans = 
X = [x0 x1 x2];
figure
fcontour(F, [-2 2 -2 2])
hold on
plot(X(1,:),X(2,:))

Problema P9.8

P9_8.png
Conjunto conjugado
combinación lineal
Necesariamente el coeficiente es cero porque A es definida positiva (existencia de mínimo)

App

Steepest Descent for Quadratic
nnd9sdq
Comparison of Methods
nnd9mc
Newton's Method
nnd9nm
Steepest Descent
nnd9sd

Referencias

El material se toma del libro de Martin Hagan et. al. enlace
function X = grades(F,x0,alpha,n)
syms x1 x2
X = [x0];
GF(x1,x2) = gradient(F);
for i = 1:n
x = x0 - alpha * GF(x0(1),x0(2));
X = [X x];
x0 = x;
end
end