Skip to contents

Construct an SSRHE-Style Local Hessian Energy Operator

Source: R/ssrhe_hessian_energy.R
ssrhe.hessian.operator.Rd

Constructs the discrete local Hessian operator used by the semi-supervised regularization framework of Kim, Steinke, and Hein (2009). The constructor is intended as a package-facing, auditable R/C++ port of the reference behavior: it returns the sparse row operator \(A\), the quadratic energy matrix \(B=A^\top A\), and optionally the supplemental stabilizer \(B_S\).

Usage

ssrhe.hessian.operator(
  X,
  k = NULL,
  tangent.dim,
  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 = FALSE,
  pinv.tol = sqrt(.Machine$double.eps),
  local.solver = c("auto", "normal.equations", "svd", "qr"),
  normal.equations.max.condition = 10000,
  return.A = TRUE,
  return.B = TRUE,
  return.BS = stabilizer,
  return.sparse = TRUE,
  return.local.diagnostics = TRUE,
  return.timing = FALSE,
  verbose = FALSE
)

Arguments

X

Numeric matrix with one observation per row. Distances for the default neighborhood search are Euclidean distances in these coordinates.

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.

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.

return.A

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

return.B

Logical. If TRUE, return \(B=A^\top A\).

return.BS

Logical. If TRUE and stabilizer = TRUE, return the supplemental stabilizer \(B_S\).

return.sparse

Logical. If TRUE, attach Matrix sparse matrices in addition to raw triplets.

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.

Value

A list of class "ssrhe.hessian.operator" containing:

  • A: sparse local Hessian row operator, if requested;

  • B: sparse quadratic energy matrix \(A^\top A\), if requested;

  • BS: optional supplemental stabilizer matrix;

  • A.triplet, BS.triplet: raw sparse triplets;

  • nn.index: neighbor index matrix for rectangular kNN neighborhoods, if used;

  • support.index: list of actual local supports used by the backend;

  • neighborhoods: support-construction diagnostics;

  • row.table: one row per Hessian component row of \(A\);

  • diagnostics: local design rank, condition, and PCA summaries;

  • local.diagnostics: if requested, one row per center vertex with support size, design rank/condition, PCA variance, chart-distortion, boundary-asymmetry, and curvature-bias proxy diagnostics;

  • timing: if requested, a phase-level elapsed-time table;

  • native.timing: if requested, a C++ backend subphase timing table for native operator construction;

  • parity: notes documenting reference-behavior conventions.

Details

For each point \(x_i\), the operator constructs a local PCA chart from a self-including local support, projects the support into \(\mathbb R^m\), and centers the local coordinates at \(x_i\). A local quadratic model is then fitted with the intercept fixed at \(f_i\): $$ f(x_j) - f(x_i) \approx \sum_{a \le b} h_{ab} z_{ja} z_{jb} + \sum_a g_a z_{ja}. $$ The local linear map from function values to the quadratic coefficients is the Matlab-compatible $$ \mathrm{RegMat} = X^+ - X^+ \mathrm{IndMat}, $$ where \(X\) is the reduced quadratic-plus-linear design and \(\mathrm{IndMat}\) copies the center value into all local rows. Diagonal Hessian components are scaled by \(\sqrt{2}\) and off-diagonal components by \(1\), so that the row energy matches the reference Hessian energy. With derivative.order = 3L, the reduced local design uses cubic monomials first, followed by quadratic and linear monomials. The returned rows estimate unique symmetric third-derivative components \(\partial_{abc} f\) with scale $$ s_{abc}=m_{abc}\sqrt{\mu_{abc}}, $$ where \(m_{abc}\) is the factorial monomial-to-derivative multiplier (\(6\) for \(aaa\), \(2\) for \(aab\), and \(1\) for three distinct indices) and \(\mu_{abc}\) is the ordered-tensor multiplicity (\(1\), \(3\), or \(6\)). Thus \(\|Af\|_2^2\) approximates the squared Frobenius norm of the full symmetric third-derivative tensor.

The returned \(A\) is the stacked sparse matrix of these local derivative coefficient rows. The \(\ell_2\) SSRHE penalty is $$ f^\top B f = \|Af\|_2^2,\quad B=A^\top A. $$ Exposing \(A\) is useful for auditability and for future \(\ell_1\)-style variants based on \(\|Af\|_1\).

Adaptive-radius graph construction uses dgraphs; fixed-k and supplied neighborhoods are package-local geosmooth paths.

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)
ssrhe.hessian.operator(X, k = 6L, tangent.dim = 1L)
#> SSRHE Hessian operator
#>   vertices: 12 
#>   operator rows: 12 
#>   neighborhood.type: knn 
#>   k: 6 
#>   tangent.dim.rule: fixed 
#>   derivative.order: 2 
#>   stabilizer: FALSE