Linear Algebra – Applying Sherman-Morrison-Woodbury for Rank 2 Update

linear algebramatricesoptimization

I was reading Nocedal and Wright, and it is stated that one may use the Sherman-Morrison-Woodbury formula
$$(A+ YGZ^*)^{-1} = A^{-1} – A^{-1}Y(G^{-1}+Z^*A^{-1}Y)^{-1}Z^*A^{-1}$$

On the hessian $H_{k+1}$ inverse approximation from the BFGS:
$$ H_{k+1}= (I – \gamma s_k y_k ^T)H_k(I-\gamma_k y_k s_k^T) + \gamma_k s_k s_k^T$$

To obtain an update formula for the Hessian $B_{k+1}$ for the BFGS:

$$B_{k+1} = B_k – \frac{B_ks_ks_k^T B_k}{s_k^TB_k s_k } + \gamma_k y_ky_k^T $$

Where the setting here is $\gamma_k = \frac{1}{y_k^Ts_k}$ and $H_k = B_k^{-1}$ and $H_k$ and $B_k$ are positive definite symmetric matrices.

My issue is that I do not see how to apply the Sherman-Morrison-Woodbury formula.
What I see is that I want
$H_{k+1} = A+ YGZ^*$ and that would give me $B_{k+1}$

Then the question is what is $A$, $Y$, $G$, and $Z$. My thoughts are that $A \not = \gamma_k y_ky_k^T$ as that is singular. Thus that leaves $A = (I – \gamma s_k y_k ^T)H_k(I-\gamma_k y_k s_k^T)$, but that is a mess to invert and when I tried it I did not get anything resembling the desired expression and also then what do I define $G,Y,Z$ to be?

Best Answer

Hint: write the RHS (I will skip the index $k$) as $$ H-\gamma sy^TH-\gamma Hys^T+\gamma^2sy^THys^T+\gamma ss^T= H+\gamma\begin{bmatrix}s & Hy\end{bmatrix}\begin{bmatrix}1+\gamma y^THy & -1\\-1 & 0\end{bmatrix}\begin{bmatrix}s^T\\y^TH\end{bmatrix}. $$ Can you finish now?

Related Solutions

[Math] Sherman–Morrison–Woodbury formula question

You don't really fix this, you should just avoid taking inverses of nearly-singular matrices like the one you have defined. Even though $B$ in this case has an exact inverse, in a numerically-approximate calculation that R does it really matters whether it is approximately singular or not. R's function rcond will tell you the reciprocal condition number of $B$, and of $C^{-1}+V A^{-1} U$; if the answer is a very small number (e.g., $10^{-9}$ or less), then the numerically-approximate computation is poorly conditioned and any answer will be rather inaccurate.

To be honest, I don't understand why you are doing this at all. Finding an eigen-decomposition of the (almost singular) matrix $B$ is hard, and you gain nothing by doing this only to insert this into the Woodbury formula. If all you know about $B$ is that it is symmetric, just calculate $A+B$ and invert that. In general, it makes very little sense to invert an arbitrary matrix by computing its eigen-decomposition. Even if a matrix is symmetric (which guarantees that it can be diagonalized), why not just invert it directly? (R will do it for you easily.) Not to mention that while $B$ is almost singular and the inverse of its diagonal form is needed to apply the formula, the matrix $A+B$ is not quite so singular, and can be inverted easily.

You are not saving any work by using Woodbury identity on a generic matrix like $A+B$. Usually, it makes a lot more sense to use it when the matrix $BCD$ (in $A+BCD$) has some very special form that makes calculating the inverses in the Woodbury identity easy: often $C$ is 1, and $B^t$ and $D$ are vectors, or $2\times n$ matrices -- then the only inverse in the problem that is hard is $A^{-1}$; for example, $DA^{-1}B$ means just solving a system of linear equations when $B$ is a vector. With the matrix $B$ that you wrote down this does not hold: $B$ is just a generic matrix, so Woodbury identity offers no help.

Finally, in general it is not advisable to calculate the inverses of matrices directly. In other words, the usual practice is not to calculate $A^{-1}$ or any other inverse explicitly, but instead whenever you need to find $A^{-1}X$, solve the set of systems of linear equations $AY=X$ instead of computing $A^{-1}$ and multiplying by $X$. So whenever in your code you have solve(A) %*% X, it should really be replaced with solve(A,X).

Are you familiar with the standard techniques for solving systems of linear equations, like Gaussian elimination and LU factorization and the like, and how they would be used to calculate matrix inverses? Your question can probably be better answered by a chapter on linear equations in an introductory numerical analysis textbook.

Optimization – Understanding DFP Updating Formula in Quasi-Newton Methods

This answer goes basically along the lines of my answer about BFGS update.

