Woodbury Matrix Inverse Identity

July 31, 2013
By

(This article was originally published at Lindons Log » Statistics, and syndicated at StatsBlogs.)

Application in Conditional Distribution of Multivariate Normal

The Sherman-Woodbury-Morrison matrix inverse identity can be regarded as a transform between Schur complements. That is, given V_{22.1}^{-1} one can obtain V_{11.2}^{-1} by using the Woodbury matrix identity and vice versa. Recall the Woodbury Identity:

V_{11.2}^{-1}=V_{11}^{-1}+V_{11}^{-1}V_{12}V_{22.1}^{-1}V_{21}V_{11}^{-1}

and

V_{22.1}^{-1}=V_{22}^{-1}+V_{22}^{-1}V_{21}V_{11.2}^{-1}V_{12}V_{22}^{-1}

I recently stumbled across a neat application of this whilst deriving full conditionals for a multivariate normal. Recall that if the data are partitioned into two blocks, Y_{1},Y_{2}, then the variance of the conditional distribution Y_{1}|Y_{2},- is the Schur complement of the block V_{22} of total variance matrix V, that is, the variance of the conditional distribution is V_{11.2}=V_{11}-V_{12}V_{22}^{-1}V_{21} which is the variance of Y_{1} subtracted by something corresponding to the reduction in uncertainty about Y_{1} gained from the knowledge about Y_{2}. If, however, V_{22} has the form of a Schur complement itself, then it may be possible to exploit the Woodbury identity above to considerably simplify the variance term. I came across this when I derived two very different-looking expressions for the conditional distribution and found them equivalent by the Woodbury identity. Consider the model

\begin{bmatrix} Y_{o}\\ Y_{a} \end{bmatrix} = \begin{bmatrix} X_{o\gamma}\\ X_{a\gamma} \end{bmatrix}\beta_{\gamma} + \varepsilon

where

\varepsilon \sim N\left( \begin{bmatrix}0\\ 0\end{bmatrix}, \sigma^{2} \begin{bmatrix}I_{no} & 0 \\ 0 & I_{na}\end{bmatrix} \right)

\beta_{\gamma}|\gamma,\sigma^{2} \sim N(0, \sigma^{2}\Lambda_{\gamma}^{-1})

.
I was seeking the distribution Y_{o}|Y_{a},\sigma^{2},\gamma and arrived there through two different paths. The distributions derived looked very different, but they turned out to be equivalent upon considering the Woodbury identity.

Method 1

This simply manipulates properties of the multivariate normal. Marginalizing over \beta one gets

Cov \begin{bmatrix} Y_{o} \\ Y_{a} \end{bmatrix} = \begin{bmatrix} X_{o\gamma} \\ X_{a\gamma} \end{bmatrix} Cov (\beta_{\gamma}) \begin{bmatrix} X_{o\gamma}^{T} & X_{a\gamma}^{T} \end{bmatrix} + Cov(\varepsilon)

Cov \begin{bmatrix} Y_{o} \\ Y_{a} \end{bmatrix} = \sigma^{2}\begin{bmatrix} X_{o\gamma}\Lambda_{\gamma}^{-1}X_{o\gamma}^{T} & X_{o\gamma}\Lambda_{\gamma}^{-1}X_{a\gamma}^{T} \\ X_{a\gamma}\Lambda_{\gamma}^{-1}X_{o\gamma}^{T} & X_{a\gamma}\Lambda_{\gamma}^{-1}X_{a\gamma}^{T} \end{bmatrix} + \sigma^{2} \begin{bmatrix} I_{no} & 0 \\ 0 & I_{na} \end{bmatrix}

.
Such that the distribution

\begin{bmatrix} Y_{o}\\ Y_{a} \end{bmatrix}|\gamma,\sigma^{2} \sim N \left( \begin{bmatrix} 0\\ 0 \end{bmatrix}, \sigma^{2} \left( \begin{bmatrix} I_{no}+ X_{o\gamma}\Lambda_{\gamma}^{-1}X_{o\gamma}^{T} & X_{o\gamma}\Lambda_{\gamma}^{-1}X_{a\gamma}^{T} \\ X_{a\gamma}\Lambda_{\gamma}^{-1}X_{o\gamma}^{T} & I_{na}+X_{a\gamma}\Lambda_{\gamma}^{-1}X_{a\gamma}^{T} \end{bmatrix} \right) \right)

