Fast Evaluation of a Double Sum

Question

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.

Answer 1

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.