Optimización¶
import numpy as np
import matplotlib.pyplot as plt
El método de la sección áurea¶
Vamos a optimizar una función $g$ de una sola variable. El método de la sección áurea es un método de intervalo.
Si el intervalo contiene un único mínimo en $[a,b]$, tendremos dos puntos interiores $x_1$ y $x_2$ de forma que $a \lt x_1 \lt x_2 \lt b.$ Teniendo en cuenta los valores de la función $g$ en $x_1$ y $x_2$, podemos descartar $[a,x_1)$ o $(x_2,b]$ manteniendo el mínimo dentro del intervalo.
Escogemos los puntos de forma que estén a una fracción $\phi$ y $1-\phi$ de uno de los extremos, con $\phi=\dfrac{\sqrt{5}-1}{2}$. $\phi$ se escoge de forma que si descartamos $(x_2,b]$ el siguiente $x_2$ sea igual que el antiguo $x_1$. Y si descartamos $[a,x_1)$ el siguiente $x_1$ sea el antiguo $x_2.$ Así, en cada iteración, solo tendremos que evaluar la función en un nuevo punto.
Al ser un método de intervalo, la longitud del intervalo final es una cota del error.
aurea
Algoritmo¶
- Sea $a_1=a$, $b_1=b$ y $\phi = \dfrac{\sqrt{5}-1}{2}$
- Para $k=1,2,\ldots,\mathrm{MaxNumIter}$
- Calcular los puntos
- $x_1=a+(1-\phi)\,(b-a)$
- $x_2=a+\phi\,(b-a)$
- si $g(x_1) \gt g(x_2)$ entonces:
$ a_{k+1}=x_1,$ $ b_{k+1}=b_k,$
- en otro caso:
$ a_{k+1}=a_{k},$ $ b_{k+1}=x_2.$
- Si se satisface el criterio de parada, parar.
Escribir una función seccionaurea(g,a,b,tol=1.e-6,maxiter=100) que tenga como argumento de entrada la función lambda g, los extremos del intervalo a, b, la cota del error absoluto tol y el número máximo de iteraciones maxiter, y como argumentos de salida el mínimo aproximado de la función y el número de iteraciones realizadas. Tomar como mínimo m el punto medio del último intervalo. Dibujar la función en el intervalo $[a,b]$ y el punto m sobre la curva.
- Calcular el mínimo de la función $g_1(x) = (2x-1.1)^2 + (3x-3.3)^2$ contenido en el intervalo $[0,3]$
- Sin modificar la función, utilizar $g_2(x) = x^5 - 3 x^2 + 1.6$ con el intervalo inicial $[0,2]$
Nota: La cota del error absoluto es la longitud del intervalo.
%run Ejercicio1.py
Una función de python para calcular el mínimo de una función es minimize:
import scipy.optimize as op
m1 = op.minimize(g1, x0=0)
print(m1.x[0])
m2 = op.minimize(g2, x0=1)
print(m2.x[0])
Representación de una función de dos variables¶
Por ejemplo, queremos representar la función $f(x,y)=x^2+y^2$ en el rectángulo $[-5,0]\times[-2,2]$
Definimos la función lambda de dos variables
f = lambda x,y: x**2 + y**2
Definimos una malla de puntos
x = np.linspace(-5,0,6) # En [-5,0] 6 puntos equiespaciados
y = np.linspace(-2,2,5) # En [-2,2] 5 puntos equiespaciados
X, Y = np.meshgrid(x,y)
Y hemos creado dos matrices
print(X)
print(Y)
Y para cada par de puntos $(x,y)$ de la malla (los que están situados en la misma posición de las matrices $X$ e $Y$) tendremos un valor de $f(x,y)$
Z = f(X,Y)
print(Z)
Podemos representar la función con curvas de nivel
plt.figure()
plt.contour(X,Y,Z,cmap='jet')
plt.colorbar()
plt.show()
plt.figure()
plt.contourf(X,Y,Z,cmap='jet')
plt.colorbar()
plt.show()
O como una superficie
from mpl_toolkits.mplot3d import Axes3D
fig1 = plt.figure(figsize=(10,5))
ax1 = plt.axes(projection ='3d')
ax1.plot_surface(X, Y, Z, cmap='jet')
plt.show()
fig2 = plt.figure(figsize=(10,5))
ax2 = plt.axes(projection ='3d')
ax2.plot_wireframe(X, Y, Z)
plt.show()
Si construimos una malla más densa tendremos una representación más suave
x = np.linspace(-5,0)
y = np.linspace(-2,2)
X, Y = np.meshgrid(x,y)
Z = f(X,Y)
plt.figure()
plt.contourf(X,Y,Z,cmap='jet')
plt.colorbar()
plt.show()
El método del gradiente con tasa de aprendizaje¶
Buscamos el mínimo de una función. Pensemos en una función de dos variables y por lo tanto, podemos pensar que estamos sobre la superficie que representa la función
- Nos situamos en un punto inicial $(x_0,y_0).$
- Buscamos la dirección del gradiente.
- Nos orientamos en sentido contrario al gradiente, puesto que el gradiente se orienta en sentido creciente de la función, hacia arriba, y nosostros buscamos un mínimo y queremos ir hacia abajo.
- Avanzamos en línea recta en la dirección contraria al gradiente con la fórmula$$ \left(\begin{array}{c} x_{1}\\ y_{1} \end{array}\right)=\left(\begin{array}{c} x_{0}\\ y_{0} \end{array}\right)-\eta\,\left(\begin{array}{c} f_x(x_0,y_0)\\ f_y(x_0,y_0) \end{array}\right) $$donde $\eta$ es una constante que tendremos que fijar. Si $\eta$ es pequeña, la sucesión convergerá lentamente. Si es demasiado grande la sucesión no convergerá.
- Buscamos de nuevo la dirección del gradiente y repetimos los pasos anteriores hasta que se cumpla una condición de parada.
- Como en el mínimo $\nabla f = \mathbf{0}$ (condición necesaria de mínimo), una condición de parada es que $\left\Vert \nabla f \right\Vert$ sea menor que una cierta tolerancia.
Escribir una función gradiente(f,eta,fx,fy,x0,y0,tol=1.e-6,maxiter=100) que tenga como argumentos de entrada la función de dos variables lambda f, la tasa de aprendizaje eta, las dos componentes del vector gradiente fx y fy de f, definidas también como funciones lambda, el punto inicial x0, y0, la tolerancia tol y el número máximo de iteraciones y como argumentos de salida un vector (array numpy) x y otro y que contengan las coordenadas $x$ e $y$ de las sucesivas iteraciones y el número de iteraciones realizadas.
Calcular el mínimo de la función $f_1(x,y) = x^2 + 3y^2$, utilizando
- Punto inicial $(2,1)$
- Tasa de aprendizaje $\eta = 0.1$,
- Tolerancia
tol = 0.001tal que $\left\Vert \nabla f \right\Vert\lt tol$
Calcular el mínimo de la función $f_2(x,y) = x^2 + y^2 - x\,y - 3\,x$, ajustando a tanteo $\eta$ y utilizando
- Punto inicial $(-2,2)$
- Tolerancia
tol = 0.001tal que $\left\Vert \nabla f \right\Vert\lt tol$ - Realizando menos de 15 iteraciones.
En ambos casos dibujar la función y la sucesión de puntos.
%run Ejercicio2.py
Aquí también podemos usar la función minimize cambiando la sintaxis de la función lambda
import scipy.optimize as op
f1 = lambda x: x[0]**2 + 3*x[1]**2
m1 = op.minimize(f1, x0 = np.array([2,1]))
print(m1.x)
f2 = lambda x: x[0]**2 + x[1]**2 - x[0]*x[1] - 3*x[0]
m1 = op.minimize(f2, x0 = np.array([-2,2]))
print(m1.x)
El método del descenso más pronunciado¶
El método es similar a anterior pero se optimiza el paso de forma que llegué hasta el punto de su trayectoria donde el valor de la función es mínimo.
Realizaremos la primera iteración por el método del gradiente usando la fórmula
$$ \left(\begin{array}{c} x_{1}\\ y_{1} \end{array}\right)=\left(\begin{array}{c} x_{0}\\ y_{0} \end{array}\right)-h\,\left(\begin{array}{c} f_x(x_0,y_0)\\ f_y(x_0,y_0) \end{array}\right)\qquad $$Cuando nos movemos en línea recta sobre la superficie estamos recorriendo una curva $g$
$$ g(h)=f\left(x_{0}-h\,f_{x}(x_0,y_0),y_{0}-h\,f_{y}(x_0,y_0)\right) $$y buscaremos el paso $h$ que hace que recorramos la curva hasta su valor mínimo. En el punto del que partimos $(x_0,y_0)$ tenemos que $h=0.$ Además $h$ ha de ser positiva porque queremos movernos en sentido contrario al gradiente y esto ya está determinado por el signo menos delante de $h.$ Por lo tanto, si buscamos con un método de intervalo como el método de la sección áurea, buscaremos el mínimo en un intervalo de la forma $[0,b].$
Escribir una función descenso(f,fx,fy,x0,y0,tol=1.e-6,maxiter=100) que tenga como argumentos de entrada la función lambda de dos variables f, las dos componentes del vector gradiente fx y fy, el punto inicial x0, y0, la tolerancia tol y el número máximo de iteraciones y como argumentos de salida un vector (array numpy) x y otro y que contengan las coordenadas $x$ e $y$ de las sucesivas iteraciones y el número de iteraciones realizadas. La función optimizará cada paso utilizando seccionaruea para hallar el h óptimo.
Calcular el mínimo de la función $f_1(x,y) = x^2 + 3y^2$, utilizando
- Punto inicial $(2,1)$
- Tolerancia
tol = 0.001tal que que $\left\Vert \nabla f \right\Vert\lt tol$
Calcular el mínimo de la función $f_2(x,y) = x^2 + y^2 - x\,y - 3\,x$
- Punto inicial $(-2,2)$
- Tolerancia
tol = 0.001tal que $\left\Vert \nabla f \right\Vert\lt tol$
En ambos casos dibujar la función y la sucesión de puntos.
%run Ejercicio3.py
El método de Newton¶
Recordamos el método de Newton para encontrar una raíz de la ecuación no lineal $f(x)=0$. Partiendo de un valor inicial $x_0$, usando la fórmula
$$ x_{k+1} = x_k-\frac{f(x_k)}{f'(x_k)} $$y obtenemos
$$ x_{1} = x_0-\frac{f(x_0)}{f'(x_0)},\qquad x_{2} = x_1-\frac{f(x_1)}{f'(x_1)}, \qquad x_{3} = x_2-\frac{f(x_2)}{f'(x_2)},\quad \ldots \quad \rightarrow \quad \alpha $$Si quisieramos encontrar un máximo o un mínimo de la función $f$, si $f$ es suficientemente regular, la condición necesaria de extremo es
$$ f'(x)=0 $$Es decir, buscamos una raíz de esta ecuación y, usando Newton, podríamos resolverlo, tomando un valor inicial $x_0$, y usando la fórmula
$$ x_{k+1} = x_k-\frac{f'(x_k)}{f''(x_k)} $$De forma análoga, para funciones de $n$ variables $f:\Omega\subset \mathbb{R}^n\to \mathbb{R}$, teniendo en cuenta que las componentes del vector gradiente de una función $\nabla f(\mathbf{x}$ son sus derivadas parciales primeras, y las componentes de la matriz Hessiana $H(\mathbf{x})$ son las derivadas parciales segundas, podemos formular Newton partiendo de un valor inicial $\mathbf{x}_0$
$$ \mathbf{x}_{k+1} = \mathbf{x}_{k}-H^{-1}(\mathbf{x}_{k}) \cdot \nabla f(\mathbf{x}_{k}) $$Esto no es una demostración. Tanto en el caso unidimensional como en el multidimensional se puede construir el método de Newton a partir de la fórmula de Taylor correspondiente, y esa sería la demostración formal.
Como en el ejercicio siguiente vamos a usar funcione de dos variables, realizaremos una iteración por el método de Newton usando la fórmula
$$ \left(\begin{array}{c} x_{1}\\ y_{1} \end{array}\right)=\left(\begin{array}{c} x_{0}\\ y_{0} \end{array}\right)-H^{-1}\left(x_{0},y_{0}\right)\cdot\nabla f\left(x_{0},y_{0}\right) $$Si consideramos que
$$ H\left(x_{0},y_{0}\right)\left(\begin{array}{c} c_{1}\\ c_{2} \end{array}\right)=\nabla f\left(x_{0},y_{0}\right)\qquad (1) $$entonces
$$ \left(\begin{array}{c} c_{1}\\ c_{2} \end{array}\right)=H^{-1}\left(x_{0},y_{0}\right)\cdot\nabla f\left(x_{0},y_{0}\right) $$y escribiremos
$$ \left(\begin{array}{c} x_{1}\\ y_{1} \end{array}\right)=\left(\begin{array}{c} x_{0}\\ y_{0} \end{array}\right)-\left(\begin{array}{c} c_{1}\\ c_{2} \end{array}\right)\qquad (2) $$donde $\left(c_{1},c_{2}\right)^{T}$ es la solución del sistema (1). En general, tiene más sentido resolver el sistema (1) que calcular la matriz inversa de $H$ y luego multiplicarla por el gradiente, porque calcular la inversa de una matriz equivale a resolver $n$ sistemas (aunque con la misma matriz de coeficientes) y de esta forma estamos resolviendo solo un sistema.
Escribir una función newton(f,gradf,Hf,x0,tol=1.e-6,maxiter=100) que tenga como argumentos de entrada la función lambda de dos variables f, la función vector gradiente de $f$ gradf, la matriz hessiana de $f$ Hf, el array numpy que contiene el punto inicial x0, la tolerancia tol y el número máximo de iteraciones y como argumentos de salida un vector (array numpy) x y otro y que contengan las coordenadas $x$ e $y$ de las sucesivas iteraciones y el número de iteraciones realizadas.
Calcular los mínimos de la función $$f(x,y) = \ln(12+2\,x^2-4\,x\,y+y^4)$$
Los elementos del vector que devuelve
gradfy de la matriz que devuelveHfse definirán con fórmulas centradas de derivación numérica.La tolerancia será
tol = 0.001tal que $\left\Vert \nabla f \right\Vert\lt tol$Intentar con distintos puntos iniciales hasta que se llegue a los dos mínimos de la función.
Dibujar la función y la sucesión de puntos.
%run Ejercicio4.py
El punto $(0,0)$ es un punto de silla porque es un máximo si cortamos la superficie en la dirección $y=x$ y un mínimo en la dirección $y=-x$
g1 = lambda h: f(h,h)
g2 = lambda h: f(h,-h)
xx = np.linspace(-1,1)
plt.figure()
plt.plot(xx,g1(xx),label='Corte de $f$ en la dirección $y = x$')
plt.plot(xx,g2(xx),label='Corte de $f$ en la dirección $y = -x$')
plt.plot(0,g1(0),'ro')
plt.legend()
plt.show()
La función de penalización¶
La idea del método de penalización es reemplazar la función objetivo $f$ por otra función
$$F(x,y)=f(x,y)+c\,P(x,y)$$y resolver el problema sin restricciones. Para ello tomamos $c$ como una constante positiva y $P$ satisfaciendo:
$P$ es continua en el dominio de $f$.
$P(x,y)\geq 0$ para todo punto del domino de $f$, y
$P(x,y)=0$ si y solo si el punto $(x,y)$ satisface las restricciones.
Minimizar la función $f(x,y)=-y$ sujeta a la condición $x^{2}+y^{2}=1$ utilizando función de penalización y el método del descenso más pronunciado
Las componentes del gradiente
fxyfyse definirán con fórmulas centradas de derivación numérica.La tolerancia será
tol = 0.01tal que $\left\Vert \nabla f \right\Vert\lt tol$Intentar con distintos puntos iniciales hasta que se llegue al mínimo de la función.
Dibujar la función y la sucesión de puntos.
Nota: Una posible función para aproximar el mínimo con restricciones sería
$$ F(x,y) = -y + 10\,(x^2+y^2-1)^2$$%run Ejercicio5.py
Minimizar la función $f(x,y)=x^{2}+y^{2}+xy-3x$ sujeta a las condiciones $x\geq0$ e $y\geq0$ utilizando la función de penalización y el método del descenso más pronunciado
Las componentes del gradiente
fxyfyse definirán con fórmulas centradas de derivación numérica.La tolerancia será
tol = 0.001tal que $\left\Vert \nabla f \right\Vert\lt tol$Intentar con distintos puntos iniciales hasta que se llegue al mínimo de la función.
Dibujar la función y la sucesión de puntos.
Nota: Una posible función para aproximar el mínimo con restricciones sería
$$F(x,y) = f(x,y) + 20 g_1(x,y) + 20 g_2(x,y)$$donde
$$ g_1(x,y)=\begin{cases} 0 & \mathrm{si} & x\ge 0\\ x^2 & \mathrm{si} & x\lt 0 \end{cases} \qquad g_2(x,y)=\begin{cases} 0 & \mathrm{si} & y\ge 0\\ y^2 & \mathrm{si} & y\lt 0 \end{cases} $$%run Ejercicio6.py
