0
$\begingroup$

This is a follow-up of How to accelerate combinations and sum calculations here?

I have tried all tricks but when N reaches 35 the time is still intimidating in a 16-core server, which is far better than my Mac 2.7 GHz Intel Core i5.

I was thinking that perhaps I could store all the binomials in a table like:

N = 80;
BinomialTable = ParallelTable[Binomial[i, j], {i, 1, N}, {j, 1, N}];

Then I changed all the Binomial[x,y] into BinomialTable[[x]][[y]] in the hope that instead of calculating the combinations each time, getting them from the table would be a better idea.

The new H is named HTable with all Binomial[x,y] replaced by BinomialTable[[x]][[y]].

However, what I got is:

NmaxCompare = 6; 
A1 = 
 ParallelTable[
   H[i, j, k, l], {i, 1, NmaxCompare}, { j, 1, 
    NmaxCompare}, {k, 1, NmaxCompare}, {l, 1, NmaxCompare}] // 
  AbsoluteTiming;

AA1 = Table[
    H[i, j, k, l], {i, 1, NmaxCompare}, { j, 1, 
     NmaxCompare}, {k, 1, NmaxCompare}, {l, 1, NmaxCompare}] // 
   AbsoluteTiming;

B1 = ParallelTable[
    HTable[i, j, k, l], {i, 1, NmaxCompare}, { j, 1, 
     NmaxCompare}, {k, 1, NmaxCompare}, {l, 1, NmaxCompare}] // 
   AbsoluteTiming;

BB1 = Table[
    HTable[i, j, k, l], {i, 1, NmaxCompare}, { j, 1, 
     NmaxCompare}, {k, 1, NmaxCompare}, {l, 1, NmaxCompare}] // 
   AbsoluteTiming;

enter image description here

It seems there is too much overhead for getting the elements from the table?

Also, ParallelTable of getting elements directly (B1) cannot even beat Table of calculating the combinations each time(AA1).

Does anyone has any idea to save the time on combinations by avoiding calculating them each time? Thanks!

All the code needed here:

N = 80;
BinomialTable = ParallelTable[Binomial[i, j], {i, 1, N}, {j, 1, N}];

H[n1_, n2_, n3_, n4_] := 
(1/(4*Sqrt[n1*(n1 + 1)*n2*(n2 + 1)*n3*(n3 + 1)*n4*(n4 + 1)]))*
   Sum[(-1)^(m1 + m2 + m3 + m4)*Binomial[n1 + 1, m1 + 2]*Binomial[n2 + 1, m2 + 2]*Binomial[n3 + 1, m3 + 2]*Binomial[n4 + 1, m4 + 2]*
     Binomial[m1 + m3, m1]*Binomial[m2 + m4, m2]*(m1 + m3 + 1)*(m2 + m4 + 1)*
     (((m1 + m3 + 2)/2^(m1 + m3))*Sum[Binomial[m1 + m3 + k + 2, k]/2^k, {k, 0, m2 + m4 + 1}] + 
      ((m2 + m4 + 2)/2^(m2 + m4))*Sum[Binomial[m2 + m4 + l + 2, l]/2^l, {l, 0, m1 + m3 + 1}]), {m1, 0, n1 - 1}, {m2, 0, n2 - 1}, {m3, 0, n3 - 1}, 
    {m4, 0, n4 - 1}]

HTable[n1_, n2_, n3_, n4_] := 
(1/(4*Sqrt[n1*(n1 + 1)*n2*(n2 + 1)*n3*(n3 + 1)*n4*(n4 + 1)]))*
   Sum[(-1)^(m1 + m2 + m3 + m4)*BinomialTable[n1 + 1, m1 + 2]*BinomialTable[n2 + 1, m2 + 2]*BinomialTable[n3 + 1, m3 + 2]*BinomialTable[n4 + 1, m4 + 2]*
     BinomialTable[m1 + m3, m1]*BinomialTable[m2 + m4, m2]*(m1 + m3 + 1)*(m2 + m4 + 1)*
     (((m1 + m3 + 2)/2^(m1 + m3))*Sum[BinomialTable[m1 + m3 + k + 2, k]/2^k, {k, 0, m2 + m4 + 1}] + 
      ((m2 + m4 + 2)/2^(m2 + m4))*Sum[BinomialTable[m2 + m4 + l + 2, l]/2^l, {l, 0, m1 + m3 + 1}]), {m1, 0, n1 - 1}, {m2, 0, n2 - 1}, {m3, 0, n3 - 1}, 
    {m4, 0, n4 - 1}]
