Speedup of KL Divergence by Monte Carlo Simulation

I started with Mathematica about just a week, so I don't know how to take all of this environment. What I want to know is if there is a way to speedup my code to compute the KL-Divergence of Gaussian Mixtures using Monte Carlo Simulation. I believe I did the correct implementation, but I was hoping it could achieve a better performance. The distributions values are just for testing purposes.

weights1 = {.3, .7};
dist11 = MultinormalDistribution[{0., 0.}, {{1., 0.}, {0., 1.}}];
dist12 = MultinormalDistribution[{3., 3.}, {{1., 0.}, {0., 1.}}];
gmm1 = MixtureDistribution[weights1, {dist11, dist12}];

weights2 = {.2, .8};
dist21 = MultinormalDistribution[{4., 4.}, {{2., 0.}, {0., 1.}}];
dist22 = MultinormalDistribution[{6., 2.}, {{1., 0.}, {0., 1.}}];
gmm2 = MixtureDistribution[weights2, {dist21, dist22}];

samplecount = 20000;
s = RandomVariate[gmm1, samplecount];

f1[x_] := f1[x] = PDF[gmm1, x];
f2[x_] := f2[x] = PDF[gmm2, x];
p1 = ParallelMap[f1, s];
p2 = ParallelMap[f2, s];
f3 = Compile[{x, y}, x * Log[2., x / y],
RuntimeAttributes -> {Listable}];
kl = ParallelSum[f3[p1[[i]], p2[[i]]], {i, 1, samplecount}] /
samplecount;
• A quick check shows that the rate-limiting lines are the p1 = ParallelMap[f1, s]; p2 = ParallelMap[f2, s]; lines. The ParallelSum line executes almost instantly. Sep 27 '14 at 2:00
• Since f1 and f2 compute PDF[gmm, x], your question is essentially equivalent to the shorter question, "What is the fastest way to evaluate PDF of a 2D distribution at a bunch of points?" Sep 27 '14 at 2:06
• However, unless I'm misreading your question, it looks like your 2D distribution has a PDF which is the sum of several multinormal Gaussian peaks. This means that it can be computed in closed-form manually, without resorting to PDF or distribution objects. Since Gaussians are compilable to "C" code, you should be able to get a speedup factor of at least a thousand. So you should be able to get much better performance. Sep 27 '14 at 2:12
• I was able to speed up your code by a factor of over 4000. I'll post a partial example in a while, but I'll have to wait until morning to polish it up as I have to go to bed soon. Sep 27 '14 at 2:29

You can make the performance-limiting parts of your code run over 4000 times as fast. The original code takes 7.71 seconds on my machine when samplecount = 2000. By correctly compiling the code, the same operation takes a few milliseconds.

First off, note that the line

p1 = ParallelMap[f1, s];

is essentially computing PDF[gmm1, x] at a large number of points in 2D space.

Second, note that PDF[gmm1, x] is the PDF associated with the weighted sum of two multinormal distributions. This means that it is expressible in closed form. Since it is expressible in closed form, it is compilable, which also means that it can be efficiently executed on large lists of datapoints. So it is natural to expect a very large speedup.

So let's derive the closed-form expression.

dist1 = MultinormalDistribution[{a, b}, {{c, d}, {e, f}}];
dist2 = MultinormalDistribution[{g, h}, {{i, j}, {k, l}}];
dist3 = MixtureDistribution[{m, n}, {dist1, dist2}];
PDF[dist3, {x, y}]

yields the PDF

$$\frac{m \exp \left(\frac{1}{2} \left(-\frac{(y-b) (-a e+b c-c y+e x)}{d e-c f}-\frac{(x-a) (-a f+b d-d y+f x)}{c f-d e}\right)\right)}{2 \pi (m+n) \sqrt{c f-d e}}+\frac{n \exp \left(\frac{1}{2} \left(-\frac{(y-h) (-g k+h i-i y+k x)}{j k-i l}-\frac{(x-g) (-g l+h j-j y+l x)}{i l-j k}\right)\right)}{2 \pi (m+n) \sqrt{i l-j k}}$$

which we can then evaluate using the coefficients used in gmm1:

(E^(1/2 (-(((-b + y) (b c - a e + e x - c y))/(
d e - c f)) - ((-a + x) (b d - a f + f x - d y))/(-d e + c f)))
m)/(2 Sqrt[-d e + c f] (m + n) \[Pi]) + (
E^(1/2 (-(((-h + y) (h i - g k + k x - i y))/(
j k - i l)) - ((-g + x) (h j - g l + l x - j y))/(-j k + i l)))
n)/(2 Sqrt[-j k + i l] (m + n) \[Pi]) /.
Thread[{a, b, c, d, e, f, g, h, i, j, k, l, m, n} -> {0, 0, 1, 0, 0,
1, 3, 3, 1, 0, 0, 1, 3/10, 7/10}]

resulting in

(7 E^(1/2 (-(-3 + x)^2 + (3 - y) (-3 + y))))/(20 \[Pi]) + (
3 E^(1/2 (-x^2 - y^2)))/(20 \[Pi])

which we can then compile as

compiledF1 =
Compile[{{x, _Real}, {y, _Real}}, (
7 E^(1/2 (-(-3 + x)^2 + (3 - y) (-3 + y))))/(20 \[Pi]) + (
3 E^(1/2 (-x^2 - y^2)))/(20 \[Pi]),
RuntimeAttributes -> {Listable}, CompilationTarget -> "C"];

We can then run the equivalent to the ParallelMap operation I highlighted in the beginning of the answer in the following way:

First@AbsoluteTiming[compiledF1 @@ (s\[Transpose])]

which executes in 1 millisecond on my laptop, a considerable speedup.

Likewise, you can verify the correctness of the computation by plotting both the original PDF and the compiled version:

Plot3D[PDF[gmm1, {x, y}], {x, -3, 6}, {y, -3, 6}, PlotRange -> All]
Plot3D[compiledF1[x, y], {x, -3, 6}, {y, -3, 6}, PlotRange -> All]

resulting in  which both look the same, confirming that the compiled version and the Distribution PDFs are in fact the same.

This naturally can be extended to the gmm2 distribution and the related ParallelMap computation, hopefully you can work out the rest of the details yourself.