Representación compacta cuasi-Newton

Descomposición matricial

La representación compacta para los métodos cuasi-Newton es una descomposición matricial , que se utiliza normalmente en algoritmos de optimización basados ​​en gradientes o para resolver sistemas no lineales . La descomposición utiliza una representación de bajo rango para el hessiano directo y/o inverso o el jacobiano de un sistema no lineal. Debido a esto, la representación compacta se utiliza a menudo para problemas grandes y optimización restringida .

La representación compacta (derecha) de una aproximación hessiana densa (izquierda) es una matriz inicial (normalmente diagonal) más una descomposición de rango bajo. Ocupa poco espacio en la memoria (áreas sombreadas) y permite realizar cálculos matriciales eficientes.

Definición

La representación compacta de una matriz cuasi-Newton para la hessiana inversa o hessiana directa de una función objetivo no lineal expresa una secuencia de actualizaciones recursivas de la matriz de rango 1 o rango 2 como una actualización de rango o rango de una matriz inicial. ^{[1] [2]} Debido a que se deriva de actualizaciones cuasi-Newton, utiliza diferencias de iteraciones y gradientes en su definición . En particular, para o las matrices rectangulares y los sistemas simétricos cuadrados dependen de la y definen las representaciones cuasi-Newton ${\displaystyle H_{k}}$ ${\displaystyle B_{k}}$ ${\displaystyle f(x):\mathbb {R} ^{n}\to \mathbb {R} }$ ${\displaystyle k}$ ${\displaystyle 2k}$ ${\displaystyle \nabla f(x_{k})=g_{k}}$ ${\displaystyle \{s_{i-1}=x_{i}-x_{i-1},y_{i-1}=g_{i}-g_{i-1}\}_{i=1}^{k}}$ ${\displaystyle r=k}$ ${\displaystyle r=2k}$ ${\displaystyle n\times r}$ ${\displaystyle U_{k},J_{k}}$ ${\displaystyle r\times r}$ ${\displaystyle M_{k},N_{k}}$ ${\displaystyle s_{i},y_{i}}$

${\displaystyle H_{k}=H_{0}+U_{k}M_{k}^{-1}U_{k}^{T},\quad {\text{ and }}\quad B_{k}=B_{0}+J_{k}N_{k}^{-1}J_{k}^{T}}$

Aplicaciones

Debido a la descomposición matricial especial, la representación compacta se implementa en software de optimización de última generación. ^{[3] [4] [5] [6]} Cuando se combina con técnicas de memoria limitada, es una técnica popular para la optimización restringida con gradientes. ^[7] Las operaciones de álgebra lineal se pueden realizar de manera eficiente, como productos matriz-vector , soluciones o descomposiciones propias . Se puede combinar con técnicas de búsqueda de línea y región de confianza , y la representación se ha desarrollado para muchas actualizaciones cuasi-Newton. Por ejemplo, el producto matriz-vector con el hessiano cuasi-Newton directo y un vector arbitrario es: ${\displaystyle g\in \mathbb {R} ^{n}}$

${\displaystyle {\begin{aligned}p_{k}^{(0)}&=J_{k}^{T}g\\{\text{solve}}\quad N_{k}p_{k}^{(1)}&=p_{k}^{(0)}\quad \quad {\text{(}}N_{k}{\text{ is small)}}\\p_{k}^{(2)}&=J_{k}p_{k}^{(1)}\\p_{k}^{(3)}&=H_{0}g\\p_{k}^{\phantom {(4)}}&=p_{k}^{(2)}+p_{k}^{(3)}\end{aligned}}}$

Fondo

