Fit SSRHE-Style Hessian L1 Regression — fit.ssrhe.hessian.l1.regression • geosmooth

Fits an $\ell_1$-adaptive SSRHE Hessian comparator using the row operator $A$ from ssrhe.hessian.operator: $$ \widehat f_\lambda = \arg\min_f \left\{ \frac{1}{2}\sum_i w_i(y_i-f_i)^2 + \lambda \|Af\|_1 \right\}. $$

Usage

fit.ssrhe.hessian.l1.regression(
  X,
  y,
  k = NULL,
  tangent.dim,
  lambda.grid = NULL,
  lambda.selection = c("cv", "fixed"),
  weights = NULL,
  n.lambda = 40L,
  nfolds = 5L,
  fold.id = NULL,
  loss = c("mse", "mae"),
  selection = c("min", "one.se"),
  nn.index = NULL,
  neighborhood.type = c("knn", "adaptive.radius", "supplied"),
  support.index = NULL,
  adaptive.k.scale = NULL,
  radius.rule = c("geomean", "max", "min"),
  radius.factor = 1.25,
  min.support = NULL,
  max.support = NULL,
  support.buffer = 2L,
  support.topup = c("nearest", "none"),
  tangent.dim.rule = c("fixed", "eigen.cumulative"),
  eigen.tolerance = 0.95,
  derivative.order = 2L,
  pinv.tol = sqrt(.Machine$double.eps),
  local.solver = c("auto", "normal.equations", "svd", "qr"),
  normal.equations.max.condition = 10000,
  solver = c("genlasso", "admm", "auto"),
  row.scaling = c("none", "l2"),
  admm.rho = NULL,
  admm.maxiter = 2000L,
  admm.abstol = 1e-06,
  admm.reltol = 1e-05,
  maxsteps = 2000L,
  minlam = 0,
  approx = FALSE,
  rtol = 1e-07,
  btol = 1e-07,
  eps = 1e-04,
  support.selection = c("rule", "cv"),
  support.grid = NULL,
  support.cv.max.candidates = 8L,
  return.local.diagnostics = FALSE,
  return.timing = FALSE,
  verbose = FALSE,
  admm.adaptive.rho = TRUE
)

Arguments

