Fit SSRHE-Style Hessian-Energy Regression

Fits the \(\ell_2\) Hessian-energy regularized estimator associated with ssrhe.hessian.operator. This is a direct SSRHE-style comparator for Hessian smoothing on point clouds; it is not a replacement for graph trend-filtering regression and is distinct from \(\ell_1\)-adaptive graph trend filtering.

Usage

fit.ssrhe.hessian.regression(
  X,
  y,
  k = NULL,
  tangent.dim,
  lambda1,
  lambda2 = 0,
  weights = NULL,
  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,
  stabilizer = lambda2 > 0,
  pinv.tol = sqrt(.Machine$double.eps),
  local.solver = c("auto", "normal.equations", "svd", "qr"),
  normal.equations.max.condition = 10000,
  ridge = 0,
  return.A = TRUE,
  return.local.diagnostics = FALSE,
  return.timing = FALSE,
  verbose = FALSE,
  lambda.selection = c("fixed", "cv", "gcv"),
  lambda1.grid = NULL,
  lambda2.grid = NULL,
  cv.control = list(),
  gcv.control = list()
)

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) or numeric matrix with nrow(X) rows. Matrix columns are fit as separate responses. NA values are allowed and are treated as unobserved by setting the corresponding observation weight to zero.

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.

lambda1

Nonnegative Hessian-energy penalty multiplier for \(f^\top B f = \|Af\|_2^2\).

lambda2

Nonnegative supplemental-stabilizer multiplier for \(f^\top B_S f\). If positive, stabilizer must be TRUE. Currently supported only for derivative.order = 2L.

weights

Optional nonnegative observation weights. May be NULL, a vector of length nrow(X), or a matrix with the same dimensions as y.

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.

stabilizer

Logical. If TRUE, also construct the supplemental stabilizer matrix described in the SSRHE supplement. Currently supported only for derivative.order = 2L.

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.

ridge

Nonnegative diagonal ridge added to the linear system for numerical stabilization.

return.A

Logical. If TRUE, return \(A\) as a sparse matrix.

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.

verbose

Logical. If TRUE, print a short native construction message.

lambda.selection

Penalty selection: "fixed" (the default) uses lambda1 and lambda2; "cv" holds out observed labels; "gcv" uses generalized cross-validation. Grid inputs require an explicit "cv" or "gcv" selection mode.

lambda1.grid, lambda2.grid

Nonnegative penalty candidates for CV or GCV. lambda1.grid is required in either selection mode; lambda2.grid = NULL means 0. Do not also supply scalar lambda1 or lambda2. Unless explicitly set, stabilizer is enabled when any supplemental penalty candidate is positive.

cv.control

Named list used only with lambda.selection = "cv". Entries: foldid = NULL (one fold label per row of X; nonpositive or missing labels are ignored), cv.folds = 5L, loss = "mse" (or "mae"), and selection = "min" (or "one.se", choosing the largest total penalty within one standard error of the minimum). Support-search entries are support.selection = "rule" (or "cv"), support.grid = NULL, and support.max.candidates = 8L. Fold count is ignored when foldid is supplied.

gcv.control

Named list used only with lambda.selection = "gcv". Entries: trace.method = "exact" (or "hutchinson"), trace.n.probes = 50L, trace.seed = NULL, support.selection = "rule" (or "gcv"), support.grid = NULL, and support.max.candidates = 8L. Unknown, unnamed, or duplicated control entries are errors. Controls for a different selection mode are rejected rather than ignored.

Value

A list of class "ssrhe.hessian.fit" containing fitted values, residuals, input response, weights, lambda parameters, objective/energy diagnostics, the reused operator, solver metadata, the call, and lambda.selection. CV adds class "ssrhe.hessian.cv.fit", cv.table, fold.id, and selection; GCV adds class "ssrhe.hessian.gcv.fit", gcv.table, gcv.trace, and selection. Support searches retain selected.support and the corresponding support.cv.table or support.gcv.table.

Details

For each response column, this function solves $$ (W + \lambda_1 B + \lambda_2 B_S + \epsilon I)\hat f = Wy, $$ where \(W\) is the diagonal matrix of observation weights and \(\epsilon\) is ridge. Fully observed data with lambda1 = lambda2 = ridge = 0 therefore reproduce the observed response exactly. Missing responses are excluded from the data-fit term by setting their weights to zero. The Matlab SSRHE semi-supervised convention is represented by a 0/1 labeled-indicator weights vector, or equivalently by setting unlabeled responses to NA: labeled vertices contribute unit diagonal data-fit terms and unlabeled vertices contribute no data-fit term.