Introduce the short notations $$ \min\|\underbrace{W^{1/2}B_kW^{1/2}}_{\hat B_k}-\underbrace{W^{1/2}BW^{1/2}}_{\hat B}\|_F $$ \begin{align} Wy_k=s_k\quad&\Leftrightarrow\quad \underbrace{\color{red}{W^{-1/2}}Wy_k}_{\hat y_k}=\underbrace{\color{red}{W^{-1/2}}s_k}_{\hat s_k}\quad\Leftrightarrow\quad \hat y_k=\hat s_k,\tag{1}\\ Bs_k=y_k\quad&\Leftrightarrow\quad \underbrace{\color{blue}{W^{1/2}}B\color{red}{W^{1/2}}}_{\hat B}\underbrace{\color{red}{W^{-1/2}}s_k}_{\hat s_k}=\underbrace{\color{blue}{W^{1/2}}y_k}_{\hat y_k}\quad\Leftrightarrow\quad \hat B\hat s_k=\hat y_k.\tag{2} \end{align} Then the problem becomes $$ \min\|\hat B_k-\hat B\|_F\quad\text{subject to }\hat B\hat s_k=\hat s_k. $$
The optimization variable $\hat B$ has the given eigenvector with the given eigenvalue, hence, it is convenient to introduce the new orthonormal basis $$ U=[u\ |\ u_\bot] $$ where $u$ is the normalized eigenvector $\hat s_k$, i.e. $$ u=\frac{\hat s_k}{\|\hat s_k\|}\tag{3}, $$ and $u_\bot$ is any orthonormal complement to $u$. Then $u^T\hat Bu=u^Tu=1$ and $u_\bot^T\hat Bu=u_\bot^Tu=0$, and the operator matrices in the new basis take the form \begin{align} U^T\hat B_kU-U^T\hat BU&=\begin{bmatrix}u^T\\ u_\bot^T\end{bmatrix}\hat B_k\begin{bmatrix}u & u_\bot\end{bmatrix}-\begin{bmatrix}u^T\\ u_\bot^T\end{bmatrix}\hat B\begin{bmatrix}u & u_\bot\end{bmatrix}=\\ &=\begin{bmatrix}\color{blue}{u^T\hat B_ku} & \color{blue}{u^T\hat B_ku_\bot}\\\color{blue}{u_\bot^T\hat B_ku} & \color{red}{u_\bot^T\hat B_ku_\bot}\end{bmatrix}-\begin{bmatrix}\color{blue}{1} & \color{blue}{0}\\\color{blue}{0} & \color{red}{u_\bot^T\hat Bu_\bot}\end{bmatrix}. \end{align} Since the Frobenius norm is unitary invariant (as it depends on the singular values only) we have \begin{align} \|\hat B_k-\hat B\|_F^2&=\|U^T(\hat B_k-\hat B)U\|_F^2= \left\|\begin{bmatrix}\color{blue}{u^T\hat B_ku-1} & \color{blue}{u^T\hat B_ku_\bot}\\\color{blue}{u_\bot^T\hat B_ku} & \color{red}{u_\bot^T\hat B_ku_\bot-u_\bot^T\hat Bu_\bot}\end{bmatrix}\right\|_F^2=\\ &=\color{blue}{(u^T\hat B_ku-1)^2+\|u^T\hat B_ku_\bot\|_F^2+\|u_\bot^T\hat B_ku\|_F^2}+\color{red}{\|u_\bot^T\hat B_ku_\bot-u_\bot^T\hat Bu_\bot\|_F^2} \end{align} The blue part cannot be affected by optimization, and to minimize the Frobenius norm, it is clear that we should make the red part zero, that is, the optimal solution satisfies $$ \color{red}{u_\bot^T\hat Bu_\bot}=\color{red}{u_\bot^T\hat B_ku_\bot}. $$
It gives the optimal solution to be \begin{align} \hat B&=U\begin{bmatrix}\color{blue}1 & \color{blue}0\\\color{blue}0 & \color{red}{u_\bot^T\hat B_ku_\bot}\end{bmatrix}U^T=\begin{bmatrix}u & u_\bot\end{bmatrix}\begin{bmatrix}1 & 0\\0 & u_\bot^T\hat B_ku_\bot\end{bmatrix}\begin{bmatrix}u^T \\ u_\bot^T\end{bmatrix}=uu^T+u_\bot u_\bot^T\hat B_ku_\bot u_\bot^T=\\ &=uu^T+(I-uu^T)\hat B_k(I-uu^T) \end{align} where we used the following representation for the projection operator to the complement of $u$ $$ I=UU^T=\begin{bmatrix}u & u_\bot\end{bmatrix}\begin{bmatrix}u^T \\ u_\bot^T\end{bmatrix}=uu^T+u_\bot u_\bot^T\quad\Leftrightarrow\quad u_\bot u_\bot^T=I-uu^T. $$
Changing variables back to the original ones is straightforward via (1), (2), (3) $$ B=W^{-1/2}\hat BW^{-1/2}=W^{-1/2}uu^TW^{-1/2}+(I-W^{-1/2}uu^TW^{1/2})B_k(I-W^{1/2}uu^TW^{-1/2}) $$ where \begin{align} \|\hat s_k\|^2&=\hat s_k^T\hat s_k=\hat y_k^T\hat s_k=y_k^Ts_k,\\ uu^T&=\frac{\hat s_k^T\hat s_k}{\|\hat s_k\|^2}=\frac{\hat y_k\hat y_k^T}{y_k^Ts_k}=W^{1/2}\frac{y_ky_k^T}{y_k^Ts_k}W^{1/2},\\ W^{-1/2}uu^TW^{-1/2}&=\frac{y_ky_k^T}{y_k^Ts_k},\\ W^{1/2}uu^TW^{-1/2}&=\frac{Wy_ky_k^T}{y_k^Ts_k}=\frac{s_ky_k^T}{y_k^Ts_k}. \end{align}

Best Answer

Related Solutions

[Math] Sherman–Morrison–Woodbury formula question

Optimization – Understanding DFP Updating Formula in Quasi-Newton Methods

Related Question