Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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@p1
ByteCount[p2 = Expand@p1]
ByteCount[CustomIntegration /@ p2]

gives

9376
1404216
13696

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?

share|improve this question

2 Answers 2

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

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

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
MaxMemoryUsed[]
29888

28313200

(Session baseline in v7 is 15MB.)

share|improve this answer

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[]};
Dynamic[Refresh[
          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];
       Do[
          temp = 1 + (m (m - 1))/2 - m ν;
          Do[
             μ = Total[Subsets[base, {ν}, {c}][[1]]];
             sum[[ν + 1]] += 1/(temp + 2 μ),
           {c, Binomial[m, ν]}],
        {ν, 0, m}]
      ];

sum == target
True

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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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