Skip to contents

The L0glm algorithm and its foundations

Source: vignettes/articles/algorithm.Rmd
algorithm.Rmd

Generalized linear models

A generalized linear model (GLM) specifies a response distribution, a linear predictor and a link relating that predictor to the conditional mean. For the standard exponential dispersion families, a unit-weight observation has density or probability mass function

f(y_i;\theta_i,\phi)= \exp\left\{\frac{y_i\theta_i-b(\theta_i)}{\phi}+c(y_i,\phi)\right\},

with

\mu_i=E(Y_i\mid x_i)=b'(\theta_i),\qquad \operatorname{Var}(Y_i\mid x_i)=\phi V(\mu_i),\qquad V(\mu_i)=b''(\theta_i).

Here \theta_i is the natural parameter and \phi the dispersion. The regression model is

\eta_i=x_i^\top\beta+o_i=g(\mu_i),\qquad \mu_i=h(\eta_i),\qquad h=g^{-1},

where o_i is a known offset. The link acts on the conditional mean. For a canonical link, \eta_i=\theta_i.

Family Variance function V(\mu) Dispersion \phi
Gaussian 1 \sigma^2
Poisson \mu 1
Bernoulli \mu(1-\mu) 1
Gamma \mu^2 Estimated or supplied
Inverse Gaussian \mu^3 Estimated or supplied

For binomial proportions, trial totals enter as observation weights. Standard references are McCullagh & Nelder (1989), Generalized Linear Models, 2nd edition and Hardin & Hilbe (2018), Generalized Linear Models and Extensions, 4th edition, Chapter 3. The illustrated IRLS account by Kane & Lewis connects these equations to R code and linear solvers.

Ordinary GLM fitting by IRLS

Iteratively reweighted least squares (IRLS) fits a GLM through a sequence of weighted least-squares problems. Let a_i be a known prior/case observation weight, and put d_i^{(t)}=h'(\eta_i^{(t)})=d\mu_i/d\eta_i. Each iteration forms the offset-adjusted working response and weights

\begin{aligned} z_i^{(t)} &=\eta_i^{(t)}-o_i+\frac{y_i-\mu_i^{(t)}}{d_i^{(t)}},\\ W_{ii}^{(t)} &=a_i\frac{\{d_i^{(t)}\}^2}{V(\mu_i^{(t)})}. \end{aligned}

Here the working response has the known offset subtracted, so the least-squares solve uses Xb directly. R’s family$linkinv(), family$mu.eta() and family$variance() supply h, h' and V, respectively. See R’s family-object reference.

Ordinary GLM IRLS

Input: design X, response y, family/link, observation weights a, offset o, a convergence tolerance and an iteration limit.

  1. Initialize valid means \mu^{(0)} from the family’s initialization rule and set \eta^{(0)}=g(\mu^{(0)}). When coefficient starts are supplied, instead use \eta^{(0)}=X\beta^{(0)}+o and check the implied means.
  2. For t=0,1,\ldots, calculate d^{(t)},z^{(t)},W^{(t)} above.
  3. Solve the weighted least-squares problem \beta^{(t+1)}=\arg\min_b \left\|W_t^{1/2}(z_t-Xb)\right\|_2^2. Its normal equations are (X^\top W_tX)\beta^{(t+1)}=X^\top W_tz_t.
  4. Update \eta^{(t+1)}=X\beta^{(t+1)}+o and \mu^{(t+1)}=h(\eta^{(t+1)}). Stop when the deviance or the available successive coefficient vectors change by less than the chosen tolerance, or when the iteration limit is reached.

A numerical solver obtains the least-squares solution without explicitly inverting the matrix. Practical implementations add domain checks and can shorten an unstable full step by backtracking or step halving. For Gaussian data with an identity link, z=y-o and W=\operatorname{diag}(a), so the unpenalized coefficient fit reduces to weighted linear regression.

IRLS is Fisher scoring, using expected information. For the weighted likelihood its score with respect to the linear predictor is

