Package {glmfsa}

glmfsa: Feature Selection with Annealing

Description

Fits sparse generalized linear models under an exact constraint on the number of non-zero coefficients, using the Feature Selection with Annealing algorithm of Barbu, She, Ding and Gramajo (2017).

Details

The algorithm alternates a gradient step on the loss with a step that removes the lowest-magnitude coefficients according to an annealing schedule. Removed variables leave the working design matrix altogether, so a fit gets cheaper as it proceeds – the property that makes the method practical when the number of features is very large.

Start with [fsa()] to fit, [cv.fsa()] to choose the model size, and 'vignette("glmfsa-intro")' for a tour. Users coming from **glmnet** will find 'vignette("glmfsa-vs-glmnet")' the fastest way in.

Author(s)

Maintainer: Adrian Barbu abarbu@fsu.edu [copyright holder]

Authors:

• Adrian Barbu abarbu@fsu.edu [copyright holder]

• Yiyuan She

Other contributors:

• Liangjing Ding [contributor]

• Gary Gramajo [contributor]

References

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. doi:10.1109/TPAMI.2016.2544315

See Also

Useful links:

• https://github.com/barbua/FSA

• Report bugs at https://github.com/barbua/FSA/issues

Extract coefficients from a fitted FSA model

Description

Extract coefficients from a fitted FSA model

Usage

## S3 method for class 'fsa'
coef(object, k = NULL, ...)

Arguments

Value

A sparse matrix with 'nvars + 1' rows (the intercept first), or a list of such matrices for '"multinomial"' and '"mgaussian"' fits.

Examples

set.seed(1)
x <- matrix(rnorm(100 * 40), 100, 40)
y <- x[, 3] - 2 * x[, 17] + rnorm(100)
coef(fsa(x, y, k = 2))

Cross-validation for FSA

Description

Chooses the model size 'k' by K-fold cross-validation. The interface and the returned object deliberately mirror [glmnet::cv.glmnet()], so existing analysis code needs only the names changed.

Usage

cv.fsa(
  x,
  y,
  family = "gaussian",
  k = NULL,
  s = 0,
  type.measure = c("default", "deviance", "mse", "mae", "class", "auc", "C"),
  nfolds = 10L,
  foldid = NULL,
  weights = NULL,
  parallel = FALSE,
  keep = FALSE,
  ...
)

Arguments

Value

An object of class '"cv.fsa"' with components 'k', 'cvm', 'cvsd', 'cvup', 'cvlo', 'nzero', 'k.min', 'k.1se', 's.min', 'index', 'name', 'fsa.fit' and, if 'keep = TRUE', 'fit.preval' and 'foldid'. When 's' was a vector, 'cvm.s' holds the full 'length(k) x length(s)' surface.

See Also

[fsa()], [plot.cv.fsa()]

Examples

set.seed(1)
n <- 200; p <- 300
x <- matrix(rnorm(n * p), n, p)
b <- numeric(p); b[c(7, 23, 61, 140)] <- 1.5
y <- as.numeric(x %*% b + rnorm(n))

cvfit <- cv.fsa(x, y, k = c(20, 10, 6, 4, 2, 1), nfolds = 5)
cvfit$k.min
coef(cvfit, k = "k.min")[1:5, ]

Methods for cross-validated FSA fits

Description

Methods for cross-validated FSA fits

Usage

## S3 method for class 'cv.fsa'
predict(object, newx, k = "k.min", ...)

## S3 method for class 'cv.fsa'
coef(object, k = "k.min", ...)

## S3 method for class 'cv.fsa'
print(x, digits = max(3, getOption("digits") - 3), ...)

## S3 method for class 'cv.fsa'
plot(x, sign.lambda = 1, ...)

Arguments

Value

'predict' returns predictions at the chosen 'k'; 'coef' returns the coefficient vector; 'plot' draws the cross-validation curve and returns 'NULL' invisibly; 'print' returns 'x' invisibly.

