Estimators

HighDimensionalModels

Lasso estimators computed by coordinate descent: lasso, sqrt_lasso, scaled_lasso, feasible_lasso!, LassoPath, the penalized generalized linear models glm_lasso and GLMLassoPath, the group lasso group_lasso and GroupLassoPath, the multi-response lasso multi_lasso and MultiLassoPath, the multitask lasso with task-specific designs multitask_lasso, MultitaskLassoPath and CVMultitaskLassoPath, the inference after selection of R's hdm (rlasso, post_lasso, double_selection), and the varying-coefficient models varying_coefficient and varying_coefficient_lasso, with the kernels of KernelSmoothing.jl. They run the coordinate descent, the penalty paths and the reweighting of ProximalSolvers.jl on the losses and penalties of ProximalBase.jl, as does coordinate_descent!, which minimizes f(x) + g(x) for a loss f (a CoordinateDifferentiableFunction or a loss that ProximalSolvers knows) and a penalty g that separates across coordinates, such as the weighted ℓ1 norm NormL1. Low-rank multivariate fits use reduced_rank, nuclear_regression, sparse_reduced_rank, sparse_svd_regression, sparse_nuclear_regression, and the additive lowrank_plus_sparse, with LowRankPath, CVLowRankPath, PenaltyGrid and CVPenaltyGrid for tuning. The lasso solvers accept real data and penalties of any type and convert them to a common floating point type: the element type of X and y promoted to a float, or that of the iterate passed to scaled_lasso! and feasible_lasso!.

Penalties on ordered coefficients and graphs are fitted by fused_lasso, trend_filter and generalized_lasso, using exact proximal steps, proximal gradient or specialized ADMM rather than coordinate descent. The dirty multitask model dirty_multitask decomposes coefficients into shared and task-specific components, with DirtyMultitaskPath and CVDirtyMultitaskPath for tuning their two penalty levels. Joint-sparse transfer from source domains to a target, TMTL(Fused) and TMTL(Debiased) of Lim and Kolar, is fitted by transfer_multitask and debias, with TransferMultitaskPath and CVTransferMultitaskPath for paths and target-only tuning.

LassoPath(X, y, λpath, options=CDOptions(); max_hat_s=Inf, standardize_x=true, intercept=false, screening=true,
          gram=:auto, weights=nothing, offset=nothing, alpha=1, penalty=nothing, lla=nothing, lower=-Inf,
          upper=Inf, rng=Random.default_rng(), kwargs...)
LassoPath(X, y, options=CDOptions(); nλ=100, λminratio, kwargs...)

Lasso solutions β0path[k] and βpath[k] of

minimize  ‖y - β₀ - Xβ‖² / (2n) + λpath[k] ∑ⱼ ωⱼ |βⱼ|

along a decreasing sequence λpath, each one warm started from the previous solution, by ProximalSolvers.solve_path. The intercept β₀ is not penalized, and it is zero unless intercept is true. The loadings are the standard deviations ωⱼ = sqrt(∑ᵢ (Xᵢⱼ - x̄ⱼ)² / n) of the columns when standardize_x is true (centered only with an intercept) and ωⱼ = 1 otherwise.

Without λpath, the path has nλ values spaced evenly on the log scale from λmax, the smallest value whose solution is zero, down to λminratio ⋅ λmax; λminratio defaults to 0.01 when X has more columns than rows and to 1e-4 otherwise. The first value is λmax times 1 + √eps, so that rounding cannot make the first solution nonzero.

With screening, each λ is solved first over the coordinates kept by the sequential strong rule of Tibshirani et al. (2012), and the optimality conditions of the others are checked afterwards; coordinates that violate them are added and the problem solved again. kkt_checks[k] is the number of checks at λpath[k]: one when the rule kept every coordinate it should have, zero without screening. When options.warm_start is false, each value is solved from zero along its own path of options.num_steps levels, without screening.

The path stops early, after the first solution with more than max_hat_s nonzero coefficients; the returned λpath then contains only the values that were solved for. The solutions are stored as SparseVectors, so a path takes memory proportional to the number of nonzero coefficients. iterations[k] and converged[k] are the passes of coordinate descent at λpath[k] and whether it met its stopping rule. See also refit_lasso_path.

The keyword arguments kwargs that are not listed are fields of CDOptions that replace those of options, and the random order of coordinate descent is drawn from rng. With gram, the path is solved on the Gram matrix of the columns, and with weights and offset the loss is the weighted one, as described for lasso; λmax and the loadings are weighted too. fit(LassoPath, X, y, ...) is the same as LassoPath(X, y, ...). The functions of StatsAPI apply to the path: coef, predict and nobs; nobs is the number of observations, also stored as nobs.

alpha and penalty give the elastic net and other penalties, as for lasso: LassoPath(X, y; alpha=0.5) is glmnet's elastic-net path at a fixed α, whose λmax is that of the lasso divided by α, and LassoPath(X, y; penalty=MCP(1, 3)) the path of ncvreg, on standardized columns with standardize_x. Along a path of a nonconvex penalty, each solution is the local one reached from the previous solution, starting from zero at λmax. Without standardize_x, the exact steps of MCP and SCAD need the condition on a described for lasso. lla computes the path by local linear approximation, as described there. With penalty = CappedL1(1, a) and intercept = true (and lla = nothing, which then means LLA(; stages = 100)) it is the path of glmtlp's "tlp" with tau = a, with standardize_x as glmtlp's standardize. With lla = LLA(; stages = 1000, start = :previous), standardize_x = false and intercept = true, the paths of MCP, SCAD and CappedL1 are those of ncpen with x.standardize = FALSE (ncpen's own standardization differs from standardize_x). Its λmax is that of the first stage, the lasso with the weights pⱼ′(0)/λ. penalty holds the penalty at level one, NormL1(1, ω) for the lasso with loadings ω; the solution at λpath[k] minimizes the loss plus that penalty at level λpath[k], applied to scales .* β when scales is not nothing (standardize_x with a penalty other than the ℓ1 norm). lower and upper bound the coefficients, as for lasso; penalty is then Constrained(g, IndBox(l, u)), and λmax is the smallest level whose solution is zero under the bounds.

LassoSolution

Result returned by lasso, sqrt_lasso, scaled_lasso and feasible_lasso!.

Fields

• x::SparseIterate: the estimated coefficients
• residuals::Vector: y - intercept - X * x
• penalty: the penalty at which x was computed, a NormL1 for the ℓ1 estimators, or the penalty of lasso's alpha or penalty at the level λ; with lower or upper, that penalty on the box of the bounds, Constrained(g, IndBox(l, u))
• σ: an estimate of the noise level (see the documentation of each solver), or nothing
• iterations::Int: passes of coordinate descent over the coordinates, summed over all the lasso fits the solver made
• converged::Bool: whether coordinate descent met its stopping rule within max_iter passes and, for scaled_lasso! and feasible_lasso!, the outer iteration met its stopping rule within its max_iter iterations
• intercept: the unpenalized intercept, zero unless the solver was called with intercept = true
• hasintercept::Bool: whether the solver fitted an intercept
• weights: the observation weights of the fit, nothing for none
• scales: nothing when penalty applies to x, or the scales s of the columns when it applies to s .* x, the coefficients of the standardized columns X[:, j] / s[j] (standardize_x with a penalty other than the ℓ1 norm)

The functions of StatsAPI apply to it: coef, predict, residuals, nobs, dof and deviance.

feasible_lasso!(x, X, y, λ0, options=IterLassoOptions(); intercept=false, rng=Random.default_rng(), kwargs...)

Heteroscedasticity-robust lasso of Belloni, Chen, Chernozhukov and Hansen (2012), which uses data-driven penalty loadings. It alternates between the lasso

β̂ = argmin  ‖y - β₀ - Xβ‖² / (2n) + λ0 ∑ⱼ Γⱼ |βⱼ|

and the update of the loadings Γⱼ = sqrt(∑ᵢ (Xᵢⱼ - x̄ⱼ)² rᵢ² / n) from the residuals r = y - β̂₀ - Xβ̂, until the relative change in Γ is below options.opt_tol (by ProximalSolvers.reweight!), or warns after options.max_iter iterations. With max_iter = 0 there is no update: x is left as the initial procedure leaves it, the loadings are those of its residuals, and converged is false, without a warning. The intercept β₀ is not penalized; it and the column means x̄ⱼ are zero unless intercept is true.

The initial residuals are chosen by options.init_procedure (see IterLassoOptions); with :InitStd they come from a lasso with penalty λ0 σinit ∑ⱼ sⱼ |βⱼ|, sⱼ = sqrt(∑ᵢ (Xᵢⱼ - x̄ⱼ)² / n). Coordinate descent starts from and overwrites x, and each lasso fit is warm started from the previous one. Returns a LassoSolution whose penalty.w holds the loadings of the last fit and whose σ is the sample standard deviation of the residuals. The keyword arguments kwargs are fields of IterLassoOptions that replace those of options, as in feasible_lasso!(x, X, y, λ0; max_iter=50), and the random order of coordinate descent is drawn from rng.

lasso(X, y, λ, options=CDOptions(); intercept=false, standardize_x=false, gram=:auto, weights=nothing,
      offset=nothing, alpha=1, penalty=nothing, lla=nothing, lower=-Inf, upper=Inf, rng=Random.default_rng(),
      kwargs...)
lasso(X, y, λ, ω, options=CDOptions(); intercept=false, gram=:auto, weights=nothing, offset=nothing,
      lower=-Inf, upper=Inf, rng=Random.default_rng(), kwargs...)

Solve the lasso

minimize  ‖y - β₀ - Xβ‖² / (2n) + λ ∑ⱼ ωⱼ |βⱼ|

by coordinate descent. The intercept β₀ is not penalized, and it is zero unless intercept is true. In the first form the loadings are ωⱼ = 1, or with standardize_x the standard deviations sqrt(∑ᵢ (Xᵢⱼ - x̄ⱼ)² / n) of the columns, with x̄ⱼ the column mean when intercept is true and zero otherwise; the coefficients are on the scale of the data either way. Returns a LassoSolution whose σ is the sample standard deviation of the residuals.

The keyword arguments kwargs are fields of CDOptions that replace those of options, as in lasso(X, y, λ; max_iter=10_000, opt_tol=1e-8), and the random order of coordinate descent is drawn from rng.

X can be any matrix and y any vector. A Matrix, or a view of a range of its rows or of some of its columns, is read in place. A SparseMatrixCSC stays sparse; coordinate descent still reads its columns entry by entry, so for now it is slower than the dense matrix. Any other matrix, such as a transpose or a view of selected rows, is copied once into a Matrix (a SparseMatrixCSC if it is sparse). The same holds for every solver of this package.

With gram = true, coordinate descent works on the Gram matrix of the (centered) columns, A = X'X / n, instead of the residuals. Each coordinate step then costs O(p) instead of O(n), after O(np²) operations and p² numbers of memory to form A. The loss, the stopping rule and so the solution are the same, to the tolerance. Forming A pays off over the many passes of a path when n ≥ p, but rarely at one λ. The default, gram = :auto, uses it for LassoPath on a dense X with n ≥ p and p ≤ 2000, and never at one λ.

With observation weights w (nonnegative, with a positive sum) and an offset o, the loss is

