The T-solvers

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 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; 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] 4.440892e-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

References

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