Given a function $f(x):R^{n} \to R$ where $x \in R^n$, and $p(x)=\mathcal{N}(x|m,V); V=diag(v)$ is a multivariate Gaussian distribution with diagonal covariance matrix $V$. We want to compute $E_{p}f(x)$ based on $m,V$ and function $f$.
One possible way is to use Taylor expansion around $m$, i.e $f(x)=f(m)+\nabla f(x)_{x=m}(x-m) $, then compute:$E_{p}f(x) \approx \int \big( f(m)+\nabla f(x)_{x=m}(x-m) \big) p(x)dx$.
However, Taylor expansion is accurate when $x$ is close to $m$. When taking expectation, each dimension of $x$ runs from $-\infty $ to $\infty$, leading $x$ is probably far from $m$.
How to estimate $E_p f(x)$ more accurate and we still obtain a parametric form (the estimation is a function of $m$ and $V$) ?
You can use Laplace's method, an introduction here: https://en.wikipedia.org/wiki/Laplace%27s_method
With your notation, $f$ is a function from $\mathbb{R}^n \mapsto \mathbb{R}$, and I will assume it is a positive function.
$$ E_p f(X) = \idotsint_{\mathbb{R}^n} e^{-g(x)} p(x) \; dx $$ where $g(x)=-\log f(x)$ (so here we do need that $f$ is positive ...). Then inserting an expression for the multivariate normal density $p(x)$, and taylor-expanding $g$, around the minimizer of $g$, $x_0$: $$ g(x)=g(x_0)+g'(x_0)'(x-x_0) + \frac12 (x-x_0)'g''(x_0)(x-x_0) + \dots $$ (and then forgetting about the dots for the approximation). This becomes $$ E_p f(X) \approx \idotsint_{\mathbb{R}^n} \exp(-g(x_0)-g'(x_0)'(x-x_0)-\frac12 (x-x_0)'g''(x_0) (x-x_0) ) \frac1{(2\pi)^{n/2} |V|^{(1/2)}} \exp(-\frac12 (x-v)'V^{-1} (x-v)) \; dx $$ here $g'$ is the gradient, $g''$ is the hessian matrix, and elsewhere ' denotes differentiation. Putting together the two exponentials, this becomes $$ \frac1{(2\pi)^{n/2} |V|^{(1/2)}} \idotsint_{\mathbb{R}^n} \exp(-g(x_0)-g'(x_0)'(x-x_0)-\frac12 (x-x_0)'g''(x_0) (x-x_0) -\frac12 (x-v)'V^{-1} (x-v)) \; dx $$ Now, to evaluate the integral we need to "complete the square": In the following, we concentrate only on the argument of the exponential in the integral. ultiply out everything in the argument, and write the resulting form as $$d+a'x -\frac12 x'(V^{-1}+g''(x_0))x $$ and then require to write it in the form $$ -\frac12 (x-\mu)'(V^{-1}+g''(x_0)) (x-\mu) $$ we find that $d=-\frac12 \mu' (V^{-1}+g''(x_0)) \mu$, $a' = -\mu' (V^{-1}+g''(x_0)) $ by solving that we find that $\mu= -a (V^{-1}+g''(x_0))^{-1}$. Then, our integral has become $$ \frac1{(2\pi)^{n/2} |V|^{(1/2)}} \idotsint_{\mathbb{R}^n} \exp( -\frac12 (x-\mu)' A^{-1} (x-\mu) ) \; dx $$ where we have introduced $A= (V^{-1}+g''(x_0))^{-1}$. Now, that has the form as the integral of a multivariate normal density, so can be found to be $$ \frac1{(2\pi)^{n/2} |V|^{(1/2)}} \cdot (2\pi)^{n/2} |A|^{1/2} $$ which I leave for you to simplify further.