Let $k\in\mathbb N$ and $y_1,\ldots,y_k\in[0,1)^d$ with $$\frac1k\sum_{i=1}^kh(y_i)\approx\int_{[0,\:1)^d}h(y)\:{\rm d}y\tag1$$ for every nice enough function $h:[0,1)^d\to\mathbb R$.
Now let $p\ge1$, $f_x:[0,1)^d\to[0,1)$ for $x\in[0,1)^d$ and $g:[0,1)^d\to[0,1)$. I need to compute the integral $$\int_{[0,\:1)^d}\left|f_x(y)-g(y)\right|^p\:{\rm d}y\approx\frac1k\sum_{i=1}^k\left|f_x(y_i)-g(y_i)\right|^p\tag2$$ for various $x\in[0,1)^d$.
The problem is: For a given $x\in[0,1)^d$, the "relevant support" of $f_x$ is very small. Given a very small $\varepsilon>0$, there is a small $\delta(x)>0$ such that if $y\in[0,1)^d$, then $0\le f_x(y)<\varepsilon$ if and only if $\|x-y\|<\delta(x)$.
I store each $q(y_i)$ in an array; say $q_i$. Now, since the "relevant support" of $f_x$ is very small, the computation of $(2)$ is wasteful and it is expensive, since my $k$ in practice is very large.
My question is: What can I do here? I had to solve a similar problem where $g=0$. In that case, I've stored $y_1,\ldots,y_k$ in an Rtree qnd queried the tree for a given $x$ to give me only the (relatively small number of) $y_i$ for which the integrand $f_x$ has a "relevant support" (i.e. does not fall below my $\varepsilon$); which are simply the $y_i$ not further apart than $\delta(x)$ from $x$.
Maybe we can somehow "reweight" the points $y_i$ in a pre-computation step such that we can fall back to this simple solution.
If it is helpful, my $f_x$ is given by $f_x(y)=h(x)\exp\left(-\left(c\frac{\|x-y\|}{s(x)}\right)^2\right)$, where $h:[0,1)^d\to[0,1)$ and $s(x):=p(x)^{-\frac1d}$ (with the convention $1/0:=\infty$).
Remark: It is not possible to compute the integrals for various $x$'s in parallel, since the next $x$ I want to compute the integral for depends on the value of the integral for the former $x$. You can imagine that the $x$'s are iterates of an optimization process.
