Maximum-likelihood estimation of mixture proportions using SQP

The mixsqp function uses a Sequential Quadratic Programming (SQP) algorithm to find the maximum likelihood estimates of mixture proportions in a (finite) mixture model. More generally, mixsqp solves the corresponding constrained, convex optimization problem, which is given below (see ‘Details’). See ‘References’ for more details about the SQP algorithm.

mixsqp(
  L,
  w = rep(1, nrow(L)),
  x0 = rep(1, ncol(L)),
  log = FALSE,
  control = list()
)

mixsqp_control_default()

Arguments

L	Matrix specifying the optimization problem to be solved. In the context of mixture-model fitting, `L[j,k]` should be the value of the kth mixture component density at the jth data point. `L` should be a numeric matrix with at least two columns, with all entries being non-negative and finite (and not missing). In some cases, it is easier or more natural to compute `log(L)`; for example, it is often easier to compute the log-likelihood rather than the likelihood. Setting `log = TRUE` will tell `mixsqp` to interpret this input as the logarithm of the data matrix. Note that, for large matrices, it is preferrable that the matrix is stored in double-precision; see `storage.mode`.
w	An optional numeric vector, with one entry for each row of `L`, specifying the "weights" associated with the rows of `L`. All weights must be finite, non-negative and not missing. Internally, the weights are normalized to sum to 1, which does not change the problem, but does change the value of the objective function reported. By default, all weights are equal.
x0	An optional numeric vector providing an initial estimate of the solution to the optimization problem. It should contain only finite, non-missing, non-negative values, and all entries of `L %*% x0` must be greater than zero (to ensure that the objective evaluates to a finite value at `x0`). The vector will be normalized to sum to 1. By default, `x0` is the vector with all equal values.
log	When `log = TRUE`, the input matrix `L` is interpreted as containing the logarithm of the data matrix.
control	A list of parameters controlling the behaviour of the optimization algorithm. See ‘Details’.

Value

A list object with the following elements:

If the SQP algorithm converges, this is the solution to the convex optimization problem. If the algorithm fails to converge, it is the best estimate of the solution achieved by the algorithm. Note that if the SQP algorithm terminates before reaching the solution, x may not satisfy the equality constraint; that is, the entries of x may not sum to 1.

value

The value of the objective function, $f(x)$, at x.

status

A character string describing the status of the algorithm upon termination.

progress

A data frame containing more detailed information about the algorithm's progress. The data frame has one row per SQP iteration. For an explanation of the columns, see the description of the verbose control parameter in ‘Details’. Missing values (NA's) in the last row indicate that these quantities were not computed because convergence was reached before computing them. Also note that the storage of these quantities in the progress data frame is slightly different than in the console output (when verbose = TRUE) as the console output shows some quantities that were computed after the convergence check in the previous iteration.

Details

mixsqp solves the following optimization problem. Let $L$ be a matrix with $n$ rows and $m$ columns containing only non-negative entries, and let $w = (w_1, \ldots, w_n)$ be a vector of non-negative "weights". mixsqp computes the value of vector $x = (x_1, \ldots, x_m)$ minimizing the following objective function, $$f(x) = -\sum_{j=1}^n w_j \log (\sum_{k=1}^m L_{jk} x_k),$$ subject to the constraint that $x$ lie within the simplex; that is, the entries of $x$ are non-negative and sum to 1. Implicitly, there is an additional constraint $L*x > 0$ in order to ensure that the objective has a finite value. In practice, this constraint only needs to be checked for the initial estimate to ensure that it holds for all subsequent iterates.

If all weights are equal, solving this optimization problem corresponds to finding the maximum-likelihood estimate of the mixture proportions $x$ given $n$ independent data points drawn from a mixture model with $m$ components. In this case, $L_{jk}$ is the likelihood for mixture component $k$ and data point $j$.

The Expectation Maximization (EM) algorithm can be used to solve this optimization problem, but it is intolerably slow in many interesting cases, and mixsqp is much faster.

A special feature of this optimization problem is that the gradient of the objective does not change with re-scaling; for example, if all the entries of matrix L are multiplied by 100, the gradient does not change. A practical benefit of this property is that the optimization algorithm will behave similarly irrespective of the scale of L; for example, the same value for the convergence tolerance convtol.sqp will have the same effect at different scales.