En el contexto del método GMRES , Walker ^[8] demostró que un producto de transformaciones de Householder (una identidad más rango 1) se puede expresar como una fórmula matricial compacta. Esto condujo a la derivación de una expresión matricial explícita para el producto de matrices identidad más rango 1. ^[7] Específicamente, para y cuando el producto de actualizaciones de rango 1 a la identidad es La actualización BFGS se puede expresar en términos de productos de los , que tienen una fórmula matricial compacta. Por lo tanto, la recursión BFGS puede explotar estas representaciones de matriz de bloques ${\displaystyle k}$ ${\textstyle S_{k}={\begin{bmatrix}s_{0}&s_{1}&\ldots s_{k-1}\end{bmatrix}},}$ ${\displaystyle ~Y_{k}={\begin{bmatrix}y_{0}&y_{1}&\ldots y_{k-1}\end{bmatrix}},}$ ${\displaystyle ~(R_{k})_{ij}=s_{i-1}^{T}y_{j-1},}$ ${\displaystyle ~\rho _{i-1}=1/s_{i-1}^{T}y_{i-1}}$ ${\textstyle ~V_{i}=I-\rho _{i-1}y_{i-1}s_{i-1}^{T}}$ ${\displaystyle 1\leq i\leq j\leq k}$ ${\displaystyle k}$ $${\displaystyle \prod _{i=1}^{k}V_{i-1}=\left(I-\rho _{0}y_{0}s_{0}^{T}\right)\cdots \left(I-\rho _{k-1}y_{k-1}s_{k-1}^{T}\right)=I-Y_{k}R_{k}^{-1}S_{k}^{T}}$$ ${\displaystyle V_{i}}$

Actualizaciones recursivas cuasi-Newton

Una familia paramétrica de actualizaciones cuasi-Newton incluye muchas de las fórmulas más conocidas. ^[9] Para vectores arbitrarios y tales que y fórmulas de actualización recursivas generales para las estimaciones hessianas inversas y directas son ${\displaystyle v_{k}}$ ${\displaystyle c_{k}}$ ${\displaystyle v_{k}^{T}y_{k}\neq 0}$ ${\displaystyle c_{k}^{T}s_{k}\neq 0}$

Al realizar elecciones específicas para los vectores de parámetros y se recuperan métodos bien conocidos. ${\displaystyle v_{k}}$ ${\displaystyle c_{k}}$

Representaciones compactas

Recopilando los vectores de actualización de las fórmulas recursivas en matrices, defina

$${\displaystyle S_{k}={\begin{bmatrix}s_{0}&s_{1}&\ldots &s_{k-1}\end{bmatrix}},}$$ $${\displaystyle Y_{k}={\begin{bmatrix}y_{0}&y_{1}&\ldots &y_{k-1}\end{bmatrix}},}$$ $${\displaystyle V_{k}={\begin{bmatrix}v_{0}&v_{1}&\ldots &v_{k-1}\end{bmatrix}},}$$ $${\displaystyle C_{k}={\begin{bmatrix}c_{0}&c_{1}&\ldots &c_{k-1}\end{bmatrix}},}$$

triangular superior

$${\displaystyle {\big (}R_{k}{\big )}_{ij}:={\big (}R_{k}^{\text{SY}}{\big )}_{ij}=s_{i-1}^{T}y_{j-1},\quad {\big (}R_{k}^{\text{VY}}{\big )}_{ij}=v_{i-1}^{T}y_{j-1},\quad {\big (}R_{k}^{\text{CS}}{\big )}_{ij}=c_{i-1}^{T}s_{j-1},\quad \quad {\text{ for }}1\leq i\leq j\leq k}$$

triangular inferior

$${\displaystyle {\big (}L_{k}{\big )}_{ij}:={\big (}L_{k}^{\text{SY}}{\big )}_{ij}=s_{i-1}^{T}y_{j-1},\quad {\big (}L_{k}^{\text{VY}}{\big )}_{ij}=v_{i-1}^{T}y_{j-1},\quad {\big (}L_{k}^{\text{CS}}{\big )}_{ij}=c_{i-1}^{T}s_{j-1},\quad \quad {\text{ for }}1\leq j<i\leq k}$$

y diagonal

$${\displaystyle (D_{k})_{ij}:={\big (}D_{k}^{\text{SY}}{\big )}_{ij}=s_{i-1}^{T}y_{j-1},\quad \quad {\text{ for }}1\leq i=j\leq k}$$

Con estas definiciones, en Brust se han desarrollado las representaciones compactas de las actualizaciones generales de rango 2 en ( 2 ) y ( 3 ) (incluidas las conocidas actualizaciones cuasi-Newton de la Tabla 1): ^[11]

$${\displaystyle U_{k}={\begin{bmatrix}V_{k}&S_{k}-H_{0}Y_{k}\end{bmatrix}}}$$

$${\displaystyle M_{k}={\begin{bmatrix}0_{k\times k}&R_{k}^{\text{VY}}\\{\big (}R_{k}^{\text{VY}}{\big )}^{T}&R_{k}+R_{k}^{T}-(D_{k}+Y_{k}^{T}H_{0}Y_{k})\end{bmatrix}}}$$

