# BFGS Formula from Kullback-Leibler Divergence

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$$?

• You're basically asking how to solve exercise 10 on page 352 of the book. It should be relatively straight forward to use the Lagrange Multipliers/KKT conditions to do this. If you expand your description of your attempt to solve this with Lagrange multipliers, then perhaps we could find your error. Jun 30, 2021 at 20:45
• Maybe you can reformulate your constrained optimization problem to the same optimization problem from this question here: math.stackexchange.com/questions/2271887/… and then use the same idea from there? This other question about DFP goes in the same direction: math.stackexchange.com/questions/2091867/… Jul 8, 2021 at 13:24

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.

• I think this answer is unnecessary patronizing. 1) The question asks about how to solve an optimization problem, and the answer points to other optimization problems. 2) The metric to optimize is known as Stein's loss and dates back to at least 1961; it is very popular in statistics and, in my opinion, is not just some "en vogue" function. Feb 22, 2023 at 16:10

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