\frac{\partial\ell}{\partial\eta_i} =\frac{a_i}{\phi} \frac{(y_i-\mu_i)h'(\eta_i)}{V(\mu_i)}, \qquad I_\beta=\frac1\phi X^\top WX.

The corresponding quadratic approximation to negative twice log likelihood has the form

Q_t(b)=C_t+\frac1\phi\left\|W_t^{1/2}(z_t-Xb)\right\|_2^2.

The common factor 1/\phi cancels from an unpenalized coefficient solve. It matters when a coefficient penalty is added. The key modification in L0glm is to replace the weighted least-squares solve by weighted ridge regression with adaptive penalties, as follows.

L0glm uses this approach to approximate best subset selection and best subset of groups (BSGS) selection, combining sparse support penalties with the likelihood appropriate to each model family. The mathematical pseudocode appendix gives the complete fitting routines, including initialization, distribution-parameter updates and SIS/ISIS screening.

A sparse likelihood target

Let \ell(\beta,\psi) denote a model’s log likelihood, including known observation weights and offsets, and let \psi collect distribution parameters such as a dispersion or negative-binomial shape. For one response, L0 selection targets

\mathcal I_\lambda(\beta,\psi) =-2\ell(\beta,\psi)+\lambda\sum_{j\in\mathcal P}\mathbf1\{\beta_j\ne0\},

where \mathcal P contains selectable predictors. Intercepts and always.include terms are unpenalized during coefficient fitting. This is the penalized form of best subset selection: lambda controls the trade-off between likelihood and the number of selected predictors. For a coefficient matrix B, independent selection counts nonzero entries; grouped selection treats each predictor row across response channels as a group. Best subset of groups (BSGS) selection then targets

\mathcal I_\lambda(B,\psi)=-2\ell(B,\psi)+\lambda\|B\|_{2,0}.

The group count is

\|B\|_{2,0}=\sum_{j\in\mathcal P}\mathbf1\{\|B_{j\cdot}\|_2>0\}.

This objective favours a compact support without the persistent linear shrinkage of an L1 penalty. The optional relaxed refit estimates effects on the selected support with negligible ridge penalization.

IRLS with adaptive ridge

Using the same GLM working response and observation weights, L0glm replaces the ordinary coefficient update with adaptive weighted ridge:

\beta^{(t+1)}\approx\arg\min_b \left\|W^{1/2}(z-Xb)\right\|_2^2+ \sum_j\Lambda_j^{(t)}(b_j-c_j)^2, \qquad \Lambda_j^{(t)}=\frac{\lambda\phi}{4|\beta_j^{(t)}|^2}.

For an unconstrained update, put D_t=\operatorname{diag}(\Lambda_j^{(t)}) and solve

(X^\top W_tX+D_t)\beta^{(t+1)}=X^\top W_tz_t+D_tc.

Multiplying the likelihood-scale quadratic by \phi puts the solve on a weighted least-squares scale, giving the dispersion factor in \Lambda_j. The factor 1/4 calibrates the inverse-square adaptive ridge update to the L0 selection threshold, giving a useful information-criterion preset. The calibration calculation below shows this first for orthogonal Gaussian models and then for the frozen IRLS quadratics.

The exponent is fixed at two, with a default zero ridge center c_j=0. Mandatory terms have \Lambda_j=0. As a coefficient shrinks, its adaptive penalty increases rapidly. Strong coefficients receive progressively less shrinkage. The design remains dense or sparse through the appropriate solver path.

The core iteration is:

  1. Initialize from training data or supplied coefficients.
  2. Update the GLM working response, observation weights and adaptive penalties.
  3. Solve on the active support, respecting the selected solver’s constraints.
  4. Check the smooth step and apply the selection cutoff.
  5. Repeat until convergence, the iteration limit or validation stopping.
  6. Optionally refit the selected support and report likelihood and recovery metrics.

For identity-link Poisson models, z=y-o and W_{ii}=a_i/\mu_i. This specialization is implemented in native C++ with OpenMP across response channels.

Foundations and developments in L0glm

Adaptive ridge and broken adaptive ridge provide the foundation for this work. Frommlet & Nuel (2016) establish the L0 calibration used here; Dai et al. (2018) and Li et al. (2021) develop BAR and scalable GLM extensions. The variational and MM identities below are established results reviewed by Abergel, Bouaziz & Nuel (2024). The likelihood scores, information matrices, IRLS linearizations and Woodbury identity are standard statistical and linear-algebra tools.

The present work brings these ideas together and develops the following adaptations, whose complete computations are given in the appendix.

Component Foundation Development in L0glm
Automatic regularization and scale GIC of Fan & Tang (2013), with adaptive-ridge calibration Scale-aware initialization, bounded dispersion updates and a common penalty convention across families
Negative binomial Standard NB2 likelihood and shape estimation Adaptive-ridge IRLS coupled to safeguarded log-shape updates, with independent or grouped channels
Multinomial Joint softmax likelihood and its quadratic approximation Class-symmetric CLR output, invariant grouped penalties, entrywise CLR selection and contrast-space ridge solves
Cox PH BAR Cox regression of Kawaguchi et al. (2020) Cached weighted Breslow risk sets, exact matrix-free curvature products, a matrix-free ridge initializer and direct/PCG dispatch
Screening SIS of Fan & Lv (2008) and likelihood ISIS of Fan, Samworth & Wu (2009) Family-aware efficient-score screening, OpenMP implementation and an offset-based ISIS/L0glm pipeline
Deconvolution L0 pulse-train deconvolution of de Rooi & Eilers (2011) Automatic regularization, Poisson and negative-binomial noise models, and shared-support multichannel extensions

These algorithmic adaptations and their application to multichannel GC/MS peak recovery form part of the methods article in preparation. The Gaussian orthogonal-design calculation below is attributed to its source; the appendix then derives the family-specific updates used in this implementation.

Why reweighted quadratics promote sparsity

Adaptive ridge has several equivalent optimization interpretations. Abergel, Bouaziz & Nuel (2024) develop the variational, majorization–minimization and local-approximation views. Two useful identities make the connection explicit.

For 0<q<2,

|b|^q=\inf_{a>0}\left\{\frac q2ab^2+ \left(1-\frac q2\right)a^{-q/(2-q)}\right\}.

For nonzero b, the minimizing auxiliary weight is a=|b|^{q-2}. Substituting this value gives |b|^q, while fixing the auxiliary weights produces a weighted quadratic problem. This explains the alternation between penalty weights and ridge solves.

For the inverse-square rule, let s=b^2. Concavity of the logarithm gives

\log s\leq\log s_t+\frac{s-s_t}{s_t},\qquad s,s_t>0.

The tangent is exact at s=s_t and has quadratic coefficient 1/b_t^2. With a positive smoothing scale \delta, the corresponding tangent weight is 1/(b_t^2+\delta^2) for \log(b^2+\delta^2). A normalized smooth log penalty approaches a nonzero indicator as \delta decreases. L0glm uses unsmoothed weights on the active support and explicitly removes numerically small coefficients.

With fixed support, fixed distribution parameters and exact minimization of a tangent majorizer, the usual MM inequality gives

F(b_{t+1})\leq Q(b_{t+1}\mid b_t) \leq Q(b_t\mid b_t)=F(b_t).

IRLS supplies a local likelihood approximation. L0glm checks the frozen-penalty smooth step by backtracking where needed, then thresholds the accepted coefficients. Dispersion updates, changing supports and relaxed refits form the surrounding sparse-model search. The reported information criterion evaluates the likelihood and parameter count of the returned model.

A hierarchical Bayesian and ARD view

Adaptive ridge also has a useful hierarchical Bayesian interpretation. A ridge step is a conditional maximum a posteriori (MAP) estimate under Gaussian priors. In the solver’s coefficient coordinates, write

\beta_j\mid\alpha_j\sim N(0,\alpha_j^{-1}),\qquad -2\log p(\beta\mid\alpha)=\sum_j\alpha_j\beta_j^2+\text{constant}.

For a GLM dispersion \phi, multiplying the local negative twice log likelihood by \phi gives the ridge system above with \Lambda_j=\phi\alpha_j. The inverse-square update corresponds to

\alpha_j^{(t)}=\frac{\lambda}{4|\beta_j^{(t)}|^2},\qquad \operatorname{Var}(\beta_j\mid\alpha_j^{(t)}) =\frac{4|\beta_j^{(t)}|^2}{\lambda}.

Thus the default prior remains zero-centered and skeptical: small fitted effects receive increasingly concentrated priors around zero, while strong effects receive wider priors. Each IRLS step updates the coefficient mode using a local quadratic likelihood; the next adaptive step updates the prior scale from that mode. The prior mean stays fixed. A supplied prior.mean changes the ridge center.

This is closely related to automatic relevance determination (ARD) and sparse Bayesian learning: each predictor has its own prior precision, and an irrelevant predictor is removed as that precision grows. A shared precision across a coefficient block gives the corresponding group-ARD interpretation. Tipping (2001) develops the sparse Bayesian learning framework using Gaussian coefficient priors with learned precisions.

The update used here is a plug-in empirical-Bayes-style scale update based on coefficient point estimates. Evidence-maximizing ARD also uses posterior uncertainty. For example, an EM variance update in a simple Gaussian hierarchy contains E(\beta_j^2\mid y)=\widehat\beta_j^2+ \operatorname{Var}(\beta_j\mid y); the inverse-square rule uses the fitted coefficient magnitude. The IC calibration sets the proportionality constant, and hard thresholding implements the limiting concentration at zero. The relaxed refit then estimates the retained effects with negligible shrinkage.

One hierarchical construction

The inverse-square rule can be obtained as the zero-rate limit of a conditional MAP update. For independent priors \alpha_j\sim\operatorname{Gamma}(a,b), with shape a>1/2 and rate b>0, the conditional log density is

\log p(\alpha_j\mid\beta_j) =\left(a-\tfrac12\right)\log\alpha_j -\left(b+\tfrac12\beta_j^2\right)\alpha_j+\text{constant}.

Its mode is (2a-1)/(\beta_j^2+2b). Choosing a=1/2+\lambda/8 gives \alpha_j=\lambda/\{4(\beta_j^2+2b)\}, which approaches the implemented active-set rule as b\downarrow0. The rate-zero limit has an improper hyperprior; a positive rate yields a smoothed inverse-square update. L0glm uses the unsmoothed rule on retained coefficients and prunes small coefficients explicitly. This construction explains the hierarchical connection while preserving the package’s likelihood-and-information-criterion formulation.

The factor of one quarter

The calibration comes from the orthogonal Gaussian analysis of Frommlet & Nuel (2016). It can also be seen directly. Let s_j=x_j^\top x_j, known variance \phi, and least-squares coefficient \widehat b_j. Exact L0 selection keeps predictor j when

\frac{s_j\widehat b_j^2}{\phi}>\lambda.

The inverse-square ridge iteration is

b_{t+1}=\frac{s_j\widehat b_j} {s_j+\lambda\phi/(4b_t^2)}.

A nonzero fixed point therefore satisfies

b^2-\widehat b_jb+\frac{\lambda\phi}{4s_j}=0.

Its discriminant is nonnegative exactly when s_j\widehat b_j^2/\phi\geq\lambda: the same selection boundary. The relaxed fit removes the remaining ridge shrinkage. This calibration also works well approximately for non-orthogonal designs, as supported by the simulations of Frommlet & Nuel (2016) and the examples here. If desired, L0glm.tune() efficiently searches a grid around the preset to select the regularization level more accurately for the data at hand, using the actual information criterion or cross-validated deviance.

Why the same factor applies to IRLS quadratics

The weighted least-squares structure gives a direct motivation for using the same calibration across GLM families. With the working response, observation weights and dispersion held fixed, every IRLS subproblem has the form

Q_t(b)=C_t+\frac1\phi\|W_t^{1/2}(z_t-Xb)\|_2^2.

For weighted-orthogonal columns, put s_j=x_j^\top W_tx_j, h_j=s_j/\phi, and let \widehat b_{j,t}=x_j^\top W_tz_t/s_j be the unpenalized working least-squares coefficient. Profile out unpenalized nuisance terms first when they are present. The coordinate contribution is h_j(b_j-\widehat b_{j,t})^2. Comparing its best nonzero value with zero gives the L0 rule

h_j\widehat b_{j,t}^{\,2}>\lambda.

The adaptive-ridge update for this same quadratic is

b_{j,\mathrm{new}}= \frac{h_j\widehat b_{j,t}}{h_j+\lambda/(4b_{j,\mathrm{old}}^2)},\qquad b_j^2-\widehat b_{j,t}b_j+\frac{\lambda}{4h_j}=0.

Its nonzero fixed-point boundary is again h_j\widehat b_{j,t}^{\,2}=\lambda. The working weights and dispersion change the curvature h_j, while the factor one quarter stays the same. This calculation holds for each frozen quadratic, whatever the family or link that generated its valid working quantities. The training-data initializer supplies the first working point; subsequent IRLS steps adapt the quadratic to the fitted means.

For correlated columns, holding the other coefficients fixed still gives a one-dimensional quadratic with this calibration. Joint support decisions also involve the off-diagonal curvature terms. This explains both the useful transfer of the quarter-factor preset and the benefit of optional tuning when the design has strong dependence.

Large-sample likelihood justification

The connection to the original likelihood becomes more accurate as its local quadratic approximation improves. On a regular, identifiable candidate model, expand around an unpenalized local likelihood maximum \widetilde\beta:

-2\ell(\widetilde\beta+d) =-2\ell(\widetilde\beta)+d^\top I_nd+R_n(d), \qquad I_n=-\nabla^2\ell(\widetilde\beta).

If the third derivatives are bounded in norm by a constant times n in the neighbourhood, Taylor’s theorem gives |R_n(d)|\leq Cn\|d\|^3. With information eigenvalues of order n, an effect at the selection boundary has scale \|d\|=O(\sqrt{\lambda/n}), so

R_n(d)=O\!\left(\frac{\lambda^{3/2}}{\sqrt n}\right).

For fixed penalty this error vanishes as n grows. For a growing IC penalty, \lambda^3/n\to0 is one sufficient condition for an absolute error tending to zero in this fixed-dimensional calculation. Fisher scoring uses expected information; under regularity, observed and expected information agree to first order. These arguments motivate the GLM extension on well-behaved active models. With growing active dimension, bounds must also control that dimension and the restricted information matrix. Increasing sample size improves local likelihood accuracy, while predictor correlations can persist. The preset and optional path tuning address these two aspects together.

Thresholds and starts

During selection the default rule removes coefficients satisfying

|b_j|<\max_k|b_k|\,\texttt{thresh} \quad\text{or}\quad |b_j|<10^{-20}, \qquad\texttt{thresh}=10^{-5}.

Thresholding uses the solver’s normalized coefficient coordinates where scaling is enabled. A zero coefficient has an infinite adaptive penalty and remains excluded from that selection pass. This shrinking active set saves substantial work. Final relaxed refitting applies no hard cutoff. Explicit zero starts restrict eligible coefficients in advance, with mandatory terms retained.

Default initialization is training-based. Estimated dispersions use a permissive ridge pilot and bounded updates, with inexpensive random-sign probes to guard against an almost-interpolating pilot. Positive supplied dispersions are held fixed; negative values initialize further estimation at their absolute value.

Grouped and special-family updates

Ordinary multichannel GLMs fit channels on a shared support and repeat shared-support selection passes. Row starts use a signed mean-absolute magnitude for signed models, avoiding cancellation across channels; nonnegative models use arithmetic row averages. Relaxed grouped fitting admits the eligible channel coefficients within each selected row.

Multinomial regression uses a joint likelihood, zero-sum CLR output and joint quadratic updates. Grouped penalties act on class-symmetric row norms; ungrouped selection acts on individual CLR entries. Negative-binomial fits update the shape parameter alongside coefficient fitting. Cox fits use a quadratic approximation to the Breslow partial likelihood, with cached risk sets and matrix-free Hessian-vector products for larger active systems.

Scaling the linear algebra

Weighted ridge has the normal equations

(X^TWX+D)b=X^TWz+Dc.

For wide systems with positive diagonal D, the Woodbury identity gives

b=c+D^{-1}X^T(W^{-1}+XD^{-1}X^T)^{-1}(z-Xc),

factoring an observation-space system instead of a predictor-space system. Zero-weight observations are omitted and unpenalized terms are profiled out. Dense direct-Gram and sparse augmentation/factorization paths exploit the appropriate storage. Required solver-boundary conversions retain formal constraints for OSQP; the native Poisson backend shares designs across threads.

For ultrahigh dimensions, SIS or offset-based ISIS adds another level of dimension reduction. Compiled screening ranks candidates, and joint refitting updates the support in small batches. This follows the screening ideas of Fan & Lv (2008) and Fan, Samworth & Wu (2009).

Information criteria and validation

The GIC preset is \lambda=\log p\log\log n, motivated by Fan & Tang (2013). L0glm.tune() can instead select candidate fits by an IC or held-out deviance. CV returns both the minimum-deviance and one-SE full-data models.

Reported GIC() uses the model’s parameter dimension: ordinary multichannel models count coefficient entries, and grouped multinomial rows contribute K-1 parameters. A row-count grouped target and a parameter-count reported criterion therefore have distinct complexity terms. The criterion reference records the conventions. Validation can also select the best IRLS checkpoint, including its fitted distribution parameters and penalties.

Appendix: mathematical pseudocode

The following algorithms describe the computations used by L0glm.fit() and L0nnpois(). They share preprocessing, ridge solves and stopping routines. Here C denotes the number of multinomial categories, m the number of ordinary response channels, and q the number of design columns after adding any intercept. The selectable predictor count is p.

Jump to initialization, single-channel GLMs, multichannel GLMs, negative binomial, multinomial, Cox PH, SIS, ISIS, or validation stopping.

A.1. Coordinates, eligibility and the ridge solve

Write M for mandatory terms: an intercept and any always.include columns. Let s_j>0 be a column scale and h_j a column shift, which may be represented implicitly when profiling a sparse design. In working coordinates,

Z_{ij}=(X_{ij}-h_j)/s_j,\qquad b_j=s_j\beta_j,\qquad b_0=\beta_0+\sum_j h_j\beta_j.

The intercept equation applies when a unit intercept is free. The same affine change of coordinates is applied to prior.mean, giving the ridge center c, or C_0 for a coefficient matrix. Predictions retain their original meaning: Zb+o=X\beta+o. Output transforms the coefficients back and restores excluded rows and entries as zeros in the original dimensions.

  • normalize="auto" uses RMS scaling for signed fits or a free intercept, and no scaling for fully nonnegative fits with a constrained intercept. RMS scaling gives \|Z_j\|_2=\sqrt n for selectable columns, using centered column norms when a free intercept is available. "l2" gives unit norm; "sd" uses sample standard deviations; "none" gives s_j=1.
  • Mandatory terms have scale one. Centering is enabled with normalization and a free unit intercept. Sparse designs retain implicit zeros; intercept profiling supplies the equivalent shift for columns whose explicit centering would fill the matrix.
  • A scalar start=1 means one in working coefficient coordinates. Other numeric starts use original input units. For ordinary GLMs and NB, explicit zero entries are excluded, except in mandatory terms. Multinomial starts exclude rows that are zero across all categories before CLR centering; category entries within an eligible row remain eligible.
  • nonnegative=TRUE constrains coefficients and any intercept to be nonnegative. Binomial and quasibinomial intercepts remain free. Multinomial CLR coefficients are unrestricted; Cox fitting omits an intercept.

Define the reusable operation

\mathcal R(Z,z,W,D,c):\quad (Z^T WZ+D)\widetilde b=Z^TWz+Dc.

Eigen solves this system and clips to the specified bounds. The constrained glmnet/bigGlm coordinate-descent and OSQP quadratic-programming backends impose bounds during optimization. For nonnegative models with positive active effects and well-chosen regularization, such as lambda="gic", ridge shrinkage often keeps genuine active effects positive, and Eigen’s simple clamping works very well. In the package’s deconvolution solver benchmarks, Eigen Cholesky LLT generally gives better support recovery and coefficient accuracy than the formally constrained alternatives, with much faster fitting.

Infinite diagonals in D remove entries from the solve; mandatory terms have diagonal zero. Free nuisance terms can be profiled by a Schur complement. For a wide system, with A=W^{1/2}Z, r=W^{1/2}z and positive penalty diagonals,

\widetilde b=c+D^{-1}A^T (I+AD^{-1}A^T)^{-1}(r-Ac).

Thus the smaller coefficient or observation system is factored. Dense and sparse designs use their respective kernels. The normalized ordinary-GLM solver may store unscaled coefficients internally; dividing its original-unit ridge diagonal by s_j^2 gives exactly the working-coordinate diagonal used below.

Unless supplied explicitly, the main iteration controls are maxit=25, tol=1e-10, thresh=1e-5 and group.maxit=6. Resolve character penalties using the selected information criterion. For GIC, \lambda=\log(p)\log\log(n_{\rm IC}): ordinary GLMs use the number of rows, NB and multinomial use the number of positive-weight rows, and Cox uses the weighted number of events. The function defaults are lambda=1 for L0glm.fit() and lambda="gic" for L0nnpois().

A.2. Ordinary-family initialization and dispersion

The ordinary GLM path evaluates family$initialize and obtains \mu^{\rm init}. The family object also supplies g, g^{-1}, \mu'=d\mu/d\eta, the variance function V, domain checks and unit deviance. These are the standard R family operations. For a numeric response in one channel:

Family V(\mu) Initial mean \mu_i^{\rm init} Scale \phi
Gaussian 1 y_i Supplied or estimated
Poisson \mu y_i+0.1 1
Binomial \mu(1-\mu) (a_i y_i+0.5)/(a_i+1) 1
Gamma \mu^2 y_i>0 Supplied or estimated
Inverse Gaussian \mu^3 y_i>0 Supplied or estimated
Quasipoisson \mu y_i+0.1 Supplied or estimated
Quasibinomial \mu(1-\mu) (a_i y_i+0.5)/(a_i+1) Supplied or estimated
quasi() Chosen variance function Defined by its initialize method Supplied or estimated

Zero-weight binomial responses are set to zero during initialization. Other ordinary family objects use their own initialization expression. For a multichannel response, this initialization is applied to each numeric channel.

Set \eta^{\rm init}=g(\mu^{\rm init}), then form

W_i^0=a_i\frac{\mu'(\eta_i^{\rm init})^2}{V(\mu_i^{\rm init})}, \qquad z_i^0=\eta_i^{\rm init}-o_i+ \frac{y_i-\mu_i^{\rm init}}{\mu'(\eta_i^{\rm init})}.

Identity links simplify to z^0=y-o. Keep rows with finite positive working weight and finite working response.

Coefficient seed with normalization enabled. Let A=(W^0)^{1/2}Z, r=(W^0)^{1/2}z^0. Put the unrestricted mandatory columns in U, and define P_U=U(U^TU)^\dagger U^T, using a rank-aware solve. For each column compute

v_j=A_j^T(I-P_U)A_j,\quad t_j=A_j^T(I-P_U)r,\quad \rho_j=\begin{cases}0,&j\in M,\\\lambda,&j\notin M,\end{cases} \qquad b_j^0=\frac{t_j+\rho_j c_j}{v_j+\rho_j}.

Uninformative columns with zero target receive zero. Set the coefficients of U aside while calculating the slopes, then recover them by regressing the remaining working response on U. For Gaussian models, first rescale the nonfree seed by

f=\frac{v^T(I-P_U)r}{v^Tv},\qquad v=(I-P_U)A b^0,

using f=0 when the denominator is zero. This calibrates its joint predicted magnitude. Finally clip the seed to the coefficient bounds.

Coefficient seed with normalize="none". Use the diagonal working-ridge seed

b_j^0=\frac{\sum_i W_i^0 X_{ij}z_i^0+\lambda c_j} {\sum_i W_i^0 X_{ij}^2+\lambda}.

Here the initialization ridge is applied to all columns; mandatory terms become unpenalized in the main iteration. For binomial/quasibinomial models with a unit intercept, compute each slope conditional on that intercept. With d_j=\sum_iW_i^0X_{ij}^2+\lambda, r_j=\sum_iW_i^0X_{ij}z_i^0+\lambda c_j and t_j=\sum_iW_i^0X_{ij}, replace its seed by (r_j-t_jr_0/d_0)/(d_j-t_j^2/d_0). For Gaussian models multiply the whole seed by \sum_iW_i^0z_i^0(Xb^0)_i/\sum_iW_i^0(Xb^0)_i^2, taking zero when its denominator vanishes. Clip to the bounds. L0nnpois() uses the Poisson identity-link specialization of these diagonal formulas.

Supplied coefficient starts replace this coefficient-seeding calculation. The first working response and weights still use the family’s initialized means. Subsequent working means come from the accepted fitted predictor.

Subroutine D: initialize and update an unknown scale

Input: the initial working model, mandatory columns, and the dispersion setting. Output: \phi and a flag recording whether it may be updated.

  1. For a fixed-scale family set \phi=1, or its declared fixed scale. Otherwise a positive supplied value fixes the scale. A negative supplied value starts estimation at its absolute value. With NULL, run the following pilot.

  2. Form A=(W^0)^{1/2}X and r=(W^0)^{1/2}z^0. With normalization enabled, profile the mandatory block and divide informative selectable columns by their conditional weighted norms. With normalize="none", divide all informative columns by their weighted norms. Fit a pilot with ridge \tau=10^{-6}, using primal or Woodbury algebra.

  3. Write R_\tau=I-H_\tau for the pilot’s residual operator, including nuisance profiling where applicable. Four reproducible Rademacher probes v_\ell, with entries \pm1, share the same factorization. Estimate

    d_2=\frac14\sum_{\ell=1}^4\|R_\tau v_\ell\|^2, \qquad \phi_{\rm pilot}=\frac{\|R_\tau r\|^2}{d_2}.

    The probe seed is 1701; the caller’s random-number state is restored.

  4. For a normalized fit, guard the residual information. Let n_w be the number of retained working rows and d_*=\min\{10,(n_w-r_M)/5\}, with r_M the nuisance rank. If d_2<d_*, also estimate d_4=\frac14\sum_\ell\|R_\tau^2v_\ell\|^2 and \nu=d_2^2/d_4. A pilot with d_2\ge d_* already passes the cheap lower bound. Retry once with \tau=0.1 if the scale is invalid, below 10^{-6} times the nuisance-only working scale, or \nu<d_*.

  5. Initialize permissively at

    \phi^0=\max\left\{\mathrm{DBL\_MIN}, \frac{\phi_{\rm pilot}}{\max(1,\lambda_{\rm GIC,ref})}\right\}.

    Replace an invalid/nonpositive pilot scale by machine epsilon first. The reference GIC uses the preprocessing’s selectable count and the good-row count for normalized fits. Unnormalized fits use the supplied design-column count and input-row count. It is independent of the requested fitted lambda.

  6. Let n_+=\#\{i:a_i>0\}. During selection, from iteration 3 onward, update only if the number of active or mandatory coefficients satisfies k_t<0.7n_+. At the accepted current mean calculate

    \phi_P=\frac{\sum_{i:a_i>0}a_i(y_i-\mu_i)^2/V(\mu_i)}{n_+-k_t}, \qquad \phi^{t+1}=\operatorname{clip}(\phi_P,\phi^t/2,2\phi^t).

    Three successive updates with unchanged support and |\phi_P/\phi^t-1|<0.01 freeze the scale. A support change unfreezes an estimated scale. Invalid/nonpositive targets leave its value unchanged.

  7. A relaxed ordinary-GLM refit holds the current scale during its solves and reports the final Pearson scale afterwards. A supplied positive scale stays fixed. The same final reporting update is used for an unpenalized or empty estimated-scale fit.

A.3. Single-channel ordinary GLMs

Algorithm 1: adaptive-ridge GLM on one response

Input: X,y,a,o,\lambda,c,M, bounds, an optional start and validation set. Output: original-scale coefficients, working quantities, dispersion and convergence diagnostics.

  1. Apply A.1, restrict to explicitly eligible columns, initialize b^0, \mu^{\rm init}, \eta^{\rm init} and \phi^0 by A.2.

  2. For t=0,\ldots,T-1: form the working model from the initialized means if t=0, or from the accepted fitted means otherwise:

    W_i^t=a_i\frac{\mu'(\eta_i^t)^2}{V(\mu_i^t)},\qquad z_i^t=\eta_i^t-o_i+\frac{y_i-\mu_i^t}{\mu'(\eta_i^t)}.

    Keep finite informative rows. For an identity link use z^t=y-o directly. Inverse-Gaussian log-link calculations use the stable equivalents W_i=a_i e^{-\eta_i} and z_i=\eta_i-o_i-1+y_i/\mu_i.

  3. Set the active set to the nonzero coefficients together with M, and form the frozen adaptive diagonals

    D_{jj}^t=\begin{cases} 0,&j\in M,\\ \lambda\phi^t/(4(b_j^t)^2),&b_j^t\ne0,\ j\notin M,\\ +\infty,&b_j^t=0,\ j\notin M. \end{cases}

    If there are no active terms, return the offset-only model. Otherwise obtain \widetilde b=\mathcal R(Z,z^t,W^t,D^t,c). A free unit intercept receives the residual weighted-least-squares correction \sum_iW_i(z_i-(Z\widetilde b)_i)/\sum_iW_i.

  4. Apply the smooth-step routine below. For an accepted selection step, remove each nonmandatory coefficient satisfying

    |b_j|<\texttt{thresh}\max_{k\notin M}|b_k| \quad\text{or}\quad |b_j|<10^{-20}.

    Keep the thresholded point when it has finite family deviance; otherwise retain the accepted smooth point. Refresh \eta=Zb+o and \mu=g^{-1}(\eta). Poisson, quasipoisson, Gamma and inverse-Gaussian working means have a machine-epsilon lower guard.

  5. Update an estimated scale by Subroutine D. If validation is supplied, score and retain the checkpoint using A.10.

  6. Stop after at least two iterations when the coefficient-change measure below is less than tol and scale updating is complete, or when scale estimation is unavailable at the current large support after iteration 3. Unpenalized and relaxed passes can stop without freezing the scale. Also stop at the iteration limit or the validation patience limit.

  7. If relaxed=TRUE, restrict to the selected terms together with M, start from their fitted coefficients and repeat with \lambda=10^{-30} and hard thresholding disabled. Restore the original coefficient dimensions.

Smooth-step routine. Freeze D, the ridge center and current nuisance parameters. Let

F_t(b)=d_f\{y,g^{-1}(Zb+o);a\}+ \sum_{j:D_{jj}<\infty}D_{jj}(b_j-c_j)^2,

where d_f is the family’s weighted unit deviance. After the first iteration, when the smooth proposal preserves support, use F_t(b^t) as reference; at initialization or a support-changing proposal require a finite trial. Try b^t+2^{-h}(\widetilde b-b^t), h=0,\ldots,50, accepting the first finite value no larger than the reference plus 10^{-10}\max(1,|F_t(b^t)|). Predictor interpolation reuses matrix products. If no trial is accepted, retain the last finite point and stop the pass. When needed, the initial feasible comparison point is found by halving the coefficient start towards a valid offset-only model, up to 60 times.

Coefficient convergence. Without normalization use \|\beta^{t+1}-\beta^t\|_2/\max(10^{-12},\|\beta^{t+1}\|_2). With normalization use the same ratio for selectable working coefficients, and, when mandatory terms exist, also the prediction-change ratio

\frac{\|\widehat f^{t+1}-\widehat f^t\|_a} {\max\{10^{-12},\sqrt{\sum_i a_i}, \|\widehat f^{t+1}-\overline f_a^{t+1}\|_a\}}, \qquad \widehat f=Zb,\quad \|v\|_a^2=\sum_i a_iv_i^2.

Take the larger change. For a Gaussian identity-link dense direct solve, the weighted Gram matrix is cached while the active support shrinks. L0nnpois() specializes Algorithm 1 to nonnegative Poisson identity-link fitting: \phi=1, W_i=a_i/\mu_i, z_i=y_i-o_i, with Eigen LLT, Woodbury where appropriate and OpenMP across channels.

A.4. Ordinary multichannel GLMs

Algorithm 2: channel fitting with shared-support passes

Input: X and Y\in\mathbb R^{n\times m}, channel weights/offsets, grouped and the controls of Algorithm 1. Output: a sparse q\times m coefficient matrix on the original scale.

  1. Prepare each channel’s starting mask and scale state. A row is eligible if any channel permits it; mandatory rows stay eligible. A shared weight or offset vector is replicated across channels.

  2. If grouped=FALSE, run Algorithm 1 once per channel, with each channel’s own seed, weights and dispersion. Proceed to the optional refit.

  3. If grouped=TRUE, run up to group.maxit shared-support passes. In the first pass use the supplied starts, or each channel’s automatic initializer. Fit the eligible channels independently by Algorithm 1, in parallel where requested, and assemble B.

  4. Before each further pass retain the row union S=M\cup\{j:\exists c,\ B_{jc}\ne0\}. For each retained row set

    r_j=\begin{cases} m^{-1}\sum_c B_{jc},&\text{nonnegative fit},\\ \operatorname{sign}(\overline B_j)\,m^{-1}\sum_c|B_{jc}|, &\text{signed fit}. \end{cases}

    Pool the returned, original-unit coefficient rows. For a signed row with zero arithmetic mean, take the sign of its largest absolute entry. Copy r_j into every eligible channel of that row, preserve the supplied zero masks, and rerun Algorithm 1 on S.

  5. Pass an estimated channel dispersion to the next pass as its negative value, allowing further optimization from the previous estimate. A fixed channel dispersion stays positive. Validation can stop these passes and restore their best checkpoint.

  6. For a relaxed refit, pool the retained rows once more and refit eligible entries using \lambda=10^{-30}, with thresholding disabled. Without grouped selection, the one-pass row union is used for this refit as well; originally excluded starting entries remain excluded. Return all q rows.

The common rows and pooled starts couple the ordinary channel fits. Each inner adaptive diagonal still uses that channel’s coefficient. Grouped NB and multinomial fitting below instead share a squared-row-norm diagonal within each IRLS iteration.

A.5. Negative binomial with estimated shape

Use family=mgcv::nb() with variance V(\mu)=\mu+\mu^2/\theta. Each response channel has its own shape \theta_c, while scale dispersion is one. The family accepts log, identity or square-root links. Its shape-setting convention is retained: nb() starts an estimated shape at one; nb(theta=-v) starts at v and estimates it; nb(theta=v) fixes it at v>0.

Algorithm 3: joint coefficient and NB-shape updates

  1. Prepare Z,C_0,M, bounds and eligible entries as in A.1. Initialize each \theta_c from the family object.

  2. With start=NULL, begin from B=0 and take one NB working-ridge step with diagonal 10^{-3} on selectable terms and zero on M. At this step calculate \eta=ZB+O, \mu=g^{-1}(\eta). If a channel has invalid/nonpositive means, use \mu_i=y_i+0.1 and \eta_i=g(\mu_i) for its working model. Supplied starts replace this pilot, retaining their entrywise eligibility restrictions.

  3. Repeat until convergence or maxit: recompute \eta,\mu and form, for every channel,

    W_{ic}=a_{ic}\frac{\mu'(\eta_{ic})^2} {\mu_{ic}+\mu_{ic}^2/\theta_c}, \qquad z_{ic}=\eta_{ic}-o_{ic} +\frac{y_{ic}-\mu_{ic}}{\mu'(\eta_{ic})}.

    Use the same positive-mean fallback when required. Set selectable adaptive diagonals to

    D_{jc}=\begin{cases} \lambda/(4\sum_hB_{jh}^2),&\text{grouped, }m>1,\\ \lambda/(4B_{jc}^2),&\text{ungrouped or single-channel}. \end{cases}

    Set mandatory diagonals to zero and excluded entries to infinity. At \lambda=0, set all permitted diagonals directly to zero. Solve each channel’s working ridge problem and apply the special-family acceptance and threshold rules below.

  4. Holding the accepted means fixed, take one bounded log-shape update in each channel whose shape is estimated, using the equations below.

  5. Stop when the larger of the relative selectable-coefficient and fitted-predictor changes is below tol, and \max_c|\log(\theta_c^{\rm new}/\theta_c)|<10^{-6}, or by the iteration/validation limit. Special-family relative changes divide by \max(1,\|B_{\rm new,selectable}\|_F) and \max(1,\|ZB_{\rm new}\|_F), respectively.

  6. For relaxed=TRUE, set \lambda=0, disable thresholding and refit. Grouped refits admit eligible entries in selected rows; ungrouped refits retain selected entries. Estimated shapes continue updating in the refit.

Log-shape step. For one channel let a=\log\theta, and write \psi, \psi_1 for the digamma and trigamma functions. At fixed means,

\begin{aligned} U_\theta={}&\sum_i a_i\left[ \psi(y_i+\theta)-\psi(\theta)+\log\theta+1 -\log(\theta+\mu_i)-\frac{y_i+\theta}{\theta+\mu_i}\right],\\ H_\theta={}&\sum_i a_i\left[ \psi_1(y_i+\theta)-\psi_1(\theta)+\frac1\theta -\frac2{\theta+\mu_i}+\frac{y_i+\theta}{(\theta+\mu_i)^2}\right],\\ G_a={}&\theta U_\theta,\qquad H_a=\theta^2H_\theta+\theta U_\theta. \end{aligned}

Use d=-G_a/H_a when H_a<0 is finite, otherwise d=\operatorname{sign}(G_a), and clip d to [-1,1]. Try a_{\rm trial}=\operatorname{clip}(a+2^{-h}d,\log10^{-6},\log10^7), up to 25 halvings. Accept the first finite NB log likelihood at least as large as its value at a; otherwise retain \theta. These bounds and the shape tolerance can be changed through control.family.

Special-family acceptance and thresholding. NB, multinomial and Cox use -2\ell+\sum D_{jc}(B_{jc}-C_{0,jc})^2 as the smooth-step objective, with distribution parameters and D frozen. Use the current objective as reference when the proposal preserves support, and a finite-value requirement when support changes; allow up to 30 halvings. Grouped NB/multinomial selection thresholds the row norm r_j=\|B_j\|_2 against thresh times the largest selectable row norm and 10^{-20}. Ungrouped and single-channel selection applies the scalar rule within each channel. Thresholding is applied during positive-\lambda selection, with the smooth point retained if the thresholded likelihood is invalid. In a grouped pass, an eligible zero entry within a nonzero row has a finite shared penalty and remains available. An all-zero selectable row has infinite penalty during positive-\lambda selection.

A.6. Joint multinomial IRLS in CLR coordinates

For mgcv::multinom(K=C-1), convert labels to a one-hot response matrix or use the supplied count matrix Y. Let N_i=\sum_cY_{ic}, with one case weight a_i per row. Store B\mathbf1_C=0, center offsets and prior centers across categories, and compute probabilities by stable softmax:

P_{ic}=\frac{\exp(\eta_{ic}-\max_h\eta_{ih})} {\sum_h\exp(\eta_{ih}-\max_k\eta_{ik})}, \qquad \eta=ZB+O.

The mgcv family specifies the category count; L0glm implements its joint fitting steps in class-symmetric CLR coordinates.

Algorithm 4: grouped or entrywise-CLR multinomial selection

  1. Define eligible rows from the original supplied starts before centering. Center supplied coefficient rows to CLR, preserving excluded rows as zero. With start=NULL, all entries are eligible and B^0=0.

  2. For an automatic start take one joint multinomial IRLS ridge step with selectable diagonal 10^{-3} and mandatory diagonal zero, from the softmax probabilities at B=0 and the supplied offsets. For a supplied start, seed any eligible zero-norm row in a grouped fit, or eligible zero entry in an ungrouped fit, from this weak-ridge step and recenter. Thus scalar start=1 centers to zero and receives the joint pilot.

  3. Repeat: use observations with a_iN_i>0 and form the current score and curvature blocks

    U_c=Z^T\{a\odot(Y_c-N\odot P_c)\},\qquad H_{cd}=Z^T\operatorname{diag} \{a_iN_iP_{ic}(\mathbf1\{c=d\}-P_{id})\}Z.

    Use D_{jc}=\lambda/(4\|B_j\|_2^2) for grouped selection and D_{jc}=\lambda/(4B_{jc}^2) for entrywise CLR selection, with the usual mandatory and excluded-entry rules.

  4. Parameterize the allowed zero-sum coefficient space by an orthonormal basis K, with v=K^T\operatorname{vec}(B), u=K^T\operatorname{vec}(U), H_K=K^THK and D_K=K^T\operatorname{diag}(\operatorname{vec}D)K. Solve

    (H_K+D_K)\widetilde v=H_Kv+u+ K^T\operatorname{diag}(\operatorname{vec}D)\operatorname{vec}(C_0), \quad \operatorname{vec}(\widetilde B)=K\widetilde v.

    Only allowed finite-penalty entries enter this equation. For shared row penalties, K=Q\otimes I, where Q^TQ=I and Q^T\mathbf1=0; its penalty stays diagonal. Entrywise selection uses a separate zero-sum basis on each row’s surviving categories, with ridge whitening to retain diagonal-ridge kernels.

  5. Apply special-family backtracking and row-norm or entrywise thresholding. Recenter surviving entries within each allowed row; a row with fewer than two allowed entries becomes zero. Keep the unthresholded smooth point if thresholding makes the likelihood invalid.

  6. Use the special-family coefficient/predictor convergence measures in Algorithm 3, or the iteration/validation limits. For a relaxed grouped fit, refit all originally eligible categories within retained rows with \lambda=0. An ungrouped refit keeps only selected CLR entries. No hard cutoff is applied during the refit. Return full-size CLR coefficients and category probabilities.

The default curvature is the exact joint multinomial IRLS curvature. The optional control.family=list(multinomial="bound") uses a_iN_i I/2 on the C-1 contrast plane. Small working systems use coefficient-space normal equations; wide systems whiten the per-observation contrast curvature and use Woodbury in n(C-1) dimensions. Those wide curvature factorizations guard probabilities at 10^{-12} and renormalize them. The default max.working=2000 guards the smaller working dimension.

A.7. Cox proportional hazards

Take Surv(time,status), with status \delta_i\in\{0,1\}, case weights a_i, and \eta_i=z_i^Tb+o_i. At each distinct event time t_r, let R_r=\{i:t_i\ge t_r\}, d_r=\sum_{i:t_i=t_r}a_i\delta_i, e_i=a_i\exp(\eta_i) and S_r=\sum_{i\in R_r}e_i. The weighted Breslow log partial likelihood is

\ell(b)=\sum_i a_i\delta_i\eta_i-\sum_r d_r\log S_r.

mgcv::cox.ph() selects this family in L0glm. Its response is the explicit time/status pair, and weights remain case weights. Sorting once by decreasing time makes risk-set sums cumulative sums. Cache the ordering, tied-time groups, group endpoints and weighted event counts; recompute only the predictor-dependent risk quantities. Subtracting a common maximum from \eta stabilizes the exponentials without changing the score or partial likelihood.

Define

q_i=e_i\sum_{r:i\in R_r}\frac{d_r}{S_r},\quad u_i=a_i\delta_i-q_i,\quad V_\eta=\operatorname{diag}(q)-\operatorname{diag}(e)\,C\,\operatorname{diag}(e), \quad C_{ih}=\sum_{r:i,h\in R_r}\frac{d_r}{S_r^2}.

Algorithm 5: sparse Cox fitting with direct or matrix-free ridge steps

  1. Omit the intercept, apply A.1, fix scale dispersion at one, and cache the risk-set bookkeeping. Supplied starts retain their eligibility masks.

  2. With start=NULL, evaluate q^0,u^0 at b=0, including offsets. Obtain a weak diagonal-working-curvature ridge pilot

    (Z^T\operatorname{diag}(q^0)Z+10^{-3}I)\widetilde b=Z^Tu^0.

    Use matrix-free preconditioned CG, including for this initializer. For q>\#\{i:q_i^0>0\}, solve its observation-space Woodbury system. Clip to bounds and backtrack from zero against -2\ell(b)+10^{-3}\|b\|^2, allowing 30 halvings and requiring an actual nonincrease. Supplied starts replace this pilot.

  3. Repeat: evaluate e,S,q,u and set D_{jj}=\lambda/(4b_j^2) on active selectable terms, zero on M, and infinity on excluded terms. Solve for an increment

    (Z^TV_\eta Z+D)\Delta=Z^Tu-D(b-c),\qquad \widetilde b=b+\Delta.

  4. For a small dense unrestricted active system, form its information matrix and use the requested direct kernel, defaulting to Eigen LLT. In automatic mode this requires k\le1024, k\le n, k\le\texttt{max.working}, and k\times\#\{r:d_r>0\}\le10^7. Sparse designs default to the matrix-free path; cox.llt.max can change these cutoffs.

  5. Otherwise use preconditioned CG on the exact operator below. Use diagonal preconditioning \operatorname{diag}(Z^T\operatorname{diag}(q)Z)+D, with a verified scaled residual tolerance of 10^{-10} and at most 2000 CG iterations by default. Automatic mode can fall back to an eligible direct system if CG fails. Coefficient bounds or a constrained solver use the direct working-system path; cox="direct" explicitly selects it.

  6. Apply special-family backtracking and scalar thresholding. Stop using the special-family convergence measure, iteration limit or validation rule. For relaxed=TRUE, repeat on selected terms with \lambda=0 and no thresholding. Return the full coefficient vector; estimate baseline cumulative hazard by \widehat H_0(t)=\sum_{r:t_r\le t}d_r/S_r.

The exact matrix-free product, requiring design products and cumulative risk-set passes, is

\begin{aligned} w&=Zv,\\ \{(Z^TV_\eta Z+D)v\} &=Z^T\left[q\odot w-e\odot \left\{\sum_{r:i\in R_r}\frac{d_r}{S_r^2} \sum_{h\in R_r}e_hw_h\right\}_{i=1}^n\right]+Dv. \end{aligned}

The initializer uses \operatorname{diag}(q^0) as its inexpensive pilot curvature. Main CG iterations use the full partial-likelihood curvature above. The direct wide path can whiten the observation information and call the same ridge kernels. The Cox help page gives the partial-likelihood and prediction conventions.

A.8. Family-aware SIS

SIS ranks each candidate by its efficient score conditional on a nuisance model. The null fit includes the supplied offset. The default screening nuisance block is an intercept for ordinary, NB and multinomial responses, and empty for Cox; standalone sis_glm() also accepts a nuisance design. This screening nuisance block has its own role alongside the joint model’s intercept and always.include terms.

Algorithm 6: native efficient-score screening

  1. Fit the small, unpenalized nuisance model with offset o. Ordinary GLMs use their family initialization and nuisance GLM fit; estimate a Pearson scale for a free-scale family. NB estimates/fixes its shape according to the family setting. Multinomial uses a joint zero-sum nuisance fit; Cox uses its weighted partial likelihood and cached risk sets.

  2. For an ordinary or NB channel form

    s_i=\frac{a_i(y_i-\mu_i)\mu'_i}{\phi V(\mu_i)},\quad W_i=\frac{a_i(\mu'_i)^2}{\phi V(\mu_i)},\quad G=(N^TWN)^\dagger,\quad u_N=N^Ts.

    Here N is the nuisance design, \phi=1 for Poisson, binomial and NB, and NB uses its fitted shape in V.

  3. For each originally eligible predictor x_j, compute

    A_j=N^TWx_j,\qquad U_j=x_j^Ts-A_j^TGu_N,\qquad I_j=x_j^TWx_j-A_j^TGA_j.

    Its score is T_j=U_j^2/I_j. A nonnegative alternative replaces U_j by \max(0,U_j). For independent response channels add their channel statistics. Treat numerically uninformative directions as zero.

  4. For multinomial use the C-1 joint score/information blocks from Algorithm 4, project out the joint nuisance block by the same Schur complement, and set T_j=U_j^TI_j^\dagger U_j. For Cox use s=u and the full observation information V_\eta, so its scalar statistic uses x_j^TV_\eta x_j and the corresponding nuisance cross-information. Remove directions below the relative 10^{-12} information tolerance before taking these inverses.

  5. Keep the largest nsis statistics. Deterministic index tie breaking gives reproducible results. OpenMP workers scan independent columns, with thread-local top-set heaps; only requested full rankings store all scores.

