Consider a $\mathcal{C}^1$ function $V:\Omega\rightarrow\mathbb{R}$ where $\Omega\subset\mathbb{R}^n$. If a random vector $X$ has a parametric density $p_\theta(\textbf{x})$ that's smooth in its parameters and uniformly supported on $\Omega$ we can use the score estimator for a stochastic derivative wrt $\theta\in\mathbb{R}^m$:
$$ \nabla \mathbb{E}\left[V(X)\right] = \mathbb{E}\left[V(X)\nabla \log p_\theta(X)\right] $$
Consider now $p_\theta$ as the density for $X\sim N(f(\theta), \epsilon I)$. Then $\nabla \log p_\theta (X)= \epsilon^{-1}J_f(\theta)^\top(f(\theta)- X)$. This then makes it look like we can get a deterministic approximate derivative for $V$ without access to a gradient oracle for $V$, for small $\epsilon$:
$$ \nabla V(f(\theta))\approx\epsilon^{-1} J_f(\theta)^\top\int d\textbf{x}\;p_\theta(\textbf{x})V(\textbf{x})(f(\theta)-\textbf{x}) $$
By continuity of $\nabla V$ equality should hold in the limit. Since the integral above looks like a Gaussian convolution, perhaps this is some kind of generalization of the standard $m$-dimensional finite difference procedure for generating an approximation to $\nabla V\circ f$.
Does this have a name? If $V$ is, say, $L$-Lipschitz at $f(\theta)$, then would using a quadrature rule for the integral yield a sensible gradient approximation (sensible being as good or better than the $O(m\epsilon^2)$ 2-norm error of symmetric finite differences)?
Edit. Perhaps the integral approximation is distracting. Let's assume we can compute the exact value of the integral and ask the same question about its accuracy for a fixed $\epsilon$ (indeed, for a Monte Carlo approximation, we can make the the integral accurate with high probability by the weak LLN). Though any thoughts on the internal integral would be much-appreciated.