Terminating forward selection

A T-solver is a forward selection algorithm on the enlarged design \([\boldsymbol X, \boldsymbol D]\) of predictors and dummies that stops as soon as T_stop dummies are active. The T-Rex selectors run one inside every random experiment, and every T-solver is available as a class of its own, for a single terminated path, for a study of the solvers, or as a building block. The terminated form of each algorithm follows the T-Rex paper [1], the greedy family the terminating greedy pursuits [2].

Family	Classes	Algorithm
LARS	`TLARS_Solver`, `TLASSO_Solver`, `TSTEPWISE_Solver`	least angle regression and its lasso and stepwise modifications [3], [4]
LARS	`TSTAGEWISE_Solver`	forward stagewise [5], [6]
LARS	`TENET_Solver`, `TENETAug_Solver`	the elastic net [7], Gram based and row augmented
LARS	`TIENET_Solver`, `TIENETAug_Solver`	the informed elastic net [8], pathwise and row augmented
LARS	`TIKNET_Solver`	the Tikhonov net, pathwise
OMP	`TOMP_Solver`, `TOOLS_Solver`	orthogonal matching pursuit and orthogonal least squares [9]–[11]
OMP	`TGP_Solver`, `TACGP_Solver`	gradient and approximate conjugate gradient pursuit [12]
OMP	`TMP_Solver`, `TNCGMP_Solver`	matching pursuit [13] and non-negative constrained generalized matching pursuit [14]
AFS	`TAFS_Solver`	adaptive forward stepwise [15]
CD	`TCDLASSO_Solver`, `TCDENET_Solver`, `TCDIENET_Solver`, `TCDTIKNET_Solver`	pathwise coordinate descent [16], [17] with exact minimizers at every dummy crossing

The path solvers (LARS and CD families) walk a regularization path, the greedy solvers (OMP and AFS families) add one variable per step. The informed solvers take a groups vector, the Tikhonov solvers a penalty matrix tikhonov_K.

One terminated path

set.seed(31)
n <- 100
p <- 50
L <- 50
true_support <- c(2, 9, 17, 33)
beta <- numeric(p)
beta[true_support] <- c(2, -2, 1.5, -1.5)
X <- matrix(rnorm(n * p), n, p)
D <- matrix(rnorm(n * L), n, L)
y <- as.numeric(X %*% beta + rnorm(n))

A solver is constructed on X, D, y and an options list, and execute_step() runs the path until T_stop dummies are active. The path can be continued by calling execute_step() again with a larger T_stop.

solver <- TLARS_Solver$new(X, D, y)
solver$execute_step(T_stop = 3)
solver$get_num_steps()
#> [1] 11
solver$get_active_predictor_indices()
#> [1]  9  2 33 17 16 23 49 47
solver$get_active_dummy_indices()
#> [1]  2  3 10

Indices are 1-based and the dummies continue after the predictors, so the active set of the enlarged design runs from 1 to \(p + L\).

solver$get_actives()
#> [1]  9  2 33 17 16 23 49 47
solver$get_dummy_start_index()
#> [1] 51
solver$get_actions()[1:4]
#> [[1]]
#> [1] 9
#> 
#> [[2]]
#> [1] 2
#> 
#> [[3]]
#> [1] 33
#> 
#> [[4]]
#> [1] 17

get_actions() records the path as signed codes per step, +j for an addition and -j for a removal. get_beta() is the coefficient vector of the last step and get_beta_path() the whole path, one column per step with the all-zero start.

dim(solver$get_beta_path())
#> [1] 100  12
round(solver$get_beta()[true_support], 3)
#> [1]  1.715 -1.872  1.079 -1.184
path_diag <- data.frame(lambda = solver$get_lambda(), rss = solver$get_rss(),
                        r2 = solver$get_r2(), cp = solver$get_cp(), dof = solver$get_dof())
round(head(path_diag, 5), 3)
#>   lambda      rss    r2      cp dof
#> 1 23.844 1066.332 0.125 776.431   2
#> 2 20.403  558.071 0.542 362.592   3
#> 3 10.417  490.974 0.597 309.696   4
#> 4  9.238  133.327 0.891  19.083   5
#> 5  2.896  130.379 0.893  18.671   6
solver$restore()    # give X and D back before the next solver claims them

The path as a picture

The coefficient paths of a terminated run show what the solver does. A predictor enters at a step and its coefficient grows from there, and a dummy that enters marks one of the T_stop crossings that end the path. One solver runs to T_stop = 1, keeps a snapshot of its path, and continues to T_stop = 3, so the second panel extends the first.

library(ggplot2)