This construction makes screening invariant to nonzero rescaling of a predictor for a two-sided alternative; a positive rescaling also preserves the one-sided alternative. Finite positive information and estimable nuisance directions are handled by rank-aware inverses. For an intercept-only identity-Poisson screening model with difficult offsets, the nuisance score \sum_i a_i\{y_i/(o_i+b_0)-1\}=0 is solved directly when required.

A.9. Combined SIS/ISIS and L0glm

Algorithm 7: offset screening followed by joint adaptive ridge

Input: the original design, family, masks and fitting controls; a batch size b; one round for SIS or up to screen.maxit=10 rounds for ISIS.

  1. Build the original eligible predictor set E, respecting starting-zero exclusions, and retain the mandatory set M. Define p_{\mathrm{eligible}}=|E\setminus M|, counting a selectable predictor when at least one response channel permits it. With screen="auto", use ISIS exactly when p_{\mathrm{eligible}}>10n; otherwise fit directly with the appropriate Algorithm 1–5 and return. screen="none" also fits directly; explicit SIS/ISIS requests apply at any design size. Here n is the number of training rows in the current fit. Resolve lambda once using the original selectable predictor count. Initialize full-size B=0 and the selected set to the mandatory predictors.
  2. For screening round r=1,\ldots,R, form the screening offset O_r=O+XB_{r-1}, including any fitted intercept. Run Algorithm 6 on the original eligible columns with this offset and retain its top b predictors, T_r.
  3. Define the candidate set C_r=S_{r-1}\cup T_r\cup M. If it equals the preceding candidate set, stop. Otherwise fit the appropriate Algorithm 1–5 on X_{C_r} with the original offset O, the resolved full-problem lambda, and screen="none". Each joint fit uses a fresh automatic start, or the user’s supplied starts restricted to these rows. Screening defaults to a relaxed joint fit and serial response-channel fitting unless overridden.
  4. Expand this fit to the original coefficient dimensions, with zero rows outside C_r, and set S_r to its nonzero row support together with the mandatory predictors. Record its likelihood, parameter dimension, GIC and optional validation score.
  5. Continue ISIS rounds until the candidate set is unchanged or the round limit is reached. Retain the round with smallest original-problem GIC, or smallest validation loss when supplied. Return its full-size sparse coefficients, support, screening history and timings.

