# 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.

## Definition

The compact representation of a quasi-Newton matrix for the inverse Hessian $H_{k}$ or direct Hessian $B_{k}$ of a nonlinear objective function $f(x):\mathbb {R} ^{n}\to \mathbb {R}$ expresses a sequence of recursive rank-1 or rank-2 matrix updates as one rank-$k$ or rank-$2k$ update of an initial matrix. Because it is derived from quasi-Newton updates, it uses differences of iterates and gradients $\nabla f(x_{k})=g_{k}$ in its definition $\{s_{i-1}=x_{i}-x_{i-1},y_{i-1}=g_{i}-g_{i-1}\}_{i=1}^{k}$. In particular, for $r=k$ or $r=2k$ the rectangular $n\times r$ matrices $U_{k},J_{k}$ and the $r\times r$ square symmetric systems $M_{k},N_{k}$ depend on the $s_{i},y_{i}$'s and define the quasi-Newton representations
    $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. When combined with limited-memory techniques it is a popular technique for constrained optimization with gradients. 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 $g\in \mathbb {R} ^{n}$ is:
    ${\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 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 $k$ identity plus rank-1 matrices. Specifically, for S$_k$=[s$_0$s$_1$…s$_{k−1}$], {\textstyle S_{k}={\begin{bmatrix}s_{0}&s_{1}&\ldots s_{k-1}\end{bmatrix}}, } $~Y_{k}={\begin{bmatrix}y_{0}&y_{1}&\ldots y_{k-1}\end{bmatrix}},$ $~(R_{k})_{ij}=s_{i-1}^{T}y_{j-1},$ $~\rho _{i-1}=1/s_{i-1}^{T}y_{i-1}$ and V$_i$=I−ρ$_{i−1}$y$_{i−1}$s$_{i−1}$$^T${\textstyle ~V_{i}=I-\rho _{i-1}y_{i-1}s_{i-1}^{T}} when $1\leq i\leq j\leq k$ the product of $k$ rank-1 updates to the identity is $\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 $V_{i}$'s, which have a compact matrix formula. Therefore, the BFGS recursion can exploit these block matrix representations

## Recursive quasi-Newton updates

A parametric family of quasi-Newton updates includes many of the most known formulas. For arbitrary vectors $v_{k}$ and $c_{k}$ such that $v_{k}^{T}y_{k}\neq 0$ and $c_{k}^{T}s_{k}\neq 0$ general recursive update formulas for the inverse and direct Hessian estimates are
By making specific choices for the parameter vectors $v_{k}$ and $c_{k}$ well known methods are recovered
Table 1: Quasi-Newton updates parametrized by vectors $v_{k}$ and $c_{k}$ :
| $v_{k}$ | ${\text{method}}$ | $c_{k}$ | ${\text{method}}$ |
| $s_{k}$ | BFGS | $s_{k}$ | PSB (Powell Symmetric Broyden) |
| $y_{k}$ | ${\text{Greenstadt's}}$ | $y_{k}$ | DFP |
| $s_{k}-H_{k}y_{k}$ | SR1 | $y_{k}-B_{k}s_{k}$ | SR1 |
| | | $P_{k}^{\text{S}}s_{k}$ | MSS (Multipoint-Symmetric-Secant) |

## Compact Representations

Collecting the updating vectors of the recursive formulas into matrices, define
$S_{k}={\begin{bmatrix}s_{0}&s_{1}&\ldots &s_{k-1}\end{bmatrix}},$ $Y_{k}={\begin{bmatrix}y_{0}&y_{1}&\ldots &y_{k-1}\end{bmatrix}},$ $V_{k}={\begin{bmatrix}v_{0}&v_{1}&\ldots &v_{k-1}\end{bmatrix}},$ $C_{k}={\begin{bmatrix}c_{0}&c_{1}&\ldots &c_{k-1}\end{bmatrix}},$
upper triangular
${\bigR_{k}{\big )}_{ij}:={\bigR_{k}^{\text{SY}}{\big )}_{ij}=s_{i-1}^{T}y_{j-1},\quad {\bigR_{k}^{\text{VY}}{\big )}_{ij}=v_{i-1}^{T}y_{j-1},\quad {\bigR_{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
${\bigL_{k}{\big )}_{ij}:={\bigL_{k}^{\text{SY}}{\big )}_{ij}=s_{i-1}^{T}y_{j-1},\quad {\bigL_{k}^{\text{VY}}{\big )}_{ij}=v_{i-1}^{T}y_{j-1},\quad {\bigL_{k}^{\text{CS}}{\big )}_{ij}=c_{i-1}^{T}s_{j-1},\quad \quad {\text{ for }}1\leq j<i\leq k$
and diagonal
$(D_{k})_{ij}:={\bigD_{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:
$U_{k}={\begin{bmatrix}V_{k}&S_{k}-H_{0}Y_{k}\end{bmatrix}}$
$M_{k}={\begin{bmatrix}0_{k\times k}&R_{k}^{\text{VY}}\\{\bigR_{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
$J_{k}={\begin{bmatrix}C_{k}&Y_{k}-B_{0}S_{k}\end{bmatrix}}$
$N_{k}={\begin{bmatrix}0_{k\times k}&R_{k}^{\text{CS}}\\{\bigR_{k}^{\text{CS}}{\big )}^{T}&R_{k}+R_{k}^{T}-(D_{k}+S_{k}^{T}B_{0}S_{k})\end{bmatrix}}$
For instance, when $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. In particular, the inverse representation is given by
$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:
$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. Using the definitions of $D_{k},L_{k}$ and $R_{k}$ from above, the inverse Hessian formula is given by
$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
$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. The compact representation for the direct Hessian was derived in
$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 $J_{k}$ in terms of $J_{k}={\begin{bmatrix}S_{k}&B_{0}Y_{k}\end{bmatrix}}$. 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 $H_{k}\leftrightarrow B_{k}$, $H_{0}\leftrightarrow B_{0}$ and $y_{k}\leftrightarrow s_{k}$ in the BFGS update), the compact representation for DFP can be immediately obtained from the one for BFGS.

#### PSB

The PSB (Powell-Symmetric-Broyden) compact representation was developed for the direct Hessian approximation. It is equivalent to substituting $C_{k}=S_{k}$ in (5)
$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 $f(x)={\widehat {k}}(x)+{\widehat {u}}(x)$, where the gradients and Hessian of ${\widehat {k}}(x)$ are known but only the gradient of ${\widehat {u}}(x)$ is known, structured BFGS formulas exist. The compact representation of these methods has the general form of (5), with specific $J_{k}$ and $N_{k}$.

#### Reduced BFGS

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

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

## Works cited

* ↑ 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.
* ↑ Brust, J. J. (2018). Large-Scale Quasi-Newton Trust-Region Methods: High-Accuracy Solvers, Dense Initializations, and Extensions (PhD thesis). University of California, Merced.
* 1 2 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.{{cite book}}: CS1 maint: location (link)
* ↑ 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.
* ↑ 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.
* ↑ 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.
* 1 2 3 4 5 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.
* ↑ 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.
* ↑ 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.{{cite journal}}: CS1 maint: multiple names: authors list (link)
* ↑ $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}$
* ↑ Brust, J. J. (2024). "Useful Compact Representations for Data-Fitting". arXiv:2403.12206 [math.OC].
* ↑ 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.
* ↑ 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: 51–70. doi:10.1023/A:1021561204463.
* ↑ 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.
* ↑ 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.
* ↑ 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.
* ↑ 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.
* 1 2 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.
* ↑ "Collected Algorithms of the ACM". calgo.acm.org.
* ↑ "TOMS Alg. 1030". calgo.acm.org/1030.zip.
* ↑ 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.