snapshot <- function(solver, label) {
  path <- solver$get_beta_path()                            # (p + L) x (steps + 1)
  entering <- abs(vapply(solver$get_actions(), function(a) a[1], numeric(1)))
  moved <- which(rowSums(path != 0) > 0)
  kind <- ifelse(moved > p, "dummy", ifelse(moved %in% true_support, "true predictor", "other predictor"))
  lines <- data.frame(panel = label, variable = rep(moved, each = ncol(path)),
                      step = rep(seq_len(ncol(path)) - 1, times = length(moved)),
                      coefficient = as.vector(t(path[moved, , drop = FALSE])),
                      kind = rep(kind, each = ncol(path)))
  marks <- data.frame(panel = label, step = seq_along(entering), dummy = entering > p,
                      label = ifelse(entering > p, paste0("D", entering - p), entering))
  list(lines = lines, marks = marks)
}

s <- TLARS_Solver$new(X, D, y)
s$execute_step(T_stop = 1)
first <- snapshot(s, "T_stop = 1")
s$execute_step(T_stop = 3)
second <- snapshot(s, "T_stop = 3, continued")
s$restore()

lines <- rbind(first$lines, second$lines)
marks <- rbind(first$marks, second$marks)
lines$panel <- factor(lines$panel, levels = c("T_stop = 1", "T_stop = 3, continued"))
marks$panel <- factor(marks$panel, levels = levels(lines$panel))
inclusions <- marks[marks$dummy, ]
inclusions$kind <- "dummy inclusion"
kinds <- c("true predictor", "other predictor", "dummy", "dummy inclusion")
ggplot(lines, aes(step, coefficient, group = variable, colour = kind, linetype = kind,
                  linewidth = kind)) +
  geom_vline(data = inclusions, aes(xintercept = step, colour = kind, linetype = kind,
                                    linewidth = kind), inherit.aes = FALSE) +
  geom_line() +
  geom_text(data = marks, aes(x = step, y = Inf, label = label), inherit.aes = FALSE,
            vjust = 1.4, size = 2.4, colour = "grey30") +
  facet_wrap(~ panel, scales = "free_x") +
  scale_colour_manual(values = c("true predictor" = "#336C68", "other predictor" = "grey45",
                                 "dummy" = "#D1495B", "dummy inclusion" = "#D1495B"),
                      breaks = kinds) +
  scale_linetype_manual(values = c("true predictor" = "solid", "other predictor" = "solid",
                                   "dummy" = "dashed", "dummy inclusion" = "dotted"),
                        breaks = kinds) +
  scale_linewidth_manual(values = c("true predictor" = 0.9, "other predictor" = 0.5,
                                    "dummy" = 0.6, "dummy inclusion" = 0.4), breaks = kinds) +
  labs(x = "Step", y = "Coefficient", colour = NULL, linetype = NULL, linewidth = NULL) +
  theme_bw(base_size = 11) +
  theme(legend.position = "right", panel.grid.minor = element_blank(),
        legend.background = element_rect(fill = "white", colour = "grey60"))

The true predictors enter first and keep growing, further predictors follow with small coefficients, and the dummies end the path. The second panel repeats the first up to its last step and continues from there, because the continuation walks the stored path on. The label over each step names the variable that enters.

The options

tsolver_options() holds every input except the data, the same list for every solver. A solver reads the fields it needs and ignores the rest.

Field	Meaning
`normalize`, `intercept`, `scaling_mode`	center and scale the columns of `X` and `D` in place (`"L2"` or `"ZSCORE"`) and center `y`
`precentered`	the caller centered the data: the active set is capped at \(n - 1\) and the intercept counts in the degrees of freedom
`tol`	the numerical tolerance of the comparisons (machine epsilon by default)
`tie_seed`	the seed of the tie-break engine (0 by default, so exact ties resolve reproducibly)
`n_threads`	the thread budget of the solver’s own passes (1 serial, at least 1, and a pass splits only above about a million entries of the design)
`exch_tie_alpha`, `exch_tie_floor`	the exchangeable-tie band of the greedy solvers
`lambda2`	the quadratic penalty of the elastic nets (0 by default)
`rho_afs`, `ncgmp_variant`, `tenet_aug_use_lars`	solver-specific switches
`verbose`	print the steps

en <- TENET_Solver$new(X, D, y, options = tsolver_options(lambda2 = 0.5))
en$execute_step(T_stop = 3)
en$get_active_predictor_indices()
#> [1]  9  2 33 17 16 29 19 25
en$get_num_removals()
#> [1] 0
en$restore()

Writers and readers

