How to calculate convolution with logarithm numerically?

535 Views Asked by At

I'm trying to compute an optimisation problem, which has a cost function involving

$$I=\int_0^1\log|x-y|\rho(y)dy$$

where $x\in[0,1]$ and $\rho$ is a probability density. Eventually, I will want to calculate the expected value of I, $\int I\rho(x)dx$.

I want to calculate this integral numerically, for example by taking the Riemann sum with spacing $h$. The problem is that the log term blows up at $x=y$. Is there a standard way to renormalise this?

I was thinking of simply ignoring the $-\infty$ terms in the Riemann sum, but that doesn't look like the right solution.

Or perhaps writing $\log[x(1-y/x)]=\log x-\sum_{n=1}^\infty \frac{y/x}{n}$ and ignoring terms n>K some K. This seems to be a bit more systematic.

Could someone help?

Thank you!

1

There are 1 best solutions below

0
On BEST ANSWER

Idea: given a small "tuning parameter" $h>0$, consider

$$I=\int_0^{x-h} \log |x-y| \rho(y) dy + \int_{x-h}^{x+h} \log |x-y| \rho(y) dy + \int_{x+h}^1 \log |x-y| \rho(y) dy.$$

Now approximate $\rho$ on $[x-h,x+h]$ by $\rho(x)$. Then the first and last integrals have no blowup problems, while the middle integral can be analytically approximated as

$$\int_{x-h}^{x+h} \log |x-y| \rho(x) dy = 2 \rho(x) \int_0^h \log(y) dy = 2 \rho(x) (h \log(h) - h).$$

Then to satisfy a tolerance of $\varepsilon$, you can choose $h$ small enough that the error in this approximation is at most $\varepsilon/2$. A simple way to approach this analysis would be:

$$\| \log(|x-y|) (\rho(y)-\rho(x) \|_{L^1([x-h,x+h])} \leq (2h - 2 h \log(h)) \| \rho(y) - \rho(x) \|_{L^\infty([x-h,x+h])}.$$

Thus everything is OK for small enough $h$ even if $\rho$ is just bounded, but you can take larger $h$ if $\rho$ has some regularity. You might get a better bound by using Cauchy-Schwarz (I haven't worked it out).

Having chosen $h$ sufficiently small, perform your ordinary quadrature sufficiently accurately on the other two intervals to have a total error of at most $\varepsilon/2$.