Compact quasi-Newton representation

The compact representation for quasi-Newton methods is a matrix decomposition, which is typically used in gradient based optimization algorithms or for solving nonlinear systems. The decomposition uses a low-rank representation for the direct and/or inverse Hessian or the Jacobian of a nonlinear system. Because of this, the compact representation is often used for large problems and constrained optimization.

The compact representation (right) of a dense Hessian approximation (left) is an initial matrix (typically diagonal) plus low rank decomposition. It has a small memory footprint (shaded areas) and enables efficient matrix computations

Definition

The compact representation of a quasi-Newton matrix for the inverse Hessian ${\displaystyle H_{k}}$ or direct Hessian ${\displaystyle B_{k}}$ of a nonlinear objective function ${\displaystyle f(x):\mathbb {R} ^{n}\to \mathbb {R} }$ expresses a sequence of recursive rank-1 or rank-2 matrix updates as one rank-${\displaystyle k}$ or rank-${\displaystyle 2k}$ update of an initial matrix. ^{[1] [2]} Because it is derived from quasi-Newton updates, it uses differences of iterates and gradients ${\displaystyle \nabla f(x_{k})=g_{k}}$ in its definition ${\displaystyle \{s_{i-1}=x_{i}-x_{i-1},y_{i-1}=g_{i}-g_{i-1}\}_{i=1}^{k}}$. In particular, for ${\displaystyle r=k}$ or ${\displaystyle r=2k}$ the rectangular ${\displaystyle n\times r}$ matrices ${\displaystyle U_{k},J_{k}}$ and the ${\displaystyle r\times r}$ square symmetric systems ${\displaystyle M_{k},N_{k}}$ depend on the ${\displaystyle s_{i},y_{i}}$'s and define the quasi-Newton representations

${\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}}$

Applications

Because of the special matrix decomposition the compact representation is implemented in state-of-the-art optimization software. ^{[3] [4] [5] [6]} When combined with limited-memory techniques it is a popular technique for constrained optimization with gradients. ^[7] Linear algebra operations can be done efficiently, like matrix-vector products, solves or eigendecompositions. It can be combined with line-search and trust region techniques, and the representation has been developed for many quasi-Newton updates. For instance, the matrix vector product with the direct quasi-Newton Hessian and an arbitrary vector ${\displaystyle g\in \mathbb {R} ^{n}}$ is:

${\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}}}$

Background

In the context of the GMRES method, Walker ^[8] showed that a product of Householder transformations (an identity plus rank-1) can be expressed as a compact matrix formula. This led to the derivation of an explicit matrix expression for the product of ${\displaystyle k}$ identity plus rank-1 matrices. ^[7] Specifically, for ${\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}}$ and ${\textstyle ~V_{i}=I-\rho _{i-1}y_{i-1}s_{i-1}^{T}}$ when ${\displaystyle 1\leq i\leq j\leq k}$ the product of ${\displaystyle k}$ rank-1 updates to the identity is $${\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}}$$ The BFGS update can be expressed in terms of products of the ${\displaystyle V_{i}}$'s, which have a compact matrix formula. Therefore, the BFGS recursion can exploit these block matrix representations

$${\displaystyle {\begin{aligned}H_{k}&=V_{k-1}H_{k-1}V_{k-1}^{T}+\rho _{k-1}s_{k-1}s_{k-1}^{T}\\&=\left(V_{k-1}\cdots V_{1}V_{0})H_{0}(V_{0}^{T}V_{1}^{T}\cdots V_{k-1}^{T}\right)+\\&{\phantom {=}}\rho _{0}\left(V_{k-1}\cdots V_{1}\right)s_{0}s_{0}^{T}\left(V_{1}^{T}\cdots V_{k-1}^{T}\right)+\\&{\phantom {=}}\quad \vdots \\&{\phantom {=}}\rho _{k-2}V_{k-1}s_{k-2}s_{k-2}^{T}V_{k-1}^{T}+\\&{\phantom {=}}\rho _{k-1}s_{k-1}s_{k-1}^{T}\end{aligned}}}$$

1

Recursive quasi-Newton updates

A parametric family of quasi-Newton updates includes many of the most known formulas. ^[9] For arbitrary vectors ${\displaystyle v_{k}}$ and ${\displaystyle c_{k}}$ such that ${\displaystyle v_{k}^{T}y_{k}\neq 0}$ and ${\displaystyle c_{k}^{T}s_{k}\neq 0}$ general recursive update formulas for the inverse and direct Hessian estimates are

