Variable selection with error control

SLOPE, knockoff+, and complementary-pairs stability selection select columns of X. Their guarantees apply under different assumptions. SLOPE with BH weights controls FDR for orthogonal Gaussian regression with known noise level. Knockoff+ controls FDR under the exchangeability and statistic assumptions below. Stability selection bounds the expected number of selected variables with low base selection probabilities; interpreting that as a false-selection bound requires a condition on the noise variables.

SLOPE

For the loss normalized by the number of observations,

\[\min_{\beta_0,\beta}\frac{\|y-\beta_0-X\beta\|_2^2}{2n} +\sum_{j=1}^p\lambda_j|\beta|_{(j)},\qquad \lambda_1\geq\cdots\geq\lambda_p\geq0.\]

using HighDimensionalModels, Random, LinearAlgebra
rng = Xoshiro(82)
n, p = 240, 100
X = sqrt(n) * Matrix(qr(randn(rng, n, p)).Q)
y = X[:, 1:20] * fill(0.3, 20) + randn(rng, n)
sol = slope(X, y; q=0.1, sigma=1)
selected = findall(!iszero, sol.x)
path = SLOPEPath(X, y; nλ=40)
coef(path), predict(sol, X)

slope_bh(p; q, sigma, n) computes $\lambda_j=\sigma\Phi^{-1}(1-qj/(2p))/\sqrt n$. The stated FDR guarantee uses $X^TX/n=I$ and independent Gaussian noise of known standard deviation $\sigma$. A correlated design, an estimated noise level, or arbitrary path selection does not inherit this guarantee. An intercept centers both X and y; with an intercept, orthogonality must hold for the centered design. With standardize_x=true, the penalty is on the coefficients of columns with mean square one; the returned coefficients are on the original scale. SparseMatrixCSC inputs stay sparse: scaling touches only stored entries. Fully stored columns are centered before scaling to avoid cancellation for large offsets; other columns are centered implicitly in matrix products. The transformed design and solver workspaces use storage proportional to nnz(X) + n + p, in addition to the coefficient path. Dense inputs retain their explicit transformation.

The solver uses accelerated proximal gradient from ProximalSolvers and the sorted ℓ1 proximal operator from ProximalBase. Paths warm start and use the strong screening rule of Larsson, Bogdan and Wallin (2020). Every reduced fit is checked by the full proximal fixed-point equation, including its zero cluster; excluded variables are added and refitted when it fails. This check is necessary because sorted penalties do not separate across coordinates. screening=false fits the full problem at every level. The path's λpath contains scalar multipliers of path.lambda. Convergence, working sizes, repairs, iterations and full residuals are retained per level. Fits that miss the requested tolerance return converged=false with a warning.

Explicit calibration

BH remains the compatibility default. slope_gaussian(p; n, q, sigma, intercept=false) exposes a Gaussian-design adjustment; single fits and paths accept calibration=:gaussian. calibration=:user requires explicit lambda weights. Positional or keyword weights always take precedence over an automatic rule; q/sigma are then unused and recorded as nothing. Unknown selectors are rejected even with supplied weights. No noise estimator is introduced: automatic weights use the caller's sigma (default 1), without verifying that it is known. The calibration field records the actual rule, noise convention, dimensions, standardization and any flat tail. Existing coefficient accessors are unchanged.

sol = slope(X, y; calibration=:gaussian, sigma=1.0, intercept=true)
path = SLOPEPath(X, y; calibration=:gaussian)
sol.calibration.method
sol.calibration.flat_from
sol.calibration.flatten_reason

For Gaussian weights let $m=n-1$ with an intercept and $m=n$ otherwise, and let $b_i=\Phi^{-1}(1-qi/(2p))$. The dimensionless recurrence is

\[a_1=b_1,\qquad a_i=b_i\sqrt{1+\frac{\sum_{j<i}a_j^2}{m-i}}.\]

The implementation stops at the first nondecrease and repeats the preceding weight; it also stops before a nonpositive denominator. Final weights are $\sigma a_i/\sqrt n$ because the loss is divided by n even after centering. This explicitly defined flat tail supports p≥n without a pole. It requires m≥2 for p>1, or m≥1 for p=1. Sigma zero yields zero weights, but an automatic path still requires a nonzero base sequence; supply path factors explicitly when using zero weights. Weights are finite, nonnegative and nonincreasing.

