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.
 A: 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.
A: 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
\begin{equation}
\hat{x}=\Sigma(\Sigma+\sigma^2)^{-1}y
\end{equation}
or even better (after seeing Alexey's answer)
\begin{equation}
\hat{x}=y - \left (I+\frac{1}{\sigma^2}\Sigma \right )^{-1}y
\end{equation}
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.
