A fast way to expand symbolic products

Question

My goal is to expand a product, where the factors are sums of symbols such as: $a_1 \times a_2 \times \ldots \times a_n$, where $a_i = b_1 + b_2 + \ldots + b_m$. The desired final expression is comprised of terms with the form $b_{i_1}b_{i_2}\ldots b_{i_n}$. The point is that I have defined an algebra between the $b_i$ symbols with upvalues and I need to fully expand the initial expression to symplify it with the defined algebra.

The final number of terms is of order $\mathcal{O}(10^5)$. The expansion usually takes 4-5 minutes (for the methods listed below). I have already tried

1. Expand,ExpandAll
2. Outer[Times,...]
3. Build a loop to successively expand the product factor by factor (expand the product of the first two factors, then multiply it with the next factor, expand it, etc.)

Is there a faster way to perform the expansion than any of the above listed options? A novel way would be a method which is an order of magnitude faster than the Expand command.

A minimal working example would look something like this:

fac[n_] := Sum[Product[Unique[x], {i, 1, RandomInteger[{1, 5}]}]  a[i], {i,1,n}]
expr = fac[#] & /@ ConstantArray[10, 6];

My first attempt (which is the fastest):

sol1 = Times @@ expr1 // Expand; // AbsoluteTiming
(*
{11.2204, Null}
*)

The second one:

sol2 = Total@
    Flatten@Outer[Times, expr[[1]], expr[[2]], expr[[3]], expr[[4]], 
      expr[[5]], expr[[6]]]; // AbsoluteTiming
(*
{31.4785, Null}
*)

And lastly:

test[list_] := With[{x = list},
  result = 1;
  For[var = 1, var <= Length@list, var++,
   result *=  list[[var]];
   result = Expand@result;
   ];
  result
  ]

sol3 = test@expr; // AbsoluteTiming
(*
{76.0036, Null}
*)

Answer 1

Not an answer, but just an attempt to analyze the problem. IMHO, what causes Mathematica to break down in this calculation is not the number of terms, but their complexity (i.e. high LeafCount). For example, we can try to substitute the products of Unique symbols with simple c[i]

fac[n_] := 
 Sum[Product[Unique[x], {i, 1, RandomInteger[{1, 5}]}] a[i], {i, 1, n}]
expr = fac[#] & /@ ConstantArray[10, 6];

counter = 0;
repFun[ex_Plus] :=
 (

  li = List @@ ex;
  Table[counter++; h[li[[i]]] -> c[counter], {i, 1, Length[li]}]
  )

The following code gives us the corresponding substitution rule

repRule = Flatten[repFun /@ expr]

which makes expr look like

exprAbbr = Map[h /@ #1 &, expr] /. repRule;

this

{c(61)+c(62)+c(63)+c(64)+c(65)+c(66)+c(67)+c(68)+c(69)+c(70),
c(71)+c(72)+c(73)+c(74)+c(75)+c(76)+c(77)+c(78)+c(79)+c(80),
c(81)+c(82)+c(83)+c(84)+c(85)+c(86)+c(87)+c(88)+c(89)+c(90),
c(91)+c(92)+c(93)+c(94)+c(95)+c(96)+c(97)+c(98)+c(99)+c(100),
c(101)+c(102)+c(103)+c(104)+c(105)+c(106)+c(107)+c(108)+c(109)+c(110),
c(111)+c(112)+c(113)+c(114)+c(115)+c(116)+c(117)+c(118)+c(119)+c(120)}

Now we can expand exprAbbr in about 2 seconds.

AbsoluteTiming[resRaw = Expand[Times @@ exprAbbr];]

However, when trying to perform backsubstitution we loose 16.5 seconds, and this is where the high LeafCount kicks in.

revRepRule = (Reverse /@ repRule) /. h -> Identity;
AbsoluteTiming[res = resRaw /. Dispatch[revRepRule];]

So at the end of the day it is true that the only sensible way to solve this problem is to switch to FORM.

In a way I find it somewhat sad, given the amount of effort put into development of Mathematica. However, it seems that good performance when handling very large expressions is not ranked very high in the priority list.

Edit:

In fact, even with FORM this expansion is not entirely trivial. Using

fac[n_] := 
 Sum[Product[Unique["x"], {i, 1, RandomInteger[{1, 5}]}] a[i], {i, 1, 
   n}]
expr = fac[#] & /@ ConstantArray[10, 6];
str = "L expr = " <> 
  StringReplace[
   ToString[Times @@ expr, InputForm], {"[" -> "(", "]" -> ")"}] <> 
  ";"

SetDirectory[NotebookDirectory[]];
file = OpenWrite["expr.frm"];
WriteString[file, str]
Close[file];

we can convert the input expression into FORM format and then evaluate it using following FORM code

Auto S x;
CF a;

#include expr.frm

.sort

Format Mathematica;

#write  <expr-expanded.m> "\"(%E)\"",expr

.end

Looking at the FORM output statistics

# time form formExpand.frm                                                                                                                                                                                              
FORM 4.1 (Aug 28 2019, v4.1-20131025-482-g0b3ab5d) 64-bits  Run: Fri Dec 27 19:31:39 2019
    Auto S x;
    CF a;

    #include expr.frm
    L expr = (x73*x74*x75*a(1) + x76*a(2) + x77*x78*x79*x80*x81*a(3) + x82*x83*x84*a
    (4) + x85*x86*x87*x88*x89*a(5) + x90*a(6) + x91*x92*a(7) + x93*x94*x95*x96*x97*a
    (8) + x100*x101*x98*x99*a(9) + x102*x103*x104*x105*x106*a(10))*(x107*x108*a(1) +
     x109*x110*a(2) + x111*x112*a(3) + x113*a(4) + x114*x115*x116*x117*a(5) + x118*x
    119*a(6) + x120*x121*a(7) + x122*x123*x124*a(8) + x125*a(9) + x126*x127*x128*x12
    9*a(10))*(x130*x131*x132*a(1) + x133*a(2) + x134*x135*a(3) + x136*x137*x138*x139
    *x140*a(4) + x141*x142*x143*x144*x145*a(5) + x146*x147*x148*x149*x150*a(6) + x15
    1*x152*x153*x154*x155*a(7) + x156*a(8) + x157*x158*a(9) + x159*x160*x161*a(10))*
    (x162*a(1) + x163*a(2) + x164*x165*a(3) + x166*a(4) + x167*a(5) + x168*x169*x170
    *x171*a(6) + x172*x173*x174*x175*a(7) + x176*a(8) + x177*x178*x179*x180*x181*a(9
    ) + x182*a(10))*(x11*x12*x13*a(1) + x14*x15*a(2) + x16*x17*x18*x19*x20*a(3) + x2
    1*x22*x23*a(4) + x24*x25*x26*x27*x28*a(5) + x29*x30*a(6) + x31*x32*x33*a(7) + x3
    4*x35*x36*x37*x38*a(8) + x39*a(9) + x40*x41*x42*a(10))*(x43*x44*x45*x46*a(1) + x
    47*x48*x49*x50*x51*a(2) + x52*x53*a(3) + x54*x55*x56*x57*a(4) + x58*x59*x60*x61*
    x62*a(5) + x63*x64*a(6) + x65*a(7) + x66*x67*a(8) + x68*x69*x70*x71*a(9) + x72*a
    (10));
    .sort

Time =       1.08 sec    Generated terms =     519379
            expr       1 Terms left      =     519379
                         Bytes used      =   84983048

Time =       2.05 sec    Generated terms =    1000000
            expr       1 Terms left      =    1000000
                         Bytes used      =  156406260

Time =       2.42 sec    Generated terms =    1000000
            expr         Terms in output =    1000000
                         Bytes used      =  154359092

    Format Mathematica;

    #write  <expr-expanded.m> "\"(%E)\"",expr

    .end

Time =       6.43 sec    Generated terms =     537298
            expr  537299 Terms left      =     537298
                         Bytes used      =   81653732

Time =       6.77 sec    Generated terms =    1000000
            expr 1000000 Terms left      =    1000000
                         Bytes used      =  154359216

Time =       7.00 sec    Generated terms =    1000000
            expr         Terms in output =    1000000
                         Bytes used      =  154359092
  7.00 sec out of 7.01 sec
form formExpand.frm  4,58s user 2,48s system 99% cpu 7,071 total

we see that FORM needs only 2.4 seconds for the expansion. But since you probably want to continue processing the expressions with Mathematica and not FORM, we need to export it into Mathematica format and then parse it from Mathematica. Here we loose about 4.2 seconds for writing the output to disk (the file is over 100 MB large!)

Then loading the result

AbsoluteTiming[
 res = Import["expr-expanded.m", "String"];
 res = ToExpression[res];]

still requires about 3 seconds. So there is a caveat here: If you do everything with FORM you can gain performance, but exchanging results between FORM and Mathematica costs time and should be, in general, minimized.

Of course, FORM has its own binary format for saving and loading results in an efficient way. Unfortunately, Mathematica cannot read that.

Notice also that if you try to naively parse the full expression directly without the String/ToExpression detour, you will loose more time. With

#write  <expr-expanded.m> "(%E)",expr

doing

AbsoluteTiming[res = Get["expr-expanded.m"];]

requires about 20 seconds!

Regarding FORM and Windows, you definitely need to install Cygwin to get it compiled. If you completely fail I can have a look at it later.

Answer 2

The question probably should be marked as a duplicate since it appears periodically, because Expand is so fundamental operation. And as it was noted by many peoples (including vsht in the comment above) we hardly can hope to beat it. However after implementing it in non commutative algebra (and after noting this question) I became myself curious how close we could get it for commutative algebra. The approach uses ListConvolve procedure from Carl Woll mentioned in my comment above.

First are the original data

fac[n_] := 
 Sum[Product[Unique[x], {i, 1, RandomInteger[{1, 5}]}] a[i], {i, 1, n}]
exprInput = fac[#] & /@ ConstantArray[10, 6];

Next I define very restrictive functions for conversion to/from new Association structure

toAssoc[expr_Plus] := 
 Map[ToString, 
  Association[
   MapAt[ToString, 
    Thread[((List @@ expr) /. a[_] :> 1) -> ((List @@ expr) /. 
        Thread[Rule[DeleteCases[Variables[expr], a[_]], 1]])], {All, 
     1}]]]


fromAssoc[expr_Association] := 
 Plus @@ KeyValueMap[ToExpression[StringJoin[#1, " ", #2]] &, expr]

fromAssoc[toAssoc[exprInput[[1]]]] - exprInput[[1]]
(* 0*)

modifiedInput = toAssoc /@ exprInput;

Next comes multiplication of single pair

myFastProductPairExpand[a_Association, b_Association] := 
 Merge[ListConvolve[Keys[a], Keys[b], {1, -1}, 0, 
   StringJoin[#1, " ", #2] -> StringJoin[a[#1], " ", b[#2]] &, 
   DeleteCases[List[##], 0] &], Total]

And formal expansion:

answerAssoc = 
   myFastProductPairExpand[modifiedInput[[1]], 
    myFastProductPairExpand[modifiedInput[[2]], 
     myFastProductPairExpand[modifiedInput[[3]], 
      myFastProductPairExpand[modifiedInput[[4]], 
       myFastProductPairExpand[modifiedInput[[5]], 
        modifiedInput[[6]]]]]]]; // AbsoluteTiming
(* 11.9 *)

To compare against

origInput = Times @@ exprInput;
answerExpand = origInput // Expand; // AbsoluteTiming
(* 10. 6 *)

and check

 fromAssoc[answerAssoc] - answerExpand
(*0*)

Of course, without taking into account of fromAssoc time.

Note that in ListConvolve we keep the convolution kernel as small as possible. If one reverses arguments, the computation time increases drastically. In order to get this speed I also used the fact that there are no powers of variables in the input and output. Otherwise we would need to implement some ordering of product keys, which would take additional time. Of course, variant with List instead of String in keys and values can be implemented (probably slower).

A slower version of similar procedure, which uses ordinary Times both for keys and values is given below. It is also more universal and do not require sorting of keys, since the job is done by Times itself.

toAssoc2[expr_Plus] := 
 Association[
  MapAt[Times, 
   Thread[((List @@ expr) /. a[_] :> 1) -> ((List @@ expr) /. 
       Thread[Rule[DeleteCases[Variables[expr], a[_]], 1]])], {All, 
    1}]]

fromAssoc2[expr_Association] := Plus @@ Times @@@ (Normal[expr])

modifiedInput2 = toAssoc2 /@ exprInput

myFastProductPairExpand2[a_Association, b_Association] := 
 Merge[ListConvolve[Keys[a], Keys[b], {1, -1}, 
   0, (#1*#2 -> a[#1]*b[#2]) &, DeleteCases[List[##], 0] &], Total]

answerAssoc2 = 
   myFastProductPairExpand2[modifiedInput2[[1]], 
    myFastProductPairExpand2[modifiedInput2[[2]], 
     myFastProductPairExpand2[modifiedInput2[[3]], 
      myFastProductPairExpand2[modifiedInput2[[4]], 
       myFastProductPairExpand2[modifiedInput2[[5]], 
        modifiedInput2[[6]]]]]]]; // AbsoluteTiming
(* 27.1 *)

fromAssoc2[answerAssoc2] - answerExpand
(* 0*)

So with ListConvolve + Associations one can implement very fast purely functional "expansion" procedure. It is still, however, slower than Expand. May be somebody could optimize the above idea further.