# Handling cases of cross terms for multi-sums

I have expressions consisting of many multi-sums and I would like to extract cross terms out of them. Consider a simple example:

$$\sum_{m_1=1}^M \sum_{m_2=1}^M \sum_{m_3=1}^M \sum_{m_4=1}^M (x_{m_1}^2 x_{m_2}^{2p} - x_{m_1}^{p+1} x_{m_2}^{p+1}) (x_{m_3}^2 x_{m_4}^{2p} - x_{m_3}^{p+1} x_{m_4}^{p+1})$$

In the end, I want to wrap the expectation operator around it and simplify the result with the usual fact that $$\mathbb{E}(m_i m_j)=\mathbb{E}(m_i)\mathbb{E}(m_j)=0$$ if $$i\neq j$$ and the random variables have zero mean.

Hence I need to find all possible cases for which $$m_1, \cdots, m_4$$ are equal. For the example above it took me two hours to extract 15 cases by hand! For example the two trivial cases are $$m_1=m_2=m_3=m_4$$ and $$m_1\neq m_2\neq m_3\neq m_4$$ but another is also $$(m_1=m_4) \neq (m_2=m_3)$$ and two out of the 15 terms would be:

$$(x_{m}^2 x_{m}^{2p} - x_{m}^{p+1} x_{m}^{p+1}) (x_{m}^2 x_{m}^{2p} - x_{m}^{p+1} x_{m}^{p+1})\,\quad 1\times \\ (x_{m_1}^2 x_{m_1}^{2p} - x_{m_1}^{p+1} x_{m_1}^{p+1}) (x_{m_3}^2 x_{m_3}^{2p} - x_{m_3}^{p+1} x_{m_3}^{p+1}) \,\quad 3M(M-1)\times$$

The suggestions given in the answers work already nicely with this sum (and slightly more complicated ones). However, the actual case is rather more complicated. It results from the determinant of a $$2\times 2$$ Gramian matrix which I can write as:

$$D = \left(\sum_{m=1}^M \frac{1}{R} \left[\sum _{r=1}^R x^{p_1}_{(m-1) R+r}\right]^2\right) \frac{1}{R}\sum_{m=1}^M \left[\sum _{r=1}^R x^{p_2}_{(m-1) R+r}\right]^2 - \left(\frac{1}{R}\sum_{m=1}^M \left(\sum_{r=1}^R x^{p_1}_{(m-1) R+r}\right) \sum_{r=1}^R x^{p_2}_{(m-1) R+r}\right)^2$$

Calculating $$\mathbb{E} D$$ is still tractable but I need $$\mathbb{E}(D^2)$$. Assuming $$M$$ and $$R$$ have the same order the number of sums is 12! For just $$R=M=5$$ this results in $$5^12=244140625$$ terms! You may note that $$D^2$$ reduces to my simple example above when $$R=1$$. The problem is that the additional sum with $$R>1$$ makes the terms asymmetric. With the suggested solutions I found already that I would end up with 70 terms for $$\mathbb{E}(D^2)$$ but since I can only use very small $$R$$ and $$M$$ I am not able to find any patterns.

If there is a method with which this works I would be very happy if I could also get $$\mathbb{E}(D^q), q>2$$ but for now I am happy with $$q=2$$.

– JimB
Commented Jan 27, 2019 at 6:52
• While this can certainly be done in general in Mathematica, I bet that @wolfies MathStatica would do this in a pretty straight forward manner.
– JimB
Commented Jan 27, 2019 at 19:23
• Is there a subscript missing that goes with $x^{p_1}$ and $x^{p_2}$ ?
– JimB
Commented Jan 30, 2019 at 1:01
• No, I just use a different notation since the index gets long: $x^{p_1}[(m-1)R+r] = x^{p_1}_{(m-1)R+r}$ (this notation is used in signal processing). With $R=1$, $x^{p_1}[(m-1)1+1] = x^{p_1}[m] = x_m^{p_1}$. EDIT: I changed the notation for consistency.
– divB
Commented Jan 30, 2019 at 1:37
• FYI: I haven't forgotten to supply what I hope is a general solution but other work has come up and I'll be able to get to this next week.
– JimB
Commented Feb 1, 2019 at 14:59

Correction

As pointed out by @divB I forgot to expand the resulting sum prior to substituting in the expectations. Here is a correction.

If the $$x_i$$ with $$i=1,2,\ldots,m$$ are all independent random variables with finite moments, then one can figure out what the expectation is in terms of the raw moments: $$\mu_k$$.

(And there certainly must be a more direct way to do this.)

