Algoritmo de Broyden-Fletcher-Goldfarb-Shanno

Método de optimización

En optimización numérica , el algoritmo Broyden–Fletcher–Goldfarb–Shanno ( BFGS ) es un método iterativo para resolver problemas de optimización no lineal sin restricciones . ^[1] Al igual que el método relacionado Davidon–Fletcher–Powell , BFGS determina la dirección de descenso preacondicionando el gradiente con información de curvatura. Lo hace mejorando gradualmente una aproximación a la matriz hessiana de la función de pérdida , obtenida solo a partir de evaluaciones de gradiente (o evaluaciones de gradiente aproximadas) a través de un método secante generalizado . ^[2]

Dado que las actualizaciones de la matriz de curvatura BFGS no requieren inversión de matriz , su complejidad computacional es solo , en comparación con el método de Newton . También se usa comúnmente L-BFGS , que es una versión de memoria limitada de BFGS que es particularmente adecuada para problemas con cantidades muy grandes de variables (por ejemplo, >1000). La variante BFGS-B maneja restricciones de caja simples. ^[3] . La matriz BFGS también admite una representación compacta , lo que la hace más adecuada para problemas grandes con restricciones. ${\displaystyle {\mathcal {O}}(n^{2})}$ ${\displaystyle {\mathcal {O}}(n^{3})}$

El algoritmo lleva el nombre de Charles George Broyden , Roger Fletcher , Donald Goldfarb y David Shanno . ^{[4] [5] [6] [7]}

Razón fundamental

El problema de optimización consiste en minimizar , donde es un vector en , y es una función escalar diferenciable. No existen restricciones sobre los valores que puede tomar. ${\displaystyle f(\mathbf {x} )}$ ${\displaystyle \mathbf {x} }$ ${\displaystyle \mathbb {R} ^{n}}$ ${\displaystyle f}$ ${\displaystyle \mathbf {x} }$

El algoritmo comienza con una estimación inicial del valor óptimo y avanza iterativamente para obtener una mejor estimación en cada etapa. ${\displaystyle \mathbf {x} _{0}}$

La dirección de búsqueda p _k en la etapa k viene dada por la solución del análogo de la ecuación de Newton:

${\displaystyle B_{k}\mathbf {p} _{k}=-\nabla f(\mathbf {x} _{k}),}$

donde es una aproximación a la matriz hessiana en , que se actualiza iterativamente en cada etapa, y es el gradiente de la función evaluada en x _k . Luego se utiliza una búsqueda de línea en la dirección p _k para encontrar el siguiente punto x _{k +1} minimizando sobre el escalar ${\displaystyle B_{k}}$ ${\displaystyle \mathbf {x} _{k}}$ ${\displaystyle \nabla f(\mathbf {x} _{k})}$ ${\displaystyle f(\mathbf {x} _{k}+\gamma \mathbf {p} _{k})}$ ${\displaystyle \gamma >0.}$

La condición cuasi-Newton impuesta a la actualización de es ${\displaystyle B_{k}}$

${\displaystyle B_{k+1}(\mathbf {x} _{k+1}-\mathbf {x} _{k})=\nabla f(\mathbf {x} _{k+1})-\nabla f(\mathbf {x} _{k}).}$

Sea y , entonces satisface ${\displaystyle \mathbf {y} _{k}=\nabla f(\mathbf {x} _{k+1})-\nabla f(\mathbf {x} _{k})}$ ${\displaystyle \mathbf {s} _{k}=\mathbf {x} _{k+1}-\mathbf {x} _{k}}$ ${\displaystyle B_{k+1}}$

${\displaystyle B_{k+1}\mathbf {s} _{k}=\mathbf {y} _{k}}$,

cual es la ecuación secante.

La condición de curvatura debe cumplirse para que sea definida positiva, lo que se puede verificar multiplicando previamente la ecuación secante por . Si la función no es fuertemente convexa , entonces la condición debe cumplirse explícitamente, por ejemplo, encontrando un punto x _{k +1} que satisfaga las condiciones de Wolfe , que implican la condición de curvatura, mediante búsqueda de línea. ${\displaystyle \mathbf {s} _{k}^{\top }\mathbf {y} _{k}>0}$ ${\displaystyle B_{k+1}}$ ${\displaystyle \mathbf {s} _{k}^{T}}$

