Clase (12/04)

Métodos de optimización irrestricta

Algoritmos de descenso

Idea: a partir de un punto obtenido, escoger una dirección para dar el próximo paso

Definición (dirección de descenso): sean $f:\mathbb{R}^n \rightarrow \mathbb{R}$, $\bar{x}\in\mathbb{R}^n$ y $d\in\mathbb{R}^n-\{0\}$, diremos que $d$ es una dirección de descenso para $f$ a partir de $\bar{x}$ si existe $\delta>0$ tal que $f(\bar{x}+td)<f(\bar{x}) \quad \forall t\in(0,\delta)$

Teorema: Si $\nabla f(\bar{x})d < 0$, entonces $d$ es dirección de descenso para $f$ a partir de $\bar{x}$

Algoritmo de descenso básico

Dados: $f,\; x^0 \in \mathbb{R}^n,\; \varepsilon>0,\; k_{MAX}>0$
k = 0
REPETIR mientras $\nabla f(x^k) > \varepsilon$ y $k<k_{MAX}$ :
     Calcular $d^k$ tal que $\nabla f(x^k)^Td^k < 0$
     Escoger $t_k>0$ tal que $f(x^k+t_kd^k)<f(x^k)$
     Hacer $x^{k+1}=x^k + t_kd^k$
     $k = k+1$

Métodos de búsqueda unidireccional

Lo que nosotros querríamos hacer es resolver el siguiente problema:

minimizar $f(x+td)$

sujeto a: $t>0$
Naturalmente, es un objetivo ambiocoso pues no resulta una tarea fácil salvo que $f$ cumpla con características muy específicas. Veremos algoritmos que permiten aproximar a la solución de ese problema.

Búsqueda Exacta - Sección Áurea

De aquí en adelante definimos $\varphi\colon [0,+\infty)\rightarrow \mathbb{R}$ de la siguiente manera: $$\varphi(t)=f(x+td)$$ Definición (Función unimodal): una función continua $\varphi\colon [0,+\infty)\rightarrow \mathbb{R}$ se dice unimodal si admite un conjunto de minimizadores $[t_1, t_2]$ y es estrictamente decreciente en $[0,t_1]$ y estrictamente creciente en $[t_2, +\infty)$

unimodal.png

Cómo encontrar un minimizador de una función unimodal

Supongamos que sabemos que un minimizador se encuentra en $[a,b]$.

  • Considerar $a<u<v<b$
  • Si $\varphi(u) < \varphi(v)$, entonces el intervalo $[v,b]$ no puede contener un minimizador y es descartado
  • Si $\varphi(u) \geq \varphi(v)$, entonces el intervalo $[a,u]$ puede ser descartado
  • Particionar el intervalo resultante y repetir el proceso.

El método de la sección áurea consiste en tomar: $$u = a + \frac{3-\sqrt{5}}{2}(b-a) \qquad v = a + \dfrac{\sqrt{5}-1}{2}(b-a)$$

Lema: en el método de la sección áurea, si $[v,b]$ es descartado, entonces $v^+=u$. Análogamente, si $[a,u]$ es descartado, entonces $u^+=v$. seccion_aurea.png

El algoritmo de la sección áurea consiste en 2 fases: en la primera, se busca un intervalo $[a,b]$ que contenga un minimizador de $\varphi$ ; en la segunda, se reduce el intervalo $[a,b]$ hasta que la precisión deseada $\varepsilon$ es alcanzada.

Algoritmo de sección áurea

Dados: $f, \; x\in\mathbb{R}^n,\; d\in\mathbb{R}^n,\; \varepsilon > 0,\; \rho>0$
Definir $\theta_1 = \frac{3-\sqrt{5}}{2}$, $\theta_2 = 1 - \theta_1$
Definir $a = 0$, $s = \rho$, $b=2\rho$
Definir $\varphi_b = f(x+b*d), \; \varphi_s = f(x+s*d)$
REPETIR mientras $\varphi_b<\varphi_s$:
     $a=s$, $s=b$, $b=2b$
     $\varphi_s = \varphi_b, \; \varphi_b = f(x+b*d)$
Definir $u = a+\theta_1(b-a), \; v=a+\theta_2(b-a)$
Definir $\varphi_u = f(x+u*d), \; \varphi_v = f(x+v*d)$
REPETIR mientras $(b-a)>\varepsilon$:
     Si $\varphi_u<\varphi_v$:
         $b=v, \; v=u,\; u=a+\theta_1(b-a)$
         $\varphi_v = \varphi_u,\; \varphi_u = f(x+u*d)$
     Si no:
         $a=u,\; u=v,\; v= a+\theta_2(b-a)$
         $\varphi_u = \varphi_v,\; \varphi_v = f(x+v*d)$
