Algorithm terminates early when x0 has some coordinates that are zero

[jwgarber]

I have a biological dataset where the algorithm incorrectly declares convergence after a single iteration when x0 = c(0, 0, 0, 0, 0, 0, 0, 0, 1).

Running mix-SQP algorithm 0.3-54 on 1000000 x 9 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):     10
number of EM iterations:    0
Computing SVD of 1000000 x 9 matrix.
Matrix is not low-rank; falling back to full matrix.
iter        objective max(rdual) nnz stepsize max.diff nqp nls
   1 +3.721806325e+00 -5.989e-07   1  ------   ------   --  --
Optimization took 0.07 seconds.
Convergence criteria met---optimal solution found.
V0 V1 V2 V3 V4 V5 V6 V7 V8 
 0  0  0  0  0  0  0  0  1 

Using x0 = rep(1/9, 9) converges to the correct solution.

Running mix-SQP algorithm 0.3-54 on 1000000 x 9 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):     10
number of EM iterations:    0
Computing SVD of 1000000 x 9 matrix.
Matrix is not low-rank; falling back to full matrix.
iter        objective max(rdual) nnz stepsize max.diff nqp nls
   1 +3.174361296e+00 +1.803e-01   9  ------   ------   --  --
   2 +3.138708847e+00 +1.778e-01   9 7.50e-01 1.59e-01   5   2
   3 +3.136557982e+00 +2.068e-02   9 1.00e+00 4.18e-02   2   1
   4 +3.136535232e+00 +3.823e-04   9 1.00e+00 3.29e-03   2   1
   5 +3.136535224e+00 +7.005e-08   9 1.00e+00 6.60e-05   2   1
   6 +3.136535224e+00 -3.489e-08   9 7.50e-01 2.26e-08   2   2
Optimization took 0.43 seconds.
Convergence criteria met---optimal solution found.
        V0         V1         V2         V3         V4         V5         V6 
0.06678421 0.03496157 0.12561704 0.03255271 0.12555548 0.02334131 0.20496663 
        V7         V8 
0.31307638 0.07314466

I believe this is because the termination condition in the main loop only checks the gradient at non-zero x_i, which on the first iteration is incorrect for x0 when some values are already zero.

    // Report on the algorithm's progress. Here we compute: the
    // smallest gradient value, or, equivalently, dual residual,
    // corresponding to the nonzero co-ordinates (gmin), which is used
    // as a convergence criterion; and the number of nonzeros in the
    // solution (nnz). Note that only the dual residuals (gmin's)
    // corresponding to the nonzero co-ordinates are relevant.
    gmin(i) = 1 + g(j1).min();

    // Check convergence. Convergence is reached with the maximum dual
    // residual is small. The negative of "gmin" is also the maximum
    // dual residual (denoted as "rdual" on p. 609 of Boyd &
    // Vandenberghe, "Convex Optimization", 2009).
    if (gmin(i) >= -convtolsqp) {
      status = 0;
      i++;
      break;
    }

[pcarbo]

@jwgarber By default there should be an initial "EM phase", and it appears that your mixsqp call has skipped that? This might at least provide a workaround for the issue you are experiencing.

If you can share your data and code, I can try to reproduce the error you are getting. It is quite possible this is a bug. Thanks.

[jwgarber]

Thanks, yes I've disabled the initial EM runs since they don't affect the number of SQP iterations for my data (though enabling it does fix the problem). It seems this can be reproduced with the default simulated data.

library(mixsqp)

set.seed(1)

L <- simulatemixdata(1000,20)$L
fit.sqp <- mixsqp(L, x0=c(rep(0.0, 19), 1.0), control=list(numiter.em=0))

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:    0
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 +2.664071413e+00 -1.934e-07   1  ------   ------   --  --
Optimization took 0.00 seconds.
Convergence criteria met---optimal solution found.

[pcarbo]

@jwgarber Yes you are right. I think there is a particular difficulty starting with solutions that are on the boundary of the constraint set. Even nudging the solution slightly away from the boundary can help with this problem, e.g.,

set.seed(1)
e <- 1e-5
L <- simulatemixdata(1000,20)$L
fit <- mixsqp(L,x0 = c(rep(e,19),1),control = list(numiter.em = 0))

I'm not sure this is strictly speaking a "bug", but it is a limitation of the method. If you leave this Issue open, it will act as a reminder to look into it and see if it can be addressed.