---
title: "Getting started with glmfsa"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting started with glmfsa}
  %\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)
library(glmfsa)
set.seed(1)
```

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

```{r first}
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
which(coef(fit)[-1] != 0)
```

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.

```{r coefs}
round(coef(fit)[c(1, 1 + c(10, 50, 100, 200, 300)), ], 3)
```

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.

```{r trace, fig.height = 4.5}
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`.

```{r path}
path <- fsa(x, y, k = c(50, 30, 20, 12, 8, 5, 3, 2, 1))
path
plot(path, label = TRUE)
```

## Choosing k

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

```{r cv, fig.height = 4.2}
cvfit <- cv.fsa(x, y, k = c(30, 20, 12, 8, 5, 3, 2, 1), nfolds = 5)
cvfit
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.

```{r cvpred}
head(predict(cvfit, newx = x, k = "k.1se"))
```

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

```{r sched, fig.height = 4.2}
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.

### You usually do not need to touch `eta`

By default `fsa()` caps the step size at the largest value for which the
gradient step is guaranteed to descend, estimated from the largest eigenvalue
of the scaled Gram matrix. This matters on correlated designs, where a fixed
step large enough to make progress in most directions is far too large along
the dominant one. The step actually used is reported:

```{r eta}
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
sort(fc$selected[[1]])
```

## Families

| `family` | Response | Notes |
|---|---|---|
| `"gaussian"` | numeric | squared-error loss |
| `"binomial"` | 2-level factor | four losses, see below |
| `"poisson"` | counts | |
| `"multinomial"` | K-level factor | variables selected for all classes jointly |
| `"mgaussian"` | numeric matrix | variables selected for all responses jointly |
| `"cox"` | `Surv(time, status)` | Breslow handling of ties, no intercept |

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

```{r losses, fig.height = 4}
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"))
```

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

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

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

## Sparse designs

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

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

## Coming from glmnet

| glmnet | fsa |
|---|---|
| `glmnet(x, y, family = "binomial")` | `fsa(x, y, family = "binomial")` |
| `lambda` | `k` -- the model size itself |
| `cv.glmnet(x, y)` | `cv.fsa(x, y)` |
| `lambda.min`, `lambda.1se` | `k.min`, `k.1se` |
| `coef(fit, s = "lambda.min")` | `coef(cvfit, k = "k.min")` |
| `predict(fit, newx, type = "response")` | same |
| `penalty.factor = 0` for variable `j` | `include = j` |
| `exclude = j` | same |
| `relax = TRUE` to debias | not needed |

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