En lugar de requerir que la matriz hessiana completa en el punto se calcule como , la matriz hessiana aproximada en la etapa k se actualiza mediante la adición de dos matrices: ${\displaystyle \mathbf {x} _{k+1}}$ ${\displaystyle B_{k+1}}$

${\displaystyle B_{k+1}=B_{k}+U_{k}+V_{k}.}$

Tanto y son matrices simétricas de rango uno, pero su suma es una matriz de actualización de rango dos. Las matrices de actualización de BFGS y DFP se diferencian de su predecesora por una matriz de rango dos. Otro método de rango uno más simple se conoce como método de rango uno simétrico , que no garantiza la definibilidad positiva . Para mantener la simetría y la definibilidad positiva de , la forma de actualización puede elegirse como . Imponiendo la condición secante, . Eligiendo y , podemos obtener: ^[8] ${\displaystyle U_{k}}$ ${\displaystyle V_{k}}$ ${\displaystyle B_{k+1}}$ ${\displaystyle B_{k+1}=B_{k}+\alpha \mathbf {u} \mathbf {u} ^{\top }+\beta \mathbf {v} \mathbf {v} ^{\top }}$ ${\displaystyle B_{k+1}\mathbf {s} _{k}=\mathbf {y} _{k}}$ ${\displaystyle \mathbf {u} =\mathbf {y} _{k}}$ ${\displaystyle \mathbf {v} =B_{k}\mathbf {s} _{k}}$

${\displaystyle \alpha ={\frac {1}{\mathbf {y} _{k}^{T}\mathbf {s} _{k}}},}$

${\displaystyle \beta =-{\frac {1}{\mathbf {s} _{k}^{T}B_{k}\mathbf {s} _{k}}}.}$

Finalmente, sustituimos y en y obtenemos la ecuación de actualización de : ${\displaystyle \alpha }$ ${\displaystyle \beta }$ ${\displaystyle B_{k+1}=B_{k}+\alpha \mathbf {u} \mathbf {u} ^{\top }+\beta \mathbf {v} \mathbf {v} ^{\top }}$ ${\displaystyle B_{k+1}}$

${\displaystyle B_{k+1}=B_{k}+{\frac {\mathbf {y} _{k}\mathbf {y} _{k}^{\mathrm {T} }}{\mathbf {y} _{k}^{\mathrm {T} }\mathbf {s} _{k}}}-{\frac {B_{k}\mathbf {s} _{k}\mathbf {s} _{k}^{\mathrm {T} }B_{k}^{\mathrm {T} }}{\mathbf {s} _{k}^{\mathrm {T} }B_{k}\mathbf {s} _{k}}}.}$

Algoritmo

Considere el siguiente problema de optimización sin restricciones donde es una función objetivo no lineal. $${\displaystyle {\begin{aligned}{\underset {\mathbf {x} \in \mathbb {R} ^{n}}{\text{minimize}}}\quad &f(\mathbf {x} ),\end{aligned}}}$$ ${\displaystyle f:\mathbb {R} ^{n}\to \mathbb {R} }$

A partir de una estimación inicial y una estimación inicial de la matriz hessiana, se repiten los siguientes pasos a medida que converge a la solución: ${\displaystyle \mathbf {x} _{0}\in \mathbb {R} ^{n}}$ ${\displaystyle B_{0}\in \mathbb {R} ^{n\times n}}$ ${\displaystyle \mathbf {x} _{k}}$

