Let $X_i$ be independent Bernoulli random variables with success probabilities $p_i$. That is, $\mathbb{P}(X_i = 1) = p_i, \mathbb{P}(X_i = 0) = 1-p_i$. Is there a closed expression or an approximate formula for the following expectation? $$\mathbb{E}\left(\frac{\sum_i \sum_j X_i X_j B_{ij}}{(\sum_i X_i)^2}\right)$$ where $B_{ij}$ are coefficients and $B_{ij} = B_{ji}, B_{ii} = 0, \forall i, j$.
In matrix form, this is a ratio of two quadratic forms (while the latter one has a power of 2)
$$\mathbb{E}\left(\frac{\mathbf{X}^T \mathbf{B} \mathbf{X}}{(\mathbf{X}^T \mathbf{X})^2}\right)$$
where $\mathbf{B}$ is a diagonal free symmetric matrix.
Intuitively, this value is the expected average value of a submatrix of $\mathbf{B}$ where the rows (and columns) are picked randomly.
It can be evaluated in exponential time by enumerating all the scenarios. But I am wondering if there are any methods that can evaluate this expectation in polynomial time.
Note that if we can evaluate $$\mathbb{E}\left(\sum_i \sum_j X_i X_j B_{ij} \Big| \sum_i X_i = m\right)$$ then we are done since the distribution of $\sum_i X_i$ is known (Poisson Binomial Disribution)
Note too that I have already figured out the following expectation $$\mathbb{E}\left(\sum_i C_i X_i \Big| \sum_i X_i = m\right) = \frac{1}{P_k}\sum_{n=0}^N \sum_{j=1}^N \frac{A^{-nm+n}}{(N+1)}p_j C_j \prod_{\substack{k=1\\k\neq j}}^N (p_k A^n + (1-p_k))$$ $$P_k = \mathbb{P}\left(\sum_i X_i = m\right) = \frac{1}{N+1}\sum_{n=0}^N A^{-nm}\prod_{j=1}^N (p_j A^n + (1-p_j)),\quad A = e^{2i\pi/(N+1)}$$ using the method introduced in Fernandez (2010).