The recurrence follows Bogdan et al. (2015), equations (26)–(27), with a first-nondecrease stopping rule and an explicit intercept dimension adjustment. It is motivated by independent Gaussian columns. Columns should have mean square one, either in the supplied data or through standardize_x; the penalty then acts on standardized coefficients and the returned coefficients remain on the input scale. Correlated designs, estimated sigma, and choosing a path level have no general finite-sample FDR guarantee from this adjustment. Here q specifies the sequence; it must not be read as a proven FDR bound for those settings. The orthogonal known-noise BH theorem stated above is separate.

SLOPE.jl 1.3.1 uses denominator max(1,n-i-1) with one-based i and applies a running minimum after the raw recurrence. R SLOPE 2.1.1 uses max(1,n-i) after converting its zero-based loop to one-based indexing. These conventions are intentionally compared separately from common-weight solver agreement; exact native-weight identity is not an acceptance condition.

Evidence and settings

The replicated study (archive pr115-gaussian-slope-20261026, report results/cluster-20261026/RESULTS.md) provides historical evidence measured on an earlier dependency stack; it was not rerun at this revision. The study compared 1,008 independent datasets per scenario, with separate estimated-noise and covariance stress cases. Its summary tables reported that, at known sigma=1/q=0.1 in the iid alternating-sign example, Gaussian adjustment reduced mean FDP from 0.1818 to 0.0733 and power from 0.9464 to 0.8570. Correlated and wide alternating-sign examples still had FDR above q=0.1 after adjustment. In the exactly orthogonal control it sacrificed power relative to BH. These historical findings support an explicit Gaussian option and retaining BH as the default; they do not extend the orthogonal BH theorem.

The study supported retaining tol=1e-7, max_iter=10000 and screening=true. Its independent checks passed every default probe on 84 fresh runtime inputs; screening improved the declared 20-level path comparisons. Stricter research tolerances could hit an iteration cap and did not establish better defaults. The study's coefficient cutoff 1e-6 defined its reported selections, not a new package selection guarantee. Its sigma=1 setting was a supplied-scale convention: the package does not establish that it is the true noise level or estimate it from std(y).

Knockoff+

X = randn(rng, n, p)
y = X[:, 1:20] * fill(0.3, 20) + randn(rng, n)
mx = knockoff_filter(X, y, 0.04; mu=zeros(p), Sigma=Matrix{Float64}(I,p,p), q=0.1, rng)
fx = knockoff_filter(X, y, 0.04; method=:fixed, q=0.1, rng)
mx.selected, mx.threshold, mx.W

gaussian_knockoffs uses the conditional Gaussian construction. For exact model-X control, the rows of X must follow the Gaussian distribution specified by the supplied population mean and covariance. Substituting estimated moments provides second-order matching without establishing finite-sample exchangeability. The filter deliberately requires explicit parameters.

fixed_knockoffs requires full column rank and $n\geq2p$; it preserves the Gram identities $\widetilde X^T\widetilde X=X^TX$ and $X^T\widetilde X=X^TX-\operatorname{diag}(s)$. Columns may have arbitrary positive norms. Both constructions use an equicorrelated solution computed on the correlation scale by default; no optimizer dependency is needed. A user-supplied feasible s is also accepted by the constructors and filter. Semidefinite conditional covariances use an eigenvalue square root, clipping only roundoff-sized negative values.

Diagonal construction choices

# Optional dependency: install Knockoffs separately, then load it explicitly.
using Knockoffs
plan = knockoff_diagonal(Sigma; strategy=:mvr)
K = gaussian_knockoffs(X, mu, Sigma; s=plan.s, noise=Z)
fit = knockoff_filter(X, y, 0.04; mu, Sigma, noise=Z, strategy=:maxent,
    construction_options=(niter=1000, tol=1e-10, λmin=1e-10),
    order=:cyclic, opt_tol=1e-9)
fit.method, fit.construction.strategy, fit.construction.diagnostics

method identifies the model (:gaussian, :fixed, :provided); strategy identifies the diagonal (:equicorrelated, :mvr, :maxent, :user). Omitting strategy selects equicorrelation, or user construction when s is supplied. Explicit :user requires s; the other strategies reject a supplied s. Provided knockoffs reject strategy, diagonal, innovations and construction options. These conflicts are checked before fitting. Legacy Gaussian parameters unused by the fixed/provided routes remain accepted for compatibility.

