Fit generalized linear model with $L_0$ or grouped row-$L_0$ regularization subject to optional nonnegativity constraints

Fits a generalized linear model with an $L_0$ penalty or, for multichannel responses, a grouped row-$L_0$ penalty, equivalent to an $L_{2,0}$ penalty, to approximate Best Subset Selection (for univariate data) or Best Subset of Groups Selection (for multivariate data) subject to optional nonnegativity constraints on model coefficients. The objective functions that can be approximated are :

for univariate data, using an $L_0$ norm penalty to achieve Best Subset Selection:

$$-2\,LL(\beta) + \lambda \|\beta\|_0$$

where $LL$ is the log-likelihood of the GLM model and $\|\beta\|_0$ is the number of non-zero elements in coefficient vector beta, i.e. the $L_0$ norm of beta.
for multivariate data, using independent elementwise $L_0$ penalized fits on each outcome:

$$-2\,LL(B) + \lambda \sum_{j=1}^{p}\sum_{k=1}^{m} I(B_{jk} \neq 0)$$
for multivariate data, using grouped variable selection in multitask learning settings to achieve Best Subset of Groups Selection:

$$-2\,LL(B) + \lambda \sum_{j=1}^{p} I(\|B_{j\cdot}\|_2 > 0)$$

Here $B$ is the $p \times m$ multichannel coefficient matrix and the grouped penalty counts active rows, i.e. predictors with a non-zero coefficient in at least one channel. Any positive within-row norm gives the same row support, so this is equivalently group-$L_0$, row-$L_0$, or $L_{2,0}$ regularization.

In all three cases, nonnegative = TRUE restricts the coefficients to be nonnegative, including any intercept, except that binomial and quasibinomial intercepts remain unrestricted:

$$\beta_j \geq 0 \quad\hbox{or}\quad B_{jk} \geq 0.$$ With nonnegative = FALSE, coefficients may have either sign.

These objectives are approximated via an iterative adaptive ridge regression procedure, see Details. The univariate and independent multichannel elementwise $L_0$ penalties can be used to approximate Best Subset Selection, while for multivariate (multitask/multichannel) data, the group-$L_0$/row-$L_0$/$L_{2,0}$ penalty can be used to approximate Best Subset of Groups Selection (via option grouped=TRUE). Fits on multichannel data can be parallelized. Option relaxed allows one to finish with a (nearly) unpenalized fit on the L0glm selected features.

The level of regularisation can be manually chosen (if a numeric value is passed to lambda) or preset to approximate an information criterion. Criterion names are case-insensitive. See information_criteria for the penalty formulas, guidance on choosing a criterion and original references.

The L0glm.fit function has various advanced features, including support for both dense and sparse design matrices X, high dimensional data (where p >> n), single or multichannel outcomes y (with either independent or grouped variable selection), support for any of the standard R GLM families & link functions, prior weights (e.g. total nr. of trials for binomial proportions), model offsets, starting coefficient values or dispersion (if prior estimate is available, otherwise optimised), optional nonnegativity constraints, unpenalized terms (e.g. intercept) and shrinkage towards given prior means (via targeted nonzero-centered ridge regression updates).

A choice of several fast C++ coded least squares solvers can be used, which means this function scales well to solve fairly large tall or wide problems with millions of observations or variables. As solvers, we currently implemented the Eigen sparse or dense Cholesky, QR, PCG & LSCG plus BDCSVD and Jacobi SVD solvers, as well as two solvers that natively support nonnegativity constraints: bigGlm from glmnet, and osqp, a fast quadratic programming solver.

Usage

L0glm.fit(
  X,
  y,
  weights = NULL,
  family = gaussian(identity),
  dispersion = NULL,
  start = NULL,
  offset = NULL,
  lambda = 1,
  intercept = (family$family == "binomial"),
  always.include = NULL,
  nonnegative = FALSE,
  normalize = "auto",
  solver = "eigen",
  method = 1L,
  parallelize = TRUE,
  ncores = L0glm_default_ncores(parallelize),
  nthreads = NULL,
  maxit = 25,
  grouped = TRUE,
  group.maxit = 6,
  relaxed = grouped & NCOL(y) > 1,
  prior.mean = NULL,
  thresh = 1e-05,
  tol = 1e-10,
  control.solver = do.call("control.solver.gen", list(solver = solver, method = method)),
  verbose = FALSE,
  openmp = TRUE,
  control.family = list(),
  screen = "auto",
  nsis = NULL,
  screen.maxit = 10L,
  screen.nthreads = NULL,
  validation = NULL,
  precluster = FALSE,
  cluster.refit = TRUE,
  control.precluster = list()
)

Arguments

X

Design matrix with predictor variables of dimension n x p, with n corresponding to the number of observations and p being the number of predictor columns. Use intercept = TRUE to add an intercept, or include a column of ones in X and set intercept = FALSE.

y

Response data with n observations, aligned with the rows of X. For ordinary GLMs, a vector of length n or an n by m matrix of response variables (tasks/channels). Special-family response formats are:

mgcv::nb(): a vector or n by m matrix of nonnegative integer counts, with one column per response channel.
mgcv::multinom(K): a factor/class vector of length n with K+1 categories (one trial per observation), or an n by K+1 nonnegative integer count matrix with one column per category. Row sums give the numbers of trials.
mgcv::cox.ph(): a right-censored survival::Surv(time, status) response, or a numeric two-column matrix with time in the first column and status in the second. Use status=1 for an observed event and status=0 for right censoring. Event indicators belong in y, not in weights.

weights

Optional finite, nonnegative prior/case observation weights. For ordinary GLMs, supply a vector of length n or an n by m matrix; for binomial proportions, weights give trial totals. Special families accept a scalar, a vector of length n, or an n by 1 matrix; NULL uses unit weights.

mgcv::nb(): these formats share weights across response channels; an n by m matrix supplies channel-specific weights.
mgcv::multinom(K): one weight per observation, shared across categories. For a count-matrix response, trial totals come from row sums of y; weights multiply each observation's log-likelihood contribution.
mgcv::cox.ph(): case weights enter both weighted event counts and risk sets. Event indicators belong in y, not in weights.

At least one positive weight is required.

family

Model family and link function. All standard GLM families and their supported link functions are accepted; see family. Special families are mgcv::nb() for negative binomial regression with fixed or estimated shape, mgcv::multinom(K) for multinomial regression with K+1 categories, and mgcv::cox.ph() for right-censored Cox proportional hazards models. See Special families under Details. Fixed-noise scalar models nb1(), Gamma_scale() and Gamma_censored() use the ordinary single- or multichannel response formats with exact likelihood-based IRLS; see signal_families.

dispersion

NULL (default) automatically initializes and estimates dispersion. A positive number fixes it at that value. A negative number initializes it at its absolute value, then estimates it; the supplied magnitude is used directly, without the automatic permissive scaling. For multichannel data, supply one number or a vector with one number per channel; fixed and estimated channels may be mixed. Zero and non-finite values are invalid. Poisson, binomial and negative binomial with known $\theta$ have dispersion fixed at 1; other fixed-scale families respect family$dispersion. See Details for initialization and updating.

start

NULL (default) requests automatic initialization. Positive-mean identity-link models (Poisson, Gamma, inverse Gaussian and negative binomial) with nonnegative designs and a compulsory background or known offset use a positive signal backprojection in original coefficient units; see compulsory_terms. This pilot respects observation weights and individual predictor scales. Other designs use the family-specific ridge pilot. Otherwise supply finite starting coefficients as a scalar, a length-p vector repeated across response columns, or a p by m matrix. Here p includes any added intercept; m is the number of response channels for ordinary GLMs and negative binomial models, K+1 for mgcv::multinom(K), and one for Cox models. Cox models do not add an intercept. Scalars are repeated over coefficients. Numeric scalar start = 1 initializes all coefficients to one in normalized working coordinates, including the intercept. With a free intercept this refers to the centered design. All other scalars and all explicit vectors/matrices use original input coefficient units. To request original-coordinate ones, pass a vector or a p x m matrix; use a matrix to distinguish this from the scalar shorthand in a one-term model. Starting values determine the initial adaptive ridge penalties, subject to the coefficient constraints and the following support restrictions.

