Can one Map a simplifying function to the summands of an Expand before fully expanding the expression in question first? Consider the following example:

I want to integrate the following product over [0,1] in the variable t.

imax = 30;
p1 = Product[a t^i + b t^(imax - i) , {i, imax}];

For that I define a custom integration (much faster than the general Integrate)

CustomIntegration[exp_] := exp /. t^k_. -> 1/(k + 1)

for which I have to use Expand first. The problem is that the expansion really eats a ton of memory, even though the final result is actually much smaller.

ByteCount[p2 = Expand@p1]
ByteCount[CustomIntegration /@ p2]



So, the expansion requires around 100-times as much memory as needed. My actual script currently needs around 70 GB of memory. So, how could I reduce that?


2 Answers 2


In this specific case I believe that this will do the same operation in one command:

Collect[p1, a, CustomIntegration] // ByteCount
Collect[p1, a, CustomIntegration] === (CustomIntegration /@ Expand[p1])

However, this apparent success is false, as one can see using MaxMemoryUsed that much more memory is allocated:

(* in a fresh kernel *)

CustomIntegration[exp_] := exp /. t^k_. -> 1/(k + 1)
imax = 50;
p1 = Product[a t^i + b t^(imax - i), {i, imax}];

Collect[p1, a, CustomIntegration] // ByteCount


(Session baseline in v7 is 15MB.)


I can't find a way to do it with both low memory and reasonable time consumption, but here is what I did when I encountered a similar problem. Basically you need to do some math for MMA.

The math part

Your p1 can be expanded by the so called multi-binomial theorem:

$$\begin{split} p_1=&\prod_{k=1}^m (a\,t^k + b\,t^{m-k} )\\ =&\sum_{j_1=0}^1 \sum_{j_2=0}^1 \cdots \sum_{j_m=0}^1 \left(\prod_{k=1}^m \binom{1}{j_k}a^{j_k} b^{1-j_k}\right)\left(\prod_{k=1}^m t^{j_k k + (m-k)(1-j_k)} \right) \end{split}$$

So we can see that there is an one-to-one map between $\{0,1,\dots,2^m\}$ and the set of all possible $\boldsymbol{j}:=(j_1,j_2,\dots,j_m)$, and we can program a loop on it. But let's go on for the moment.

Let $\nu_{\boldsymbol{j}}:=\sum_{k=1}^m j_k$, $\mu_{\boldsymbol{j}}:=\sum_{k=1}^m j_k k$,

$$ p_1=\sum_{\forall \boldsymbol{j}} a^{\nu_{\boldsymbol{j}}}b^{m-\nu_{\boldsymbol{j}}}\,t^{2\mu_{\boldsymbol{j}}-m\nu_{\boldsymbol{j}}+(m-1)m/2} $$

Now we can substitute your integration rule into $p_1$, so CustomIntegration /@ p2 should be

$$\begin{split} I=&\sum_{\forall \boldsymbol{j}}\frac{a^{\nu_{\boldsymbol{j}}}b^{m-\nu_{\boldsymbol{j}}}}{2\mu_{\boldsymbol{j}}-m\nu_{\boldsymbol{j}}+(m-1)m/2+1}\\ =&\sum_{\nu=0}^m a^\nu b^{m-\nu} \sum_{\forall \boldsymbol{j}(\nu)}\frac{1}{2\mu_{\boldsymbol{j}(\nu)}-m\nu+(m-1)m/2+1} \end{split}$$

where $\mu_{\boldsymbol{j}(\nu)}$ is the $\mu_{\boldsymbol{j}}$ given $\nu$ out of $m$ elements of $\boldsymbol{j}$ equal to $1$ while others equal to $0$.

So what we need to do is to compute the coefficients of $a^\nu b^{m-\nu}$. And that is the place I failed.

The program part

Here is what I tried.

First we setup some monitors:

miu = {MemoryInUse[]};
          AppendTo[miu, MemoryInUse[]];
          If[Length[miu] > 3000, miu = miu[[-3000 ;; -1]]];
          ListLinePlot[miu, PlotRange -> All],
        UpdateInterval -> 0]]

Dynamic[ListLinePlot[{target, sum}, PlotRange -> All]]

Then the result by your original method:

imax = 20; 
p1 = Product[a*t^i + b*t^(imax - i),  {i, imax}]; 
target = p1 // ExpandAll // CustomIntegration //
                  CoefficientList[#, {a}] & // # /. b -> 1 &;

Now my VERY slow (but memory-economical) way:

sum = ConstantArray[0, imax + 1];
Module[{m = imax, temp, base, μ},
       base = Range[m];
          temp = 1 + (m (m - 1))/2 - m ν;
             μ = Total[Subsets[base, {ν}, {c}][[1]]];
             sum[[ν + 1]] += 1/(temp + 2 μ),
           {c, Binomial[m, ν]}],
        {ν, 0, m}]

sum == target

For imax$>20$ the time consumption will be unbearable. I believe a proper parallelization and compilation will make it a bit faster, but the gain might still be limited. Also, with a naive compilation, you loose the infinite precision.


Your Answer

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

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