Thank you for your questions! To answer **3**: Given a data point $x^{(i)}$, the ELBO-optimal posterior $q(z|x^{(i)})$ depends on both the prior $p(z|\theta)$ and the likelihood $p_\gamma(x^{(i)}|z)$. Because $x^{(i)}$ is a fixed value pulled from the data set, we need to think of the likelihood as a function over $z$, writing $\ell(z|x^{(i)})=p_\gamma(x^{(i)}|z)$.
To find the optimal posterior, we must integrate $p(z|\theta)$ (which is easy) and $\ell(z|x^{(i)})$ (which is impossible due to the neural network decoder). The central innovation of the SVAE is that, instead of approximating the posterior directly as in amortized VI, we approximate the likelihood function $\ell(z|x^{(i)})$ and then "combine" it with the prior. By combine, we mean: we imagine a world where our prior is still $p(z|\theta)$ but our likelihood is now $\hat{\ell}_\phi(z|x^{(i)})$, and optimize the posterior in that setting. This is the surrogate objective: the ELBO in this hypothetical world.
Let's consider the standard normal case with latent dimension 1. The prior is a $\mathcal{N}(0,1)$. Because we cannot integrate the decoder likelihood $\ell(z|x^{(i)})$ over $z$, we approximate it with a normal distribution whose mean is $\mu$ and precision (i.e. $\frac{1}{\sigma^2}$) is $\tau$. We write $\hat{\ell}_\phi(z|x^{(i)}) = \mathcal{N}(z;\mu,\tau^{-1})$, where $\mu$ and $\tau$ depend on $x^{(i)}$. In this setting, the surrogate loss is the ELBO given our $\mathcal{N}(0,1)$ prior and our $\mathcal{N}(\mu,\tau^{-1})$ (surrogate) likelihood. If the prior and likelihood are both Gaussian, then the optimal posterior (which minimizes the surrogate loss) is also Gaussian, with mean $\frac{\tau}{\tau+1}\mu$ and precision $\tau+1$. Note that in training we would never have to actually evaluate the surrogate loss: it is set up so we can find its optimizer in closed form.
We now proceed to the Gaussian mixture model, where we will assume all model parameters ($\pi,\mu,\Sigma$) are point estimates for simplicity. The *full* generative model for an SVAE with a latent Gaussian mixture (and fixed output variances) is:
$$k\sim\text{Cat}(\pi),\quad z \sim\mathcal{N}(\mu_k,\Sigma_k),\quad x\sim\mathcal{N}(\mu_\gamma(z),I)$$
Where $\mu_\gamma(z)$ is the *decoder* network. The ELBO for this model is:
$$\mathcal{L}= E_{q(k, z)}\left[ \log \frac{p(k)p(z\mid k) p_\gamma(x \mid z)}{q(k, z)}\right]\leq \log p(x)$$
The amortized inference approach would estimate $q(k,z)$ via $q_\phi(k,z\mid x)$. However, in this case computing the loss requires sampling from $q_\phi(k, z\mid x)$, which is problematic with respect to gradients as discussed elsewhere (enumeration over $k$ is possible in the GMM case, but is infeasible for time-series models).
An alternative used by several models is to make the (structured) mean-field assumption and assume that $q(k, z)$ factorizes as $q(k, z)=q(k)q(z)$.
$$\mathcal{L}=E_{q(k)q(z)}\left[ \log\frac{p(k)p(z\mid k)p_\gamma(x \mid z)}{q(k)q( z)}\right]\leq\log p(x)$$
The advantage of this factorization is that both $E_{q(k)}\left[\log p(k) \right]$ and $E_{q(k)q(z)}\left[\log p(z\mid k) \right]$ can be computed in closed-form and thus don't require sampling, while the reconstruction term: $E_{q(z)}\left[\log p_\gamma(x\mid z) \right]$ *no longer depends on $k$* and thus can be estimated with the standard reparameterization trick. As previously discussed, if the decoder took $k$ as input we would be left with $E_{q(z)q(k)}\left[\log p_\gamma(x\mid z, k) \right]$ and would be back to needing to sample $k$, thus to answer **2**, this is how we avoid sampling $k$ and it can be seen as a modeling choice. In the case of GMM, this choice does not hurt us: the continuous latent variable contains strictly more information than its cluster assignment so should be sufficient as an input to the decoder.
If we again use amortized inference for both of these distributions, so that $q(k)$ becomes $q_\phi(k\mid x)$ and $q(z)$ becomes $q_\phi(z\mid x)$, we get the *stochastic inference network* (SIN) approach. This approach performs poorly in empirical tests due to the strong decoupling of $q(k)$ and $q(x)$ ($q(k)$ is also not directly informed by the decoder).
The SVAE approach is to consider a world where the likelihood function isn't parameterized by neural networks, but is instead a simple Gaussian (the likelihood is only a function of $z$, not of $k$, so we don't need to define an approximate likelihood on $k$). GMMs with Gaussian likelihoods are well studied: in this case, optimizing the surrogate loss cannot be done in close form but can be done via *coordinate ascent variational inference* (CAVI). To answer **1**, CAVI alternates updates of $q(z)$ and $q(k)$ to optimize both using expectations, not samples. As a result, $q(k)$ can be optimized without sampling and its contributions to the loss can be evaluated without sampling, but it informs reconstruction: the parameters of $q(z)$ depend on the expectations of $q(k)$.