With normalize = TRUE or intercept = TRUE (the defaults) a solver is a writer. It centers and scales the columns of X and D in place, holds the statistics, and writes the original values back when restore() is called or when the object is garbage collected. A second writer on a claimed matrix is an error, because a second in-place transformation would move the data by rounding. y is never written.

X_before <- X + 0
w <- TLASSO_Solver$new(X, D, y)
max(abs(X - X_before))      # the scaled columns
#> [1] 4.137846
w$restore()
max(abs(X - X_before))      # the original values
#> [1] 1.110223e-16

A pipeline that runs several solvers on one design scales the data once in front with a scaler and constructs readers, solvers with normalize = FALSE and intercept = FALSE that never write.

sx <- MLLpNormScaler$new(norm_type = 2, center = TRUE)
invisible(sx$fit_transform_inplace(X))
sd <- MLLpNormScaler$new(norm_type = 2, center = TRUE)
invisible(sd$fit_transform_inplace(D))
y_c <- y - mean(y)
reader_opts <- tsolver_options(normalize = FALSE, intercept = FALSE, precentered = TRUE)

With the data pre-scaled, every family runs on the same design.

constructors <- list(
  TLARS = TLARS_Solver, TLASSO = TLASSO_Solver, TSTEPWISE = TSTEPWISE_Solver,
  TSTAGEWISE = TSTAGEWISE_Solver, TENET = TENET_Solver, TENETAug = TENETAug_Solver,
  TOMP = TOMP_Solver, TOOLS = TOOLS_Solver, TGP = TGP_Solver, TACGP = TACGP_Solver,
  TMP = TMP_Solver, TNCGMP = TNCGMP_Solver, TAFS = TAFS_Solver,
  TCDLASSO = TCDLASSO_Solver, TCDENET = TCDENET_Solver)
rows <- lapply(names(constructors), function(nm) {
  s <- constructors[[nm]]$new(X, D, y_c, options = reader_opts)
  s$execute_step(T_stop = 3)
  data.frame(solver = nm, steps = s$get_num_steps(),
             predictors = paste(s$get_active_predictor_indices(), collapse = " "))
})
do.call(rbind, rows)
#>        solver steps            predictors
#> 1       TLARS    11 9 2 33 17 16 23 49 47
#> 2      TLASSO    11 9 2 33 17 16 23 49 47
#> 3   TSTEPWISE     7             9 2 33 17
#> 4  TSTAGEWISE    11 9 2 33 17 16 23 49 47
#> 5       TENET    11 9 2 33 17 16 23 49 47
#> 6    TENETAug    11 9 2 33 17 16 23 49 47
#> 7        TOMP     7             9 2 33 17
#> 8       TOOLS     7             9 2 33 17
#> 9         TGP    10        9 2 33 17 47 3
#> 10      TACGP     9        9 2 33 17 47 3
#> 11        TMP     9             9 2 33 17
#> 12     TNCGMP     9             9 2 33 17
#> 13       TAFS     7             9 2 33 17
#> 14   TCDLASSO     3 2 9 16 17 23 33 47 49
#> 15    TCDENET     3 2 9 16 17 23 33 47 49

The informed solvers and the Tikhonov solvers take their structure as an extra argument.

groups <- rep(1:5, each = 10)
ien <- TIENET_Solver$new(X, D, y_c, groups = groups,
                         options = tsolver_options(normalize = FALSE, intercept = FALSE,
                                                   precentered = TRUE, lambda2 = 0.5))
ien$execute_step(T_stop = 3)
ien$get_active_predictor_indices()
#> [1]  9  2 33 17 16 23
tik <- TCDTIKNET_Solver$new(X, D, y_c, tikhonov_K = diag(p),
                            options = tsolver_options(normalize = FALSE, intercept = FALSE,
                                                      precentered = TRUE, lambda2 = 0.5))
tik$execute_step(T_stop = 3)
tik$get_active_predictor_indices()
#> [1]  2  9 16 17 19 25 29 33

Saving and loading

save() serializes the state of a solver to a file and load() restores it into a solver constructed on the same data in the same state. The loaded solver is a reader and continues the path where it was saved.

f <- tempfile(fileext = ".bin")
a <- TLARS_Solver$new(X, D, y_c, options = reader_opts)
a$execute_step(T_stop = 2)
a$save(f)
b <- TLARS_Solver$new(X, D, y_c, options = reader_opts)
b$load(f)
b$get_num_steps()
#> [1] 8
b$execute_step(T_stop = 4)
c <- TLARS_Solver$new(X, D, y_c, options = reader_opts)
c$execute_step(T_stop = 4)
identical(b$get_actives(), c$get_actives())
#> [1] TRUE
unlink(f)

