stephenslab / mixsqp

R package for fast maximum-likelihood estimation of mixture proportions using SQP.
https://bit.ly/2NtYHWT
Other
11 stars 7 forks source link

NPMLE example where mixsqp fails to converge #42

Closed pcarbo closed 4 years ago

pcarbo commented 4 years ago

Thanks to @willwerscheid for this example. The mixsqp solution is only very slightly worse than the KWDual solution, but fails to satisfy the convergence criterion.

library(mixsqp)
load("npmle_mixsqp.RData")
out1 <- mixkwdual(L)
out2 <- mixsqp(L,control = list(maxiter.sqp = 100))
cat(sprintf("objective value at KWDual solution: %0.12f\n",out1$value))
cat(sprintf("objective value at mixsqp solution: %0.12f\n",out2$value))
i <- which(out2$x > 0)
y <- eigen(out2$hessian[i,i])$values
plot(seq(1,length(i)),y,pch = 20,col = "dodgerblue",log = "y",
     xlab = "nonzero co-ordinate",ylab = "eigenvalue")

npmle_mixsqp.RData.gz

pcarbo commented 4 years ago

With commit 98f1fa0 (version 0.3-37), the output should look something like this:

library(mixsqp)
load("npmle_mixsqp.RData")
out1 <- mixkwdual(L)
out2 <- mixsqp(L)
cat(sprintf("KWDual solution: %0.12f\n",mixobjective(L,out1$x)))
cat(sprintf("mixsqp solution: %0.12f\n",mixobjective(L,out2$x)))
# Running mix-SQP algorithm 0.3-37 on 10000 x 100 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 10000 x 100 matrix.
# Matrix is not low-rank; falling back to full matrix.
# iter        objective max(rdual) nnz stepsize max.diff nqp nls
#    1 +2.022268957e+00  -- EM --  100 1.00e+00 2.03e-02  --  --
#    2 +2.003455962e+00  -- EM --  100 1.00e+00 6.47e-03  --  --
#    3 +1.999928739e+00  -- EM --  100 1.00e+00 3.27e-03  --  --
#    4 +1.998632776e+00  -- EM --  100 1.00e+00 2.04e-03  --  --
#    5 +1.997957461e+00  -- EM --  100 1.00e+00 1.46e-03  --  --
#    6 +1.997546459e+00  -- EM --  100 1.00e+00 1.13e-03  --  --
#    7 +1.997277008e+00  -- EM --  100 1.00e+00 9.24e-04  --  --
#    8 +1.997091797e+00  -- EM --  100 1.00e+00 7.78e-04  --  --
#    9 +1.996959463e+00  -- EM --  100 1.00e+00 6.68e-04  --  --
#   10 +1.996861483e+00  -- EM --  100 1.00e+00 5.81e-04  --  --
#    1 +1.996786460e+00 +2.844e-02 100  ------   ------   --  --
#    2 +1.996726676e+00 +2.290e-02  80 1.00e+00 2.72e-03  20   1
#    3 +1.996678436e+00 +2.105e-02  60 1.00e+00 2.14e-02  20   1
#    4 +1.996604077e+00 +1.313e-02  40 1.00e+00 9.67e-02  20   1
#    5 +1.995947165e+00 +4.678e-03  26 1.00e+00 4.41e-02  20   1
#    6 +1.995896017e+00 +1.397e-03  24 1.00e+00 7.86e-03  20   1
#    7 +1.995880170e+00 +1.145e-05  25 1.00e+00 2.07e-02  20   1
#    8 +1.995875775e+00 +2.871e-06  27 1.00e+00 2.58e-02  20   1
#    9 +1.995873237e+00 +4.126e-06  27 1.00e+00 4.86e-02  20   1
#   10 +1.995871719e+00 +1.707e-06  28 1.00e+00 6.42e-02  20   1
#   11 +1.995867351e+00 +4.008e-06  28 1.00e+00 1.33e-01  20   1
#   12 +1.995867351e+00 -5.928e-08  28 1.00e+00 2.57e-05  20   1
# Optimization took 1.10 seconds.
# Convergence criteria met---optimal solution found.
# KWDual solution: 1.995865487812
# mixsqp solution: 1.995867491013