Algoritmo de Broyden-Fletcher-Goldfarb-Shanno

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 precondicionando 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 únicamente a partir de evaluaciones del gradiente (o evaluaciones aproximadas del gradiente) mediante un método de secante generalizado . ^{[ 2 ]}

Dado que las actualizaciones de la matriz de curvatura BFGS no requieren inversión de matriz , su complejidad computacional es solo${\displaystyle {\mathcal {O}}(n^{2})}$, en comparación con${\displaystyle {\mathcal {O}}(n^{3})}$en 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 un número muy grande 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 con grandes restricciones.

El algoritmo recibe su nombre de Charles George Broyden , Roger Fletcher , Donald Goldfarb y David Shanno . ^{[ 4 ] [ 5 ] [ 6 ] [ 7 ]} Es una instancia de un algoritmo más general de John Greenstadt. ^{[ 8 ]}

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

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

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}),}$

dónde${\displaystyle B_{k}}$es una aproximación a la matriz hessiana en${\displaystyle \mathbf {x} _{k}}$, que se actualiza iterativamente en cada etapa, y${\displaystyle \nabla f(\mathbf {x} _ {k})}$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${\displaystyle f(\mathbf {x} _{k}+\gamma \mathbf {p} _{k})}$sobre el escalar${\displaystyle \gamma >0.}$

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

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

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

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

que es la ecuación de la secante.

La condición de curvatura${\displaystyle \mathbf {s} _{k}^{\mathsf {T}}\mathbf {y} _{k}>0}$debería quedar satisfecho${\displaystyle B_{k+1}}$ser definida positiva, lo cual se puede verificar premultiplicando la ecuación de la secante con${\displaystyle \mathbf {s} _{k}^{\mathsf {T}}}$. Si la función no es fuertemente convexa , entonces la condición debe imponerse explícitamente, por ejemplo, encontrando un punto x _{k +1} que satisfaga las condiciones de Wolfe , que implican la condición de curvatura, utilizando una búsqueda lineal.

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

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

Ambos${\displaystyle U_{k}}$y${\displaystyle V_{k}}$son matrices simétricas de rango uno, pero su suma es una matriz de actualización de rango dos. Las matrices de actualización BFGS y DFP difieren 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 positividad definida . Para mantener la simetría y la positividad definida de${\displaystyle B_{k+1}}$, el formulario de actualización se puede elegir como${\displaystyle B_{k+1}=B_{k}+\alpha \mathbf {u} \mathbf {u} ^{\mathsf {T}}+\beta \mathbf {v} \mathbf {v} ^{\mathsf {T}}}$. Imponiendo la condición de secante,${\displaystyle B_{k+1}\mathbf {s} _{k}=\mathbf {y} _{k}}$. Elegir${\displaystyle \mathbf {u} =\mathbf {y} _ {k}}$y${\displaystyle \mathbf {v} =B_ {k}\mathbf {s} _ {k}}$, podemos obtener: ^{[ 9 ]}

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

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

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

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

Consideremos el siguiente problema de optimización sin restricciones. $${\displaystyle {\begin{aligned}{\underset {\mathbf {x} \in \mathbb {R} ^{n}}{\text{minimizar}}}\quad &f(\mathbf {x} ),\end{aligned}}}$$ dónde${\displaystyle f:\mathbb {R} ^{n}\to \mathbb {R} }$es una función objetivo no lineal y dos veces diferenciable .

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

1. Obtén una dirección${\displaystyle \mathbf {p} _{k}}$resolviendo${\displaystyle B_{k}\mathbf {p} _{k}=-\nabla f(\mathbf {x} _{k})}$.
2. Realizar una optimización unidimensional ( búsqueda lineal ) para encontrar un tamaño de paso aceptable.${\displaystyle \alpha _{k}}$en la dirección encontrada en el primer paso. Si se realiza una búsqueda lineal exacta, entonces${\displaystyle \alpha _{k}=\arg \min _{\alpha }f(\mathbf {x} _{k}+\alpha \mathbf {p} _{k})}$En la práctica, una búsqueda lineal inexacta suele ser suficiente, con un resultado aceptable.${\displaystyle \alpha _{k}}$que cumplen las condiciones de Wolfe .
3. Colocar${\displaystyle \mathbf {s} _{k}=\alpha _{k}\mathbf {p} _{k}}$y actualización${\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}^{\mathsf {T}}}{\mathbf {y} _{k}^{\mathsf {T}}\mathbf {s} _{k}}}-{\frac {B_{k}\mathbf {s} _{k}\mathbf {s} _{k}^{\mathsf {T}}B_{k}^{\mathsf {T}}}{\mathbf {s} _{k}^{\mathsf {T}}B_{k}\mathbf {s} _{k}}}}$.

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

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

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

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

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

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

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