Ordinary GLMs and mgcv::nb(): exact zeros exclude individual coefficients separately in each channel, throughout selection, grouped passes and relaxed refits. Rows zero in every channel are screened out.
mgcv::cox.ph(): use a scalar, length-p vector or p by 1 matrix; exact zeros exclude the corresponding predictors.
mgcv::multinom(K): a starting matrix has one column per category, in factor-level or response-matrix column order. Rows are centered across categories to give CLR coefficients. Rows supplied as all zero are excluded before centering and remain zero during selection and relaxed refits. Individual zero entries in an otherwise nonzero row remain eligible. Common nonzero rows, including scalar start = 1, center to zero but remain eligible; a weak-ridge step seeds eligible zero starts.

Intercepts and always.include terms remain eligible even with zero starts. Thus start = 0 fits only mandatory terms, if any, whereas NULL selects automatic initialization. Returned sparse coefficient matrices retain all original predictor rows and response/category columns, including zero rows excluded by the start. For binomial and quasibinomial models with an intercept, automatic slope starts account for its cross-products with the predictors. During selection, additional small coefficients can be excluded according to thresh.

offset

Optional finite, known additive offsets on the linear-predictor (link) scale. For ordinary GLMs, supply a vector of length n or an n by m matrix. For special families, NULL uses zero offsets; a scalar is replicated over observations and response columns, while a length-n vector or n by 1 matrix supplies one offset per observation, shared across response columns.

mgcv::nb(): an n by m matrix supplies channel-specific offsets on the chosen link scale.
mgcv::multinom(K): an n by K+1 matrix supplies category-specific offsets, including every category. Offsets are centered across categories within each row. A scalar or shared observation offset therefore cancels from the multinomial probabilities.
mgcv::cox.ph(): use a scalar, length-n vector or n by 1 matrix of offsets on the log-relative-hazard scale.

lambda

Either a scalar with the $L_0$ penalty for univariate or independent multichannel fits (when grouped=FALSE) or the grouped row-$L_0$/$L_{2,0}$ penalty (when grouped=TRUE), or a case-insensitive criterion name: "aic", "bic", "gic", "ebic", or "mbic". For example, "GIC" and "gic" are equivalent. The preset is a per-coefficient penalty on the $-2\log L$ scale: AIC uses 2, BIC uses $\log n$, and GIC uses $\log p\,\log\log n$. Use preset.lambda to obtain the numeric value directly, for example preset.lambda("gic", p = ncol(X), n = nrow(X)). Here $p$ counts the original selectable predictors, excluding intercepts and always.include terms, before screening by explicit starts. A pool with no selectable predictors uses $p=1$ in logarithms. See information_criteria for all formulas, intended uses and references. Set lambda = 0 to disable penalization. Presets calibrate the adaptive ridge updates to the chosen selection criterion. In L0glm.fit(), an omitted penalty is "gic" with screening and 1 otherwise. L0nnpois() defaults to "gic".

intercept

Logical: Whether to include an intercept term in the model. If `TRUE`, a column of ones will be added to the design matrix `X`, and this term will **not** be penalized. If `FALSE`, no intercept column will be added, and the user must ensure that `X` already contains an intercept if needed. Defaults to `TRUE` for binomial models and `FALSE` otherwise in `L0glm.fit()`, and to `FALSE` in `L0nnpois()`. With `nonnegative = TRUE`, the intercept is also constrained to be nonnegative, except for binomial and quasibinomial models, where it may have either sign.

always.include

A vector of coefficient indices that should always be estimated & included in the model. These terms are unpenalized and remain eligible when their current estimate or explicit starting value is zero. Coefficients with a name that includes "Intercept" in their name are always included. Use NULL or 0 to request no additional mandatory terms; recognized intercepts remain mandatory. If `intercept = TRUE`, the first column (intercept) will always be estimated and included in the model. Indices refer to columns of the supplied `X`. Nonnegative compulsory terms are refined in a small constrained block; see [compulsory_terms] for smooth-background modelling.

nonnegative

A single logical value. If TRUE, all coefficients, including any intercept and always.include terms, must be nonnegative. For binomial and quasibinomial models only, intercept columns (columns of ones, whether supplied in X or added by intercept = TRUE) remain unrestricted. If FALSE, all coefficients may have either sign. Defaults to FALSE in L0glm.fit(). The specialized L0nnpois() function always enforces nonnegativity and has no nonnegative argument. Unpenalized terms use a constrained conditional solve, while selectable coefficients use the chosen ridge solver. See [compulsory_terms].

normalize

Preprocessing convention: "auto" (default) uses "rms" for signed fits or a free intercept, and "none" for fully nonnegative fits. Other choices are FALSE or "none", TRUE or "l2" for unit L2 norm, "rms" for unit root-mean-square ($\sum_i X_{ij}^2/n=1$), and "sd" for unit sample standard deviation. With a free intercept, column means are accounted for and the norms use centered columns. The original X and y are unchanged. Dense working columns are centered once; sparse working columns are shifted only when this adds no stored entries, with the other means handled by the intercept. A constrained intercept prevents this centering, including in L0nnpois(). Mandatory terms and zero-scale columns use scale one. Penalties, starts, ridge centers and bounds are transformed consistently. Returned coefficients and explicit vector/matrix starts retain their input units; numeric scalar start=1 is the normalized-coordinate shorthand described above. Non-"none" modes also use these scales for automatic initialization, thresholding and convergence. Automatic signed fits are therefore invariant to changes of predictor units, up to numerical precision. Sparse designs remain sparse. The training centers and scales are returned in preprocessing; predictions use the original columns directly.

solver

Solver used for the ridge regression step in the IRLS algorithm (see details). It should be one of "osqp" (osqp quadratic programming solver), "bigGlm" (using bigGlm coordinate descent; alias "glmnet") or "eigen" (one of the RcppEigen Eigen solvers). "osqp" and "bigGlm" support nonnegativity constraints and scale well to larger problems, with especially "bigGlm" being suitable for very large problems. The "eigen" solvers apply nonnegativity constraints by clamping negative proposals to zero. This works well in practice when the true active effects are positive and regularization is well chosen, for example lambda="gic": ridge shrinkage stabilizes the updates and often keeps those effects positive. In the package's deconvolution solver benchmarks, Eigen Cholesky LLT generally gives better support recovery and coefficient accuracy than the formally constrained alternatives, together with much faster fits. See the measured solver comparisons. Use method=1L for Eigen Cholesky LLT.

method

When solver="eigen", method used for the eigen solver (see details). It should be an integer from 1 to 6 :

1L = Cholesky LLT solver, with sparse AMD ordering by default

2L = Cholesky LDLT solver, with sparse AMD ordering by default

3L = Preconditioned Conjugate Gradient (PCG) with incomplete Cholesky preconditioner

4L = Least Squares Conjugate Gradient (LSCG) with least squares diagonal preconditioner

5L = Column-pivoted QR solver. Sparse QR uses numerical column pivoting and COLAMD fill-reducing ordering by default.

6L = BDCSVD solver for dense systems, or Jacobi SVD for sparse systems.

Defaults to 1L (Cholesky LLT), which is often the fastest.

parallelize

A logical value indicating whether in the multivariate case, the columns of the response matrix y should be fitted in parallel. If TRUE, each column of y is fitted concurrently, potentially speeding up computation. If FALSE, the columns are fitted sequentially.

