# How to simplify a complicated Sum in terms of power Sums?

For example, I have:

$a=\sum _{r=1}^n x_r \left(\left(\sum _{i=1}^n x_i-x_r\right){}^2-\sum _{i=1}^n x_i^2\right)$

a =
Sum[Subscript[x, r]*
((Sum[Subscript[x, i], {i, 1, n}] - Subscript[x, r])^2 -
Sum[Subscript[x, i]^2,  {i, 1, n}]),
{r, 1, n}]


and $\sum _{i=1}^n x_i=s_1, \sum _{i=1}^n x_i^2=s_2, \sum _{i=1}^n x_i^3=s_3$

 c2 = Sum[Subscript[x, i]^2, {i, 1, n}] == s2


I would like to represent $a$ by $s_1, s_2, s_3$, e.g. $a=s_1^3+s_3-3s_1s_2$ . How should I do it? I tried Solve or Eliminate, but couldn't find a way.

-
Could you provide the Mathematica code of your definitions and by a bit clearer about what kind of output you expect? –  Sjoerd C. de Vries Mar 10 '13 at 23:29
Related question on SO –  Michael E2 Mar 11 '13 at 12:57
@colinfang I changed the title. The question didn't seem to me to have to do with subscripts, but with manipulating sums. If you don't like the new title, you can change it back, or change it to something better. –  Michael E2 Mar 11 '13 at 13:30

As bound variables, $r$ and $i$ must play no role in the expansion and so they shouldn't even appear in our solution. It should be equally evident that $n$ is just along for the ride as a placeholder for the upper limit; we could call it anything, and therefore we may call it nothing and ignore it except in expressions where it will necessarily appear in the output. Notice, too, that $x$ disappears under the summation: it really just serves as a placeholder for arguments of other expressions being summed.

Consequently, expressions like "$\sum_{i=1}^n f(x_i)$" should be understood as operations on $f$ itself. What properties should they have? Only the obvious ones implied by the ring operations $+$ and $\times$ in the algebra of polynomials in the $x_i$. This gives just four simple rules, after we establish the unique roles of the symbols "$x$", "$s$", and "$n$":

ClearAll[sum, x, s, n];

1. Converting sums of powers of $x$ into the power sum variables $s_i$:

sum[x] := Subscript[s, 1]; sum[x^p_] := Subscript[s, p];


sum[a_ + b_] := sum[a] + sum[b];

3. Linearity (with respect to "scalars" which do not have any $x_i$ in them):

sum[Times[a_, y__]] /; FreeQ[a, x] := a sum[Times[y]];

4. The effect of summing constant values (this is the only place $n$ need appear):

sum[a_] /; FreeQ[a, x] := n a;


That should do it. But to perform the algebra, we need to expand algebraic combinations of everything possible into powers of $x$ (and otherwise leave everything else alone):