$${\displaystyle H_{k+1}=H_{k}+{\frac {(s_{k}-H_{k}y_{k})v_{k}^{T}+v_{k}(s_{k}-H_{k}y_{k})^{T}}{v_{k}^{T}y_{k}}}-{\frac {(s_{k}-H_{k}y_{k})^{T}y_{k}}{(v_{k}^{T}y_{k})^{2}}}v_{k}v_{k}^{T}}$$

2

$${\displaystyle B_{k+1}=B_{k}+{\frac {(y_{k}-B_{k}s_{k})c_{k}^{T}+c_{k}(y_{k}-B_{k}s_{k})^{T}}{c_{k}^{T}s_{k}}}-{\frac {(y_{k}-B_{k}s_{k})^{T}s_{k}}{(c_{k}^{T}s_{k})^{2}}}c_{k}c_{k}^{T}}$$

3

By making specific choices for the parameter vectors ${\displaystyle v_{k}}$ and ${\displaystyle c_{k}}$ well known methods are recovered

Table 1: Quasi-Newton updates parametrized by vectors ${\displaystyle v_{k}}$ and ${\displaystyle c_{k}}$

${\displaystyle v_{k}}$	${\displaystyle {\text{method}}}$	${\displaystyle c_{k}}$	${\displaystyle {\text{method}}}$
${\displaystyle s_{k}}$	BFGS	${\displaystyle s_{k}}$	PSB (Powell Symmetric Broyden)
${\displaystyle y_{k}}$	${\displaystyle {\text{Greenstadt's}}}$	${\displaystyle y_{k}}$	DFP
${\displaystyle s_{k}-H_{k}y_{k}}$	SR1	${\displaystyle y_{k}-B_{k}s_{k}}$	SR1
		${\displaystyle P_{k}^{\text{S}}s_{k}}$ [10]	MSS (Multipoint-Symmetric-Secant)

Compact Representations

Collecting the updating vectors of the recursive formulas into matrices, define

$${\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}},}$$

upper triangular

$${\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}$$

lower triangular

$${\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}$$

and 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}$$

With these definitions the compact representations of general rank-2 updates in ( 2 ) and ( 3 ) (including the well known quasi-Newton updates in Table 1) have been developed in Brust: ^[11]

$${\displaystyle H_{k}=H_{0}+U_{k}M_{k}^{-1}U_{k}^{T},}$$

4

$${\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}}}$$

and the formula for the direct Hessian is

$${\displaystyle B_{k}=B_{0}+J_{k}N_{k}^{-1}J_{k}^{T},}$$

5

$${\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}}}$$

For instance, when ${\displaystyle V_{k}=S_{k}}$ the representation in ( 4 ) is the compact formula for the BFGS recursion in ( 1 ).

Specific Representations

Prior to the development of the compact representations of ( 2 ) and ( 3 ), equivalent representations have been discovered for most known updates (see Table 1).

BFGS

Along with the SR1 representation, the BFGS (Broyden-Fletcher-Goldfarb-Shanno) compact representation was the first compact formula known. ^[7] In particular, the inverse representation is given by

${\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]}$ The direct Hessian approximation can be found by applying the Sherman-Morrison-Woodbury identity to the inverse Hessian:

${\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

The SR1 (Symmetric Rank-1) compact representation was first proposed in. ^[7] Using the definitions of ${\displaystyle D_{k},L_{k}}$ and ${\displaystyle R_{k}}$ from above, the inverse Hessian formula is given by

${\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}}$

The direct Hessian is obtained by the Sherman-Morrison-Woodbury identity and has the form

${\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}}$

MSS

The multipoint symmetric secant (MSS) method is a method that aims to satisfy multiple secant equations. The recursive update formula was originally developed by Burdakov. ^[12] The compact representation for the direct Hessian was derived in ^[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}}$

Another equivalent compact representation for the MSS matrix is derived by rewriting ${\displaystyle J_{k}}$ in terms of ${\displaystyle J_{k}={\begin{bmatrix}S_{k}&B_{0}Y_{k}\end{bmatrix}}}$. ^[14] The inverse representation can be obtained by application for the Sherman-Morrison-Woodbury identity.

DFP

Since the DFP (Davidon Fletcher Powell) update is the dual of the BFGS formula (i.e., swapping ${\displaystyle H_{k}\leftrightarrow B_{k}}$, ${\displaystyle H_{0}\leftrightarrow B_{0}}$ and ${\displaystyle y_{k}\leftrightarrow s_{k}}$ in the BFGS update), the compact representation for DFP can be immediately obtained from the one for BFGS. ^[15]

