the actual computation factorizes into $\nabla_z p(x, z_m)$ that depends only on $x$ and $z_m$ and into $\nabla_\lambda f(\epsilon, \lambda)$ and $\nabla_{\lambda}|\text{det}_{J_f}(\epsilon_m, \lambda)|$ that depend only on $\epsilon_m$ and $\lambda$. The former part is typically considerably more computationally expensive. The observation that the slower part does not directly depend on the parameters hints that it should be possible to avoid re-computing the term even if the approximation changes, and this indeed is the case as explained next.

Assume we have already estimated the gradient at some parameters $\lambda$, implying that we have also evaluated $\nabla_z p(x, z_m)$ for some set of $z_m$. The question now is how to estimate the gradient at parameters $\lambda' = \lambda + \delta$ that are (typically only slightly) different. It turns out this can be done using the well-known concept of *importance sampling* originally designed for approximating expectations when we cannot directly draw samples from the density of interest. In our case, however, we *could* draw samples directly from the new approximation, but *choose not to* since the gradient can be estimated also using the old approximation as a proposal distribution. That is, we are using importance sampling for an unusual reason but can still 