ncores

Optional integer specifying the number of cores to use when parallelize=TRUE. If parallelize=TRUE, defaults to half the detected number of physical cores when that can be determined, falling back to half the number of logical cores. The requested value is capped at the detected number of physical cores when that can be determined. With the L0nnpois backend, ncores instead controls OpenMP threads sharing the design matrix and is capped at the logical core count.

nthreads

Optional integer specifying the number of BLAS threads to use per fit process. The default NULL leaves the main R session's thread settings unchanged and uses one BLAS thread per parallel worker. Setting nthreads can be useful to benchmark serial fits or tune per-worker threading on systems with a multithreaded BLAS. The L0nnpois backend uses ncores; it ignores nthreads.

maxit

Integer specifying the maximum number of IRLS iterations to perform per channel and fitting pass (default 25). Fits stop when the convergence criterion is met; inspect convergence to distinguish this from reaching the iteration cap.

grouped

Logical (default TRUE) choosing the unit of sparse selection. For multichannel ordinary GLMs and negative binomial models, TRUE uses a grouped row-$L_0$/$L_{2,0}$ penalty: a predictor is counted once when any response-channel coefficient in its row is nonzero. FALSE uses elementwise $L_0$ selection separately in each channel. For special families:

mgcv::nb(): grouped selection shares adaptive ridge weights across each predictor row and thresholds its row norm; ungrouped selection uses channel-specific coefficients and weights. Each channel retains its own shape parameter, whether grouping is enabled or not.
mgcv::multinom(K): TRUE selects predictor rows using the norm of their K+1 CLR category coefficients and counts each nonzero row once. FALSE selects and counts individual nonzero CLR entries. Both fit the same joint multinomial likelihood with zero-sum coefficient rows; ungrouped fitting is not independent binomial regression and does not select entries of an arbitrary ILR contrast basis. A nonzero CLR row needs at least two nonzero category entries.
mgcv::cox.ph(): there is one coefficient per predictor, so grouped and ungrouped selection coincide.

Mandatory terms are unpenalized. Starting-value exclusions remain in force; the returned sparse coefficient matrix keeps all original dimensions.

group.maxit

Integer specifying the maximum number of iterations to perform for the group penalized fitting procedure when grouped = TRUE. In each iteration, the columns of the response matrix y will be refitted on the active row support using shared row-level starting values derived from the current coefficients, which then imposes a group penalty. This argument applies to ordinary GLMs, not the special-family loops; those update grouping within IRLS and use maxit per pass.

relaxed

Logical controlling a final refit to reduce shrinkage bias. FALSE returns the selection fit. For ordinary GLMs, TRUE uses a nearly unpenalized refit with lambda = 1e-30 and unsmoothed adaptive ridge weights. Special families instead use an unpenalized IRLS pass with lambda = 0, on the following support:

mgcv::nb(): grouped refits admit all eligible channel coefficients within selected predictor rows; ungrouped refits admit only entries nonzero at the end of selection. Estimated shape parameters continue to be updated during refitting; shapes fixed in the constructor remain fixed.
mgcv::multinom(K): grouped refits jointly estimate all category coefficients within retained predictor rows, allowing category entries that were zero during selection to become nonzero. Ungrouped refits jointly estimate only the selected nonzero CLR entries and hold all other category entries at zero. Both retain the joint multinomial likelihood and zero-sum row constraints; a row with $s$ admitted category entries has $s-1$ free directions. No new predictor rows are introduced.
mgcv::cox.ph(): refit the selected log-relative-hazard coefficients using the Cox partial likelihood, without an intercept.

Hard thresholding is disabled during refitting. Relaxed nonnegative fits finish with a conditional boundary pass: a positive coefficient may be set to exactly zero when the endpoint derivative supports zero and the actual deviance or negative log likelihood does not increase. This uses no coefficient-size cutoff. It supports ordinary GLMs, negative binomial at its fitted shape and Cox partial likelihood; multinomial CLR coefficients remain unrestricted. A validation checkpoint is retained if polishing would increase its held-out deviance. Explicit starting-value exclusions are respected; intercepts and always.include terms remain eligible and unpenalized. Sparse output retains the original dimensions. The default is TRUE for grouped multichannel ordinary GLMs and NB, and for grouped multinomial models with either factor or count-matrix responses. It is FALSE for Cox, ungrouped fits and other single-channel fits. Set relaxed = TRUE explicitly to refit an ungrouped model. Screened L0glm.fit() models instead default to TRUE, including single-channel models, unless explicitly overridden. Each special-family pass has its own maxit limit; inspect convergence to assess whether the refit converged.

prior.mean

Optional finite centers for the ridge penalty in each IRLS update. NULL (default) uses zero centers. Supply values in original input coefficient units; normalization and any eligible predictor centering are handled internally. Unlike start = 1, prior.mean = 1 has no special normalized-coordinate meaning. For ordinary GLMs, use a length-p vector shared across channels or a p by m matrix, including any added intercept row. Special families also accept a scalar repeated over predictors and response/category columns:

mgcv::nb(): use a length-p vector or p by m matrix for regression-coefficient centers. These do not specify a center for the shape parameter $\theta$ or dispersion.
mgcv::cox.ph(): use a length-p vector or p by 1 matrix for log-relative-hazard coefficients, without an intercept. These are not centers for the baseline hazard.
mgcv::multinom(K): use a p by K+1 matrix in response-category order. Targets are centered across categories within each predictor row to give CLR targets, $a_{jc} \leftarrow a_{jc}-(K+1)^{-1}\sum_h a_{jh}$. Adding a common constant to all category targets in a row has no effect. A scalar or length-p vector is shared across categories and therefore gives zero CLR targets. Supply category contrasts for nonzero CLR targets.

The ridge penalty is centered on these targets, but adaptive weights and support counts are still based on coefficients relative to zero, not on deviations from prior.mean. These targets are not starting coefficients and do not change the eligibility rules of start: a nonzero target cannot reactivate an excluded starting coefficient or an all-zero multinomial starting row. Such rows remain zero in the full-size sparse output. Unpenalized intercepts and always.include terms are not shrunk toward their targets. For special families, the final relaxed refit removes the ridge penalty on retained support, so these centers no longer exert shrinkage during that refit.

thresh

Finite nonnegative scalar, default 1e-5, controlling hard thresholding during selection. Cutoffs use normalized working coefficients $b$; for selectable slopes, abs(b) = abs(beta)*scale, where beta is in original input units and scale = 1 with normalize = "none". Ordinary GLM updates use individual coefficient magnitudes, comparing each with thresh times the largest selectable magnitude in its response channel, or the absolute floor 1e-20. For special families:

Grouped multichannel NB and grouped multinomial models threshold $r_j=\|b_j\|_2$. An entire predictor row is removed when $r_j < \texttt{thresh}\max_{k\in S} r_k$ or $r_j < 10^{-20}$, where $S$ contains selectable, nonmandatory rows. Individual entries within retained rows are not hard thresholded, even if very small.
Ungrouped NB and multinomial models threshold individual entries: $|b_{jc}| < \texttt{thresh}\max_{k\in S}|b_{kc}|$ or $|b_{jc}| < 10^{-20}$. Each response/category column has its own relative cutoff. Cox and single-channel NB use this scalar rule.
After multinomial thresholding, surviving entries are recentered within their allowed category set to maintain zero-sum CLR rows, leaving excluded entries at zero. Rows with fewer than two allowed entries are zero. Row norms, not arithmetic row means, are used for grouped thresholding.