PSB

The PSB (Powell-Symmetric-Broyden) compact representation was developed for the direct Hessian approximation. ^[16] It is equivalent to substituting ${\displaystyle C_{k}=S_{k}}$ in ( 5 )

${\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]}$

Structured BFGS

For structured optimization problems in which the objective function can be decomposed into two parts ${\displaystyle f(x)={\widehat {k}}(x)+{\widehat {u}}(x)}$, where the gradients and Hessian of ${\displaystyle {\widehat {k}}(x)}$ are known but only the gradient of ${\displaystyle {\widehat {u}}(x)}$ is known, structured BFGS formulas exist. The compact representation of these methods has the general form of ( 5 ), with specific ${\displaystyle J_{k}}$ and ${\displaystyle N_{k}}$. ^[17]

Reduced BFGS

The reduced compact representation (RCR) of BFGS is for linear equality constrained optimization ${\displaystyle {\text{ minimize }}f(x){\text{ subject to: }}Ax=b}$, where ${\displaystyle A}$ is underdetermined. In addition to the matrices ${\displaystyle S_{k},Y_{k}}$ the RCR also stores the projections of the ${\displaystyle y_{i}}$'s onto the nullspace of ${\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}$$

For ${\displaystyle B_{k}}$ the compact representation of the BFGS matrix (with a multiple of the identity ${\displaystyle B_{0}}$) the (1,1) block of the inverse KKT matrix has the compact representation ^[18]

${\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}}$

Limited Memory

Pattern in a limited-memory updating strategy. With a memory parameter of ${\displaystyle m=7}$, the first ${\displaystyle k\leq m}$ iterations fill the matrix (e.g., an upper triangular for the compact representation, ${\displaystyle R_{k}={\text{triu}}(S_{k}^{T}Y_{k})}$). For ${\displaystyle k>m}$ the limited-memory techniques discards the oldest information and add a new column to the end.

