Fast nonnegative identity-link Poisson L0 or grouped L0 penalized fit

`L0nnpois()` is a fast RcppEigen/OpenMP implementation of [L0glm.fit()] for the nonnegative Poisson model with identity link. It shares the main fitting arguments with `L0glm.fit()` for use in deconvolution workflows, and internally uses the sparse or dense Eigen Cholesky LLT ridge solver. `L0glm.fit(openmp = TRUE)` selects this backend automatically for eligible fits. Compared with the general R fitting path (`openmp = FALSE`), it keeps the design matrix `X` in shared memory across outcome channels. Dense inputs retain dense storage and sparse inputs retain sparse storage. The dense solver forms weighted crossproducts using R's BLAS library. Eligible wide active systems use the Woodbury identity to factor an observation-space system rather than a predictor-space system. The default penalty is `lambda = "gic"`. Optional `precluster=TRUE` screens summed channel groups and refines on the original responses; see [preclustering] for aggregate-only output and clustering controls. Smooth IRLS updates use the same fixed-penalty backtracking safeguard as `L0glm.fit()`, followed by the selection cutoff. Trial steps reuse predictor differences and require no additional ridge solve. Every coefficient, including any intercept, is constrained to be nonnegative. The LLT sparse-selection update clamps negative proposals to zero. Unpenalized intercepts and `always.include` terms are refined using a small constrained Poisson solve; see [compulsory_terms]. With positive active effects and well-chosen regularization such as `lambda="gic"`, ridge shrinkage often keeps genuine active effects positive. See the solver discussion in [L0glm.fit()] for measured comparisons with the formally constrained alternatives. This constraint cannot be disabled. Use `L0glm.fit(nonnegative = FALSE)` if coefficients of either sign are required.

Usage

