Evaluate high-dimensional Gaussian with variance matrix $\sigma^{2}I_{n_{t}\times n_{t}}+\boldsymbol{\Sigma}_{t}\boldsymbol{\Sigma}_{t}^{'}$

Question

I need to compute the log-likelihood function in a high-dimensional Gaussian time-series. I have the following model:

$\mathbf{y}_{t}\left|\mathcal{F}_{t-1}\sim\mathcal{N}\left(\mathbf{\boldsymbol{\mu}}_{t},\sigma^{2}I_{n_{t}\times n_{t}}+\boldsymbol{\Sigma}_{t}\boldsymbol{\Sigma}_{t}^{'}\right)\right.$

here $\sigma^{2}$ is a constant scaler. $\boldsymbol{\Sigma}_{t}$ is a time-varying $n_{t}\times1$ vector where $n_{t}$ is big. $I_{n_{t}\times n_{t}}$ is the identity matrix of dimension $n_{t}$.

The log-likelihood will then be given by:

$\log p\left(\mathbf{y}_{1},...,\mathbf{y}_{T}\right)=\sum_{t=1}^{T}-\frac{1}{2}\left(n_{t}\log\left(2\pi\right)+\log\left(\left|\sigma^{2}I_{n_{t}\times n_{t}}+\boldsymbol{\Sigma}_{t}\boldsymbol{\Sigma}_{t}^{'}\right|\right)+\left(\mathbf{y}_{t}-\mathbf{\boldsymbol{\mu}}_{t}\right)^{'}\left(\sigma^{2}I_{n_{t}\times n_{t}}+\boldsymbol{\Sigma}_{t}\boldsymbol{\Sigma}_{t}^{'}\right)^{-1}\left(\mathbf{y}_{t}-\mathbf{\boldsymbol{\mu}}_{t}\right)\right)$

Given the special structure of the variance matrix

$\sigma^{2}I_{n_{t}\times n_{t}}+\boldsymbol{\Sigma}_{t}\boldsymbol{\Sigma}_{t}^{'}$

is there anyway to exploit this structure to evaluate the inverse, determinant and in turn the full log-likelihood computationally fast?

Answer 1

To calculate likelihood you need to calculate the determinant $|\sigma^2 I + \Sigma \Sigma^T|$ and the quadratic function $(\mathbf{y} - \boldsymbol{\mu}_t)^{\mathrm T} (\sigma^2 I + \Sigma \Sigma^{\mathrm T})^{-1} (\mathbf{y} - \boldsymbol{\mu}_t)$.

It holds that $|A| = \prod_{i}^{n_t} a_i$, where $a_i$ are eigenvalues of $A$. Matrix $\sigma^2 I + \Sigma \Sigma^{\mathrm T}$ has $n_t - 1$ eigenvalues $\sigma^2$ and one eigenvalue $\sigma^2 + \Sigma^{\mathrm T} \Sigma$ (see for example this question). So, $$|\sigma^2 I + \Sigma \Sigma^T| = \sigma^{2(n_t - 1)} (\sigma^2 + \Sigma^{\mathrm T} \Sigma).$$ Complexity of determinant calculation is $O(n_t)$.

To get inverse of the matrix $\sigma^2 I + \Sigma \Sigma^{\mathrm T}$ you should use Sherman–Morrison formula: $$(\sigma^2 I + \Sigma \Sigma^{\mathrm T})^{-1} = \frac{1}{\sigma^2} I - \frac{1}{\sigma^4} \frac{\Sigma \Sigma^{\mathrm T}}{1 + \frac{1}{\sigma^2} \Sigma^{\mathrm T} \Sigma}.$$ To further reduce computation cost you can calculate $\xi = \Sigma^T (\mathbf{y} - \boldsymbol{\mu}_t)$, as $$(\mathbf{y} - \boldsymbol{\mu}_t)^{\mathrm T} (\sigma^2 I + \Sigma \Sigma^{\mathrm T})^{-1} (\mathbf{y} - \boldsymbol{\mu}_t) = \frac{1}{\sigma^2} (\mathbf{y} - \boldsymbol{\mu}_t)^{\mathrm T} (\mathbf{y} - \boldsymbol{\mu}_t) - \frac{1}{\sigma^4} \frac{\xi^2}{1 + \frac{1}{\sigma^2} \Sigma^{\mathrm T} \Sigma}.$$ As we need only to multiply vectors of size $n_t$ computational complexity is again $O(n_t)$.

So, we get the results for one $\boldsymbol{\mu}_t$ using $O(n_t)$ operations. For different $t$ we then get total computational cost $O(\sum_{t = 1}^T n_t)$.