construction_options is a named tuple separate from lasso CDOptions. Equicorrelation accepts safety; user diagonals accept no optimizer options. Optimized construction accepts niter, tol, λmin, robust, verbose, and optimality_tol. MVR and maximum entropy require the optional Knockoffs.jl extension. It supplies only the established diagonal optimizer; Gaussian and fixed-X sampling stay native. Knockoffs.jl is MIT-licensed; the extension was validated with Knockoffs.jl 3.1.4. Only its diagonal optimizer is used. Its full sampler is not used because the pinned version's conditional sampler applies the noise with the wrong factor. Loading the extension does not change the equicorrelated default or existing seeded sampling arithmetic.

Optimized diagonals minimize reconstructability or maximize entropy on the correlation scale, then return Float32/Float64 covariance/Gram units. The actual rounded diagonal must pass the original, unperturbed objective's stationarity check and both native scale-relative feasibility checks. Backend iterate-change tolerance alone is insufficient. Failure raises an error without switching strategy or projecting an infeasible diagonal. The backend returns neither an actual iteration count nor a convergence status; diagnostics mark those fields unavailable and report independent stationarity separately. The optimized defaults are niter=1000, tol=1e-10, λmin=1e-10, robust=false, with optimality_tol=max(1e-6,100eps(input precision)). The recorded matrix study supports these controls; the backend's looser package defaults can return feasible diagonals that miss original-objective stationarity.

KnockoffDiagonal owns its diagonal and metadata. Passing s=plan.s explicitly records user construction; it validates that vector against the supplied covariance and does not reuse cached factors. KnockoffResult.construction records the actual strategy, owned diagonal and diagnostics. The existing seven-argument result constructors remain available with unspecified metadata.

An optimized diagonal does not repair a misspecified covariance, a statistic without swap antisymmetry, or a fixed-X procedure without sufficiency. The population-model assumptions and fixed-X restrictions above still apply.

The replicated construction study (archive pr117-knockoff-strategies-20261026, report results/cluster-20261026/RESULTS.md) records historical evidence. It was measured on an earlier dependency stack and was not rerun at this revision. The study uses 1,008 independent datasets per scenario and common innovations, lasso statistic and penalty across constructions. For the recorded AR(.4) design with alternating signals, power rises from 0.64% with equicorrelation to 39.90% with MVR, at estimated FDR 6.92%. Same-sign cases lose power, and optimized construction costs substantially more time: MVR full filters take 1.08–9.28 times the equicorrelation time on the tested cases. Equicorrelation, safety 1e-8, and knockoff+ target q=0.1 therefore remain the compatibility defaults. These results concern the listed simulation laws only; they do not supply a data-driven strategy or penalty selector, or a misspecified-model guarantee.

The statistic is $W_j=|\widehat\beta_j|-|\widehat\beta_{p+j}|$ from the joint lasso on $[X\ \widetilde X]$. The filter fits without an intercept or standardization after constructing knockoffs. Scaling must precede fixed-X construction. This fixed-X interface assumes a model without an intercept. For an intercept, project both X and y onto an orthonormal basis orthogonal to the constant vector, then construct knockoffs in those n−1 coordinates (requiring n−1 ≥ 2p). Subtracting column means alone does not ensure that the knockoffs are orthogonal to the removed intercept. A prespecified positive penalty, independent of y, preserves fixed-X sufficiency. Cross-validating this penalty does not have the same fixed-X guarantee. Model-X also requires a swap-antisymmetric statistic; nonunique lasso solutions and exact ties need care. A failure of the joint lasso to converge raises an error.

knockoff_threshold returns the smallest observed positive magnitude t with