Mandatory terms are excluded from cutoff maxima and never hard thresholded. Removed rows or entries stay excluded during the selection pass; grouped relaxed refits may recover eligible entries within retained rows as described in relaxed. Starting-value exclusions remain binding, and sparse output keeps its original dimensions. Setting thresh = 0 disables only the relative cutoff; the absolute floor remains. Final relaxed refits apply neither cutoff. For special families, cutoffs are applied only in selection passes with lambda > 0. Coefficient constraints still apply. thresh does not smooth or otherwise alter the adaptive ridge weights.

tol

Convergence tolerance (default 1e-10). Iterations stop when the relative change in the $L_2$ norm of the coefficients drops below this threshold. With preprocessing, the test uses scaled selectable slopes and, when mandatory terms exist, changes in predictions, so that the predictor or response origin cannot hide slope changes. When dispersion is estimated, selection also waits until its update has stabilized as described in Details, subject to maxit.

control.solver

List of parameters controlling the optimization solvers. Auto-initialised with control.solver.gen when NULL. See ?control.solver.gen for more details on the available parameters for the different solvers. For sparse Eigen Cholesky LLT/LDLT, control.solver$NaturalOrdering selects Eigen natural ordering (TRUE) or AMD ordering (FALSE, default). Sparse QR defaults to Eigen COLAMDOrdering (ordering=1L); ordering=0L uses AMD and ordering=2L uses NaturalOrdering.

verbose

Logical indicating if the algorithm progression should be printed to the console.

openmp

Logical; defaults to TRUE. For nonnegative identity-link Poisson fits with solver="eigen" and method=1L, use the compiled Rcpp OpenMP L0nnpois backend. Validation-based stopping is supported directly by this backend. Other families, links and solvers use the general R fitting path. Set FALSE to use that path explicitly. Solver controls requesting another factorization or rowaug=FALSE also retain the R path. This option selects the backend; parallelize controls whether its outcome channels run concurrently.

control.family

Family-specific IRLS controls. A named list with multinomial="irls" (joint block IRLS) or "bound" (upper curvature bound), theta.min=1e-6, theta.max=1e7, theta.tol=1e-6, max.working=2000L, cox="auto", cox.llt.max=c(dense=1024L,sparse=0L), cox.cg.tol=1e-10 and cox.cg.maxit=2000L. Cox also accepts "direct" or "pcg"; a scalar LLT threshold applies to both storage types. The working-size guard applies to dense direct Cox and multinomial systems, not matrix-free Cox products. See the dedicated family help pages.

screen

Screening strategy in L0glm.fit(): "auto" (default), "none", "SIS" or "ISIS", ignoring case. Auto enables ISIS exactly when $p_{\mathrm{eligible}}>10n$, where $n$ is the number of rows of X. The eligible count excludes mandatory terms and predictors excluded by explicit zero starts in every response channel. No separate minimum sample size or predictor count is imposed. Explicit SIS/ISIS works at any size; "none" disables it. SIS screens once; ISIS repeatedly screens using the fitted predictor as an offset and refits the original response on the expanded candidate set. See Screening under Details and sis_glm.

nsis

Positive integer screening batch size, capped at ncol(X). By default SIS uses $2\lfloor n/\log n\rfloor$; ISIS uses $\operatorname{round}\{n/(\log p\,\log\log n)\}$, at least one. Here $p$ is the original selectable predictor count. ISIS also retains currently selected and mandatory rows. Ignored when screening is off. With validation, a vector compares several batch sizes using validation error; NULL retains the GIC-based ISIS default.

screen.maxit

Maximum ISIS rounds, default 10. A repeated candidate set stops the search early; SIS always uses one round. This is separate from maxit, the inner IRLS limit. Ignored when screening is off.

screen.nthreads

Positive integer OpenMP thread count for screening. NULL (default) uses the number of physical cores, falling back to one if detection is unavailable. Independent of nthreads (BLAS) and ncores (outcome fitting). See sis_glm for caps and serial fallback. Ignored when screening is off.

validation

Optional list for validation-based IRLS stopping and batch/start selection, supported for all families and links. Supply X and y for held-out observations, using the same predictor columns and response format as the training data. Named predictors are matched automatically. Optional weights and offset accept scalars, observation vectors or matrices following the training formats. Controls are patience = 3L, minit = 3L, tol = 1e-8, and a scoring-only mu_floor = 1e-8 for positive-mean families. Gaussian identity models compare starts = c("default", "marginal"), with conditional marginal slopes as a second initializer; other models use starts = "default". An explicit start restricts the comparison to that start and preserves its exclusions. rule = "min" chooses minimum validation error; rule = "1se" chooses the sparsest candidate within one empirical standard error of the minimum. The default NULL fits to the usual numerical stopping criterion.

precluster

Enable response-channel preclustering, default FALSE. Currently supports nonnegative identity-link Poisson models. See preclustering for the fitting sequence and mathematical model.

cluster.refit

With preclustering, refine selection on every original channel (default TRUE). FALSE returns a model of summed cluster outcomes at the requested lambda. This is separate from relaxed.

control.precluster

Named advanced settings for clustering and screening. Defaults include at most 32 groups, Ward Poisson clustering, GIC cut selection, unit row bins and a screening factor of 0.1. See preclustering for all entries.

Value

An object of class L0glm_fit, with list components below. Special-family models additionally inherit from L0glm_extended_fit. Use coef(), logLik(), AIC(), BIC(), GIC(), EBIC() and MBIC() to extract coefficients and scores; see model_criteria and L0glm.metrics. print() shows a compact overview; summary() adds fit metrics and a bounded coefficient display; see model_summary. predict() and fitted() provide original-scale predictions; see model_predictions for new data, offsets and special-family types.

coefficients: a sparse dgCmatrix with the fitted coefficients (single column for the univariate case). If `intercept = TRUE`, the first row is for the intercept.
weights: final IRLS working weights for ordinary GLMs and native Poisson; prior observation weights for special-family models.
prior.weights, offset: original prior observation weights and offsets retained for ordinary-model scoring. Special-family models retain prior weights in weights.
dispersion, dispersion.estimated: fitted dispersions and flags distinguishing estimated from fixed dispersion for ordinary GLMs. Special families retain their nuisance-parameter information separately.
init.lambda: initial lambda.
lambdas: a vector or matrix with the adaptive penalties lambdas.
family: The error structure used to fit the data.
z: a vector or matrix containing the working response in the final iteration of L0glm.fit
X: design matrix X, including the intercept column if added
y: response vector or matrix y.
p: original number of selectable predictors.
screening: screening method, batch size, candidate/retained rows, outer iteration history, best round and stopping information when screening is used; otherwise absent (NULL).
validation: when requested, validation deviance and separate selection and relaxed diagnostics, including the best iteration, number of attempted iterations, stopping flag and per-iteration history. Multichannel ordinary models also retain pass/channel diagnostics. The starts table records each start/batch candidate's validation deviance, fitting time and nonzero count; selected_start and selected_nsis identify the returned model. minimum, se, rule and scoring describe the selection criterion.
always.include: indices of unpenalized mandatory coefficients, including recognized intercepts, in the returned coefficient matrix.
preprocessing: training column centers and scales, resolved normalization mode, centering eligibility, and original selectable count. Coefficients are already on the original input scale.
convergence: one row per channel and selection or relaxed pass, containing its stage, pass, channel, iteration count, convergence flag, number of dispersion updates and number of reactivations after a freeze. The flag records whether the iterative stopping criterion was met.

Details

L0glm approximates the $L_0$, independent multichannel elementwise $L_0$, or grouped row-$L_0$/$L_{2,0}$ penalized GLM objectives via an iterative adaptive ridge regression procedure.

Algorithm for univariate, single-channel data

