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

Slow performance using 0.1-119 #31

Open jean997 opened 5 years ago

jean997 commented 5 years ago

I am having an issue with the new version of mixsqp. I have found it to be mostly faster however I seem to have certain instances where it takes a very very long time. I can file a github issue but I thought I would post here as well. I will attach my test data that replicates this problem. To test it I installed two different versions of mixsqp in different directories. Here is the code I used with output:

library(mixsqp, lib.loc="~/mixsqp/0.1-97")
#library(mixsqp, lib.loc="~/mixsqp/master")
library(ashr)

dat <- readRDS("mixsqp_debug.RDS")

system.time(
 w_res <- with(dat, ashr::mixSQP(matrix_lik =matrix_lik,
                                            prior=c(null_wt, rep(1, K-1)),
                                            weights=rep(1, nrow(matrix_lik))))
)

The top two lines are for switching versions. Using 0.1-97 it completes in 11-15 seconds. Using version 0.1-119 I have not let it run long enough to finish but it is more than 10 minutes.

The data file seems to be too big to submit via github.

pcarbo commented 5 years ago

@jean997 I can confirm that mix-SQP is getting stuck:

library(mixsqp)
load("mixsqp_jean.RData")
out1 <- mixsqp(L,w,x0,control = list(maxiter.sqp = 24))
out2 <- mixkwdual(L,w)

Here's the console output:

Running mix-SQP algorithm 0.1-119 on 842860 x 42 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-15
l.s. sufficient decrease:   1.0e-02
step size reduction factor: 7.5e-01
minimum step size:          1.0e-08
max. iter (SQP):            24
max. iter (active-set):     43
number of EM iterations:    4
iter        objective max(rdual) nnz stepsize max.diff nqp nls
   1 +5.628443405e-01  -- EM --   42 1.00e+00 2.33e-02  --  --
   2 +4.460648143e-01  -- EM --   42 1.00e+00 1.94e-02  --  --
   3 +3.955883104e-01  -- EM --   42 1.00e+00 1.53e-02  --  --
   4 +3.702957467e-01  -- EM --   42 1.00e+00 1.26e-02  --  --
   1 +3.702957467e-01 +1.162e-01  42  ------   ------   --  --
   2 +3.972471946e-01 +2.754e+00  42 9.90e-01 6.15e-01  43   1
   3 +3.310772356e-01 +1.266e+00  35 9.90e-01 1.47e-01  43   1
   4 +3.183761190e-01 +5.023e-01  23 9.90e-01 1.44e-01  43   1
   5 +3.142609116e-01 +1.141e-01  13 9.90e-01 2.57e-01  43   1
   6 +3.134638969e-01 +1.667e-02  11 9.90e-01 1.17e-01  43   1
   7 +3.133306492e-01 +2.563e-03   9 9.90e-01 6.78e-02  43   1
   8 +3.133266812e-01 +7.871e-04   9 9.90e-01 2.98e-02  43   1
   9 +3.133183983e-01 +2.767e-04   9 9.90e-01 4.68e-03  43   1
  10 +3.133180167e-01 +2.715e-04   9 9.90e-01 1.87e-04  43   1
  11 +3.133180124e-01 +2.714e-04   9 9.90e-01 2.10e-06  43   1
  12 +3.133180124e-01 +2.714e-04   9 9.90e-01 2.11e-08  43   1
  13 +3.133180124e-01 +2.714e-04   9 9.90e-01 5.42e-10  43   1
  14 +3.133180124e-01 +2.714e-04   9 2.35e-01 1.26e-10  43   6
  15 +3.133180124e-01 +2.714e-04   9 9.91e-02 5.31e-13  43   9
  16 +3.133180124e-01 +2.714e-04   9 3.13e-01 5.66e-12  43   5
  17 +3.133180124e-01 +2.714e-04   9 9.91e-02 1.42e-12  43   9
  18 +3.133180124e-01 +2.714e-04   9 3.14e-02 8.28e-13  43  13
  19 +3.133180124e-01 +2.714e-04   9 7.45e-04 2.86e-14  43  26
  20 +3.133180124e-01 +2.714e-04   9 1.33e-04 5.86e-15  43  32
  21 +3.133180124e-01 +2.714e-04   9 7.44e-03 2.33e-13  43  18
  22 +3.133180124e-01 +2.714e-04   9 3.15e-07 6.94e-18  43  53
  23 +3.133180124e-01 +2.714e-04   9 1.77e-04 5.41e-15  43  31
  24 +3.133180124e-01 +2.714e-04   9 5.61e-08 5.96e-27  43  59
Failed to converge within iterations limit.

See here for the data used to run this example.

pcarbo commented 5 years ago

@jean997 However, I should note the two solutions are very close to each other in terms of objective value:

> max(abs(out1$x - out2$x))
[1] 0.01287
> print(out1$value,digits = 12)
[1] 0.313318012633
> print(out2$value,digits = 12)
[1] 0.313315465836

I don't know exactly why mix-SQP isn't converging successfully, and it may take some time to pinpoint the exact issue---it does seem like a poorly behaved objective surface. As a workaround, I suggest providing a safer upper limit on the number of SQP iterations.