function construction - Improve timings for calculating coeficients

Is there a quicker way of finding ther coefficients of something like

$$(x_{1}^{2}+x_{2}^{2}+x_{3}^{2}+x_{4}^{2}+x_{5}^{2}+x_{6}^{2}+x_{7}^{2}+x_{8}^{2}+x_{9}^{2}+x_{10}^{2})\times\\ (x_{1}+x_{2}+x_{3}+x_{4}+x_{5}+x_{6}+x_{7}+x_{8}+x_{9}+x_{10})^{8}$$

I have used

  s10 = Table[ToExpression["Subscript[x, " <> ToString[i] <> "]"], {i, 10}];
  s10a = CoefficientRules@Expand[Plus @@ (s10^2)*(Plus @@ s10)^8];
  Table[s10a[[n]], {n, Flatten[Position[s10a, #] & /@ (PadRight[#, 10, 0] & /@ 
   IntegerPartitions@10), 1][[All, 1]]}]

to find all of the coefficients of the distinct monomial terms using Mr.WIzard's suggestion below. Can I improve on this?

Answer

If interested in the evolution of this, read on, otherwise, skip to near the bottom for latest iteration...

Give this a whirl, ran fine on a junky old netbook with 1 Gig that I keep at the local cigar shop:

var = 10
pow = 12