For univariate data, the approach builds on L0 penalized regression by de Rooi & Eilers (2011), Frommlet & Nuel (2016), Liu et al. (2017), Dai et al. (2018 & 2020) & Li et al. (2021) and is based in a simple modification of the IRLS (Iteratively Reweighted Least Squares) algorithm used to fit GLMs (cf. McCullagh & Nelder 1989, Hastie & Pregibon 1992, Hardin & Hilbe 2007). In this algorithm, each iteration of the IRLS algorithm updates the coefficient estimates using a quadratic approximation of the log likelihood of the GLM objective. This is done using a weighted least squares step. In our adapted algorithm, this step is replaced by a (zero or nonzero centered) adaptive ridge regularised regression with optional nonnegativity constraints on the coefficients. In each iteration, this uses the following ridge penalties for the penalized terms (excluding the intercept, if included): $$\Lambda = (1/4)\lambda\sigma^2 w$$ with adaptive penalty weights $$w = 1 / |\beta|^2$$ where $\beta$ are the coefficient estimates in each iteration. During selection, the relative cutoff and absolute numerical floor specified by thresh remove coefficients from later solves. The final relaxed refit does not apply this hard cutoff. The inverse-square weights follow Frommlet & Nuel (2016). Nonnegativity constraints are natively supported if the bigGlm or the osqp solver is used, or are obtained via clipping if the eigen solvers are used. The factor of 1/4 was shown to be required to approximate the L0 penalty objective and has been proven for the Gaussian case and under orthogonality (Frommlet & Nuel 2016), but in practice also works well for non-orthogonal designs and other noise families. In an IRLS step, fixed working weights and dispersion give a weighted least-squares quadratic with curvature $h_j=x_j^T W x_j/\phi$. For weighted-orthogonal columns its L0 selection boundary is $h_j\widehat b_j^2=\lambda$; the adaptive-ridge fixed-point equation $b_j^2-\widehat b_j b_j+\lambda/(4h_j)=0$ has the same boundary. Thus the quarter factor carries over directly to each working quadratic. Standard likelihood regularity conditions motivate this approximation increasingly well as sample size grows on an identifiable active model. Correlated predictors couple the joint support decisions. If desired, L0glm.tune() searches around the preset using the information criterion or cross-validated deviance. The Algorithm and theory vignette details the quadratic and large-sample arguments. Dispersion (the variance for Gaussian models) can be fixed or estimated during IRLS using dispersion = NULL or a negative initial value.

Adaptive ridge weights

The exponent is fixed at two. Weights use normalized working coefficients $b$ relative to zero. Ordinary single-channel ridge diagonals are $D_j=(\lambda/4)\phi/|b_j|^2$, with scale dispersion $\phi$. Special families have scale dispersion one. Grouped multichannel NB and multinomial fits use $D_{jc}=\lambda/(4\|b_j\|_2^2)$, shared across each row. For multinomial models the norm includes all K+1 CLR category coefficients. Ungrouped NB and multinomial fits use $D_{jc}=\lambda/(4|b_{jc}|^2)$; Cox and single-channel NB use the same scalar rule. Mandatory terms have zero ridge penalty and excluded coefficients remain excluded. Hard thresholding is controlled separately by thresh. Special-family relaxed refits use lambda = 0.

Optimization interpretation

Adaptive ridge uses successive quadratic penalties to promote sparsity. Abergel, Bouaziz & Nuel (2024) derive a variational representation for $L_q$ penalties and explain their majorization–minimization interpretation. With zero ridge centers and fixed dispersion, inverse-square weights give tangent quadratic penalties for a log penalty on the current nonzero coefficients. The inequality $$\log s \leq \log s_t + (s-s_t)/s_t,\qquad s,s_t>0,$$ with $s=b_j^2$, yields a quadratic coefficient proportional to $1/(b_j^{(t)})^2$. For grouped norms the same argument uses $s=\|b_j\|_2^2$. These inverse-square weights give the broken adaptive ridge rule. IRLS approximates the likelihood, and hard thresholding reduces subsequent systems to the active support. Information criteria evaluate the resulting models using their likelihood and support size. Exact surrogate minimization with a fixed likelihood and penalty gives the usual MM descent argument; the package also accommodates dispersion updates, constraints and relaxed refits as described below.

Hierarchical Bayesian and ARD interpretation

A ridge update is a conditional maximum a posteriori (MAP) estimate under Gaussian coefficient priors. With the default zero ridge center, the inverse-square update corresponds to prior precision $\alpha_j=\lambda/(4\beta_j^2)$ in the solver's coefficient coordinates; its ridge coefficient is $\phi\alpha_j$ for GLM dispersion $\phi$. The prior remains zero-centered and skeptical: small fitted effects receive concentrated priors and strong effects receive wider priors. The next IRLS step updates the coefficient mode using a local quadratic likelihood. This gives a plug-in empirical-Bayes-style connection to automatic relevance determination (ARD), with shared precisions providing a group-ARD view. Evidence-maximizing ARD also accounts for posterior uncertainty; the adaptive ridge rule here updates scales from coefficient point estimates. The ridge center stays fixed; prior.mean specifies its value. See Tipping (2001) and the Algorithm and theory vignette for this connection.

Automatic initialization uses one weak ridge fit to the initial working response, with columns scaled to unit weighted norm and ridge constant 1e-6. Four reproducible random-sign probes share this factorization to estimate the residual noise degrees of freedom, $\operatorname{tr}\{(I-H)^T(I-H)\}$, where $H$ is the pilot hat matrix. The caller's random-number state is preserved. The pilot residual scale is divided by max(1, preset.lambda("gic", p, n)) for a permissive start independent of the fitted lambda; its support is never imposed on the main fit. A supplied negative value bypasses this pilot. With preprocessing enabled, mandatory terms are profiled out and column information is conditional on them, preserving changes of predictor origin when an unrestricted intercept is present. The normalized pilot also checks the effective number of residual directions, $\nu=\operatorname{tr}(R^2)^2/\operatorname{tr}(R^4)$, with $R=I-H$. If $\nu<\min(10,(N-q)/5)$, where $q$ is the rank of the mandatory terms, or the estimated variance is below 1e-6 times the nuisance-adjusted working-response variance, one retry uses pilot ridge 0.1. This stabilizes initialization near interpolation. The probes reuse the same factorization, and ordinary pilots pass a trace bound without additional solves.

Starting at iteration 3, when the number $k$ of active coefficients is below $0.7N$, dispersion is updated toward $$\widetilde\phi = \frac{\sum_{i:w_i>0}w_i(y_i-\mu_i)^2/V(\mu_i)}{N-k}.$$ Here $N$ counts positive-weight observations, $k$ includes the intercept, and $V$ is family$variance, evaluated at the newly fitted mean. Each update is bounded to $[\phi/2,2\phi]$. After three consecutive updates with unchanged support and a target within 1% of the current dispersion, it is held fixed until the support changes again. A later support loss reactivates these bounded updates in the same IRLS loop. Each subsequent grouped selection pass may update it again on the new support. A relaxed refit holds dispersion fixed during its coefficient solves and reports the final Pearson residual scale for estimated channels; user-supplied positive dispersions remain fixed throughout.

These inexpensive Pearson moment updates use the active coefficient count to approximate residual degrees of freedom. In nearly interpolating models, little residual information is available to estimate noise; if the active model remains too large for an update, the starting scale is retained. The pilot normally requires one factorization (two with the rare retry); subsequent scale updates require only $O(N)$ arithmetic, without nested fitting or further factorizations. Supplying a dispersion value avoids the pilot factorization, which can be appreciable for very large dense designs even when the main solver is iterative.

Shared numerical safeguards

