In this paper (bottom of p. 4), the authors state the following (I've added all definitions below):
Given the deterministic mapping $z=g_{\phi}(\epsilon, x)$ we know that: $q_\phi (z|x)\Pi_i dz_i = p(\epsilon)\Pi_i d\epsilon_i$. Therefore, $$\int q_\phi(z|x)f(z)dz = \int p(\epsilon)f(z)d\epsilon=\int p(\epsilon)f(g_\theta(\epsilon, x))d\epsilon$$
where:
$z$ is a continuous random variable from some conditional distribution: $z \sim q_\phi(z|x)$
$z$ can be expressed using a deterministic, vector-valued function $g_\phi(\epsilon, x)$ parameterized by $\phi$ using $\epsilon \sim p(\epsilon)$. For example, for $z \sim N(\mu(x), \sigma(x)^2)$, we could write $z = \mu(x) + \sigma(x)\epsilon$ where $\epsilon \sim N(0,1)$.
$f$ is some function.
I can't derive this formally (I can intuitively see how they plugged in the equality before the "Therefore," and that it 'kind of makes sense'). I was trying to see if I could somehow prove this using change of variables, but since $\epsilon$ is itself sampled, I don't see how that applies.
What is going on here, formally? (Ideally - I would like to know which theorems are being used and how).
This question was addressed here, but not rigorously.
I used some better notation. These are assumed:
$$f_Z(z)= q_\phi(z|x)$$ $$Z \sim g_{\phi}(\mathcal E, x)$$ $$f_{\mathcal E}(\epsilon)= p(\epsilon)$$
where $f_Z(z)$ and $f_{\mathcal E}(\epsilon)$ denote the pdfs of random variables $Z$ and $\mathcal E$, respectively. To avoid confusion, I here replace $f$ in the OP by $H$.
Then,
$$\int q_\phi(z|x)H(z)dz=\int f_Z(z)H(z)dz=\mathbb E [H(Z)]\\ =\mathbb E \big [H(\color{blue}{g_{\phi}(\mathcal E, x)}) \big ]=\int f_{\mathcal E}(\epsilon)H(g_\phi(\epsilon, x))d\epsilon=\int p(\epsilon)H(g_\phi(\epsilon, x))d\epsilon.$$