A related feature is that the solution to the optimization problem is invariant to rescaling the rows of L; for example, the solution will remain the same after all the entries in a row of L are multiplied by 10. A simple normalization scheme divides each row by the largest entry in the row so that all entries of L are at most 1: L <- L / apply(L,1,max) Occasionally, it can be helpful to normalize the rows when some of the entries are unusually large or unusually small. This can help to avoid numerical overflow or underflow errors.

The SQP algorithm is implemented using the Armadillo C++ linear algebra library, which can automatically take advantage of multithreaded matrix computations to speed up mixsqp for large L matrices, but only when R has been configured with a multithreaded BLAS/LAPACK library (e.g., OpenBLAS).

The control argument is a list in which any of the following named components will override the default optimization algorithm settings (as they are defined by mixsqp_control_default):

normalize.rows: When normalize.rows = TRUE, the rows of the data matrix L are automatically scaled so that the largest entry in each row is 1. This is the recommended setting for better stability of the optimization. When log = TRUE, this setting is ignored becase the rows are already normalized. Note that the objective is computed on the original (unnormalized) matrix to make the results easier to interpret.
force.sparse.init: The computational complexity of the SQP updates grows rapidly by the number of nonzeros in the current solution estimate. When force.sparse.init = TRUE, and there are more than 20 nonzeros in the solution before beginning the SQP updates, only the 20 largest entries are kept, the rest are forced to zero.
tol.svd: Setting used to determine rank of truncated SVD approximation for L. The rank of the truncated singular value decomposition is determined by the number of singular values surpassing tol.svd. When tol.svd = 0 or when L has 4 or fewer columns, all computations are performed using full L matrix.
convtol.sqp: A small, non-negative number specifying the convergence tolerance for SQP algorithm; convergence is reached when the maximum dual residual in the Karush-Kuhn-Tucker (KKT) optimality conditions is less than or equal to convtol.sqp. Smaller values will result in more stringent convergence criteria and more accurate solutions, at the expense of greater computation time. Note that changes to this tolerance parameter may require respective changes to convtol.activeset and/or zero.threshold.searchdir to obtain reliable convergence.
convtol.activeset: A small, non-negative number specifying the convergence tolerance for the active-set step. Smaller values will result in higher quality search directions for the SQP algorithm but possibly a greater per-iteration computational cost. Note that changes to this tolerance parameter can affect how reliably the SQP convergence criterion is satisfied, as determined by convtol.sqp.
zero.threshold.solution: A small, non-negative number used to determine the "active set"; that is, it determines which entries of the solution are exactly zero. Any entries that are less than or equal to zero.threshold.solution are considered to be exactly zero. Larger values of zero.threshold.solution may lead to speedups for matrices with many columns, at the (slight) risk of prematurely zeroing some co-ordinates.
zero.threshold.searchdir: A small, non-negative number used to determine when the search direction in the active-set step is considered "small enough". Note that changes to this tolerance parameter can affect how reliably the SQP convergence criterion is satisfied, as determined by convtol.sqp, so choose this parameter carefully.
suffdecr.linesearch: This parameter determines how stringent the "sufficient decrease" condition is for accepting a step size in the backtracking line search. Larger values will make the condition more stringent. This should be a positive number less than 1.
stepsizereduce: The multiplicative factor for decreasing the step size in the backtracking line search. Smaller values will yield a faster backtracking line search at the expense of a less fine-grained search. Should be a positive number less than 1.
minstepsize: The smallest step size accepted by the line search step. Should be a number greater than 0 and at most 1.
identity.contrib.increase: When the Hessian is not positive definite, a multiple of the identity is added to the Hessian to ensure a unique search direction. The factor for increasing the identity contribution in this modified Hessian is determined by this control parameter.
eps: A small, non-negative number that is added to the terms inside the logarithms to sidestep computing logarithms of zero. This prevents numerical problems at the cost of introducing a small inaccuracy in the solution. Increasing this number may lead to faster convergence but possibly a less accurate solution.
maxiter.sqp: Maximum number of SQP iterations to run before reporting a convergence failure; that is, the maximum number of quadratic subproblems that will be solved by the active-set method.
maxiter.activeset: Maximum number of active-set iterations taken to solve each of the quadratic subproblems. If NULL, the maximum number of active-set iterations is set to min(20,1 + ncol(L)).
numiter.em: Number of expectation maximization (EM) updates to perform prior to running mix-SQP. This can help ensure convergence of mix-SQP when the initial solution is very poor.
verbose: If verbose = TRUE, the algorithm's progress and a summary of the optimization settings are printed to the console. The algorithm's progress is displayed in a table with one row per SQP (outer loop) iteration, and with the following columns: "iter", the (outer loop) SQP iteration; "objective", the value of the objective function (see $f(x)$) at the current estimate of the solution, $x$; "max(rdual)", the maximum "dual residual" in the Karush-Kuhn-Tucker (KKT) conditions, which is used to monitor convergence (see convtol.sqp); "nnz", the number of non-zero co-ordinates in the current estimate, as determined by zero.threshold.solution; "max.diff", the maximum difference in the estimates between two successive iterations; "nqp", the number of (inner loop) active-set iterations taken to solve the quadratic subproblem; "nls", the number of iterations in the backtracking line search.