1. Obtenga una dirección resolviendo . ${\displaystyle \mathbf {p} _{k}}$ ${\displaystyle B_{k}\mathbf {p} _{k}=-\nabla f(\mathbf {x} _{k})}$
2. Realice una optimización unidimensional ( búsqueda lineal ) para encontrar un tamaño de paso aceptable en la dirección encontrada en el primer paso. Si se realiza una búsqueda lineal exacta, entonces . En la práctica, una búsqueda lineal inexacta suele ser suficiente, con un tamaño de paso aceptable que satisfaga las condiciones de Wolfe . ${\displaystyle \alpha _{k}}$ ${\displaystyle \alpha _{k}=\arg \min f(\mathbf {x} _{k}+\alpha \mathbf {p} _{k})}$ ${\displaystyle \alpha _{k}}$
3. Establecer y actualizar . ${\displaystyle \mathbf {s} _{k}=\alpha _{k}\mathbf {p} _{k}}$ ${\displaystyle \mathbf {x} _{k+1}=\mathbf {x} _{k}+\mathbf {s} _{k}}$
4. ${\displaystyle \mathbf {y} _{k}={\nabla f(\mathbf {x} _{k+1})-\nabla f(\mathbf {x} _{k})}}$.
5. ${\displaystyle B_{k+1}=B_{k}+{\frac {\mathbf {y} _{k}\mathbf {y} _{k}^{\mathrm {T} }}{\mathbf {y} _{k}^{\mathrm {T} }\mathbf {s} _{k}}}-{\frac {B_{k}\mathbf {s} _{k}\mathbf {s} _{k}^{\mathrm {T} }B_{k}^{\mathrm {T} }}{\mathbf {s} _{k}^{\mathrm {T} }B_{k}\mathbf {s} _{k}}}}$.

La convergencia se puede determinar observando la norma del gradiente; dado algún , se puede detener el algoritmo cuando Si se inicializa con , el primer paso será equivalente a un descenso del gradiente , pero los pasos posteriores se refinan cada vez más mediante , la aproximación al hessiano. ${\displaystyle \epsilon >0}$ ${\displaystyle ||\nabla f(\mathbf {x} _{k})||\leq \epsilon .}$ ${\displaystyle B_{0}}$ ${\displaystyle B_{0}=I}$ ${\displaystyle B_{k}}$

El primer paso del algoritmo se lleva a cabo utilizando la inversa de la matriz , que se puede obtener de manera eficiente aplicando la fórmula de Sherman-Morrison al paso 5 del algoritmo, dando ${\displaystyle B_{k}}$

${\displaystyle B_{k+1}^{-1}=\left(I-{\frac {\mathbf {s} _{k}\mathbf {y} _{k}^{T}}{\mathbf {y} _{k}^{T}\mathbf {s} _{k}}}\right)B_{k}^{-1}\left(I-{\frac {\mathbf {y} _{k}\mathbf {s} _{k}^{T}}{\mathbf {y} _{k}^{T}\mathbf {s} _{k}}}\right)+{\frac {\mathbf {s} _{k}\mathbf {s} _{k}^{T}}{\mathbf {y} _{k}^{T}\mathbf {s} _{k}}}.}$

Esto se puede calcular de manera eficiente sin matrices temporales, reconociendo que es simétrico y que y son escalares, utilizando una expansión como ${\displaystyle B_{k}^{-1}}$ ${\displaystyle \mathbf {y} _{k}^{\mathrm {T} }B_{k}^{-1}\mathbf {y} _{k}}$ ${\displaystyle \mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k}}$

${\displaystyle B_{k+1}^{-1}=B_{k}^{-1}+{\frac {(\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k}+\mathbf {y} _{k}^{\mathrm {T} }B_{k}^{-1}\mathbf {y} _{k})(\mathbf {s} _{k}\mathbf {s} _{k}^{\mathrm {T} })}{(\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k})^{2}}}-{\frac {B_{k}^{-1}\mathbf {y} _{k}\mathbf {s} _{k}^{\mathrm {T} }+\mathbf {s} _{k}\mathbf {y} _{k}^{\mathrm {T} }B_{k}^{-1}}{\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k}}}.}$

Por lo tanto, para evitar cualquier inversión de la matriz, se puede aproximar la inversa de la hessiana en lugar de la hessiana misma: ^[9] ${\displaystyle H_{k}{\overset {\operatorname {def} }{=}}B_{k}^{-1}.}$

A partir de una estimación inicial y una matriz hessiana invertida aproximada , se repiten los siguientes pasos a medida que converge a la solución: ${\displaystyle \mathbf {x} _{0}}$ ${\displaystyle H_{0}}$ ${\displaystyle \mathbf {x} _{k}}$

