Variational Autoencoder and Covariance Matrix Why does the encoder from a variational autoencoder map to a vector of means and a vector of standard deviations? Why does it not instead map to a vector of means and a covariance matrix?
Is it because we want our latent vectors to have zero covariance between components?
 A: The diagonal covariance matrix is an explicit statement about the kind of latent representation the researcher wants the model to learn: a representation that can be modeled as independent Gaussians.
Additionally, @Firebug points out in comments that a symmetric, PD matrix can be diagonalized without any loss of information. In other words, for some symmetric, PD matrix $A$, we can write $A=PDP^\top$ for $D$ some diagonal matrix and $P$ can be chosen to be orthonormal. This retains the same information in the sense that $A$ is rotated to have orthogonal coordinates.
Purely from the perspective of abstraction, there is no reason you must be limited to learning a latent representation that is composed of independent Gaussians. However, the computational side seems challenging.
The standard VAE encoder for a single sample produces latent parameters $(\mu, \sigma)$ its input. Then it uses the re-parameterization trick to draw random samples from that distribution. There are $d$ elements in each of $\mu$ and $\sigma$ so the total number of latent parameters is $2d$.
An alternative model that includes a covariance matrix would need some method to produce a covariance matrix, so the output of the encoder is $(\mu, \Sigma)$.
If your latent space has dimension $d$, you're doing inference on each of the $d$ elements of $\mu$ and each of the $\frac{d(d+1)}{2}$ elements of $\Sigma$ (because $\Sigma$ is symmetric by definition), for a total of $\frac{d(d+3)}{2}$ elements. Any time you have more than 1 latent dimension, the covariance matrix model will have more latent parameters to learn compared to the diagonal model.
Furthermore, the multivariate normal distribution requires that $\Sigma$ be positive definite, so we must somehow guarantee that, for each sample, we generate a PD matrix. (Using an alternative strategy, such as factorizing into standard deviations and correlation matrix $\Omega$, i.e. $\Sigma = (\sigma I) \Omega (\sigma I)$, will increase the number of effective parameters without solving the PD problem, since now we must guarantee that $\Omega$ is PD.)
Additionally, we must also be able to backprop that procedure so that the encoder weights can be updated. This may or may not be possible, depending on the strategy used to generate $\Sigma$ and draw a deviate from the multivariate Gaussian.
These three issues -- more parameters, assuring differentiability, positive definiteness -- are challenging.
If you're contemplating undertaking research to overcome these challenges, that's great! But one must ask, why is this a good model? What problems does it solve which are not solved by the diagonal Gaussian VAE model, or an alternative non-Gaussian VAE model (e.g. a Dirichlet VAE)?
A: I'm working on something like this now. 
Let's say the researcher wants a full covariance latent variable $Z ~ N(mu, \Sigma)$. Like the multiplication trick from VAE, we can have a matrix multiplication trick. Still sample from a high dimensional unit Gaussian, $e \sim N(0, I)$. 
There is this property of linear transformations of high D Gaussians: https://www.statlect.com/probability-distributions/normal-distribution-linear-combinations#hid2. If we want the transformed covariance to be $\Sigma$, we have to hit $e$ with the matrix square root of the proposal covariance $\Sigma^{1/2}e \sim N(0, \Sigma)$. 
Someone above pointed out that $\Sigma$ must be symmetric positive definite. This can be modeled as a sum of rank 1 matrices, formed as the outer product of each component vector. This is a sum of quadratic forms, so we know it will at least be positive semidefinite. To get positive definite, we have to choose a sufficient number of component vectors. By CLT (> 30 components), the off diagonal entries will be zero mean Gaussian, while the main diagonal remains positive. The number of components needed in practice might depend on other factors. 
I tried this with MNIST and got decent results. There might be a way to enforce PD by parametrizing $\Sigma$ with Inverse Wishart prior, but that also requires one of the parameters to be a symmetric PD matrix. 
A: A variational autoencoder is based on variational inference. The zero covariance is an assumption, mean-field variational family, which makes optimization easier, since the latent variables are independent.
Check chapter 2.3 of this paper: https://arxiv.org/abs/1601.00670
They also briefly dive into other families:

One way to expand the family is to add dependencies
  between the variables (Saul and Jordan, 1996; Barber and Wiegerinck, 1999); this is called structured variational inference.

