Make code faster and in C-compatible manner?

Questions
, , , ,
1

Hello,

I need to make the following code very fast:

B[a_, m_, t_] := Module[{n, p, q, seq, s, f, b},
   n = 1 - m; p = m/a; q = m (1 - 1/a);
   seq[i_, j_] := 
    Sequence @@ 
     Flatten[Table[{i == K || j == t - K, s[i, j, K]}, {K, 1, 
        Round[t/2]}], 1];
   s[i_, j_, K_] := 
    Sum[Binomial[i, i - k] Binomial[t - i, j - i + k] (p q)^
       k ((n + p) (n + q))^(K - k), {k, 0, K}];
   f[i_, j_] := 
    If[(t - i - j) >= 0, (n + p)^(t - i - j), (n + q)^-(t - i - j)];
   b[i_, j_] := Which[
     j == t, p^(j - i) f[i, j] Binomial[t, i],
     i == t, q^(i - j) f[i, j] Binomial[t, j],
     i <= j, 
     p^(j - i) f[i, j] Binomial[t, i] Which[Evaluate[seq[i, j]]],
     i > j, 
     q^(i - j) f[i, j] Binomial[t, j] Which[Evaluate[seq[j, i]]]];
   Return[Table[b[i, j], {i, 1, t}, {j, 1, t}]]];

It generates a t x t numerical matrix for an integer a and a real number m between 0 and 1. For low t it is fast enough but I need it for t=2000. I estimated that this would take years on a desktop computer… So I guess I need to compile the whole thing to C-code and possibly parallelize it. The problem is the seq function. It says that it cannot be compiled because the expression “Sequence” is not allowed. It generates a sequence of conditions of length t/2 and also uses the function s inside it. How can it be written in C-compatible manner? And is it possible to compile a function that uses another inside it? So I have a lot of questions. Unfortunately, my programming skills are only basic so any help would greatly be appreciated.

2

In case anyone is interested: the main brake in the above code is the seq function. It creates useless redundancy because the list of conditions it generates can be integrated directly into Do-loops as follows:

B[a_, m_, t_] := Module[{n, p, q, s, f, b},
   n = 1 - m; p = m/a; q = m (1 - 1/a);
   s[i_, j_, K_] := 
    Sum[Binomial[i, i - k] Binomial[t - i, j - i + k] (p q)^
       k ((n + p) (n + q))^(K - k), {k, 0, K}];
   f[i_, j_] := 
    If[(t - i - j) >= 0, (n + p)^(t - i - j), (n + q)^-(t - i - j)];
   Do[b[i, j] = 
     p^(j - i) f[i, j] Binomial[t, i]/Binomial[t, j] s[i, j, i], {i, 
     1, Round[t/2] + 1}, {j, i, t - i}];
   Do[b[i, j] = 
     p^(j - i)
       f[i, j] Binomial[t, i]/Binomial[t, j] s[i, j, t - j], {j, 
     Round[t/2], t - 1}, {i, t - j + 1, j}];
   Do[b[i, t] = p^(t - i) f[i, t] Binomial[t, i], {i, 1, t}];
   Do[b[j, i] = (q/p)^(j - i)
       b[i, j] Binomial[t, j]/Binomial[t, i], {i, 1, t}, {j, i + 1, 
     t}];
   Return[Table[b[i, j], {i, 1, t}, {j, 1, t}]]];

This way it is possible to reduce the time for t=2000 from thousands of years to hours, still non-compiled…

3

After 8 years, I am again struggling with my matrix B to get it still faster… The problem is that the code involves big numbers generated by Binomial[t, i]. For instance, Binomial[t, t/2]=2.04815162699310^600 for t=2000. So this is well above double floating of 2^1023=8.9884710^307. In the code, there are also factorials t!, Anyone has an idea how such numbers can be handled to prevent the compiled code to generate errors and reverse back to the slow uncompiled code?

4

Maybe post your question here. They will definitely answer your question and you will be satisfied

I’m not a good programmer who knows tricks to speed up the code and I won’t help you in your situation.

Regards M.I.

5

I rewrote a bit. This still uses seq but not as a Sequence object. It has a final clause that is never actually needed.