expand[a_] /; ! FreeQ[a, sum] := Map[Expand[#, x] &, a, Infinity];
expand[a_] := a;


### Examples

The question:

sum[x ((sum[x] - x)^2 - sum[x^2])] // expand


$s_1^3-3 s_1 s_2+s_3$

Variance:

sum[(x - sum[x]/n)^2] // expand


$-\frac{s_1^2}{n}+s_2$

Skewness:

sum[(x - sum[x]/n)^3] / sum[(x - sum[x]/n)^2] ^(3/2) // expand


$\frac{\frac{2 s_1^3}{n^2}-\frac{3 s_1 s_2}{n}+s_3}{\left(-\frac{s_1^2}{n}+s_2\right){}^{3/2}}$

etc.

We might worry about misapplications, such as to non-polynomial functions. Not to fear:

sum[Exp[(x + 1)^2]] // expand


$\text{sum}\left[e^{1+2 x+x^2}\right]$

The expansion proceeds insofar as it can, but sum does not know how to go any further, and so stops. sum is also ignorant of other non-polynomial objects, but when they can be converted to polynomials, it works:

sum[Series[Exp[x], {x, 0, 4}]  // Normal] // expand


$n+s_1+\frac{s_2}{2}+\frac{s_3}{6}+\frac{s_4}{24}$

(This is a series expansion through order four of $\sum_{i=1}^n e^{x_i}$.)

### Edit

Michael E2, in an extension to his answer, reminds us of the value of simplifying multiple sums. His generalization further clarifies the nature of these operations. Because the variable name "$x$" does not explicitly appear in the $s_i$ notation, the variable name is immaterial. What matters are the index names, only insofar as they are used to connect an $x$ with its enclosing summation: as I remarked at the outset, as bound variables they must disappear at the end.

Thus, for example, we could write sum[i^3, i] for $\sum_{i=1}^n x_i^3$: the first argument to sum is a polynomial (or rational function, even) and the second one is a symbol to indicate what is varying over the summation.

A comparable extension is trivial to make: just include the second argument explicitly in the definition of sum. So that we can reproduce the previous solution (where x was the only variable), we can have the second argument default to x if it's missing. Here is the entire (generalized) solution:

ClearAll[sum, s, n];
sum[a_] := sum[a, x];
sum[x_, x_] := Subscript[s, 1]; sum[x_^p_, x_] := Subscript[s, p];
sum[Times[a_, y__], x_] /; FreeQ[a, x] := a sum[Times[y], x];
sum[a_ + b_, x_] := sum[a, x] + sum[b, x];
sum[a_, x_] /; FreeQ[a, x] := n a;
expand[a_] /; ! FreeQ[a, sum] := Map[Expand, a, Infinity]; expand[a_] := a;


For example, compute $\sum_{k=1}^n\sum_{j=1}^n\sum_{i=1}^n (x_i+x_j+x_k - 3\bar{x})^2$, where $\bar{x} = \sum_{i=1}^n x_i / n$ is the average:

sum[sum[sum[(i + j + k - 3 sum[i, i]/n)^3, i], j], k]


$6 s_1^3-9 n s_1 s_2+3 n^2 s_3$

(Notice that it was no problem for i to appear in two summations.) You do have to be a little careful not to go overboard and stretch the use of sum too far:

sum[x^2, x^2] // expand


$s_1$

This surprising result is correct: the symbol x^2 in the second argument is being used to represent the variable summed over; in terms of this pattern, the first argument x^2 is the first power, whence the result is indeed $s_1$.

One can rather easily go further with generalizations, depending on the need. The next big step would be to implement multivariate sums, such as $\sum_{i=1}^n\sum_{j=1}^n x_i y_j^2$, perhaps to be written $s_1(x)s_2(y)$ or--assuming a specific ordering $(x,y)$ of the variables--even as $s_{1,2}$. As we appear to have gone beyond the scope of the present question, I leave these generalizations to interested readers. I hope that the examples given here have shown how stripping the notation down to its essentials can lead to simple and clear solutions.

-
+1 Here I was trying to formulate it as a general problem, but wasn't able to articulate it coherently –  belisarius Mar 12 '13 at 18:20
+1 Nice and clear. Glad if any hints were helpful -- oops, you meant belisarius. That's ok :) –  Michael E2 Mar 12 '13 at 18:22
One thing I am stuck on: in principle it is possible to format these sum expressions into exactly the original notation used in the question. For instance, Format[sum[x^p_]] := \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$i = 1$$, $$n$$] \*SuperscriptBox[ SubscriptBox[$$x$$, $$i$$], $$p$$]\) handles simple cases like sum[x^3]. But how to extend this to the general case? (Among other things, we need a way to introduce different bound variable names such as "i" and "r", or perhaps $i_1, i_2,$ etc.) –  whuber Mar 12 '13 at 18:28

I guess you have to teach Mathematica some rules of algebra for sums:

a=Sum[x[r] ((Sum[x[i], {i, n}] - x[r])^2 - Sum[x[i]^2, {i, n}]), {r, n}];

rules ={Sum[x[i_], {i_, n}] -> s[1], Sum[x[i_]^(k_), {i_, n}] :> s[k], (* OP's terms *)
(* some rules of algebra: *)
Sum[(a_Plus)*(b_), {i_, n}] :> Sum[Expand[a*b], {i, n}],
Sum[(a_)*(b_), {i_, n}] :>
(Sum[a*b, {i, n}] /. Sum[(c_)?(FreeQ[#1, i] & )*(d_), {i, n}] :> c*Sum[d, {i, n}]),
Sum[(a_) + (b_), {i_, n}] :> Sum[a, {i, n}] + Sum[b, {i, n}]};

a //. rules
(* -> s[1]^3 - 3 s[1] s[2] + s[3] *)


Using Simplify

[Edit: I added some more rules to handle some multiple sums and powers]

One can add some algebraic transformation rules and a special complexity function to Simplify and Mathematica will expand the sums as far as possible. The complexity function (normally just LeafCount) is altered to strong favor more sums and updated to favor evaluating constant sums to n times the summand. The coefficients -10 and -20 need to be large enough but seem to be otherwise arbitrary based on experimentation; however, I suspect that with more summation limits (in double sums etc.), the more one has to compensate for the greater number of leaves in the Sum.

sRules = {Sum[x[i_], {i_, n}] -> s[1], Sum[x[i_]^(k_), {i_, n}] :> s[k]};
(* factorRule factors out constant factors *)
factorRule = HoldPattern[Sum][summand_Times, lims__] :>
Select[summand, ! InternalDependsOnQ[#, First /@ {lims}] &]*
Sum[Select[summand, InternalDependsOnQ[#, First /@ {lims}] &], lims];
(* expandSum expands the summand and distributes the Sum *)
expandRule = HoldPattern[Sum][summand : _Times | _Power, lims__] :>
Sum[Expand[summand], lims];
expandSum[e_] := Distribute[e /. expandRule];  (* break sums apart *)
simplifySum[Sum[summand_, lims__]] :=  (* for multiple sums *)
Fold[Sum, mySimplify[summand], {lims}];

complexity[e_] := -10 Count[e, _Sum, {0, Infinity}] +  (* favors more Sums *)
-20 Count[e, n, {0, Infinity}] +   (* favors more sums and evaluated constant sums *)
Count[e, Sum[_Times, __], {0, Infinity}] +   (* favors factoring out constants *)
LeafCount[e];   (* favors smaller expression *)

Simplify[a,
TransformationFunctions -> {Automatic, # /. factorRule &, expandSum},
ComplexityFunction -> complexity] /. sRules


Compared to first: plus -- uses Simplify; minus -- undocumented internal function, explicit algebra rules.

More examples

Summand is a function of two indices:

b = Sum[x[r] Sum[(x[i] - x[r])^2, {i, n}], {r, n}];
mySimplify[b] /. sRules
(* -> -2 s[1]^2 + 2 n s[2] *)


Double sum:

c = Sum[(x[i] - x[r])^2, {r, n}, {i, n}];
mySimplify[c] /. sRules
(* -> -2 s[1]^2 + 2 n s[2] *)


As whuber alluded to in a comment to his answer, you can format the s[k] with subscripts as follows:

Format[s[k_]] := Subscript[s, k]


By the way, the rule for factorRule was found on an archive of the newsgroup sci.math.symbolic (19 Nov 2007). I don't know about Internal functions except what I've learned on this site.

-
+1 But I really don't like none of our solutions. There must be an easier way –  belisarius Mar 11 '13 at 3:06
It would be nice if someone found one. None of the usual functions seem to simplify indefinite sums. –  Michael E2 Mar 11 '13 at 3:11
@belisarius Of course, you can get the answer like this: Eliminate[{y == a, s[3] == Sum[x[i]^3, {i, n}], s[2] == Sum[x[i]^2, {i, n}], s[1] == Sum[x[i], {i, n}]} /. n -> 3, Table[x[i], {i, 3}]], for any n >= 3. –  Michael E2 Mar 11 '13 at 3:23
Yep, but the whole thing is about keeping n symbolic ... –  belisarius Mar 11 '13 at 3:25
@belisarius You give good hints :-). I have proposed a solution based on them. –  whuber Mar 12 '13 at 17:40
a = Sum[x[r] ((Sum[x[i], {i, n}] - x[r])^2 - Sum[x[i]^2, {i, n}]), {r, n}];

Module[{t0, t1},
rule = Sum[p1___  x[a_]^i__  p2___, {a_, n}] :> p1 p2 s[i];
a /. x[k_] -> t0 x[k]^t1
/. Sum[a__, l_List] :> (Sum[#, l] & /@ Expand@a)
//. rule /.  t0 -> 1 /. t1 -> 1]

(*s[1]^3 - 3 s[1] s[2] + s[3]*)

-
More succinct than mine. And now you have s[1] instead of s1. I guess I'll post mine. But there's one part I'm not proud of...at least. –  Michael E2 Mar 11 '13 at 2:56
@MichaelE2 It took me six iterations to get here, and I think it isn't robust enough –  belisarius Mar 11 '13 at 2:59

If you are willing to work with specific integer values of n then you can use SymmetricReduction. Here is an example with n=4.

poly[n_] :=
Sum[Subscript[x,
r]*((Sum[Subscript[x, i], {i, 1, n}] - Subscript[x, r])^2 -
Sum[Subscript[x, i]^2, {i, 1, n}]), {r, 1, n}]

First[
SymmetricReduction[poly[4], Table[Subscript[x, i], {i, 4}],
Table[s[j], {j, 4}]]]

(* Out[146]= -s[1]^3 + 3 s[1] s[2] + 3 s[3] *)

-
However, it resolves into elementary symmetric polynomials, not the ones I expect. –  colinfang Mar 11 '13 at 15:18

Daniel's approach can be saved, if what you want is to express your expression entirely in terms of power sums. Since the power sums are always expressible in terms of the elementary symmetric polynomials (by Newton-Girard), you should be able to always process the output of SymmetricReduction[] using GroebnerBasis[]. To wit,

poly[n_Integer] := Sum[Subscript[x, r]*((Sum[Subscript[x, i], {i, 1, n}] -
Subscript[x, r])^2 - Sum[Subscript[x, i]^2, {i, 1, n}]), {r, 1, n}]

With[{n = 3},
First @ GroebnerBasis[
Prepend[Table[
p[k] - First[SymmetricReduction[Sum[Subscript[x, i]^k, {i, n}],
Table[Subscript[x, i], {i, n}], Table[s[j], {j, n}]]],
{k, n}],
First[SymmetricReduction[poly[n],
Table[Subscript[x, i], {i, n}], Table[s[j], {j, n}]]]],
Table[p[j], {j, n}], Table[s[j], {j, n}]]]
p[1]^3 - 3 p[1] p[2] + p[3]

-

Not entirely clear what you mean, but here is a solution that works if the structure of your equation fulfills certain constraints, but it will not be looking for mathematical equivalence in general.

Sum[Subscript[x, r] ((Sum[Subscript[x, i], {i, 1, n}] - Subscript[x, r])^2 -
Sum[Subscript[x, i]^2, {i, 1, n}]), {r, 1, n}]
/. {Sum[Subscript[_, i], {i, 1, n}] -> s1, Sum[Subscript[_, i]^2, {i, 1, n}] -> s2}
`

-
plz see my updates, I would like to have an equation without any x –  colinfang Mar 10 '13 at 23:53
@colinfang You provided a substitution rule for the literal sums over i only. You can't use the same substitution for the sum over the r part, as that can't be equal to s1. –  Sjoerd C. de Vries Mar 10 '13 at 23:59
hmm, ye, So I am asking how I can use mathmatica to help me solve/eliminate this maths equations, rather than myself doing it on paper –  colinfang Mar 11 '13 at 0:09