L0nnpois(
  X,
  y,
  weights = NULL,
  family = poisson(identity),
  dispersion = NULL,
  start = NULL,
  offset = NULL,
  lambda = "gic",
  intercept = FALSE,
  always.include = NULL,
  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,
  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

The identity-link Poisson family, poisson("identity").

dispersion

Poisson dispersion is fixed at 1. Accepts NULL or 1; other finite nonzero numeric values are ignored with a warning. Negative values cannot request dispersion estimation for this family.

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; if TRUE, add an intercept column to X. The intercept is constrained to be nonnegative, just like every other coefficient. Defaults to FALSE.

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.

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 OpenMP threads to use across outcome channels when parallelize=TRUE. For L0nnpois, the requested value is capped only at the detected number of logical cores, not at the number of physical cores.

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

Optional list of solver controls. L0nnpois always uses Eigen Cholesky LLT internally, but for sparse matrices control.solver$NaturalOrdering can be set to TRUE or FALSE to choose Eigen NaturalOrdering or AMDOrdering. The default is FALSE. Integer control.solver$ordering uses 1L for natural ordering and 0L meaning AMD ordering.

verbose

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

validation

Optional held-out data and stopping controls as in L0glm.fit(). Held-out deviance is scored directly in C++, retaining OpenMP threading across channels and restoring the best selection or relaxed-refit checkpoint. Validation predictors are shared across threads.

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 the same main components as [L0glm.fit()]: `coefficients`, `weights`, `dispersion`, `init.lambda`, `lambdas`, `family`, `z`, `X`, and `y`. Original prior weights and offsets are retained separately from the working weights for scoring. Methods coef(), logLik(), AIC(), BIC(), GIC(), EBIC() and MBIC() are available; see model_criteria. print() and summary() provide compact model displays; see model_summary. predict() and fitted() return original-scale Poisson means; see model_predictions for new data and prediction offsets.

Examples

## Fast nonnegative identity-link Poisson L0 or grouped L0 penalized fit
## All coefficients, including an optional intercept, are always nonnegative.
## No nonnegative, lower or upper argument is needed or accepted.
n <- p <- 1000
npeaks <- 100

sim.pois.mc <- simul(method = 4, n = n, p = p, s = npeaks / p, seed = 2)
dim(sim.pois.mc$beta_true) # 1000 x 176 channels

## Multichannel fit. OpenMP is used across outcome channels when available.
system.time(fit.mc <- L0nnpois(
  sim.pois.mc$X,
  sim.pois.mc$y,
  lambda = "gic",
  grouped = TRUE,
  relaxed = TRUE,
  parallelize = TRUE,
  ncores = round(parallel::detectCores(logical = FALSE)/2)
)) 
# with n <- p <- 1000 and npeaks <- 100   0.56s with 4 cores on my machine
# with n <- p <- 10000 and npeaks <- 1000 9.41s with 4 cores on my machine

dim(fit.mc$coefficients)
sum(fit.mc$coefficients != 0)

## A six-coefficient nonnegative spline background, refined on every channel.
B <- splines::bs(seq(0, 1, length.out=n), df=6, intercept=TRUE)
set.seed(3)
alpha <- outer(c(2, 5, 3, 7, 4, 2), rep(1, ncol(sim.pois.mc$y)))
Y_background <- sim.pois.mc$y + matrix(rpois(length(sim.pois.mc$y), B %*% alpha),
  n, ncol(sim.pois.mc$y))
compulsory <- p + seq_len(ncol(B))
fit_background <- L0nnpois(cbind(sim.pois.mc$X, Matrix::Matrix(B, sparse=TRUE)),
  Y_background, always.include=compulsory, lambda="gic", ncores=4L,
  precluster=TRUE, control.precluster=list(k.max=32L, rowbins=5L))
background_hat <- B %*% as.matrix(coef(fit_background)[compulsory, , drop=FALSE])
fit_background$nuisance$seconds

if (FALSE) { # \dontrun{
## Compare with the general R fitting path (process-based parallelism):
system.time(fit.mc <- L0glm.fit(
  X = sim.pois.mc$X,
  y = sim.pois.mc$y,
  family = poisson(identity),
  lambda = "gic",
  nonnegative = TRUE,
  openmp = FALSE,
  grouped = TRUE,
  relaxed = TRUE,
  parallelize = TRUE,
  ncores = round(parallel::detectCores(logical = FALSE)/2)
)) 
# with n <- p <- 1000 and npeaks <- 100   11.84s with 4 cores on my machine
# with n <- p <- 10000 and npeaks <- 1000 38.55s with 4 cores on my machine

dim(fit.mc$coefficients)
sum(fit.mc$coefficients != 0)
} # }

## Single-channel fit uses the same nonnegative model.
system.time(fit.sc <- L0nnpois(
  sim.pois.mc$X,
  sim.pois.mc$y[, 1L],
  lambda = "gic",
  grouped = FALSE,
  relaxed = FALSE,
  parallelize = FALSE
)) 
# with n <- p <- 1000 and npeaks <- 100   0.03s on my machine
# with n <- p <- 10000 and npeaks <- 1000 0.13s on my machine

sum(fit.sc$coefficients != 0)

## Native validation stopping with an independent repeat measurement.
small <- simul(n = 300L, p = 300L, n_nonzero = 20L, m = 4L, seed = 7L)
set.seed(8)
y.validation <- matrix(rpois(length(small$y_true), small$y_true), 300L, 4L)
validated <- L0nnpois(small$X, small$y, lambda = "gic", ncores = 2L,
  validation = list(X = small$X, y = y.validation, patience = 3L))
validated$validation$deviance
validated$validation$selection$history

## Compare with the general R path for a single channel:
system.time(fit.sc <- L0glm.fit(
  X = sim.pois.mc$X,
  y = sim.pois.mc$y[, 1L],
  family = poisson(identity),
  lambda = "gic",
  nonnegative = TRUE,
  openmp = FALSE,
  grouped = FALSE,
  relaxed = FALSE,
  parallelize = FALSE
)) 
# with n <- p <- 1000 and npeaks <- 100   0.04s on my machine
# with n <- p <- 10000 and npeaks <- 1000 0.25s on my machine

sum(fit.sc$coefficients != 0)

Fast nonnegative identity-link Poisson L0 or grouped L0 penalized fit

Usage

Arguments

Value

See also

Examples