# Subtraction of Monte Carlo integrals - Catastrophic cancellation

I am attempting to estimate a quantity $$Q$$ which is given by the difference between two functions of Monte Carlo integrals over some set of points $$\{x_i\}_{i=1}^N$$, call the estimator $$\hat{Q}$$:

$$\hat{Q} = \frac{1}{N} \sum_i f(x_i) - \left(\frac{1}{N}\sum_j g(x_j)\right)^2~.$$

Unfortunately in the cases where both terms are near-equal, catastrophic cancellation results and the final approximation is not trustworthy. Is there any way possible to avoid catastrophic cancellation in this setting?

Edit: Note this is not an attempt to compute variance, despite it's similarity!

• If this is an attempt to compute the variance of something you may be better not "simplifying" the variance (i.e. to have some quadratic in the first summation). Failing that you may be able to have a rough approximation of $\bar{g}$ (the thing inside the $^2$ in the second term) and adjust both terms accordingly, in a manner akin to $\overline{(y_i-m)^2} - (\bar{y}-m)^2$ (I hope I got that right, I didn't check the algebra) where $m$ is some very rough approximation to $\bar{y}$. Commented Jul 11 at 1:44
• Unfortunately this is not a variance computation, although it does look similar! $Q$ is not shift invariant so I suppose the proposal wouldn't work in this case... Commented Jul 11 at 1:58
• Are the $f_i = f(x_i)$ values non-negative? Commented Jul 11 at 5:55
• Yes, the $f(x_i)$ values are non-negative. The full expression looks something like $$-(A,A) + \vert(B,C)\vert^2~,$$ where $(\cdot, \cdot)$ denotes a Hermitian inner product, which is estimated via Monte Carlo. Commented Jul 11 at 14:56
• Thanks! Then my initial comment, from "failing that" on applies. While it was not a variance, the same strategies should work with the little adaptation given. If $f_i$ is a square already, its even more direct (you'd already the $y_i$ values before squaring). The posted answer gives the details. I've removed the unneeded second part and added a few additional comments including brief discussion of whuber's suggestion Commented Jul 12 at 0:45

Catastrophic cancellation will occur when $$\bar{f}=\frac{1}{n} \sum_i f(x_i)$$ is large, positive and nearly equal to $$\bar{g}^2 =\left( \frac{1}{n} \sum_i g(x_i) \right)^2$$.

If all the $$f(x_i)$$ are non-negative (which was confirmed in comments):

Let $$y_i = f(x_i)^\frac12$$. Let $$g_i=g(x_i)$$ and similarly for $$f$$.

$$\begin{eqnarray} \frac{1}{n} \sum_i (y_i-m)^2 &=&\frac{1}{n} \sum_i (y_i^2-2y_i m + m^2)\\ &=&\frac{1}{n}\sum_i y_i^2-2m\sum_i (y_i - m)\\ &=&\frac{1}{n}\sum_i f_i-2m \bar{e}_y \text{ ... }\text{ where } \bar{e}_y=\frac{1}{n}\sum_i (y_i - m)\\ \end{eqnarray}$$

$$\begin{eqnarray} \left( \frac{1}{n} \sum_i (g_i-m) \right)^2 &=&\left( \frac{1}{n} (\sum_i g_i)-m \right)^2\\ &=& (\bar{g}-m)^2 \\ &=& \bar{g}^2-2\bar{g}m+m^2\\ &=& \bar{g}^2-2m(\bar{g}-m) \\ &=& \bar{g}^2-2m\bar{e}_g \end{eqnarray}$$ where $$\bar{e}_g=\frac{1}{n}\sum_i (g_i - m)$$

Hence

$$\begin{eqnarray} \frac{1}{n} \sum_i f_i - \left( \frac{1}{n} \sum_i g_i \right)^2 &=& \frac{1}{n}\sum_i f_i-2m \bar{e}_y + 2m \bar{e}_y - (\bar{g}^2-2m\bar{e}_g) - 2m\bar{e}_g\\ &=&\frac{1}{n} \sum_i (y_i-m)^2-\left( \frac{1}{n} \sum_i (g_i-m) \right)^2 \\ & & \quad +2m (\bar{e}_y - \bar{e}_g) \end{eqnarray}$$

So $$\hat{Q}=\hat{Q}_m + 2m (\bar{e}_y - \bar{e}_g)$$

where $$\hat{Q}_m=\frac{1}{n} \sum_i (y_i-m)^2-\left( \frac{1}{n} \sum_i (g_i-m) \right)^2$$.

There's still subtraction here (both inside the summations and at the end), but if $$m$$ is reasonably well-chosen it should lead to considerably less loss of accuracy.

Advice: Choose $$m$$ in the rough ballpark of either $$\bar{y}$$ or $$\bar{g}$$, whichever is easier. They should be about the same (or the cancellation would not have been so catastrophic). If the calculation is not "on-line", average a few randomly selected $$g$$'s or average a few randomly selected $$y$$'s as is convenient; if both are convenient, you could even do both and average them.

If there's cancellation of any consequence in the difference of $$\bar{e}$$ terms a better choice of $$m$$ might be called for but you can also compute each term in the difference $$=(y_i-g_i)$$, subtract those term by term and average those instead. This is not likely to be an issue.

If you happen to have $$\bar{g}$$ before this step of calculating $$\hat{Q}$$, use that for $$m$$ and drop the remainder adjustment for the $$g$$ term (which is $$0$$).

In many cases, even just $$m=g_1$$ may suffice, so this might be used online (since $$g_1$$ will be available at the start).

Note that one could adapt Welford's algorithm to this purpose in very similar fashion and obtain a more stable online algorithm still.

whuber points out in comments an even simpler approach that doesn't rely on the individual $$f_i\geq 0$$, just $$\bar{f}$$ itself (which should always be the case if you're getting catastrophic cancellation with $$\bar{g}^2$$).

Let $$a=\bar{f}^\frac12$$. Then $$\hat{Q}=(a-\bar{g})(a+\bar{g})$$. Now that first term will still have cancellation in it, but the impact should tend to be smaller.

Having played around with it a bit, it seems that it does often help but it turns out it's not always better, if found some examples fairly quickly where sometimes it makes the error worse. It might well work fine in your application but there are numerous alternatives that might be used. I'd expect the error here mostly comes in the cancellation from $$\sqrt{\bar{f}}-\bar{g}$$ (rather than in computing the square root), and in that case you could improve the accuracy of that part if its at issue, but if you need to take that step this approach that may not offer much benefit over a more direct one.

• +1. Provided the sum of the $f_i$ is positive, it's tempting to factor $N^2\hat Q=(\sqrt{N\sum f_i})^2-\left(\sum g_i\right)^2=A^2-B^2$ as $(A+B)(A-B).$ There will be less cancellation in $A-B$ and $A+B$ is stably computed.
– whuber
Commented Jul 11 at 13:23
• Yes, that's an excellent point Commented Jul 11 at 21:58