# Variance of summation of products of all pairs of Bernoulli random variables

I have $$n$$ independent Bernoulli random variables with parameter $$p$$: $$\{X_1, \dots, X_n\}$$ and I want to find the variance of the sum of the products of all pairs of these random variables. Specifically, I want to calculate

$$Var\left[\sum_{ij \in [n]^2} (c_i X_i \cdot X_j)\right]$$ where $$c_{ij} \geq 0$$ are constants.

Even a somewhat tight upper bound (tighter than $$n^2$$) is also helpful.

For simplicity, we symmetrize $$(c_{ij})$$ by letting $$s_{ij} = \frac{1}{2}(c_{ij} + c_{ji})$$. Also, define $$Y$$ by

$$Y = \sum_{(i,j) \in [n]^2} c_{ij} X_i X_j = \sum_{(i,j) \in [n]^2} s_{ij} X_i X_j .$$

Then

\begin{align*} \mathbf{Var}(Y) = \mathbf{Cov}(Y, Y) = \sum_{(i,j), (k,l) \in [n]^2} s_{ij}s_{kl} \mathbf{Cov}(X_i X_j, X_k X_l). \end{align*}

Now by noting that $$X_i^2 = X_i$$, we have

\begin{align*} \mathbf{Cov}(X_i X_j, X_k X_l) &= \mathbf{E}[X_i X_j X_k X_l] - \mathbf{E}[X_i X_j] \mathbf{E}[X_k X_l] \\ &= p^{|\{i, j, k, l\}|} - p^{|\{i, j\}|}p^{|\{k, l\}|}. \end{align*}

Using the idea as in this answer, we can classify all the different scenarios for the values of $$|\{i, j, k, l\}|$$, $$|\{i, j\}|$$, and $$|\{k, l\}|$$. Consequently,

\begin{align*} \mathbf{Var}(Y) &= \Biggl( \sum_{i} s_{ii}^2 \Biggr) p(1-p) + 4 \Biggl( \sum_{i, j \text{ distinct}} s_{ii}s_{ij} \Biggr) p^2(1-p) \\ &\quad + 2 \Biggl( \sum_{i, j \text{ distinct}} s_{ij}^2 \Biggr) p^2(1-p^2) + 4 \Biggl( \sum_{i, j, k \text{ distinct}} s_{ij}s_{jk} \Biggr) p^3(1 - p). \end{align*}

Let $$M = \max_{(i,j) \in [n]^2} s_{ij}$$. Then an upper bound can be obtained:

\begin{align*} \mathbf{Var}(Y) &\leq M^2 \biggl[ \begin{gathered} np(1-p) + 4n(n-1) p^2(1-p) \\ + 2n(n-1)p^2(1-p^2) + 4n(n-1)(n-2) p^3(1-p) \end{gathered} \biggr] \\ &= M^2 np(1-p) [(2n-2)(2n-3)p^2 + 6(n-1)p + 1] \end{align*}

Note that this bound is also tight, in the sense that if $$(s_{ij})$$ is a constant array, then the equally holds for the above inequality. In particular, $$\mathbf{Var}(Y)$$ is in general of order $$n^3$$.

Below is a comparison of the values of $$\mathbf{Var}(Y)$$ by a direct computation and using the formula. Mathematica 13 is used to produce this.

n = 4;
(* Direct computation *)
dist = Table[x[k] \[Distributed] BernoulliDistribution[p], {k, n}];
Sum[Expectation[x[i]*x[j]*x[k]*x[l], dist] - Expectation[x[i]*x[j], dist] Expectation[x[k]*x[l], dist], {i, n}, {j, n}, {k, n}, {l, n}]
(* Using formula *)
n p (1 - p) + 4 n (n - 1) p^2 (1 - p) + 2 n (n - 1) p^2 (1 - p^2) + 4 n (n - 1) (n - 2) p^3 (1 - p) // Expand
Clear[n, dist];