glmfsa and glmnet

Selection accuracy at a matched model size

The comparison below is deliberately generous to glmnet: it fits its whole lambda path and we then pick the model whose size is closest to the true number of relevant variables – an oracle choice no real user could make. fsa is simply asked for k = 10.

compare <- function(d, k = 10) {
  f <- fsa(d$x, d$y, k = k)
  g <- glmnet(d$x, d$y)
  i <- which.min(abs(g$df - k))
  gs <- which(as.matrix(coef(g, s = g$lambda[i]))[-1, 1] != 0)
  c(fsa = length(intersect(f$selected[[1]], d$true)) / k,
    glmnet = length(intersect(gs, d$true)) / k)
}

reps <- t(replicate(10, compare(make_data(300, 2000))))
round(colMeans(reps), 3)
#>    fsa glmnet 
#>   0.99   0.51

On a strongly correlated AR(1) design with n = 300 and p = 2000, FSA recovers a larger fraction of the true variables. The gap widens as the design becomes harder; it closes on easy, well-separated problems where both methods simply get everything right.

Coefficient bias

d <- make_data(400, 1000)
f <- fsa(d$x, d$y, k = 10)
g <- glmnet(d$x, d$y)
i <- which.min(abs(g$df - 10))

data.frame(
  truth  = 1,
  fsa    = round(mean(coef(f)[d$true + 1, 1]), 3),
  glmnet = round(mean(as.matrix(coef(g, s = g$lambda[i]))[d$true + 1, 1]), 3))
#>   truth   fsa glmnet
#> 1     1 1.002  0.146

Every true coefficient equals 1. FSA’s estimates sit on it; the lasso’s are shrunk, by an amount that grows with lambda and therefore with how sparse you asked the model to be. Over a wider sweep (20 runs at n = 500, p = 2000) the averages are 1.00 for FSA against 0.24 for the lasso at the same model size.

Mislabeled data

This is the case with no glmnet equivalent at all. With 15% of the labels randomised, the Lorenz loss of the paper keeps finding the signal after the logistic loss has started chasing the noise.

noisy <- function(rep) {
  set.seed(rep)
  d <- make_data(400, 1000, task = "classification", noise = 0.15)
  g <- glmnet(d$x, d$y, family = "binomial")
  i <- which.min(abs(g$df - 10))
  gs <- which(as.matrix(coef(g, s = g$lambda[i]))[-1, 1] != 0)
  yb <- ifelse(d$y == 1, 1, -1)
  c(`fsa-logistic` = length(intersect(
      fsa(d$x, yb, family = "binomial", loss = "logistic",
          k = 10, s = 1e-4)$selected[[1]], d$true)),
    `fsa-lorenz` = length(intersect(
      fsa(d$x, yb, family = "binomial", loss = "lorenz",
          k = 10, s = 1e-4)$selected[[1]], d$true)),
    `glmnet-lasso` = length(intersect(gs, d$true))) / 10
}
round(colMeans(t(sapply(1:8, noisy))), 3)
#> fsa-logistic   fsa-lorenz glmnet-lasso 
#>        0.512        0.613        0.412

Speed

It is worth being blunt here, because the paper’s headline speed comparisons are against boosting and interior-point L1 solvers, not against glmnet’s coordinate descent, which is very fast indeed.

Matched on output – both asked for a path of 20 models:

tm <- function(p) {
  d <- make_data(1000, p)
  kg <- unique(round(exp(seq(log(100), log(2), length.out = 20))))
  c(p = p,
    fsa = system.time(fsa(d$x, d$y, k = kg))[["elapsed"]],
    glmnet = system.time(glmnet(d$x, d$y, nlambda = 20))[["elapsed"]])
}
round(as.data.frame(t(sapply(c(500, 2000, 8000), tm))), 3)
#>      p   fsa glmnet
#> 1  500 0.168  0.076
#> 2 2000 0.229  0.063
#> 3 8000 0.520  0.277

The two are of the same order; on sparse designs glmnet is several times faster. Speed is not the reason to prefer fsa. The reason is what you get for the time: at equal cost, a model that is both smaller and more accurate.

The end-to-end comparison a user actually pays for is a cross-validated fit. Across repeated runs at n = 2000 with ten true variables, cv.fsa and cv.glmnet both recovered all ten – but cv.glmnet needed roughly 55 to 62 variables to do it, against 11 or 12 for cv.fsa. That is the practical difference: an L1 path reaches full recovery only at a lambda loose enough to admit five times as many variables as it needs.