The sum that you mention is

Sum[(x[m1]^2 x[m2]^(2 p) - x[m1]^(p + 1) x[m2]^(p + 1))*
(x[m3]^2 x[m4]^(2 p) - x[m3]^(p + 1) x[m4]^(p + 1)),
{m1, 1, m}, {m2, 1, m}, {m3, 1, m}, {m4, 1, m}]


which returns

$$\sum _{m_1=1}^m \sum _{m_2=1}^m \sum _{m_3=1}^m \sum _{m_4=1}^m \left(x_{m_1}^2 x_{m_2}^{2 p}-x_{m_1}^{p+1} x_{m_2}^{p+1}\right) \left(x_{m_3}^2 x_{m_4}^{2 p}-x_{m_3}^{p+1} x_{m_4}^{p+1}\right)$$

Now we'll expand the sum for specific values of $$m$$ and substitute $$\mu_k$$ for $$x_{m_i}^k$$ (as that step is equivalent to taking the expectation) and look for a pattern:

Table[{m, Expand[Sum[(x[m1]^2 x[m2]^(2 p) - x[m1]^(p + 1) x[m2]^(p + 1)) (x[m3]^2 x[m4]^(2 p) - x[m3]^(p + 1) x[m4]^(p + 1)),
{m1, 1, m}, {m2, 1, m}, {m3, 1, m}, {m4, 1, m}]]} /. x[_]^k_ -> μ[k],
{m, 2, 10}] /. {μ[2] -> Subscript[μ, 2], μ[4] -> Subscript[μ, 4],
μ[2 p] -> Subscript[μ, 2 p], μ[4 p] -> Subscript[μ, 4 p], μ[1 + p] -> Subscript[μ, 1 + p],
μ[3 + p] -> Subscript[μ, 3 + p], μ[2 + 2 p] -> Subscript[μ, 2 (1 + p)],
μ[1 + 3 p] -> Subscript[μ, 1 + 3 p]}