Each ISIS round scans all originally eligible predictors. A predictor omitted from one joint fit remains eligible for a later screening round. The current fitted predictor enters the next screening offset; every joint fit uses the original model offset.

For L0glm.fit() the batch/retained-set defaults are

\begin{aligned} b_{\rm SIS}&=\max\{1,\ 2\lfloor n/\log(\max(n,3))\rfloor\},\\ b_{\rm ISIS}&=\max\left\{1,\ \operatorname{round} \left[\frac{n}{\log(\max(p,2))\log\log(\max(n,3))}\right]\right\}, \end{aligned}

capped at the available predictor count. Standalone sis_glm() defaults to \max\{1,\lfloor n/\log(\max(n,3))\rfloor\}. An explicit nsis overrides these sizes. The automatic rule p_{\mathrm{eligible}}>10n has no separate minimum sample size or predictor count. L0glm.tune() uses the same dispatch independently for each full-data and CV training fit, with n equal to that fit’s training row count.

GIC comparison uses the original p and the fitted parameter count, including estimated distribution parameters where appropriate. For a grouped multinomial fit, an active row contributes C-1 free category directions; for ordinary multichannel models, nonzero channel entries contribute separate parameters. This follows the reported criterion conventions. Screening scores themselves require no fit of a full p-dimensional model.

