I'm trying to find the MAP estimate for a model by gradient descent. My prior is multivariate Gaussian with a known covariance matrix.

On a conceptual level, I think I know how to do this, but I was hoping for some help with the details. In particular, if there is an easier way to approach the problem, then that would be especially useful.

Here's what I think I need to do:

• For each dimension, find the conditional distribution, given my current position in the other dimensions.
• This gives me a local univariate Gaussian in each dimension, with the correct mean and standard deviation.
• I think that the gradient should just be a vector of derivatives for each of these univariate distributions.

My question has two parts:

1. Is this the best approach to take, or is there an easier way?
2. Assuming I need to go this route, what's the best way to go about finding these conditional distributions?
-
Is there any reason why you want to do this with gradient descent? Finding the MAP of a MVN with some prior sounds like a fairly well studied problem. Since the MVN is self-conjugate, a fully Bayesian approach should even be possible. –  bayerj Mar 15 '14 at 19:05
@bayerj Good question. The prior is MVN, but the likelihood isn't. I think that limits my options. –  David J. Harris Mar 15 '14 at 19:06
Ah ok, I did not get that. –  bayerj Mar 15 '14 at 19:59

Let's see if I understand you correctly. You have a model $p(y|x, \theta)$ conditioned on some observation $x$ and a set of parameters $\theta$ and a prior $p(\theta)$ leading to a joint likelihood of $\mathcal{L} = p(y|x, \theta)p(\theta)$. The parameters are distributed according to a known multivariate normal, i.e. $\theta \sim \mathcal{N}(\mu, \Sigma)$. You want to find the MAP solution to this problem, i.e. $$\text{argmax}_{\theta} \mathcal{L}.$$ A special case of this problem is well studied in the neural networks community, known as weight decay. In that case, $\mu=\mathbf{0}$ and $\Sigma = \mathbf{I}\sigma^2$.

As you already noted, the trick is that $\text{argmax}_{\theta} \mathcal{L} = \text{argmax}_{\theta} \log \mathcal{L}$. When you take the log of the Gaussian density, many ugly terms (the exponential) vanish and you will end up with sth like $\log p(\theta) = {1 \over 2}(\theta - \mu)^T\Sigma^{-1}(\theta - \mu) + \text{const}$. If you differentiate that, Sam Roweis' matrix identities will come in handy and let you arrive at

$${1 \over 2}{\partial (\theta - \mu)^T\Sigma^{-1}(\theta - \mu) \over \partial \theta} = \Sigma^{-1}(\theta - \mu).$$

(Please verify, this was done quickly and in my head.) Together with the derivatives of your model, you can use off-the-shelf optimizers to arrive at a MAP solution.

Update: Incorporated comment by David J. Harris. Formulas should be correct now.

-
(+1) This looks like exactly what I need. I'm going to do a bit of verification this afternoon and then I'll hit the "accept" checkmark if everything works out. Thank you! –  David J. Harris Mar 15 '14 at 20:49
I forgot to add: if your model is simple (i.e. linear in the parameters), equating the derivative of the log-likelihood with zero and solving for $\theta$ might even work out. –  bayerj Mar 15 '14 at 21:06
I've played around with this a little bit numerically, and I think it's off by a factor of 2 but otherwise correct. Perhaps it cancels with the $1/2$ from the formula for the multivariate normal density? Thanks again! –  David J. Harris Mar 15 '14 at 21:41
That should be it, yes. I forgot about it! –  bayerj Mar 16 '14 at 13:11
It should be noted that a custom multivariate Gaussian prior (with non-diagonal covariance matrix) defines Tikhonov regularization –  Barmaley.exe Dec 26 '14 at 8:09

If the likelihood is not Gaussian it's not possible to say if there are analytic results. Also, the second bullet is then incorrect in general. Since Gaussian prior and general likelihood does not make for conditional gaussian distributions on the vector components.

One way to get the MAP would be to do a full Bayesian analysis, e.g. using MCMC and use the samples from the posterior to estimate it. [In which case you'd have better information available than only using the MAP.] Out of interest- why not go down this route anyhow?

Another approach could be to do (I've not seen this done generally so someone please correct me if it's nuts):

$p(\theta|x) = \frac{p(x|\theta)p(\theta)}{p(x)}$

$l(\theta|x) = l(x|\theta) + l(\theta) - l(x)$

$\frac{dl(\theta|x)}{d\theta} = \frac{dl(x|\theta)}{d\theta} + \frac{dl(\theta)}{d\theta} = 0$

Then solve for $\theta$ (probably numerically).

-
thanks for your input. I may not have been clear: right now, I'm just interested in finding the gradient for the prior. The gradient of the log-posterior is just the gradient of the log-likelihood plus the gradient of the log-prior, so finding these two gradients separately should be okay. –  David J. Harris Mar 15 '14 at 21:45
Ehm, your proposal to solve for $\theta$ numerically is... the basic approach. How else is one supposed to optimize the target function without finding a root of a gradient? –  Barmaley.exe Dec 26 '14 at 8:06