$$\begin{array}{cc} 2 & 6 \mu _{2 (p+1)}^2+2 \mu _4 \mu _{4 p}-8 \mu _{p+3} \mu _{3 p+1} \\ 3 & 6 \mu _{4 p} \mu _2^2+12 \mu _{2 p} \mu _{2 (p+1)} \mu _2-24 \mu _{p+1} \mu _{3 p+1} \mu _2+6 \mu _4 \mu _{2 p}^2+18 \mu _{2 (p+1)}^2+6 \mu _4 \mu _{4 p}+24 \mu _{p+1}^2 \mu _{2 (p+1)}-24 \mu _{2 p} \mu _{p+1} \mu _{p+3}-24 \mu _{p+3} \mu _{3 p+1} \\ 4 & 24 \mu _{p+1}^4-48 \mu _2 \mu _{2 p} \mu _{p+1}^2+96 \mu _{2 (p+1)} \mu _{p+1}^2-96 \mu _{2 p} \mu _{p+3} \mu _{p+1}-96 \mu _2 \mu _{3 p+1} \mu _{p+1}+24 \mu _2^2 \mu _{2 p}^2+24 \mu _4 \mu _{2 p}^2+36 \mu _{2 (p+1)}^2+24 \mu _2^2 \mu _{4 p}+12 \mu _4 \mu _{4 p}+48 \mu _2 \mu _{2 p} \mu _{2 (p+1)}-48 \mu _{p+3} \mu _{3 p+1} \\ 5 & 120 \mu _{p+1}^4-240 \mu _2 \mu _{2 p} \mu _{p+1}^2+240 \mu _{2 (p+1)} \mu _{p+1}^2-240 \mu _{2 p} \mu _{p+3} \mu _{p+1}-240 \mu _2 \mu _{3 p+1} \mu _{p+1}+120 \mu _2^2 \mu _{2 p}^2+60 \mu _4 \mu _{2 p}^2+60 \mu _{2 (p+1)}^2+60 \mu _2^2 \mu _{4 p}+20 \mu _4 \mu _{4 p}+120 \mu _2 \mu _{2 p} \mu _{2 (p+1)}-80 \mu _{p+3} \mu _{3 p+1} \\ 6 & 360 \mu _{p+1}^4-720 \mu _2 \mu _{2 p} \mu _{p+1}^2+480 \mu _{2 (p+1)} \mu _{p+1}^2-480 \mu _{2 p} \mu _{p+3} \mu _{p+1}-480 \mu _2 \mu _{3 p+1} \mu _{p+1}+360 \mu _2^2 \mu _{2 p}^2+120 \mu _4 \mu _{2 p}^2+90 \mu _{2 (p+1)}^2+120 \mu _2^2 \mu _{4 p}+30 \mu _4 \mu _{4 p}+240 \mu _2 \mu _{2 p} \mu _{2 (p+1)}-120 \mu _{p+3} \mu _{3 p+1} \\ 7 & 840 \mu _{p+1}^4-1680 \mu _2 \mu _{2 p} \mu _{p+1}^2+840 \mu _{2 (p+1)} \mu _{p+1}^2-840 \mu _{2 p} \mu _{p+3} \mu _{p+1}-840 \mu _2 \mu _{3 p+1} \mu _{p+1}+840 \mu _2^2 \mu _{2 p}^2+210 \mu _4 \mu _{2 p}^2+126 \mu _{2 (p+1)}^2+210 \mu _2^2 \mu _{4 p}+42 \mu _4 \mu _{4 p}+420 \mu _2 \mu _{2 p} \mu _{2 (p+1)}-168 \mu _{p+3} \mu _{3 p+1} \\ 8 & 1680 \mu _{p+1}^4-3360 \mu _2 \mu _{2 p} \mu _{p+1}^2+1344 \mu _{2 (p+1)} \mu _{p+1}^2-1344 \mu _{2 p} \mu _{p+3} \mu _{p+1}-1344 \mu _2 \mu _{3 p+1} \mu _{p+1}+1680 \mu _2^2 \mu _{2 p}^2+336 \mu _4 \mu _{2 p}^2+168 \mu _{2 (p+1)}^2+336 \mu _2^2 \mu _{4 p}+56 \mu _4 \mu _{4 p}+672 \mu _2 \mu _{2 p} \mu _{2 (p+1)}-224 \mu _{p+3} \mu _{3 p+1} \\ 9 & 3024 \mu _{p+1}^4-6048 \mu _2 \mu _{2 p} \mu _{p+1}^2+2016 \mu _{2 (p+1)} \mu _{p+1}^2-2016 \mu _{2 p} \mu _{p+3} \mu _{p+1}-2016 \mu _2 \mu _{3 p+1} \mu _{p+1}+3024 \mu _2^2 \mu _{2 p}^2+504 \mu _4 \mu _{2 p}^2+216 \mu _{2 (p+1)}^2+504 \mu _2^2 \mu _{4 p}+72 \mu _4 \mu _{4 p}+1008 \mu _2 \mu _{2 p} \mu _{2 (p+1)}-288 \mu _{p+3} \mu _{3 p+1} \\ 10 & 5040 \mu _{p+1}^4-10080 \mu _2 \mu _{2 p} \mu _{p+1}^2+2880 \mu _{2 (p+1)} \mu _{p+1}^2-2880 \mu _{2 p} \mu _{p+3} \mu _{p+1}-2880 \mu _2 \mu _{3 p+1} \mu _{p+1}+5040 \mu _2^2 \mu _{2 p}^2+720 \mu _4 \mu _{2 p}^2+270 \mu _{2 (p+1)}^2+720 \mu _2^2 \mu _{4 p}+90 \mu _4 \mu _{4 p}+1440 \mu _2 \mu _{2 p} \mu _{2 (p+1)}-360 \mu _{p+3} \mu _{3 p+1} \\ \end{array}$$

We see the exact same products of moments for all values of $$m\geq 4$$.

The general form by inspection is

λ[m_] := Expand[FactorialPower[m, 4]
\!$$\*SubsuperscriptBox[\(μ$$, $$2$$, $$2$$]\)
\!$$\*SubsuperscriptBox[\(μ$$, $$2\ p$$, $$2$$]\) +
FactorialPower[m, 2] Subscript[μ, 4] Subscript[μ, 4 p] +
FactorialPower[m, 3] (Subscript[μ, 4]
\!$$\*SubsuperscriptBox[\(μ$$, $$2\ p$$, $$2$$]\) + \!$$\*SubsuperscriptBox[\(μ$$, $$2$$, $$2$$]\
\*SubscriptBox[$$μ$$, $$4\ p$$]\)) +
FactorialPower[m, 4] (-2 Subscript[μ, 2] Subscript[μ, 2 p]
\!$$\*SubsuperscriptBox[\(μ$$, $$1 + p$$, $$2$$]\) +
\!$$\*SubsuperscriptBox[\(μ$$, $$1 + p$$, $$4$$]\)) +
FactorialPower[m,
3] (2 Subscript[μ, 2] Subscript[μ, 2 p] Subscript[μ,
2 (1 + p)] + 4
\!$$\*SubsuperscriptBox[\(μ$$, $$1 + p$$, $$2$$]\) Subscript[μ, 2 (1 + p)]) +
FactorialPower[m, 2] 3
\!$$\*SubsuperscriptBox[\(μ$$, $$2\ \((1 + p)$$\), $$2$$]\) -
4 FactorialPower[m, 3] (
Subscript[μ, 2 p] Subscript[μ, 1 + p] Subscript[μ,
3 + p] +
Subscript[μ, 2] Subscript[μ, 1 + p] Subscript[μ,
1 + 3 p]) -
4 FactorialPower[m, 2] Subscript[μ, 3 + p] Subscript[μ,
1 + 3 p]]


