# Fast Evaluation of a Double Sum

Let

• $$q$$ be a probability distribution on $$\mathcal{X}$$,
• $$w$$ be a nonnegative function from $$\mathcal{X}$$ to $$\mathbf{R}$$ which is bounded away from $$0$$ and $$\infty$$, and
• $$s$$ be a bounded function from $$\mathcal{X}$$ to $$\mathbf{R}_+$$ with mean $$0$$ under $$q$$

I am interested in unbiased estimation of the following quantity:

\begin{align} \gamma = \int_{\mathcal{X}^N} \left( \prod_{i = 1}^N q (x_i)\right) \cdot \frac{w(x_1)}{\sum_{i = 1}^N w(x_i)} \cdot s(x_1) \, dx_{1:N} \end{align}

given $$(N + 1)$$ iid samples from $$q$$, where $$N$$ is a positive integer.

Now, given that $$s$$ has mean $$0$$ under $$q$$, I can actually rewrite this expression in terms of a covariance, i.e. define

\begin{align} W ( x_1 | x_{2:N} ) = \frac{w(x_1)}{\sum_{i = 1}^N w(x_i)}, \end{align}

then \begin{align} \gamma = \int_{\mathcal{X}^{N - 1}} \left( \prod_{i = 2}^N q (x_i)\right) \cdot \text{Cov}_{q(x_1)} \left( W ( x_1 | x_{2:N} ), s(x_1) \right) \, dx_{2:N} \end{align}

Given samples $$x_{0:N}$$, one can then form an unbiased estimator of the inner covariance as

\begin{align} \hat{c} (x_0, x_1 | x_{2:N} ) \triangleq \frac{1}{2} \left\{ W ( x_0 | x_{2:N} ) - W ( x_1 | x_{2:N} ) \right\} \cdot \left\{ s(x_0) - s(x_1) \right\}, \end{align}

and so one can form an overall estimator by re-using this estimator for all splits of the data, i.e.

\begin{align} \frac{1}{N(N+1)} \sum_{i = 0}^N \sum_{j = 0}^N \hat{c} (x_i, x_j | x_{0:N \setminus \{i, j\} } ). \end{align}

Unfortunately, naive evaluation of the above sum costs $$\mathcal{O} \left( N^2 \right)$$, which is more expensive than I would like. My question is thus: can one simplify this sum such that it can be evaluated exactly in time $$\mathcal{O} \left( N \right)$$?

One reason why I suspect that this may be possible is that $$\hat{c}$$ can be re-written as

\begin{align} \hat{c} (x_i, x_j | x_{0:N \setminus \{i, j\} } ) = \frac{1}{2} \left\{ \frac{w (x_i)}{W - w (x_j)} - \frac{w (x_j)}{W - w (x_i)} \right\} \cdot \left\{ s(x_i) - s(x_j) \right\}, \end{align}

where $$W = \sum_{i = 0}^N w (x_i)$$, which may be able to simplify some things.

• The term $W(x_1|x_{2:N})\times s(x_1)$ is an unbiased estimator of $\gamma$. – Xi'an Nov 26 '20 at 18:28
• @Xi'an indeed! My aim is not just to construct an unbiased estimator, but to construct one which i) uses the knowledge that $s$ is mean-0 to reduce variance, ii) uses the full information of all of the samples which I have generated, and iii) can be evaluated quickly. The single-term estimator satisfies only iii), as it has cost $\mathcal{O} \left( N \right)$. A mid-way point would be to use $$\frac{1}{N+1} \sum_{i = 0}^N \hat{c} (x_i, x_{i+1} | x_{0:N \setminus \{i, i + 1 \} } )$$ (interpreting indices $\mod (N+1)$), which also costs $\mathcal{O} \left( N \right)$ to form. – πr8 Nov 26 '20 at 19:13

## 1 Answer

I think that the following should work: re-write

\begin{align} \hat{c} (x_i, x_j | x_{0:N \setminus \{i, j\} } ) &= \frac{1}{2} \left\{ \frac{w (x_i)}{W - w (x_j)} - \frac{w (x_j)}{W - w (x_i)} \right\} \cdot \left\{ s(x_i) - s(x_j) \right\} \\ &= \frac{1}{2} \left\{ \frac{1}{W - w (x_j)} \cdot \left[ w(x_i) \cdot s(x_i) \right] + \frac{1}{W - w (x_i)} \cdot \left[ w(x_j) \cdot s(x_j) \right] \right\} \\ &- \frac{1}{2} \left\{ \frac{s(x_j)}{W - w (x_j)} \cdot w(x_i) + \frac{s(x_i)}{W - w (x_i)} \cdot w(x_j) \right\} \end{align}

which allows for the reformulation

\begin{align} \hat{\gamma} &\triangleq \frac{1}{N(N+1)} \sum_{i = 0}^N \sum_{j = 0}^N \hat{c} (x_i, x_j | x_{0:N \setminus \{i, j\} } ) \\ &= \frac{1}{N(N+1)} \sum_{i = 0}^N \sum_{j = 0}^N \left\{ \frac{1}{W - w (x_j)} \cdot \left[ w(x_i) \cdot s(x_i) \right] \right\} \\ &- \frac{1}{N(N+1)} \sum_{i = 0}^N \sum_{j = 0}^N \left\{ \frac{s(x_j)}{W - w (x_j)} \cdot w(x_i) \right\} \\ &= \frac{N+1}{N} \cdot \left\{ \overline{\frac{1}{W - w}} \cdot \overline{ws} - \overline{\frac{s}{W-w}} \cdot \overline{w} \right\} \end{align}

where

\begin{align} \overline{\frac{1}{W - w}} &\triangleq \frac{1}{N+1} \sum_{i = 0}^N \frac{1}{W - w (x_i) } \\ \overline{ws} &\triangleq \frac{1}{N+1} \sum_{i = 0}^N w(x_i) s(x_i) \\ \overline{\frac{s}{W-w}} &\triangleq \frac{1}{N+1} \sum_{i = 0}^N \frac{s(x_i)}{W - w(x_i)} \\ \overline{w} &\triangleq \frac{1}{N+1} \sum_{i = 0}^N w(x_i) \end{align}

which suggests that the whole estimator can be computed in time $$\approx 4N = \mathcal{O} \left( N \right)$$, noting that $$\overline{w} = \frac{W}{N + 1}$$ does not need to be re-computed.

Some additional manipulations reveal that this estimator can be written as

\begin{align} \hat{\gamma} &= \left\{ \frac{1}{N+1} \sum_{i = 0}^N \frac{1}{W - w(x_i)} \right\} \cdot \left\{ \frac{1}{N} \sum_{i = 0}^N \left( s(x_i) - \overline{s} \right) \cdot \left( w(x_i) - \overline{w} \right) \right\} \\ &- \left\{ \frac{1}{N} \sum_{i = 0}^N \left( \frac{1}{W - w(x_i)} - \overline{\frac{1}{W - w}} \right) \cdot \left( s(x_i) - \overline{s} \right) \right\} \cdot \left\{ \frac{1}{N+1} \sum_{i = 0}^N w(x_i) \right\} \end{align}

which may suggest implementations of this estimator which achieve improved numerical stability.