Ordinary GLMs start their working model with family$initialize. Each ridge solve proposes a coefficient update, followed by a shared backtracking check. Dispersion and adaptive ridge penalties $d_j$ are held fixed during this check. For ordinary GLMs its target is $$Q(\beta)=D(y,\mu(\beta))+ \sum_{j:d_j<\infty}d_j(\beta_j-a_j)^2,$$ where $a_j$ is prior.mean. Special families use their joint log likelihood in the corresponding fixed-penalty target. Smooth updates on an unchanged support must decrease this target within numerical tolerance. If necessary, the coefficient step is halved until its means are valid and its objective is finite and acceptable. Predictor differences are reused, without another ridge solve or factorization. A feasible backtracking anchor leaves the requested first-step penalty weights and coefficient eligibility intact.

Hard thresholding follows the smooth update as a discrete selection step. Bound clipping can also remove coefficients; these support changes require finite valid predictions but can increase the smooth ridge target. This preserves the adaptive-ridge selection path while safeguarding continuous refits. A threshold deletion is retained only when it leaves valid means. Starting-zero exclusions and coefficient bounds remain in force; relaxed refits apply no hard threshold. Failed backtracking retains the previous finite coefficients with a warning and converged=FALSE.

Distribution-specific kernels provide valid parameter domains, working weights and stable deviance evaluations. For example, the log-link inverse Gaussian uses $D=\sum_i w_i(y_i/\mu_i-1)^2/y_i$ and $W_i=w_i/\mu_i$ to avoid intermediate overflow. The step-acceptance rule is shared across families and the compiled Poisson implementation. See Marschner (2011) for deviance backtracking in ordinary GLMs.

Fisher scoring uses the expected Hessian in each quadratic approximation. If option relaxed is set to TRUE the algorithm will refit the L0glm selected variables using a nearly unpenalized model (with lambda=1E-30) without hard thresholding. This can reduce bias and increase sensitivity. If an appropriate level of regularisation lambda is chosen, the L0 norm penalized solution approximates Best Subset selection. To achieve variable selection consistency, the "bic" lambda preset is recommended for low dimensional settings, whilst "gic", "ebic" or "mbic" may work better for high dimensional settings, subject to the conditions discussed in information_criteria. Setting it to "aic" is recommended if predictive performance is of key interest, as the AIC criterion asymptotically approximates the leave-one-out cross-validation error.

Algorithm for multivariate/multichannel data

For the multivariate ("multitask learning") case there are two options. First, if grouped=FALSE, an independent L0glm model is fit on each of the responses, which corresponds to an elementwise $L_0$ penalty on the coefficient matrix. Second, if grouped=TRUE, grouped variable selection will be performed via a row-$L_0$ penalty, conventionally denoted $L_{2,0}$ in the multitask/MMV literature. This is achieved by taking the independent elementwise $L_0$ solution, and then, for all variables that have nonzero coefficients in at least one of the channels, repeating the $L_0$ penalized fits but using a shared row-level start derived from the coefficients across all channels (response variables) so as to achieve a group penalty and perform grouped feature selection. Active rows are detected by a row-wise nonzero count, equivalently by any positive within-row norm. With signed coefficients, the shared start uses mean absolute magnitude with the sign of the row mean. If that mean is exactly zero, the sign of the first largest-magnitude entry is used. Thus opposite signs across channels cannot reduce the starting magnitude by cancellation. Nonnegative fits retain arithmetic averaging. This is then further iterated until convergence. Finally, as for single channel data, if relaxed=TRUE, the algorithm can finish with a final step whereby the selected variables are refit using a nearly unpenalized model (with lambda=1E-10). Compared to other group penalized regression methods, such as the group LASSO or multitask lasso, which use an $L_{2,1}$ penalty, this approach produces a sparser solution that better approximates the Best Subset of Groups solution if lambda is adequately chosen.

Screening

Screening reuses sis_glm for all ordinary and special families. With screening enabled, omitted lambda, relaxed and parallelize use "gic", TRUE and FALSE, respectively. Explicit arguments are preserved. Other defaults are unchanged. This reproduces the small-batch ISIS workflow used in the ultra-wide Gaussian benchmark with a single L0glm.fit(X, y) call.

Each candidate model starts afresh. The original offsets and prior observation weights are used for refitting; the fitted linear predictor is added to the offset only for the next screen. The best GIC round is returned, or the last round for quasi families without an available GIC. Penalty presets and GIC always use the original selectable predictor count, before screening. Explicit starting-zero exclusions remain in force, including per-channel exclusions, while mandatory terms stay eligible. Coefficients and predictions retain the original dimensions and units; screened-out coefficients are zero. Sparse designs stay sparse and dense designs stay dense.

screening records the resolved method, automatic-dispatch flag, batch size, candidates, retained indices, outer history and chosen round. It also records the number of OpenMP screening threads actually used. Outer stopping describes the screening search; inner convergence describes each fitted model. Use screen="none" to fit the full design directly.

Special families

The constructors mgcv::nb(), mgcv::multinom(K) and mgcv::cox.ph() select family-specific IRLS implemented using this package's ridge kernels. NB estimates one shape per channel unless fixed in the constructor. Multinomial fits return all K+1 category coefficients in CLR coordinates, with a default free intercept and invariant grouped row penalties. Cox fits use right-censored responses, Breslow ties, no intercept, cached risk sets and automatic LLT/matrix-free PCG routing. See L0glm_negative_binomial, L0glm_multinomial and L0glm_cox for the mathematics, examples and controls. These family loops are serial: parallelize, ncores, openmp and group.maxit do not govern them. Ordinary GLM and native Poisson paths use the parallelize, ncores and openmp controls.

Validation-based stopping

validation=list(X=V,y=v) tracks held-out deviance using a common stopping rule. Ordinary GLMs use family$dev.resids, averaged per weighted response; Gaussian identity therefore uses mean squared error. Multichannel models aggregate response-channel deviances. Multinomial models use joint deviance per weighted trial, and Cox models use partial deviance per weighted event with risk sets formed from the validation observations. For negative binomial models with estimated shape, the saturated reference is fixed at the Poisson saturated likelihood, so changing training theta cannot shift the reference between iterations. This gives a predictive log-likelihood score; a fixed theta uses ordinary NB deviance. After at least minit iterations, patience consecutive iterations without an improvement exceeding tol * max(1, best_deviance) stop the pass. The best checkpoint is returned, including its dispersion and penalties. A training-only nuisance fit supplies the initial checkpoint. Relaxed refitting starts from the selected checkpoint and can retain it if refitting increases validation error. SIS/ISIS score their candidate rounds on the same validation set. A vector nsis, with SIS/ISIS enabled, compares its batch sizes at the requested lambda. Each candidate starts afresh, and the nuisance-only model is included in the comparison. Automatic starts compare the default initializer and, for Gaussian identity models, uncontracted conditional marginal slopes. Starts use training data only. An explicit start restricts the comparison to that start. For rule="1se", the standard error uses per-observation squared validation losses from the minimum-error candidate; Cox uses event contributions. Among eligible models, fewer nonzero coefficients take priority, followed by a smaller batch size. With nsis=NULL, each screened candidate uses the usual default batch. Coefficient updates use training data only. Validation selects the stopping point, so use separate test observations to assess prediction performance.

References