A.10. Validation checkpoints and returned quantities

When validation is supplied, training data alone determine initializers, normalization, working weights, screening statistics and ridge updates. First fit a mandatory-only baseline. At each accepted iteration evaluate held-out mean deviance and save an improved checkpoint when

L_{\rm val}^{\rm new}<L_{\rm val}^{\rm best} -\epsilon_{\rm val}\max\{1,|L_{\rm val}^{\rm best}|\}.

Here \epsilon_{\rm val} is validation$tol. The defaults are tolerance 10^{-8}, minimum iteration 3 and patience 3. Stop when the minimum iteration has been reached and the patience count is exhausted. Restore the best coefficients and nuisance parameters, including dispersion or NB shape, plus their support/penalty state. Ordinary/native working quantities are refreshed at a restored checkpoint. Grouped ordinary models also track validation across shared-support passes. A relaxed pass can retain its selection checkpoint as its starting best fit.

Ordinary GLMs use weighted mean unit deviance; NB with estimated shape uses predictive log loss relative to a fixed Poisson saturated reference; fixed-shape NB uses the same fixed NB saturated reference throughout. Multinomial uses deviance per trial and Cox partial deviance per event, computed in validation risk sets. Positive-mean validation scoring has a default floor of 10^{-8}. These scoring conventions keep the reference comparable while fitted distribution parameters change.