X: Numeric matrix with one observation per row. Distances for the default neighborhood search are Euclidean distances in these coordinates.
y: Numeric response vector of length nrow(X). Matrix responses are not yet supported for the $\ell_1$ path.
k: Integer number of nearest neighbors per point, including the point itself, used when neighborhood.type = "knn". This follows the SSRHE Matlab convention. For adaptive-radius neighborhoods, k may be omitted when adaptive.k.scale is supplied.
tangent.dim: Integer tangent dimension used for the local PCA chart. Required when tangent.dim.rule = "fixed". If NULL and tangent.dim.rule = "eigen.cumulative", the dimension is selected locally from the PCA variance ratio.
lambda.grid: Optional nonnegative lambda grid. For lambda.selection = "fixed", this must contain exactly one value. For lambda.selection = "cv", NULL builds a default grid from the fitted full-data generalized-lasso path.
lambda.selection: "cv" for observed-label K-fold cross-validation or "fixed" for a supplied fixed lambda.
weights: Optional nonnegative observation weights. Missing responses automatically receive zero weight.
n.lambda: Number of default lambda candidates when lambda.grid = NULL and lambda.selection = "cv".
nfolds: Number of validation folds over observed positive-weight labels. Ignored when fold.id is supplied.
fold.id: Optional integer fold assignments of length nrow(X).
loss: Validation loss, currently "mse" or "mae".
selection: Selection rule. "min" chooses the smallest mean validation loss. "one.se" chooses the largest lambda within one standard error of the minimum.
nn.index: Optional integer matrix of neighbor indices, with nrow(X) rows and k columns. Each row must contain its center vertex. If NULL, Euclidean k-NN including self is computed in C++.
neighborhood.type: Local support rule. "knn" uses the original rectangular self-including kNN neighborhoods. "adaptive.radius" builds variable-size supports from dgraphs::create.rknn.graph(). "supplied" uses support.index directly.
support.index: Optional list of integer vectors, one per row of X. Each element gives a variable-size local support and must contain its center vertex.
adaptive.k.scale: Integer local-scale k used by dgraphs::create.rknn.graph() when neighborhood.type = "adaptive.radius".
radius.rule, radius.factor: Adaptive-radius graph parameters passed to dgraphs::create.rknn.graph().
min.support: Optional minimum local support size for adaptive-radius supports. Defaults to the local quadratic design size plus support.buffer, clamped to nrow(X).
max.support: Optional maximum local support size. Oversized adaptive supports are trimmed to the closest max.support vertices while keeping the center vertex.
support.buffer: Nonnegative integer added to the local quadratic design size when choosing the default min.support.
support.topup: How undersized adaptive-radius supports are enlarged. "nearest" appends ambient nearest neighbors. "none" leaves supports unchanged and lets the C++ operator reject undersized supports.
tangent.dim.rule: Either "fixed" or "eigen.cumulative".
eigen.tolerance: Cumulative local PCA variance threshold used when tangent.dim.rule = "eigen.cumulative" and tangent.dim = NULL.
derivative.order: Integer derivative order of the local SSRHE-style operator. 2L is the original local Hessian energy. 3L is an experimental third-derivative energy whose rows estimate unique symmetric third-derivative tensor components with factorial and tensor-multiplicity scaling.
pinv.tol: Nonnegative tolerance multiplier for local pseudoinverses.
local.solver: Local least-squares backend used to map local function values to derivative coefficients. "auto" is the default: it uses normal equations when the local design is full rank and has condition number no larger than normal.equations.max.condition, and otherwise falls back to SVD. "normal.equations" requests the normal-equation solve for full-rank local designs and falls back only on hard numerical failures. "svd" is the most stable reference path. "qr" uses pivoted QR and falls back to SVD for rank-deficient local designs.
normal.equations.max.condition: Positive condition-number guard used by local.solver = "auto" before accepting the normal-equation backend. The default is deliberately conservative because normal-equation solves square the effective condition number and order-3 near-minimum local supports can be numerically fragile.
solver: Solver backend. "genlasso" uses the generalized-lasso path to propose a grid when needed, then refines every candidate with ADMM on the original objective. "admm" skips the path and requires an explicit grid. "auto" also proceeds when path construction fails, provided an explicit grid is supplied.
row.scaling: Optional row scaling for the penalty matrix before solving. "l2" scales nonzero rows of $A$ to unit Euclidean norm.
admm.rho, admm.maxiter, admm.abstol, admm.reltol: ADMM controls used by all backends, including refinement of generalized-lasso candidates. admm.rho = NULL chooses an initial algorithm penalty of 1 / max(1, mean(rowSums(D^2))), where D is the penalty operator after any requested row scaling. A positive numeric value sets the initial penalty explicitly. This is separate from statistical regularization lambda.grid. The iteration limit is a hard cap.
maxsteps, minlam, approx, rtol, btol, eps, verbose: Controls passed to genlasso.
support.selection: Support-profile selection rule. "rule" uses the supplied adaptive.k.scale, min.support, and max.support. "cv" is currently supported for neighborhood.type = "adaptive.radius" with lambda.selection = "cv" and chooses among rows of support.grid by an outer response cross-validation loop.
support.grid: Optional data frame of support profiles with columns adaptive.k.scale, min.support, and optional max.support. If NULL, ssrhe.support.grid builds a compact default grid.
support.cv.max.candidates: Maximum number of support profiles to try when support.selection = "cv".
return.local.diagnostics: Logical. If TRUE, compute additional R-side local chart diagnostics such as chart distortion and boundary asymmetry. These diagnostics are useful for operator audits, but can be skipped in fitting and cross-validation paths for speed.
return.timing: Logical. If TRUE, attach a phase-level elapsed time table to the returned operator. The timings separate R-side validation, neighborhood/support construction, native local operator construction, optional local diagnostics, sparse matrix assembly, and output finalization. For adaptive-radius neighborhoods, subphase timings are also attached in neighborhoods$timing.
admm.adaptive.rho: Logical; adapt the algorithm penalty during the first 1000 iterations by balancing tolerance-normalized primal and dual residuals. Set FALSE to hold the initial penalty fixed.

Value

A list of class "ssrhe.hessian.l1.fit" containing fitted values, residuals, selected lambda, the SSRHE operator, generalized-lasso path metadata, lambda-grid fitted values, and CV diagnostics when requested.

Details