Memory

get_memory_bytes() reports the heap bytes a solver holds (its copy of y, the residual, the correlations, the Cholesky factor, the coefficient path). X and D are borrowed and not counted.

a$get_memory_bytes()
#> [1] 19696

Examples

The TRexSelector_Examples repository holds runnable demos and Monte Carlo studies of every T-solver and five performance studies, with their results and a README per demo, in R, Python and C++.

Bibliography

[1]

J. Machkour, M. Muma, and D. P. Palomar, “The terminating-random experiments selector: Fast high-dimensional variable selection with false discovery rate control,” Signal Process., vol. 231, 2025, doi: 10.1016/j.sigpro.2025.109894.

[2]

F. Scheidt and M. Muma, “Terminating greedy pursuits for high-dimensional false discovery rate control,” in Proc. Eur. Signal process. Conf. (EUSIPCO), 2026, pp. 2511–2515.

[3]

B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani, “Least angle regression,” Ann. Statist., vol. 32, no. 2, pp. 407–499, 2004, doi: 10.1214/009053604000000067.

[4]

R. Tibshirani, “Regression shrinkage and selection via the lasso,” J. R. Stat. Soc. Ser. B Stat. Methodol., vol. 58, no. 1, pp. 267–288, 1996, doi: 10.1111/j.2517-6161.1996.tb02080.x.

[5]

T. Hastie, J. Taylor, R. Tibshirani, and G. Walther, “Forward stagewise regression and the monotone lasso,” Electron. J. Stat., vol. 1, pp. 1–29, 2007, doi: 10.1214/07-EJS004.

[6]

R. J. Tibshirani, “A general framework for fast stagewise algorithms,” J. Mach. Learn. Res., vol. 16, pp. 2543–2588, 2015.

[7]

H. Zou and T. Hastie, “Regularization and variable selection via the elastic net,” J. R. Stat. Soc. Ser. B Stat. Methodol., vol. 67, no. 2, pp. 301–320, 2005, doi: 10.1111/j.1467-9868.2005.00503.x.

[8]

J. Machkour, M. Muma, and D. P. Palomar, “The informed elastic net for fast grouped variable selection and FDR control in genomics research,” in Proc. IEEE 9th int. Workshop comput. Adv. Multi-sensor adapt. Process. (CAMSAP), 2023, pp. 466–470. doi: 10.1109/CAMSAP58249.2023.10403489.

[9]

Y. C. Pati, R. Rezaiifar, and P. S. Krishnaprasad, “Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition,” in Proc. 27th asilomar conf. Signals, syst., comput., 1993, pp. 40–44. doi: 10.1109/ACSSC.1993.342465.

[10]

S. Chen, S. A. Billings, and W. Luo, “Orthogonal least squares methods and their application to non-linear system identification,” Int. J. Control, vol. 50, no. 5, pp. 1873–1896, 1989, doi: 10.1080/00207178908953472.

[11]

T. Blumensath and M. E. Davies, “On the difference between orthogonal matching pursuit and orthogonal least squares,” Mar. 2007.

[12]

T. Blumensath and M. E. Davies, “Gradient pursuits,” IEEE Trans. Signal Process., vol. 56, no. 6, pp. 2370–2382, 2008, doi: 10.1109/TSP.2007.916124.

[13]

S. G. Mallat and Z. Zhang, “Matching pursuits with time-frequency dictionaries,” IEEE Trans. Signal Process., vol. 41, no. 12, pp. 3397–3415, 1993, doi: 10.1109/78.258082.

[14]

F. Locatello, R. Khanna, M. Tschannen, and M. Jaggi, “A unified optimization view on generalized matching pursuit and Frank-Wolfe,” in Proc. 20th int. Conf. Artif. Intell. Statist. (AISTATS), 2017, vol. 54.

[15]

I. Zhang and R. Tibshirani, “Adaptive forward stepwise: A method for high sparsity regression,” J. Mach. Learn. Res., vol. 27, pp. 1–24, 2026, Available: http://jmlr.org/papers/v27/25-0151.html

[16]

J. Friedman, T. Hastie, H. Höfling, and R. Tibshirani, “Pathwise coordinate optimization,” Ann. Appl. Stat., vol. 1, no. 2, pp. 302–332, 2007, doi: 10.1214/07-AOAS131.

[17]

J. Friedman, T. Hastie, and R. Tibshirani, “Regularization paths for generalized linear models via coordinate descent,” J. Stat. Softw., vol. 33, no. 1, pp. 1–22, 2010, doi: 10.18637/jss.v033.i01.

The T-solvers