Let $d\in\mathbb N$, $p:[0,1)^d\to[0,\infty)$ be Lebesgue integrable, $\sigma>0$, $$\varphi(x,y):=p(x)e^{-\frac{\|x-y\|^2}{\sigma^2}}\;\;\;\text{for }x,y\in[0,1)^d$$ and $$A(x,y):=\sum_{i=1}^k\varphi(x_i,y)\;\;\;\text{for }x\in([0,1)^d)^k\text{ and }y\in[0,1)^d$$ for some $k\in\mathbb N$. Now, let $\lambda$ denote the Lebesgue measure on $\mathcal B(\mathbb R)$, $$E(x):=\int_{[0,\:1)^d}\left|A(x,\;\cdot\;)-p\right|\:{\rm d}\lambda^{\otimes d}.\tag1$$ Note that $$\nabla_{x_i}E(x)=\int_{[0,\:1)^d}\operatorname{sgn}(A(x,y)-p(y))e^{-\frac{\|x_i-y\|^2}{\sigma^2}}\left(\nabla p(x_i)+2p(x_i)\frac{y-x_i}{\sigma^2}\right)\:\lambda^{\otimes d}({\rm d}y)\tag2.$$
Question: I want to numerically compute $(1)$ and $(2)$. At the moment, I do this on a regular grid, but both the error$^1$ and the computation time is very bad when $k$ is large.
What is the "state-of-the-art" approach to compute this?
$^1$ Assume for example that each $x_i$ on the grid used for the numerical computation; then $E(x)$ will be computed as $0$; which is clearly dramatically wrong.