For Gaussian identity-link validation, the default candidate starts are the automatic seed and an uncalibrated conditional-marginal seed. Supplied starts override this comparison. The marginal seed uses b_j=t_j/v_j from A.2, without the initialization ridge or joint Gaussian calibration, and recovers the mandatory coefficients from the residual weighted fit. An nsis vector adds a screening-batch grid; include the mandatory-only model among the candidates. Select minimum validation loss, or the sparsest candidate within one SE when validation$rule="1se". For per-observation validation losses l_i and normalized response masses \omega_i, the candidate SE is \{\sum_i\omega_i^2(l_i-\overline l)^2/(1-\sum_i\omega_i^2)\}^{1/2}, taking zero when only one positive mass is available. L0glm.tune(tune="cv") wraps the same fitting routines in training folds; each fold computes its own initializers and screening.

Returned coefficients always use the original predictor units and original dimensions. Model objects also retain the distribution/link, offsets, prior weights, normalization, selected lambda, convergence information, fitted distribution parameters and validation/screening diagnostics when requested.

Nonnegative compulsory terms and spline backgrounds

Separate the unpenalized block from the sparse slopes, \eta_j=X\beta_j+B\alpha_j+o_j,\qquad\alpha_j\in\mathcal C_j. Here \mathcal C_j imposes nonnegativity on compulsory terms, with the usual unrestricted binomial intercept. Conditional on \beta_j, solve the small constrained likelihood problem for \alpha_j. An active-set quadratic solver uses exact Gaussian curvature, observed Poisson curvature for nonnegative identity-link signals, or a safeguarded Fisher-scoring quadratic for other scalar GLMs. Cox uses its partial-likelihood score and information. Negative-binomial shape updates remain conditional on the current fitted means.

