---
title: "glmfsa and glmnet"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{glmfsa and glmnet}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7,
                      fig.height = 4.2, dpi = 110)
has_glmnet <- requireNamespace("glmnet", quietly = TRUE)
knitr::opts_chunk$set(eval = has_glmnet)
library(glmfsa)
set.seed(20)
```

```{r nog, eval = !has_glmnet, echo = FALSE, results = "asis"}
cat("**glmnet is not installed, so the comparisons below were not run.**\n")
```

Both packages fit sparse linear models. They disagree about what you should be
asked to specify, and that disagreement is the whole story.

`glmnet` minimizes `loss + lambda * penalty(beta)`. Sparsity is a *consequence*
of the penalty: you choose `lambda`, and the model size is whatever falls out.
`fsa` minimizes the loss subject to `||beta||_0 <= k`. Sparsity is the
*constraint*: you choose `k`, and get exactly that many variables.

Three things follow.

**You specify what you care about.** If your collaborator needs a 20-gene
panel, `fsa(x, y, k = 20)` gives one. With `glmnet` you search the `lambda`
path for whichever value happens to produce 20 non-zeros, and depending on the
grid there may not be one.

**No shrinkage bias.** The L1 penalty that performs the selection also pulls the
selected coefficients toward zero. This is why refitting an unpenalized model on
the selected support ("the relaxed lasso") is standard practice. FSA's
constraint does no shrinking, so its coefficients are already unbiased.

**The problem shrinks as you fit.** Eliminated variables leave the design
matrix, so late iterations operate on a small problem. Coordinate descent must
keep cycling over all the variables still in its active set.

```{r setup2}
library(glmnet)

make_data <- function(n, p, rho = 0.9, k = 10, task = "regression",
                      noise = 0) {
  z <- matrix(rnorm(n * p), n, p); x <- z
  for (j in 2:p) x[, j] <- rho * x[, j - 1] + sqrt(1 - rho^2) * z[, j]
  ti <- seq_len(k) * 10L
  b <- numeric(p); b[ti] <- 1
  if (task == "regression") {
    y <- as.numeric(x %*% b + rnorm(n))
  } else {
    y <- as.integer(x %*% b > 0)
    if (noise > 0) { fl <- sample(n, round(noise * n))
                     y[fl] <- sample(0:1, length(fl), TRUE) }
  }
  list(x = x, y = y, true = ti)
}
```

## 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`.

```{r accuracy}
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)
```

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

```{r 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))
```

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.

```{r 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)
```

## 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:

```{r timing}
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)
```

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:

```{r fast}
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])
c(default = length(intersect(f1$selected[[1]], d$true)),
  fast    = length(intersect(f2$selected[[1]], d$true)))
```

## 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.

```{r real, eval = has_glmnet}
data(dexter)
dim(dexter$x)
sum(Matrix::colSums(dexter$x != 0) == 0)   # columns with no support at all
```

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

```{r realcv, eval = has_glmnet}
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)

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]))
})
```

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

| Dataset | k | glmfsa (logistic) | glmfsa (Lorenz) | lasso | elastic net |
|---|---|---|---|---|---|
| Gisette | 10 | 9.5 | **7.9** | 15.7 | 16.3 |
| Gisette | 20 | 7.0 | **5.8** | 13.4 | 14.7 |
| Gisette | 50 | **3.5** | 3.7 | 7.4 | 10.5 |
| Dexter | 5 | 25.3 | **18.3** | 23.0 | 26.3 |
| Dexter | 10 | 18.3 | **13.7** | 16.3 | 20.3 |
| Dexter | 50 | 11.7 | 14.7 | **8.7** | 9.7 |

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.