$$\mu _2^2 m^{(4)} \mu _{2 p}^2+3 m^{(2)} \mu _{2 (p+1)}^2+\mu _4 m^{(2)} \mu _{4 p}+m^{(3)} \left(\mu _2^2 \mu _{4 p}+\mu _4 \mu _{2 p}^2\right)+m^{(4)} \left(\mu _{p+1}^4-2 \mu _2 \mu _{2 p} \mu _{p+1}^2\right)+m^{(3)} \left(4 \mu _{2 (p+1)} \mu _{p+1}^2+2 \mu _2 \mu _{2 p} \mu _{2 (p+1)}\right)-4 m^{(2)} \mu _{p+3} \mu _{3 p+1}-4 m^{(3)} \left(\mu _{2 p} \mu _{p+1} \mu _{p+3}+\mu _2 \mu _{p+1} \mu _{3 p+1}\right)$$

This general function also works for $$m>1$$. There are probably Mathematica functions to make this more automatic.

• Hmm that looked promising at first but I think you are not allowed to just substitute the moments. The expectation operator is linear and hence requires the terms to be a sum (which is the reason why I am even doing this). I think you should wrap the sum into an Expand. And even then it still does not solve the root issue: E(xi xj)=E(xi)E(xj) only if i!=j. So there must be a procedure that figures out the cross terms (i.e. when the different ms are equal and when not).
– divB
Commented Jan 27, 2019 at 10:58
• You are correct. When I looked at an individual sum for a specific value of $m$ I mistakenly thought that it was expanded. But it was not. I'll fix that shortly.
– JimB
Commented Jan 27, 2019 at 16:58
• Thank you for fixing - I get the idea and aplied to my scenario. Unfortunately I said I have many multi-sums - this was just a simple example. My actual issue arises from a Taylor expansion of a determinant - many sums! My current term is a 6-fold sum (over different counting variables) and THAT squared, i.e. 12-fold sum! If I take the largest counting variable ($m=5$) I am already at $5^{12}=244140625$ terms. They SHOULD drop to just 72 terms. But even for small numbers I had to cancel computation after hours. Worse than that, to FIND a pattern, I would need $m=10$ or so.
– divB
Commented Jan 28, 2019 at 5:34
• Hence this does not seem to be feasible...
– divB
Commented Jan 28, 2019 at 5:34
• You might want to look into symmetric functions: F. N. David, M. G. Kendall & D. E. Barton (1966) Symmetric Function and Allied Tables, Cambridge University Press. That's what I used in the second answer I gave. Also, the second answer is certainly feasible for the example you gave. If possible, giving an example that has the complexity of your real issue might generate more efficient approaches.
– JimB
Commented Jan 28, 2019 at 5:56

