Here we explore the use of active-set and interior-point methods (the latter implemented by the commercial software MOSEK) for solving the quadratic subproblem inside SQP.
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, Mosek and JuMP Julia packages, as well as some function definitions used in the code chunks below.
using Pkg
using Random
using Printf
using Distributions
using SparseArrays
using LinearAlgebra
using Mosek
using JuMP
include("../code/datasim.jl");
include("../code/likelihood.jl");
include("../code/mixSQP.jl");
Next, initialize the sequence of pseudorandom numbers.
Random.seed!(1);
Let's begin with a smaller example with 50,000 samples.
n = round(Int,5e4);
z = normtmixdatasim(n);
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 the mix-SQP algorithm a couple of times to precompile the relevant functions.
out = mixSQP(L,qpsubprob = "activeset",lowrank = "none",maxiter = 10,verbose = false);
out = mixSQP(L,qpsubprob = "mosek",lowrank = "none",maxiter = 10,verbose = false);
Fit the model using the SQP algorithm, with an active-set method to find the solution to the quadratic program at each SQP iteration.
@time out1 = mixSQP(L,qpsubprob = "activeset",lowrank = "none");
Next fit the model again using the same SQP algorithm, with the active-set method replaced by MOSEK.
@time out2 = mixSQP(L,qpsubprob = "mosek",lowrank = "none");
Both runs converged to a solution in a small number of iterations. The solutions are very similar:
maximum(abs.(out1["x"] - out2["x"]))
Solving the quadratic programs is only a small fraction of the total effort. Nonetheless, the effort of the active-set computations is much less than solving the quadratic subproblems with MOSEK.
@printf "Total runtime of active set method: %0.3f s.\n" sum(out1["qptiming"])
@printf "Total runtime of interior point method: %0.3f s.\n" sum(out2["qptiming"])
Let's now explore the accuracy and runtime of the active-set and MOSEK solvers in a larger data set.
z = normtmixdatasim(round(Int,1e5));
As before, we compute the $n \times k$ conditional likelihood matrix for a mixture of zero-centered normals. This time, we use a finer grid of $k = 40$ normal densities to compute this matrix.
k = 40;
sd = autoselectmixsd(z,nv = k);
L = normlikmatrix(z,sd = sd);
size(L)
Now we fit the model using the two variants of the SQP algorithm.
@time out1 = mixSQP(L,qpsubprob = "activeset",lowrank = "none",verbose = false);
@time out2 = mixSQP(L,qpsubprob = "mosek",lowrank = "none",verbose = false);
The first SQP run with the active-set method is only slightly faster. And, as before, the solutions are very similar:
maximum(abs.(out1["x"] - out2["x"]))
The amount of time spent solving the quadratic programs is again only a small proportion of the total:
@printf "Total runtime of active set method: %0.3f s.\n" sum(out1["qptiming"])
@printf "Total runtime of interior point method: %0.3f s.\n" sum(out2["qptiming"])
Therefore, although the active-set method is much faster than MOSEK, the overall impact on performance is relatively small.
The section gives information about the computing environment used to generate the results contained in this notebook, including the version of Julia, and the versions of the Julia packages used here.
Pkg.status("Distributions");
Pkg.status("Mosek");
Pkg.status("JuMP");
versioninfo()