y la fórmula para el hessiano directo es

$${\displaystyle J_{k}={\begin{bmatrix}C_{k}&Y_{k}-B_{0}S_{k}\end{bmatrix}}}$$

$${\displaystyle N_{k}={\begin{bmatrix}0_{k\times k}&R_{k}^{\text{CS}}\\{\big (}R_{k}^{\text{CS}}{\big )}^{T}&R_{k}+R_{k}^{T}-(D_{k}+S_{k}^{T}B_{0}S_{k})\end{bmatrix}}}$$

Por ejemplo, cuando la representación en ( 4 ) es la fórmula compacta para la recursión BFGS en ( 1 ). ${\displaystyle V_{k}=S_{k}}$

Representaciones específicas

Antes del desarrollo de las representaciones compactas de ( 2 ) y ( 3 ), se descubrieron representaciones equivalentes para la mayoría de las actualizaciones conocidas (ver Tabla 1).

BFGS

Junto con la representación SR1, la representación compacta BFGS (Broyden-Fletcher-Goldfarb-Shanno) fue la primera fórmula compacta conocida. ^[7] En particular, la representación inversa está dada por

${\displaystyle H_{k}=H_{0}+U_{k}M_{k}^{-1}U_{k}^{T},\quad U_{k}={\begin{bmatrix}S_{k}&H_{0}Y_{k}\end{bmatrix}},\quad M_{k}^{-1}=\left[{\begin{smallmatrix}R_{k}^{-T}(D_{k}+Y_{k}^{T}H_{0}Y_{k})R_{k}^{-1}&-R_{k}^{-T}\\-R_{k}^{-1}&0\end{smallmatrix}}\right]}$ La aproximación hessiana directa se puede encontrar aplicando la identidad de Sherman-Morrison-Woodbury a la hessiana inversa:

${\displaystyle B_{k}=B_{0}+J_{k}N_{k}^{-1}J_{k}^{T},\quad J_{k}={\begin{bmatrix}B_{0}S_{k}&Y_{k}\end{bmatrix}},\quad N_{k}=\left[{\begin{smallmatrix}S^{T}B_{0}S_{k}&L_{k}\\L_{k}^{T}&-D_{k}\end{smallmatrix}}\right]}$

SR1

La representación compacta SR1 (Rango simétrico 1) se propuso por primera vez en ^[7] . Utilizando las definiciones de y de arriba, la fórmula hessiana inversa se da por ${\displaystyle D_{k},L_{k}}$ ${\displaystyle R_{k}}$

${\displaystyle H_{k}=H_{0}+U_{k}M_{k}^{-1}U_{k}^{T},\quad U_{k}=S_{k}-H_{0}Y_{k},\quad M_{k}=R_{k}+R_{k}^{T}-D_{k}-Y_{k}^{T}H_{0}Y_{k}}$

La hessiana directa se obtiene mediante la identidad de Sherman-Morrison-Woodbury y tiene la forma

${\displaystyle B_{k}=B_{0}+J_{k}N_{k}^{-1}J_{k}^{T},\quad J_{k}=Y_{k}-B_{0}S_{k},\quad N_{k}=D_{k}+L_{k}+L_{k}^{T}-S_{k}^{T}B_{0}S_{k}}$

Manuscrito

El método de la secante simétrica multipunto (MSS) es un método que tiene como objetivo satisfacer múltiples ecuaciones secantes. La fórmula de actualización recursiva fue desarrollada originalmente por Burdakov. ^[12] La representación compacta para la hessiana directa se derivó en ^[13]

${\displaystyle B_{k}=B_{0}+J_{k}N_{k}^{-1}J_{k}^{T},\quad J_{k}={\begin{bmatrix}S_{k}&Y_{k}-B_{0}S_{k}\end{bmatrix}},\quad N_{k}=\left[{\begin{smallmatrix}W_{k}(S_{k}^{T}B_{0}S_{k}-(R_{k}-D_{k}+R_{k}^{T}))W_{k}&W_{k}\\W_{k}&0\end{smallmatrix}}\right]^{-1},\quad W_{k}=(S_{k}^{T}S_{k})^{-1}}$