B2[a_, m_, t_] := Module[{n, p, q, seq, s, f, b},
  n = 1 - m;
  p = m/a;
  q = m - p;
  s[i_, j_, K_] := 
   Sum[Binomial[i, i - k]  Binomial[t - i, j - i + k]  (p  q)^
      k  ((n + p)  (n + q))^(K - k), {k, 0, K}];
  f[i_, j_] := 
   If[(t - i - j) >= 0, (n + p)^(t - i - j), (n + q)^-(t - i - j)];
  seq[i_, j_] := Which[
    i <= Round[t/2] && t - j <= Round[t/2],
    s[i, j, Min[i, t - j]]
    ,
    i <= Round[t/2],
    s[i, j, i]
    ,
    t - j <= Round[t/2],
    s[i, j, t - j]
    ,
    True, 0];
  b[i_, j_] := Which[
    j == t,
    p^(j - i)  f[i, j]  Binomial[t, i]
    ,
    i == t,
    q^(i - j)  f[i, j]  Binomial[t, j]
    ,
    i <= j,
    p^(j - i)  f[i, j]  Binomial[t, i]  seq[i, j]
    ,
    i > j,
    q^(i - j)  f[i, j]  Binomial[t, j] seq[j, i]];
  Table[b[i, j], {i, 1, t}, {j, 1, t}]]

As you note, this can get into trouble with machine numbers.

Timing[B2[3, .22, 400];]

(* During evaluation of In[127]:= General::munfl: 3.75211*10^-299 5.02373*10^-10 is too small to represent as a normalized machine number; precision may be lost.

During evaluation of In[127]:= General::munfl: 2.75154*10^-300 5.88718*10^-10 is too small to represent as a normalized machine number; precision may be lost.

During evaluation of In[127]:= General::munfl: 2.0178*10^-301 6.89904*10^-10 is too small to represent as a normalized machine number; precision may be lost.

During evaluation of In[127]:= General::stop: Further output of General::munfl will be suppressed during this calculation.

Out[127]= {30.3564, Null} *)

It handles higher precision without signs of distress though, and this is not hugely slower.

Timing[b2big = B2[3, .22`20, 400];]

(* Out[129]= {51.2323, Null} *)
6

If you change Sum to: Total@Table code is 20% faster.

7

Thank you Daniel and Mariusz for your suggestions. Meanwhile I found that the following code is the fastest version:

  Bmatrix[a_, t_, m_] := Block[{n, p, q, B},
   n = 1 - m; p = m/a; q = m (1 - 1/a); B = ConstantArray[0., {t, t}];
   Do[B[[i, j]] = 
     j! /i! p^(-i + j) (n + p)^(-j + t) (n + q)^
      i Hypergeometric2F1Regularized[-i, j - t, 1 - i + j, (
       p q)/((n + p) (n + q))], {i, 1, Floor[t/2]}, {j, i, t - i}];
   Do[B[[i, j]] = 
     j! /i! p^(-i + j) (n + p)^(-j + t) (n + q)^
      i Hypergeometric2F1Regularized[-i, j - t, 1 - i + j, (
       p q)/((n + p) (n + q))], {j, Ceiling[t/2], t - 1}, {i, 
     t - j + 1, j}];
   Do[B[[i, t]] = p^(t - i) (n + q)^i Binomial[t, i], {i, 1, t}];
   Do[B[[i, j]] = (q/p)^(i - j) Binomial[t, i]/
      Binomial[t, j] B[[j, i]], {i, 2, t}, {j, 1, i - 1}];
   Return[B]];

It can be C-compiled but j! /i! and similar expressions produce too big numbers. To prevent this, I could use

((p^j) j! )/(p^i i!) = FactorialPower[p j, j - i, p]

in the first Do loop because there j>=i and 0<p<<1. But this doesn’t help because strangely the thus compiled code is slower than the uncompiled one…

Edit: I just found out that Product[p (j - k), {k, 0, j - i - 1}] is much faster than FactorialPower[p j, j - i, p] if compiled. Shouldn’t built-in functions be more rapid?

8

You might experiment with taking logs, PowerExpanding using LogGamma instead of Log[Factorial[...]], and exponentiating at the end. If there is intermediate coefficient growth followed by cancelation, the LogGamma usage should ameliorate the issue numerically.