Definir $\bar{t} = \dfrac{u+v}{2}$

Valores usuales para los parámetros: $\varepsilon = 10^{-5},\; \rho=1$

Cuando $\varphi$ es unimodal, la sección áurea funciona muy bien para encontrar una aproximación al minimizador. Si $\varphi$ no es unimodal, este algoritmo puede no ser eficiente.

Búsqueda inexacta - Condición de Armijo

A diferencia de la sección áurea, este algoritmo no busca minimizar $\varphi(t)$ sino encontrar $t$ tal que haya una buena reducción en la función objetivo. Más específicamente, dados $\bar{x}\in\mathbb{R}^n$, $d$ dirección de descenso y $\eta\in(0,1)$ busca $\bar{t}$ tal que: $$f(\bar{x}+\bar{t}d) \leq f(\bar{x})+\eta\bar{t}\nabla f(\bar{x})^Td$$ es decir, la reducción debe ser proporcional al tamaño del paso.

Interpretación gráfica (Armijo)

Armijo.PNG

Algoritmo de Búsqueda de Armijo

Dados: $f,\;\bar{x}\in\mathbb{R}^n$, $d\in \mathbb{R^n}$ dirección de descenso, $\gamma,\;\eta\in(0,1)$
t=1
REPETIR mientras $f(\bar{x}+td) > f(\bar{x})+\eta t\nabla f(\bar{x})^Td$
     $t=\gamma t$

Valores usuales para los parámetros: $\gamma = 0.7,\; \eta=0.45$

Búsqueda inexacta - Condiciones de Wolfe

La condicion de Armijo impone cotas a cuan grande puede ser el paso en la direccion $d$. Sin embargo, puede ocurrir que el paso sea tan pequeño que $x^k$ no converja a un minimizador local.

Por ejemplo, si $f(x)=x^2$, si comenzamos en $x^0=2$, $d=-1$ es dirección de descenso en $x^k = 1+2^{-k}$ y $t_k = 2^{-k-1}$ cumple con la condición y efectivamente se logra un descenso en $f$. Sin embargo, $x_k\rightarrow 1$ que no es el minimizador de $f$. No se converge al mínimo pues los pasos son muy pequeños.

Entonces, necesitamos otra condición que acote inferiormente a $t$.

Wolfe agrega otra condición sobre $t$, la condición de curvatura: $$\nabla f (\bar{x} + \bar{t}d)^Td \geq \zeta \nabla f(\bar{x})^T d$$ donde $\zeta\in(\eta,1)$ con $\eta$ siendo la constante de la condición de Armijo. El lado izquierdo es $\varphi'(\bar{t})$, por lo que la condición impone que la pendiente de $\varphi$ en $\bar{t}$ sea mayor que $\zeta$ veces la pendiente inicial.

Si $\varphi'$ es muy negativa $\Rightarrow$ se puede decrecer mucho en esta dirección
Si $\varphi'$ no es muy negativa o es positiva $\Rightarrow$ terminar la búsqueda lineal, no se pueden lograr (muchas) mejoras

Condiciones de Wolfe:

$$\begin{array}{rcl} f(\bar{x}+\bar{t}d) &\leq& f(\bar{x})+c_1\bar{t}\nabla f(\bar{x})^Td \\ \nabla f (\bar{x} + \bar{t}d)^Td &\geq& c_2 \nabla f(\bar{x})^T d \end{array}$$

Con $0<c_1<c_2<1$.

Algoritmo de Búsqueda de Wolfe

Dados: $f,\; \bar{x}\in\mathbb{R}^n,\; d$ dirección de descenso,$\; 0<c_1<c_2<1$
Definir $\alpha = 0,\; t=1,\; \beta = +\infty$ (beta=np.inf)
REPETIR
     Si $f(\bar{x}+td) > f(\bar{x})+c_1t\nabla f(\bar{x})^Td$ :
         Definir $\beta=t,\; t=\frac{1}{2}(\alpha+\beta)$
     De lo contrario, si $\nabla f (\bar{x} + td)^Td < c_2 \nabla f(\bar{x})^T d$:
         Definir $\alpha=t ,\; t=\begin{cases} 2\alpha \quad \text{si } \beta=+\infty \\ \frac{1}{2}(\alpha+\beta) \quad \text{c.c.} \end{cases} $
     Si no: PARAR

Valores usuales para parámetros: $c_1= 0.5,\; c_2=0.75$
El comando break de Python permite salir de un ciclo while o for