Examples

set.seed(1)
x <- matrix(rnorm(150 * 60), 150, 60)
y <- as.numeric(x[, c(3, 9)] %*% c(2, -2) + rnorm(150))
cvfit <- cv.fsa(x, y, k = c(10, 5, 2, 1), nfolds = 5)
plot(cvfit)
head(predict(cvfit, x, k = "k.1se"))

Deviance of a fitted FSA model

Description

Deviance of a fitted FSA model

Usage

## S3 method for class 'fsa'
deviance(object, ...)

Arguments

Details

Returns '2 * sum(weights) * loss', where 'loss' is the weighted mean training loss. For 'family = "gaussian"' this is the residual sum of squares; for '"binomial"' with the logistic loss, '"poisson"' and '"multinomial"' it is the usual deviance up to the saturated-model constant; for '"cox"' it is minus twice the Breslow partial log-likelihood.

Value

A numeric vector with one deviance per fitted value of 'k'.

Examples

set.seed(1)
x <- matrix(rnorm(100 * 30), 100, 30)
deviance(fsa(x, x[, 1] + rnorm(100), k = c(5, 2)))

Dexter: text categorization with 20,000 features and 300 observations

Description

One of the five datasets of the NIPS 2003 Feature Selection Challenge. The task is to decide whether a Reuters newswire article is about corporate acquisitions. Documents are represented as a bag of words, and the feature set was deliberately padded: of the 20,000 columns, 10,053 are **probes** – randomly permuted features with no relationship to the label, included so that a method's false-discovery behavior can be measured. The design is very sparse (0.47 which makes it a natural test for cardinality-constrained fitting.

Usage

dexter

Format

A list with four elements:

x

300 x 20000 sparse 'dgCMatrix' of word counts (training).

y

Factor of length 300 with levels '"neg"' and '"pos"', balanced 150/150.

xvalid

300 x 20000 sparse 'dgCMatrix' (validation).

yvalid

Factor of length 300, balanced 150/150.

Details

Test-set labels were never released, so the published validation split is supplied here as the held-out set.

Source

UCI Machine Learning Repository, https://archive.ics.uci.edu/dataset/168/dexter. Licensed CC BY 4.0.

References

Guyon, I., Gunn, S., Ben-Hur, A. and Dror, G. (2004). Result analysis of the NIPS 2003 feature selection challenge. Advances in Neural Information Processing Systems 17.

Examples

data(dexter)
dim(dexter$x)
fit <- fsa(dexter$x, dexter$y, family = "binomial", k = 20)
mean(predict(fit, dexter$xvalid, k = 20, type = "class") !=
       as.character(dexter$yvalid))

Download the Gisette dataset

Description

Gisette is another dataset from the NIPS 2003 Feature Selection Challenge: separating the handwritten digits 4 and 9, with 5,000 features of which 2,500 are probes (random features carrying no signal). At 6,000 training observations it is too large to ship inside an R package, so it is downloaded from the UCI Machine Learning Repository on first use and cached.

Usage

fetch_gisette(cache_dir = tools::R_user_dir("glmfsa", "cache"), quiet = FALSE)

Arguments

Value

A list shaped like [dexter]: 'x' (6000 x 5000 numeric matrix), 'y' (factor, 3000/3000), 'xvalid' (1000 x 5000), 'yvalid' (factor, 500/500).

Source

UCI Machine Learning Repository, https://archive.ics.uci.edu/dataset/170/gisette. Licensed CC BY 4.0.

References

Guyon, I., Gunn, S., Ben-Hur, A. and Dror, G. (2004). Result analysis of the NIPS 2003 feature selection challenge. Advances in Neural Information Processing Systems 17.

Examples

## Not run during package checks: the first call downloads ~21 MB from UCI.
## Not run: 
g <- fetch_gisette()
fit <- fsa(g$x, g$y, family = "binomial", k = 50)
mean(predict(fit, g$xvalid, k = 50, type = "class") !=
       as.character(g$yvalid))

