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 de 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 de Hesse de la función de pérdida , obtenida únicamente a partir de evaluaciones de gradiente (o evaluaciones de gradiente aproximadas) mediante un método secante generalizado . ^[2]

Dado que las actualizaciones de la matriz de curvatura BFGS no requieren inversión de la matriz , su complejidad computacional es única , en comparación con el método de Newton . También es de uso común L-BFGS , que es una versión de BFGS con memoria limitada que es particularmente adecuada para problemas con un gran número de variables (p. ej., >1000). La variante BFGS-B maneja restricciones de caja simples. ^[3] ${\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 es minimizar , donde es un vector en , y es una función escalar diferenciable. No hay restricciones sobre los valores que pueden 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 procede de forma iterativa 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 de Hesse 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 lineal en la dirección p _k para encontrar el siguiente punto x _{k +1} minimizando 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}).}$

Deja y , luego 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}}$,

que 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 aplicarse explícitamente, por ejemplo, encontrando un punto x _{k +1} que satisfaga las condiciones de Wolfe , que implican la condición de curvatura, utilizando la búsqueda lineal. ${\displaystyle \mathbf {s} _{k}^{\top }\mathbf {y} _{k}>0}$ ${\displaystyle B_{k+1}}$ ${\displaystyle \mathbf {s} _{k}^{T}}$

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

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

Ambas y son matrices simétricas de rango uno, pero su suma es una matriz de actualización de rango dos. La matriz de actualización de BFGS y DFP se diferencia 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 precisión positiva . Para mantener la simetría y la precisión positiva de , la forma de actualización se puede elegir 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 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

$${\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 de Hesse, se repiten los siguientes pasos hasta llegar 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 de líneas ) para encontrar un tamaño de paso aceptable en la dirección encontrada en el primer paso. Si se realiza una búsqueda de línea exacta, entonces . En la práctica, una búsqueda de línea inexacta suele ser suficiente, con condiciones de Wolfe aceptablemente satisfactorias . ${\displaystyle \alpha _{k}}$ ${\displaystyle \alpha _{k}=\arg \min f(\mathbf {x} _{k}+\alpha \mathbf {p} _{k})}$ ${\displaystyle \alpha _{k}}$
3. Configurar 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 algunos , uno puede detener el algoritmo cuando se inicializa con , el primer paso será equivalente a un descenso de 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, usando 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 matriz, se puede aproximar la inversa del hessiano en lugar del propio hessiano: ^[9] ${\displaystyle H_{k}{\overset {\operatorname {def} }{=}}B_{k}^{-1}.}$

A partir de una suposición inicial y una matriz de Hesse invertida aproximada , se repiten los siguientes pasos hasta llegar 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 de líneas ) para encontrar un tamaño de paso aceptable en la dirección encontrada en el primer paso. Si se realiza una búsqueda de línea exacta, entonces . En la práctica, una búsqueda de línea inexacta suele ser suficiente, con condiciones de Wolfe aceptablemente satisfactorias . ${\displaystyle \alpha _{k}}$ ${\displaystyle \alpha _{k}=\arg \min f(\mathbf {x} _{k}+\alpha \mathbf {p} _{k})}$ ${\displaystyle \alpha _{k}}$
3. Configurar 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 máxima verosimilitud o inferencia bayesiana), se pueden estimar intervalos creíbles o intervalos de confianza para la solución a partir de la inversa de la matriz de Hesse final ^{[ cita necesaria ]} . Sin embargo, estas cantidades están técnicamente definidas por la verdadera matriz de Hesse, y la aproximación BFGS puede no converger a la verdadera matriz de Hesse. ^[10]

Nuevos desarrollos

La fórmula de actualización de BFGS se basa en gran medida en que la curvatura sea estrictamente positiva y acotada desde 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 habitualmente curvaturas negativas o casi nulas. Esto puede ocurrir al optimizar un objetivo no convexo o al emplear un enfoque de región de confianza en lugar de una búsqueda de líneas. 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 llamadas 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 cuadro) 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 MATLAB Optimization Toolbox , 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.

Ver 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-Newton
• Rango uno simétrico

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, Nueva Jersey: 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 una optimización limitada", 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étricas variables", 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 por medios variacionales", Matemáticas de la Computación , 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", Matemáticas de la Computación , 24 (111): 647–656, doi : 10.1090/S0025-5718-1970-0274029-X , Señor  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 de 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 GSL 2.6". www.gnu.org . Consultado el 22 de noviembre de 2020 .
13. "R: Optimización de uso 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 .

Otras lecturas

• 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 ed.), 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 , Serie internacional en investigación de operaciones y ciencia de la gestión, vol. 116 (Tercera ed.), Nueva York: Springer, págs. xiv+546, ISBN 978-0-387-74502-2, señor  2423726
• Kelley, CT (1999), Métodos iterativos de optimización , Filadelfia: Sociedad de Matemáticas Industriales y Aplicadas, págs. 71–86, ISBN 0-89871-433-8