References

Y. Kim, P. Carbonetto, M. Stephens and M. Anitescu (2018). A fast algorithm for maximum likelihood estimation of mixture proportions using sequential quadratic programming. arXiv:1806.01412 https://arxiv.org/abs/1806.01412.

Examples

set.seed(1)
n <- 1e5
m <- 10
w <- rep(1,n)/n
L <- simulatemixdata(n,m)$L
out.mixsqp <- mixsqp(L,w)
#> Running mix-SQP algorithm 0.3-17 on 100000 x 10 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-10
#> 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):     11
#> number of EM iterations:    10
#> Computing SVD of 100000 x 10 matrix.
#> Matrix is not low-rank; falling back to full matrix.
#> iter        objective max(rdual) nnz stepsize max.diff nqp nls
#>    1 +6.649406794e-01  -- EM --   10 1.00e+00 7.79e-02  --  --
#>    2 +6.336404241e-01  -- EM --   10 1.00e+00 2.61e-02  --  --
#>    3 +6.201431305e-01  -- EM --   10 1.00e+00 2.41e-02  --  --
#>    4 +6.131770842e-01  -- EM --   10 1.00e+00 2.06e-02  --  --
#>    5 +6.091788165e-01  -- EM --   10 1.00e+00 1.71e-02  --  --
#>    6 +6.066948308e-01  -- EM --   10 1.00e+00 1.42e-02  --  --
#>    7 +6.050509636e-01  -- EM --   10 1.00e+00 1.18e-02  --  --
#>    8 +6.039062025e-01  -- EM --   10 1.00e+00 9.95e-03  --  --
#>    9 +6.030755992e-01  -- EM --   10 1.00e+00 8.45e-03  --  --
#>   10 +6.024525768e-01  -- EM --   10 1.00e+00 7.25e-03  --  --
#>    1 +6.024525768e-01 +2.365e-02  10  ------   ------   --  --
#>    2 +6.058450302e-01 +4.073e+00   4 1.00e+00 3.91e-01  11   1
#>    3 +5.999690521e-01 +1.894e+00   5 1.00e+00 3.21e-02  10   1
#>    4 +5.997818561e-01 +7.929e-01   5 1.00e+00 2.95e-04   2   1
#>    5 +5.997500584e-01 +2.749e-01   5 1.00e+00 1.98e-03   2   1
#>    6 +5.996811464e-01 +1.482e-01   5 1.00e+00 5.64e-03   2   1
#>    7 +5.994718234e-01 +6.121e-02   5 1.00e+00 9.54e-03   2   1
#>    8 +5.991459366e-01 +1.505e-02   5 1.00e+00 8.23e-03   2   1
#>    9 +5.989768704e-01 +1.193e-03   5 1.00e+00 2.88e-03   2   1
#>   10 +5.989591648e-01 +8.215e-06   5 1.00e+00 2.42e-04   2   1
#>   11 +5.989590364e-01 +3.904e-10   5 1.00e+00 1.63e-06   2   1
#> Optimization took 0.54 seconds.
#> Convergence criteria met---optimal solution found.
f <- mixobjective(L,out.mixsqp$x,w)
print(f,digits = 16)
#> [1] 0.5989590363728419

# We can also compare this result with solution found from an 
# interior-point approach called via the "KWDual" function from the 
# REBayes package. (This requires installation of the MOSEK 
# optimization library as well as the REBayes package, so we have
# made this step optional.)
if (FALSE) {
out.kwdual <- mixkwdual(L,w)
print(mixobjective(L,out.kwdual$x,w),digits = 16)
}