The most common use of the compact representations is for the limited-memory setting where ${\displaystyle m\ll n}$ denotes the memory parameter, with typical values around ${\displaystyle m\in [5,12]}$ (see e.g., ^{[18] [7]} ). Then, instead of storing the history of all vectors one limits this to the ${\displaystyle m}$ most recent vectors ${\displaystyle \{(s_{i},y_{i}\}_{i=k-m}^{k-1}}$ and possibly ${\displaystyle \{v_{i}\}_{i=k-m}^{k-1}}$ or ${\displaystyle \{c_{i}\}_{i=k-m}^{k-1}}$. Further, typically the initialization is chosen as an adaptive multiple of the identity ${\displaystyle H_{k}^{(0)}=\gamma _{k}I}$, with ${\displaystyle \gamma _{k}=y_{k-1}^{T}s_{k-1}/y_{k-1}^{T}y_{k-1}}$ and ${\displaystyle B_{k}^{(0)}={\frac {1}{\gamma _{k}}}I}$. Limited-memory methods are frequently used for large-scale problems with many variables (i.e., ${\displaystyle n}$ can be large), in which the limited-memory matrices ${\displaystyle S_{k}\in \mathbb {R} ^{n\times m}}$ and ${\displaystyle Y_{k}\in \mathbb {R} ^{n\times m}}$ (and possibly ${\displaystyle V_{k},C_{k}}$) are tall and very skinny: ${\displaystyle S_{k}={\begin{bmatrix}s_{k-l-1}&\ldots &s_{k-1}\end{bmatrix}}}$ and ${\displaystyle Y_{k}={\begin{bmatrix}y_{k-l-1}&\ldots &y_{k-1}\end{bmatrix}}}$.

Implementations

Open source implementations include:

• ACM TOMS algorithm 1030 implements a L-SR1 solver ^{[19] [20]}
• R's optim general-purpose optimizer routine uses the L-BFGS-B method.
• SciPy's optimization module's minimize method also includes an option to use L-BFGS-B.
• IPOPT with first order information

Non open source implementations include:

• Artelys Knitro nonlinear programming (NLP) solvers use compact quasi-Newton matrices ^[3]
• L-BFGS-B (ACM TOMS algorithm 778) ^[21]

Works cited

1. Nocedal, J.; Wright, S.J. (2006). Numerical Optimization. Springer Series in Operations Research and Financial Engineering. Springer New York, NY. doi:10.1007/978-0-387-40065-5. ISBN   978-0-387-30303-1.
2. Brust, J. J. (2018). Large-Scale Quasi-Newton Trust-Region Methods: High-Accuracy Solvers, Dense Initializations, and Extensions (PhD thesis). University of California, Merced.
3. Byrd, R. H.; Nocedal, J; Waltz, R. A. (2006). "KNITRO: An integrated package for nonlinear optimization". Large-Scale Nonlinear Optimization. Nonconvex Optimization and Its Applications. Vol. 83. In: Di Pillo, G., Roma, M. (eds) Large-Scale Nonlinear Optimization. Nonconvex Optimization and Its Applications, vol 83.: Springer, Boston, MA. pp. 35–59. doi:10.1007/0-387-30065-1_4. ISBN   978-0-387-30063-4.
4. Zhu, C.; Byrd, R. H.; Lu, P.; Nocedal, J. (1997). "Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization". ACM Transactions on Mathematical Software. 23 (4): 550–560. doi:10.1145/279232.279236.
5. Brust, J.; Burdakov, O.; Erway, J.; Marcia, R. (2022). "Algorithm 1030: SC-SR1: MATLAB software for limited-memory SR1 trust-region methods". ACM Transactions on Mathematical Software. 48 (4): 1–33. doi:10.1145/3550269.
6. Wächter, A.; Biegler, L. T. (2006). "On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming". Mathematical Programming. 106: 25–57. doi:10.1007/s10107-004-0559-y.
7. Byrd, R. H.; Nocedal, J.; Schnabel, R. B. (1994). "Representations of Quasi-Newton Matrices and their use in Limited Memory Methods". Mathematical Programming. 63 (4): 129–156. doi:10.1007/BF01582063. S2CID   5581219.
8. Walker, H. F. (1988). "Implementation of the GMRES Method Using Householder Transformations". SIAM Journal on Scientific and Statistical Computing. 9 (1): 152–163. doi:10.1137/0909010.
9. Dennis, Jr, J. E.; Moré, J. J. (1977). "Quasi-Newton methods, motivation and theory" (PDF). SIAM Review. 19 (1): 46–89. doi:10.1137/1019005. hdl: 1813/6056 .
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, J. J. (2024). "Useful Compact Representations for Data-Fitting". arXiv: 2403.12206 [math.OC].
12. Burdakov, O. P. (1983). "Methods of the secant type for systems of equations with symmetric Jacobian matrix". Numerical Functional Analysis and Optimization. 6 (2): 1–18. doi:10.1080/01630568308816160.
13. Burdakov, O. P.; Martínez, J. M.; Pilotta, E. A. (2002). "A limited-memory multipoint symmetric secant method for bound constrained optimization". Annals of Operations Research. 117 (1–4): 51–70. doi:10.1023/A:1021561204463.
14. Brust, J. J.; Erway, J. B.; Marcia, R. F. (2024). "Shape-changing trust-region methods using multipoint symmetric secant matrices". Optimization Methods and Software. 39 (5): 990–1007. arXiv: 2209.12057 . doi:10.1080/10556788.2023.2296441.
15. Erway, J. B.; Jain, V.; Marcia, R. F. (2013). Shifted limited-memory DFP systems. In 2013 Asilomar Conference on Signals, Systems and Computers. IEEE. pp. 1033–1037.
16. Kanzow, C.; Steck, D. (2023). "Regularization of limited memory quasi-Newton methods for large-scale nonconvex minimization". Mathematical Programming Computation. 15 (3): 417–444. arXiv: 1911.04584 . doi: 10.1007/s12532-023-00238-4 .
17. Brust, J. J; Di, Z.; Leyffer, S.; Petra, C. G. (2021). "Compact representations of structured BFGS matrices". Computational Optimization and Applications. 80 (1): 55–88. arXiv: 2208.00057 . doi:10.1007/s10589-021-00297-0.
18. Brust, J. J; Marcia, R.F.; Petra, C.G.; Saunders, M. A. (2022). "Large-scale optimization with linear equality constraints using reduced compact representation". SIAM Journal on Scientific Computing. 44 (1): A103 –A127. arXiv: 2101.11048 . Bibcode:2022SJSC...44A.103B. doi:10.1137/21M1393819.
19. "Collected Algorithms of the ACM". 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: Algorithm 778: L-BFGS-B, FORTRAN routines for large scale bound constrained optimization". ACM Transactions on Mathematical Software. 23 (4): 550–560. doi: 10.1145/279232.279236 . S2CID   207228122.