The lower-level operator is matched to the Kim–Steinke–Hein SSRHE Matlab construction: self-including kNN neighborhoods, local PCA charts, a fixed-intercept local quadratic fit, doubled diagonal Hessian components, and optional supplemental stabilizer. The package exposes both \(A\) and \(B=A^\top A\); the fitted \(\ell_2\) estimator uses \(B\).

Penalty selection

Fixed fitting supports response vectors or matrices. CV and GCV currently select penalties for one response vector at a time. Label CV removes each validation fold from the data-fit term and scores its held-out observed positive-weight labels. It supports missing responses. GCV requires a fully observed response and strictly positive observation weights, and minimizes $$\frac{n^{-1}\sum_i w_i(y_i-\hat f_i)^2} {(1-\mathrm{tr}(S_\lambda)/n)^2},$$ where \(S_\lambda=(W+\lambda_1 B+\lambda_2 B_S+\epsilon I)^{-1}W\). The exact trace is deterministic. The Hutchinson option uses Rademacher probes, reports trace standard errors, and is reproducible with gcv.control = list(trace.method = "hutchinson", trace.seed = ...).

Both selection modes can also search adaptive-radius support profiles: set support.selection to the same mode in its control list. Each row of support.grid contains adaptive.k.scale, min.support, and optionally max.support. A NULL grid is constructed by ssrhe.support.grid. Each profile requires a new operator and a full penalty search, so this can substantially increase runtime.

Migration

The former fit.ssrhe.hessian.regression.cv() and fit.ssrhe.hessian.regression.gcv() entry points have been retired. Use this function with lambda.selection = "cv" or "gcv". Move fold.id to cv.control$foldid, nfolds to cv.control$cv.folds, and loss/selection to cv.control. Move gcv.trace.method, gcv.trace.n.probes, and gcv.trace.seed to trace.method, trace.n.probes, and trace.seed in gcv.control. Support search settings also belong in the active control list; both former support candidate limits are now named support.max.candidates.

Choosing controls

First fit: supply X, y, neighborhood size k, tangent.dim, and lambda1; use lambda.selection = "fixed" for a specified penalty.

Selection: lambda.selection, lambda1.grid, and lambda2.grid select penalties. Fold settings belong in cv.control; trace settings belong in gcv.control. Label CV and GCV are different criteria. Support selection rebuilds the operator for each candidate and multiplies this work.

Geometry: neighborhood, adaptive-radius, tangent and derivative controls determine the operator. Reuse it when fitting multiple responses.

Numerics: ridge and local conditioning controls stabilize solves. Exact GCV traces require more solves than stochastic trace estimates; the latter report sampling uncertainty and need a reproducible probe seed.

Diagnostics: inspect the selection table, energies, observed-label mask, and solver metadata. fitted() and residuals() preserve response matrices.

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.regression(
  X, sin(2 * pi * X[, 1]), k = 6L, tangent.dim = 1L, lambda1 = 0.1
)
#> SSRHE Hessian regression fit
#>   responses: 1 
#>   lambda1: 0.1 
#>   lambda2: 0 
#>   objective: 1.509 
fit.ssrhe.hessian.regression(
  X, sin(2 * pi * X[, 1]), k = 6L, tangent.dim = 1L,
  lambda.selection = "cv", lambda1.grid = c(0.01, 0.1),
  cv.control = list(cv.folds = 2L)
)
#> SSRHE Hessian regression CV fit
#>   responses: 1 
#>   selected lambda1: 0.1 
#>   selected lambda2: 0 
#>   CV mse : 0.2803 
fit.ssrhe.hessian.regression(
  X, sin(2 * pi * X[, 1]), k = 6L, tangent.dim = 1L,
  lambda.selection = "gcv", lambda1.grid = c(0.01, 0.1)
)
#> SSRHE Hessian regression GCV fit
#>   responses: 1 
#>   selected lambda1: 0.1 
#>   selected lambda2: 0 
#>   GCV: 0.3712 
#>   trace(S): 2.167