# Can we compile this function/distribution?

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

• 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. – Henrik Schumacher Aug 7 '19 at 16:28
• 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. – Henrik Schumacher Aug 7 '19 at 16:30
• 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). – Arnoud Buzing Aug 7 '19 at 16:33
• @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 – An old man in the sea. Aug 7 '19 at 16:48
• @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. – An old man in the sea. Aug 7 '19 at 17:56

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[CompileGetElement[X, i], {i, 1, 3}];
MM = Table[CompileGetElement[M, i], {i, 1, 3}];
ΣΣ = Table[CompileGetElement[Σ, i, j], {i, 1, 3}, {j, 1, 3}];
det = Det[ΣΣ];
c = 3. Log[2. Pi];

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.

• 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. ;) – An old man in the sea. Aug 7 '19 at 18:02
• You're welcome. Don't hesitate to ask me if anything is unclear. – Henrik Schumacher Aug 7 '19 at 18:28
• Hi Henrik, I've just had some time, and I have some questions. 1 - What does CompileGetElement do exactly? I've tried searching through the mathematica documentation, but neither ' nor GetElement are in it. (to be continued) – An old man in the sea. Aug 9 '19 at 8:27
• 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... – An old man in the sea. Aug 9 '19 at 8:33
• Sorry for the delay. Many thanks Henrik! +1 ;) – An old man in the sea. Sep 15 '19 at 18:02