1. Obtenga una dirección resolviendo . ${\displaystyle \mathbf {p} _{k}}$ ${\displaystyle \mathbf {p} _{k}=-H_{k}\nabla f(\mathbf {x} _{k})}$
2. Realice una optimización unidimensional ( búsqueda lineal ) para encontrar un tamaño de paso aceptable en la dirección encontrada en el primer paso. Si se realiza una búsqueda lineal exacta, entonces . En la práctica, una búsqueda lineal inexacta suele ser suficiente, con un tamaño de paso aceptable que satisfaga las condiciones de Wolfe . ${\displaystyle \alpha _{k}}$ ${\displaystyle \alpha _{k}=\arg \min f(\mathbf {x} _{k}+\alpha \mathbf {p} _{k})}$ ${\displaystyle \alpha _{k}}$
3. Establecer y actualizar . ${\displaystyle \mathbf {s} _{k}=\alpha _{k}\mathbf {p} _{k}}$ ${\displaystyle \mathbf {x} _{k+1}=\mathbf {x} _{k}+\mathbf {s} _{k}}$
4. ${\displaystyle \mathbf {y} _{k}={\nabla f(\mathbf {x} _{k+1})-\nabla f(\mathbf {x} _{k})}}$.
5. ${\displaystyle H_{k+1}=H_{k}+{\frac {(\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k}+\mathbf {y} _{k}^{\mathrm {T} }H_{k}\mathbf {y} _{k})(\mathbf {s} _{k}\mathbf {s} _{k}^{\mathrm {T} })}{(\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k})^{2}}}-{\frac {H_{k}\mathbf {y} _{k}\mathbf {s} _{k}^{\mathrm {T} }+\mathbf {s} _{k}\mathbf {y} _{k}^{\mathrm {T} }H_{k}}{\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k}}}}$.

En problemas de estimación estadística (como la máxima verosimilitud o la inferencia bayesiana), los intervalos creíbles o los intervalos de confianza para la solución se pueden estimar a partir de la inversa de la matriz hessiana final ^{[ cita requerida ]} . Sin embargo, estas cantidades están definidas técnicamente por la verdadera matriz hessiana, y la aproximación BFGS puede no converger a la verdadera matriz hessiana. ^[10]

Desarrollos futuros

La fórmula de actualización de BFGS depende en gran medida de que la curvatura sea estrictamente positiva y esté acotada a partir de cero. Esta condición se cumple cuando realizamos una búsqueda lineal con condiciones de Wolfe en un objetivo convexo. Sin embargo, algunas aplicaciones de la vida real (como los métodos de programación cuadrática secuencial) producen rutinariamente curvaturas negativas o casi nulas. Esto puede ocurrir cuando se optimiza un objetivo no convexo o cuando se emplea un enfoque de región de confianza en lugar de una búsqueda lineal. También es posible producir valores espurios debido al ruido en el objetivo. ${\displaystyle \mathbf {s} _{k}^{\top }\mathbf {y} _{k}}$

En tales casos, se puede utilizar una de las denominadas actualizaciones BFGS amortiguadas (ver ^[11] ) que modifican y/o para obtener una actualización más robusta. ${\displaystyle \mathbf {s} _{k}}$ ${\displaystyle \mathbf {y} _{k}}$

Implementaciones notables

Las implementaciones de código abierto notables son:

• ALGLIB implementa BFGS y su versión de memoria limitada en C++ y C#
• GNU Octave utiliza una forma de BFGS en su fsolvefunción, con extensiones de región de confianza .
• El GSL implementa BFGS como gsl_multimin_fdfminimizer_vector_bfgs2. ^[12]
• En R , el algoritmo BFGS (y la versión L-BFGS-B que permite restricciones de caja) se implementa como una opción de la función base optim(). ^[13]
• En SciPy , la función scipy.optimize.fmin_bfgs implementa BFGS. ^[14] También es posible ejecutar BFGS utilizando cualquiera de los algoritmos L-BFGS configurando el parámetro L en un número muy grande.
• En Julia , el paquete Optim.jl implementa BFGS y L-BFGS como una opción de solución para la función optimizar() (entre otras opciones). ^[15]

Las implementaciones propietarias notables incluyen:

• El software de optimización no lineal a gran escala Artelys Knitro implementa, entre otros, los algoritmos BFGS y L-BFGS.
• En la caja de herramientas de optimización de MATLAB , la función fminunc utiliza BFGS con búsqueda de línea cúbica cuando el tamaño del problema se establece en "escala media".
• Mathematica incluye BFGS.
• LS-DYNA también utiliza BFGS para resolver problemas implícitos.

Véase también

