BFGS Formula from Kullback-Leibler Divergence

Question

On page 411 in this book, the authors give the following BFGS formula $$ \boxed{\boldsymbol C_{\textrm{BFGS}} = \boldsymbol C + \underbrace{\frac{\boldsymbol g^\top\boldsymbol\delta+\boldsymbol g^\top\boldsymbol C\boldsymbol g}{(\boldsymbol g^\top\boldsymbol\delta)^2}\boldsymbol\delta\boldsymbol\delta^\top-\frac{1}{\boldsymbol g^\top\boldsymbol\delta}\left(\boldsymbol\delta\boldsymbol g^\top\boldsymbol C + (\boldsymbol\delta\boldsymbol g^\top\boldsymbol C)^\top\right)}_{\textrm{BFGS update}}}\tag{1} $$ by considering the constrained optimisation problem: $$\begin{array}{rlcll} \min_{\boldsymbol A} & \mathcal D(\boldsymbol 0, \boldsymbol C \mid \boldsymbol 0, \boldsymbol A) \\ \textrm{subject to}: & \boldsymbol A\boldsymbol g = \boldsymbol\delta, \\ & \boldsymbol A ~~= \boldsymbol A^\top, \end{array}\tag{2}$$ where $$ \mathcal D(\boldsymbol 0, \boldsymbol C \mid \boldsymbol 0, \boldsymbol A) := \frac12\left(\mathrm{tr}\left(\boldsymbol A^{-1}\boldsymbol C\right) - \log\left(\det(\boldsymbol A^{-1}\boldsymbol C)\right) - n\right)$$ is the Kullback-Leibler divergence between the normal distributions $\mathcal N(\boldsymbol 0, \boldsymbol C)$ and $\mathcal N(\boldsymbol 0, \boldsymbol A)$.

Using Lagrange multipliers with the Lagrangian $$ \mathcal L(\boldsymbol A, \boldsymbol \beta) := \mathcal D(\boldsymbol 0, \boldsymbol C \mid \boldsymbol 0, \boldsymbol A) + \boldsymbol \beta^\top(\boldsymbol A\boldsymbol g - \boldsymbol \delta)$$ gives (also using symmetry of $\boldsymbol A$ and $\boldsymbol C$) $$ \frac{\partial}{\partial\boldsymbol A}\mathcal L(\boldsymbol A, \boldsymbol \beta) = \frac12\left(-\boldsymbol A^{-1}\boldsymbol C\boldsymbol A^{-1} + \boldsymbol A^{-1}\right) + \boldsymbol \beta\boldsymbol g^\top, \quad \boldsymbol{\nabla}_{\boldsymbol \beta}\mathcal L(\boldsymbol A, \boldsymbol \beta) = \boldsymbol A\boldsymbol g - \boldsymbol \delta,$$ and equating these both equal to zero gives $\boldsymbol A\boldsymbol g = \boldsymbol \delta$ and $\boldsymbol A^{-1}(\boldsymbol I - \boldsymbol C\boldsymbol A^{-1}) = -2\boldsymbol \beta \boldsymbol g^\top$. It is not obvious to me how this last expression could be used to identify the corresponding values of $\boldsymbol \beta$ and $\boldsymbol A$. I tried multiplying through by $\boldsymbol \delta$ to use the $\boldsymbol A\boldsymbol g = \boldsymbol \delta$ constraint, yielding $\boldsymbol{\delta}\boldsymbol \beta^\top\boldsymbol \delta = \frac12\left(\boldsymbol C\boldsymbol g - \boldsymbol \delta\right)$, but it is still unclear from this (1) how $\boldsymbol \beta$ is specified and (2) how $\boldsymbol{A}$ would be recovered.

What is the correct way to solve the constrained optimisation problem (2) to obtain the BFGS formula (1) with the method of Lagrange multipliers that gives values for $\boldsymbol A$ and $\boldsymbol \beta$?