Otra representación compacta equivalente para la matriz MSS se deriva reescribiéndola en términos de . ^[14] La representación inversa se puede obtener mediante la aplicación de la identidad de Sherman-Morrison-Woodbury. ${\displaystyle J_{k}}$ ${\displaystyle J_{k}={\begin{bmatrix}S_{k}&B_{0}Y_{k}\end{bmatrix}}}$

DFP

Dado que la actualización de DFP (Davidon Fletcher Powell) es dual de la fórmula BFGS (es decir, intercambiando y en la actualización de BFGS ), la representación compacta para DFP se puede obtener inmediatamente a partir de la de BFGS. ^[15] ${\displaystyle H_{k}\leftrightarrow B_{k}}$ ${\displaystyle H_{0}\leftrightarrow B_{0}}$ ${\displaystyle y_{k}\leftrightarrow s_{k}}$

BSP

La representación compacta PSB (Powell-Symmetric-Broyden) fue desarrollada para la aproximación hessiana directa. ^[16] Es equivalente a sustituir en ( 5 ) ${\displaystyle C_{k}=S_{k}}$

${\displaystyle B_{k}=B_{0}+J_{k}N_{k}^{-1}J_{k}^{T},\quad J_{k}={\begin{bmatrix}S_{k}&Y_{k}-B_{0}S_{k}\end{bmatrix}},\quad N_{k}=\left[{\begin{smallmatrix}0&R_{k}^{\text{SS}}\\(R_{k}^{\text{SS}})^{T}&R_{k}+R_{k}^{T}-(D_{k}+S_{k}^{T}B_{0}S_{k})\end{smallmatrix}}\right]}$

BFGS estructurado

Para los problemas de optimización estructurada en los que la función objetivo se puede descomponer en dos partes , donde se conocen los gradientes y el hessiano de pero solo se conoce el gradiente de, existen fórmulas BFGS estructuradas. La representación compacta de estos métodos tiene la forma general de ( 5 ), con y específicos . ^[17] ${\displaystyle f(x)={\widehat {k}}(x)+{\widehat {u}}(x)}$ ${\displaystyle {\widehat {k}}(x)}$ ${\displaystyle {\widehat {u}}(x)}$ ${\displaystyle J_{k}}$ ${\displaystyle N_{k}}$

BFGS reducido

La representación compacta reducida (RCR) de BFGS es para optimización lineal con restricciones de igualdad , donde está subdeterminado . Además de las matrices, la RCR también almacena las proyecciones de los sobre el espacio nulo de ${\displaystyle {\text{ minimize }}f(x){\text{ subject to: }}Ax=b}$ ${\displaystyle A}$ ${\displaystyle S_{k},Y_{k}}$ ${\displaystyle y_{i}}$ ${\displaystyle A}$

$${\displaystyle Z_{k}={\begin{bmatrix}z_{0}&z_{1}&\cdots z_{k-1}\end{bmatrix}},\quad z_{i}=Py_{i},\quad P=I-A(A^{T}A)^{-1}A^{T},\quad 0\leq i\leq k-1}$$

Para la representación compacta de la matriz BFGS (con un múltiplo de la identidad ) el bloque (1,1) de la matriz KKT inversa tiene la representación compacta ^[18] ${\displaystyle B_{k}}$ ${\displaystyle B_{0}}$

${\displaystyle K_{k}={\begin{bmatrix}B_{k}&A^{T}\\A&0\end{bmatrix}},\quad B_{0}={\frac {1}{\gamma _{k}}}I,\quad H_{0}=\gamma _{k}I,\quad \gamma _{k}>0}$

${\displaystyle {\big (}K_{k}^{-1}{\big )}_{11}=H_{0}+U_{k}M_{k}^{-1}U_{k}^{T},\quad U_{k}={\begin{bmatrix}A^{T}&S_{k}&Z_{k}\end{bmatrix}},\quad M_{k}=\left[{\begin{smallmatrix}-AA^{T}/\gamma _{k}&\\&G_{k}\end{smallmatrix}}\right],\quad G_{k}=\left[{\begin{smallmatrix}R_{k}^{-T}(D_{k}+Y_{k}^{T}H_{0}Y_{k})R_{k}^{-1}&-H_{0}R_{k}^{-T}\\-H_{0}R_{k}^{-1}&0\end{smallmatrix}}\right]^{-1}}$

Memoria limitada