\[\frac{1+\#\{j:W_j\leq-t\}}{\max(1,\#\{j:W_j\geq t\})}\leq q,\]

or Inf (no selection). Ties are included. knockoff_filter always uses this knockoff+ rule. The standalone threshold also supports offset=0, the original knockoff rule with a modified-FDR guarantee. Supplied knockoffs can be used with method=:provided, Xk=..., leaving their validity to the caller.

Complementary-pairs stability selection

st = stability_selection(X, y, 0.12; q=10, cutoff=0.8, npairs=50, rng)
st.selected, st.frequencies, st.pfer_bound

Each pair consists of two disjoint half samples. For odd n, one randomly chosen row per pair is omitted. On each half the base selector fits the lasso at the supplied penalty and keeps at most q nonzero coefficients, ranked by magnitude on the fitted scale. This is the fixed-penalty, capped-coefficient selector, not a lasso-path stopping rule. standardize_x=true and intercept=true are the defaults; these transformations are fitted separately on each half.

The frequency of a variable is its average selection indicator over all 2B halves; variables at or above cutoff are retained. pair_frequencies records simultaneous selection in both halves. The result also retains each selection and the pairs, so the calculation can be inspected. Explicit pairs can be shared with another implementation. Draws and per-half RNG seeds are made before fitting. julia --threads=N enables concurrent fits; parallel=false uses the same draws serially. Failed fits raise an error.

For $\tau=\mathtt{cutoff}>1/2$, Shah and Samworth's distribution-free bound is

\[E|\widehat S_{\rm CPSS}\cap L_\theta| \leq\frac{\theta q}{2\tau-1},\qquad L_\theta=\{j:P(j\in\widehat S_{\rm half})\leq\theta\}.\]

The deterministic cap makes q an upper bound on the expected base selection size. The reported pfer_bound = q^2/(p*(2cutoff-1)) uses $\theta=q/p$, as R's stabs does with assumption="none". This bounds expected false selections only when all noise variables belong to $L_{q/p}$; exchangeability and a no-worse-than-random base selector are sufficient conditions. It is not an unconditional FDR guarantee. stability_bound(...; theta=...) also exposes the general bound. The observed average selection size is descriptive and does not replace q in the guarantee. No unimodality or r-concavity improvement is assumed.

References and validation

• Bogdan et al. (2015), SLOPE.
• Su and Candès (2016), minimax adaptation.
• Larsson, Bogdan and Wallin (2020), strong screening.
• Barber and Candès (2015), fixed-X knockoffs.
• Candès et al. (2018), model-X knockoffs.
• Meinshausen and Bühlmann (2010), stability selection.
• Shah and Samworth (2013), complementary pairs.

The R golden-data generator is test/data/selection/golden.R. The replicated benchmark and profile harness and the Monte Carlo recipe are in the selection study (snapshot); the selection group of the benchmark harness measures the current fits. Monte Carlo evidence checks the listed finite simulation designs; it does not prove a guarantee.

SLOPEPath(X, y, factors=nothing; lambda=nothing, calibration=:bh, q=0.1, sigma=1,
          nλ=50, λmin_ratio=0.01, screening=true, intercept=false,
          standardize_x=false, max_iter=10000, tol=1e-7)

SLOPESolution

slope(X, y, lambda; intercept=false, standardize_x=false, max_iter=10000, tol=1e-7)
slope(X, y; lambda=nothing, calibration=:bh, q=0.1, sigma=1, kwargs...)

slope_bh(p; q=0.1, sigma=1, n=1)

slope_gaussian(p; n, q=0.1, sigma=1, intercept=false)

KnockoffDiagonal

knockoff_diagonal(Sigma; strategy=nothing, s=nothing, construction_options=(;))

KnockoffResult

equicorrelated_s(Sigma; safety=1e-8)

fixed_knockoffs(X; s=nothing, strategy=nothing, construction_options=(;))

gaussian_knockoffs(X, mu, Sigma; rng=Random.default_rng(), s=nothing, noise=nothing,
                  strategy=nothing, construction_options=(;))

knockoff_filter(X, y, lambda; method=:gaussian, mu=nothing, Sigma=nothing,
                Xk=nothing, q=0.1, rng=Random.default_rng(), options=CDOptions(),
                strategy=nothing, s=nothing, noise=nothing, construction_options=(;), kwargs...)

knockoff_threshold(W; q=0.1, offset=1)

lasso_coefdiff(X, Xk, y, lambda; options=CDOptions(), kwargs...)

StabilitySelectionResult

complementary_pairs(n, npairs; rng=Random.default_rng())

stability_bound(p, q, cutoff; theta=q/p)

stability_selection(X, y, lambda; q, cutoff=0.9, npairs=50, pairs=nothing,
                    parallel=true, intercept=true, standardize_x=true,
                    rng=Random.default_rng(), options=CDOptions(), kwargs...)