Illustration of mixsqp applied to a small data set, and a large one

Youngseok Kim, Peter Carbonetto and Matthew Stephens

2023-12-20

In this vignette, we illustrate the use of the sequential quadratic programming (SQP) algorithm implemented in mixsqp.

Environment set-up

Load the mixsqp package.

library(mixsqp)

Next, initialize the sequence of pseudorandom numbers.

set.seed(1)

Generate a small data set

We begin with a small example to show how mixsqp works.

L <- simulatemixdata(1000,20)$L
dim(L)
# [1] 1000   20

This call to simulatemixdata created an \(n \times m\) conditional likelihood matrix for a mixture of zero-centered normals, with \(n = 1000\) and \(m = 20\). By default, simulatemixdata normalizes the rows of the likelihood matrix so that the maximum entry in each row is 1.

Fit mixture model

Now we fit the mixture model using the SQP algorithm:

fit.sqp <- mixsqp(L)
# Running mix-SQP algorithm 0.3-54 on 1000 x 20 matrix
# convergence tol. (SQP):     1.0e-08
# conv. tol. (active-set):    1.0e-10
# zero threshold (solution):  1.0e-08
# zero thresh. (search dir.): 1.0e-14
# l.s. sufficient decrease:   1.0e-02
# step size reduction factor: 7.5e-01
# minimum step size:          1.0e-08
# max. iter (SQP):            1000
# max. iter (active-set):     20
# number of EM iterations:    10
# Computing SVD of 1000 x 20 matrix.
# Matrix is not low-rank; falling back to full matrix.
# iter        objective max(rdual) nnz stepsize max.diff nqp nls
#    1 +6.825854400e-01  -- EM --   20 1.00e+00 3.43e-02  --  --
#    2 +6.608901094e-01  -- EM --   20 1.00e+00 1.12e-02  --  --
#    3 +6.501637569e-01  -- EM --   20 1.00e+00 8.83e-03  --  --
#    4 +6.441429345e-01  -- EM --   20 1.00e+00 7.64e-03  --  --
#    5 +6.405379612e-01  -- EM --   20 1.00e+00 6.44e-03  --  --
#    6 +6.382623445e-01  -- EM --   20 1.00e+00 5.36e-03  --  --
#    7 +6.367520429e-01  -- EM --   20 1.00e+00 4.46e-03  --  --
#    8 +6.357009493e-01  -- EM --   20 1.00e+00 3.75e-03  --  --
#    9 +6.349366492e-01  -- EM --   20 1.00e+00 3.18e-03  --  --
#   10 +6.343584376e-01  -- EM --   20 1.00e+00 2.73e-03  --  --
#    1 +6.339053898e-01 +1.856e-02  20  ------   ------   --  --
#    2 +6.281996199e-01 +1.384e-03   4 1.00e+00 4.36e-01  20   1
#    3 +6.281978243e-01 +8.849e-07   4 1.00e+00 3.56e-03   2   1
#    4 +6.281978243e-01 -1.816e-08   4 1.00e+00 6.59e-06   2   1
# Optimization took 0.00 seconds.
# Convergence criteria met---optimal solution found.

In this example, the SQP algorithm converged to a solution in a small number of iterations.

By default, mixsqp outputs information on its progress. It begins by summarizing the optimization problem and the algorithm settings used. (Since we did not change these settings in the mixsqp call, all the settings shown here are the default settings.)

After that, it outputs, at each iteration, information about the current solution, such as the value of the objective (“objective”) and the number of nonzeros (“nnz”).

The “max(rdual)” column shows the quantity used to assess convergence. It reports the maximum value of the “dual residual”; the SQP solver terminates when the maximum dual residual is less than conv.tol, which by default is \(10^{-8}\). In this example, we see that the dual residual shrinks rapidly toward zero.

Another useful indicator of convergence is the “max.diff” column—it reports the maximum difference between the solution estimates at two successive iterations. We normally expect these differences to shrink as we approach the solution, which is precisely what we see in this example.

This information is also provided in the return value, which we can use, for example, to create a plot of the objective value at each iteration of the SQP algorithm:

numiter <- nrow(fit.sqp$progress)
plot(1:numiter,fit.sqp$progress$objective,type = "b",
     pch = 20,lwd = 2,xlab = "SQP iteration",
     ylab = "objective",xaxp = c(1,numiter,numiter - 1))

Session information

This next code chunk gives information about the computing environment used to generate the results contained in this vignette, including the version of R and the packages used.

sessionInfo()
# R version 3.6.2 (2019-12-12)
# Platform: x86_64-apple-darwin15.6.0 (64-bit)
# Running under: macOS Catalina 10.15.7
# 
# Matrix products: default
# BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
# LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
# 
# locale:
# [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
# 
# attached base packages:
# [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
# [1] mixsqp_0.3-54
# 
# loaded via a namespace (and not attached):
#  [1] Rcpp_1.0.8      lattice_0.20-38 digest_0.6.23   grid_3.6.2     
#  [5] R6_2.4.1        jsonlite_1.7.2  magrittr_2.0.1  evaluate_0.14  
#  [9] highr_0.8       stringi_1.4.3   rlang_1.0.6     cli_3.5.0      
# [13] irlba_2.3.3     jquerylib_0.1.4 Matrix_1.3-4    bslib_0.3.1    
# [17] rmarkdown_2.21  tools_3.6.2     stringr_1.4.0   xfun_0.36      
# [21] yaml_2.2.0      fastmap_1.1.0   compiler_3.6.2  htmltools_0.5.4
# [25] knitr_1.37      sass_0.4.0