M.E. Tipping, Sparse Bayesian Learning and the Relevance Vector Machine. Journal of Machine Learning Research 1, 211-244 (2001). https://www.jmlr.org/papers/v1/tipping01a.html.
J. de Rooi & P. Eilers, Deconvolution of pulse trains with the L0 penalty. Analytica chimica acta 705, 218-226 (2011).
F. Frommlet & G. Nuel, An adaptive ridge procedure for L0 regularization. PloS one 11, e0148620 (2016).
I. Marschner, glm2: Fitting Generalized Linear Models with Convergence Problems. The R Journal (2011). doi:10.32614/RJ-2011-012 .
Z. Liu, F. Sun & D.P. McGovern, Sparse generalized linear model with L0 approximation for feature selection and prediction with big omics data. BioData mining 10, 1-12 (2017).
L. Dai, K. Chen, Z. Sun, Z. Liu & G. Li, Broken adaptive ridge regression and its asymptotic properties. Journal of multivariate analysis 168, 334-351 (2018).
L. Dai, K. Chen & G. Li, The broken adaptive ridge procedure and its applications. Statistica Sinica 30, 1069-1094 (2020).
N. Li, X. Peng, E. Kawaguchi, M.A. Suchard & G. Li, A scalable surrogate L0 sparse regression method for generalized linear models with applications to large scale data. Journal of Statistical Planning and Inference 213, 262-281 (2021).
R. Abergel, O. Bouaziz & G. Nuel, A review on the Adaptive-Ridge Algorithm with several extensions. Statistics and Computing 34, 140 (2024). doi:10.1007/s11222-024-10440-6 .
W. Guo, S. Ma & Z. Liu, l0ara: Sparse Generalized Linear Model with L0 Approximation for Feature Selection. R package version 0.1.6 (2022). Available at https://CRAN.R-project.org/package=l0ara
J.W. Hardin & J.M. Hilbe, Generalized linear models and extensions. (Stata press, 2007).
T.J. Hastie & D. Pregibon, Generalized linear models. Chapter 6 of Statistical Models in S eds J. M. Chambers and T. J. Hastie, Wadsworth & Brooks/Cole (1992).
P. McCullagh & J.A. Nelder, Generalized Linear Models. London: Chapman and Hall (1989).

Author

Tom Wenseleers, Laura Lanzas & Christophe Vanderaa

Tom Wenseleers tom.wenseleers@kuleuven.be \ Laboratory of Socioecology and Social Evolution, KU Leuven \ https://bio.kuleuven.be/eeb/tw/research \ ORCID: 0000-0002-1434-861X

Examples

# Single-channel Poisson deconvolution: shifted copies of a known PSF.
single <- simul(method=2, n=500, p=500, s=.1, seed=2)
fit_single <- L0glm.fit(single$X, single$y, family=poisson(identity),
  lambda="gic", nonnegative=TRUE, relaxed=TRUE, parallelize=FALSE)
L0glm.metrics(fit_single, beta_true=single$beta_true)[
  c("gic", "TP", "FP", "FN", "TN", "F1_score", "MCC")]

old_par <- par(no.readonly=TRUE)
par(mfrow=c(1, 1), mgp=c(2, .7, 0))
# The default PSF has peak one. Scale both spike trains by the same PSF peak
# to display them in signal units; negative values indicate mirroring only.
spike_scale <- max(single$X)
beta_single <- as.numeric(coef(fit_single))
mu_single <- as.numeric(fitted(fit_single))
height <- max(as.numeric(single$y), as.numeric(single$y_true), mu_single,
  as.numeric(single$beta_true)*spike_scale, beta_single*spike_scale)
plot(single$x, single$y, type="l", col="grey65", ylim=c(-1, 1)*height*1.15,
  xlab="Time", ylab="Signal / spike height (fit mirrored)",
  main="Single-channel Poisson deconvolution")
abline(h=0, col="grey75", lty=3)
lines(single$x, single$y_true, col="grey30", lty=2)
segments(single$x_beta, 0, single$x_beta,
  as.numeric(single$beta_true)*spike_scale, col="navy")
lines(single$x, -mu_single, col="orange", lwd=1.5)
segments(single$x_beta, 0, single$x_beta, -beta_single*spike_scale, col="red")
legend("topright", c("Observed signal", "True signal", "True spikes",
  "Fitted signal (mirrored)", "Fitted spikes (mirrored)"),
  col=c("grey65", "grey30", "navy", "orange", "red"),
  lty=c(1, 2, 1, 1, 1), bty="n", cex=.75)

# Multichannel Poisson deconvolution: select rows shared across channels.
# Retain every mass channel generated by the simulation.
multi <- simul(method=4, n=500, p=500, s=.05, seed=2)
fit_multi <- L0glm.fit(multi$X, multi$y, family=poisson(identity),
  lambda="gic", nonnegative=TRUE, grouped=TRUE, relaxed=TRUE,
  parallelize=FALSE)
L0glm.metrics(fit_multi, beta_true=multi$beta_true)[
  c("gic", "TP", "FP", "FN", "TN", "F1_score", "MCC",
    "within_row_F1_score", "within_row_MCC")]
# The within-row scores assess channel sparsity among correctly recovered rows.
dim(coef(fit_multi))                  # original predictors by response channels
par(mfrow=c(2, 1))
plotmat(multi$beta_true, xlab="Spike location", ylab="Channel",
  main="True multichannel coefficients")
plotmat(coef(fit_multi), xlab="Spike location", ylab="Channel",
  main="Estimated multichannel coefficients")
par(old_par)

# Precluster into at most 32 channel groups, select on their summed counts,
# then refine and return coefficients for every original mass channel.
system.time(fit_multi_clustered <- L0glm.fit(multi$X, multi$y,
  family=poisson(identity), lambda="gic", nonnegative=TRUE,
  grouped=TRUE, relaxed=TRUE, parallelize=TRUE, ncores=4L,
  precluster=TRUE, cluster.refit=TRUE,
  control.precluster=list(k.max=32L, rowbins=5L)))
fit_multi_clustered$preclustering$clusters
fit_multi_clustered$preclustering$seconds
L0glm.metrics(fit_multi_clustered, beta_true=multi$beta_true)[
  c("gic", "TP", "FP", "FN", "TN", "F1_score", "MCC",
    "within_row_F1_score", "within_row_MCC")]
dim(coef(fit_multi_clustered))         # all original predictors and channels

# A smooth nonnegative background: append a low-df spline basis and keep
# those columns unpenalized. Poisson backgrounds add directly to the signal.
B <- splines::bs(seq(0, 1, length.out=nrow(multi$X)), df=6, intercept=TRUE)
alpha <- outer(c(2, 5, 3, 7, 4, 2), seq(.5, 1.5, length.out=ncol(multi$y)))
set.seed(3)
Y_background <- multi$y + matrix(rpois(length(multi$y), B %*% alpha),
  nrow(multi$y), ncol(multi$y))
X_background <- cbind(multi$X, Matrix::Matrix(B, sparse=TRUE))
compulsory <- ncol(multi$X) + seq_len(ncol(B))
fit_background <- L0glm.fit(X_background, Y_background,
  family=poisson(identity), lambda="gic", nonnegative=TRUE,
  always.include=compulsory, grouped=TRUE, relaxed=TRUE,
  precluster=TRUE, ncores=4L,
  control.precluster=list(k.max=32L, rowbins=5L))
background_hat <- B %*% as.matrix(coef(fit_background)[compulsory, , drop=FALSE])
# Every original mass channel is retained in both the peak and background fit.
dim(coef(fit_background))
fit_background$nuisance$seconds
L0glm.metrics(fit_background,
  beta_true=rbind(multi$beta_true, Matrix::Matrix(alpha, sparse=TRUE)))[
  c("gic", "TP", "FP", "FN", "TN", "F1_score", "MCC",
    "within_row_F1_score", "within_row_MCC")]

# High-dimensional Gaussian regression with positive and negative effects.
# Run manually: the dense design alone occupies about 4 GB.
if (FALSE) { # \dontrun{
n <- 500L
p <- 1000000L
set.seed(123)
X_hd <- matrix(rnorm(n*p), n, p)
support <- sample.int(p, 60L)
beta_hd <- numeric(p)
beta_hd[support] <- runif(60L, -.8, .8)
mean_hd <- as.numeric(X_hd[, support, drop=FALSE] %*% beta_hd[support])
set.seed(1)
y_hd <- mean_hd + rnorm(n)

# screen="auto" activates ISIS exactly when p_eligible > 10*n. Screened fits
# default to the original-p GIC preset and a relaxed refit.
# Screening threads are independent of BLAS threads (nthreads).
system.time(fit_hd <- L0glm.fit(X_hd, y_hd, screen.nthreads=8, nthreads=1))
fit_hd$screening[c("method", "screen_size", "iterations", "best_iteration", "nthreads")]
fit_hd$screening$history
L0glm.metrics(fit_hd, beta_true=beta_hd)[
  c("gic", "TP", "FP", "FN", "TN", "F1_score", "MCC")]
dim(coef(fit_hd))                    # all p rows, including screened-out zeros
# Explicit alternatives: screen="SIS", "ISIS", or "none".
} # }

