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$
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.

| improve this answer | |
$\endgroup$

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.