Lema: sean $f:\mathbb{R}^n\rightarrow \mathbb{R}$, $f\in C^1$, $\bar{x}\in\mathbb{R}^n$ y $d$ una dirección de descenso, entonces una de las siguientes dos situaciones pueden ocurrir con el método antes expuesto para las condiciones de Wolfe:
i) El procedimiento temrina en una cantidad finita de pasos, devolviendo un valor $\bar{t}$ que satisface las condiciones de Wolfe
ii) El procedimiento no termina: el parámetro $\beta$ nunca toma un valor finito, $\alpha$ se vuelve positivo en la primera iteración y se duplica con las iteraciones siguientes, y $f(x+td)\rightarrow -\infty$

Interpretación gráfica (2da. Condición de Wolfe)

wolfe.PNG

Interpretación gráfica (Condiciones de Wolfe)

AW.PNG

Métodos de Optimización Irrestricta

Método del Gradiente (o de Cauchy)

En este método, se toma como dirección de descenso: $$d = -\nabla f(x)$$ Observar que, si $\nabla f(x) \neq 0$, efectivamente $d$ es dirección de descenso: $$\nabla f(x)^T d = -\left\Vert \nabla f(x) \right\Vert ^2 < 0$$

Algoritmo del Método del Gradiente

Dados: $f, x^0 \in \mathbb{R}^n,\; \varepsilon>0,\; k_{MAX}>0 $
$k = 0$
REPETIR mientras $\lVert\nabla f(x^{k})\rVert > \varepsilon$ y $k<k_{MAX}$ :
     Definir $d^k = -\nabla f(x^k)$
     Obtener $t_k>0$ tal que $f(x^k+t_kd^k)<f(x^k)$ (con sección áurea, Armijo o Wolfe)
     Hacer $x^{k+1} = x^k + t_kd^k$
     $k = k+1$

Lema: en el algoritmo anterior, sea $t_k$ el obtenido por una minimiación local de $f(x^k+td^k)$, entonces $<d^{k+1}, d^k> = 0$

Demo: como $t_k$ es minimizador local de $\varphi$ tenemos que: $$0 =\varphi'(t_k) = \nabla f(x^k +t_kd^k)^Td^k=\nabla f(x^{k+1})^Td^k= (-d^{k+1})^Td^k$$ metodo_grad.png

Método de Newton

En este método, se elige como dirección de descenso: $$d = -Hf(x)^{-1}\nabla f(x)$$ Notar que, si $Hf(x)$ es definida positiva, $d$ es dirección descenso: $$\nabla f(x)^T d= -\nabla f(x)^THf(x)^{-1}\nabla f(x) < 0$$

Observación: $d$ puede no estar bien definido si $Hf(x)$ es singular o puede no ser una dirección de descenso si $Hf(x)$ no es definida positiva.

Teorema: si $Hf(x)$ es definida positiva para todo $x\in\mathbb{R}^n$, el Método de Newton es globalmente convergente.

Algoritmo del Método de Newton

Dados: $f, x^0 \in \mathbb{R}^n,\; \varepsilon>0,\; k_{MAX}>0$
$k=0$
REPETIR mientras $\lVert\nabla f(x^{(k)})\rVert > \varepsilon$ y $k<k_{MAX}$:
     Definir $d^{(k)}$ como la solución de $Hf(x^{(k)})d^{(k)} = -\nabla f(x^{(k)})$
     Determinar el tamaño del paso $t_k > 0$
     Hacer $x^{(k+1)} = x^{(k)} + t_kd^{(k)}$
     k = k + 1

¿Qué pasa si $Hf(x)$ no es definida positiva o si es singular?

thinking.gif

Modificación de Levenberg-Marquardt

Considerar $$x^{(k+1)} = x^{(k)} - t_k(Hf(x^{(k)})+\mu I)^{-1}\nabla f(x^{(k)})$$ donde $\mu>0$ es lo suficientemente grande como para que $B^{(k)}= Hf(x^{(k)})+\mu I$ sea definida positiva.

Obs: si $\mu$ es demasiado grande, $d^{(k)}$ tiende a estar en la dirección de $-\nabla f(x^{(k)})$, pues $(Hf(x^{(k)})+\mu I)^{-1}\nabla f(x^{(k)}) \approx \frac{1}{\mu}\nabla f(x^{(k)})$ . En cambio, si $\mu$ es pequeño, $d^{(k)}$ se parece más a la dirección de descenso del Método de Newton

Con la modificación de L-M, se obtiene un método que utiliza las ventajas del Método del Gradiente y del Método de Newton. Pensándolo con ideas generales:

  • cerca del mínimo $\approx$ $Hf(x)$ definida positiva $\approx$ $\mu=0$ o $\mu$ pequeño $\approx$ Método de Newton $\approx$ convergencia rápida cerca de un mínimo
  • lejos del mínimo $\approx$ $Hf(x)$ no definida positiva $\approx$ $\mu$ grande $\approx$ Método del gradiente $\approx$ más velocidad de descenso

