Variational inference with dependent variables Classical variational inference uses mean field theory because of its computational benefits, i.e. assumes that the latent variables are independent. For gaussian distribution, it wants to find a variational normal distribution with diagonal covariance matrix which is close enough to the standard normal distribution as the prior density. 
What is the the form of evidence lower bound (ELBO) for an especial case with dependent latent variables? As we know, in mean  field family, ELBO contains two parts, a KL divergence (generally between a normal distribution with diagonal covariance matrix and the standard normal distribution) and the reconstruction error. If we assume that our latent variables have a multivariate normal distribution with a full covariance matrix(that shows they are dependent), how does the model change? (I know what is the formula for KL divergence between normal distributions, but it seems the variational ELBO should involve some additional terms).  Thanks beforehand.
 A: The ELBO is, as usual, given by 
\begin{align}
\mathcal{L} ( q ) &= \log p ( x ) - KL ( q ( z ) || p ( z | x) ) \\
&= \int q ( z ) \log \frac{ p ( x, z ) }{ q ( z )} dz \\
&= \mathbf{E}_q \left[ \log p ( x, z ) \right] + \mathcal{H} ( q )
\end{align}
where $\mathcal{H}$ is the entropy.
Now, if $q ( z ) = \mathcal{N} ( \mu, \Sigma )$, one has that $ \mathcal{H} ( q ) = \frac{1}{2} \log \det \left( 2 \pi e \Sigma \right)$. 
Moreover, if $\xi \sim \mathcal{N} ( 0, I )$, then $\mu + \Sigma^{1/2} \xi$ is distributed as $q$. As a result, one can write 
\begin{align}
\mathbf{E}_q \left[ \log p ( x, z ) \right] = \mathbf{E}_{\xi \sim \mathcal{N} ( 0, I )} \left[ \log p ( x, \mu + \Sigma^{1/2} \xi ) \right].
\end{align}
This allows us to write the ELBO as a function of $( \mu, \Sigma )$:
\begin{align}
\mathcal{L} ( \mu, \Sigma ) &= \mathbf{E}_{\xi \sim \mathcal{N} ( 0, I )} \left[ \log p ( x, \mu + \Sigma^{1/2} \xi ) \right] + \frac{1}{2} \log \det \left( 2 \pi e \Sigma \right).
\end{align}
One can estimate gradients of the first term by noting that
\begin{align}
\nabla_{ ( \mu, \Sigma) } \mathbf{E}_{\xi \sim \mathcal{N} ( 0, I )} \left[ \log p ( x, \mu + \Sigma^{1/2} \xi ) \right] = \mathbf{E}_{\xi \sim \mathcal{N} ( 0, I )} \left[  \nabla_{ ( \mu, \Sigma) } \log p ( x, \mu + \Sigma^{1/2} \xi ) \right]
\end{align}
and thus gradients can be estimated unbiasedly by Monte Carlo. Gradients of the second term can be obtained analytically.
In closing, one can compute unbiased estimates of the gradient of $\mathcal{L} ( \mu, \Sigma )$, which allows for efficient gradient-based optimisation of the ELBO. This is roughly the idea behind Automatic Differentiation Variational Inference.
