Can a variational autoencoder be interpreted as a mixture of Gaussians? In a variational autoencoder (VAE), we have an encoder network $E_{\phi}$ that maps inputs $x$ to the distribution parameters of the approximate posterior $q_{\phi}(z \vert x)$. Most commonly we model this distribution over latent as a diagonal-covariance Gaussian, so we have
$$ (\mu_{\phi}(x), \Sigma_{\phi}(x)) = E_{\phi}(x) $$
Where $\Sigma_{\phi}(x)$ are the diagonal elements of the covariance matrix corresponding to datapoint $x$. This results in a different set of distribution parameters for each datapoint $x$, where the computational cost of learning the parameters of the posterior are 'amortized' through learning the model parameters of the encoder $E_{\phi}$.
In view of this, can the distribution over the latent space produced by the VAE be thought of as a very large mixture of Gaussians with a number of components equal to the number of data points we have? i.e. for $N$ data points $\{x_1, x_2, \ldots x_N\}$, do we have
$$ q_{\phi}(z \vert x) = \frac{1}{N} \sum_{n=1}^N \delta(x,x_n) \mathcal{N}\left(\mu_{\phi}(x_n), \Sigma_{\phi}(x_n)\right)$$
In other words, can the latent variable $Z$ be modelled as a mixture of $N$ equally weighted Gaussian components, with component $n$ having distribution $\mathcal{N}\left(\mu_{\phi}(x_n), \Sigma_{\phi}(x_n)\right)$?
 A: In the VAE scenario, the encoder $q(z|x)$ already works for any $x$ (any $x$ that the encoder network has seen during training, of course), so there's no gain in representing it as a mixture of trainset distributions $q(z|x_n)$ with binary weights.
However, look at the aggregated posterior, $$ q^\text{agg}(z) = \frac{1}{N} \sum_{n=1}^N q(z|x_n) $$It is a finite (and non-degenerate) mixture of Gaussians, representing average encoding distribution. In a sense, this is the distribution you should use to sample $z$ from after you trained the model (and not the prior $p(z)$). This is because the decoder network was effectively trained on samples from $q^\text{agg}(z)$ and works best on them. In theory, Variational Inference tries to make the aggregated posterior $q^\text{agg}(z)$ and the prior $p(z)$ as close as possible, but in practice it might not succeed.
For more, see the paper on VampPrior.

Leaving VAE's inference model aside, you can claim that the generative part of the VAE is actually a (possibly infinite) mixture. Indeed, the most popular choice for the decoding distribution $p(x|z)$ is Gaussian, and $p(z)$ is also typically Gaussian. Thus,$$ p(x) = \int p(x|z) p(z) dz = \int \mathcal{N}(x \mid \mu_x(z), \Sigma_x(z)) \mathcal{N}(z \mid 0, I) dz $$So effectively here we construct a distribution $\mathcal{N}(x \mid \mu_x(z), \Sigma_x(z))$ for every $z \in \mathbb{R}^d$ and then mix all these distributions (uncountably many!) with weights of the standard multivariate Gaussian distribution. Such a mixture is too complicated to work with directly, therefore we resort to Variational Inference.
A: You seem to be mixing up two components together.
The latent $Z$ of a standard VAE is effectively a mixture of Gaussians, but it has nothing to do with the number of datapoints.
The diagonal covariance is effectively a notation trick; a VAE latent space comprises $N$ instances of independent Gaussian 'particles', where $N$ is a model hyperparameter. 
Each of those particles has a mean and variance encoded by the output of the... well, encoder. The encoder can have any number of layers in its own right and tries to figure out an embedding that can be described nicely using a bunch of Gaussians.