Here is an alternative approach with a bit less investigation and additional typing involved. (And I'm sure, too, that this can be automated much more efficiently.)

First define some helper functions to process the sums (i.e., take the expectations) in the form that you have:

$$\sum _{m_1=1}^m x_{m_1}^{e_1}$$ $$\sum _{m_1=1}^m \sum _{m_2=1}^m x_{m_1}^{e_1} x_{m_2}^{e_2}$$ $$\sum _{m_1=1}^m \sum _{m_2=1}^m \sum _{m_3=1}^m x_{m_1}^{e_1} x_{m_2}^{e_2} x_{m_3}^{e_3}$$ $$\sum _{m_1=1}^m \sum _{m_2=1}^m \sum _{m_3=1}^m \sum _{m_4=1}^m x_{m_1}^{e_1} x_{m_2}^{e_2} x_{m_3}^{e_3} x_{m_4}^{e_4}$$

sum1[e_, m_] := m μ[e]
sum2[e1_, e2_, m_] := m μ[e1 + e2] + m (m - 1) μ[e1] μ[e2]
sum3[e1_, e2_, e3_, m_] :=
m (m - 1) (m - 2) μ[e1] μ[e2] μ[e3] +
m (m - 1) μ[e1] μ[e2 + e3] +
m (m - 1) μ[e2] μ[e1 + e3] +
m (m - 1) μ[e3] μ[e1 + e2] + m μ[e1 + e2 + e3]
sum4[e1_, e2_, e3_, e4_, m_] :=
m (m - 1) (m - 2) (m - 3) μ[e1] μ[e2] μ[e3] μ[e4] +
m (m - 1) (m - 2) μ[e1 + e2] μ[e3] μ[e4] +
m (m - 1) (m - 2) μ[e2] μ[e1 + e3] μ[e4] +
m (m - 1) (m - 2) μ[e1] μ[e2 + e3] μ[e4] +
m (m - 1) μ[e1 + e2 + e3] μ[e4] +
m (m - 1) (m - 2) μ[e2] μ[e3] μ[e1 + e4] +
m (m - 1) μ[e2 + e3] μ[e1 + e4] +
m (m - 1) (m - 2) μ[e1] μ[e3] μ[e2 + e4] +
m (m - 1) μ[e1 + e3] μ[e2 + e4] +
m (m - 1) μ[e3] μ[e1 + e2 + e4] +
m (m - 1) (m - 2) μ[e1] μ[e2] μ[e3 + e4] +
m (m - 1) μ[e1 + e2] μ[e3 + e4] +
m (m - 1) μ[e2] μ[e1 + e3 + e4] +
m (m - 1) μ[e1] μ[e2 + e3 + e4] + m μ[e1 + e2 + e3 + e4]


For the example you gave the expanded value is

Expand[(x[m1]^2 x[m2]^(2 p) - x[m1]^(p + 1) x[m2]^(p + 1)) (x[m3]^2 x[m4]^(2 p) -
x[m3]^(p + 1) x[m4]^(p + 1))]
(* x[m1]^2       x[m2]^(2 p)   x[m3]^2       x[m4]^(2 p) -
x[m1]^(1 + p) x[m2]^(1 + p) x[m3]^2       x[m4]^(2 p) -
x[m1]^2       x[m2]^(2 p)   x[m3]^(1 + p) x[m4]^(1 + p) +
x[m1]^(1 + p) x[m2]^(1 + p) x[m3]^(1 + p) x[m4]^(1 + p) *)


So using the helper functions the desired expected value for your expression is

sum4[2, 2 p, 2, 2 p, m] - sum4[1 + p, 1 + p, 2, 2 p, m] -
sum4[2, 2 p, 1 + p, 1 + p, m] + sum4[1 + p, 1 + p, 1 + p, 1 + p, m]

(* (-3 + m) (-2 + m) (-1 + m) m μ[2]^2 μ[
2 p]^2 + (-2 + m) (-1 + m) m μ[4] μ[
2 p]^2 + (-2 + m) (-1 + m) m μ[2]^2 μ[
4 p] + (-1 + m) m μ[4] μ[4 p] -
2 (-3 + m) (-2 + m) (-1 + m) m μ[2] μ[2 p] μ[
1 + p]^2 + (-3 + m) (-2 + m) (-1 + m) m μ[1 + p]^4 -
4 (-2 + m) (-1 + m) m μ[2 p] μ[1 + p] μ[3 + p] +
2 (-2 + m) (-1 + m) m μ[2] μ[2 p] μ[2 + 2 p] +
4 (-2 + m) (-1 + m) m μ[1 + p]^2 μ[2 + 2 p] +
3 (-1 + m) m μ[2 + 2 p]^2 -
4 (-2 + m) (-1 + m) m μ[2] μ[1 + p] μ[1 + 3 p] -
4 (-1 + m) m μ[3 + p] μ[1 + 3 p] *)


When $$m=8$$

sum4[2, 2 p, 2, 2 p, m] - sum4[1 + p, 1 + p, 2, 2 p, m] -
sum4[2, 2 p, 1 + p, 1 + p, m] + sum4[1 + p, 1 + p, 1 + p, 1 + p, m] /. m -> 8

(* 1680 μ[2]^2 μ[2 p]^2 + 336 μ[4] μ[2 p]^2 + 336 μ[2]^2 μ[4 p] + 56 μ[4] μ[4 p] -
3360 μ[2] μ[2 p] μ[1 + p]^2 + 1680 μ[1 + p]^4 - 1344 μ[2 p] μ[1 + p] μ[3 + p] +
672 μ[2] μ[2 p] μ[2 + 2 p] + 1344 μ[1 + p]^2 μ[2 + 2 p] + 168 μ[2 + 2 p]^2 -
1344 μ[2] μ[1 + p] μ[1 + 3 p] - 224 μ[3 + p] μ[1 + 3 p] *)


For the more complicated function $$D$$ obtaining a compact formula for the expectation of $$D^2$$ for general $$M$$ and $$R$$ is a lot more difficult. For the case of $$M>0$$ and $$R=1$$ the expected value of $$D^2$$ is given by the following function:

fD2R1[M_] :=
FactorialPower[M, 4] μ[2 p1]^2 μ[2 p2]^2 +
FactorialPower[M, 3] μ[4 p1] μ[2 p2]^2 +
FactorialPower[M, 3] μ[2 p1]^2 μ[4 p2] +
FactorialPower[M, 2] μ[4 p1] μ[4 p2] +
FactorialPower[M, 3] 4 μ[p1 + p2]^2 μ[2 p1 + 2 p2] -
FactorialPower[M, 4] 2 μ[2 p1] μ[2 p2] μ[p1 + p2]^2 +
FactorialPower[M, 4] μ[p1 + p2]^4 -
FactorialPower[M, 3] 4 μ[2 p2] μ[p1 + p2] μ[3 p1 + p2] +
FactorialPower[M, 3] 2 μ[2 p1] μ[2 p2] μ[2 p1 + 2 p2] +
FactorialPower[M, 2] 3 μ[2 p1 + 2 p2]^2 -
FactorialPower[M, 3] 4 μ[2 p1] μ[p1 + p2] μ[p1 + 3 p2] -
FactorialPower[M, 2] 4 μ[3 p1 + p2] μ[p1 + 3 p2]


For $$M=50$$ this results in

fD2R1[50]
(* 5527200 μ[2 p1]^2 μ[2 p2]^2 + 117600 μ[4 p1] μ[2 p2]^2 + 117600 μ[2 p1]^2 μ[4 p2] +
2450 μ[4 p1] μ[4 p2] - 11054400 μ[2 p1] μ[2 p2] μ[p1 + p2]^2 +
5527200 μ[p1 + p2]^4 - 470400 μ[2 p2] μ[p1 + p2] μ[3 p1 + p2] +
235200 μ[2 p1] μ[2 p2] μ[2 p1 + 2 p2] + 470400 μ[p1 + p2]^2 μ[2 p1 + 2 p2] +
7350 μ[2 p1 + 2 p2]^2 - 470400 μ[2 p1] μ[p1 + p2] μ[p1 + 3 p2] -
9800 μ[3 p1 + p2] μ[p1 + 3 p2] *)


Still working on the general case...

Here is a more formal approach that I think will be feasible even for your more complex sums. (Essentially it's a formalizing of the "alternative approach" answer I gave earlier.)

One turns the sums into power sums, then into a sum of augmented symmetric sums for which the expectations are readily evident. A general reference for constructing augmented symmetric sums is Merca, 2015 and tables of the associated coefficients can be found at David, et. al., 1967. Alternatively, the software add-on for Mathematica called mathStatica has built-in functions to make this process much, much easier.

Here is a simple example of the process (which will be followed sometime later with your original request).

Consider the need to find the expectation of a sum

$$\sum _{m_1=1}^m \sum _{m_2=1}^m x_{m_1}^2 x_{m_2}^2$$

which we can rewrite:

$$\sum _{m_1=1}^m \sum _{m_2=1}^m x_{m_1}^2 x_{m_2}^2=\left(\sum _{m_1=1}^m x_{m_1}^2\right) \sum _{m_2=1}^m x_{m_2}^2=\sum _{m_1=1}^m x_{m_1}^4+\sum _{m_1\neq m_2} x_{m_1}^2 x_{m_2}^2$$

With the last equation we see that because there are no more overlaps in indices, we can immediately write down the expected values:

$$m \mu^\prime_4 + m(m-1)\mu_2^{\prime 2}$$

where $$\mu^\prime_k$$ is the $$k$$-th raw moment.

To use the power sum and augmented symmetric function notation we have

$$s_k=\sum _{i=1}^m x_i^k$$

$$a_{\left\{k_1,k_2\right\}}=\sum _{m_1\neq m_2} x_{m_1}^{k_1} x_{m_2}^{k_2}$$

For the simple example we have

$$\sum _{m_1=1}^m \sum _{m_2=1}^m x_{m_1}^2 x_{m_2}^2=s_2^2=a_{\{4\}}+a_{\{2,2\}}$$

with

$$E(a_{\{4\}})=m \mu_4^\prime$$ $$E(a_{\{2,2\}})=m(m-1) \mu_2^{\prime 2}$$

For your original request here are the commands to find the expectation using functions from mathStatica:

$$\left(\sum _{m_1=1}^m \sum _{m_2=1}^m \left(x_{m_1}^2 x_{m_2}^{2 p}-x_{m_1}^{p+1} x_{m_2}^{p+1}\right)\right){}^2=\left(\left(\sum _{m_1=1}^m x_{m_1}^2\right) \sum _{m_2=1}^m x_{m_2}^{2 p}-\left(\sum _{m_1=1}^m x_{m_1}^{p+1}\right) \sum _{m_2=1}^m x_{m_2}^{p+1}\right){}^2$$ $$=\left(s_2 s_{2 p}-s_{p+1}^2\right){}^2=s_{p+1}^4-2 s_2 s_{2 p} s_{p+1}^2+s_2^2 s_{2 p}^2$$

sum = PowerSumToAug[{2, 2, 2 p, 2 p}, s, a][[2]] -
2 PowerSumToAug[{2, 2 p, 1 + p, 1 + p}, s, a][[2]] +
PowerSumToAug[{(1 + p)^4}, s, a][[2]]


$$2 a_{\{2,4 p+2\}}+a_{\{4,4 p\}}+2 a_{\{2 p,2 p+4\}}+4 a_{\{p+1,3 p+3\}}+5 a_{\{2 p+2,2 p+2\}}+a_{\{2,2,4 p\}}+4 a_{\{2,2 p,2 p+2\}}+a_{\{4,2 p,2 p\}}+6 a_{\{p+1,p+1,2 p+2\}}+a_{\{2,2,2 p,2 p\}}-2 \left(a_{\{2,4 p+2\}}+a_{\{2 p,2 p+4\}}+2 a_{\{p+1,3 p+3\}}+2 a_{\{p+3,3 p+1\}}+a_{\{2 p+2,2 p+2\}}+a_{\{2,2 p,2 p+2\}}+2 a_{\{2,p+1,3 p+1\}}+2 a_{\{2 p,p+1,p+3\}}+a_{\{p+1,p+1,2 p+2\}}+a_{\{2,2 p,p+1,p+1\}}+a_{\{4 p+4\}}\right)+a_{\{p+1,p+1,p+1,p+1\}}+2 a_{\{4 p+4\}}$$

expectAug[t_] := (Thread[μ[t]] /. List -> Times) Product[m - i, {i, 0, Length[t] - 1}]
expectedSum = sum /. Subscript[a, t_] :> expectAug[t]

(* (-3 + m) (-2 + m) (-1 + m) m μ[2]^2 μ[2 p]^2 + (-2 + m) (-1 + m) m μ[4] μ[2 p]^2 +
(-2 + m) (-1 + m) m μ[2]^2 μ[4 p] + (-1 + m) m μ[4] μ[4 p] +
(-3 + m) (-2 + m) (-1 + m) m μ[1 + p]^4 +
4 (-2 + m) (-1 + m) m μ[2] μ[2 p] μ[2 + 2 p] +
6 (-2 + m) (-1 + m) m μ[1 + p]^2 μ[2 + 2 p] +
5 (-1 + m) m μ[2 + 2 p]^2 + 2 (-1 + m) m μ[2 p] μ[4 + 2 p] +
4 (-1 + m) m μ[1 + p] μ[3 + 3 p] + 2 (-1 + m) m μ[2] μ[2 + 4 p] + 2 m μ[4 + 4 p] -
2 ((-3 + m) (-2 + m) (-1 + m) m μ[2] μ[2 p] μ[1 + p]^2 +
2 (-2 + m) (-1 + m) m μ[2 p] μ[1 + p] μ[3 + p] +
(-2 + m) (-1 + m) m μ[2] μ[2 p] μ[2 + 2 p] + (-2 + m) (-1 + m) m μ[1 + p]^2 μ[2 + 2 p] +
(-1 + m) m μ[2 + 2 p]^2 + (-1 + m) m μ[2 p] μ[4 + 2 p] +
2 (-2 + m) (-1 + m) m μ[2] μ[1 + p] μ[1 + 3 p] +
2 (-1 + m) m μ[3 + p] μ[1 + 3 p] +
2 (-1 + m) m μ[1 + p] μ[3 + 3 p] +
(-1 + m) m μ[2] μ[2 + 4 p] + m μ[4 + 4 p]) *)


I think this approach can certainly handle the 12-deep sums you mention.

Update

The expectation for $$D$$ and $$D^2$$ (and even \$D^3 and higher powers) can be constructed using the functions from mathStatica. Here's how to do that:

d = -(\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$m = 1$$, $$M$$]$$\(( \*UnderoverscriptBox[\(\[Sum]$$, $$r = 1$$, $$R$$]
\*SuperscriptBox[$$x[r + \((\(-1$$ + m)\)\ R]\), $$p1$$])\)\ $$\*UnderoverscriptBox[\(\[Sum]$$, $$r = 1$$, $$R$$]
\*SuperscriptBox[$$x[ r + \((\(-1$$ + m)\)\ R]\), $$p2$$]\)\)\))^2 + (\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$m = 1$$, $$M$$]
\*SuperscriptBox[$$( \*UnderoverscriptBox[\(\[Sum]$$, $$r = 1$$, $$R$$]
\*SuperscriptBox[$$x[ r + \((\(-1$$ + m)\)\ R]\), $$p1$$])\), $$2$$]\)) \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$m = 1$$, $$M$$]
\*SuperscriptBox[$$( \*UnderoverscriptBox[\(\[Sum]$$, $$r = 1$$, $$R$$]
\*SuperscriptBox[$$x[r + \((\(-1$$ + m)\)\ R]\), $$p2$$])\), $$2$$]\)
(* Rewrite the function so that *)
d = d /. \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$r = 1$$, $$R$$]
\*SuperscriptBox[$$x[r + \((\(-1$$ + m)\)\ R]\), $$p_$$]\) ->
s[-1 + m, p]