1. Obtén una dirección${\displaystyle \mathbf {p} _{k}}$resolviendo${\displaystyle \mathbf {p} _{k}=-H_{k}\nabla f(\mathbf {x} _{k})}$.
2. Realizar una optimización unidimensional ( búsqueda lineal ) para encontrar un tamaño de paso aceptable.${\displaystyle \alpha _{k}}$en la dirección encontrada en el primer paso. Si se realiza una búsqueda lineal exacta, entonces${\displaystyle \alpha _{k}=\arg \min _{\alpha }f(\mathbf {x} _{k}+\alpha \mathbf {p} _{k})}$En la práctica, una búsqueda lineal inexacta suele ser suficiente, con un resultado aceptable.${\displaystyle \alpha _{k}}$que cumplen las condiciones de Wolfe .
3. Colocar${\displaystyle \mathbf {s} _{k}=\alpha _{k}\mathbf {p} _{k}}$y actualización${\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}^{\mathsf {T}}\mathbf {y} _{k}+\mathbf {y} _{k}^{\mathsf {T}}H_{k}\mathbf {y} _{k})(\mathbf {s} _{k}\mathbf {s} _{k}^{\mathsf {T}})}{(\mathbf {s} _{k}^{\mathsf {T}}\mathbf {y} _{k})^{2}}}-{\frac {H_{k}\mathbf {y} _{k}\mathbf {s} _{k}^{\mathsf {T}}+\mathbf {s} _{k}\mathbf {y} _{k}^{\mathsf {T}}H_{k}}{\mathbf {s} _{k}^{\mathsf {T}}\mathbf {y} _{k}}}}$.

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

La fórmula de actualización de BFGS depende en gran medida de la curvatura.${\displaystyle \mathbf {s} _{k}^{\mathsf {T}}\mathbf {y} _{k}}$Ser estrictamente positivo y acotado lejos de cero. Esta condición se cumple al realizar una búsqueda lineal con condiciones de Wolfe en un objetivo convexo. Sin embargo, algunas aplicaciones prácticas (como los métodos de programación cuadrática secuencial) suelen producir 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 lineal. También es posible obtener valores espurios debido al ruido en el objetivo.

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

Algunas implementaciones de código abierto destacadas son:

• ALGLIB implementa BFGS y su versión de memoria limitada en C++ y C#.
• GNU Octave utiliza una forma de BFGS en su fsolvefuncionamiento, con extensiones de región de confianza .
• GSL implementa BFGS como gsl_multimin_fdfminimizer_vector_bfgs2 . ^{[ 13 ]}
• 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(). ^{[ 14 ]}
• En SciPy , la función scipy.optimize.fmin_bfgs implementa BFGS. ^{[ 15 ]} También es posible ejecutar BFGS utilizando cualquiera de los algoritmos L-BFGS estableciendo el parámetro L en un número muy grande. Este es también uno de los métodos predeterminados que se utilizan al ejecutar scipy.optimize.minimize sin restricciones. ^{[ 16 ]}
• En Julia , el paquete Optim.jl implementa BFGS y L-BFGS como una opción de solucionador para la función optimize() (entre otras opciones). ^{[ 17 ]}
• Stan implementa BFGS junto con la diferenciación automática como una opción para resolver problemas de estimación de máxima verosimilitud y estimación de máxima probabilidad a posteriori .

Entre las implementaciones propietarias más destacadas se 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 lineal 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.

• 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, pp. 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 operativa y ciencias de la gestión, vol.  116 (Tercera  ed.), Nueva York: Springer, pp.  xiv+546, ISBN 978-0-387-74502-2, MR 2423726 
• Kelley, CT (1999), Métodos iterativos para la optimización , Filadelfia: Society for Industrial and Applied Mathematics, págs. 71–86 , ISBN  0-89871-433-8