In this example, we compare the runtime and accuracy of the EM algorithm, the mix-SQP algorithm, and the interior-point method implemented by the MOSEK commercial solver (and called via the KWDual function in the R package REBayes).
Before attempting to run this Julia code, make sure your computer is properly set up to run this code by following the setup instructions in the README of the git repository.
We begin by loading the Distributions, LowRankApprox, RCall and other packages, as well as some function definitions used in the code chunks below.
using Pkg
using Printf
using Random
using Distributions
using LinearAlgebra
using SparseArrays
using LowRankApprox
using RCall
include("../code/datasim.jl");
include("../code/likelihood.jl");
include("../code/mixEM.jl");
include("../code/mixSQP.jl");
include("../code/REBayes.jl");
Next, initialize the sequence of pseudorandom numbers.
Random.seed!(1);
Let's begin with a smaller example with 50,000 samples.
z = normtmixdatasim(round(Int,5e4));
Compute the $n \times k$ likelihood matrix for a mixture of zero-centered normals, with $k = 20$. Note that the rows of the likelihood matrix are normalized by default.
sd = autoselectmixsd(z,nv = 20);
L = normlikmatrix(z,sd = sd);
size(L)
First we run each of the optimization algorithms once to precompile the relevant functions.
outem = mixEM(L,maxiter = 10);
outsqp = mixSQP(L,maxiter = 10,verbose = false);
outip = REBayes(L);
Next, let's fit the model using the three algorithms.
@time xem, tem = mixEM(L,tol = 1e-4,maxiter = 1000);
@time xip, tip = REBayes(L);
@time outsqp = mixSQP(L,verbose = false);
The mix-SQP algorithm algorithm is much faster than the other two methods, with EM being the slowest.
The quality of the IP and SQP solutions is very similar, whereas the EM solution is much worse:
fem = mixobjective(L,xem);
fip = mixobjective(L,xip);
fsqp = mixobjective(L,outsqp["x"]);
fbest = minimum([fem fip fsqp]);
@printf "Difference between EM solution and best solution: %0.2e\n" fem - fbest
@printf "Difference between IP solution and best solution: %0.2e\n" fip - fbest
@printf "Difference between SQP solution and best solution: %0.2e\n" fsqp - fbest
Next, let's see what happens when we apply these three algorithms to a larger data set.
z = normtmixdatasim(round(Int,1e5));
As before, we compute the $n \times k$ likelihood matrix for a mixture of zero-centered normals. This time, we use a finer grid of $k = 100$ normal densities.
sd = autoselectmixsd(z,nv = 100);
L = normlikmatrix(z,sd = sd);
size(L)
Now we fit the model using the three approaches.
@time xem, tem = mixEM(L,tol = 1e-4,maxiter = 1000);
@time xip, tip = REBayes(L);
@time outsqp = mixSQP(L,verbose = false);
In this example, the mix-SQP algorithm reaches a solution much faster than the both EM and IP approaches.
As before, the quality of the IP and SQP solutions is similar, whereas the EM solution is much worse.
fem = mixobjective(L,xem);
fip = mixobjective(L,xip);
fsqp = mixobjective(L,outsqp["x"]);
fbest = minimum([fem fip fsqp]);
@printf "Difference between EM and best solutions: %0.2e\n" fem - fbest
@printf "Difference between IP and best solutions: %0.2e\n" fip - fbest
@printf "Difference between SQP and best solutions: %0.2e\n" fsqp - fbest
The section gives information about the computing environment used to generate the results contained in this notebook, including the version of Julia, R and the packages used.
Pkg.status("Distributions")
Pkg.status("LowRankApprox")
Pkg.status("RCall")
versioninfo()
Since we called the KWDual function in R, it is also useful to record information about R, and the R packages used.
R"sessionInfo()"