6
$\begingroup$

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}
*)
$\endgroup$
8
  • $\begingroup$ If you could give a minimal working example, I think you'd get more responses (other than responses to close the question). $\endgroup$
    – JimB
    Dec 24, 2019 at 18:05
  • 3
    $\begingroup$ Since the calculation is not very complicated, I don't think there's much we can do by using Mathematica. I suggest you to use FORM in this case. $\endgroup$
    – Wen Chern
    Dec 24, 2019 at 18:50
  • 1
    $\begingroup$ @WenChen: Didn't know about FORM. Thanks. nikhef.nl/~form $\endgroup$
    – JimB
    Dec 24, 2019 at 19:32
  • $\begingroup$ Thanks, I'll give FORM a shot! I wait a day though and mark the question as closed if no miracle method shows up until tomorrow. $\endgroup$
    – dzsoga
    Dec 24, 2019 at 20:11
  • 2
    $\begingroup$ Your might be interesting in the approach mathematica.stackexchange.com/questions/148230/… and Carl Woll in particular. My suggestions: 1) instead of using "Plus" as container to collect different terms use Association. 2) Split a product of n terms into pairwise products. 3) use the above method of Carl Woll to Keys of association. I implemented all this to Clifford product of multivectors. The speed gain was 100x when compared to Plus and implementation with nonpaired Expand. $\endgroup$
    – Acus
    Dec 27, 2019 at 8:46

2 Answers 2

6
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ Thank you, very insightful answer! I have found a mathematica package, using FORM. Unfortunately nothing ever is that easy, I was not able to install FORM on my win 10 OS yet, but working on it. $\endgroup$
    – dzsoga
    Dec 27, 2019 at 17:20
  • $\begingroup$ I updated my answer to show how this expansion can be done using FORM. $\endgroup$
    – vsht
    Dec 27, 2019 at 18:39
  • $\begingroup$ Thanks! I looked into it and the problem is, that I have a 64 bit OS and and the FORM installer executable is 16 bit. There is an emulator for such things but only for 32 bit OS :( If you are interested here is the doc for the mathematica package for FORM: arxiv.org/abs/1610.09331 $\endgroup$
    – dzsoga
    Dec 27, 2019 at 18:46
  • 1
    $\begingroup$ Hmm, I also tried to compile FORM on W10 using Cygwin and it seems that it doesn't work out of the box, although I remember having succeeded with W7 few years ago. I guess the reason is that in our community people mostly use Linux or macOS, so nobody really cares if a code compiles on Windows. Perhaps it is worth creating a bug report here github.com/vermaseren/form/issues. $\endgroup$
    – vsht
    Dec 27, 2019 at 21:27
  • 1
    $\begingroup$ @dzsoga In fact, one can compile a 64-bit FORM on a 64-bit version of Windows 10 without much efforts, if one uses the right tools. See here: github.com/vermaseren/form/issues/342 $\endgroup$
    – vsht
    Jan 14, 2020 at 12:53
3
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.