Answer 1

The history of Quasi-Newton methods is replete with examples of optimizing the step-wise matrix update with respect to various merit functions.

In the 1970s, Broyden (and the other initials of BFGS) used the Frobenius norm: $$\big\|C_{k+1}-C_k\big\|_F$$

In the 1980s, Michael Todd used the condition number: $$\kappa(C)=\frac{\lambda_{max}}{\lambda_{min}}$$

In the 1990s, Byrd & Nocedal introduced the merit function: $$\psi(C) = {\rm Tr}(C) - \log(\det(C))$$

and Wolkowicz introduced his "sizing" metric: $$\omega(C) = \frac{\frac 1n{\rm Tr}(C)}{(\det C)^{1/n}} = \frac{\lambda_{arithmetic\_mean}}{\lambda_{geometric\_mean}}$$

This new KL-metric is obviously inspired by recent developments in Machine Learning, but it's not anything radically new, or even all that interesting.

If you wish to learn how to optimize with respect to some novel en vogue merit function, go back and read those classic papers from the 70s, 80s, and 90s.

Answer 2

What we want to minimize is the Stein loss, aka the KL divergence between two Gaussians of same mean, aka the Bregman divergence induced by the negative logdet. I'll write in slightly different notation because those are the ones I have already typed in my class.

$$\min_H - \mathrm{log\,det} (H) + \langle H_k^{-1}, H \rangle \quad \text{subject to} \quad H = H^\top, Hs = y \enspace.$$ Introduce the Lagrangian: $$\mathcal{L}(H, \lambda, \Theta) = - \log \det H + \langle H_k^{-1}, H \rangle + \langle \lambda, Hs - y \rangle + \langle \Theta, H - H^\top \rangle .$$

Notice that $\langle \lambda, Hs \rangle$ is a scalar, equal to its trace $\mathrm{tr}\, \lambda^\top Hs = \mathrm{tr} \, s \lambda^\top H = \langle \lambda s^\top, H \rangle$ where this time the scalar product is over matrices.

Similarly, $\langle \Theta, H - H^\top \rangle = \langle \Theta, H \rangle - \langle \Theta, H^\top \rangle = \langle \Theta, H \rangle - \langle \Theta^\top, H \rangle = \langle \Theta - \Theta^\top, H \rangle$.

This gives us the two gradients we need to write the first order optimality condition over $H$, $\nabla_H \mathcal{L} = 0$. Using that the gradient of logdet is the inverse, we obtain: $$H^{-1} = H_k^{-1} + \lambda s^\top + \Theta^\top - \Theta$$

First we get rid of $\Theta$: $H$ must be symmetric hence $H^{-1}$ too; $H_k$ is symmetric (it's the approximate Hessian at the previous iteration) and this gives: $$\Theta^\top - \Theta = \frac12 ( s \lambda^\top - \lambda s^\top)$$ so $$H^{-1} = H_k^{-1} + \frac{1}{2} (\lambda s^\top + s \lambda^\top)$$ To find $\lambda$, we use the other condition, $Hs = y$ aka $s = H^{-1} y$. This gives $$\frac12 s^\top y \lambda = (s - H_k^{-1} y) - \frac12 \lambda^\top y s$$ This may seem impossible to solve as $\lambda$ appears twice, but it nevertheless tells us something: there exists a scalar $a$ such that $$ \lambda = \frac{2}{s^\top y} (s - H_k^{-1} y) + a s$$ Plug back this expression to obtain $H^{-1}$ as a function of $a$, then use $H^{-1} y = s$ to solve in $a$ (it's tedious, but a good exercise).

Then plugging back $a$ into $\lambda$ we obtain $$\lambda = \frac{2}{s^\top y} (s - H_k^{-1} y) + \left(\frac{1}{(s^\top y)^2} y^\top H_k^{-1} y - \frac{1}{s^\top y} \right) s$$ and this gives the desired solution for $H^{-1}$.