## End(Not run)

Fit a sparse generalized linear model by Feature Selection with Annealing

Description

'fsa()' fits a linear model under an exact constraint on the number of non-zero coefficients, using the Feature Selection with Annealing algorithm of Barbu et al. (2017). It alternates a gradient step on the loss with a step that removes the lowest-magnitude coefficients according to an annealing schedule, so that the working design matrix physically shrinks as the fit proceeds.

Usage

fsa(
  x,
  y,
  family = c("gaussian", "binomial", "poisson", "multinomial", "mgaussian", "cox"),
  k = NULL,
  nk = 10L,
  loss = NULL,
  s = 0,
  niter = 500L,
  mu = 300,
  eta = NULL,
  mom = 0.9,
  nesterov = FALSE,
  minibatch = 0L,
  schedule = c("inverse", "exponential", "linear"),
  perc_exp = 0.8,
  standardize = TRUE,
  intercept = TRUE,
  weights = NULL,
  offset = NULL,
  include = NULL,
  exclude = NULL,
  nwarm = 10L,
  adapt = TRUE,
  cap = TRUE,
  balance = TRUE,
  keep_trace = FALSE,
  verbose = 0L
)

Arguments

Details

Where [glmnet::glmnet()] asks for a penalty 'lambda' and gives you whatever model size that implies, 'fsa()' asks for the model size 'k' directly. The surviving coefficients are not shrunk toward zero by a sparsity penalty, so the "fit the lasso to select, then refit on the selected variables to de-bias" two-step is unnecessary.

Value

An object of class '"fsa"', a list with components:

Choosing the tuning parameters

The defaults follow the paper and are deliberately not aggressive. 'niter' controls how long the anneal runs, 'mu' how fast variables are dropped (large 'mu' drops most variables almost immediately and is much cheaper; 'mu = 0' drops them linearly), and 'eta' is the gradient step size. Figure 6 of the paper shows prediction performance is flat over roughly an order of magnitude in each of the three, so tuning them is rarely the best use of effort – tune 'k' instead, with [cv.fsa()].

The 'k' path

Passing a decreasing vector to 'k' returns one model per entry, computed in a single annealed pass: the run anneals down to 'k[1]', records the fit, then continues annealing to 'k[2]', and so on. The supports are therefore nested and the whole path costs little more than the single fit at 'min(k)'.

References

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. doi:10.1109/TPAMI.2016.2544315

See Also

[cv.fsa()] to choose 'k', [predict.fsa()], [plot.fsa()], [plot_trace()]

Examples

set.seed(1)
n <- 200; p <- 500
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)
which(coef(fit)[-1] != 0)        # the five true variables

## a whole path of model sizes in one pass
fit2 <- fsa(x, y, k = c(20, 10, 5, 2))
fit2

Per-observation classification losses used by FSA

Description

Evaluates the margin-based losses available for 'family = "binomial"', as defined in equations (6)-(8) of Barbu et al. (2017). Useful for seeing why the Lorenz loss is the robust choice: it is the only one of the four whose value stays bounded in growth as the margin becomes very negative, so a handful of mislabeled points cannot dominate the fit.

Usage

fsa_losses(u, loss = c("logistic", "lorenz", "svm", "hinge"), h = 0.1)

Arguments

Value

A numeric vector of loss values.

Examples

u <- seq(-6, 3, length.out = 200)
matplot(u, sapply(c("logistic", "lorenz", "svm", "hinge"),
                  function(l) fsa_losses(u, l)),
        type = "l", lty = 1, ylab = "loss", xlab = "margin")

Annealing schedule (internal)

Description

Exposed so the schedule used by the C++ engine can be plotted and tested from R without re-implementing it.

Usage

fsa_schedule_cpp(e, niter, kStar, M0, mu, sched, perc_exp)

Arguments

Value

integer vector of retained-variable counts.

