Intractability in Variational Autoencoders
I believe that you got bogged down by this thought:


  I understand that it is not possible to analytically compute this integral to optimize for $\theta$.


From the paper I understood that:

$$
p_\theta (x) = \int p_{\theta} (z) p_{\theta} (x|z) dz
$$

Can be analytically solved if we simplify $p_{\theta} (x|z)$.  i.e. in the case where the Bayes' Rule for its posterior is easy.  This integral is instead intractable if the Bayes Rule for $p_{\theta} (x|z)$ turn out to be a very complex function, e.g. $(z|x)$ is very complex.  In general integrating a simple distribution is easy, conditional distributions can be (very) complex to integrate.

The reparametization trick on the other hand is what makes the paper's point.  They do not try to compute the the integrals as you derived them but they attempt to substitute all $(z|x)$ with a $g_{\phi}$ function (parametrized by $\phi$) which will build $z$ from $x$ and $\epsilon$.

The idea is that we choose $p(\epsilon)$ and instead of estimating $q_{\theta}(z|x)$ we estimate $\phi$ (in $g_{\phi}(z_{parametrized}, x$).  Estimating (SGD) a value for a parametrization of a vector function is not difficult.



I understood (analogically and very lousily) the reparametrization they do in not a very different way than what we do with, say, integrals of sphere or cylinder-like surfaces: we reparametrize the variables and recalculate the interval. Here we want the interval to not have the two sums of the KL (including both distributions) but one sum.