# Optimising coefficient finding functions

I'm seeking a generalisation to this question, where I am looking to list the coefficients (without repetition) to monomials of the form

$$(x_{1}^{3}+x_{2}^{3}+x_{3}^{3}+x_{4}^{3})^4\cdot(x_{1}^{2}+x_{2}^{2}+x_{3}^{2}+x_{4}^{2})\cdot(x_{1}+x_{2}+x_{3}+x_{4})^{2}$$

for example.

Motivation:

I am currently using the following functions:

ip[n_] := Times @@@ Flatten[Table[{Table[Prime@ k, {k,
Length[IntegerPartitions[n][[s]]]}]}^{IntegerPartitions[n][[s]]},
{s, Length[IntegerPartitions[n]]}], 1]

alpha1[omega_] :=  Join[{a^omega},   Times @@ # & /@ Rest@(List @@ # & /@
(If[# == 1, 1, If[PrimePowerQ[#] == True, Power @@ {a, FactorInteger[#][[1, 2]]},
(Power @@@Transpose@{Take[ToExpression[CharacterRange["a", "z"]],
Length@FactorInteger@#], (Reverse@Sort@FactorInteger[#][[All, 2]])})]] & /@ip@omega))]

timesABC[omega_] := Module[{
aa = Flatten[#] & /@ (DeleteDuplicates[#] & /@ Part[#] & /@ (Split[#] & /@
(Reverse@Flatten[#] & /@ (Split[#] & /@ IntegerPartitions@omega)))),
bb = Table[Length[#] & /@ ((Split[#] & /@ (Reverse@Flatten[#] & /@ (Split[#] & /@
IntegerPartitions@omega)))[[nn]]), {nn,Length@IntegerPartitions@omega}],
cc = Table[Plus @@ (Take[ToExpression[CharacterRange["a", "z"]], omega]^nn), {nn,omega}]},
Times @@ # & /@ Table[(cc[[#]] & /@ aa[[#]] & /@ Range@Length@aa)[[nn, #]]^
bb[[nn, #]] & /@ Range@Length@bb[[nn]], {nn, Length@bb}]]


There are no problems with these functions - they all work fine. The problem is with this:

partialON[omega_] := Flatten[#] & /@ ((Part[#, 1] & /@ Part[#] & /@
Module[{ff = Expand[#] & /@ timesABC[omega]},
Table[Transpose@{Table[DeleteCases[Coefficient[List @@ ff[[m]],
alpha1[omega][[n]]], 0], {n, Length@alpha1@omega}], alpha1[omega]},
{m, Length@IntegerPartitions@omega}]]) /. {} -> {0})

partialON@4

(*out

{{1, 0, 0, 0, 0}, {1, 1, 0, 0, 0}, {1, 0, 2, 0, 0}, {1, 2, 2, 2, 0}, {1, 4, 6, 12, 24}}

*)


What this function does is probably best illustrated in table form: where the top row are the product inputs, and the rows are the coefficients of the right column.

The function partialON works as desired, but is not optimised, hence the question. I can apply rasher's solution to the top column of the table, but only for two groups of parentheses at present. I would ideally like to be able to go as high as partialON@20. Currently, 12 is stretching things on my little machine.

Just for clarity, the full code for generating above type tables:

aZ[omega_] := Module[{aa = Table["(" <> ToString[Superscript[Subscript[x, 1],
If[n == 1, "", n]] // TraditionalForm] <> "+.. " <> ToString[Superscript[Subscript[
x, omega], If[n == 1, "", n]] // TraditionalForm] <> ")", {n, 40}],
bb = Table[Plus @@ (Take[ToExpression[CharacterRange["a", "z"]], omega]^n),
{n, omega}],
cc = Reverse@(Flatten[#] & /@ (DeleteDuplicates[#] & /@ Part[#] & /@ (Split[#] & /@
(Reverse@Flatten[#] & /@ (Split[#] &  IntegerPartitions@omega))))),
dd = Reverse@Table[Length[#] & /@ ((Split[#] & /@ (Reverse@Flatten[#] & /@ (Split[#]
& /@ IntegerPartitions@omega)))[[n]]), {n,Length@IntegerPartitions@omega}]},
Column@# & /@ Table[aa[[Part[#] & /@ cc[[n]]]]^(Part[#] & /@ dd[[n]]),
{n, Length@dd}]]
alphax[n_] := ToString[If[n == 1, 1, If[PrimePowerQ[n] == True,Power @@
{Subscript[x, 1], FactorInteger[n][[1, 2]]}, Times @@ (Power @@@ Transpose@{Table[
ToExpression["Subscript[x, " <> ToString[i] <> "]"], {i, Length@FactorInteger@n}],
(Reverse@ Sort@FactorInteger[n][[All, 2]] )})]] // TraditionalForm]
oG3[omega_] := Text@Grid[Join[{Join[Rotate[#, 0] & /@ (Reverse@Flatten@(aZ[omega])),
{""}]}, Transpose@ Join[partialON[omega] /. 0 -> "", {alphax[#] & /@ ip@omega}]],
Frame -> All, Background -> {None, {{White,  Lighter@Lighter@Lighter@Lighter@Gray}}},
Alignment -> {Left, Center},   ItemSize -> {Join[
Table[7, {x, Length@IntegerPartitions@omega}], {Automatic}]}]

oG3@5

• Nice question, +1. I recall this work deals with some partition experimentation, perhaps add info. on what these structures are (e.g., "I'm looking at counting how many ways foo in sets of bar can be partitioned into baz sets...") might lead to more compact solutions.
– ciao
Feb 16, 2015 at 22:24
• OK, a quick stab at this netted 30-100X faster than CoefficientRules on 5 variables, "internal" powers between 2 and 5, "external" powers between 2 and 5. I'm verifying correctness (using a different method, I don't think the mathematical trickery used earlier applies here). Rather than decode the code and the paper involved (interesting article, BTW), can you specify some useful to you upper bounds on number of vars, maximum "internal" and "external" powers, and if the powers are disjoint, or disjoint for internal/external, or.... this will help with deciding on which routes to pursue.
– ciao
Feb 17, 2015 at 7:40
• Also, might be worthwhile to query on Math.StackExchange - see if someone there is expert in the area and might have a short-cut method for this kind of coefficient gathering. I've also got a sneaking suspicion the group theoretical functions in Mathematica might have a use here.
– ciao
Feb 17, 2015 at 7:43
• @rasher that sounds great - 30-100X sounds good to me :) re MSE, I have asked previously, but will try to word better. Meanwhile I will look up groupfunctions on MMA. I don't know re variables - I have just gone with that safe (and time consuming) option where var=ipow+pow using your definitions from the last question. Feb 17, 2015 at 9:35
• See the non-answer answer for now, need to hit the hay, will ponder the core problem and the rest of your code this week.
– ciao
Feb 17, 2015 at 10:06

This is not (at this time) a direct answer to the core question, perhaps better thought of as a placeholder. Nonetheless, it will serve as an example of how to speed up this problem, and like things in general.

You need to start thinking in "Mathematica" terms. That is, whenever possible, think of how you might manipulate things en masse, as operations on vectors, matrices, lists, etc. If you can vectorize things in Mathematica, it;s nearly always a win.

By way of example, I started to decode your code in the OP to get a better understanding of what you're after. Here's the first function, Mathematica-ized:

newip[n_] := Times @@@ (Prime@Range[Length /@ #]^#) &[IntegerPartitions[n]];


How much of a difference does that make? We have to view it in Log scale to even see the faster way: By the time N is 30, this is over 250X faster.

I'll update this as I go through your code, and if the current ideas for the core problem flesh out, I'll add them.

Update 24/02/2015:

OK, here's the coefficient finder I've been tooling with.

Some caveats: It's obviously a work-in-progress, parts within that could be combined are separated to allow profiling of them individually, some parts certainly could be further optimized (e.g. using outer vs mapping over sets for low cardinality partition sets), etc., but perhaps there's some useful stuff here for you.

I've not tested it on huge arguments compared to native Mathematica methods - I do my play puzzling on leisure netbook & laptop, so no testing on real workstations was done...

This takes arguments similar to the earlier answer in the linked question, e.g.

(x + x + x + x + x + x)^10 (x^2 + x^2 +
x^2 + x^2 + x^2 + x^2)^10 (x^3 + x^3 + x^3 +
x^3 + x^3 + x^3)^2


is represented with the first two arguments being the "internal" and "external" powers:

{1, 2, 3, 2, 2} , {10, 1, 2, 4, 5}


the third argument is the exponent set you're after, e.g.:

{17, 9, 7, 1, 1, 1}


representing an exponent set in six variables. You can use things like

{17, 9, 10} or  {17, 9, 10, 0, 0, 0}


for terms with zero exponents, they are treated equivalently.

A quick test on the above case gives this an ~2500X speed advantage over your coeff, and about a 10X advantage over Coefficient in Mathematica. I'd expect the gap to grow as problem size grows. It is also much gentler on memory resources.

As stated in my comment, getting all is still faster using something like CoefficientArray (I've yet to find the Kryptonite for this problem, any luck for you over at Math.StackExchange?), but I'm guessing that's a no-go for you due to memory or other resource pressure?

getCoeff[ipow_, epow_, expset_] :=
Module[{f, fastRF, myIP, px = {{}}, target = DeleteCases[expset, 0],
set = Join @@ MapThread[ConstantArray, {ipow, epow}], nres, x, bin,
tset, breaks, accs, splits, seqs, res},

(* fugly patch for now *)
If[(Tr[ipow*epow] < Tr@target) || (Min@target < Min@ipow), Return];

fastRF[a_List, b_List] := Module[{c, o, x}, c = Join[b, a];
o = Ordering[c];
x = 1 - 2 UnitStep[-1 - Length[b] + o];
x = FoldList[Max[#, 0] + #2 &, x];
x[[o]] = x;
Pick[c, x, -1]];

myIP[set_, used_, left_, n_] :=
With[{tmp = fastRF[set, used]},
(* Another fugly patch for non-conformance *)
If[Length@tmp < left, Return[{{}}]];
IntegerPartitions[n, Length@tmp - left, DeleteDuplicates@tmp]];

tset = Tally[set];

Module[{cnt = #[], tmp, reset = #[]},
bin[#[], x_] := (tmp = Binomial[cnt, x]; cnt -= x;
If[cnt == 0, cnt = reset]; tmp)] & /@ tset;

(f[#[], z_] = z <= #[]) & /@ Tally[set];

nres = Module[{z = 0},
Nest[With[{zz = #}, z++;
Join @@ (Module[{mip = myIP[set, #, Length@target - z, target[[z]]], tmp = #},
mip = Join[tmp, #] & /@ mip;

Pick[mip, Apply[And, Apply[f, Tally /@ mip, {2}], {1}]]] & /@
zz)] &, px, Length@target]];

breaks = Accumulate@target;
accs = Accumulate /@ nres;

Reap[Do[
x = Join @@ Position[accs[[z]], Alternatives @@ breaks];
x = Differences[Prepend[x, 0]];
splits =
Transpose@{Span @@@ Transpose[{Accumulate[Prepend[Most@x, 0]] + 1, Accumulate[x]}]};
segs = Tally /@ Rest@Extract[nres[[z]], Prepend[splits, {{}}]];
res = Times @@ Join @@ Apply[bin, segs, {2}];
Sow[res];, {z, Length@nres}]] //
If[#[] == {}, 0, Tr@#[[2, 1]]] &

]


A quick run-through:

Sub-function myIP is a restricted partition generator. It takes size and member restrictions and is used to build the candidates for each of the target exponent set components. It uses Mr. Wizard's fastRF to cull the restrictions. I may have a faster method, but it's an inconsequential part of the process, and his code is so pretty, so I kept it.

I the build a set of custom binomial coefficient functions bin[...]. These return the coefficient given the target value but keep track internally of how many of that value are available, resetting to the initial tally when exhausted for each candidate partition set.

Following, I build a set of f[...] functions that allow quick validation of the tallies of candidate partition sets (i.e., does the set have a valid amount of each possible piece).

The next step fills nres with valid partition sets using the above.

Finally, I parse the candidate partition sets, splitting them at the points where the running "stuttering" totals equal the target coefficients, tally these, extract the results using an undocumented form of extract I discovered for speed, and applying the adaptive bin binomial coefficient extractors. Lastly, I check if the result is empty (coefficient set was invalid) and return 0 (invalid) or the total of the coefficients (valid).

I hope at least some piece of this finds use, and will continue pondering this interesting puzzle (I'm convinced there's a better way, mathematically justified, that perhaps someone with deep multinomial/poly-product knowledge can answer with over at your Math SE query).

Looking at your test harness, there's an obvious issue - you can't just Quiet something and expect that errors sans messages means the same performance. Even though using Quiet does usually reduce timings when errors cause messages, there's still overhead in Mathematica handling the issue, and of course the likelihood that there's a cascade of problems presented to the code... also on large result lists, using rules can be quite inefficient for replacements. Look into Dispatch and using Replace directly with a limited level specification if you have tasks that need to do such things in the future.

Since the original question, comments, Math.SE question, etc. all implied the feeding of well-conditioned exponent lists, this was built with that in mind. I've added another ugly patch for now, but be advised, it just short-ciruits the problem, but it also means there's work, sometimes a lot, getting done where it need not be.

That said, I took your code, and just changed one line to be:

getCoeff[Sequence @@ Transpose@Tally@list1, #] & /@ list2


in your partialONk function. No other changes made, so there's extra work converting the argument format, and the extra work being done in my function because bad sets are passed to it. Nonetheless, the following is the result (netbook timings, I ran out of patience at 12): Both grow quickly (expected from the problem), but the original coeff explodes.

There appears to be much overhead from the rest of your code, based on just timing the coefficient functions. I'd venture re-thinking it, and in particular generating only valid exponent sets, will go a long way in improving timing. As it stands, much of the boost in finding the coefficients is getting swamped in other overhead. By example, the two functions, over a range of 1 to 5 components, each with random powers between 2 and 5 for internal/external power, and a random exponent set of appropriate size. I gave up at six components - getCoeff finished in under a second, I aborted coeff after twenty minutes. At five components, getCoeff was ~2800 times faster... N.B. This had to be done with a Log scale so getCoeff would even show up on the plot: • yes, much cleaner (and faster!) :) I'll continue to look through the group theory stuff too, though not turned anything up yet that's too useful ... Feb 17, 2015 at 20:37
• @martin Spent some more time on this, frankly drawing blanks at this point - method I thought might be useful turns out not to be, so back to drawing board. Highly recommend pursuing the core question (fast coefficient determination given an exponent set since the latter can be generated easily). I'd venture there's someone there that might have the specific background to provide info or advise there's no real short-cut. I'll continue pondering...
– ciao
Feb 19, 2015 at 6:00
• ok - will follow up on MSE - will let you know if I get any joy over there. I have been looking at it & there certainly seems to be some sort of pattern, but it is eluding me for now too. Feb 19, 2015 at 9:36
• posted a question here, hope I haven't missed anything. Feb 19, 2015 at 10:20
• @martin See chat here, an idea that so far seems solid, not sure how practical, might lead to some revelation perhaps...
– ciao
Feb 20, 2015 at 6:00

Have managed to get sn improvement in the method given in the OP using rasher's reveltation in this chat session & Mr.Wizard's removeFrom

subs[num_, list_] :=
DeleteDuplicates[Sort[#] & /@ Select[Subsets[list], Total@# == num &]]

removeFrom[b_List, a_List] := Module[{f}, f[_] = 0;
(f[#] = -#2) & @@@ Tally[a];
Pick[b, UnitStep[f[#]++ & /@ b], 1]]

coeff[list1_, list2_] := With[{cc = subs[#, list1] & /@ list2},
With[{df =
DeleteCases[
If[removeFrom[Sort@list1, Sort@Flatten@#] == {}, #,
0] & /@ (Tuples[cc]), 0]},
Total[Times @@ # & /@
Table[With[{ff =
Join[{Flatten[Split[#] & /@ df[[nn]], 1],
Join[{list1},
Most@(removeFrom[list1, #] & /@
Table[Flatten[
Take[Flatten[Split[#] & /@ df[[nn]], 1], n]], {n,
Length@Flatten[Split[#] & /@ df[[nn]], 1]}])]}]},
Table[
Binomial[Count[ff[[2, m]], First@ff[[1, m]]],
Length[ff[[1, m]]]], {m,
Length@Flatten[Split[#] & /@ df[[nn]], 1]}]], {nn,
Length@df}]]]]

aZ2[omega_] :=
Module[{aa = Table[n, {n, 100}],
bb = Table[
Plus @@ (Take[ToExpression[CharacterRange["a", "z"]], omega]^
n), {n, omega}],
cc = Reverse@(Flatten[#] & /@ (DeleteDuplicates[#] & /@
Part[#] & /@ (Split[#] & /@ (Reverse@
Flatten[#] & /@ (Split[#] & /@
IntegerPartitions@omega))))),
dd = Reverse@
Table[Length[#] & /@ ((Split[#] & /@ (Reverse@
Flatten[#] & /@ (Split[#] & /@
IntegerPartitions@omega)))[[n]]), {n,
Length@IntegerPartitions@omega}]},
Table[{aa[[Part[#] & /@ cc[[n]]]], (Part[#] & /@ dd[[n]])}, {n,
Length@dd}]]

partialONk[omega_] :=
With[{list2 =
Table[Flatten[
Join[Table[#[], {xx, #[]}]] & /@
Transpose@(Reverse@aZ2[omega])[[yy]]], {yy,
Length@aZ2[omega]}]},
Table[With[{list1 = list2[[zz]]},
coeff[list1, #] & /@ list2], {zz, Length@list2}]]

partialONk


which can calculate partialONk without too much trouble (though timing is not great). The code could certainly benefit from some optimizing though!