Gradientes Conjugados¶
Definición: sea $A\in\mathbb{R}^{n\times n}$ definida positiva, diremos que los vectores $d^0, d^1,\dots, d^{n-1}\in\mathbb{R}^n-\{0\}$ son A-conjugados si: $$ (d^i)^TAd^j = 0$$
Lema: sea $A\in\mathbb{R}^{n\times n}$ definida positiva, cualquier conjunto de vectores A-conjugados es linealmente independiente.
Sea $f\colon \mathbb{R}^n \rightarrow \mathbb{R}$ dada por: $$\frac{1}{2}x^TAx + bx^T +c$$ con $A\in\mathbb{R}^{n\times n}$ definida positiva, $b\in\mathbb{R}^n$, $c\in \mathbb{R}$.
Sea $\{d^0, \dots, d^k\}$ un conjunto de vectores A-conjugados cualquiera, la idea consiste en tomar: $$x^{k+1} = x^{k} + t_kd^k \quad \text{donde} \quad t_k = \underset{t\in\mathbb{R}}{\mathrm{argmin}}\{f(x^k + td^k)\}$$ Notar que, como no sabemos si los $d^k$ son direcciones de descenso, no pedimos que $t_k$ sea necesariamente positivo. Nuevamente, definimos $\varphi(t) = f(x^k + t d^k)$ y, como $f$ es cuadrática, podemos hallar el $t$ óptimo: $$\begin{array}{l} 0 = \varphi'(t_k) = \nabla f(x^{k} + t_kd^k)^Td^k = \nabla f(x^{k+1})^Td^k \\ \nabla f(x^{k+1}) = Ax^{k+1} + b =A(x^k + t_kd^k) + b = \nabla f(x^k) + t_kAd^k \\ \Rightarrow t_k = - \dfrac{\nabla f(x^k)^T d^k}{(d^k)^TAd^k}\end{array}$$
Teorema: Sea $f\colon \mathbb{R}^n \rightarrow \mathbb{R}$ dada por $\frac{1}{2}x^TAx + bx^T +c$ con $A\in\mathbb{R}^{n\times n}$ definida positiva, $b\in\mathbb{R}^n$, $c\in \mathbb{R}$, entonces dado $x^0\in\mathbb{R}^n$, la secuencia definida anteriormente alcanza el minimizador $x^\ast$ en $n$ pasos (i.e. $x^n = x^\ast$).
Lema: sea $x^0\in \mathbb{R}^n$, la secuencia definida anteriormente cumple que: $$ \nabla f(x^k)^Td^j = 0 \quad \forall j=0,1,\dots,k-1 $$
Gradientes conjugados¶
La idea consiste en que las direcciones conjugadas $\{d^0, \dots, d^k\}$ sean definidas a partir del gradiente de $f$. Sea $x_0\in\mathbb{R}^n$, se define $d^0 = -\nabla f(x_0)$ y, para $k=0,1,\dots,n-2$ se define: $$d^{k+1} = -\nabla f(x^{k+1}) + \beta_kd^k$$ donde $\beta_k$ es calculado de manera tal que $d^k$ y $d^{k+1}$ sean A-conjugados: $$\beta_k = \dfrac{(d^k)^T A \nabla f(x^{k+1})}{(d^k)^TAd^k}$$
Obs: $d^k$ son direcciones de descenso: $$\nabla f(x^k)^T d^k = \nabla f(x^k)^T (-\nabla f(x^k)+\beta_{k-1}d^{k-1}) = -\lVert \nabla f(x^k) \rVert ^2 $$
Algoritmo de Gradientes Conjugados (para cuadráticas)¶
Dados: $f,\; x^0 \in \mathbb{R}^n, \; \varepsilon > 0, \; k_{MAX} > 0$
$d^0 = -\nabla f(x^0)$
$k = 0$
REPETIR mientras $\lVert \nabla f(x^k) \rVert > \varepsilon$ y $k<k_{MAX}$
$t_k = - \dfrac{\nabla f(x^k)^T d^k}{(d^k)^TAd^k}$
$x^{k+1} = x^k + t_kd^k$
$\beta_k = \dfrac{(d^k)^T A \nabla f(x^{k+1})}{(d^k)^TAd^k}$
$d^{k+1} = -\nabla f(x^{k+1}) + \beta_kd^k$
$k = k+1$
Extensión para funciones no cuadráticas¶
Se deben encontrar otras maneras de obtener la longitud del paso $t_k$ y el coeficiente $\beta_k$
Obs: si $t_k$ minimiza $\varphi$ entonces las $d^k$ son direcciones de descenso. Sin embargo, la búsqueda de Armijo no alcanza para asegurar tal propiedad.
$$t_k \rightarrow \color{green}{\text{sección áurea}}\text{, Armijo, }\color{green}{\text{Wolfe}}$$Para definir el coeficiente $\beta_k$ se puede utilizar la fórmula de Fletcher y Reeves:
$$\beta^{FR}_k = \dfrac{\nabla f(x^{k+1})^T \nabla f(x^{k+1})}{\nabla f(x^{k})^T \nabla f(x^{k})} $$Para el caso de funciones cuadráticas, se puede demostrar que $\beta_k = \beta^{FR}_k$
Gradientes Conjugados para funciones no cuadráticas no necesariamente terminan en $n$ pasos. Una práctica que da buenos resultados es reiniciar el $\beta_k$ cada $n$ iteraciones.
Algoritmo de Gradientes Conjugados¶
Dados: $f,\; x^0 \in \mathbb{R}^n, \; \varepsilon > 0, \; k_{MAX} > 0$
$d^0 = -\nabla f(x^0)$
$k = 0$
REPETIR mientras $\lVert \nabla f(x^k) \rVert > \varepsilon$ y $k<k_{MAX}$
Calcular $t_k$
Hacer $x^{k+1} = x^k + t_kd^k$
Si $k + 1 \neq 0 \mod(n)$:
$\beta^{FR}_k = \dfrac{\nabla f(x^{k+1})^T \nabla f(x^{k+1})}{\nabla f(x^{k})^T \nabla f(x^{k})} $
Si no:
$\beta_k = 0$
Definir $d^{k+1} = -\nabla f(x^{k+1}) + \beta_kd^k$
$k=k+1$
Métodos de Cuasi-Newton¶
La idea, como en el caso del Método de Newton, es aproximar $f$ mediante una expresión cuadrática. Sin embargo, en los métodos cuasi-Newton, se utiliza, en vez de $Hf(x)$, una matriz $B_k$ que sea simétrica definida positiva: $$m_k(d) = f(x^k) + \nabla f(x^{k})^T d+\frac{1}{2}d^TB_kd$$
Como hemos visto, el minimizador del modelo cuadrático viene dado por: $$-B_k^{-1}\nabla f(x^{k})$$
Notaremos $H_k = B_k^{-1}$.
Algoritmo general de un método Cuasi-Newton¶
Dados: $f,\; x^0 \in \mathbb{R}^n, \; H_0\in \mathbb{R}^{n\times n}\;\;\text{definida positiva},\; \varepsilon > 0, \; k_{MAX} > 0$
$k = 0$
REPETIR mientras $\lVert \nabla f(x^k) \rVert > \varepsilon$, $\lVert x^{k+1} - x^k \rVert > \varepsilon$ y $k<k_{MAX}$:
Definir $d^k = -H_k\nabla f(x^k)$
Obtener $t_k$ (sección áurea, Armijo o Wolfe)
Hacer $x^{k+1} = x^k + t_kd^k$
Determinar $H_{k+1}$ definida positiva
$k=k+1$
Obs: si $H_k = I \rightarrow$ Método del Gradiente
si $H_k = Hf(x^k)^{-1} \rightarrow$ Método de Newton
Obs: como $H_0$ podemos tomar la identidad, ya que es definida positiva. También podrÃa tomarse $Hf(x_0)$ si es definida positiva.
Sean $q^k = \nabla f(x^{k+1}) - \nabla f(x^k)$ y $p^k = x^{k+1} - x^k$, se demuestra que una propiedad que debe satisfacer $H_{k+1}$ definida por el algoritmo es que: $$H_{k+1}q^j = p^j \quad \forall j=0,1,\dots,k$$
Método DFP (Davindon, Fletcher y Powell)¶
$$H_{k+1} = H_k + \dfrac{p^k(p^k)^T}{(p^k)^Tq^k} - \dfrac{H_kq^k(q^k)^TH_k}{(q^k)^TH_kq^k}$$Método BFGS (Broyden, Fletcher, Goldfarb y Shanno)¶
$$H_{k+1} = H_k + \left(1+\dfrac{(q^k)^TH_kq^k}{(p^k)^Tq^k}\right)\dfrac{p^k(p^k)^T}{(p^k)^Tq^k} - \dfrac{p^k(q^k)^TH_k + H_kq^k(p^k)^T}{(p^k)^Tq^k} $$Producto externo en Python: numpy.outer(a,b)
Se puede demostrar que si $t_k$ es mÃnimo local de $\varphi$ entonces $H_{k+1}$ de ambos métodos es definida positiva.
Región de confianza¶
Idea: se establece un modelo de la función objetivo en un entorno del punto actual. Se calcula un minimizador aproximado del modelo en la región de confianza establecida. Si el minimizador ofrece una reducción razonable de la función objetivo, nos movemos a ese punto y repetimos el proceso. De lo contrario, se achica la región de confianza y se busca un nuevo minimizador.
Métodos antes vistos: establecer dirección $\rightarrow$ establecer longitud de paso
Región de confianza: establecer longitud de paso $\rightarrow$ establecer dirección
Dado $\Delta_k>0$ consideramos la región de confianza: $$\{ x\in\mathbb{R}^n \colon \lVert x - x^k \rVert \leq \Delta_k\} $$ Modelo cuadrático para aproximar $f$ en dicha región: $$q_k(x) = f(x^k) + \nabla f(x^k)^T (x-x^k) + \frac{1}{2}(x-x^k)^TB_k(x-x^k)$$ siendo $B_k$ una matriz lo suficientemente buena para aproximar $f$ (por ejemplo, $B_k = Hf(x^k)$)
Sean: $$d = x - x^k \qquad m_k(d) = q_k(x^k+d)$$
En la primera etapa del método se resuelve (de manera aproximada) el siguiente subproblema: $$\begin{array}{rl} \min & m_k(d) = f(x^k) + \nabla f(x^k)^Td + \frac{1}{2}d^TB_kd \\ \text{s.a.} & \lVert d \rVert \leq \Delta_k\end{array}$$ Obteniendo de esta manera el valor de $d^k$.
La siguiente etapa es evaluar cuánto mejora la función en $x^k + d^k$. Resulta razonable pedir que la mejora en la función objetivo sea comparable a la mejora obtenida en el modelo (función cuadrática). Se definen: $$ared = f(x^k) - f(x^k+d^k) \qquad pred = m_k(0)-m_k(d^k) \qquad \rho_k = \dfrac{ared}{pred}$$
Se aceptará el nuevo punto $x^k + d^k$ si $\rho_k > \eta$ para cierto $\eta>0$ constante. Si no se acepta el nuevo punto, se reduce $\Delta_k$ y se vuelve a resolver el subproblema.
Solución aproximada del subproblema: Punto de Cauchy¶
Notación: $f_k = f(x^k) \qquad g_k = \nabla f(x^k)$
Idea:
Primero, hallar la solución de la versión lineal del subproblema: $$d_k^S = \underset{d\in\mathbb{R}^n}{\mathrm{argmin}} f_k + g_k^Td \qquad \text{s.a } \lVert d \rVert \leq \Delta_k $$ La solución es: $$d_k^S = -\frac{\Delta_k}{\lVert g_k \rVert} g_k $$
El segundo paso es encontrar el escalar $\tau_k$ tal que: $$\tau_k = \underset{\tau\geq 0}{\mathrm{argmin}}\;m_k(\tau d_k^S) \qquad \text{s.a }\lVert \tau d_k^S \rVert \leq \Delta_k$$
Dos casos:
$\mathbf{g_k^TB_kg_k \leq 0}$ : la función $m_k(\tau d_k^S)$ decrece monótonamente cuando $g_k\neq 0$, luego $\tau_k$ puede ser tan grande como lo permita la restricción, es decir, $\tau_k = 1$
$\mathbf{g_k^TB_kg_k > 0}$ : $m_k(\tau d_k^S)$ resulta cuadrática convexa sobre $\tau$, luego $\tau_k$ es el minimizador de la cuadrática irrestricta ( o sea, $\lVert g_k\rVert^3/(\Delta_kg_k^TB_kg_k)$) o tan grande como lo permita la restricción, lo que ocurra primero.
Se tiene entonces que, una aproximación a la solución del subproblema es: $$d_k^C = -\tau_k\frac{\Delta_k}{\lVert g_k \rVert} g_k$$ donde: $$\tau_k = \begin{cases} 1 \quad \text{si} \quad g_k^TB_kg_k \leq 0 \\ \min\left(\frac{\lVert g_k\rVert^3}{\Delta_kg_k^TB_kg_k},1\right) \quad \text{c.c.}\end{cases}$$
Algoritmo de Región de confianza con Punto de Cauchy¶
Utilizaremos $B_k = Hf(x^k)$
Dados: $f,\; x^0 \in \mathbb{R}^n,\; \Delta_0 > 0,\; \eta\in [0,\frac{1}{4}), \; \varepsilon > 0, \; k_{MAX} > 0$
$k = 0$
REPETIR mientras $\lVert \nabla f(x^k) \rVert > \varepsilon$ y $k<k_{MAX}$:
Calcular $d_k$ con Punto de Cauchy
Calcular $\rho_k$
Si $\rho_k>\eta$:
$x^{k+1} = x^k + d^k$
Si no:
$x^{k+1} = x^k$
Si $\rho_k<\frac{1}{4}$:
$\Delta_{k+1} = \frac{\Delta_k}{2}$
Si no:
Si $\rho_k > \frac{3}{4}$ y $\lVert d_k \rVert = \Delta_k$:
$\Delta_{k+1} = 2\Delta_k$
Si no:
$\Delta_{k+1} = \Delta_k$
$k = k+1$
Posibles valores por defecto: $\Delta_0 = 1, \; \eta=0.2$
Consignas¶
- Implementar el método de gradientes conjugados para funciones en general. Permitir elegir qué método utilizar para encontrar la longitud del paso, el método por defecto debe ser Wolfe.
- Implementar el método de Cuasi-Newton que considera como $H_0$ a la matriz identidad y utiliza el método DFP para actualizar $H_k$.
- Implementar el método de región de confianza con punto de Cauchy y $B_k = Hf(x_k)$
Perfil de desempeño¶
Diseñado por Dolan y Moré en 2001 y es una herramienta que se emplea para comparar la aptitud de distintos algoritmos para resolver un conjunto de problemas.
Dados un conjunto de problemas $P$ y un conjunto de algoritmos $S$, definimos como $c_{s,p}$ el costo de resolver el problema $p \in P$ con el algoritmo $s \in S$. Si el algoritmo $s$ no puede resolver el problema $p$, se define $c_{s,p}=M$ con $M\in\mathbb{R}$ un valor adecuado de manera tal que: $$c_{s,p}=M \Leftrightarrow \text{el algoritmo $s$ no resuelve el problema $p$}$$
Se asume que al menos un algoritmo resuelve el problema $p$. Definimos el Ãndice de desempeño relativo como: $$r_{s,p}= \frac{c_{s,p}}{\min \{c_{j,p}\colon j\in S\}}$$
Observar que $r_{s,p}\geq 1$ y, si $r_{s,p}=1$, entonces es algoritmo $s$ es uno de los mejores (o el mejor) para resolver el problema $p$. Mientras mayor sea $r_{s,p}$, peor es el desempeño de $s$ al resolver $p$. Finalmente, definimos la función de desempeño $\rho_s:[1,+\infty)\rightarrow [0,1]$ para el algoritmo $s$ como: $$\rho_s(\tau) = \frac{\#\{p\in P\colon r_{s,p}\leq \tau\}}{\# P}$$
AsÃ, $\rho_s(\tau)$ es el porcentaje de problemas que el algoritmo $s$ resuelve con hasta $\tau$ veces el costo del algoritmo más eficiente. Por ejemplo, $\rho_s(1)$ es la proporción de problemas que el algoritmo $s$ resuelve con el menor costo.
Si se define: $$r_{max} = \underset{s\in S, p\in P\;\colon\; c_{s,p}<M}{\max} r_{s,p}$$ se tiene que $\rho_s(r_{max})$ es el número de problemas resueltos por el algoritmo $s$. El valor $\rho_s(1)$ se denomina la eficiencia de $s$ y $\rho_s(r_{max})$ es su robustez.
La definición del costo está relacionada a la medida de desempeño que se quiera estudiar. Se puede considerar como costo al tiempo que necesita el algoritmo para resolver el problema, a la cantidad de iteraciones, a la cantidad de veces que es evaluada la función objetivo, etc.
Para analizar el perfil de desempeño se suelen graficar en conjunto las funciones $\rho_s$ para cada $s\in S$.
Ejemplo¶
Supongamos que tenemos cuatro algoritmos $S = \{A_0,\; A_1,\; A_2,\; A_3\}$ y queremos comparar su desempeño en base a seis problemas $P =\{P_0,\; P_1,\; P_2,\;P_3,\;P_4,\;P_5,\;\}$. Consideraremos como costo $c_{s,p}$ a la cantidad de iteraciones que necesita el algoritmo $s$ para resolver el problema $p$. Si no puede resolverlo, asignamos $c_{s,p}=10^8$. Por "no puede resolverlo" entendemos que el algoritmo $s$ supera el lÃmite de iteraciones.
Obviamente, el primer paso consiste en ejecutar cada uno de los cuatro algoritmos sobre cada problema y registrar el costo (iteraciones) de cada uno. AsÃ, podemos formar la matriz $C$ de costos:
Ejemplo¶
Sobre un conjunto de problemas, se realiza el perfil de desempeño de cuatro algoritmos: LANCELOT, LOQO, MINOS, SNOPT. Se considera como medida de desempeño al tiempo de ejecución.
Observando el valor de $\rho_s(1)$ para cada $s$, observamos que LOQO resuelve aproximadamente el $61\%$ de los problemas en el menor tiempo. En otras palabras, hay un $61\%$ de probabilidades de que LOQO sea el más rápido al resolver un problema.
El perfil de desempeño también muestra, por ejemplo, que LOQO o MINOS pueden resolver alrededor del $70\%$ de los problemas tardando hasta 4 veces más que el algoritmo más rápido. Por otro lado MINOS puede resolver casi el $80\%$ de los problemas tardando hasta 8 veces más que el algoritmo más rápido, mientras que en ese tiempo LOQO sigue resolviendo alrededor del $70\%$ de los problemas.
Si realizamos el mismo gráfico para $\tau\in[1,100]$:
Si nos interesa encontrar un algoritmo que resuelve alrededor del $75\%$ de los problemas con el menor tiempo, MINOS es el mejor candidato. Sin embargo, si buscamos un algoritmo con mayor probabilidad de resolver los problemas, se observa que SNOPT logra resolver el $90\%$.