In a Monte Carlo approximation of a product of expectations, can the same samples be used for both expectations? I am trying to approximate a product of expectations:
$\operatorname{E}[f(x)]\operatorname{E}[g(x)]=\sum_x P(x) f(x) \sum_x P(x) g(x)$
with $N$ Monte Carlo samples $(x_1,x_2,...,x_N)$ from $P(x)$:
$\operatorname{E}[f(x)]\operatorname{E}[g(x)] \approx \left( \frac{1}{N} \sum_i f(x_i) \right) \left( \frac{1}{N} \sum_i g(x_i) \right)$
The Monte Carlo samples are very expensive to generate. What are the consequences of using only $N$ samples total, reusing the same samples for each expectation (vs. using $2N$ samples so that each expectation uses a different set of $N$ samples)? Will that produce an incorrect/biased estimate, higher variance, etc.?
More Context (added Aug 2, 2013)
Specifically, the application that inspired this question concerns gradient descent optimization of entropy in Boltzmann machines. Assume $M$ visible/observable binary neurons (Bernoulli variables), and no hidden neurons for simplicity. The state of each neuron is $X_i \in \{0,1\}$, and the joint state of all neurons is $X \in \{0,1\}^M$. These neurons are fully connected with a set of ${M \choose 2}$ weights $\theta_{ij}$ between neurons $i$ and $j$ with $i \neq j$ (a symmetric $M \times M$ parameter matrix with $0$s on the diagonal). Assume the following relationships:
$$ p(\left. X \right| \theta) = \frac{1}{Z_{\theta}} \exp \left(-\frac{E(X,\theta)}{T} \right)$$
$$ Z_{\theta} = \sum_X \exp \left( -\frac{E(X,\theta)}{T} \right) $$
$$ E(X,\theta) = -\sum_{ij} \theta_{ij} X_i X_j$$
where $Z_{\theta}$ is the partition function, $E(X,\theta)$ is the energy function, and $T$ is the temperature (a global weight scalar which is usually, but not always, assumed to be $1$ for Boltzmann machines; note that Boltzmann's constant $k_B$ is assumed to be $1$ here and is not shown).
For this model, say we want to maximize the entropy of $X$ by adjusting the parameters $\theta_{ij}$. We can compute the gradient of the entropy expression and apply a gradient descent procedure to maximize it.
The entropy expression is:
$$ S_{\theta} = -\sum_X p(\left. X \right| \theta) \log p(\left. X \right| \theta)$$
The partial derivative of the entropy wrt each weight $\theta_{ij}$ is (derivation not shown):
$$ \frac{\partial}{\partial \theta_{ij}} S_{\theta}
= \frac{1}{T^2} \left[ \sum_X p(\left. X \right| \theta) \, E(X,\theta) X_i X_j - \sum_X p(\left. X \right| \theta) \, E(X,\theta) \sum_X p(\left. X \right| \theta) \, X_i X_j \right] $$
These expectations cannot generally be computed exactly if $M$ is large because they require summing over $2^M$ possible states. Instead, we just approximate them with $N$ Monte Carlo samples from $p(\left. X \right| \theta)$:
$$ \frac{\partial}{\partial \theta_{ij}} S_{\theta}
\approx \frac{1}{T^2} \left[ \left( \frac{1}{N} \sum_k^N E(X^k,\theta) X_i^k X_j^k \right) - \left( \frac{1}{N} \sum_k^N E(X^k,\theta) \right) \left( \frac{1}{N} \sum_k^N X_i^k X_j^k \right) \right] $$
where $X^k$ and $X_i^k$ refer to the $k$th Monte Carlo sample. Note that optimizing entropy alone would simply make all $\theta_{ij}=0$ because that maximizes $S_{\theta}$; however, in practice we may want to combine this objective with something else, like a maximum likelihood objective, and optimize a weighted combination of the two objectives. In that case, the max likelihood objective would make the model fit the data, and the max entropy objective would regularize the model to avoid overfitting.
Regarding the original question posted here, I was curious if this gradient approximation could use only $N$ Monte Carlo samples total (i.e. can the product of the two expectations use the same $N$ samples?), or whether it requires $2N$ samples (one set of $N$ samples for each expectation in the product) to remain an unbiased estimator of the gradient.
Note that, in any case, the expectation on the left side of the gradient expression can reuse the $N$ samples from either one of the two expectations on the right side without bias because we can simply rearrange things to put a summation outside the outer brackets. Specifically, we can do either this:
$$ \frac{\partial}{\partial \theta_{ij}} S_{\theta}
\approx \frac{1}{T^2} \frac{1}{N} \sum_k^N E(X^k,\theta) \left( X_i^k X_j^k - \frac{1}{N} \sum_l^N X_i^l X_j^l \right) $$
or this:
$$ \frac{\partial}{\partial \theta_{ij}} S_{\theta}
\approx \frac{1}{T^2} \frac{1}{N} \sum_k^N X_i^k X_j^k \left( E(X^k,\theta) - \frac{1}{N} \sum_l^N E(X^l,\theta) \right) $$
...reducing the problem from 3 to 2 expectations. Again, the question is (essentially) whether it is unwise to try to reuse the same set of $N$ Monte Carlo samples for these two expectations.
 A: Expanding on comments a little bit, I believe that you should use all the data for both expectations.
Consistency: This is completely unaffected by using the same samples by the continuous mapping theorem. Sensible enough that, in the big, it doesn't matter if you reuse samples, so there is some assurance that if $N$ is big enough you should be okay.
Variance: Let $\sigma^2_f = \mbox{Var}(f(X))$, $\sigma^2_g = \mbox{Var}(g(X))$ and $\sigma_{fg} = \mbox{Cov}(f(X), g(X))$, and $\mu_f$ and $\mu_g$ defined similarly. The variance of your estimator with reused samples can be shown to be
$$
\frac{\mu_g^2\sigma_f^2 + 2\mu_f\mu_g\sigma_{fg} + \mu_f^2\sigma_g^2}{N} + O(N^{-2}),
$$
where $N$ is the number observations. If we run $2N$ iterations and use $N$ to estimate one and the other we get 
$$
\frac{\mu_g^2\sigma^2_f + \mu^2_f\sigma^2_g}{N} + O(N^{-2}),
$$
instead. This is a simple application of the multivariate delta method, and while you can perhaps get some insight using exact calculations, this seems good enough for intuition. What we see is that, in fact, it is possible for reusing samples to decrease the variance if $\mu_f\mu_g\sigma_{fg} < 0$. That is, we might do half the work reusing data (only $N$ versus $2N$ samples) and get a better result.
To see more clearly we should use all the data, consider the following. The variance using $2N$ samples for both is to $O(N^{-1})$
$$
\frac{\mu_g^2\sigma_f^2 + 2\mu_f\mu_g\sigma_{fg} + \mu_f^2\sigma_g^2}{2N} \le \frac{\mu_g^2\sigma_f^2 + 2|\mu_f\mu_g\sigma_f\sigma_g| + \mu_f^2\sigma_g^2}{2N},
$$
by Cauchy-Schwartz ($|\sigma_{fg}| \le \sigma_f\sigma_g$). But
$$
\frac{\mu_g^2\sigma^2_f + \mu^2_f\sigma^2_g}{N} - \frac{\mu_g^2\sigma_f^2 + 2|\mu_f\mu_g\sigma_f\sigma_g| + \mu_f^2\sigma_g^2}{2N} = \frac{(|\mu_f\sigma_g| - |\mu_g\sigma_f|)^2}{2N} \ge 0.
$$
This says the following: we can either draw $2N$ samples and use all $2N$ in both calculations or we can use $N$ for $\mu_f$ and $N$ for $\mu_g$. We have just proved that - at least to $O(N^{-1})$ - the second strategy cannot possibly be better than the first. 
Bias: As a general rule, the bias of many reasonable-looking estimators is $O(N^{-1})$ where as the standard error is typically $O(N^{-1/2})$. So the bias term is not the important term, for large enough $N$. Of course, there are ways of mucking this up. Doing a little bit more, it can be shown that
$$
E(\bar f \bar g) = \frac{N-1}{N}\mu_f \mu_g + \frac{1}{N^2}\sum_{i = 1} ^ N E(f(X_i)g(X_i)) = \mu_f\mu_g + \frac{\sigma_{fg}}{N}.
$$
So the bias is $O(N^{-1})$ and contributes $O(N^{-2})$ to the mean squared error, which we already have up to $O(N^{-2})$ to begin with. But, at any rate, the magnitude of the bias is determined by $\sigma_{fg}$, which you can estimate using your sample and satisfy yourself that it isn't too large. 