Plot the coefficient path of a fitted FSA model

Description

Draws each selected coefficient against the model size 'k'. Because FSA supports along a path are nested, a coefficient enters the plot at the largest 'k' where it survives and continues to the sparsest model that retains it.

Usage

## S3 method for class 'fsa'
plot(x, xvar = c("k", "df"), label = FALSE, nlabel = 10L, response = 1L, ...)

Arguments

Value

'NULL', invisibly. Called for its side effect of drawing the plot.

Examples

set.seed(1)
x <- matrix(rnorm(200 * 100), 200, 100)
y <- as.numeric(x[, c(5, 20, 60)] %*% c(2, -1.5, 1) + rnorm(200))
plot(fsa(x, y, k = c(30, 20, 15, 10, 6, 4, 3, 2, 1)))

Plot an annealing schedule

Description

Draws the number of surviving variables against the iteration index for the schedules implemented by the package – the figure that makes the cost of a fit obvious, since total work is proportional to the area under the curve (paper Fig. 3 and Table 1).

Usage

plot_schedule(
  M = 1000,
  k = 10,
  niter = 500,
  mu = c(0, 10, 100, 300),
  schedule = c("inverse", "exponential", "linear"),
  perc_exp = 0.8,
  ...
)

Arguments

Value

Invisibly, a matrix of retained-variable counts with one column per value of 'mu'.

Examples

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

Plot the annealing trajectory

Description

Reproduces Figure 2 of Barbu et al. (2017): every coefficient's value against the iteration number, with curves terminating as the annealing schedule eliminates their variables. It is the most direct picture of what the algorithm does, and a quick diagnostic – if coefficients are still moving sharply at the final iteration, 'niter' is too small or 'eta' too large.

Usage

plot_trace(object, highlight = NULL, show_size = TRUE, ...)

Arguments

Value

'NULL', invisibly. Called for its side effect of drawing the plot.

Examples

set.seed(1)
x <- matrix(rnorm(200 * 200), 200, 200)
y <- as.numeric(x[, c(10, 50, 120)] %*% c(2, -2, 1.5) + rnorm(200))
fit <- fsa(x, y, k = 3, niter = 120, keep_trace = TRUE)
plot_trace(fit, highlight = c(10, 50, 120))

Predict from a fitted FSA model

Description

Predict from a fitted FSA model

Usage

## S3 method for class 'fsa'
predict(
  object,
  newx,
  k = NULL,
  type = c("link", "response", "class", "coefficients", "nonzero"),
  newoffset = NULL,
  ...
)

## S3 method for class 'fsa'
fitted(object, ...)

Arguments

Details

For 'family = "binomial"' with a margin-based 'loss' ('"lorenz"', '"svm"', '"hinge"'), 'type = "response"' returns 'plogis()' of the margin. That is monotone in the decision value – so it ranks observations correctly and is fine for ROC or AUC – but it is not a calibrated probability, because those losses do not estimate one. Use 'loss = "logistic"' if you need probabilities.

Value

A matrix of predictions with one column per value of 'k', or an array 'nobs x K x length(k)' for multi-response families.

Examples

set.seed(1)
x <- matrix(rnorm(200 * 50), 200, 50)
y <- factor(x[, 5] + x[, 20] > 0)
fit <- fsa(x, y, family = "binomial", k = c(5, 2))
head(predict(fit, x, k = 2, type = "response"))
predict(fit, type = "nonzero")

Print a fitted FSA model

Description

Print a fitted FSA model

Usage

## S3 method for class 'fsa'
print(x, digits = max(3, getOption("digits") - 3), ...)

## S3 method for class 'fsa'
summary(object, digits = max(3, getOption("digits") - 3), ...)

Arguments

Value

'x', invisibly. Called for the printed summary of the 'k' path.

Examples

set.seed(1)
x <- matrix(rnorm(100 * 30), 100, 30)
print(fsa(x, x[, 1] + rnorm(100), k = c(5, 2)))