# Computationally efficient Gaussian MAP estimation algorithm in MATLAB

I have a MAP estimation model for a Gaussian prior and i.i.d Gaussian noise:

$$y=x+n$$

where $x\sim\mathcal{N}(0,\Sigma)$ and $n\sim \mathcal{N}(0,\sigma^2I)$. The MAP estimate is given by

$$\hat{x} = \arg\max P(X|Y) = (I+\sigma^2\Sigma^{-1})^{-1} y$$

Everything works fine, but the direct application of this, including two matrix inversions, is not computationally efficient. Is there any general purpose algorithm to make this work faster?

Thanks.

• You could look to approximate the outer inversion using a Taylor expansion and truncate at a suitable point. Commented Mar 2, 2015 at 11:18
• Having said that, I don't think you should be doing a second inversion to solve a problem of the form $x = A^{-1}y$ as it is an inefficient way to do it. Commented Mar 2, 2015 at 11:28
• @rwolst , thank you, I am indeed using MATLAB's LS solution (the backslash operator) instead of the last inversion.
– yoki
Commented Mar 2, 2015 at 12:14

While this method doesn't do anything clever in terms of the structure or algorithm, it is quicker in Matlab to do the following:

Notice that $$\hat{x}=\Sigma(\Sigma+\sigma^2)^{-1}y$$

or even better (after seeing Alexey's answer) $$\hat{x}=y - \left (I+\frac{1}{\sigma^2}\Sigma \right )^{-1}y$$

We can compare these in Matlab to the initial naive implementation using the code

n = 2000;

A = rand(n);
y = rand(n,1);

% Naive
tic
inv(eye(n) + inv(A))*y;
toc

% Mine
tic
B = A+eye(n);
A*(B\y);
toc

% Alexey
tic
B = A+eye(n);
y - B\y;
toc


I get answers of

Elapsed time is 1.126651 seconds.
Elapsed time is 0.246517 seconds.
Elapsed time is 0.202166 seconds.


In order to get any more really significant speedup, I'm guessing you would have to start taking advantage of the structure of your matrices if possible.

I think the first thing you should try is to use Woodbury matrix identity. Using it you get from your initial formula: $$(I + \sigma^2 \Sigma^{-1})^{-1} = I - I \left(\frac{1}{\sigma^2} \Sigma + I \right)^{-1} I = I - \left(\frac{1}{\sigma^2} \Sigma + I \right)^{-1}.$$

Using this formula you require only one matrix inversion instead of two. Another benefit is that you have a nonnegative matrix, and as you add an identity matrix to it, you get a matrix to inverse with mininal eigenvalue at least 1.

Another thing to try is to explore structure of your $\Sigma$ matrix.