VAE: why we do not sample again after decoding and before reconstruction loss? In many of the VAE schematics and in the original paper, a sampling step is present after decoding and before the reconstruction loss as shown in the image below. The image comes from Stanford CS321n.

In many of the code implementations though, this step is not present. For example in the Keras implementation available here: https://keras.io/examples/variational_autoencoder/
In the latent space z they sample with the Lambda layer, but at the end of the decoder there is just a Dense layer with a sigmoid activation.
Is the sigmoid doing something I don't understand mathematically? Is the VAE math still valid without this sampling step?
It is not only in code implementations, in some other schematics and textual material it seems to be ignored (see next image).

 A: The most important point stems from the confusion that the tilde $\sim$ implies a sampling operation. But $\sim$ does not imply that something is sampled, which is an algorithmic/computational operation. It indicates that something is distributed according to some distribution.
Now, when we train a VAE, we want to get gradients of the ELBO. The form of the ELBO used in VAEs is typically 
$$\mathcal{L} = \mathbb{E}_{z \sim q}\left[ \log p(x|z) \right] - \mathop{KL}\left[ q(z|x) || p(z)\right].$$
In its vanilla form, the KL of the VAE's ELBO can be computed efficiently with Monte Carlo estimates from $q$. 
The first term, the reconstruction term or likelihood term, can often be computed in closed form if $z$ is given. Especially in the two most prevalent cases–that of a Bernoulli and that of a Gaussian log likelihood.
Hence, if $x|z \sim D$ with $D$ being some tractable distribution, there is no need to sample from it, as what we are interested in is $\log p(x|z)$, which is often tractable by itself.
A: In the case of image generation, the network will usually output a reconstructed image $\hat{\mathbf{x}}$ and you will compute the reconstruction term as
\begin{equation}
    \mathcal{L}_{\text{reconstr}}(\theta) = {\lVert \mathbf{x} - f_\theta(\mathbf{x}_i,\mathbf{z}) \rVert}^2 = {\lVert \mathbf{x} - \hat{\mathbf{x}} \rVert}^2
\end{equation}
Were $f_\theta(\mathbf{x}_i,\mathbf{z})$ is the decoder.
This is the same thing that happens when you do linear regression. In normal linear regression, you estimate the mean of the likelihood with $\mathbf w^\top\mathbf x_i=\hat y_i$, then evaluate the normal likelihood in the true point $y_i$. And that's what you try to maximise. If you work out the math, and you assume $\sigma$ to be fixed, then maximizing the likelihood is equivalent to minimizing the squared error.
\begin{equation}
    \max\;\mathcal{N}\left(y_i|\hat{y}_i,\sigma\right) \Leftrightarrow \min\;{\lVert y_i - \hat{y}_i \rVert}^2
\end{equation}
Here you are doing the same thing. The decoder is generating the mean of a $N\times M$ normal (where $N\times M$ is the size of the image), and you are essentially doing
\begin{equation}
    \max\;\mathcal{N}\left(\mathbf x_i|f_\theta(\mathbf{x}_i,\mathbf{z}),\mathbf{\sigma}\right) \Leftrightarrow \min\;{\lVert \mathbf x_i - \hat{\mathbf x}_i \rVert}^2
\end{equation}
while assuming that $\mathbf{\sigma}$ is fixed. Technically, you could also make the network generate $\mathbf{\sigma}$, and then optimize a more refined L2 loss that also includes $\mathbf{\sigma}$ somewhere. That way, your network would generate all parameters of a normal distribution, not just the mean, and for each latent vector $\mathbf z$ you could actually sample several synthetic images $\hat{\mathbf x}$.
However, we are often happy with generating one image for each $\mathbf z$. Thus, it is okay if we only take the most likely one, i.e. the mean of the resulting predictive posterior. After all, this is also what we do with point estimation in linear regression.