result = With[{ar = Array[x, var]}, CoefficientArrays[Tr[ar^2]*Tr[ar]^pow][[-1]] //
              (#["NonzeroValues"][[Ordering[#["NonzeroPositions"]]]]) &];

The results are quite sparse (e.g., for 10 variables and power of 12, only about 0.00000082% of the array is non-zero), so this should get you somewhere with memory pressure. It is quite a bit faster in my tests on the turdbook, YMMV.

If the order of results is not relevant (that is, you only care about the coefficient values), drop the ordering part, that will save time and some memory. The results with ordering applied follow the permutations of the partitions of the degree of the result (that is, the order that say CoefficientRules poops it out), so if you want to decorate the results with the coefficient vector, just generate them (trivial) and transpose...

Update: Here's a direct way to generate the coefficient for a specific exponent configuration. Quite quick, so for huge arguments, you can just iterate the list of valid exponent configurations and use it to get the coefficient. It is specific to your OP (squared sum of vars multiplied by some power of the sum of same vars), can easily be changed to do other setups. No sanity checking done - if you use it, either add some, or just be sure the exponent lists make sense...

The code:

coeff = With[{c = Clip[#, {2, 0}, {0, -2}], t = #, l = Length@#},
             Tr@(Multinomial @@@ Map[(t + UnitVector[l, #]*c) &, Pick[Range@l, c, -2]])] &;

An example of use (just for example, I generate some valid example coefficient list for vars variables and pow power, check it for validity, then use the function):

var = 50;

pow = 200;

exponents = RandomSample@
             Nest[Append[#, RandomInteger[{0, Ceiling[(pow + 2 - Tr[#])/3]}]] &, {}, var]

Through[{Length, Tr}[exponents]] == {var, pow + 2}

coeff[exponents] // Timing

Resulting in:

{30, 0, 1, 0, 0, 4, 1, 0, 1, 18, 2, 0, 0, 0, 0, 1, 0, 1, 11, 8, 0, 0, 0, 2, 0, 3, 2, 0, 10, 1, 24, 0, 0, 4, 1, 3, 17, 0, 1, 0, 0, 10, 1, 1, 0, 33, 8, 0, 2, 1}

True

{0., 2993477628373446010341511351677305932386459268334635334294479632840326025742434118385910676527613671812233157327614068574441257146217469230067380916523101753985716500952966641450100671258052610349229670400000000000000000}

So, below timer threshold on my lounging netbook. Near nil memory needs (basically just the storage needed for results if you iterate coefficients, or storage of results + list of coefficient lists if you generate the latter in bulk) and wildly parallelizable, this could handily beat the earlier method on many-core and/or distributed kernels.

Here's a flash I had re: using the above update. Since I don't know from your OP what you're doing with the coefficients/exponents, not sure if it will be of use... but:

var = 20;
pow = 25;


result = With[{ip = IntegerPartitions[pow + 2, {var}, Range[0, pow + 2]] // 
                 Pick[#, Unitize@Subtract[Max /@ #, 1], 1] &, r=Range@var},
           {#, With[{c = Clip[#, {2, 0}, {0, -2}], t = #},
            Tr@(Multinomial @@@ Map[(t + UnitVector[var, #]*c) &, Pick[r, c, -2]])],
            var!/Times @@ (Tally[#][[All, 2]]!)} & /@ ip];

Length@result

Total@result[[All, 3]]


RandomSample[result, 5]


(*

2980

4154246671960


{{{10,4,4,2,2,2,1,1,1,0,0,0,0,0,0,0,0,0,0,0},111314418415200000,846518400},
{{13,4,3,3,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0},518948337600000,105814800},
{{10,4,3,2,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0},1632611470089600000,1496523600},
{{7,7,6,3,3,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0},2968384491072000,6976800},
{{10,4,4,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0},845989579955520000,83140200}}

*)

So, we see that for the case of 20 variables, power of 25, there are 2980 distinct coefficients over 4,154,246,671,960 monomials. The snippet of the output shows the form of the exponents, the coefficient for that form, and the number of monomials with that form.

A "form" is just a permutation of the form in the output, e.g., say a simpler case had an output including {{2,3,4}, x, y}. That means all monomials with coefficients that are a permutation of {2,3,4} have coefficient x, and there are y of them.

This also means of course you can just take some arbitrary exponent form, sort it, and look up its sorted form in the output to get the details (or not, if it is an impossible exponent combination).

The 20/25 example took under a second on the goofy loungebook, so I'd guess a few hundredths of a second on a workstation...

The explanation:

(* Here's what we'll deconstruct *)
var = 4
ipow = 2
pow = 2

result = With[{ip = IntegerPartitions[pow + ipow, {var}, Range[0, pow + ipow]] // 
                Pick[#, UnitStep@Subtract[Max /@ #, ipow], 1] &, r = Range@var},

          {#, With[{c = Clip[#, {ipow, 0}, {0, -ipow}], t = #}, 
           Tr@(Multinomial @@@ Map[(t + UnitVector[var, #]*c) &, Pick[r, c, -ipow]])], 
           var!/Times @@ (Tally[#][[All, 2]]!)} & /@ ip]

(* and the direct constructions for comparisions *)
ar = Array[x, var];
argl = Tr[ar^ipow];
argr = Tr[ar]^pow;
arg = argl*argr;
cr = CoefficientRules[arg];

Note there's a couple of changes, one to generalize it, the other was an unneeded assignment. Throughout, var is the number of variables for each side of your construction, pow is the exponent for the right side, and ipow is the "internal" exponents for the left side.

Using the values above, let's look at the coefficient rules and see if there's anything interesting:

{{4, 0, 0, 0} -> 1, {3, 1, 0, 0} -> 2, {3, 0, 1, 0} -> 2, {3, 0, 0, 1} -> 2, 
 {2, 2, 0, 0} -> 2, {2, 1, 1, 0} -> 2, {2, 1, 0, 1} -> 2, {2, 0, 2, 0} -> 2, 
 {2, 0, 1, 1} -> 2, {2, 0, 0, 2} -> 2, {1, 3, 0, 0} -> 2, {1, 2, 1, 0} -> 2, 
 {1, 2, 0, 1} -> 2, {1, 1, 2, 0} -> 2, {1, 1, 0, 2} -> 2, {1, 0, 3, 0} -> 2, 
 {1, 0, 2, 1} -> 2, {1, 0, 1, 2} -> 2, {1, 0, 0, 3} -> 2, {0, 4, 0, 0} -> 1, 
 {0, 3, 1, 0} -> 2, {0, 3, 0, 1} -> 2, {0, 2, 2, 0} -> 2, {0, 2, 1, 1} -> 2, 
 {0, 2, 0, 2} -> 2, {0, 1, 3, 0} -> 2, {0, 1, 2, 1} -> 2, {0, 1, 1, 2} -> 2, 

 {0, 1, 0, 3} -> 2, {0, 0, 4, 0} -> 1, {0, 0, 3, 1} -> 2, {0, 0, 2, 2} -> 2, 
 {0, 0, 1, 3} -> 2, {0, 0, 0, 4} -> 1}

Let's group them by exponent set, regardless of a set's order:

gb = Gather[cr, Sort[#1[[1]]] == Sort[#2[[1]]] &];

Column[gb, Automatic, 1]

(*


{
 {{{4, 0, 0, 0} -> 1, {0, 4, 0, 0} -> 1, {0, 0, 4, 0} -> 1, 
   {0, 0, 0, 4} -> 1}},

 {{{3, 1, 0, 0} -> 2, {3, 0, 1, 0} -> 2, {3, 0, 0, 1} -> 2, 
   {1, 3, 0, 0} -> 2, {1, 0, 3, 0} -> 2, {1, 0, 0, 3} -> 2, 
   {0, 3, 1, 0} -> 2, {0, 3, 0, 1} -> 2, {0, 1, 3, 0} -> 2, 
   {0, 1, 0, 3} -> 2, {0, 0, 3, 1} -> 2, {0, 0, 1, 3} -> 2}},

 {{{2, 2, 0, 0} -> 2, {2, 0, 2, 0} -> 2, {2, 0, 0, 2} -> 2, 

   {0, 2, 2, 0} -> 2, {0, 2, 0, 2} -> 2, {0, 0, 2, 2} -> 2}},

 {{{2, 1, 1, 0} -> 2, {2, 1, 0, 1} -> 2, {2, 0, 1, 1} -> 2, 
   {1, 2, 1, 0} -> 2, {1, 2, 0, 1} -> 2, {1, 1, 2, 0} -> 2, 
   {1, 1, 0, 2} -> 2, {1, 0, 2, 1} -> 2, {1, 0, 1, 2} -> 2, 
   {0, 2, 1, 1} -> 2,{0, 1, 2, 1} -> 2, {0, 1, 1, 2} -> 2}}
}

*)

We can see that for a given distinct set of exponent components, the coefficients are the same for all permutations of that set. Makes sense algebraically. So, let's chop them to just the root permutation of each distinct set:

{4,0,0,0}->1
{3,1,0,0}->2
{2,2,0,0}->2
{2,1,1,0}->2

So, if we can construct a black box that given one of the valid distinct exponent sets it returns the coefficient, we've cut down the work hugely already - instead of feeding it all possible exponent sets, we need only distinct sets, an exponentially smaller group.

Now, that pattern looks kind of familiar:

ip = IntegerPartitions[pow + ipow, {var}, Range[0, pow + ipow]];


ip // Column

(*
{4,0,0,0}
{3,1,0,0}
{2,2,0,0}
{2,1,1,0}
{1,1,1,1}
*)

Hmm, close... we know that the sum of exponents of each term of your construction must be pow+ipow. But we also know from algebra that the maximum exponent must be at least ipow. So, chop the ones not meeting that requirement:

ip = IntegerPartitions[pow + ipow, {var}, Range[0, pow + ipow]] // 
   Pick[#, UnitStep@Subtract[Max /@ #, ipow], 1] &;

ip // Column
ip == gb[[All, 1, 1]]

(*
{4,0,0,0}
{3,1,0,0}

{2,2,0,0}
{2,1,1,0}

True
*)

Looks good, and matches the results from the MMA CoefficientRules.

A quick side-trip re: Pick, etc. What is being done there is I get a list of the maximum value from each integer partition, subtract ipow from it, and do a UnitStep on the result. This then allows me to pick the integer partitions that meet the criteria - UnitStep gives a 1 for hits, 0 for misses, and pick returns only the partitions with 1 in the corresponding pick mask. You could use position and other methods, but when lists get larger, this is much, much faster.

OK, we've got the exponent sets that are valid, now what? Well, it turns out our friend algebra helps again. It turns out that the coefficient for an exponent set is nothing more than the sum of the multinomial coefficients for all possible combinations of the set where ipow can be subtracted from an element and leave a non-negative number.

So, how to get at this? Clip. It gets us not only the places in a set where the condition is met, but in the same swoop generates pieces we'll need later:

clipped = Clip[#, {ipow, 0}, {0, -ipow}] & /@ ip;

clipped // Column

(*
{-2,0,0,0}
{-2,0,0,0}
{-2,-2,0,0}
{-2,0,0,0}
*)

So, we basically say "Anything that's not at least ipow, make it zero, otherwise, make it the negative of ipow..."

Putting the result side-by-side with the valid exponent sets might make it clearer:

{-2,0,0,0}  {4,0,0,0}
{-2,0,0,0}  {3,1,0,0}
{-2,-2,0,0} {2,2,0,0}
{-2,0,0,0}  {2,1,1,0}

See how the elements of each exponent set that are at least ipow are "marked" with -2 in the result?

So, we want to (for each exponent set) get the result(s) from each possible valid subtraction of ipow:

r = Range@var;

tmp = (t = #; c = Clip[#, {ipow, 0}, {0, -ipow}];
    Map[(t + UnitVector[var, #]*c) &, Pick[r, c, -ipow]]) & /@ ip;

Row[{clipped // Column, ip // Column, tmp // Column}, "  "]

(*
{-2,0,0,0}   {4,0,0,0}  {{2,0,0,0}}
{-2,0,0,0}   {3,1,0,0}  {{1,1,0,0}}

{-2,-2,0,0}  {2,2,0,0}  {{0,2,0,0},{2,0,0,0}}
{-2,0,0,0}   {2,1,1,0}  {{0,1,1,0}}
*)

So we see the mask, the set, and the result(s) of doing the subtraction on all valid places.

A note on Pick and UnitVector: Again, this formed to be as quick as possible. In the pick, I have the list of possible positions {1,2,3...var} in r. So it returns only the rs that have a negative ipow in the corresponding position in c, the clipped partition. As before, position/etc. could be used, but this is generally much more efficient for this kind of thing.

The results of the pick are mapped over, so I take the partition (in t), use UnitVector to "mask the mask" to each of the valid positions, and add the terms (since this is just a -ipow, that accomplishes the subtraction).

Again, this could be done many ways, this is just a pretty quick way. Might be room for improvement - honestly, I just barfed this out while lounging, there may well be a way to do this "bottom-up", or a more direct way...

OK, so recall we want the sums of the corresponding multinomial coefficients:

mults = (t = #; c = Clip[#, {ipow, 0}, {0, -ipow}];

     Multinomial @@@ Map[(t + UnitVector[var, #]*c) &, Pick[r, c, -ipow]]) & /@ ip;

gets us what we want - it applies Multinomial to all results from the mapping.

Let's look at the above all together (sorry, as image, formatting got weird):

Row[Column /@ {clipped, ip, tmp, mults, Total /@ mults, validExp}, "    "]

So we see indeed, as the math would tell us, the sum of the multinomials of the candidates for a given exponent set is the coefficient for that set.

In the function, all of this is done in one place:

{#, With[{c = Clip[#, {ipow, 0}, {0, -ipow}], t = #}, 

    Tr@(Multinomial @@@ Map[(t + UnitVector[var, #]*c) &, Pick[r, c, -ipow]])], 
     var!/Times @@ (Tally[#][[All, 2]]!)} & /@ ip

So we get a list of the exponent set, the resulting coefficient, and the count of monomial terms matching (the latter is just the combinatorial math - (number of slots)!/(Product of factorials of numbers of each distinct element). Tally is used to count how many of each, just the count is kept, factorial done on each, and then multiplication of all together.

The speed comes from minimizing the work to be done up front, and then using efficient mechanisms within wherever "listy" operations occur on potentially large lists.

The next "bottleneck" would probably be the generation of the partitions - given large enough arguments, there will be a correspondingly large number of them. No shortcuts around that - you need them all.

Should you get the that point, you can simply generate the partitions iteratively, either with roll-your-own code, or using something like the NextPartition function in the built-in Combinatorica package. Past that, you'd probably run into limits of how large of a result MMA can handle, or the universe will be in heat death...

Lastly, here's the raw function adapted to arbitrary "internal" (left side) powers - what you'd use if you're iterating valid coefficient sets or want singlet results:

coeff = With[{c = Clip[#, {#2, 0}, {0, -#2}], t = #, l = Length@#, d = #2},
    Tr@(Multinomial @@@ Map[(t + UnitVector[l, #]*-d) &, Pick[Range@l, c, -d]])] &;

It takes the same first argument (the exponent set) and a second argument, that is the power all of the left-side variables are raised to, e.g,. cl here is a random list I made representing a 500 variable exponent set, with the "internal" exponents 100 and the "outer" exponent 25000...

IntegerLength[coeff[cl, 100]] // Timing

(* {0.246003, 64032} *)

So, a quarter second on my loungebook for a result that has ~2^16 digits - I'd venture a few thousandths of a sec. on one of my workstations.

Well, I think that's it, feel free to comment if further clarification is needed.

Update: A cleaner method, ~2-6X faster on cases I tested. N.b.: I removed the "count" column, you can add it back in if desired, speed comparisons were done with the same removal from the prior method to have a level playing field...

var=10

pow=15
ipow=2

newresult = With[{ip = IntegerPartitions[pow + ipow, {var}, Range[0, pow + ipow]] // 
                  Pick[#, UnitStep@Subtract[Max /@ #, ipow], 1] &},
               Transpose[{ip, Tr/@Divide[Multinomial @@@ ip, Divide[FactorialPower[pow + ipow, ipow], 
                          FactorialPower[Pick[ip, UnitStep@Subtract[ip, ipow], 1], ipow]]]}]]

Here's a quick benchie I ran while raiding the humidor. var set to 10, ipow to 2, pow ranging from 2 to 12. On the loungebook, so figure timings on workstation with advantages of outright speed and many more cores s/b order(s) of magnitude faster for all.

We can see the exponential growth (expected) from the work required to get the coefficients, compounded by the repeated lookups in your OP method. My OP method shows the benefit of using the more efficient means in MMA, but still has exponential growth, just with a big head start. Lastly, the new method (last of the ones in this post) shows... well, I'm not sure. Most of the times were below timer resolution, the last one being a couple of a hundredth of a second (~25,000X faster than OP).

Finally (I think...), here's an optimization I cooked up while optimizing the tummy at dinner. We can take advantage of the structure of the valid exponent sets, and when the number of variables is large, drop exponents that contribute nothing to the results:

newerresult = 
 With[{ip = IntegerPartitions[pow + ipow, {Min[var, pow + 1]}, Range[0, pow + ipow]] // 
            Pick[#, UnitStep@Subtract[Max/@#, ipow], 1] &},
      Transpose[{ip, Tr/@Divide[Multinomial@@@ip, 
                 Divide[FactorialPower[pow + ipow, ipow], 
                        FactorialPower[Pick[ip, UnitStep@Subtract[ip, ipow], 1], ipow]]]}]]

With, for example, var=500, pow=40, ipow=10, this carries about ~1/20 the memory footprint for the exponent sets as the last, and (due to chopped calculation, and probably more cache hits) is ~32X faster yet again over the last. I'm pretty pleased with that!

(Do note, the returned coefficient forms will be truncated to the length of the highest cardinality of non-zero exponents, so just pad it with zeros to length=number of variables if you need it fleshed out. I did it this way so memory is conserved even in the result.)

I'd say we've made some solid progress on this interesting puzzle!