This way, for example s[0,p1] has no common random variables with s[1,p1] and because of that statistical independence, the expectations will be simpler to construct.

Now we can expand d for a specific $$M$$ and general $$R$$ into power sums which will be converted to augmented symmetric sums and then finally to the expectations. We do that by using two functions (using mathStatica and copied nearly identical from the associated textbook) and a replacement rule.

expectAug[t_, r_] := (Thread[Subscript[
\!$$\*OverscriptBox[\(\[Mu]$$, $$\[Prime]$$]\), t]] /.
List -> Times) Product[r - i, {i, 0, Length[t] - 1}]

w[p_, r_] := PowerSumToAug[p, , a][[2]] /. Subscript[a, t_] :> expectAug[t, r]

rule = {s[i_, p1]^k1_ s[i_, p2]^k2_ :> w[{p1^k1, p2^k2}, R],
s[i_, p1]^k1_ s[i_, p2] :> w[{p1^k1, p2}, R],
s[i_, p1] s[i_, p2]^k2_ :> w[{p1, p2^k2}, R],
s[i_, p1] s[i_, p2] :> w[{p1, p2}, R],
s[i_, p_]^k_ :> w[{p^k}, R],
s[i_, p_] :> w[{p}, R]};


For the expectation of $$D$$ we have