It follows that the conditional distribution is
Y_{o}|Y_{a},\gamma,\sigma^{2} \sim N \left( X_{o\gamma}\Lambda_{\gamma}^{-1}X_{a\gamma}^{T} \left[ X_{a\gamma}\Lambda_{\gamma}^{-1}X_{a\gamma}^{T} + I_{na}\right]^{-1}Y_{a}, I_{no} + X_{o\gamma}\Lambda_{\gamma}^{-1}X_{o\gamma}^{T} - X_{o\gamma}\Lambda_{\gamma}^{-1}X_{a\gamma}^{T} \left[ I_{na} + X_{a\gamma}\Lambda_{\gamma}^{-1}X_{a\gamma}^{T} \right]^{-1} X_{a\gamma}\Lambda_{\gamma}^{-1}X_{o\gamma}^{T}\right).
This looks a bit nasty, but notice that V_{22}^{-1} looks like it too could be a Schur complement of some matrix.

Method 2

An alternative route to this distribution is

f(Y_{o}|Y_{a},\sigma^{2},\gamma)=\int f(Y_{o}|\sigma^{2},\gamma,\beta_{\gamma})\pi(\beta_{\gamma}|Y_{a},\sigma^{2},\gamma)d\beta_{\gamma}

where

\beta_{\gamma}|Y_{a},\gamma,\sigma^{2}\sim N \left( (X_{a\gamma}^{T}X_{a\gamma}+\Lambda_{\gamma})^{-1}X_{a\gamma}^{T}Y_{a}, \sigma^{2}(X_{a\gamma}^{T}X_{a\gamma}+\Lambda_{\gamma})^{-1} \right).

It follows that

Y_{o}|Y_{a},\gamma,\sigma^{2} \sim N\left( X_{o\gamma}(X_{a\gamma}^{T}X_{a\gamma}+\Lambda_{\gamma})^{-1}X_{a\gamma}^{T}Y_{a}, \sigma^{2} (I_{no} + X_{o\gamma} (X_{a\gamma}^{T}X_{a\gamma}+\Lambda_{\gamma})^{-1}X_{o\gamma}^{T}) \right)

which looks different from the distribution obtained through method 1. The expression for the variance is a lot neater. They are in fact identical by the Woodbury identity.

Comparison

Mean

X_{o\gamma}\Lambda_{\gamma}^{-1}X_{a\gamma}^{T} \left[ X_{a\gamma}\Lambda_{\gamma}^{-1}X_{a\gamma} + I_{na}\right]^{-1}Y_{a}\\ =X_{o\gamma} \left[ X_{a\gamma} + X_{a\gamma}^{-T}\Lambda_{\gamma} \right]^{-1}Y_{a}\\ =X_{o\gamma} \left[ X_{a\gamma}^{-T}X_{a\gamma}^{T}X_{a\gamma} + X_{a\gamma}^{-T}\Lambda_{\gamma}\right]^{-1}Y_{a}\\ =X_{o\gamma} \left[ X_{a\gamma}^{T}X_{a\gamma} + \Lambda_{\gamma}\right]^{-1} X_{a\gamma}^{T}Y_{a}\

So mean1=mean2.

Variance

By the Woodbury Identity it follows that

\Lambda_{\gamma}^{-1} - \Lambda_{\gamma}^{-1}X_{a\gamma}^{T} \left[ I_{na} + X_{a\gamma}\Lambda_{\gamma}^{-1}X_{a\gamma}^{T} \right]^{-1} X_{a\gamma}\Lambda_{\gamma}^{-1} = (X_{a\gamma}^{T}I_{na}X_{a\gamma}+\Lambda_{\gamma})^{-1}.

Therefore

X_{o\gamma}\Lambda_{\gamma}^{-1}X_{o\gamma}^{T}-X_{o\gamma}\Lambda_{\gamma}^{-1}X_{a\gamma}^{T} \left[ I_{na}+X_{a\gamma}\Lambda_{\gamma}^{-1}X_{a\gamma}^{T} \right]^{-1} X_{a\gamma}\Lambda_{\gamma}^{-1}X_{o\gamma}^{T}={X_{o\gamma} (X_{a\gamma}^{T}X_{a\gamma}+\Lambda_{\gamma})^{-1}X_{o\gamma}^{T}}\\

and so variance1=variance2. The trick is recognizing the form of the formulas at the top of the page, then one can write the variance as a much neater expression.

The post Woodbury Matrix Inverse Identity appeared first on Lindons Log.



Please comment on the article here: Lindons Log » Statistics

Tags: , ,

Subscribe

Email:

  Subscribe