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.
 A: 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.