$\endgroup$
  • $\begingroup$ You expect readers to bounce to other questions to get code for a working example? No, just.... no. Please provide a complete minimal working example with all needed code. $\endgroup$ – ciao May 17 '16 at 4:15
  • $\begingroup$ "perhaps I could store all the binomials in a table" - Mathematica already caches internally. $\endgroup$ – J. M. will be back soon May 17 '16 at 4:18
  • $\begingroup$ @ciao done edited $\endgroup$ – James May 17 '16 at 4:23
  • $\begingroup$ @J.M. but it seems from input[359] and [360] getting from the table is quicker? $\endgroup$ – James May 17 '16 at 4:24
  • $\begingroup$ AIn't much of a difference, AFAICT. Try using RepeatedTiming[] as a check. $\endgroup$ – J. M. will be back soon May 17 '16 at 4:44
3
$\begingroup$

This is a very quick-and-dirty, but gives 4X speed-up (7X with tweak for symmetry) on my crappy netbook for the n=8 case, don't have time or patience to test bigger cases to see scaling differences.

Perhaps a description of what you're trying to calculate? It appears to be some combinatorial problem, there may well be a much more efficient scheme to do what you want.

In any case, give this a whirl (and if you're going to iterate cases through various n, don't clear the helper functions between runs):

Clear[s1, s2, bc, mx, sq]
Hx[i_, j_, k_, l_] /; (i > j || k > l) = Missing[];
Hx[n1_, n2_, n3_, n4_] := Hx[n3, n4, n1, n2] =
  mx[n1, n2, n3, n4]*
   Sum[(-1)^(m1 + m2 + m3 + m4)*bc[n1 + 1, m1 + 2]*bc[n2 + 1, m2 + 2]*
     bc[n3 + 1, m3 + 2]*bc[n4 + 1, m4 + 2]*bc[m1 + m3, m1]*
     bc[m2 + m4, m2]*(m1 + m3 + 1)*(m2 + m4 + 1)*
     (sq[m1, m3]*s1[m1, m2, m3, m4] + (sq[m2, m4])*
        s2[m1, m2, m3, m4]), 
      {m1, 0, n1 - 1}, {m2, 0, n2 - 1}, {m3, 0,n3 - 1}, {m4, 0, n4 - 1}];

s1[m1_, m2_, m3_, m4_] := s1[m1, m2, m3, m4] = 
   Sum[bc[m1 + m3 + k + 2, k]/2^k, {k, 0, m2 + m4 + 1}];

s2[m1_, m2_, m3_, m4_] := s2[m1, m2, m3, m4] = 
   Sum[bc[m2 + m4 + l + 2, l]/2^l, {l, 0, m1 + m3 + 1}];

bc[a_, b_] := bc[a, b] = Binomial[a, b];

sq[a_, b_] := sq[a, b] = ((a + b + 2)/2^(a + b));

SetAttributes[mx, Orderless];
mx[n1_, n2_, n3_, n4_] := mx[n1, n2, n3, n4] = 
     (1/(4*Sqrt[n1*(n1 + 1)*n2*(n2 + 1)*n3*(n3 + 1)*n4*(n4 + 1)]));

(* to make a table, e.g.: *)
nn = 8
res2 = Table[Hx[a, b, c, d], {a, nn}, {b, nn}, {c, nn}, {d, nn}]; 
$\endgroup$
  • $\begingroup$ Cool! I also find this Define a function recursively. I was wondering could we do the same tweak for symmetry like s2[m1_, m2_, m3_, m4_] := s2[m3, m2, m1, m4] like Hx? $\endgroup$ – James May 17 '16 at 15:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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