∑ᵢ wᵢ (yᵢ - oᵢ - β₀ - Xᵢ'β)² / (2 ∑ᵢ wᵢ).

Dividing by ∑ᵢ wᵢ follows glmnet: λ means the same with and without weights, unit weights give the unweighted fit, integer weights the fit of replicated rows, and zero weights leave observations out. The intercept, the loadings of standardize_x (sqrt(∑ᵢ wᵢ (Xᵢⱼ - x̄ⱼ)² / ∑ᵢ wᵢ)) and the Gram matrix are weighted as well. The residuals are y - o - β₀ - Xβ, and with weights σ is their weighted root mean square, sqrt(∑ᵢ wᵢ rᵢ² / ∑ᵢ wᵢ).

Other penalties

alpha < 1 gives glmnet's elastic net, λ (α ∑ⱼ |βⱼ| + (1 - α) ‖β‖² / 2), for 0 < α ≤ 1. penalty gives another penalty of ProximalBase whose level λ replaces: NormL1(1, w) and NormL1Nonneg(1, w) (the lasso and the nonnegative lasso with weights w), ElasticNet(1, λ2, w) (λ ∑ⱼ wⱼ|βⱼ| + λ2 ‖β‖²/2, with λ2 fixed), Scaled(ElasticNet(α, 1 - α, w), 1) (the elastic net with weights on its ℓ1 part), and the nonconvex MCP(1, a), SCAD(1, a), CappedL1(1, a) (λ min(|βⱼ|, a), the capped ℓ1 penalty of T. Zhang, 2010, and the truncated ℓ1 penalty of Shen, Pan and Zhu, 2012) and LogPenalty(1, ε) (λ log(1 + |βⱼ|/ε)), whose parameters are kept. The level of the penalty given is not used: MCP(1, 3) means "MCP with a = 3".

alpha solves the objective of glmnet's documentation. glmnet's gaussian family scales y to unit standard deviation internally, which divides the ridge part by sd(y), the population standard deviation of y (weighted with weights). So glmnet's solutions are those of penalty = Scaled(ElasticNet(α, (1 - α) / sd(y)), 1), and they equal those of alpha when sd(y) = 1.

With standardize_x, the penalty applies to the coefficients of the standardized columns X[:, j] / sⱼ, sⱼ their standard deviations, as in glmnet and ncvreg; the coefficients returned are those of the columns of X. For the ℓ1 norm this is the same as the loadings ωⱼ = sⱼ, which is how it is computed. For the other penalties X is divided by s once (a copy), and s is stored as the solution's scales.

A nonconvex penalty has local solutions, and the one found depends on the algorithm and the starting point. As in ncvreg, MCP and SCAD are computed by exact coordinate steps along options.num_steps levels from λmax, the smallest level with a zero solution, down to λ, each solved to the tolerance and warm started from the previous one; num_steps = 0 solves at λ from zero. Use LassoPath for the solutions at many levels. With options.warm_start false, a convex penalty other than the ℓ1 norm is also solved along such a path.

The exact steps need each coordinate problem to be convex: a vⱼ > 1 for MCP and (a - 1) vⱼ > 1 for SCAD, where vⱼ is the curvature of the loss along column j, its mean square ∑ᵢ (Xᵢⱼ - x̄ⱼ)² / n (centered with an intercept, weighted with weights, and of the standardized column with standardize_x). The path throws an ArgumentError when a column fails it, since the step could then leave zero where the test of λmax and screening says that zero is optimal. With standardize_x, vⱼ = 1, so ncvreg's conditions a > 1 and a > 2 suffice; without it, a column of variance at most 1/a fails. Take standardize_x = true, a larger a, or lla.

lla computes a nonconvex penalty by local linear approximation (LLA, or reweighted ℓ1) instead: lla = LLA(; stages = 2, tol = 1e-12, start = :lasso) (ProximalSolvers' LLA, exported) solves a sequence of weighted lassos with the weights pⱼ′(|βⱼ|)/λ at the previous solution, each a convex problem, by the coordinate descent of options; the inner solver of the LLA given is not used. From zero the first stage is the lasso. The stages stop after stages, or once the weights change by at most tol, relative to the largest one.

• LLA(), two stages, is the two-step estimator of Fan, Xue and Zou (2014); with MCP and SCAD it has no condition on a.
• LLA(; stages = 100) iterates until the weights stop changing, as the multi-stage convex relaxation of T. Zhang (2010), the difference of convex functions of Shen, Pan and Zhu (2012, glmtlp's "tlp") and the reweighted ℓ1 of Candès, Wakin and Boyd (2008). CappedL1 and LogPenalty have no exact path and use it when lla is nothing.
• With start = :lasso the stages at λ start from the lasso at λ, as in the papers; with start = :previous they start from the solution at the previous level of a path of options.num_steps levels from λmax, as ncpen does.

lla does not take bounds (lower, upper), and it applies to these four penalties only. With lla, converged says whether coordinate descent converged in every stage.

Bounds

lower and upper bound the coefficients, lower ≤ β ≤ upper, as glmnet's lower.limits and upper.limits do: each is a number or a vector with a bound per coefficient, and zero must be feasible, lower ≤ 0 ≤ upper. For example, lower = 0 gives the nonnegative lasso, and upper = 0 with lower = 0 for a coefficient keeps it at zero. The intercept is not bounded. The penalty becomes Constrained(g, IndBox(l, u)), which keeps coordinate descent separable, and λmax and screening use only the directions the bounds allow, so that λmax of the nonnegative lasso is maxⱼ (Xⱼ'y / n)₊ / ωⱼ. With standardize_x and a penalty other than the ℓ1 norm, which applies to s .* β, the bounds of the stored penalty are s .* lower and s .* upper. Bounds work with every penalty above, in lasso, LassoPath, adaptive_lasso (through its keywords), CVLassoPath and glm_lasso.

refit_lasso_path(path::LassoPath, X, y)

Least squares refits on the supports along a LassoPath. Returns a Dict that maps each distinct support S (a vector of column indices) to X[:, S] \ y, or, when the path has an intercept, to [1 X[:, S]] \ y, whose first entry is the intercept.

Consecutive supports along a path differ by a few columns, so the Cholesky factor of X[:, S]'X[:, S] is updated from one support to the next (ProximalSolvers' RefitCache) instead of computed again. A support with at least n columns, or one whose columns are nearly collinear, is refit by \, which gives the minimum-norm solution.

scaled_lasso!(x, X, y, λ, ω, options=IterLassoOptions(); intercept=false, weights=nothing, offset=nothing,
             rng=Random.default_rng(), kwargs...)
scaled_lasso(X, y, λ, ω, options=IterLassoOptions(); intercept=false, weights=nothing, offset=nothing,
            rng=Random.default_rng(), kwargs...)
scaled_lasso!(x, f::CDSqrtQuadraticLoss, λ, ω, options=CDOptions(); rng=Random.default_rng(), kwargs...)

Scaled lasso of Sun and Zhang (2012), which estimates the coefficients and the noise level jointly:

minimize  ‖y - β₀ - Xβ‖² / (2nσ) + σ/2 + λ ∑ⱼ ωⱼ |βⱼ|   over β and σ > 0.

The intercept β₀ is not penalized, and it is zero unless intercept is true. The minimum over σ is at σ̂ = ‖y - β₀ - Xβ‖ / √n, so β̂ solves the square-root lasso of Belloni, Chernozhukov and Wang (2011),

minimize  ‖y - β₀ - Xβ‖ / √n + λ ∑ⱼ ωⱼ |βⱼ|,

which one run of coordinate descent solves (on CDSqrtLassoLoss, with the penalty level √n λ), with the options options.options_cd; the other fields of IterLassoOptions are not used. Coordinate descent starts from and overwrites x (unless options.options_cd.warm_start is false). Returns a LassoSolution whose σ is σ̂, and whose penalty is that of the lasso that β̂ solves, NormL1(λσ̂, ω).

The third form solves the same problem from sufficient statistics: with the loss CDSqrtQuadraticLoss(A, b, c), A = X'X/n, b = -X'y/n and c = y'y/n, it needs neither X nor y, and a coordinate update costs O(p) instead of O(n). With A = S, b = -S[:, j], c = S[j, j] and ω[j] = Inf it is the scaled lasso of variable j on the others from their covariance matrix S, as in the precision matrix estimator of Sun and Zhang (2013); for all the variables, check S once and build the losses with check = false. Its solution has no residuals (an empty vector) and no intercept.

The keyword arguments kwargs are fields of the options that replace those of options (of IterLassoOptions in the first two forms, of CDOptions in the third), and the random order of coordinate descent is drawn from rng. With observation weights w and an offset o, the loss is ∑ᵢ wᵢ rᵢ² / (2σ ∑ᵢ wᵢ) + σ/2 for the residuals r = y - o - β₀ - Xβ, so that σ̂ = sqrt(∑ᵢ wᵢ rᵢ² / ∑ᵢ wᵢ).

Before version 0.6.1, scaled_lasso! alternated lasso fits with updates of σ̂ until the relative change of σ̂ was below options.opt_tol (by default 1e-2), which stopped short of the solution and solved a lasso at every iteration.

scaled_lasso(X, y, λ, ω, options=IterLassoOptions())

Same as scaled_lasso!, starting from a zero iterate.

sqrt_lasso(X, y, λ, options=CDOptions(); intercept=false, standardize_x=true, gram=:auto, weights=nothing,
          offset=nothing, rng=Random.default_rng(), kwargs...)
sqrt_lasso(X, y, λ, ω, options=CDOptions(); intercept=false, gram=:auto, weights=nothing, offset=nothing,
          rng=Random.default_rng(), kwargs...)

Solve the square-root lasso of Belloni, Chernozhukov and Wang (2011)

minimize  ‖y - β₀ - Xβ‖ + λ ∑ⱼ ωⱼ |βⱼ|

by coordinate descent, with an unpenalized intercept β₀ that is zero unless intercept is true. In the first form the loadings are the standard deviations ωⱼ = sqrt(∑ᵢ (Xᵢⱼ - x̄ⱼ)² / n) of the columns when standardize_x is true (centered only with an intercept) and ωⱼ = 1 otherwise. The optimal λ does not depend on the noise level. Returns a LassoSolution whose σ is the sample standard deviation of the residuals. The keyword arguments kwargs, gram and rng are as for lasso; on the Gram matrix the loss is the square root of a quadratic (as in CDSqrtQuadraticLoss). With observation weights w and an offset o, the loss is the weighted norm sqrt(∑ᵢ wᵢ (yᵢ - oᵢ - β₀ - Xᵢ'β)²), the loadings are weighted as in lasso, and σ is the weighted root mean square of the residuals.

adaptive_lasso(X, y, λ, loss::ProximalBase.ScalarLoss, options=CDOptions(); init=:lasso, initλ=λ, γ=1,
              intercept=true, standardize_x=true, kwargs...)

The adaptive lasso of a generalized linear model: as above, with the loss of glm_lasso (ProximalBase.Logistic(), ProximalBase.Poisson(), ...) in place of least squares. The first estimate β̃ is the glm_lasso fit at initλ (init = :lasso), the ridge-penalized fit ∑ᵢ wᵢ ℓ(β₀ + xᵢ'β + oᵢ, yᵢ) / ∑ᵢ wᵢ + initλ ‖β‖² / 2 on the standardized columns with standardize_x (init = :ridge, glmnet's alpha = 0), or a given vector; the second fit is glm_lasso with the loadings ωⱼ = 1/|β̃ⱼ|^γ and no standardization. The keywords kwargs (weights, offset, lower, upper, newton_tol, newton_max_iter) apply to both fits. Returns the GLMLassoSolution of the second fit.

adaptive_lasso(X, y, λ, options=CDOptions(); init=:lasso, initλ=λ, γ=1, intercept=false, standardize_x=false,
              kwargs...)

The adaptive lasso of Zou (2006),

minimize  ‖y - β₀ - Xβ‖² / (2n) + λ ∑ⱼ |βⱼ| / |β̃ⱼ|^γ,

whose loadings come from a first estimate β̃: large coefficients are penalized less, so that they are less biased, and small ones more, so that fewer variables that do not matter are selected. β̃ is

• init = :lasso: the lasso at initλ, by default λ, with intercept, standardize_x and the keywords kwargs;
• init = :ridge: ridge regression at initλ, ‖y - β₀ - Xβ‖² / (2n) + initλ ‖β‖² / 2, on the standardized columns with standardize_x (glmnet's alpha = 0);
• init = β̃, a vector: for example least squares, X \ y when n > p, or the coefficients of a lasso whose level was chosen by cross-validation.

A coefficient with β̃ⱼ = 0 has an infinite loading and stays at zero. The second fit has no standardization, as with γ = 1 the loadings already make it invariant to the scales of the columns. The keywords kwargs (weights, offset, gram, rng and the fields of CDOptions) apply to both fits. Returns the LassoSolution of the second fit, whose penalty is NormL1(λ, ω) with the loadings ωⱼ = 1/|β̃ⱼ|^γ.

Reference: H. Zou (2006). The adaptive lasso and its oracle properties. Journal of the American Statistical Association 101(476), 1418–1429.

GLMLassoPath(X, y, loss = ProximalBase.Logistic(), options = CDOptions(); nλ = 100, λminratio, kwargs...)
GLMLassoPath(X, y, λpath, loss = ProximalBase.Logistic(), options = CDOptions();
             intercept = true, standardize_x = true, weights = nothing, offset = nothing,
             alpha = 1, penalty = nothing, lower = -Inf, upper = Inf, screening = true, max_hat_s = Inf,
             devmax = 0.999, newton_tol = 1e-8, newton_max_iter = 100, rng = Random.default_rng())

Solutions β0path[k] and βpath[k] of the penalized generalized linear model of glm_lasso along a decreasing sequence λpath, each one warm started from the previous solution and its intercept, by proximal Newton with the inner coordinate descent of options (ProximalSolvers' solve_path with ProxNewton).

Without λpath, the path has nλ values spaced evenly on the log scale from λmax, the smallest value whose solution is zero (at the intercept of the fit with β = 0), down to λminratio ⋅ λmax; λminratio defaults to 0.01 when X has more columns than rows and to 1e-4 otherwise. The first value is λmax times 1 + √eps, so that rounding cannot make the first solution nonzero.

With screening, each λ is solved first over the coordinates kept by the sequential strong rule of Tibshirani et al. (2012), and the optimality conditions of the others are checked afterwards (kkt_checks[k] times at λpath[k]), as in LassoPath. The path stops early after the first solution with more than max_hat_s nonzero coefficients, or once the fit explains a fraction devmax of the null deviance, as glmnet's does: along a logistic path on separable data the coefficients grow without bound, and the explained fraction devratio[k] = 1 - deviance / nulldeviance tends to one. devmax = 1 turns this off. The returned λpath then contains only the values that were solved for. iterations[k] and converged[k] are the passes of coordinate descent at λpath[k], summed over the proximal Newton steps, and whether the solve met its stopping rules. nobs is the number of observations.

alpha, penalty, lower and upper are as for glm_lasso; along a path of a nonconvex penalty, each solution is the local one reached from the previous solution. penalty holds the penalty at level one (NormL1(1, ω) for the lasso), applied to scales .* β when scales is not nothing. coef, predict and nobs apply to the path, as for LassoPath, and fit(GLMLassoPath, X, y, ...) is GLMLassoPath(X, y, ...).

GLMLassoSolution

Result returned by glm_lasso.

Fields

• x::SparseIterate: the estimated coefficients β
• intercept: the unpenalized intercept β₀, zero unless the solver was called with intercept = true
• penalty: the penalty at which x was computed: NormL1(λ, ω), the penalty of alpha or penalty at the level λ, or either on the box of lower and upper, Constrained(g, IndBox(l, u))
• loss: the scalar loss of the model, such as ProximalBase.Logistic()
• linearpredictor::Vector: η = β₀ + Xβ + offset
• deviance: 2 ∑ᵢ wᵢ (ℓ(ηᵢ, yᵢ) - minᵤ ℓ(u, yᵢ)), twice the negative log-likelihood ratio against the saturated model for the logistic and Poisson losses, with the observation weights wᵢ
• nulldeviance: the deviance of the fit with β = 0 (the intercept and the offset only); 1 - deviance / nulldeviance is the fraction of the deviance that the fit explains
• iterations::Int: passes of coordinate descent, summed over the proximal Newton steps
• converged::Bool: whether proximal Newton, and coordinate descent at each of its steps, met their stopping rules
• hasintercept::Bool: whether the solver fitted an intercept
• scales: nothing when penalty applies to x, or the scales s of the columns when it applies to s .* x (standardize_x with a penalty other than the ℓ1 norm), as for LassoSolution

The functions of StatsAPI apply to it: coef, predict (on the scale of the response or of the linear predictor), fitted, nobs, dof, deviance and nulldeviance.

glm_lasso(X, y, λ, loss = ProximalBase.Logistic(), options = CDOptions();
         intercept = true, standardize_x = true, weights = nothing, offset = nothing,
         alpha = 1, penalty = nothing, lower = -Inf, upper = Inf, newton_tol = 1e-8, newton_max_iter = 100,
         rng = Random.default_rng())

Solve the penalized generalized linear model

minimize  ∑ᵢ wᵢ ℓ(β₀ + xᵢ'β + oᵢ, yᵢ) / ∑ᵢ wᵢ + λ ∑ⱼ ωⱼ |βⱼ|

for a scalar loss ℓ of ProximalBase: ProximalBase.Logistic(), the negative log-likelihood of logistic regression (yᵢ ∈ [0, 1], usually {0, 1}); ProximalBase.Poisson(), Poisson regression with the log link (yᵢ ≥ 0); ProximalBase.Huber(δ); or a loss of your own. The observation weights wᵢ (nonnegative, one when omitted) and the offset oᵢ (zero when omitted) are those of weights and offset. The intercept β₀ is not penalized, and it is zero unless intercept is true. With standardize_x the loadings ωⱼ are the standard deviations of the columns, weighted by w and centered at their weighted means with an intercept, so that the coefficients are on the scale of the data; otherwise ωⱼ = 1. These are glmnet's conventions and defaults. lower and upper bound the coefficients, lower ≤ β ≤ upper, as for lasso.

alpha and penalty give the elastic net and the other penalties, as for lasso: alpha < 1 is glmnet's elastic net λ (α ∑ⱼ |βⱼ| + (1 - α) ‖β‖² / 2), and penalty a penalty of ProximalBase whose level λ replaces (ElasticNet(1, λ2), MCP(1, a), SCAD(1, a), NormL1Nonneg(1), or NormL1(1, w) for weights). With standardize_x, such a penalty applies to the coefficients of the standardized columns, as in glmnet; the coefficients returned are those of the columns of X, and the scales are stored in the solution's scales. The elastic net equals glmnet's binomial and poisson fits. A nonconvex penalty (MCP, SCAD) is computed along options.num_steps levels from λmax whatever options.warm_start, as ncvreg does; num_steps = 0 solves at λ from zero. The solutions are stationary points of the objective above. ncvreg's binomial and poisson families do not solve it: their coordinate steps scale the concavity of the penalty by the weighted curvature of each column at the current fit (they divide by v (1 - 1/a) where the exact step divides by v - 1/a), so their solutions differ (in its gaussian family v = 1 and the two agree).

The problem is solved by proximal Newton (ProximalSolvers' ProxNewton, the method of glmnet): each step minimizes the weighted least-squares model of the loss at the current fit by coordinate descent with options, then backtracks, until the objective changes by less than newton_tol relative to its value or after newton_max_iter steps. With options.warm_start false (or a nonconvex penalty), the solution is computed along options.num_steps levels from λmax (see GLMLassoPath), each warm started from the previous one; otherwise it is computed directly, starting from β = 0 and the intercept of the fit with β = 0. The random order of coordinate descent is drawn from rng. Returns a GLMLassoSolution.

y = Float64.(rand(n) .< 1 ./ (1 .+ exp.(-X * β)))
sol = glm_lasso(X, y, 0.05)                               # logistic
sol = glm_lasso(X, counts, 0.05, ProximalBase.Poisson(); offset = log.(exposure))

QuantileLassoPath(X, y, λpath; tau=0.5, method=:exact, max_hat_s=Inf, kwargs...)
QuantileLassoPath(X, y; nλ=100, λminratio=0.01, kwargs...)

Fit mean(rho_tau(y - intercept - X*β - offset)) + λ*sum(w .* abs.(β)) along a nonincreasing, finite, nonnegative grid. Every saved level retains both coefficients and raw dual scores, including at zero residuals. path[k] is a QuantileLassoSolution; path.λpath contains the returned levels and requested_λpath the requested grid. complete is false when a limit stops the path. There is no coefficient or dual interpolation: coef(path, λ) requires a saved level.

max_hat_s=Inf limits the number of nonzero slopes. It must be a nonnegative real number other than NaN. After saving the first fit with count(!iszero, β) > max_hat_s, later levels are skipped. The intercept does not count; nonzero slopes with zero penalty factors do. The exceeding fit is retained, as in LassoPath and GLMLassoPath; QR.jl dropped it. If levels are skipped, complete=false and status=:max_hat_s. Exceeding the limit only at the last requested level leaves complete=true and the usual status.

method=:exact uses HiGHS dual simplex with an analytic quantile starting basis and reuses that basis/model along the path. Load HiGHS to activate the optional extension. initial_basis=:default disables the analytic start; simplex=:primal selects primal simplex. time_limit limits the entire exact path between levels and passes the remaining time to HiGHS. Preparation and diagnostics can overrun that limit, and automatic-grid preparation precedes the path timer. A stopped path retains available primal–dual results and the termination status; it never fills missing levels. The support limit is checked only after an OPTIMAL level with a saved pair. Non-OPTIMAL termination or a missing pair takes precedence. The path time budget is checked before each level; a support stop after an optimal level takes precedence over a budget exhausted during that level.

method=:smooth uses the existing proximal Newton path with the Epanechnikov convolution loss at explicit bandwidth (default 0.2). Its raw dual scores are the smoothed loss scores, not exact check-loss duals. time_limit is currently supported only by the exact backend. options, newton_tol, newton_max_iter, screening and rng control proximal Newton.

The intercept is unpenalized (intercept=true). penalty_factors are finite nonnegative slope weights, one by default. With standardize_x=true (the default), they are multiplied by column population standard deviations, centered when fitting an intercept. Coefficients and duals remain on the original data scale. Zero factors and λ=0 are supported on explicit grids; the automatic grid requires positive effective factors. offset is fixed and unpenalized. Observation weights and formula fits are not implemented here. Inputs are converted to Float64; sparse designs remain sparse.

The automatic grid starts at the intercept-only score maximum, inflated for roundoff, and ends at λminratio times that maximum. This is a numerical path, not pivotal calibration or cross-validation; use quantile_lambda or CVQuantileLassoPath for those choices. tol=2e-5 is the external KKT and relative raw-gap tolerance. Backends supporting convolution certificates also use it during smoothed fitting, honoring a smaller positive options.gaptol. An explicit options.gaptol=0 keeps the legacy objective-change rule; older backends retain their original stopping and are checked afterward.

For positive λ and factors, the certificate projects scores onto their box and intercept constraint, then contracts to meet the slope constraints. With zero λ or any zero factor it uses the feasible zero vector, whose bound may be loose; numerical convergence still checks the raw pair. No inferential validity is asserted merely from solver convergence.

QuantileLassoSolution

A quantile-lasso fit on the original columns of X. β, intercept, residuals, λ and the effective penalty_factors describe the fitted problem. dual contains the raw observation-level scores; see dual_scores. certificate_dual is a separate feasible vector used only for lower bounds. No rounding or projection is applied to dual or β.

objective is the fitted (exact or smoothed) penalized objective; check_objective always uses the check loss. duality_gap uses the raw dual, whereas certificate_gap uses the feasible certificate. kkt_residual, dual_feasibility and complementarity check the returned primal–dual pair. converged uses those numerical checks at the requested tolerance; native_converged separately records the backend's stopping flag. status retains its termination reason. Duals can be nonunique at zero residuals.

QuantileRefit

An unpenalized quantile fit on the sorted column indices in support. β has one entry per original column and is zero off that support; intercept is the fitted intercept. fit is the restricted QuantileLassoSolution at λ=0, including its raw dual scores, objectives, numerical diagnostics and native termination status.

coef puts the intercept first when fitted. predict takes the full design and an optional offset. residuals, nobs and dual_scores use fit. dof counts length(support) + fit.hasintercept, including selected columns whose refitted coefficient is zero. This result provides computation, not post-selection inference.

dual_scores(sol::QuantileLassoSolution)
dual_scores(path::QuantileLassoPath)
dual_scores(r::QuantileRefit)

Raw observation scores u. At a valid solution, tau-1 ≤ uᵢ ≤ tau. A stopped fit can retain an infeasible iterate; check converged and its numerical diagnostics before using it. At an exact solution, uᵢ=tau for a positive residual and tau-1 for a negative residual; scores at zero residuals come from the optimizer, not from a residual-sign rule. With an intercept, sum(u)=0; slope stationarity is X'u/n ∈ λ*w .* ∂abs(β). For a smooth fit these are convolution-loss scores. A whole path returns an n-by-level matrix; dual_scores(path[k]) returns one vector. The returned vectors are owned by the saved fits, not by a live optimizer.

QR.jl's historical getXi convention is -dual_scores(sol)/nobs(sol), the dual of y - a - X*β == rp - rn with the mean check-loss objective. Apply the same conversion to a path; transpose it for QR.jl's level-by-observation layout. This changes only sign and normalization, and always uses the raw scores.

quantile_lasso(X, y, λ; kwargs...)

Fit one level of QuantileLassoPath, returning a QuantileLassoSolution. A backend limit with no available primal–dual pair throws an error; an available but inaccurate pair has converged=false.

quantile_refit(X, y, support; tau=0.5, method=:exact, intercept=true,
               offset=nothing, bandwidth=0.2, tol=2e-5, time_limit=Inf,
               initial_basis=:quantile, simplex=:dual, options=CDOptions(...),
               newton_tol=1e-10, newton_max_iter=100, screening=true,
               rng=Random.default_rng())
quantile_refit(X, y, sol::QuantileLassoSolution; tau=sol.tau, offset=nothing, kwargs...)

Refit without a penalty on a vector of integer column indices, possibly empty. The support is copied, sorted and checked for duplicates and indices outside 1:size(X,2). Use findall to convert a Bool mask. A scalar tau returns a QuantileRefit; a nonempty vector returns one refit per quantile in its supplied order, including repeats. Quantiles must lie in (0,1). Inputs are converted to Float64; sparse designs remain sparse.

The solution form takes the exact nonzero support of sol.β and inherits its method, bandwidth and intercept setting. Dimensions are checked, then y - offset - sol.intercept - X*sol.β must reproduce sol.residuals within sqrt(eps(Float64))*(1 + maximum(abs, y - offset)) in maximum norm. This detects incompatible responses, offsets or scaling but does not establish that X is the original design: a change with ΔX*sol.β == 0 passes.

Columns outside the support are fixed at zero, including zero-loading columns with zero fitted coefficients. List them explicitly to keep them. No support threshold, penalty factors, standardization or support limit is accepted. tol, options, Newton and screening controls have the same meaning as for QuantileLassoPath.

Exact refits require using HiGHS. One model is reused across quantiles by changing only score bounds. Every requested refit is attempted, including after a non-optimal native termination with an available primal–dual pair. Inspect r.fit.status, native_converged and converged separately. A missing pair throws an error. time_limit is one deadline for the entire call; an expired deadline before a requested solve throws an error. Preparation and diagnostics can overrun it. Finite limits are rejected for method=:smooth, which fits each quantile independently without warm starts.

Diagnostics use only the restricted design. At λ=0 the feasible certificate is zero and its gap can be loose; convergence checks the raw pair. Scores and coefficients can be nonunique in degenerate problems; the optimizer's values are retained. These refits do not provide post-selection inference.

refit_quantile_path(path::QuantileLassoPath, X, y; tau=nothing,
                    offset=nothing, kwargs...)

Refit each distinct exact-nonzero support of the saved fits once. Return a Dict{Vector{Int},QuantileRefit}, indexed by refits[findall(!iszero, path[k].β)]. An empty path returns an empty dictionary; a truncated path uses only its saved fits.

The method and intercept setting come from path. Each support uses the quantile and bandwidth of its first saved fit, unless a scalar tau is supplied. Each such source fit passes the dimension and residual consistency checks of quantile_refit, with the same limits on what they establish. Solver controls and the interruption contract are those of quantile_refit. For exact refits, time_limit is one deadline across all supports. Models are not reused across different supports. No inference is provided.

CVQuantileLassoPath(X, y[, λpath]; nfolds=10, folds=nothing,
                   rng=Random.default_rng(), threads=true, tau=0.5,
                   nλ=100, λminratio=0.01, offset=nothing, kwargs...)

Select a quantile-lasso penalty by held-out check loss, for either exact or convolution-smoothed fits. kwargs go to QuantileLassoPath, including method, bandwidth, loadings, standardization and solver controls. max_hat_s is forwarded to the full-data path and every fold path. A support stop on the full data truncates the shared grid. A fold that stops earlier has NaN errors at its unreached levels; these are ineligible, with cvm=Inf and cvse=NaN, even when no fold reaches them. An explicit decreasing λpath is shared by all folds; otherwise the full-data path supplies the common grid. Scaling and centering are computed using each training fold. The automatic common grid itself uses the full data.

cvm and cvse use fold-size weighting and the same standard-error convention as CVLassoPath. λmin minimizes cvm; λ1se is the largest eligible penalty within one standard error of that minimum. A level is eligible only when the full-data fit and every fold's fit pass numerical convergence. Ineligible levels have cvm=Inf; fold_errors is a folds-by-level matrix with NaN for failed/missing fold fits. If no level is eligible, an error is raised. time_limit, when supplied, applies separately to each exact path.

The result retains the full-data path, λpath, cvm, cvse, nnz, λmin, λ1se, folds and fold_errors. coef(cv), predict(cv, X) and dual_scores(cv) select the saved full-data fit at λ1se, or at λmin with λ=:λmin. Duals are never averaged across folds or interpolated. Use dual_scores(cv.path) to access the whole full-data path. Offsets follow observations into training and held-out evaluation; pass new prediction offsets explicitly. Observation weights are not supported.

quantile_lambda(X; tau=0.5, c=2.0, gamma=0.1, nsim=1000,
                penalty_factors=nothing, standardize_x=true,
                intercept=true, rng=Random.default_rng())

Simulate a fixed-quantile pivotal penalty in the normalization of quantile_lasso: mean check loss plus λ*sum(w .* abs.(β)). For independent uniforms Uᵢ, simulate max_j abs(sum_i Zᵢⱼ*(tau - (Uᵢ < tau)))/(n*wⱼ) and return c times its empirical (1-gamma) quantile (linear interpolation, R's type 7). Z is X centered by column when fitting an unpenalized intercept, and X otherwise. w contains the effective loadings used by the fit, including its population-standard-deviation scaling when standardize_x=true.

This is the single-quantile score simulation of Belloni and Chernozhukov (2011), translated to this objective's scale. It depends on X, not y, estimates no noise scale, and does not calibrate a simultaneous quantile process. The paper's defaults here are c=2 and gamma=0.1; nsim controls Monte Carlo precision. Use a seeded rng for reproducibility and the same loading/intercept settings in the subsequent fit. All effective slope loadings must be positive; unpenalized slopes require a separately justified calibration or an explicit λ. The rule calibrates exact check-loss scores; it is not a bandwidth calibration or a finite-sample guarantee for a convolution-smoothed fit.

GroupLassoPath(X, y, groups, options = CDOptions(); nλ = 100, λminratio, kwargs...)
GroupLassoPath(X, y, λpath, groups, options = CDOptions(); intercept = true, standardize = :orthonormal,
               groupweights = sqrt.(length.(groups)), screening = true, max_groups = Inf)

Solutions β0path[k] and βpath[k] of the group lasso of group_lasso along a decreasing sequence λpath, each one warm started from the previous solution, by ProximalSolvers' solve_path.

Without λpath, the path has nλ values spaced evenly on the log scale from λmax, the smallest value whose solution is zero, down to λminratio ⋅ λmax; λminratio defaults to 0.01 when X has more columns than rows and to 1e-4 otherwise. The first value is λmax times 1 + √eps.

With screening, each λ is solved first over the groups kept by the sequential strong rule, and the optimality conditions of the others are checked afterwards (kkt_checks[k] times at λpath[k]). The path stops early after the first solution with more than max_groups nonzero groups; the returned λpath then contains only the values that were solved for. activegroups[k] lists the nonzero groups at λpath[k], and iterations[k] and converged[k] are the passes of coordinate descent there and whether it met its stopping rule. coef, predict and nobs apply to the path, as for LassoPath, and fit(GroupLassoPath, X, y, ...) is GroupLassoPath(X, y, ...).

GroupLassoSolution

Result returned by group_lasso.

Fields

• x::SparseVector: the coefficients β, on the scale of the data
• intercept: the unpenalized intercept β₀, zero unless intercept = true
• λ: the penalty level
• groups::Vector{Vector{Int}}, groupweights::Vector: the groups and their weights wₖ
• activegroups::Vector{Int}: the groups with nonzero coefficients
• iterations::Int, converged::Bool: the passes of coordinate descent and whether it met its stopping rule
• hasintercept::Bool, nobs::Int: whether the solver fitted an intercept, and the number of observations

coef, predict, nobs and dof apply to it.

group_lasso(X, y, λ, groups, options = CDOptions(); intercept = true, standardize = :orthonormal,
           groupweights = sqrt.(length.(groups)))

The group lasso

minimize  ‖y - β₀ - Xβ‖² / (2n) + λ ∑ₖ wₖ ‖X̃_Gₖ βGₖ‖₂ / √n

for disjoint groups Gₖ that cover the columns of X (give a column its own group to treat it alone, and a zero weight to leave a group unpenalized), with an unpenalized intercept β₀ that is zero unless intercept is true, and X̃_G the columns of G, centered with an intercept. The weights default to wₖ = √|Gₖ|.

With standardize = :orthonormal, the default, the penalty is on the fitted values of each group, as above, which makes the estimate invariant to a change of basis within a group (Simon and Tibshirani, 2012; R's grpreg). It is fitted on orthonormalized groups: each X̃_G is factored by a pivoted QR as Z_G R̃, with Z_G'Z_G / n = I, the group lasso is solved in the coefficients θ_G = R̃ β_G by block coordinate descent, whose steps are then exact, and mapped back. The columns of a group that are linear combinations of its other columns get zero coefficients. With standardize = :none, the penalty is λ ∑ₖ wₖ ‖βGₖ‖₂ on the columns as given.

Coordinate descent runs with options; with options.warm_start false the solution is computed along options.num_steps levels from λmax. Returns a GroupLassoSolution. See also GroupLassoPath.

SparseGroupLassoPath(X, y, groups, options = CDOptions(); nλ = 100, λminratio, kwargs...)
SparseGroupLassoPath(X, y, λpath, groups, options = CDOptions(); alpha = 0.95, intercept = true,
                     standardize_x = true, groupweights = sqrt.(length.(groups)), screening = true,
                     max_groups = Inf, rng = Random.default_rng())

Solutions β0path[k] and βpath[k] of the sparse group lasso of sparse_group_lasso along a decreasing sequence λpath, each one warm started from the previous solution, by ProximalSolvers' solve_path, at the fixed mix alpha.

Without λpath, the path has nλ values spaced evenly on the log scale from λmax, the smallest value whose solution is zero, down to λminratio ⋅ λmax; λminratio defaults to 0.01 when X has more columns than rows and to 1e-4 otherwise. The first value is λmax times 1 + √eps. A group at zero is optimal while ‖soft(∇_G, λα)‖₂ ≤ λ(1 - α) w_G (soft the soft-thresholding of its gradient); λmax is the largest root of these conditions at zero, and the sequential strong rule (with screening) and its checks compare each group's root with λ.

The path stops early after the first solution with more than max_groups nonzero groups; the returned λpath then contains only the values that were solved for. activegroups[k] lists the nonzero groups at λpath[k]. coef, predict and nobs apply to the path, as for LassoPath, and fit(SparseGroupLassoPath, X, y, ...) is SparseGroupLassoPath(X, y, ...).

SparseGroupLassoSolution

Result returned by sparse_group_lasso.

Fields

• x::SparseVector: the coefficients β, on the scale of the data
• intercept: the unpenalized intercept β₀, zero unless intercept = true
• λ, alpha: the penalty level and the mix of the ℓ1 and group terms
• groups::Vector{Vector{Int}}, groupweights::Vector: the groups and their weights wₖ
• activegroups::Vector{Int}: the groups with nonzero coefficients
• scales: nothing, or the scales s of the columns when the penalty applies to s .* x (standardize_x)
• iterations::Int, converged::Bool: the passes of coordinate descent and whether it met its stopping rule
• hasintercept::Bool, nobs::Int: whether the solver fitted an intercept, and the number of observations

coef, predict, nobs and dof apply to it.

sparse_group_lasso(X, y, λ, groups, options = CDOptions(); alpha = 0.95, intercept = true, standardize_x = true,
                 groupweights = sqrt.(length.(groups)), rng = Random.default_rng())

The sparse group lasso of Simon, Friedman, Hastie and Tibshirani (2013),

minimize  ‖y - β₀ - Xβ‖² / (2n) + λ (α ∑ⱼ |β̃ⱼ| + (1 - α) ∑ₖ wₖ ‖β̃_Gₖ‖₂),

which selects groups, as the group lasso does, and also variables within the selected groups, for disjoint groups Gₖ that cover the columns of X. With standardize_x, β̃ = s .* β are the coefficients of the columns divided by their standard deviations s (centered with an intercept); otherwise β̃ = β. R's sparsegl and SGL scale the columns to unit Euclidean norm instead, √n times the standard deviation, so their level λ is √n λ here. The coefficients returned are those of the columns of X. alpha mixes the two terms: alpha = 1 is the lasso and alpha = 0 the group lasso on the columns as given (group_lasso(...; standardize = :none)); the default 0.95 is SGL's. The weights default to wₖ = √|Gₖ|, and a zero weight with alpha = 0 leaves a group unpenalized. The intercept β₀ is not penalized, and it is zero unless intercept is true.

It is solved by block coordinate descent on the groups (ProximalSolvers' CD on an AtomIterate, whose block steps are the exact proximal operator of the penalty plus a quadratic bound), with options; with options.warm_start false the solution is computed along options.num_steps levels from λmax. The random order is drawn from rng. Unlike group_lasso, the groups are not orthonormalized: the ℓ1 term is not invariant to a change of basis within a group. Returns a SparseGroupLassoSolution. See also SparseGroupLassoPath.

Reference: N. Simon, J. Friedman, T. Hastie and R. Tibshirani (2013). A sparse-group lasso. Journal of Computational and Graphical Statistics 22(2), 231–245.

CVMultiLassoPath(X, Y, options = CDOptions(); nfolds = 10, folds = nothing,
                 rng = Random.default_rng(), threads = true, nλ = 100, λminratio, kwargs...)

Cross-validation of MultiLassoPath, which minimizes

‖Y - 1β₀' - XB‖²_F / (2n) + λ ∑ⱼ hⱼ((sⱼ B[j, t])ₜ).

The row penalty is selected by row_norm and alpha as in multi_lasso. Observations are assigned to nfolds balanced folds, shared by all responses, from rng, or to the labels in folds. Each training fold refits the centering, scales and coefficients at the full-data path's levels. The held-out error of fold k is the mean squared error over its observations and responses.

Fields

• path::MultiLassoPath, λpath::Vector: the full-data path and its levels
• cvm::Vector, cvse::Vector: the fold-size-weighted error and its standard error, as in CVLassoPath
• cvm_response::Matrix: the held-out mean squared error of each response (a row) at each level (a column)
• nrows::Vector{Int}: the number of selected rows of the full-data path
• λmin, λ1se: the minimum-error and one-standard-error levels, among those reached and converged in every fold
• folds::Vector{Int}: the fold labels of the observations

Unreached or nonconverged fold levels have no error (NaN); they are never selected. Each fold's generator is seeded from rng, independently of threads. coef(cv) and predict(cv, X) use λ1se, or λmin with λ = :λmin.

MultiLassoPath(X, Y, options = CDOptions(); nλ = 100, λminratio, kwargs...)
MultiLassoPath(X, Y, λpath, options = CDOptions(); intercept = true, standardize_x = true,
               screening = true, max_rows = Inf, row_norm = :l2, alpha = 0,
               rng = Random.default_rng())

Solutions β0path[k] and Bpath[k] of the multi-response lasso of multi_lasso along a decreasing sequence λpath, each one warm started from the previous solution, by ProximalSolvers' solve_path for the default penalty and by row coordinate descent for the row ℓ1/ℓ∞ and sparse-group penalties. row_norm and alpha record the penalty on the path; the old positional constructor fills :l2 and zero.

Without λpath, the path has nλ values spaced evenly on the log scale from λmax, the smallest value whose solution is zero, down to λminratio ⋅ λmax; λminratio defaults to 0.01 when X has more columns than rows and to 1e-4 otherwise. The first value is λmax times 1 + √eps.

With screening, each λ is solved first over the rows kept by the sequential strong rule, and the optimality conditions of the others are checked afterwards (kkt_checks[k] times at λpath[k]). The path stops early after the first solution with more than max_rows nonzero rows; the returned λpath then contains only the values that were solved for. activerows[k] lists the selected predictors at λpath[k]. coef and predict give a vector with a matrix per level (coef(path, λ) one matrix, interpolated in λ), nobs the number of observations, and fit(MultiLassoPath, X, Y, ...) is MultiLassoPath(X, Y, ...).

MultiLassoSolution

Result returned by multi_lasso.

Fields

• B::SparseMatrixCSC: the p × q matrix of coefficients, a column per response
• intercept::Vector: the unpenalized intercepts of the responses, zero unless intercept = true
• λ: the penalty level
• activerows::Vector{Int}: the rows with a nonzero coefficient; entries of a selected row can still be zero
• iterations::Int, converged::Bool: the passes of coordinate descent over the rows and whether it met its stopping rule
• hasintercept::Bool, nobs::Int: whether the solver fitted intercepts, and the number of observations
• row_norm::Symbol, alpha: the row norm and share of the entrywise penalty

coef (a (p + 1) × q matrix, the intercepts in the first row), predict (an n × q matrix), nobs and dof apply to it.

multi_lasso(X, Y, λ, options = CDOptions(); intercept = true, standardize_x = true,
            row_norm = :l2, alpha = 0, rng = Random.default_rng())

The multi-response lasso of the n × q responses Y on X,

minimize  ‖Y - 1β₀' - XB‖²_F / (2n) + λ ∑ⱼ h((sⱼ B[j, t])ₜ)

where h(u) = alpha * ‖u‖₁ + (1 - alpha) * ‖u‖₂ with row_norm = :l2, or h(u) = ‖u‖∞ with row_norm = :linf (which requires alpha = 0). alpha must lie in [0, 1]. The default row ℓ1/ℓ2 penalty is glmnet's mgaussian family. Row selection allows zero entries within a selected row. The intercepts β₀ are not penalized, and they are zero unless intercept is true. The scales sⱼ are the population standard deviations of the columns, centered with an intercept, with standardize_x, and one otherwise. Coefficients are returned on the original scale. With one response all three penalties reduce to lasso with matching standardization.

It is solved by block coordinate descent on the rows of B, with options. The default uses ProximalSolvers' multi-response least-squares model; the new penalties use row coordinate descent on the shared design, without stacking. Its steps are exact, and orders are drawn from rng; with options.warm_start false the solution is computed along options.num_steps levels from λmax. Returns a MultiLassoSolution. See also MultiLassoPath and CVMultiLassoPath.

CVMultitaskLassoPath(Xs, ys, options = CDOptions(); nfolds = 10, folds = nothing,
                     rng = Random.default_rng(), threads = true, nλ = 100, λminratio, kwargs...)

K-fold cross-validation of the path of MultitaskLassoPath (with the same options and keywords, such as intercept, standardize_x, task_weights, penalty_factors, row_norm, alpha and max_rows), as CVLassoPath does for one task: the path on all the data gives the levels, each fold's path is computed on the other folds at those levels, and the error on the held-out fold estimates the prediction error at each level. The automatic grid needs finite λmax; sparse group with alpha > 0 permits zero factors.

The folds are drawn within each task: the observations of every task are assigned to the K = nfolds folds in a balanced random way, drawn from rng (K must not exceed the size of any task). folds gives them instead, a vector of labels 1, …, K per task, with every label used in every task. Each fold's path refits the whole procedure on the observations of the other folds, with the same task_weights argument: the centering, the scales, the weights of :observations and the design caches are computed on them, so that nothing leaks from the held-out fold. Every task must keep two training observations in each fold with an intercept, and one without.

The error of fold k at a level is

eₖ = ∑ₜ (cₜ / nₜ) ∑ᵢ (y_ti - ŷ_ti)² / wₖ,   wₖ = ∑ₜ cₜ mₜₖ / nₜ,

where the inner sum runs over the mₜₖ observations of task t in fold k, and cₜ and nₜ are those of all the data. cvm is the average of the fold errors weighted by the wₖ, and cvse its standard error, sqrt(∑ₖ wₖ (eₖ - cvm)² / ∑ₖ wₖ / (K - 1)); with one task they are those of CVLassoPath. cvm_task holds the held-out mean squared errors MSEₜ of the tasks, whose folds are weighted by the mₜₖ; at the levels that every fold reaches, cvm = ∑ₜ cₜ MSEₜ / ∑ₜ cₜ, the held-out counterpart of the loss. λmin minimizes cvm, and λ1se is the largest level whose cvm is within cvse of the minimum (the one-standard-error rule). A fold path that stops early (max_rows) leaves the remaining levels out of the averages; nonconverged fold levels are treated the same way and are never selected. Levels that no fold reaches have NaN in cvm, cvse and cvm_task. λmin and λ1se are chosen among the levels that every fold reaches and converges at.

With threads, the folds run in parallel on the available threads; each fold draws its random coordinate orders from a generator seeded from rng, so that the result does not depend on the number of threads.

Returns a CVMultitaskLassoPath; coef(cv) and predict(cv, Xs) use the solution at λ1se, or at λmin with λ = :λmin, and fit(CVMultitaskLassoPath, Xs, ys, ...) is CVMultitaskLassoPath(Xs, ys, ...).

CVMultitaskLassoPath

Result returned by CVMultitaskLassoPath(Xs, ys).

Fields

• path::MultitaskLassoPath: the path on all the data
• λpath::Vector: its levels, the grid of every fold
• cvm::Vector, cvse::Vector: the cross-validated prediction error at each level, and its standard error (NaN at the levels that no fold reaches, and cvse also where only one does)
• cvm_task::Matrix: the held-out mean squared error of each task (a row) at each level (a column), NaN where no fold reaches the level
• nrows::Vector{Int}: the number of nonzero rows of path at each level
• λmin, λ1se: among the levels reached and converged in every fold, the level with the smallest cvm, and the largest level whose cvm is within one standard error of it
• folds::Vector{Vector{Int}}: the fold of each observation, a vector per task
• task_weights::Vector: the task weights cₜ of the path on all the data

MultitaskLassoPath(Xs, ys, options = CDOptions(); nλ = 100, λminratio, kwargs...)
MultitaskLassoPath(Xs, ys, λpath, options = CDOptions(); intercept = true, standardize_x = true,
                   task_weights = :equal, penalty_factors = nothing, rng = Random.default_rng(),
                   screening = true, max_rows = Inf, row_norm = :l2, alpha = 0)

Solutions β0path[k] and Bpath[k] of the multitask lasso of multitask_lasso (with the same keywords) along a decreasing sequence λpath, each one warm started from the previous solution. The default uses ProximalSolvers' solve_path, with momentum on the stacked problem unless momentum = false or order = :cyclic, and the options unchanged on a shared design, as in multitask_lasso; the row ℓ1/ℓ∞ and sparse-group penalties use row coordinate descent on the tasks' own designs.

Without λpath, the path has nλ values spaced evenly on the log scale from λmax, the smallest value whose solution is zero, down to λminratio ⋅ λmax; λminratio defaults to 0.01 when there are fewer observations than coefficients, ∑ₜ nₜ < pq, and to 1e-4 otherwise, the rule of LassoPath for one task and of MultiLassoPath for a shared design. The first value is λmax times 1 + √eps. This grid needs positive penalty factors for row ℓ1/ℓ2 and ℓ1/ℓ∞. Sparse group with alpha > 0 also allows zero factors, since the entries remain penalized. The grid needs finite λmax; otherwise give λpath. For ℓ∞, λmax is the largest row gradient ℓ1 norm divided by its factor; for sparse group it is the largest root of ‖soft(gⱼ(0), alpha * λ)‖₂ = (1 - alpha) * ωⱼ * λ.

With screening, each λ is solved first over the rows kept by the sequential strong rule, and the optimality conditions of the others are checked afterwards (kkt_checks[k] times at λpath[k]). The path stops early after the first solution with more than max_rows nonzero rows; the returned λpath then contains only the values that were solved for.

Fields

• λpath::Vector: the levels solved for
• Bpath::Vector{SparseMatrixCSC}, β0path::Vector{Vector}: the p × q coefficients and the intercepts of the tasks at each level
• activerows::Vector{Vector{Int}}: the selected predictors at each level
• iterations::Vector{Int}, converged::Vector{Bool}, kkt_checks::Vector{Int}: the passes of coordinate descent at each level, whether it met its stopping rule, and the checks of the screening
• intercept::Bool: whether the path has intercepts
• nobs::Int, task_nobs::Vector{Int}, task_weights::Vector, penalty_factors::Vector, row_norm::Symbol, alpha: as in MultitaskLassoSolution

coef gives a vector with a coefficient matrix per level, with the intercepts in the first row when the path has them (coef(path, λ) one matrix, interpolated linearly in λ), predict(path, Xs) a vector with an nₜ × L matrix of fitted values per task, a column per level, nobs the number of observations, and fit(MultitaskLassoPath, Xs, ys, ...) is MultitaskLassoPath(Xs, ys, ...).

MultitaskLassoSolution

Result returned by multitask_lasso.

Fields

• B::SparseMatrixCSC: the p × q matrix of coefficients, a column per task, on the scale of the data
• intercept::Vector: the unpenalized intercepts of the tasks, zero unless intercept = true
• λ: the penalty level
• activerows::Vector{Int}: the selected predictors, the rows of B with a nonzero coefficient; entries of a selected row can still be zero
• iterations::Int, converged::Bool: the passes of coordinate descent over the rows and whether it met its stopping rule
• hasintercept::Bool: whether the solver fitted intercepts
• nobs::Int, task_nobs::Vector{Int}: the number of observations of all the tasks, N = ∑ₜ nₜ, and those of each task, nₜ
• task_weights::Vector, penalty_factors::Vector: the task weights cₜ and the penalty factors ωⱼ of the fit
• row_norm::Symbol, alpha: the row norm and share of the entrywise penalty

coef (a (p + 1) × q matrix, the intercepts in the first row), predict (a vector of fitted values per task), nobs and dof apply to it.

multitask_lasso(Xs, ys, λ, options = CDOptions(); intercept = true, standardize_x = true,
                task_weights = :equal, penalty_factors = nothing, row_norm = :l2, alpha = 0,
                rng = Random.default_rng())

The multitask lasso of the responses ys[t] on the designs Xs[t], t = 1, …, q, with the same p columns in every task (the same predictors, in the same order) and nₜ rows, which may differ between the tasks,

minimize  ∑ₜ cₜ ‖y_t - β₀ₜ - X_t b_t‖² / (2nₜ) + λ ∑ⱼ hⱼ((s₁ⱼ b₁ⱼ, …, s_qⱼ b_qⱼ))

over the intercepts β₀ₜ and the coefficients B = [b_1 … b_q]. The penalty on the rows of B selects predictors across tasks; a selected row can have zero entries. With row_norm = :l2, hⱼ(u) = alpha * ‖u‖₁ + (1 - alpha) * ωⱼ * ‖u‖₂, the row ℓ1/ℓ2 penalty by default (alpha = 0) and the multitask sparse group for alpha > 0. With row_norm = :linf, hⱼ(u) = ωⱼ * ‖u‖∞, and alpha must be zero. alpha must lie in [0, 1].

• The task weights cₜ are given by task_weights: :equal, the default, cₜ = 1, so that the mean squared error of each task counts once; :observations, cₜ = q nₜ / N with N = ∑ₜ nₜ, so that each observation counts once (the loss is q times the pooled loss ∑ₜ ‖y_t - β₀ₜ - X_t b_t‖² / (2N)); or a vector of positive weights. With equal sizes the first two are the same.
• With standardize_x, sₜⱼ is the standard deviation sqrt(∑ᵢ (X_t[i, j] - x̄ₜⱼ)² / nₜ) of column j within task t, centered at its mean x̄ₜⱼ with an intercept and at zero otherwise: the penalty is on the coefficients of the columns standardized within each task. Otherwise sₜⱼ = 1. Either way, a column that is zero in task t, or with an intercept constant there, has a zero coefficient in that task. The coefficients are on the scale of the data.
• The penalty factors ωⱼ are given by penalty_factors, finite, nonnegative and not all zero, one by default. They weight the row term only: a zero factor leaves the entrywise term when alpha > 0.
• The intercepts are not penalized, and they are zero unless intercept is true.

With unit penalty factors, one design X for all the tasks and equal task weights c, this is multi_lasso(X, [y_1 … y_q], λ / c), and with one task it is lasso(X_1, y_1, λ / c₁), both with the same intercept and standardize_x and unit factors; these reductions hold for both new penalties too.

For the default penalty, the tasks are stacked into one least-squares problem over vec(B), whose design has the centered and scaled designs of the tasks on its diagonal and is stored as a sparse matrix; a sparse design stays sparse without an intercept, and with an intercept its centered columns are dense. It is solved by block coordinate descent on the rows of B (ProximalSolvers' group coordinate descent), with options; with options.warm_start false the solution is computed along options.num_steps levels from λmax. The random order of coordinate descent is drawn from rng. On this problem momentum is on by default (see CDOptions): each pass over the active rows starts from a point extrapolated from the last two, which saves passes at small penalties; momentum = false or order = :cyclic turns it off. When the tasks share one dense design (the same matrix, or equal ones) and have equal weights, the model of multi_lasso is used instead, on that design and with the options as multi_lasso takes them, without momentum unless momentum = true. Real data of any types are converted to a common floating point type, and options are passed as CDOptions.

The row ℓ1/ℓ∞ and sparse-group penalties use row coordinate descent on the tasks' own designs, without stacking or copying a shared design. Sparse columns are centered on the fly. Full passes alternate with active-row passes. With relative_tol, the change rule measures each change Δ of U[j, t] as |Δ| * sqrt(aₜⱼ), where aₜⱼ = cₜ ‖x̃ₜⱼ‖² / (nₜ sₜⱼ²), against opt_tol * r, with response scale r = sqrt(∑ₜ cₜ ‖ỹₜ‖² / nₜ). The columns and responses are centered with an intercept. The scale has a positive numerical floor; at r = 0 the zero coefficients are exact. Without relative_tol, it uses absolute changes of U. Single fits require a duality gap only with a positive gaptol; continuation uses the path target. There is no gap at λ = 0, with completely unpenalized rows, or at zero objective.

Returns a MultitaskLassoSolution. See also MultitaskLassoPath and CVMultitaskLassoPath.

CVGLMLassoPath

Result returned by CVGLMLassoPath(X, y, loss).

Fields

• path::GLMLassoPath: the path on all the data
• λpath::Vector: its levels, the grid of every fold
• cvm::Vector, cvse::Vector: the cross-validated measure at each level, and its standard error (cvm is NaN at the levels that no fold reaches, cvse where fewer than two do)
• nnz::Vector{Int}: the number of nonzero coefficients of path at each level
• λmin, λ1se: among the levels that every fold reaches, the level with the smallest cvm, and the largest level whose cvm is within one standard error of it
• folds::Vector{Int}: the fold of each observation
• measure::Symbol: the measure of the held-out error

CVGLMLassoPath(X, y, loss = ProximalBase.Logistic(), options = CDOptions(); measure = :deviance,
               nfolds = 10, folds = nothing, rng = Random.default_rng(), threads = true, nλ = 100,
               λminratio, weights = nothing, offset = nothing, kwargs...)

K-fold cross-validation of the path of GLMLassoPath (with the same loss, options and keywords, such as intercept, standardize_x, alpha or penalty), as CVLassoPath does for least squares, with glmnet's measures of the held-out error (type.measure of cv.glmnet):

• :deviance, the default: 2 ∑ᵢ wᵢ (ℓ(ηᵢ, yᵢ) - minᵤ ℓ(u, yᵢ)) / ∑ᵢ wᵢ, for the logistic and Poisson losses;
• :mse and :mae: the mean squared and absolute errors of the predicted means μᵢ (probabilities for logistic regression);
• :class: the misclassification rate of μᵢ > 1/2, for logistic regression.

The path on all the data gives the levels; each fold's path is computed on the other folds at those levels, with their observation weights and offset, and the measure is computed on the held-out fold, weighted by its weights and with its offset in the linear predictor. cvm averages the folds, weighted by their sums of weights (their sizes without weights), and cvse is its standard error, λmin minimizes cvm and λ1se is the largest level within cvse of the minimum. A fold path that stops early (max_hat_s, devmax) leaves its remaining levels out of the averages: cvm is NaN at the levels that no fold reaches and cvse where fewer than two do, and λmin and λ1se are chosen among the levels that every fold reaches. The folds and threads are as for CVLassoPath: the result does not depend on the number of threads.

Returns a CVGLMLassoPath; coef(cv) and predict(cv, X) use the solution at λ1se, or at λmin with λ = :λmin.

CVLassoPath

Result returned by CVLassoPath(X, y).

Fields

• path::LassoPath: the path on all the data
• λpath::Vector: its levels, the grid of every fold
• cvm::Vector, cvse::Vector: the cross-validated mean squared prediction error at each level, and its standard error (cvm is NaN at the levels that no fold reaches, cvse where fewer than two do)
• nnz::Vector{Int}: the number of nonzero coefficients of path at each level
• λmin, λ1se: among the levels that every fold reaches, the level with the smallest cvm, and the largest level whose cvm is within one standard error of it
• folds::Vector{Int}: the fold of each observation

CVLassoPath(X, y, options = CDOptions(); nfolds = 10, folds = nothing, rng = Random.default_rng(),
            threads = true, nλ = 100, λminratio, intercept = false, standardize_x = true, kwargs...)

K-fold cross-validation of the lasso path of LassoPath (with the same options and keywords): the path on all the data gives the levels, each fold's path is computed on the other folds at those levels, and the mean squared error on the held-out fold estimates the prediction error at each level. Standardization and the intercept are computed on each training fold, so nothing leaks from the held-out one.

cvm is the average of the fold errors weighted by the fold sizes, and cvse its standard error, sqrt(∑ₖ nₖ (eₖ - cvm)² / n / (K - 1)), as in glmnet. λmin minimizes cvm, and λ1se is the largest level whose cvm is within cvse of the minimum (the one-standard-error rule). A fold path that stops early (max_hat_s) leaves the remaining levels out of the averages: cvm is NaN at the levels that no fold reaches and cvse where fewer than two do, and λmin and λ1se are chosen among the levels that every fold reaches.

The folds are a balanced random assignment drawn from rng, or the vector folds of fold labels 1, …, K given, for example grouped or blocked folds for dependent data. With threads, the folds run as tasks on the available threads; each fold draws its random coordinate orders from a generator seeded from rng, so that the result does not depend on the number of threads.

Observation weights and an offset (see lasso) follow the observations into the folds. The fold errors are then the weighted mean squared errors of y - offset, and the folds are weighted by their sums of weights. With gram = :auto, the choice for the full data is used for every training fold, so that holding out observations does not switch a nearly square design from the Gram form to the residual form. gram = true or false selects the form explicitly for all fits.

Returns a CVLassoPath; coef(cv) and predict(cv, X) use the solution at λ1se, or at λmin with λ = :λmin.

ebic(path::LassoPath, X, y; γ = 0.5)
aic(path::LassoPath, X, y)
bic(path::LassoPath, X, y)

Information criteria along the lasso path fitted to X and y, one per level of path.λpath:

• aic = n log(RSS/n) + 2 df;
• bic = n log(RSS/n) + log(n) df;
• ebic = bic + 2γ log C(p, s), the extended BIC of Chen and Chen (2008) for large model spaces, with s nonzero coefficients among p.

RSS is the residual sum of squares and df the number of nonzero coefficients plus one for an intercept, an unbiased estimate of the degrees of freedom of the lasso (Zou, Hastie and Tibshirani, 2007). A level with df ≥ n gets NaN. The level that minimizes a criterion c is path.λpath[argmin(c)], with the NaNs left out.

aic(path::LassoPath, X, y)

The Akaike criterion n log(RSS/n) + 2 df at each level of the lasso path; see ebic.

bic(path::LassoPath, X, y)

The Bayesian criterion n log(RSS/n) + log(n) df at each level of the lasso path; see ebic.

coef(cv::CVGLMLassoPath; λ = :λ1se)
predict(cv::CVGLMLassoPath, X; λ = :λ1se, type = :response, offset = nothing)

The coefficients ([β₀; β] with an intercept), and the predictions at the rows of X on the scale of the response or of the linear predictor (see glm_lasso), of the path on all the data at cv.λ1se, at cv.λmin with λ = :λmin, or at a level λ within the path.

coef(cv::CVLassoPath; λ = :λ1se)
predict(cv::CVLassoPath, X; λ = :λ1se)

The coefficients ([β₀; β] with an intercept), and the predictions at the rows of X, of the path on all the data at cv.λ1se, at cv.λmin with λ = :λmin, or at a level λ within the path.

CVGeneralizedLassoPath(X, y, D, levels; nfolds = 10, folds = nothing,
                       rng = Random.default_rng(), threads = true, kwargs...)

Prediction-error cross-validation on a common explicit penalty grid, using GeneralizedLassoPath. The same fixed coefficient penalty is used in each training fold; centering is learned only on its training observations. User fold labels allow blocked or grouped validation. Fields path, λpath, cvm, cvse, λmin, λ1se and folds follow CVLassoPath. Nonconverged fold fits raise an error instead of entering the CV averages.

GeneralizedLassoPath(X, y, D, levels; kwargs...)
GeneralizedLassoPath(y, D, levels; kwargs...)

Generalized-lasso fits at the supplied decreasing nonnegative penalty levels. ADMM warm starts its primal and scaled dual variables and reuses the normal equations. Consecutive fits with the same rho also reuse the factorization, including an automatically chosen value. The keywords are those of generalized_lasso. The field solutions holds each fit; λpath holds its penalty level. coef, predict, dof, deviance, aic and bic return one column or value per level. coef(path, λ) requires an exact grid level: this ADMM path does not compute or interpolate homotopy knots. Information criteria return NaN at unconverged fits or when df >= n.

StructuredLassoSolution

Fit from fused_lasso, generalized_lasso or trend_filter. x contains the coefficients, intercept is the unpenalized intercept, residuals are on the original data, and objective is sum(abs2, residuals)/(2n) + penalty. iterations, converged, primal_residual and dual_residual describe the solver termination. df estimates the generalized-lasso degrees of freedom. Use coef, predict, fitted, residuals, nobs, deviance and dof as for other fits.

cv_fused_lasso(X, y, levels; l1_ratio = 0.5, nfolds = 10, folds = nothing,
               rng = Random.default_rng(), threads = true, kwargs...)
cv_trend_filter(X, y, levels; order = 1, x = nothing, kwargs...)
cv_trend_filter(y, levels; kwargs...)

Cross-validation using fused_lasso_path or trend_filter_path. For signals, masked training observations keep the complete location grid. Fold work streams prediction errors without computing degrees of freedom or retaining a training path. The full-data returned path uses direct signal dispatch. max_concurrent_folds = 2 bounds simultaneous signal-fold work when threads = true; use one for serial fold memory. No held-out response enters a training fit or changes its difference penalty. order is the polynomial degree, and x must be strictly increasing. Returns a CVGeneralizedLassoPath.

cv_trend_filter(X, y, levels; order = 1, x = nothing, kwargs...)
cv_trend_filter(y, levels; kwargs...)

Cross-validate a trend-filtering path. Signal folds retain the full location grid and use diagonal observation Grams and eligible banded systems. Positive penalties require enough training locations to identify the degree-order polynomial component; unidentified folds fail without adding a ridge. See cv_fused_lasso for the fold, grid and tuning conventions.

difference_matrix(n, order = 1; x = nothing)

Sparse difference matrix with n - order rows. order = 0 gives the identity; order = 1 gives beta[j+1] - beta[j]. For strictly increasing locations x, higher differences use the divided-difference convention of trend filtering: D[k+1] = D[1] * diag(k / (x[j+k]-x[j])) * D[k]. Unit-spaced locations give ordinary successive differences.

fused_lasso(X, y, λ1, λ2; intercept = false, max_iter = 20_000,
            abstol = 1e-8, reltol = 1e-7, df_tol = 1e-4)
fused_lasso(y, λ1, λ2; kwargs...)

Minimize norm(y - intercept - X*beta)^2/(2n) + λ1*norm(beta,1) + λ2*sum(abs, diff(beta)). Regression uses accelerated proximal gradient with backtracking and ProximalBase's exact sparse fused-lasso proximal operator. The signal form is a single proximal step. Columns retain their original scales and order. The stopping test is the proximal-gradient residual in coefficient units, norm(beta - prox(beta - step*gradient))/step, against sqrt(p)*abstol + reltol*max(norm(X'y/n), 1). Before reporting convergence, the residual is checked at the returned coefficients using the original design; dual_residual reports that value. The defaults above apply to Float64 data; other floating types retain abstol = 1e-7, reltol = 1e-6. This is an optimality-residual test, not a uniform relative coefficient-error guarantee.

fused_lasso_path(X, y, levels; l1_ratio = 0.5, kwargs...)
fused_lasso_path(y, levels; kwargs...)

Warm-started fused regression at λ1 = l1_ratio*λ, λ2 = (1-l1_ratio)*λ for each decreasing value in levels. Keywords follow fused_lasso. Returns a GeneralizedLassoPath.

generalized_lasso(X, y, D, λ; intercept = false, rho = nothing,
                  max_iter = 20_000, abstol = 1e-9, reltol = 1e-8,
                  df_tol = 1e-4)
generalized_lasso(y, D, λ; kwargs...)

Minimize norm(y - intercept - X*beta)^2/(2n) + λ*norm(D*beta, 1) by ADMM. D can be dense or sparse, for example graph_difference. The signal form uses the identity design without allocating it. ADMM factors X'X + rho*D'D once; [X; D] must have full column rank. rho is on the unnormalized-loss scale and defaults to the larger of n*λ and the average column squared norm of the centered design (one for a signal fit).

No standardization is applied, since it would change the difference penalty. At zero penalty the minimum-norm least-squares solution is returned. The result reports nonconvergence if the residual stopping criteria are unmet. dof is rank(X * null(D[inactive, :])), including an intercept when fitted; inactive rows satisfy abs(D*beta) <= df_tol. At transition points the degrees-of-freedom estimate is sensitive to this numerical tolerance. Compared with R's genlasso, use lambda_R = n*λ. The default residual tolerances above apply to Float64 data; other floating types retain abstol = 1e-7, reltol = 1e-6. Residual convergence does not certify a uniform relative coefficient error on ill-conditioned designs.

graph_difference(p, edges; weights = nothing)

Sparse oriented incidence matrix of a graph on vertices 1:p. Each pair (i, j) in edges contributes the row weight * (beta[j] - beta[i]). Weights must be positive and finite; orientation does not affect the penalty.

trend_filter(y, λ; order = 1, x = nothing, kwargs...)
trend_filter(X, y, λ; order = 1, x = nothing, kwargs...)

Trend filtering with λ * norm(difference_matrix(p, order+1; x)*beta, 1) and loss norm(y-X*beta)^2/(2n). order = k gives piecewise polynomials of degree k (order zero is total-variation denoising). x gives the ordered coefficient locations, not the rows of the design.

For order > 0, the specialized ADMM of Ramdas and Tibshirani (2016) splits the order-k divided differences and applies exact TV denoising in each iteration. Sparse banded normal equations avoid a dense identity in signal fits. Keywords and termination diagnostics follow generalized_lasso. For a signal, default rho = n*λ*spacing^order follows the paper's spacing rule. With a design, it is at least the average column squared norm times spacing^(2order). R's glmgen::trendfilter uses lambda = n*λ.

Float64 linear signal fits (order = 1, positive reltol) stop using a feasible primal-dual gap at most 0.3reltol times the primal objective, with primal and dual residuals at most five times their usual bounds. The gap is checked every ten iterations and at the first and final iteration. Near a zero objective the original residual test is used. Other fits and absolute-only tolerances retain residual stopping. Other Float64 trend fits default to abstol = 1e-8, reltol = 1e-7. Linear signals and other floating types retain 1e-7, 1e-6. Explicit tolerances override these defaults, including on paths and CV folds.

trend_filter_path(X, y, levels; order = 1, x = nothing, kwargs...)
trend_filter_path(y, levels; kwargs...)

Warm-started specialized ADMM path for trend_filter, at supplied decreasing levels. Returns a GeneralizedLassoPath. For order zero, each signal fit is a single exact proximal step.

aic(path::GeneralizedLassoPath)
bic(path::GeneralizedLassoPath)

Gaussian information criteria n*log(RSS/n) + c*df, where c = 2 for AIC and c = log(n) for BIC. df = rank(X*null(D_inactive)), rather than a count of nonzero coefficients. Values are NaN for unconverged fits or df >= n. Omitted Gaussian constants are common to all levels.

PostLassoResult

Result returned by post_lasso: the least-squares coefficients x (a SparseVector, zero off the selected columns) and intercept, their heteroskedasticity-robust (HC1) standard errors se and se_intercept, the selected columns, the residuals, and vcov, the HC1 covariance of [β₀; β_S] (of β_S without an intercept), and hasintercept.

The functions of StatsAPI apply to it, with [β₀; β] (β without an intercept) as the coefficients: coef, stderror, vcov (a sparse matrix, zero off the selected columns), confint (normal intervals, zero width off the selected columns), predict, residuals, nobs, dof and deviance.

RLassoSolution

Result returned by rlasso. Its positional constructor is internal; obtain a solution by calling rlasso.

Fields

• x::SparseVector: the coefficients, refitted by least squares on the selected columns when post is true (as hdm's rlasso returns them), the lasso coefficients otherwise
• intercept: the intercept, zero without one
• selected::Vector{Int}: the selected columns, equal to findall(!iszero, lasso)
• λ0, loadings::Vector: the nominal penalty level on the scale of feasible_lasso!, and the loadings updated from the final residuals (unchanged on an empty-selection return; hdm's lambda is 2n λ0 Γ)
• residuals::Vector, σ: the residuals of the final coefficients and their standard deviation
• iterations::Int: the lasso fits of the procedure
• post::Bool
• hasintercept::Bool: whether the procedure fitted an intercept
• lasso::SparseVector, lasso_intercept: a copy of the penalized coefficients of the final coordinate-descent solve, and their intercept (zero without one)
• lasso_λ0, lasso_loadings::Vector: the penalty level and a copy of the loadings actually used by that solve, before the residual-based update; lasso_λ0 == λ0/2 if the final solve was the first one and post is true, and lasso_λ0 == λ0 otherwise
• converged::Bool: the outer standard-deviation stopping criterion was met; an empty-selection return also counts as converged
• solver_converged::Bool: every coordinate-descent solve converged, including the solve preceding an empty-selection return
• rank::Int: the numerical rank of the selected columns after centering when fitting an intercept, zero for no selected columns. For BLAS floating types this is the rank returned by the pivoted-QR least-squares solve used by _ols: incremental estimates of the smallest/largest singular values of leading R blocks, with rcond = min(n, length(selected)) * eps(T). A block is rejected when smax * rcond > smin. This is not the diagonal threshold used by rank(::QRPivoted). Other floating types follow Julia's generic pivoted-QR solve, which assumes full rank on a successful solve.
• tuning::NamedTuple: (c, γ, max_iter, tol, options) as used by the fit. max_iter limits the outer loadings loop; options.max_iter limits the coordinate-descent passes of each solve

coef, predict, residuals, nobs, dof and deviance apply to it.

TreatmentEffect

Result returned by double_selection for the prespecified treatment coefficient in original units. The positional constructor is internal.

estimate, se, t, pvalue, and ci give the estimate, standard error, normal statistic, two-sided normal p-value, and interval at the requested level. selected is the union of the nuisance supports and, for double selection, the forced controls. method, n, vcov_type, and post record the procedure. refit records whether partialling used OLS nuisance residuals and is nothing for double selection. selected_y, selected_d, and forced record the separate supports and the controls forced into double selection (empty for partialling).

rank is the rank of the intercept and selected-control basis (1 for partialling); dof_residual is the final regression's residual degrees of freedom. target_variation is the fraction of centered treatment variation remaining after projection. max_leverage and the observation-level influence refer to the final regression. nuisance_y and nuisance_d retain the RLassoSolutions, including their outer and solver convergence flags. Non-convergence is reported, without rejecting the effect estimate.

coef, stderror, vcov and confint (at any level) give the estimate as a model with one coefficient; nobs is n.

double_selection(y, d, X; method = :doubleselection, vcov = :hdm, post = true,
                 refit = nothing, forced = Int[], level = 0.95, options = CDOptions(), rng = default_rng(), kwargs...)
double_selection(y, d, X, fit_y::RLassoSolution, fit_d::RLassoSolution;
                 method = :doubleselection, vcov = :hdm, refit = nothing, forced = Int[], level = 0.95)

Inference for the prespecified coefficient α in y = μ + αd + Xβ + ε, in original units, with an intercept throughout. This procedure does not justify choosing the target after inspecting estimated effects. kwargs (c, γ, max_iter, tol) go to rlasso; options controls its coordinate descent and rng supplies both fits' random orders. Double selection fits d then y; partialling fits y then d. intercept=false, weights, and offsets are unsupported.

Double selection uses the sorted union of the d-on-X and y-on-X supports (the latter is the reduced form) and forced. Forced controls enter only the final regression, as hdm's I3, without changing either nuisance fit. forced must be an AbstractVector with an integer element type, excluding Bool, unique indices in 1:p, and no Bool elements. Invalid values raise ArgumentError before fitting; pass column indices, not a Bool mask. The control basis W contains an intercept and centered selected controls. With τ = max(1e-7, sqrt(eps(T))), a centered column of norm at most τ times its raw norm is aliased with the intercept. Remaining columns are scaled to unit norm; pivoted QR retains the r directions with abs(R[j,j]) > τ. Projection gives v = M_W d, w = M_W y, base leverages h_W, k = r+2, and nominal control count s, including forced and aliased columns. Rank-deficient controls are allowed when the treatment is identified; the target remains the original coefficient of d, and selection does not establish true-support recovery. This rule differs from R's dqrdc2 and from the nuisance QR solve. Partialling centers the chosen nuisance residuals to obtain v and w, with h_W = 1/n, k = 2, and s = 0. refit::Union{Nothing,Bool}=nothing defaults to the fits' post. If refit == post, their native residuals are used without recomputation. With post=true, refit=false, residuals use the stored final penalized solve: response - lasso_intercept - X*lasso. With post=false, refit=true, OLS residuals use the same pinned projection rule as double selection, separately on each nuisance's selected controls. These modes keep the fits, supports, penalty levels and loadings fixed. Passing precomputed fits permits controlled comparisons; changing native post also changes the default penalty constant from 1.1 to 0.5.

Partialling rejects nonempty forced (native hdm ignores I3 here). When effective refit=true, each nuisance must have length(selected)+1<n, even if some columns are aliases. Double selection only accepts refit=nothing; its final OLS is part of the estimator. HC corrections for partialling always use k=2; they do not count the nuisance controls as final-stage parameters.

In both layouts, α̂ = dot(v,w)/dot(v,v), e = w - α̂v, φᵢ = n*vᵢ*eᵢ/dot(v,v), and hᵢ = h_Wᵢ + vᵢ²/dot(v,v). The variance is:

The default preserves hdm 0.3.2's variance conventions. Intervals are α̂ ± z*se; t statistics and two-sided p-values use the standard normal.

Inputs must be finite, with n≥3, p≥1, 0<level<1, and nonconstant y and d. All variances require k<n; hdm double selection also requires n-s-1≥1. Reject norm(v) ≤ τ*norm(d-mean(d)), norm(e) ≤ τ*norm(y-mean(y)), and, for HC2/HC3, maximum(h) > 1-τ. Invalid quantities raise ArgumentError with their values; incompatible dimensions raise DimensionMismatch.

The precomputed method consumes no RNG and does not mutate either fit. It checks dimensions, promoted element type, intercepts, matching post flags, sorted unique lasso supports, and residual consistency with y/d and X to a budget sqrt(eps(T)) * (norm(response,Inf) + abs(b0) + maximum(abs.(X)*abs.(x))) for norm(residuals - (response - b0 - X*x), Inf). These checks establish consistency, not provenance: they cannot establish which data or algorithm produced a record. Nuisance non-convergence is retained in the result and displayed by show.

feasible_lasso!(x, X, y, options = IterLassoOptions(); λ0 = :theory, c = 1.1, γ = 0.1 / log(n), intercept = false)

feasible_lasso! at the penalty level of theory_lambda(n, p; c, γ). R's hdm uses c = 0.5 instead when it does not refit by least squares (rlasso(…, post = FALSE)); see rlasso for hdm's procedure.

post_lasso(X, y, sol; intercept = true)

Least squares of y on the columns that the lasso sol selected (a LassoSolution, an RLassoSolution, or a vector of column indices), with an intercept unless intercept is false, and the heteroskedasticity-robust standard errors HC1, n/(n - k) ⋅ (Z'Z)⁻¹ Z' Diagonal(e²) Z (Z'Z)⁻¹ for the design Z = [1 X_S] with k columns. Returns a PostLassoResult.

rlasso(X, y, options = CDOptions(); post = true, intercept = true, c = post ? 1.1 : 0.5,
       γ = 0.1 / log(n), max_iter = 15, tol = 1e-5, rng = default_rng())

The lasso with data-driven penalty loadings as R's hdm computes it (rlasso with its defaults: heteroskedastic and X-independent penalty), the procedure of Belloni, Chen, Chernozhukov and Hansen (2012):

1. Start from the residuals of least squares on the five columns most correlated with y, and compute the loadings Γⱼ from them (see feasible_lasso!).
2. Solve the lasso ‖y - β₀ - Xβ‖² / (2n) + λ0 ∑ⱼ Γⱼ |βⱼ| with λ0 = theory_lambda(n, p; c, γ), at half the level in the first step when post is true, as hdm does.
3. With post, refit least squares on the selected columns; update the loadings from the residuals; repeat 2–3 until the standard deviation of the residuals changes by less than tol, or max_iter times.

The coefficients returned are the refitted ones when post is true. With intercept, the columns and y are centered and β₀ = ȳ - x̄'β. Each lasso is solved by coordinate descent with options, warm started from the previous one. Random coordinate orders are drawn from rng::AbstractRNG; a private generator leaves the default RNG unchanged. The default preserves the stream of calls seeded with Random.seed!.

Returns an RLassoSolution, retaining both the final penalized solve (lasso, lasso_intercept, lasso_λ0, lasso_loadings) and the returned coefficients/residuals, together with outer/solver convergence, selected-design rank, and effective tuning. The outer max_iter limits loadings iterations; options.max_iter limits coordinate-descent passes per solve. See also post_lasso and double_selection.

theory_lambda(n, p; c = 1.1, γ = 0.1 / log(n))

The penalty level λ0 = c Φ⁻¹(1 - γ/(2p)) / √n of Belloni, Chen, Chernozhukov and Hansen (2012) for the lasso with data-driven loadings,

minimize  ‖y - β₀ - Xβ‖² / (2n) + λ0 ∑ⱼ Γⱼ |βⱼ|,   Γⱼ = sqrt(∑ᵢ (Xᵢⱼ - x̄ⱼ)² εᵢ² / n),

the scale of feasible_lasso!. It is R's hdm penalty 2c √n Φ⁻¹(1 - γ/(2p)) for ‖y - Xβ‖² + λ ∑ⱼ Γⱼ |βⱼ|, divided by 2n, with hdm's defaults for c and γ. Φ⁻¹ is the standard normal quantile.

CVLowRankPath(X, Y; nfolds = 5, folds = nothing, rng = Random.default_rng(), kwargs...)

K-fold cross-validation of LowRankPath, on a common penalty or rank grid. Centering, adaptive weights and initialization are recomputed on each training fold. Errors are mean squared errors per response, weighted by fold size. cvm, cvse, folds, path and index (minimum-error solution) are returned. coef and predict use that solution. A failed fold solve is an error; it is never silently selected. For rank paths choose ranks estimable in every training fold (an invalid rank raises an error).

LowRankPath(X, Y; method = :nuclear, λpath = nothing, ranks = nothing,
            rank = 1, nλ = 50, λminratio = 0.01, tol = nothing,
            inner_tol = nothing, kwargs...)

Fit a path with method = :nuclear, :sparse, :rank, or :sparse_svd. Nuclear and sparse paths use decreasing λpath (or an automatic log grid); :sparse holds rank fixed. :rank fits the supplied ranks, by default zero through the maximum estimable rank. :sparse_svd fits the supplied rank bounds with component penalties chosen by BIC. solutions contains independent LowRankSolutions. Nuclear paths reuse the gradient setup and warm start the coefficient matrix. solver = :panoc selects the optional ProximalAlgorithms backend for nuclear paths. block_solver = :fista selects block FISTA for sparse paths; defaults are solver = :apg and block_solver = :cd. Sparse paths warm start nonzero fits, restarting at the RRR initialization after a zero fit to avoid an absorbing zero block. Adaptive weights are fixed from the training-data RRR fit over the whole path. Select using argmin(lowrank_ic(path)), or CVLowRankPath. For sparse paths, tol and inner_tol default to max(1e-7, 100eps(T)) and max(1e-9, 100eps(T)) (1e-7 and 1e-9 in Float64); other methods use 1e-7 and 1e-9. Explicit tolerances are unchanged. Each sparse block finishes with guarded Newton steps on its nonzero rows when the first-order solver needs more than 100 first-order KKT checks (100 CD sweeps or 1000 FISTA iterations), and thereafter in that fit. inner_max_iter bounds first-order iterations; Newton steps are additional. Paths do not warn on non-convergence; sparse_reduced_rank does.

LowRankSolution

A multivariate regression fit. λ is its penalty level, B the coefficient matrix, intercept the response intercepts, rank the numerical rank of B, and activerows the selected predictors. iterations, converged and objective describe the optimization (:sparse_svd reports the combined residual loss); a nonconvex fit reports convergence to a block stationary point, not a global minimum. coef, predict, residuals, nobs, deviance, dof, aic and bic apply. dof counts the parameters of the selected rank/support model; it is not a Stein unbiased degrees-of-freedom estimate.

lowrank_ic(solution; criterion = :bic, df = :naive)
lowrank_ic(path; criterion = :bic)

Gaussian model-selection scores log(RSS) + penalty * dof/(n*q); :aic uses 2, :bic uses log(n*q), and :bicp uses 2log(n*q) for row-sparse regression and 2log(p*q) otherwise, as in rrpack. Degrees of freedom count the selected rank/support model. Smaller is better. Penalized fits use a parameter-count approximation rather than exact Stein degrees of freedom. df = :exact uses the Stein formula for rank-constrained fits (Mukherjee et al., 2015); a cutoff through tied singular values has infinite score.

nuclear_regression(X, Y, λ; solver = :apg, intercept = true, max_iter = 10_000, tol = 1e-7)

Minimize ‖Y - 1β₀' - XB‖²_F / (2n) + λ ‖B‖_* by accelerated proximal gradient, using ProximalBase.NuclearNorm for singular-value thresholding. Columns remain on their original scales. tol bounds the Frobenius norm of the proximal gradient mapping relative to max(1, ‖Xc'Yc/n‖_F). See LowRankPath for warm starts over several penalties. solver = :panoc uses ProximalAlgorithms' PANOC with the same loss, nuclear proximal operator, and stopping certificate; first load ProximalAlgorithms.

reduced_rank(X, Y, rank; intercept = true, rtol = nothing)

Minimize ‖Y - 1β₀' - XB‖²_F subject to rank(B) ≤ rank. The truncated SVD of the least-squares fitted values gives the solution (Anderson, 1951; Izenman, 1975). A compact SVD of X handles wide and rank-deficient designs and chooses the minimum-norm coefficient matrix. rtol is the relative singular-value cutoff for X, by default max(size(X)...) * eps(T). The intercept is unpenalized; columns are kept on their original scales.

select_rank(X, Y; sigma2 = nothing, multiplier = 2, μ = nothing, intercept = true)

The rank selection criterion of Bunea, She and Wegkamp (2011): minimize RSS + μ rank(B). Keep fitted singular values whose squares exceed μ. By default μ = multiplier * sigma2 * (√rank(Xc) + √q)^2, with centered Xc when fitting an intercept. If omitted, sigma2 is estimated from the full least-squares residuals using q * (n - rank(Xc) - intercept) degrees of freedom. A saturated design requires a supplied variance or μ. Returns (solution, rank, μ, sigma2, singularvalues).

sparse_reduced_rank(X, Y, rank, λ; intercept = true, weights = nothing,
                    adaptive = false, gamma = 2, max_iter = 1000, tol = nothing,
                    inner_max_iter = 10_000, inner_tol = nothing, block_solver = :cd)

The row-sparse reduced-rank regression of Chen and Huang (2012): minimize ‖Y - 1β₀' - XB‖²_F/(2n) + λ ∑ⱼ weights[j] ‖B[j,:]‖₂, with rank at most rank. Alternates a group-lasso fit of A and an orthogonal Procrustes update of V in B = AV'. It starts from reduced-rank least squares. This problem is nonconvex; convergence is local. adaptive = true uses inverse row norms of the initial RRR fit to power gamma, with a sqrt(eps(T)) floor. Explicit positive weights override adaptive weights. No column rescaling is performed. For rrpack's unscaled loss use λ_rrpack = n * λ. block_solver = :fista replaces row coordinate descent with restarted block FISTA, using the same weighted group proximal map and KKT stopping rule. When a block needs more than 100 first-order KKT checks (100 CD sweeps or 1000 FISTA iterations), that block and later blocks in the fit finish with guarded Newton steps on the nonzero rows. inner_max_iter bounds first-order iterations; Newton steps are additional. The defaults for tol and inner_tol are max(1e-7, 100eps(T)) and max(1e-9, 100eps(T)), respectively (1e-7 and 1e-9 in Float64). Explicit tolerances are unchanged. A fit that does not converge emits a warning.

sparse_svd_regression(X, Y, rank; intercept = true, nλ = 50,
                      λminratio = 1e-5, adaptive = true, gamma = 2,
                      criterion = :bic, max_iter = 1000, tol = 1e-7,
                      inner_max_iter = 10_000, inner_tol = 1e-9)

Sparse singular-vector regression by one round of iterative exclusive extraction (Chen, Chan and Stenseth, 2012; rrpack's rssvd). From the SVD of the rank-constrained least-squares coefficient matrix, form the exclusive response layer for each component. Fit each layer with a weighted rank-one lasso, selecting its penalty by :bic, :bicp or :aic, including the zero component. Both predictor and response loadings can be sparse. The rank is an upper bound; selected layers need not remain mutually orthogonal.

Adaptive inverse weights use the absolute initial singular vectors to power gamma, with zero entries replaced by eps(T). Each component's grid runs from 0.9λmax to λminratio*λmax, followed by zero, where λmax = maximum(abs.(wu .* (Xc'Ylayer/n) .* wv')). The loss is divided by 2n. Optimization is nonconvex and convergence is local. Returns (solution, components, λ, scores, convergedpath, selected); the scores include every grid value and the zero model. Nonconverged candidates have infinite scores and their flags are preserved in convergedpath. Component scores use nnz(u)+nnz(v)-1 parameters.

dirty_multitask(Xs, ys, λ_shared, λ_specific, options = CDOptions(); shared_norm = :linf,
                intercept = true, standardize_x = true, task_weights = :equal,
                penalty_factors = nothing, rng = Random.default_rng())
dirty_multitask(X, Y, λ_shared, λ_specific, options = CDOptions(); kwargs...)

The dirty model of Jalali et al. (2010), with total coefficients B = A + S,

minimize  ∑ₜ cₜ ‖yₜ - β₀ₜ - Xₜ(aₜ + sₜ)‖² / (2nₜ)
            + ∑ⱼ ωⱼ (λ_shared ‖Ã[j, :]‖∞ + λ_specific ‖S̃[j, :]‖₁).

Here Ã[j, t] = scale[j, t] * A[j, t], and likewise for S̃, with scales, task weights and unpenalized intercepts as in multitask_lasso. Xs[t] and ys[t] may have unequal task sizes; X, Y share one design, a response per column. shared_norm = :l2 replaces the shared norm by ℓ2, a variant of the paper's model. Levels must be finite and nonnegative. Penalty factors multiply both terms; a zero factor leaves a row unpenalized.

The solver minimizes over B using the infimal-convolution row penalty. For ℓ∞ this is SLOPE with weights ωⱼ min(λ_specific, max(λ_shared - (i - 1)λ_specific, 0)) on the fitted coordinates. No stacked design is built; sparse designs are centered on the fly. options controls the row solver, including continuation when warm_start = false.

Components are split in standardized coordinates, then rescaled. A zero shared level gives A = B; a zero specific level gives S = B; both zero give A = B. At positive levels, the lasso endpoint is handled first (A = 0), then the row-penalty endpoint (S = 0). Interior splits clip the standardized row at its ceil(λ_shared / λ_specific)-th magnitude for ℓ∞ (the most-shared choice at integer ratios), or at the positive root described in the guide for ℓ2. The split is canonical and need not be identified by the data.

Returns DirtyMultitaskSolution. See DirtyMultitaskPath and CVDirtyMultitaskPath for warm starts and tuning.

DirtyMultitaskPath(Xs, ys[, λpath], options = CDOptions(); ratio, shared_norm = :linf,
                   nλ = 100, λminratio = nothing, screening = true, max_rows = Inf, kwargs...)
DirtyMultitaskPath(X, Y[, λpath], options = CDOptions(); ratio, kwargs...)

A decreasing path of specific levels λ, with shared levels ratio * λ. ratio must be positive and finite. Total coefficients warm start each level; components are recomputed. The automatic grid has nλ log-spaced levels from λmax * (1 + √eps) to λminratio times that value. λminratio defaults to 0.01 if sum(nₜ) < pq, and 1e-4 otherwise.

For gradient gⱼ on fitted coordinates, λmax = maxⱼ dⱼ(gⱼ(0)), where dⱼ(g) = max(‖g‖₁ / ratio, ‖g‖∞) / ωⱼ for ℓ∞, and replaces ℓ1 by ℓ2 for the variant. A nonzero gradient on an unpenalized row requires explicit levels. Screening uses this gauge in the sequential strong rule and checks all excluded rows. max_rows stops after the first fit exceeding that many nonzero rows.

λpath, Bpath, Apath, Spath, β0path, activerows, sharedrows, objective, gap, kkt, iterations and converged hold each level's results; kkt_checks counts screening repair rounds. ratio, shared_norm, intercept, nobs, task_nobs, task_weights and penalty_factors record the fit settings. coef(path) gives a matrix per level; coef(path, λ) requires an exact path level, preserving its decomposition and certificate. predict(path, Xs) gives an nₜ × L matrix per task; predict(path, X) gives an n × q matrix per level.

CVDirtyMultitaskPath(Xs, ys, options = CDOptions(); ratios = nothing, shared_norm = :linf,
                     nfolds = 10, folds = nothing, rng = Random.default_rng(), threads = true,
                     nλ = 100, kwargs...)
CVDirtyMultitaskPath(X, Y, options = CDOptions(); kwargs...)

Tune the ratio and specific level by held-out prediction error. The default uses five geometric ratios on [1, q] for ℓ∞, or [1, √q] for ℓ2. All ratios use the same per-task folds. The shared-design method holds out whole rows across responses. Each training fit refits its preprocessing and task weights; its grid is the corresponding full-data path's grid.

The score is that of CVMultitaskLassoPath, weighted by full-data task weights. Nonconverged fold levels are unavailable. Selection compares only levels reached and converged in every fold and on the full data. Ties go to the larger ratio, then the larger specific level; there is no 1se rule.

paths holds the full-data paths and ratios their ratios. cvm and cvse are length(ratios) × nλ matrices, padded with NaN for shorter paths; cvm_task[t, i, j] is task t's held-out MSE at ratio i, level j. Means use the folds that reach and converge at a level. selected is the selected CartesianIndex(ratio_index, level_index), ratio and λmin its settings, path its full-data path, and solution its converged full-data fit. folds stores the labels per task, common across responses for a shared design. coef(cv) and predict(cv, Xnew) use solution.

sparse_nuclear_regression(X, Y, λ_rank, λ_sparse; sparsity = :entry,
                          weights = nothing, intercept = true,
                          max_iter = 10_000, tol = 1e-7)

Minimize ‖Y - 1β₀' - XB‖²_F/(2n) + λ_rank ‖B‖_* + λ_sparse P(B) by ProximalSolvers' Condat–Vũ splitting. Both penalties act on the same B. For sparsity = :entry, P(B) = ∑ⱼ weights[j] ∑ₜ |B[j,t]|; for :row, P(B) = ∑ⱼ weights[j] ‖B[j,:]‖₂. Row weights are positive and finite, one by default. Levels are finite and nonnegative. Columns keep their original scales, and intercepts are unpenalized.

Returns a SparseNuclearSolution. tol bounds its relative duality gap at the returned coefficients. With both levels zero, the solver's KKT residual rule is used and gap = NaN. A constant centered design gives the intercept-only fit. See PenaltyGrid for warm starts and CVPenaltyGrid for tuning both levels.

lowrank_plus_sparse(X, Y, λ_rank, λ_sparse; sparsity = :entry,
                    weights = nothing, intercept = true,
                    max_iter = 10_000, tol = 1e-7)

Minimize ‖Y - 1β₀' - X(L + S)‖²_F/(2n) + λ_rank ‖L‖_* + λ_sparse P(S) by restarted accelerated proximal gradient on (L, S), with step 1/(2‖Xc‖²_op/n) and ProximalBase's nuclear and sparsity proximal maps. For sparsity = :entry, P(S) = ∑ⱼ weights[j] ∑ₜ |S[j,t]|; for :row, P(S) = ∑ⱼ weights[j] ‖S[j,:]‖₂. Weights are positive finite predictor weights, one by default. Observation, task and domain weights are unsupported. Levels are finite and nonnegative. Intercepts are unpenalized, columns keep their original scales, and the loss sums over responses without dividing by their number. Only a shared design matrix is supported.

Returns a LowRankPlusSparseSolution. With positive levels, tol bounds the relative intersection-dual gap at the returned pair, up to floating-point accuracy. A small gap does not bound component error. If λ_rank = 0, the minimum-norm least-squares coefficient is assigned to L and S = 0, including when both levels are zero. If only λ_sparse = 0, it is assigned to S and L = 0. These axes have gap = NaN and use the least-squares normal equations. A zero centered design gives an intercept-only fit. Components need not be unique, and no identifiability constraints are imposed. See PenaltyGrid for warm starts and CVPenaltyGrid for tuning both levels; automatic axes use the full-data zero-pair references.

PenaltyGrid(f, data...; λ1 = nothing, λ2 = nothing, nλ = (10, 10),
            λminratio = (0.01, 0.01), kwargs...)

Fit an estimator implementing the internal two-penalty protocol, currently sparse_nuclear_regression. solutions[i,j] uses levels λ1[i] and λ2[j], whose meanings are recorded in names ((:rank, :sparse)). Explicit axes are nonempty, decreasing, finite and nonnegative; zeros and repeated levels are allowed. Automatic axes descend logarithmically from the respective marginal zero thresholds times 1 + √eps.

Fits proceed down columns, starting each column from the head of the previous one. Warm states and returned solutions are independent. coef and predict return matrices of results over the grid.

CVPenaltyGrid(f, data...; nfolds = 5, folds = nothing,
              rng = Random.default_rng(), threads = true, kwargs...)

K-fold validation of PenaltyGrid. Every training fold refits preprocessing on its own observations, using the numerical axes of the full-data grid. Each fold has its own seeded generator and workspaces. For sparse_nuclear_regression, cvm is mean squared prediction error per response, weighted by held-out row counts. cvse is sqrt(∑ₖ wₖ (errorₖ - cvm)^2 / (∑ₖ wₖ (K-1))).

A pair is eligible only if the full-data fit and every fold fit converge; otherwise cvm and cvse are NaN. An empty eligible set raises an error. index selects the minimum cvm, breaking ties in column-major order; λmin = (grid.λ1[index[1]], grid.λ2[index[2]]). There is no two-dimensional one-standard-error rule. coef and predict use this full-data solution.

transfer_multitask(target, sources, λ₀, options = CDOptions(); contrast_weights = :sqrt,
                   intercept = true, solver = :cd, admm = (;), rng = Random.default_rng())

TMTL(Fused) of Lim and Kolar (2026): joint-sparse transfer from source domains ℓ = 1, …, L to a target domain ℓ = 0. Domain ℓ has nₗ observations of the same K responses on the same p predictors,

minimize  ∑ₗ ∑ₖ ‖yˡᵏ - β₀ˡᵏ - Xˡᵏ Bˡ[:, k]‖² / (2NK)
            + λ₀ ∑ⱼ ‖B⁰[j, :]‖₂ + λ₀ ∑ₗ≥₁ aₗ ∑ⱼ ‖Bˡ[j, :] - B⁰[j, :]‖₂

over the p × K domain blocks B⁰, …, Bᴸ and unpenalized intercepts β₀ˡᵏ per domain and response, with N = ∑ₗ nₗ over all domains, the target included. The target block is row sparse; the contrasts Bˡ - B⁰ are row sparse, the sources' blocks need not be. The reported estimate is the fused Ŵᶠ = ∑ₗ (nₗ/N) B̂ˡ, not B̂⁰; target predictions use the intercepts ȳ⁰ - x̄⁰'Ŵᶠ. debias adds the target correction (TMTL(Debiased)).

target = (X⁰, Y⁰) and sources is a vector of (Xˡ, Yˡ), possibly empty. Yˡ is nₗ × K; Xˡ is one nₗ × p design for all responses, or a vector of K designs, one per response (task-specific designs, Appendix B). Sizes may differ across domains; the responses of a domain share its nₗ rows. Designs are not standardized, since the contrasts compare coefficients across domains on the predictors' own scales. With intercept, columns and responses are centered within each domain and response. contrast_weights is :sqrt (aₗ = 8√(nₗ/N)), :linear (aₗ = 8nₗ/N) or a vector of L positive weights. λ₀ must be finite and nonnegative.

The default solver is coordinate descent in the reparametrization θ⁰ = B⁰, θˡ = Bˡ - B⁰ (Appendix C.1), whose penalty λ₀‖θ⁰‖₂,₁ + ∑ₗ λ₀aₗ‖θˡ‖₂,₁ separates over components and rows. A row of θ⁰ enters every domain, and its gradient and curvature are summed over them; a row of θˡ enters domain ℓ. Each row step uses the largest coordinate curvature of the row, exact for a shared design and a majorization otherwise. A coefficient whose column is constant in one domain can be nonzero, since the others determine it; only coefficients that enter no loss stay zero. options are those of CDOptions, as for the row solver of multitask_lasso, including continuation without warm starts; a positive gaptol requires the global relative duality gap described in the guide.

solver = :admm uses the paper's two-block ADMM (Appendix A) with the exact pair update of eq. (37), cached least-squares factorizations, residual balancing for the first adapt iterations, and its constraint-wise primal and dual residual stopping rule. admm overrides the settings (; max_iter = 10_000, abstol = 1e-8, reltol = 1e-6, rho = 1.0, rho_max = 1e4, balance = 10.0, increase = 2.0, decrease = 1.5, adapt = 1_000); constraint c stops when its primal residual is at most √(pK) abstol + reltol ⋅ max(‖B‖, ‖Γ‖) and its dual residual at most √(pK) abstol + reltol ⋅ ρ‖U‖. ADMM reports θ⁰ = Γ⁰⁰ and the contrasts Γˡ - Γ⁰ˡ, so zero rows are exact. options and rng are not used by ADMM.

Returns a TransferMultitaskSolution. See TransferMultitaskPath and CVTransferMultitaskPath for paths and tuning.

debias(fit::TransferMultitaskSolution, target, λ̃, options = CDOptions();
       rng = Random.default_rng())

TMTL(Debiased): the target correction of the fused fit Ŵᶠ = fit.fused,

D̂ = argmin  ∑ₖ ‖y⁰ᵏ - d₀ᵏ - X⁰ᵏ(Ŵᶠ[:, k] + D[:, k])‖² / (2n₀K) + λ̃ ∑ⱼ ‖D[j, :]‖₂,

with unpenalized intercepts d₀ᵏ when the fit has them, and Ŵᴰ = Ŵᶠ + D̂. It is a group lasso of the target residuals Y⁰ - X⁰Ŵᶠ on the unstandardized target design, the same as multi_lasso(X⁰, Y⁰ - X⁰Ŵᶠ, K λ̃; standardize_x = false) for a shared design. "Debiased" refers to the aggregation bias of Ŵᶠ in point estimation, not to inference. The correction uses the same target data as the fused fit, which target = (X⁰, Y⁰) must give again; λ̃ must be finite and nonnegative. A column constant in the target has a zero correction.

The row solver of multitask_lasso computes D̂, with options and its relative duality gap. The fused diagnostics are kept; the correction's are reported separately. A corrected fit is corrected again from fit.fused. Returns a TransferMultitaskSolution with W = Ŵᴰ and target intercepts ȳ⁰ - x̄⁰'Ŵᴰ.

TransferMultitaskPath(target, sources, options = CDOptions(); nλ = 100, λminratio = nothing, kwargs...)
TransferMultitaskPath(target, sources, λpath, options = CDOptions(); contrast_weights = :sqrt,
                      intercept = true, screening = true, rng = Random.default_rng())

Fits of transfer_multitask by coordinate descent along a decreasing sequence of levels λ₀, each warm started from the previous one, with fixed contrast weights. The automatic grid has nλ log-spaced levels from λmax ⋅ (1 + √eps) down to λminratio times that value; λminratio defaults to 0.01 when N < (L + 1)p and to 1e-4 otherwise.

λmax = max(maxⱼ ‖∑ₗ Sˡⱼ‖₂, maxₗ≥₁,ⱼ ‖Sˡⱼ‖₂ / aₗ), with scores Sˡⱼₖ = ⟨x̃ˡᵏⱼ, ỹˡᵏ⟩/(NK) at zero coefficients, is the smallest level at which every block, and not only Ŵᶠ, is zero. With screening, each level first solves the blocks kept by the sequential strong rule for each component (θ⁰ with weight one, θˡ with weight aₗ), then checks the optimality of every excluded block row and repairs (kkt_checks counts the rounds).

λpath, solutions (a TransferMultitaskSolution per level), kkt_checks, λmax, contrast_weights, nobs and domain_nobs are its fields. coef(path) and predict(path, X; domain) give a result per level; coef(path, λ) requires an exact path level.

CVTransferMultitaskPath(target, sources, options = CDOptions(); families = (:sqrt, :linear),
                        nfolds = 5, folds = nothing, λpath = nothing, nλ = 100, λminratio = nothing,
                        λ̃path = nothing, nλ̃ = 100, λ̃minratio = nothing, intercept = true,
                        screening = true, rng = Random.default_rng(), threads = true)

Tune TMTL(Fused) and then its correction by held-out target prediction, with coordinate descent and the options and screening of TransferMultitaskPath. families are contrast-weight families (:sqrt, :linear or a weight vector; one symbol is one family), and λpath gives levels common to all families. Folds are over the target rows, the same for every response, family, level and correction level; every source row stays in every training fit. Each fold refits on its target training rows and all sources, at the full-data grid of each family, recomputing sizes, N, the contrast weights of :sqrt and :linear, and the centering from that training data.

The held-out score of a fold is the target sum of squared errors divided by its held-out rows times K, using target intercepts from the training rows; cvm averages the folds weighted by their sizes, and cvse is its standard error. (family, λ₀) minimizes cvm among the candidates converged on the full data and in every fold; ties go to the earlier family, then the larger level. The correction level λ̃ is then selected with the same folds: each fold corrects its own fused fit at the selected (family, λ₀) on its training residuals, so the full-data fused fit is never used inside a fold. The correction grid λ̃path defaults to nλ̃ log-spaced levels from the full-data correction's λ̃max ⋅ (1 + √eps), the smallest level giving D̂ = 0, down to λ̃minratio times that value (default 0.01 when n₀ < p, else 1e-4), so a zero correction is a candidate. Ties go to the larger λ̃.

Grids come from the full data, and the folds are reused for the sequential selection, so cvm and cvm_correction are not unbiased estimates of the selected fit's error; independent test data are needed for that.

Fields: paths (the full-data path of each family), families, cvm and cvse (families × levels), selected (CartesianIndex(family, level)), family, λ₀, λ̃path, cvm_correction, cvse_correction, λ̃, fused (the full-data fused fit at the selection), solution (that fit with the full-data correction at λ̃) and folds. coef(cv) and predict(cv, X) use solution. With threads, folds run in parallel with generators seeded from rng, so the result does not depend on the number of threads.