This function exposes the SSRHE local polynomial derivative estimator as an $\ell_1$ penalty, rather than as the original SSRHE $\ell_2$ Hessian-energy penalty used by fit.ssrhe.hessian.regression. With derivative.order = 2L, $A$ contains local quadratic/Hessian rows. With derivative.order = 3L, $A$ contains local cubic third-derivative tensor rows with the same symmetric component scaling used by ssrhe.hessian.operator. The existing $\ell_2$ fit solves a linear system involving $B=A^\top A$. This function instead keeps the rows of $A$ and penalizes their absolute values. It is therefore closer in spirit to graph trend filtering, while retaining SSRHE's local-PCA derivative construction.

The default solver backend is genlasso. For larger or numerically fragile third-order operators, solver = "admm" fits the requested lambda values directly without computing a generalized-lasso path, and all backends use ADMM for the returned coefficients. solver = "auto" can additionally recover from path construction errors with an explicit grid. Cross-validation removes held-out labels from the data-fit term by setting their weights to zero, then scores predictions on those held-out labels.

ADMM uses proximal numerical stabilization, which does not add a ridge penalty to the requested objective. Convergence requires primal feasibility, the dual stopping residual, and stationarity of the original weighted L1 objective to meet their reported tolerances. A finite fixed-lambda iterate that reaches its cap is returned with a warning and solver$converged = FALSE; it is not a completed fit. Inspect solver$admm before using it, and consider a larger iteration limit or an adaptive penalty. Non-finite selected fits are errors.

CV candidates are eligible only when every fold has an acceptable solve. ADMM iterates that do not converge are not scored. The CV result retains fold.status, fold.converged, fold.diagnostics, n.successful.folds, and eligible. If no candidate qualifies, an error of class ssrhe_l1_cv_error retains fold diagnostics for inspection with tryCatch. A selected full-data ADMM fit that fails to converge is also an error in CV mode.

A generalized-lasso path alone can be inaccurate for nearly rank-deficient operators. It is retained in path for diagnostics, but its coefficients are not the returned fit. solver$backend = "genlasso_admm_refined" identifies this route; solver$path.objective records the unrefined objectives at the grid values. Backend warnings, including its numerical ridge for underdetermined designs, are retained in solver$warnings. CV uses the original objective without that additional ridge. The default ADMM tolerances are absolute 1e-6 and relative 1e-5; convergence describes numerical stopping at those tolerances, not predictive accuracy. With row.scaling = "l2", the reported objective uses the scaled penalty; energies$penalty.l1 is that penalty norm, while energies$hessian.l1 retains the unscaled operator norm.

With support.selection = "cv", each adaptive-radius support candidate builds a fresh SSRHE operator and runs the usual lambda CV with shared fold assignments. The selected fit is the candidate with the smallest selected CV error. This is useful when the default support-size rule is too rigid, but it can be much slower than tuning lambda for one fixed operator.

Choosing controls

First fit: specify X, y, k, tangent.dim, and lambda.grid. For a single supplied penalty set lambda.selection = "fixed" and solver = "admm" to avoid constructing a path.

Selection: fold.id, nfolds, loss, and selection control label CV. A grid requires one solve per penalty per fold; support CV also rebuilds geometry. Incomplete candidates are excluded, not silently scored.

Geometry: neighborhood, tangent and derivative controls define A. row.scaling changes the penalty itself and thus the statistical objective.

Numerics: admm.* controls stopping and algorithmic penalties for all backends. Path controls (maxsteps, rtol, btol, eps) affect the optional proposal path. They do not certify the final original-objective solve.

Diagnostics: solver$admm, candidate/fold status and objective components describe numerical acceptance. path and beta.grid can use substantial memory; use a compact explicit grid with ADMM when no full path is needed.

References

Kim, K. I., Steinke, F., and Hein, M. (2009). Semi-supervised regression using Hessian energy with an application to semi-supervised dimensionality reduction. Advances in Neural Information Processing Systems 22. https://papers.nips.cc/paper_files/paper/2009/hash/f4552671f8909587cf485ea990207f3b-Abstract.html

Examples

X <- matrix(seq(0, 1, length.out = 12), ncol = 1)
fit.ssrhe.hessian.l1.regression(
  X, sin(2 * pi * X[, 1]), k = 6L, tangent.dim = 1L,
  lambda.grid = 0.05, lambda.selection = "fixed", solver = "admm"
)
#> SSRHE Hessian L1 regression fit
#>   solver: admm 
#>   status: converged  
#>   responses: 1 
#>   lambda.selection: fixed 
#>   lambda: 0.05 
#>   objective: 1.523