3
$\begingroup$

I need to maximize a function numerically, thus the speed of this function is very important. Here is the most time consuming part of this function.

pPoisson[lambda_, mu_, x_, 
  y_] := (E^(-lambda - mu)*lambda^x*mu^y)/(x!*y!)
list[a1_, b1_, a2_, b2_, c_, upper_ : 15] := 
 Module[{lambda = Exp[a1 - b2 + c], mu = Exp[a2 - b1], m},
  m = Table[pPoisson[lambda, mu, i, j], {j, 0, upper}, {i, 0, upper}];
  {Total[UpperTriangularize[m, 1], Infinity], Total[Diagonal[m]], 
   Total[LowerTriangularize[m, -1], Infinity]}]
Do[list[1.3, 0.6, 0.2, 0.2, 0.17, 15], {i, 1000}]; // AbsoluteTiming
(* {1.06158, Null} *)

The pPoisson is the product of 2 Poisson Distribution probability mass function. The function list generates a matrix of pPoisson and calculates the sum of the matrix upper, diagnal and lower part. Running list 1000 times takes 1.06 s. This is not the speed I would like to have. Therefore, I tried to compile it, but the compiled version seems to be 2x slower.

clist = Compile[{{a1, _Real}, {b1, _Real}, {a2, _Real}, {b2, _Real}, \
{c, _Real}, {upper, _Integer}},
  Module[{lambda = Exp[a1 - b2 + c], mu = Exp[a2 - b1], i, j},
   {Sum[(Exp[(-lambda - mu)]*lambda^i*mu^j)/(Product[x, {x, 1, i}]*
        Product[x, {x, 1, j}]), {i, 1, upper}, {j, 0, i - 1}],
    Sum[(Exp[(-lambda - mu)]*lambda^i*
        mu^i)/(Product[x, {x, 1, i}]^2), {i, 0, upper}],
    Sum[(Exp[(-lambda - mu)]*lambda^i*mu^j)/(Product[x, {x, 1, i}]*
        Product[x, {x, 1, j}]), {j, 1, upper}, {i, 0, j - 1}]}],
  CompilationTarget -> "C",
  Parallelization -> True, RuntimeOptions -> "Speed"]
Do[clist[1.3, 0.6, 0.2, 0.2, 0.17, 15], {i, 1000}]; // AbsoluteTiming
(* {1.9038489`, Null} *)

Any idea why and how to accelarate it?

$\endgroup$
2
  • $\begingroup$ By using immediate assignments (=) instead of delayed assignments (:=) you may speed up your code quite a bit even without compilation. wolfram.com/language/elementary-introduction/2nd-ed/… $\endgroup$
    – Roman
    Commented Jul 26, 2020 at 11:02
  • 2
    $\begingroup$ Look up CompilePrint, MainEvaluate, and other pointers here and here $\endgroup$
    – Michael E2
    Commented Jul 26, 2020 at 16:06

1 Answer 1

3
$\begingroup$

The slowness is due to several instances of MainEvaluate. I replaced the Product factorials with Gamma, which turns out to be compilable.

clist = Compile[{{a1, _Real}, {b1, _Real}, {a2, _Real}, {b2, _Real}, 
   {c, _Real}, {upper, _Integer}},
  Module[{
    lambda = Exp[a1 - b2 + c],
    mu = Exp[a2 - b1],
    i, j},
   {Sum[(Exp[(-lambda - mu)]*lambda^i*mu^j)/(Gamma[1 + i]*
        Gamma[1 + j]), {i, 1, upper}, {j, 0, i - 1}]
    , Sum[(Exp[(-lambda - mu)]*lambda^i*mu^i)/(Gamma[1 + i]^2), {i, 0,
       upper}]
    , Sum[(Exp[(-lambda - mu)]*lambda^i*mu^j)/(Gamma[1 + i]*
        Gamma[1 + j]), {j, 1, upper}, {i, 0, j - 1}]
    }]
  , CompilationTarget -> "C"
  (*,Parallelization\[Rule]True*) (* useless without Listable *)
  , RuntimeOptions -> "Speed"]

Do[clist[1.3, 0.6, 0.2, 0.2, 0.17, 15], {i, 1000}]; // RepeatedTiming
(*  {0.0067, Null}  *)

The OP's clist took about 2 sec.

$\endgroup$

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.