Effective nested summation

My problem is the following: given a list of probabilities $$p_0, p_1, \dots, p_M$$ and an integer $$n$$, I need to compute a function of a real variable $$x$$

$$W(x) = \sum_{m_1 = 0}^M \sum_{m_2 = 0}^M \dots \sum_{m_n = 0}^M \left( \prod_{j=1}^n \frac{ p_{m_j}}{ m_j !}\right) \frac{ \left(\sum_{i=1}^n m_i\right)! }{ ( - n)^{\sum_{i=1}^n m_i} } L_{\sum_{i=1}^n m_i } (x),$$ where $$L_i (x)$$ is the Laguerre polynomial LaguerreL[i,x].

My current (inefficient) solution is this:

w[x_, n_, probs_] :=   (* probs = {p0, p1, p2, .., pM} *)
Module[ {mM = Length@probs -1 , mj },
With[
{ sumIt = Sequence@@Table[ { mj[i], 0 , mM}, {i, 1, n}]},
Sum[
Module[ { sumMJ = Sum[ mj[j], {j,1, n}]},
Product[ probs[[ mj[j]+1]]/ ((mj[j])!), { j , 1, n}]
sumMJ! ( -1)^sumMJ / n^sumMJ LaguerreL[ sumMJ , x ]],
sumIt]
]]


It takes quite a long time to compute, e.g. with $$n=9$$ and $$M = 3$$ it returns the result after about 2 minutes. I guess this is because the approach is quite procedural. I would like to derive this function for the values $$n \leq 40$$ and $$M \leq 5$$, which I am afraid would take a few weeks.

Therefore, my question is, is it possible to compute this function more efficiently in Mathematica?

1 Answer

2nd Update: I've added the rationale I used to obtain the more compact sum.

Update: This is now (hopefully) an answer. At end of the case where $$M=1$$ is the general case.

Original "answer":

This is an extended comment. By setting $$M$$ to some value and then looking for a pattern for various values of $$n$$, one might find a quicker formula.

For example, consider $$M=1$$. Looking at the coefficients for a few values of $$n$$ and using FindSequenceFunction, I get the following:

c[n_] := DifferenceRoot[
Function[{y, i}, {(-(n + 1) + i) y[i] + n y[1 + i] == 0,  y[1] == 1}]][#] & /@ Range[n + 1]

w2[x_, n_] := Table[p[0]^(n - i) (-1)^i p[1]^i LaguerreL[i, x], {i, 0, n}] . c[n]

AbsoluteTiming[w2[x, 15];]
(* {0.013407, Null} *)
AbsoluteTiming[w[x, 15, {p[0], p[1]}];]
(* {8.23932, Null} *)

AbsoluteTiming[w2[x, 40];]
(* {0.0418239, Null} *)
AbsoluteTiming[w[x, 40, {p[1], p[2]}];]
(* Still waiting after several minutes *)


General case:

w3[x_, n_, m_] := Module[{t, r},
t = Flatten[Permutations[#] & /@ (PadRight[#, m + 1] & /@IntegerPartitions[n, m + 1]), 1];
r = Range[0, m];
c1 = (n!/Times @@ (#!)) (1/Times @@ (r!^#)) (r . #)!/(-n)^(r . #) & /@ t;
c2 = Product[p[i]^#[[i + 1]], {i, 0, m}] & /@ t;
c3 = LaguerreL[r . #, x] & /@ t;
Total[c1 c2 c3]
]

AbsoluteTiming[w[x, 3, {p[0], p[1], p[2]}]]


AbsoluteTiming[w3[x, 3, 2]]


Check on if the results are the same:

w[x, 3, {p[0], p[1], p[2]}] - w3[x, 3, 2]
(* 0 *)


A more lengthy calculation:

AbsoluteTiming[w3[x, 40, 4];]
(* {251.448, Null} *)


Rationale for the more compact summation:

We have the following sum:

$$W(x) = \sum_{m_1 = 0}^M \sum_{m_2 = 0}^M \dots \sum_{m_n = 0}^M \left( \prod_{j=1}^n \frac{ p_{m_j}}{ m_j !}\right) \frac{ \left(\sum_{i=1}^n m_i\right)! }{ ( - n)^{\sum_{i=1}^n m_i} } L_{\sum_{i=1}^n m_i } (x)$$

Because we can write

$$\prod_{j=1}^n \frac{ p_{m_j}}{ m_j !}=\prod_{i=0}^M p_i^{k_i}/\prod_{j=1}^n m_j !$$

where $$\sum_{i=0}^M k_i=n$$, for each combination of the $$m_i$$ indices the $$n(M+1)$$ terms are of the form

$$\left(\prod_{i=0}^M p_i^{k_i}/\prod_{j=1}^n m_j ! \right) \frac{ \left(\sum_{i=1}^n m_i\right)! }{ ( - n)^{\sum_{i=1}^n m_i} } L_{\sum_{i=1}^n m_i } (x)$$

For example, suppose $$M=2$$, $$n=5$$, $$m_1=0$$, $$m_2=1$$, $$m_3=1$$, $$m_4=2$$, and $$m_5=1$$. The associated term is

$$\left(\prod_{i=0}^M p_i^{k_i}/\prod_{j=1}^n m_j ! \right) \frac{ \left(\sum_{i=1}^n m_i\right)! }{ ( - n)^{\sum_{i=1}^n m_i} } L_{\sum_{i=1}^n m_i } (x)=p_0 p_1^3 p_2 \frac{5!}{0! (1!)^3 2! (-5)^5}L_5 (x)$$

which simplifies to

$$-\frac{12}{625}p_0 p_1^3 p_2 L_5 (x)$$

There are 19 other combinations of the $$m_i$$ values that result in the same exact term. How to calculate that multiplier? Of the 5 values of $$m_i$$ there is one 0, three 1's, and one 2. So 1, 3, and 1 are the powers of $$p_0$$, $$p_1$$, and $$p_2$$, respectively. The desired count is given by

$$\frac{5}{1!3!1!}=20$$

The coefficient associated with $$p_0 p_1^3 p_2 L_5 (x)$$ is then

$$20*(-\frac{12}{625})=-\frac{48}{125}$$

So (again, for this example) we just need to generate all of the combinations of $$p_0^{k_0} p_1^{k_1} p_2^{k_2}$$ where $$k_i \geq 0$$ and $$k_0+k_1+k_2=n$$, construct the associated term, and then multiply by the number of combinations of the $$m_i$$ that will generate that term (which is $$\frac{n!}{k_0!k_1!k_2!}$$).

• Thank you very much, this is very impressive. I would absolutely love to learn more about the way you were able to come up with this solution. Dec 7, 2023 at 8:48