Comparison with other packages
This benchmark was inspired by this post. The full benchmark code is benchmark/benchmark_crosslang.jl, with the figures drawn by benchmark/plot_crosslang.jl. The original univariate-only version is kept as a Jupyter notebook.
Scope and limitations of competing packages
A key distinction of ExpectationMaximization.jl is its genericity: it works with any mixture of distributions supported by Distributions.jl (univariate, multivariate, continuous, discrete, or custom), without any modification to the core algorithm. The competing packages benchmarked here are, in contrast, largely restricted to Gaussian mixtures:
| Package | Language | Gaussian only? | Notes |
|---|---|---|---|
Sklearn | Python | Yes | Hardcoded Gaussian[2] |
mixtools | R | Mostly | Supports some other families but not extensible |
mixem | Python | Mostly | Numerically fragile[3] |
GaussianMixtures.jl | Julia | Yes | Hardcoded Gaussian |
ExpectationMaximization.jl | Julia | No | Any Distributions.jl distribution. Possibilities for custom distributions via the Distributions.jl interface. |
The benchmark below only tests the Gaussian mixture case (the most common), which is deliberately the strongest case for the specialized packages. Despite this, ExpectationMaximization.jl remains highly competitive.
Why is ExpectationMaximization.jl fast?
No heavy programming tricks are used. The performance comes from standard Julia best practices:
- E-step: memory allocated once and reused at every iteration (the posteriors overwrite the log-likelihood matrix in place),
@views, type-stable code, and a fused allocation-free row-wise log-sum-exp/softmax. - M-step: delegates to
fit_mlefromDistributions.jl, which is well-optimized for each distribution (e.g., see the Multivariate Normal implementation).
Many more optimizations are possible, however, I'd like to keep the code as simple and readable as possible. Note that as of v0.3.5, I am testing LLM suggestions to improve performance without sacrificing readability (too much). If you have suggestions, please open an issue or PR.
Results
The published run is the one produced by the Benchmark (Julia + R + Python) GitHub Actions workflow, on an ubuntu-latest runner (AMD EPYC 7763) with a single Julia thread and a single BLAS thread. A CI runner is slower and noisier than a workstation, but it is the same machine for all four backends and anybody can re-run it, which matters more here than absolute speed.
All benchmark cases are shown in a single figure: one panel per (K, D), with the ratio of each backend's fit time to ExpectationMaximization.jl's against the sample size N. Above the dashed line means slower than ExpectationMaximization.jl; the vertical axis is logarithmic because mixtools is several orders of magnitude slower on the multivariate cases.
Only the ratio is shown. Absolute times are not comparable across cases — each case fixes its own number of EM iterations — whereas a ratio between two backends at the same N of the same case always compares equal amounts of work.
Every backend gets the same data, the same initial conditions, the same full-covariance model and the same, fixed number of EM iterations, and the benchmark verifies all of that as it runs rather than assuming it: after each timed fit it checks that the backends agree on the fitted weights, means and variances, and that each one really performed the number of iterations it was asked for. The full transcript of the published run — cross-check results, iteration counts, per-case totals — is kept alongside the figures:
benchmark/results/benchmark_log_latest.txt
Three caveats that transcript will show you.
mixtoolsperforms a multicycle ECM step — two E-steps per iteration, the second one after the means have been updated — so it does more work per iteration than the other three, and its intermediate iterates differ from theirs. It is therefore compared with them at its fixed point, with everything run to convergence. Part of its distance from the other backends is this extra work, not R being slow.GaussianMixtures.jlends its loop on an M-step, while the other three evaluate the likelihood of the parameters they finish on and so pay one further E-step. At the same iteration count it therefore does about 6% less work atiters = 8and 2.5% atiters = 20. Asking it for one more iteration would buy an extra M-step too, so the difference is reported rather than papered over.- A backend whose fit exceeds 100 s is dropped from the larger
Nof that case, so its curve stops early rather than costing hours. On the published run this happened tomixtoolsin both multivariate cases.
On univariate Gaussian mixtures, ExpectationMaximization.jl is 3-5× faster than mixtools (R) and Sklearn (Python) — up to 7× and 16× respectively at N = 500, where those two still pay a fixed per-call overhead. GaussianMixtures.jl is the one backend ahead of it here, by 1.6-2× over most of the range; the next section says why. On multivariate Gaussian mixtures it leads everything tested: 1.7-3.4× faster than GaussianMixtures.jl, 2.6-19× faster than Sklearn, and three orders of magnitude faster than mixtools, whose mvnormalmixEM is the outlier of the whole benchmark. And unlike every other package here, it does all this while accepting an arbitrary Distributions.jl mixture.
Reproducing these figures
julia --threads=1 --project=benchmark benchmark/benchmark_crosslang.jl # writes benchmark/results/*
julia --project=benchmark benchmark/plot_crosslang.jl # writes benchmark/timing_crosslang*.svgThe benchmark and the plotting are separate scripts. The benchmark writes, under benchmark/results/, both a dated copy and a _latest one of:
| file | contents |
|---|---|
benchmark_timings_<date>.csv | one row per (case, backend, N): case,backend,K,D,N,time_s,iters_asked,iters_run |
benchmark_log_<date>.txt | the run transcript — the file linked above |
system_info_<date>.txt | Julia, R and Python versions, package versions, thread counts |
Both are rewritten after every case, so an interrupted run still leaves the cases that did finish.
If you have comments to improve these benchmarks, they are welcome.
Cross-language comparisons are inherently imperfect[1]. PythonCall.jl and RCall.jl introduce a small overhead (~few milliseconds), which was verified to be negligible here.
There is also a Julia-only suite in benchmark/benchmarks.jl, run on every pull request by AirspeedVelocity.jl, which covers univariate, multivariate, weighted and Classical and Stochastic EM. It allows tracking performance regressions and improvements over time. These results are displayed in each PR as a comment.
- 1
@btimewithRCall.jlandPythonCall.jlmay add a small overhead; see this discussion. Timings were cross-checked againstRmicrobenchmarkand Pythontimeit, which gave consistent results.BenchmarkTools.jlautomatically determines the number of repetitions needed for a reliable estimate. - 2
Sklearn'sGaussianMixtureused to run K-means initialization even when initial conditions are explicitly provided — see this discussion, issue, and PR. I should be fix by now. - 3
mixemoverflows for $n \gtrsim 500$ due to a fragilelogsumexpimplementation and was excluded from the benchmark.