Getting started with glmfsa

The idea in one paragraph

Feature Selection with Annealing fits a linear model under an exact constraint on how many coefficients may be non-zero. It starts from beta = 0 and repeats two steps: take a gradient step on the loss, then throw away the coefficients with the smallest magnitude, according to a schedule that tightens from all M variables down to the k you asked for. Variables that are removed leave the working design matrix altogether, so each iteration is cheaper than the last.

That is the whole algorithm. It comes from Barbu, She, Ding and Gramajo (2017), along with a theorem guaranteeing that it converges and that the surviving support contains the true one with probability tending to one.

A first fit

n <- 300; p <- 1000
x <- matrix(rnorm(n * p), n, p)
beta <- numeric(p); beta[c(10, 50, 100, 200, 300)] <- 1.5
y <- as.numeric(x %*% beta + rnorm(n))

fit <- fsa(x, y, family = "gaussian", k = 5)
fit
#> 
#> Call: fsa(x = x, y = y, family = "gaussian", k = 5)
#> 
#> Feature Selection with Annealing -- family: gaussian
#> 300 observations, 1000 variables; fitted in 0.023 s
#> 
#>  k Df  %Dev   Loss
#>  5  5 91.84 0.4597
which(coef(fit)[-1] != 0)
#> [1]  10  50 100 200 300

Five variables requested, five returned, and they are the five that generated the data. Note what you did not have to do: there was no penalty parameter to choose, no path to search for the model size you wanted, and no refitting step to remove shrinkage bias from the surviving coefficients.

round(coef(fit)[c(1, 1 + c(10, 50, 100, 200, 300)), ], 3)
#> (Intercept)         V10         V50        V100        V200        V300 
#>       0.113       1.548       1.503       1.370       1.427       1.486

The estimates sit on the true value of 1.5. A lasso fit at the same model size would return coefficients pulled systematically toward zero, which is why L1 users routinely refit an unpenalized model on the selected variables.

Watching the algorithm run

keep_trace = TRUE records every coefficient at every iteration, which plot_trace() then draws. This reproduces Figure 2 of the paper and is the quickest diagnostic there is: curves should have flattened well before the last iteration.

tr <- fsa(x, y, k = 5, niter = 200, keep_trace = TRUE)
plot_trace(tr, highlight = c(10, 50, 100, 200, 300))

The gray curves are nuisance variables being eliminated; the dashed red line (right axis) is the number of variables still in the model. By iteration 100 the problem has shrunk from 1000 variables to 5, and the remaining coefficients simply converge.

A whole path of model sizes

Pass a decreasing vector to k and you get one model per entry from a single annealed pass. The supports are nested, and the path costs barely more than the single fit at the smallest k.

path <- fsa(x, y, k = c(50, 30, 20, 12, 8, 5, 3, 2, 1))
path
#> 
#> Call: fsa(x = x, y = y, k = c(50, 30, 20, 12, 8, 5, 3, 2, 1))
#> 
#> Feature Selection with Annealing -- family: gaussian
#> 300 observations, 1000 variables; fitted in 0.028 s
#> 
#>   k Df  %Dev   Loss
#>  50 50 95.48 0.2546
#>  30 30 94.36 0.3178
#>  20 20 93.68 0.3560
#>  12 12 92.90 0.3999
#>   8  8 92.40 0.4284
#>   5  5 91.84 0.4597
#>   3  3 60.86 2.2060
#>   2  2 47.07 2.9830
#>   1  1 27.08 4.1090
plot(path, label = TRUE)

Choosing k

cv.fsa() works like cv.glmnet(), with k in place of lambda.

cvfit <- cv.fsa(x, y, k = c(30, 20, 12, 8, 5, 3, 2, 1), nfolds = 5)
cvfit
#> 
#> Call: cv.fsa(x = x, y = y, k = c(30, 20, 12, 8, 5, 3, 2, 1), nfolds = 5)
#> 
#> Measure: Deviance  
#> 
#>     k Measure      SE Nonzero
#> min 5  0.9703 0.05045       5
#> 1se 5  0.9703 0.05045       5
plot(cvfit)