• Algoritmo BHHH
• Fórmula de Davidon-Fletcher-Powell
• Descenso de gradiente
• L-BFGS
• Algoritmo de Levenberg-Marquardt
• Método Nelder-Mead
• Búsqueda de patrones (optimización)
• Métodos cuasi-newtonianos
• Rango uno simétrico
• Representación compacta cuasi-Newton

Referencias

1. Fletcher, Roger (1987), Métodos prácticos de optimización (2.ª ed.), Nueva York: John Wiley & Sons , ISBN 978-0-471-91547-8
2. Dennis, JE Jr. ; Schnabel, Robert B. (1983), "Métodos secantes para minimización sin restricciones", Métodos numéricos para optimización sin restricciones y ecuaciones no lineales , Englewood Cliffs, NJ: Prentice-Hall, págs. 194-215, ISBN 0-13-627216-9
3. Byrd, Richard H.; Lu, Peihuang; Nocedal, Jorge; Zhu, Ciyou (1995), "Un algoritmo de memoria limitada para optimización con restricciones limitadas", SIAM Journal on Scientific Computing , 16 (5): 1190–1208, CiteSeerX 10.1.1.645.5814 , doi :10.1137/0916069 
4. Broyden, CG (1970), "La convergencia de una clase de algoritmos de minimización de doble rango", Revista del Instituto de Matemáticas y sus Aplicaciones , 6 : 76–90, doi :10.1093/imamat/6.1.76
5. Fletcher, R. (1970), "Un nuevo enfoque para los algoritmos de métrica variable", Computer Journal , 13 (3): 317–322, doi : 10.1093/comjnl/13.3.317
6. Goldfarb, D. (1970), "Una familia de actualizaciones de métricas variables derivadas de medias variacionales", Mathematics of Computation , 24 (109): 23–26, doi : 10.1090/S0025-5718-1970-0258249-6
7. Shanno, David F. (julio de 1970), "Condicionamiento de métodos cuasi-Newton para la minimización de funciones", Mathematics of Computation , 24 (111): 647–656, doi : 10.1090/S0025-5718-1970-0274029-X , MR  0274029
8. Fletcher, Roger (1987), Métodos prácticos de optimización (2.ª ed.), Nueva York: John Wiley & Sons , ISBN 978-0-471-91547-8
9. Nocedal, Jorge; Wright, Stephen J. (2006), Optimización numérica (2.ª ed.), Berlín, Nueva York: Springer-Verlag , ISBN 978-0-387-30303-1
10. Ge, Ren-pu; Powell, MJD (1983). "La convergencia de matrices métricas variables en optimización sin restricciones". Programación matemática . 27 (2). 123. doi :10.1007/BF02591941. S2CID  8113073.
11. Jorge Nocedal; Stephen J. Wright (2006), Optimización numérica
12. "Biblioteca científica GNU — documentación de GSL 2.6". www.gnu.org . Consultado el 22 de noviembre de 2020 .
13. "R: Optimización de propósito general". stat.ethz.ch . Consultado el 22 de noviembre de 2020 .
14. "scipy.optimize.fmin_bfgs — Guía de referencia de SciPy v1.5.4". docs.scipy.org . Consultado el 22 de noviembre de 2020 .
15. "Opciones configurables de Optim.jl". julianlsolvers .

Lectura adicional

• Avriel, Mordecai (2003), Programación no lineal: análisis y métodos , Dover Publishing, ISBN 978-0-486-43227-4
• Bonnans, J. Frédéric; Gilbert, J. Charles; Lemaréchal, Claude ; Sagastizábal, Claudia A. (2006), "Métodos newtonianos", Optimización numérica: aspectos teóricos y prácticos (segunda edición), Berlín: Springer, págs. 51–66, ISBN 3-540-35445-X
• Fletcher, Roger (1987), Métodos prácticos de optimización (2.ª ed.), Nueva York: John Wiley & Sons , ISBN 978-0-471-91547-8
• Luenberger, David G. ; Ye, Yinyu (2008), Programación lineal y no lineal , International Series in Operations Research & Management Science, vol. 116 (tercera edición), Nueva York: Springer, pp. xiv+546, ISBN 978-0-387-74502-2, Sr.  2423726
• Kelley, CT (1999), Métodos iterativos para optimización , Filadelfia: Sociedad de Matemáticas Industriales y Aplicadas, págs. 71–86, ISBN 0-89871-433-8