Obs: Como pedíamos que $f\in C^3$, en particular resulta que $Hf(x)$ es simétrica $\forall x \in \mathbb{R}^n$, entonces los autovalores de $Hf(x)$ son reales $\forall x \in \mathbb{R}^n$.

Obs: Si $\lambda_1, \dots, \lambda_n$ son los autovalores de $Hf(x)$, entonces $\lambda_1+\mu, \dots, \lambda_n+\mu$ son los autovalores de $B=Hf(x)+\mu I$, y $B$ tiene los mismos autovectores que $Hf(x)$. Sea $v_i$ autovector de $Hf(x)$: $$\begin{array}{rcl} Bv_i &=& (Hf(x)+\mu I)v_i \\ &=& Hf(x)v_i+\mu Iv_i \\ &=& \lambda_iv_i+\mu v_i \\ &=& (\lambda_i+\mu)v_i \end{array}$$ Entonces, en caso de que $Hf(x)$ no sea definida positiva, se puede utilizar su mínimo autovalor como una aproximación a $-\mu$

Algoritmo del Método de Newton con modificación de Levenberg-Marquardt

Dados: $f, x^0 \in \mathbb{R}^n,\; k_{MAX}>0,\; \varepsilon>0,\; \gamma>0$
$k=0$
REPETIR mientras $\lVert\nabla f(x^{(k)})\rVert > \varepsilon$ y $k<k_{MAX}$:
     $B = Hf(x^{(k)})$
     $\mu = \min \{ \lambda \colon \lambda \text{ es autovalor de } B\}$
     Si $\mu \leq 0$:
         $B = B+(-\mu+\gamma)I$
     Definir $d^{(k)}$ como la solución de $Bd^{(k)} = -\nabla f(x^{(k)})$
     Determinar el tamaño del paso $t_k > 0$
     Hacer $x^{(k+1)} = x^{(k)} + t_kd^{(k)}$
     k = k + 1

Utilizar la función eigvals de numpy.linalg para calcular los autovalores.

Comparación de los métodos, aplicados a: $$f(x, y) = 2(e^{-x^2-y^2}- e^{-(x-1)^2 - (y-1)^2})$$

comp_1.png

comp_2.png

comp_3.png

Consignas

Implementar una función para cada una de las siguientes tareas:

  • devolver la longitud del paso utilizando sección áurea
  • devolver la longitud de un paso que cumpla con la Condición de Armijo
  • devolver la longitud de un paso que cumpla con las Condiciones de Wolfe
  • encontrar mínimo utilizando el Método del Gradiente
  • encontrar mínimo utilizando el Método de Newton
  • encontrar mínimo utilizando el Método de Newton con modificación de Levenberg-Marquardt

Las últimas tres deben tener metodo como un input opcional que permita elegir qué metodo utilizar (sección áurea, Armijo o Wolfe) para encontrar la longitud del paso. Deben imprimir la cantidad de iteraciones que realizaron, el minimizador y el valor de la función en el mismo.

Aplicar los tres métodos sobre la función de Rosenbrock en 2 dimensiones y comparar los resultados. El minimizador de esta función es $x^*=(1,1)$ y $f(x^*)=0$

Se pueden aplicar a otras funciones del archivo fx_test.py. Se pueden importar todas las funciones de ese archivo con el siguiente comando

In [ ]:
from fx_test import *

Para utilizar la librería plotly hay que instalar también jupyter. Para instalar ambos paquetes, desde la terminal se ejecuta, para Python 2.7 : pip install jupyter y luego pip install plotly.

Para Python 3+: pip3 install jupyter y luego pip3 install plotly

Dependiendo de su configuración, podría ser necesario anteponer el comando sudo

Para graficar las funciones del archivo fx_test pueden usar la función plot_function que se encuentra en ese archivo. Por ejemplo, si se desea graficar la función de ackley:

In [ ]:
plot_function(ackley)

Pueden modificar los límites del gráfico con el input opcional limits e ingresando una tupla $(x_1, x_2, y_1, y_2)$ para que se grafique la función en el dominio $[x_1, x_2]\times [y_1, y_2]$. Por ejemplo:

In [ ]:
plot_function(ackley, limits=(-1,1,-3,2))

grafica la función de Ackley en $[-1,1]\times [-3,2]$. El argumento opcional points toma una lista de puntos $(x,y)$ y los grafica sobre la superficie de la función (se ve mejor utilizando plotly