(ed1 = Table[{M0, Expand[d /. M -> M0] //. rule /. M -> M0}, {M0, 1, 10}]) // TableForm


Looking at the expectation for different values of M, we a pattern and the general formula is

expectationD[M_, R_] := M (M - 1) (((-1 + R) R
\!$$\*SubsuperscriptBox[ OverscriptBox[\(\[Mu]$$, $$\[Prime]$$], $$p1$$, $$2$$]\) + R Subscript[
\!$$\*OverscriptBox[\(\[Mu]$$, $$\[Prime]$$]\), 2 p1]) ((-1 + R) R
\!$$\*SubsuperscriptBox[ OverscriptBox[\(\[Mu]$$, $$\[Prime]$$], $$p2$$, $$2$$]\) + R Subscript[
\!$$\*OverscriptBox[\(\[Mu]$$, $$\[Prime]$$]\),
2 p2]) - ((-1 + R) R Subscript[
\!$$\*OverscriptBox[\(\[Mu]$$, $$\[Prime]$$]\), p1] Subscript[
\!$$\*OverscriptBox[\(\[Mu]$$, $$\[Prime]$$]\), p2] + R Subscript[
\!$$\*OverscriptBox[\(\[Mu]$$, $$\[Prime]$$]\), p1 + p2])^2)


For the expectation of $$D^2$$ we have

(ed2 = Table[{M0, Expand[d^2 /. M -> M0] //. rule /. M -> M0}, {M0, 1, 10}]) // TableForm


This can be done for any other moment however, the computation time increases dramatically.

(ed3 = Table[{M0, Expand[d^3 /. M -> M0] //. rule /. M -> M0}, {M0, 1, 10}]) // TableForm