2
$\begingroup$

I would like to compile this function to make it faster:

logMultinormalDens2[x_, mean_, var2_] := Module[{},        
   t = 0.5*
     Flatten[(x - mean)].LinearSolve[var2, Flatten[(x - mean)], 
       Method -> "Cholesky"];
   (*Flatten - The LinearSolve receives only {} not {{}}*)

   det = Det[2*Pi*var2]^(0.5);

   -t - Log[det]

   ];

I've tried writing :=Compile[{},Module[(...)];]; but it takes the same time to evaluate as before. I've also tried using Block, and it didn't improve much either...

Edit: This function is called at least 300000 times (for each Time period. I want to use it several periods). Each time this distribution is called, we use a different mean and var2 (using some sum and matrix multiplication) from some other lists. Ex: mean_i = A_i+B_i, where A_i and B_i are in some other lists. var2 is a $3\times 3$ matrix

$\endgroup$
7
  • 1
    $\begingroup$ No, there is in general no point in compiling LinearSolve since its backend consists already of compiled libraries. Maybe the calling overhead is reduces a bit if you call the compiled function many times, but usually the difference is negligible. $\endgroup$ Aug 7, 2019 at 16:28
  • $\begingroup$ But tell me, how large is the matrix var2? Also, it is relevant to know how the data looks like onto which you want to apply the function. Maybe there is some structure (e.g., redundance) that can be exploited to achieve a speed-up. $\endgroup$ Aug 7, 2019 at 16:30
  • 1
    $\begingroup$ A sample usage of the function is always helpful to me, to understand how you intend to use it (and how to answer your answer). $\endgroup$ Aug 7, 2019 at 16:33
  • $\begingroup$ @HenrikSchumacher Hi Henrik, the var2 is a $3\times 3$ matrix. This function is called at least 15000 times(*20 by parallelization * number of periods). Each time this distribution is called, we use a different mean and var2 (using some sum and matrix multiplication) from some other lists. Ex: mean_i = A_i+B_i, where A_i and B_i are in some other lists $\endgroup$ Aug 7, 2019 at 16:48
  • 1
    $\begingroup$ @HenrikSchumacher I'm not sure. So far, for the several smaller simulations I've done (10 000 calls), I haven't received a single warning of ill-conditioned matrices. $\endgroup$ Aug 7, 2019 at 17:56

1 Answer 1

6
$\begingroup$

One may simply write out the matrix inverses symbolically (by Cramer's formula) and compile the resulting expression.

For dimension 3, Cramer's formula is not that inefficient. And if we know that the matrices are not too badly conditioned, that the instability (for which Cramer's formula is also infamous) are not that severe.

clogMultinormalDens2 = 
  Block[{XX, X, MM, M, ΣΣ, Σ, adjunct, det, c},
   XX = Table[Compile`GetElement[X, i], {i, 1, 3}];
   MM = Table[Compile`GetElement[M, i], {i, 1, 3}];
   ΣΣ = Table[Compile`GetElement[Σ, i, j], {i, 1, 3}, {j, 1, 3}];
   det = Det[ΣΣ];
   c = 3. Log[2. Pi];
   adjunct = Inverse[ΣΣ] det;

   With[{code = -0.5 (((XX - MM).adjunct.(XX - MM))/det + c + Log[det])}, 
    Compile[{{X, _Real, 1}, {M, _Real, 1}, {Σ, _Real, 2}},
     code,
     CompilationTarget -> "C",
     RuntimeAttributes -> {Listable},
     Parallelization -> True,
     RuntimeOptions -> "Speed"
     ]
    ]
   ];

Createing some test data. Adding a small multiple of the identity matrix ensures that the matrices are well-conditioned.

n = 100000;
x = RandomReal[{-1, 1}, {n, 3}];
μ = RandomReal[{-1, 1}, {n, 3}];
A = Plus[
   0.1 ConstantArray[IdentityMatrix[3, WorkingPrecision -> MachinePrecision], n],
   Map[#\[Transpose].# &, RandomReal[{-1, 1}, {n, 3, 3}]]
   ];

Now a test run:

a = MapThread[logMultinormalDens2, {x, μ, A}]; // AbsoluteTiming // First
b = clogMultinormalDens2[x, μ, A]; // AbsoluteTiming // First
Max[Abs[1 - b/a]]
Max[Abs[b - a]]

1.76667

0.005171

1.24345*10^-14

2.16716*10^-13

More than 300 times faster and, at least in this case, with acceptable relative and absolute errors.

$\endgroup$
11
  • 1
    $\begingroup$ Henrik, many thanks. I'm going to have to study your answer =D it's full of things I don't know. I'll get in touch later on. ;) $\endgroup$ Aug 7, 2019 at 18:02
  • $\begingroup$ You're welcome. Don't hesitate to ask me if anything is unclear. $\endgroup$ Aug 7, 2019 at 18:28
  • $\begingroup$ Hi Henrik, I've just had some time, and I have some questions. 1 - What does Compile`GetElement do exactly? I've tried searching through the mathematica documentation, but neither ' nor GetElement are in it. (to be continued) $\endgroup$ Aug 9, 2019 at 8:27
  • $\begingroup$ 2 - From a previous question of mine, we(it was also your answer, incidentally =D mathematica.stackexchange.com/questions/191421/…) concluded that if I had to use a different x, mean and var2, we should use something similar to MapThread. How would we do that for this answer? This should increase the expected computing time for your procedure... $\endgroup$ Aug 9, 2019 at 8:33
  • 1
    $\begingroup$ Sorry for the delay. Many thanks Henrik! +1 ;) $\endgroup$ Sep 15, 2019 at 18:02

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.