If you do want fsa to be fast, the annealing schedule is the dial. The cost is proportional to the area under the schedule curve, so a larger mu and a smaller niter cut it directly, and the paper’s Figure 6 shows that prediction performance is insensitive to both over a wide range:

d <- make_data(1000, 8000)
rbind(
  default = system.time(f1 <- fsa(d$x, d$y, k = 10))[1:3],
  fast    = system.time(f2 <- fsa(d$x, d$y, k = 10, mu = 1000,
                                  niter = 150))[1:3])
#>         user.self sys.self elapsed
#> default     1.661    0.095   0.343
#> fast        1.071    0.066   0.268
c(default = length(intersect(f1$selected[[1]], d$true)),
  fast    = length(intersect(f2$selected[[1]], d$true)))
#> default    fast 
#>      10      10

Real data

The package ships dexter and can fetch its larger companion gisette – both from the NIPS 2003 Feature Selection Challenge, and both deliberately padded with probes, random features inserted by the organizers so that false discoveries are measurable.

data(dexter)
dim(dexter$x)
#> [1]   300 20000
sum(Matrix::colSums(dexter$x != 0) == 0)   # columns with no support at all
#> [1] 12249

On a 300 x 20,000 design the ridge is not optional. Cross-validate it jointly with k:

set.seed(11)
fid <- sample(rep(1:5, length.out = nrow(dexter$x)))
cv <- cv.fsa(dexter$x, dexter$y, family = "binomial", loss = "lorenz",
             k = c(50, 20, 10, 5), s = c(0, 1e-3, 1e-2, 1e-1, 1),
             foldid = fid, type.measure = "class")
c(k = cv$k.min, s = cv$s.min)
#>  k  s 
#> 50  1

err <- function(p) mean((p > 0) != (as.character(dexter$yvalid) == "pos"))
g <- glmnet(dexter$x, dexter$y, family = "binomial")
pg <- predict(g, dexter$xvalid, type = "link")
sapply(c(5, 10, 20, 50), function(k) {
  f <- fsa(dexter$x, dexter$y, family = "binomial", loss = "lorenz",
           k = k, s = cv$s.min)
  i <- which.min(abs(g$df - k))
  c(glmfsa = err(drop(predict(f, dexter$xvalid, type = "link"))),
    lasso  = err(pg[, i]))
})
#>        [,1]      [,2]      [,3]       [,4]
#> glmfsa 0.28 0.1466667 0.1466667 0.13333333
#> lasso  0.23 0.1633333 0.1333333 0.08666667

The full picture across both datasets, from the benchmark in the package repository:

Validation error, percent. Gisette goes to glmfsa at every model size – a 10-variable glmfsa model beats a 20-variable lasso. Dexter changes hands with k: the Lorenz loss wins for very small models, the lasso for larger ones. With 300 observations and 7,751 usable columns, Dexter is exactly the regime where shrinkage pays, which is the next section’s point.

Where glmnet is the better tool

It would be a poor recommendation that did not say this plainly.

• The lasso is convex. Its solution path is unique and does not depend on initialization or step size. FSA’s constraint set is not convex; the theory guarantees convergence to a local optimum and selection consistency under a slow enough schedule, which is weaker.
• Small, well-conditioned problems. With p in the dozens and a design that is not badly correlated, coordinate descent is extremely hard to beat and the elastic net’s grouping behavior on correlated predictors is often exactly what you want.
• Maturity. glmnet has been in production use for over fifteen years and its edge cases are very well explored.
• You genuinely want shrinkage. If the goal is prediction under heavy noise rather than identifying a support, a shrunken estimator is frequently the better one. The Dexter result above is this case: n = 300 against 7,751 usable columns, where the lasso’s bias buys more than it costs once the model is allowed more than about 20 variables.

A reasonable division: use fsa when the model size is itself the thing you need to control, when p is large, when the labels may be wrong, or when you want unbiased coefficients without a refitting step. Use glmnet otherwise.

Reference

Barbu, A., She, Y., Ding, L. and Gramajo, G. (2017). Feature Selection with Annealing for Computer Vision and Big Data Learning. IEEE Transactions on Pattern Analysis and Machine Intelligence 39(2), 272-286.