# Basic controls: intercept constraints, dispersion and binomial fits.
x <- seq(-1, 1, length.out=30)
X <- matrix(x, ncol=1)

# Nonnegativity includes the Gaussian intercept.
constrained <- L0glm.fit(X, -2+3*x, family=gaussian(),
  intercept=TRUE, nonnegative=TRUE, lambda=1e-8,
  dispersion=1, parallelize=FALSE)
as.numeric(constrained$coefficients)

# Estimate Gaussian variance from an automatic or user-supplied start.
set.seed(1)
y_gaussian <- 2*x + rnorm(length(x), sd=.5)
automatic <- L0glm.fit(X, y_gaussian, family=gaussian(),
  lambda="gic", dispersion=NULL, maxit=100, parallelize=FALSE)
initialized <- L0glm.fit(X, y_gaussian, family=gaussian(),
  lambda="gic", dispersion=-.25, maxit=100, parallelize=FALSE)
# dispersion=.25 instead keeps the variance fixed at .25.

# A binomial intercept may remain negative with a nonnegative slope.
y <- round(1000*plogis(-2+x))/1000
logistic <- L0glm.fit(X, y, family=binomial(),
  weights=rep(1000, length(y)), nonnegative=TRUE,
  lambda=1e-8, parallelize=FALSE)
as.numeric(logistic$coefficients)

# Gaussian validation-based stopping: keep a separate test set for assessment.
set.seed(42)
X_validation <- matrix(rnorm(200*40), 200, 40)
y_validation <- 3 + 2*X_validation[, 1] - X_validation[, 2] + rnorm(200)
train <- 1:120
fit_stopped <- L0glm.fit(X_validation[train, ], y_validation[train],
  intercept=TRUE, lambda="gic", relaxed=TRUE, maxit=100,
  validation=list(X=X_validation[-train, ], y=y_validation[-train]),
  parallelize=FALSE)
fit_stopped$validation$selection[c("best_iteration", "iterations", "stopped")]

# Optionally select the ISIS batch on the same validation set.
# Omit nsis to keep the GIC-based default batch.
fit_batch <- L0glm.fit(X_validation[train, ], y_validation[train],
  intercept=TRUE, lambda="gic", screen="ISIS", nsis=c(1, 3, 5, 10),
  relaxed=TRUE, maxit=100, validation=list(
    X=X_validation[-train, ], y=y_validation[-train], rule="1se"),
  parallelize=FALSE)
fit_batch$validation[c("selected_start", "selected_nsis", "deviance")]
fit_batch$validation$starts
# Sparse special families: GIC preset fits, GIC tuning and three-fold CV.
if (requireNamespace("mgcv", quietly=TRUE)) {
  set.seed(204)
  n <- 300L
  X_special <- matrix(rnorm(n*20), n, 20)
  folds_special <- sample(rep(1:3, length.out=n))
  family_par <- par(mgp=c(2, .7, 0))

  # Negative binomial: two nonzero slopes among 20 predictors.
  mu_nb <- exp(.4 + .8*X_special[, 1] - .6*X_special[, 2])
  y_nb <- rnbinom(n, mu=mu_nb, size=2)
  # nb() estimates theta; nb(theta=2) would keep it fixed at two.
  nb_fit <- L0glm.fit(X_special, y_nb, family=mgcv::nb(),
    intercept=TRUE, lambda="gic", relaxed=TRUE, maxit=100, tol=1e-8)
  nb_fit$theta
  coef(nb_fit)
  nb_gic <- L0glm.tune(X_special, y_nb, family=mgcv::nb(),
    intercept=TRUE, tune="gic", nlambda=5,
    relaxed=TRUE, maxit=100, tol=1e-8)
  nb_cv <- L0glm.tune(X_special, y_nb, family=mgcv::nb(),
    intercept=TRUE, tune="cv", nlambda=5, nfolds=3, foldid=folds_special,
    relaxed=TRUE, maxit=100, tol=1e-8)
  c(GIC=nb_gic$lambda_best, CV=nb_cv$lambda_best, CV_1se=nb_cv$lambda_1se)
  predict(nb_cv$fit_1se, X_special[1:5, , drop=FALSE], type="response")

  # Three-category multinomial: only two predictor rows affect probabilities.
  eta <- cbind(0, .9*X_special[, 1] - .7*X_special[, 2],
    -.8*X_special[, 1] + .6*X_special[, 2])
  probabilities <- exp(eta - apply(eta, 1, max))
  probabilities <- probabilities/rowSums(probabilities)
  y_multinomial <- factor(vapply(seq_len(n), function(i)
    sample.int(3, 1, prob=probabilities[i, ]), integer(1)), levels=1:3)
  multinomial_fit <- L0glm.fit(X_special, y_multinomial,
    family=mgcv::multinom(K=2), lambda="gic", grouped=TRUE,
    relaxed=TRUE, maxit=200, tol=1e-8)
  coef(multinomial_fit)               # original-scale, zero-sum CLR rows
  multinomial_gic <- L0glm.tune(X_special, y_multinomial,
    family=mgcv::multinom(K=2), tune="gic", nlambda=5, grouped=TRUE,
    relaxed=TRUE, maxit=200, tol=1e-8)
  multinomial_cv <- L0glm.tune(X_special, y_multinomial,
    family=mgcv::multinom(K=2), tune="cv", nlambda=5,
    nfolds=3, foldid=folds_special, grouped=TRUE,
    relaxed=TRUE, maxit=200, tol=1e-8)
  c(GIC=multinomial_gic$lambda_best, CV=multinomial_cv$lambda_best,
    CV_1se=multinomial_cv$lambda_1se)
  predict(multinomial_cv$fit_1se, X_special[1:5, , drop=FALSE], type="response")

  # Right-censored Cox model: two nonzero log-hazard effects, no intercept.
  if (requireNamespace("survival", quietly=TRUE)) {
    event_time <- rexp(n, rate=exp(.7*X_special[, 1] - .5*X_special[, 2]))
    censor_time <- rexp(n, rate=.4)
    y_cox <- survival::Surv(pmin(event_time, censor_time),
      event_time <= censor_time)
    cox_fit <- L0glm.fit(X_special, y_cox, family=mgcv::cox.ph(),
      lambda="gic", relaxed=TRUE, maxit=100, tol=1e-8)
    coef(cox_fit)
    cox_gic <- L0glm.tune(X_special, y_cox, family=mgcv::cox.ph(),
      tune="gic", nlambda=5, relaxed=TRUE, maxit=100, tol=1e-8)
    cox_cv <- L0glm.tune(X_special, y_cox, family=mgcv::cox.ph(),
      tune="cv", nlambda=5, nfolds=3, foldid=folds_special,
      relaxed=TRUE, maxit=100, tol=1e-8)
    c(GIC=cox_gic$lambda_best, CV=cox_cv$lambda_best, CV_1se=cox_cv$lambda_1se)
    predict(cox_cv$fit_1se, X_special[1:5, , drop=FALSE], type="risk")
  }
  par(family_par)
}

Fit generalized linear model with \(L_0\) or grouped row-\(L_0\) regularization subject to optional nonnegativity constraints