Patrón en una estrategia de actualización de memoria limitada. Con un parámetro de memoria de , las primeras iteraciones llenan la matriz (por ejemplo, un triángulo superior para la representación compacta, ). Para las técnicas de memoria limitada se descarta la información más antigua y se agrega una nueva columna al final. ${\displaystyle m=7}$ ${\displaystyle k\leq m}$ ${\displaystyle R_{k}={\text{triu}}(S_{k}^{T}Y_{k})}$ ${\displaystyle k>m}$

El uso más común de las representaciones compactas es para la configuración de memoria limitada donde denota el parámetro de memoria, con valores típicos alrededor de (ver, por ejemplo, ^{[18] [7]} ). Entonces, en lugar de almacenar el historial de todos los vectores, uno limita esto a los vectores más recientes y posiblemente o . Además, típicamente la inicialización se elige como un múltiplo adaptativo de la identidad , con y . Los métodos de memoria limitada se utilizan con frecuencia para problemas a gran escala con muchas variables (es decir, pueden ser grandes), en los que las matrices de memoria limitada y (y posiblemente ) son altas y muy delgadas: y . ${\displaystyle m\ll n}$ ${\displaystyle m\in [5,12]}$ ${\displaystyle m}$ ${\displaystyle \{(s_{i},y_{i}\}_{i=k-m}^{k-1}}$ ${\displaystyle \{v_{i}\}_{i=k-m}^{k-1}}$ ${\displaystyle \{c_{i}\}_{i=k-m}^{k-1}}$ ${\displaystyle H_{k}^{(0)}=\gamma _{k}I}$ ${\displaystyle \gamma _{k}=y_{k-1}^{T}s_{k-1}/y_{k-1}^{T}y_{k-1}}$ ${\displaystyle B_{k}^{(0)}={\frac {1}{\gamma _{k}}}I}$ ${\displaystyle n}$ ${\displaystyle S_{k}\in \mathbb {R} ^{n\times m}}$ ${\displaystyle Y_{k}\in \mathbb {R} ^{n\times m}}$ ${\displaystyle V_{k},C_{k}}$ ${\displaystyle S_{k}={\begin{bmatrix}s_{k-l-1}&\ldots &s_{k-1}\end{bmatrix}}}$ ${\displaystyle Y_{k}={\begin{bmatrix}y_{k-l-1}&\ldots &y_{k-1}\end{bmatrix}}}$

Implementaciones

Las implementaciones de código abierto incluyen:

• El algoritmo 1030 de ACM TOMS implementa un solucionador L-SR1 ^{[19] [20]}
• optimLa rutina optimizadora de propósito general de R utiliza el método L-BFGS-B.
• El método de minimización del módulo de optimización de SciPy también incluye una opción para utilizar L-BFGS-B.
• IPOPT con información de primer orden

Las implementaciones que no son de código abierto incluyen:

• Los solucionadores de programación no lineal (PLN) de Artelys Knitro utilizan matrices cuasi-Newton compactas ^[3]
• L-BFGS-B (algoritmo ACM TOMS 778) ^[21]

Obras citadas

1. Nocedal, J.; Wright, SJ (2006). Optimización numérica. Springer Series in Operations Research and Financial Engineering. Springer Nueva York, NY. doi :10.1007/978-0-387-40065-5. ISBN 978-0-387-30303-1.
2. Brust, JJ (2018). Métodos de regiones de confianza cuasi-Newton a gran escala: solucionadores de alta precisión, inicializaciones densas y extensiones (tesis doctoral). Universidad de California, Merced.
3. Byrd, RH; Nocedal, J; Waltz, RA (2006). "KNITRO: Un paquete integrado para optimización no lineal". Optimización no lineal a gran escala . Optimización no convexa y sus aplicaciones. Vol. 83. En: Di Pillo, G., Roma, M. (eds) Optimización no lineal a gran escala. Optimización no convexa y sus aplicaciones, vol. 83.: Springer, Boston, MA. págs. 35-59. doi :10.1007/0-387-30065-1_4. ISBN 978-0-387-30063-4.{{cite book}}: CS1 maint: location (link)
4. Zhu, C.; Byrd, RH; Lu, P.; Nocedal, J. (1997). "Algoritmo 778: L-BFGS-B: subrutinas Fortran para optimización a gran escala con restricciones de límites". ACM Transactions on Mathematical Software (TOMS) . 23 (4): 550-560. doi :10.1145/279232.279236.
5. Brust, J.; Burdakov, O.; Erway, J.; Marcia, R. (2022). "Algoritmo 1030: SC-SR1: software MATLAB para métodos de región de confianza SR1 de memoria limitada". ACM Transactions on Mathematical Software (TOMS) . 48 (4): 1-33. doi :10.1145/3550269.
6. Wächter, A.; Biegler, LT (2006). "Sobre la implementación de un algoritmo de búsqueda de línea con filtro de punto interior para programación no lineal a gran escala". Programación matemática . 106 : 25-57. doi :10.1007/s10107-004-0559-y.
7. Byrd, RH; Nocedal, J.; Schnabel, RB (1994). "Representaciones de matrices cuasi-Newton y su uso en métodos de memoria limitada". Programación matemática . 63 (4): 129–156. doi :10.1007/BF01582063. S2CID  5581219.
8. Walker, HF (1988). "Implementación del método GMRES utilizando transformaciones de Householder". Revista SIAM sobre computación científica y estadística . 9 (1): 152–163. doi :10.1137/0909010.
9. Dennis, Jr, JE; Moré, JJ (1977). "Métodos, motivación y teoría quasi-newtoniana". SIAM Review . 19 (1): 46-89. doi :10.1137/1019005. hdl :1813/6056.{{cite journal}}: CS1 maint: multiple names: authors list (link)
10. ${\displaystyle S_{k+1}={\begin{bmatrix}s_{0}&\ldots &s_{k}\end{bmatrix}},~P_{k}^{\text{S}}=I-S_{k+1}(S_{k+1}^{T}S_{k+1})^{-1}S_{k+1}^{T}}$
11. Brust, JJ (2024). "Representaciones compactas útiles para el ajuste de datos". arXiv : 2403.12206 [math.OC].
12. Burdakov, OP (1983). "Métodos del tipo secante para sistemas de ecuaciones con matriz jacobiana simétrica". Análisis funcional numérico y optimización . 6 (2): 1–18. doi :10.1080/01630568308816160.
13. Burdakov, OP; Martínez, JM; Pilotta, EA (2002). "Un método secante simétrico multipunto de memoria limitada para optimización con restricciones de límites". Annals of Operations Research . 117 : 51–70. doi :10.1023/A:1021561204463.
14. Brust, JJ; Erway, JB; Marcia, RF (2024). "Métodos de regiones de confianza que cambian de forma utilizando matrices secantes simétricas multipunto". Optimization Methods and Software : 1–18. arXiv : 2209.12057 . doi :10.1080/10556788.2023.2296441.
15. Erway, JB; Jain, V.; Marcia, RF (2013). Sistemas DFP de memoria limitada desplazada . En 2013, Conferencia Asilomar sobre señales, sistemas y computadoras. IEEE. págs. 1033–1037.
16. Kanzow, C.; Steck, D. (2023). "Regularización de métodos cuasi-Newton de memoria limitada para minimización no convexa a gran escala". Computación de programación matemática . 15 (3): 417–444. doi :10.1007/s12532-023-00238-4.
17. Brust, J. J; Di, Z.; Leyffer, S.; Petra, CG (2021). "Representaciones compactas de matrices BFGS estructuradas". Optimización computacional y aplicaciones . 80 (1): 55–88. doi :10.1007/s10589-021-00297-0.
18. Brust, J. J; Marcia, RF; Petra, CG; Saunders, MA (2022). "Optimización a gran escala con restricciones de igualdad lineal utilizando representación compacta reducida". Revista SIAM de Computación Científica . 44 (1): A103–A127. arXiv : 2101.11048 . Código Bibliográfico :2022SJSC...44A.103B. doi :10.1137/21M1393819.
19. "Algoritmos recopilados del ACM (CALGO)". calgo.acm.org .
20. "TOMS Alg. 1030". calgo.acm.org/1030.zip .
21. Zhu, C.; Byrd, Richard H.; Lu, Peihuang; Nocedal, Jorge (1997). "L-BFGS-B: Algoritmo 778: L-BFGS-B, rutinas FORTRAN para optimización restringida a gran escala". ACM Transactions on Mathematical Software . 23 (4): 550–560. doi : 10.1145/279232.279236 . S2CID  207228122.