k.min minimizes the cross-validated deviance; k.1se is the smallest model within one standard error of it – the same conservative choice lambda.1se makes, oriented the way model size runs.

head(predict(cvfit, newx = x, k = "k.1se"))
#>              k5
#> [1,]  2.1478843
#> [2,]  6.9527092
#> [3,]  1.8817696
#> [4,] -1.3607630
#> [5,]  0.6662497
#> [6,] -4.1652220

The tuning parameters

Three parameters control the anneal, and the paper’s Figure 6 shows performance is flat across roughly an order of magnitude in each. They are worth understanding but rarely worth tuning.

• niter – how long the anneal runs. The schedule reaches k at niter/2, so the second half is pure refinement of the final model.
• mu – how aggressively variables are dropped. Larger is faster.
• eta – the gradient step size.

plot_schedule() makes the cost explicit: total work is proportional to the area under these curves.

plot_schedule(M = 1000, k = 10, niter = 500, mu = c(0, 1, 10, 100, 300, 1000))

With mu = 300, almost all of the 1000 variables are gone within 50 iterations, so the bulk of the run operates on a tiny problem.

xc <- matrix(rnorm(n * 300), n, 300) + rnorm(n)   # every column shares a factor
bc <- numeric(300); bc[c(10, 60, 150, 220, 280)] <- 1
yc <- as.numeric(xc %*% bc + rnorm(n) * 0.5)

fc <- fsa(xc, yc, k = 5, eta = 1)        # asked for 1 ...
fc$eta_final                             # ... capped to what is actually safe
#> [1] 0.48259
sort(fc$selected[[1]])
#> [1]  10  60 150 220 280

Families

Robust classification losses

For family = "binomial", loss selects the per-observation loss. The Lorenz loss is introduced by this paper and is the reason to reach for fsa() when labels are unreliable: it grows only logarithmically as the margin becomes very negative, so a handful of mislabeled points cannot dominate the fit.

u <- seq(-20, 3, length.out = 400)
L <- sapply(c("logistic", "lorenz", "svm", "hinge"),
            function(l) fsa_losses(u, l))
matplot(u, L, type = "l", lty = c(1, 1, 2, 1), lwd = 2,
        col = hcl.colors(4, "Dark 3"), xlab = "margin", ylab = "loss")
legend("topright", colnames(L), lty = c(1, 1, 2, 1), lwd = 2, bty = "n",
       col = hcl.colors(4, "Dark 3"))

yb <- as.integer(x %*% beta > 0)
flip <- sample(n, 0.15 * n); yb[flip] <- 1 - yb[flip]   # 15% wrong labels
sapply(c("logistic", "lorenz"), function(l) {
  s <- fsa(x, factor(yb), family = "binomial", loss = l, k = 5)$selected[[1]]
  length(intersect(s, c(10, 50, 100, 200, 300)))
})
#> logistic   lorenz 
#>        5        5

The ridge prior, and when it stops being optional

s adds s * ||beta||^2 to the objective (paper eq. 4). The default is 0, which is faithful to the paper and the right choice when n comfortably exceeds p. When p is much larger than n it is not: there is not enough data to pin down the surviving coefficients, and the anneal makes its keep-or-kill decisions on noisy estimates. cv.fsa() will tune s jointly with k if you pass it a vector.

data(dexter)                      # 300 observations, 20,000 features
set.seed(11)
cv <- cv.fsa(dexter$x, dexter$y, family = "binomial", loss = "lorenz",
             k = c(20, 10, 5), s = c(0, 1e-2, 1), nfolds = 5,
             type.measure = "class")
c(k = cv$k.min, s = cv$s.min)
#>  k  s 
#> 20  1

Sparse designs

Sparse matrices are handled natively and are never densified; the sparsity structure is preserved through standardization by carrying the centering analytically.

xs <- Matrix::rsparsematrix(500, 2000, density = 0.01)
bs <- numeric(2000); bs[c(7, 300, 1500)] <- 4
ys <- as.numeric(xs %*% bs + rnorm(500) * 0.2)
sort(fsa(xs, ys, k = 3)$selected[[1]])
#> [1]    7  300 1500

Coming from glmnet

vignette("glmfsa-vs-glmnet") works through where each package wins.

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.