Comparing Laplace Approximation and Variational Inference Does anyone know of any references that look at the relationship between the Laplace approximation and variational inference (with normal approximating distributions)? Namely I'm looking for something like conditions on the distribution being approximated for the two approximations to coincide.
Edit: To give some clarification, suppose you want to approximate some distribution with density $f(\theta)$ which you only know up to proportionality. When using the Laplace approximation, you approximate it with the density of a Normal distribution with mean $\hat{\mu}_1$ and covariance $\hat{\Sigma}_1$ where $\hat{\mu}_1=\arg \max_{\theta}f(\theta)$ and $\hat{\Sigma}_1=[-\nabla\nabla \log f(\theta)\mid_{\theta=\hat{\mu}}]^{-1}$. When using variational inference with a normal approximating distribution, you approximate it with the density of a Normal distribution with mean $\hat{\mu}_2$ and covariance $\hat{\Sigma}_2$, where $(\hat{\mu}_2,\hat{\Sigma}_2)=\arg \max_{(\mu,\Sigma)}KL(\phi_{(\mu,\Sigma)}||f)$, $KL$ is the KL-Divergence, and $\phi_{(\mu,\Sigma)}$ denotes a Normal density with mean and covariance $(\mu,\Sigma)$. Under what conditions do we have $(\hat{\mu}_1,\hat{\Sigma}_1)=(\hat{\mu}_2,\hat{\Sigma}_2)$?
 A: I am not aware of any general results, but in this paper the authors have some thoughts for Gaussian variational approximations (GVAs) for generalized linear mixed model (GLMMs). Let $\vec y$ be the observed outcomes, $X$ be a fixed effect design matrix, $Z$ be a random effect design, denote an unknown random effect $\vec U$, and consider a GLMM with densities:
$$
\begin{align*}
f_{\vec Y\mid\vec U} (\vec y;\vec u) &= 
  \exp\left(\vec y^\top(X\vec\beta + Z\vec u) 
  - \vec 1^\top b(X\vec\beta + Z\vec u) 
  + \vec 1^\top c(\vec y)\right) \\
f_{\vec U}(\vec u) &= \phi^{(K)}(\vec u;\vec 0, \Sigma) \\
f(\vec y,\vec u) &= f_{\vec Y\mid\vec U} (\vec y;\vec u)f_{\vec U}(\vec u)
\end{align*}
$$
where I use the same notation as in the paper and $\phi^{(K)}$ is a $K$-dimensional multivariate normal distribution density function.
Using a Laplace Approximation
Let
$$
g(\vec u) = \log f(\vec y,\vec u).
$$
Then we use the approximation
$$
\log\int \exp(g(\vec u)) d\vec u \approx 
  \frac K2\log{2\pi  - \frac 12\log\lvert-g''(\widehat u)\rvert}
  + g(\widehat u)
$$
where
$$
\widehat u = \text{argmax}_{\vec u} g(\vec u).
$$
Using a Gaussian Variational Approximation
The lower bound in the GVA with a mean
$\vec\mu$ and covariance matrix $\Lambda$ is:
$$
\begin{align*}
\int \exp(g(\vec u)) d\vec u &\approx
   \vec y^\top(X\vec\beta + Z\vec\mu) 
  - \vec 1^\top B(X\vec\beta + Z\vec\mu, \text{diag}(Z\Lambda Z^\top)) \\ 
&\hspace{25pt}+ \vec 1^\top c(\vec y) + \frac 12 \Big( 
  \log\lvert\Sigma^{-1}\rvert + \log\lvert\Lambda\rvert
  -\vec\mu^\top\Sigma^{-1}\vec\mu \\
&\hspace{25pt} - \text{trace}(\Sigma^{-1}\Lambda)
  + K \Big)  \\
B(\mu,\sigma^2) &= \int b(\sigma x + \mu)\phi(x) d x
\end{align*}
$$
where $\text{diag}(\cdot)$ returns a diagonal matrix.
Comparing the Two
Suppose that we can show that $\Lambda\rightarrow 0$ (the estimated conditional covariance matrix of the random effects tends towards zero). Then the lower bound (disregarding a determinant) tends towards:
$$
\begin{align*}
\int \exp(g(\vec u)) d\vec u &\approx
   \vec y^\top(X\vec\beta + Z\vec\mu) 
  - \vec 1^\top b(X\vec\beta + Z\vec\mu) \\ 
&\hspace{25pt}+ \vec 1^\top c(\vec y) + \frac 12 \Big( 
  \log\lvert\Sigma^{-1}\rvert 
  -\vec\mu^\top\Sigma^{-1}\vec\mu + K\Big) \\
&= g(\vec\mu) + \dots
\end{align*}
$$
where the dots do not depend on the model parameters, $\vec\beta$ and $\Sigma$. Thus, maximizing over $\vec\mu$ yields $\vec\mu\rightarrow \widehat u$. Then the only difference between the Laplace approximation and the GVA is a
$$
- \frac 12\log\lvert -g''(\widehat u)\rvert
$$
term. We have that
$$
-g''(\widehat u) = \Sigma^{-1} + Z^\top b''(X\vec\beta + Z\vec u)Z
$$
where the derivatives are with respect to $\vec\eta = X\vec\beta + Z\vec u$. This does not tend towards zero as the conditional distribution of the random effects becomes more peaked. However, still very hand wavy, it may cancel out with the
$$
\frac 12\log\lvert\Lambda\rvert = -\frac 12\log\lvert\Lambda^{-1}\rvert
$$
term we disregarded in the lower bound. The first order condition for $\Lambda$ is:
$$
\Lambda^{-1} = \Sigma^{-1} + Z^\top B^{(2)}(X\vec\beta + Z\vec\mu, \text{diag}(Z\Lambda Z^\top)Z
$$
where
$$
B^{(2)}(\mu,\sigma^2) = \int b''(\sigma x+ \mu)\phi(x) dx.
$$
Thus, if $\vec\mu \approx \widehat u$ and $\Lambda \approx 0$ then:
$$
\Lambda^{-1} \approx \Sigma^{-1} + Z^\top b''(X\vec\beta + Z\vec u)Z
$$
and the Laplace approximation and the GVA yield the same approximation of the log marginal likelihood.
Notes
Do also see the annals paper Ryan Warnick mentions.
A: There's a nice old Neural Computation paper on the relationship between the Laplace approximation and variational inference with a Gaussian proxy posterior:
http://www0.cs.ucl.ac.uk/staff/c.archambeau/publ/neco_mo09_web.pdf
In fine, the variational approximation is equivalent to requiring the Laplace approximation to hold on average, where the average is taken under the proxy posterior, as opposed to just "locally."  Thus, the mean of the proxy posterior under a Laplace approximation is the point (assuming there's only one) where the gradient of the true log-posterior is zero; whereas the mean of the proxy posterior under the variational Gaussian approximation is the point that renders the average of the gradient of the true log-posterior zero.  Similarly for the covariance matrix.