The outer procedure is:

  1. Initialize a full model. For automatic starts in positive-mean identity models, a nonnegative partition-of-unity spline basis is initialized through a constant fit, then expanded into its spline coefficients.
  2. Solve the small compulsory block, keeping the sparse signal fixed.
  3. Refit candidate slopes with B\alpha included in the link-scale offset, using the original penalty and predictor universe. Use two grouped passes unless group.maxit is supplied.
  4. Accept an update when the original-data likelihood plus complexity improves; validation stopping uses held-out deviance instead. Repeat at most twice and finish with a compulsory-block polish.

Compulsory parameters are included in dispersion degrees of freedom even when their contribution is temporarily represented as a fixed offset. Explicit zero starts continue to exclude selectable coefficients. A one-iteration ridge pilot (maxit=1) retains its single-pass definition. Preclustering supplies Poisson candidates; all outer scores and final coefficients use the original channels. See compulsory terms for details.

Implementation map. The equations above can be followed directly in the source:

Routine Source
Coordinates, scaling and mandatory terms normalization.R, mandatory.R
Ordinary single-channel IRLS L0glm_1_fit_internal.R
Multichannel passes and refitting L0glm_2_fit.R, group_start.R
Coefficient seeds and dispersion pilot initialization.R, dispersion.R
Constrained compulsory blocks nuisance.R, nuisance.cpp
Native Poisson backend L0nnpois.cpp
Special-family passes families_common.R, family_nb.R, family_multinomial.R, family_cox.R, family_cox_pcg.R
Screening and outer ISIS sis_context.R, sis.cpp, isis_glm.R, screening.R
Backtracking and validation irls.R, validation.R

References