Package 'geosmooth'

Description

Geometric smoothing and conditional expectation methods for coordinate data, point clouds, and weighted graphs. The package provides local polynomial, model-averaged, trend-filtering, graph low-pass, occupation-density, and Hessian-energy methods.

Start here

Use browseVignettes("geosmooth") to open the installed HTML guides, including vignette("function-guide", package = "geosmooth"). The online guide links to the full reference. The task links below work in installed help without an internet connection.

Finding a function

Read vignette("function-guide", package = "geosmooth") for a task-based guide to all exported functions, prediction limits, and examples.

Local regression

Start with fit.lps or fit.malps.

Lifting trend filters

Use fit.lpl.tf or fit.slpl.tf for residual penalties with optional prediction synchronization.

Graph smoothing

Use fit.metric.graph.lowpass for spectral smoothing or fit.graph.trend.filtering for graph-difference penalties.

Hessian penalties

Use fit.ssrhe.hessian.regression for quadratic energy or fit.ssrhe.hessian.l1.regression for an absolute-value penalty.

Probability mass

Use fit.density for support weights, fit.subject.od for repeated visits, and normalize.density to convert an existing fitted field.

Synthetic examples

Combine components with synthetic.spec and draw data with materialize.synthetic.

Author(s)

Maintainer: Pawel Gajer [email protected]

Authors:

• Pawel Gajer [email protected]

Other contributors:

• Gael Guennebaud (Eigen) [contributor]

• Benoit Jacob (Eigen) [contributor]

• Authors of Eigen (Authorship and copyright in the included Eigen library) [copyright holder]

• Yixuan Qiu (Spectra) [contributor, copyright holder]

• Contributors to Spectra (Copyright in the included Spectra library) [copyright holder]

• Sunil Arya (ANN) [contributor, copyright holder]

• David M. Mount (ANN) [contributor, copyright holder]

• University of Maryland (ANN) [copyright holder]

See Also

vignette("function-guide", package = "geosmooth"), vignette("synthetic-datasets", package = "geosmooth")

Examples

X <- matrix(seq(0, 1, length.out = 30), ncol = 1)
y <- sin(2 * pi * X[, 1])
fit <- fit.lps(X, y, support.grid = 8L, degree.grid = 1L,
               kernel.grid = "tricube", foldid = rep(1:3, 10))
head(predict(fit))

Description

Applies every requested low-pass parameter to one response or a matrix of responses using a reusable spectral basis.

Usage

apply.metric.graph.lowpass.path(
  basis,
  y,
  eta.grid,
  filter.type = c("heat_kernel", "tikhonov", "cubic_spline", "gaussian", "exponential",
    "butterworth"),
  block.size = NULL,
  truncation.tol = 1e-04,
  unresolved.action = c("warn", "error", "allow"),
  exact.zero = TRUE
)

Arguments

Value

A list of class "metric.graph.lowpass.path". For one response, fitted.values is an $N$ by $J$ matrix. For multiple responses it is an $N$ by $J$ by $S$ array.

Examples

adj <- list(2L, c(1L, 3L), c(2L, 4L), 3L)
lengths <- list(1, c(1, 1), c(1, 1), 1)
basis <- metric.graph.lowpass.basis(adj, lengths, n.eigenpairs = 4L,
                                    eigen.solver = "dense")
apply.metric.graph.lowpass.path(basis, 1:4, eta.grid = c(0, 0.1, 1))

Description

Summarize a synthetic dataset as one metadata row

Usage

## S3 method for class 'synthetic_dataset'
as.data.frame(x, row.names = NULL, optional = FALSE, ...)

Arguments

Value

A one-row data frame.

Examples

x <- materialize.synthetic(synthetic.registry.spec("G1"), n = 20L, seed = 1L)
as.data.frame(x)

Description

Reuses the supports, local charts, prediction supports, averaging weights, selected degree, and local-solver controls from an existing fit.malps object, then repeatedly calls refit with bootstrap-style case weights. This is a fixed-selection uncertainty diagnostic: it measures refit variability conditional on the fitted MALPS support profile and does not rerun cross-validation or rebuild supports in each replicate.

Usage

bootstrap.malps(
  object,
  B = 200L,
  weight.type = c("bayesian", "multinomial"),
  y = NULL,
  probs = NULL,
  conf.level = 0.95,
  seed = NULL,
  max.failures = max(10L, B),
  keep.weights = FALSE,
  verbose = FALSE,
  ...
)

Arguments

Details

Two weight generators are available. "bayesian" draws positive exponential weights and rescales them to have mean one. This is the default because it preserves positive local support weights and therefore avoids many local-support degeneracies. "multinomial" draws ordinary nonparametric bootstrap counts with total count nrow(object$X); this can fail when a replicate assigns zero weight to all observations in a local support, and such failures are recorded.

Value

A list of class "malps_bootstrap" containing:

• replicate.fitted.values: an $n \times B$ matrix of successful replicate fitted values;

• replicate.weights: an optional $n \times B$ matrix of bootstrap weights when keep.weights = TRUE;

• failures: a data frame with columns attempt, failure, and message;

• summary: pointwise fitted values, bootstrap summaries, and empirical lower/upper intervals;

• bootstrap metadata such as B.completed, n.failures, conf.level, and weight.type.

Examples

X <- matrix(seq(0, 1, length.out = 20), ncol = 1)
fit <- fit.malps(X, X[, 1]^2, degree = 1L,
                 support.type = "knn", support.size = 8L)
bootstrap.malps(fit, B = 5L, seed = 1L)

Description

Runs one candidate graph per process, with dynamically scheduled jobs and no nested fold parallelism. On Unix, forked workers inherit the loaded package sources. Run from Rscript, not a GUI or multithreaded parent process.

Usage

compare.masked.graph.fields(operators, studies, workers = 1L, ...)

Arguments

Value

Named graph results, each containing named study results or an explicit error record. Retains per-observation predictions for common-support scoring.

Description

Compare two synthetic datasets scientifically

Usage

compare.synthetic.dataset(x, y, tolerance = c(1e-12, 1e-10))

Arguments

Value

A list with equality and mismatch details.

Examples

spec <- synthetic.registry.spec("G1")
x <- materialize.synthetic(spec, n = 20L, seed = 1L)
y <- materialize.synthetic(spec, n = 20L, seed = 1L)
compare.synthetic.dataset(x, y)

Description

Select smoothing within outer training folds, predict held-out groups, and finally tune and refit on all observations. Graph construction is label-free and fixed (transductive prediction). Unsupported predictions remain NA.

Usage

cv.masked.graph.field(
  operator,
  vertex,
  y,
  group,
  lambdas = 10^seq(-6, 2, length.out = 17),
  outer.folds = 5L,
  inner.folds = 4L,
  seed = 1L,
  fold.id = NULL,
  progress = NULL
)

Arguments

Value

Outer predictions, baseline predictions, scores on supported cases, coverage, tuning records, fold assignments and final full-data field. Scores do not account for selecting a graph after comparing graphs.

Description

Fits a local chart-kernel field by evaluating a Nadaraya–Watson style smoother at each evaluation point. The function is a general fitted-field model, not an occupation-density wrapper: density workflows should call normalize.density on the returned fit.

Usage

fit.chart.kernel(
  X,
  y,
  X.eval = NULL,
  support.size = min(15L, nrow(X)),
  kernel = c("gaussian", "tricube", "epanechnikov", "triangular"),
  bandwidth.multiplier = 1,
  support.grid = NULL,
  kernel.grid = NULL,
  bandwidth.multiplier.grid = NULL,
  foldid = NULL,
  cv.folds = 5L,
  cv.seed = 1L,
  coordinate.method = c("coordinates", "local.pca"),
  chart.dim = NULL,
  chart.dim.grid = NULL,
  selection.strategy = c("grid", "sparse_kd", "plateau_kd"),
  chart.dim.max = NULL,
  geometry.margin = 0L,
  auto.chart.support.metric = c("coordinates", "operator", "both"),
  auto.chart.selection.metric = c("coordinates", "operator"),
  quadrature.weights = NULL,
  denominator.floor = sqrt(.Machine$double.eps),
  return.details = TRUE
)

Arguments

Details

For an evaluation point $x_u$, the method chooses a local support $U_u$, builds either centered ambient coordinates or a local PCA chart, and computes

$\widehat f(x_u) = \frac{\sum_{i\in U_u} y_i K_h(z_{ui})} {\sum_{i\in U_u} q_i K_h(z_{ui})}.$

Here $z_{ui}$ is the local chart coordinate of $x_i-x_u$, $q_i$ is an optional quadrature weight, and $K_h$ is the selected kernel.

Value

A list with class "chart_kernel" containing fitted.values, source/evaluation supports, selected controls, and denominator diagnostics.

Examples

X <- matrix(seq(0, 1, length.out = 21), ncol = 1)
y <- exp(-20 * (X[, 1] - 0.5)^2)
fit.chart.kernel(X, y, support.size = 7L, kernel = "gaussian")

Description

Dispatches to a dedicated density estimator. Subject-occupation density estimation is one application: construct a sparse mass vector over a fixed support set and call this generic density layer.

Usage

fit.density(
  X,
  weights = NULL,
  method = c("empirical", "graph_random_walk"),
  graph = NULL,
  graph.control = list(),
  density.control = list(),
  return.details = TRUE,
  ...
)

Arguments

Value

A list of class "density_fit" with fields method.id, status, rho, empirical.rho, fitted.raw, theta, accounting, smoothness, timing, diagnostics, and warnings.

Examples

X <- matrix(seq(0, 1, length.out = 6), ncol = 1)
fit.density(X, weights = c(1, 0, 2, 0, 0, 1), method = "empirical")

Description

Fits graph trend filtering on a supplied undirected graph. Supported orders are 0L, 1L, and 2L:

$\widehat\beta_\lambda = \arg\min_{\beta\in\mathbb R^n} \left\{ \frac{1}{2}\|y-\beta\|_2^2 + \lambda\|\Delta_w^{(k+1)}\beta\|_1 \right\}.$

Usage

fit.graph.trend.filtering(
  adj.list,
  weight.list = NULL,
  y,
  order = 0L,
  lambda.grid = NULL,
  lambda.selection = c("cv", "fixed"),
  weight.rule = c("conductance", "sqrt.conductance", "unit"),
  operator.family = c("graph.laplacian.recursive", "path.divided.difference"),
  path.family = c("branch.continuation", "all.simple"),
  path.weighting = c("unit", "anchor.mean"),
  n.lambda = 40L,
  nfolds = 5L,
  foldid = NULL,
  maxsteps = 2000L,
  minlam = 0,
  approx = FALSE,
  rtol = 1e-07,
  btol = 1e-07,
  eps = 1e-04,
  verbose = FALSE
)

Arguments

Details

Graph trend filtering is intended as an $\ell_1$-adaptive comparator to the $\ell_2$ graph smoothers now housed in geosmooth and the legacy rdgraph comparator. The low-pass smoother penalizes broad quadratic graph roughness through a term like

$\eta f^\top L f,$

which shrinks high-frequency variation everywhere. Graph trend filtering instead penalizes absolute graph differences. For order = 0L, this is

$\lambda \sum_e \omega_e |\beta_j-\beta_i|,$

so many edge differences can become exactly zero while selected edges carry jumps. In this sense order = 0L is locally adaptive: fitted values can be piecewise constant over graph regions. Higher orders use

$\Delta_w^{(1)} = D_w,\qquad \Delta_w^{(2)} = L_w,\qquad \Delta_w^{(3)} = D_w L_w,$

with $L_w = D_w^\top D_w$.

The supplied weight.list, when used, is interpreted as a graph trend-filtering edge weight or conductance for operator.family = "graph.laplacian.recursive". For operator.family = "path.divided.difference", weight.list is interpreted as metric edge length and is used to form local path coordinates. Use weight.rule = "sqrt.conductance" when the recursive graph square-root convention is desired, and weight.rule = "unit" for an unweighted recursive graph fused-lasso penalty.

The solver backend is currently genlasso. Cross-validation refits the generalized-lasso problem on vertex-held-out training sets using a selection matrix $X$; this is useful for small and moderate diagnostic problems but is not yet a scalable production graph-trend-filtering solver.

Value

A list of class "graph.trend.filtering.fit" containing fitted values, residuals, selected lambda, the operator, path metadata, and CV diagnostics when requested.

References

Wang, Y.-X., Sharpnack, J., Smola, A., and Tibshirani, R. J. (2016). Trend filtering on graphs. Journal of Machine Learning Research, 17(105), 1–41. https://www.jmlr.org/papers/v17/15-147.html

Examples

if (requireNamespace("genlasso", quietly = TRUE)) {
  adj <- list(2L, c(1L, 3L), c(2L, 4L), 3L)
  fit.graph.trend.filtering(
    adj, y = c(0, 0.2, 1.8, 2), lambda.grid = 0.3,
    lambda.selection = "fixed", weight.rule = "unit"
  )
}

Description

Fits a local likelihood model at each evaluation point and reads off one raw fitted value at that evaluation anchor. Density and Bernoulli workflows should convert the returned fitted field with normalize.density when the target is a subject-occupation density.

Usage

fit.local.likelihood(
  X,
  y,
  X.eval = NULL,
  likelihood.family = c("density", "bernoulli"),
  support.size = min(15L, nrow(X)),
  degree = 1L,
  kernel = c("gaussian", "tricube", "epanechnikov", "triangular"),
  bandwidth.multiplier = 1,
  support.grid = NULL,
  degree.grid = NULL,
  kernel.grid = NULL,
  bandwidth.multiplier.grid = NULL,
  lambda.ridge.grid = NULL,
  foldid = NULL,
  cv.folds = 5L,
  cv.seed = 1L,
  coordinate.method = c("coordinates", "local.pca"),
  chart.dim = NULL,
  chart.dim.grid = NULL,
  selection.strategy = c("grid", "sparse_kd", "plateau_kd"),
  chart.dim.max = NULL,
  design.margin = 2L,
  auto.chart.support.metric = c("coordinates", "operator", "both"),
  auto.chart.selection.metric = c("coordinates", "operator"),
  quadrature.weights = NULL,
  lambda.ridge = 1e-08,
  min.local.mass = sqrt(.Machine$double.eps),
  min.nonzero.mass = 1L,
  fallback = c("degree0", "zero", "chart_kernel", "na"),
  optimizer = c("newton", "optim"),
  max.iter = 50L,
  tol = 1e-08,
  return.details = TRUE
)

Arguments

Details

The likelihood.family = "density" branch uses a local exponential tilt of the chart reference measure. The "bernoulli" branch uses a weighted local logistic likelihood and returns fitted probabilities at the evaluation anchors.

Value

A list with class "local_likelihood" containing fitted.values, selected controls, and local solver diagnostics.

Examples

X <- matrix(seq(0, 1, length.out = 21), ncol = 1)
y <- exp(-20 * (X[, 1] - 0.5)^2)
y <- y / sum(y)
fit.local.likelihood(X, y, likelihood.family = "density",
                     support.size = 7L, degree = 0L)

Description

Fits an LPL-TF estimator with a fixed operator:

$\hat f = \arg\min_f \frac12\|y-f\|_2^2 + \lambda\|A_{\mathrm{LPL}}f\|_1.$

Usage

fit.lpl.tf(
  X = NULL,
  y,
  operator = NULL,
  adj.list = NULL,
  weight.list = NULL,
  graph = NULL,
  degree = 2L,
  lambda.grid = NULL,
  lambda = NULL,
  lambda.selection = c("cv", "fixed"),
  operator.grid = NULL,
  foldid = NULL,
  cv.folds = 5L,
  cv.loss = c("mse", "rmse", "mae"),
  cv.seed = NULL,
  cv.repeats = 1L,
  solver = c("genlasso"),
  selection = c("min", "one.se"),
  n.lambda = 80L,
  maxsteps = 2000L,
  minlam = 0,
  approx = FALSE,
  rtol = 1e-07,
  btol = 1e-07,
  eps = 1e-04,
  verbose = FALSE,
  ...
)

Arguments

Value

A list of class "lpl_tf".

Examples

if (requireNamespace("genlasso", quietly = TRUE)) {
  X <- matrix(seq(0, 1, length.out = 18), ncol = 1)
  fit.lpl.tf(X, sin(2 * pi * X[, 1]), degree = 1L,
             support.type = "knn", support.size = 7L,
             lambda = 0.1, lambda.selection = "fixed")
}

Description

Fits a local polynomial smoother (LPS) and selects its support size, polynomial degree, and kernel by cross-validation. By default, the smoother works in the observed ambient coordinates: each prediction point uses its nearest training points in Euclidean distance, centers the support at the prediction point, fits a weighted local polynomial, and uses the fitted intercept as the prediction.

Usage

fit.lps(
  X,
  y,
  foldid = NULL,
  support.grid = c(10L, 15L, 20L),
  degree.grid = 0:2,
  kernel.grid = c("gaussian", "tricube"),
  cv.folds = 5L,
  cv.seed = 1L,
  X.eval = NULL,
  coordinate.method = c("coordinates", "local.pca"),
  chart.dim = NULL,
  chart.dim.grid = NULL,
  local.chart.method = c("pca", "second.order.svd"),
  auto.chart.support.metric = c("coordinates", "operator", "both"),
  auto.chart.selection.metric = c("coordinates", "operator"),
  backend = c("auto", "R", "cpp", "cpp.local.pca"),
  design.basis = c("orthogonal.polynomial.drop", "monomial", "weighted.qr",
    "weighted.qr.drop"),
  design.drop.tol = 1e-08,
  ridge.multiplier.grid = c(0, 1e-10, 1e-08),
  ridge.condition.max = 1e+12,
  unstable.action = c("na", "mean"),
  outcome.family = c("gaussian", "bernoulli", "binomial"),
  bandwidth.multiplier.grid = 1,
  keep.cv.predictions = FALSE,
  ridge.shrinkage.target = c("zero", "local.mean"),
  selection.strategy = c("grid", "sparse_kd", "plateau_kd"),
  chart.activation = c("none", "subject.od"),
  chart.activation.response = NULL,
  chart.activation.control = list(),
  chart.dim.max = NULL,
  design.margin = 2L
)

Arguments

Details

The optional coordinate.method = "local.pca" mode keeps the same support and kernel weighting rule, but builds the local polynomial in a local PCA chart centered at each prediction point. With chart.dim = "auto", the chart dimension is estimated as one global scalar from observed X only, using the same shared local-PCA dimension helper used by LPL-TF and S-LPL-TF. The experimental chart.dim = "local.auto" mode estimates an input-only local chart dimension separately for each prediction anchor.

Value

A list of class "lps" with response-scale fitted.values, unmodified local least-squares fitted.values.raw, selected parameters, a candidate CV table, the requested local.chart.method, and the effective chart method used for reporting. In outcome.family = "bernoulli" or "binomial" mode, fitted.values are response-scale probabilities in [0,1], fitted.values.raw are the un-clipped conditional-expectation estimates for "bernoulli" and the fitted probabilities for "binomial", cv.table$cv.brier.observed is the observed CV Brier score of the response-scale (clipped) predictions with Inf for candidates having any non-finite prediction, and probability.diagnostics records raw/clipped probability ranges, out-of-range fractions, and Brier/log-loss diagnostics. The Brier and log-loss diagnostics are defined only when the fitted predictions have the same length as the training response, which is the default X.eval = X path. In "binomial" mode, logistic.diagnostics records local logistic solve attempts, convergence statuses, fallback-path counts, event-rate fallback counts, and NA failure counts separately for CV and final fitting.

Choosing controls

First fit: supply X, y, foldid, and small support.grid, degree.grid, and kernel.grid candidates. See the installed function guide.

Selection: cv.folds, cv.seed, and selection.strategy control the search. Candidate combinations and folds multiply the fitting work; use lps.nested.cv() for an outer assessment of the tuning procedure.

Geometry: coordinate.method, chart.dim, and local chart controls describe neighborhoods and coordinate systems, not a different response model.

Numerics: backend, design.basis, rank thresholds and ridge controls govern stable local solves. Inspect lps.backend.diagnostics() before changing them. A numerical ridge is distinct from neighborhood selection.

Diagnostics: selected, cv.table, and diagnostics retain selection and local-fit information. A dense smoother.matrix needs one entry per evaluation/training pair; it is restricted to supported configurations.

References

Fan, J. and Gijbels, I. (1996). Local Polynomial Modelling and Its Applications. Chapman and Hall/CRC. ISBN 9780412983214.

Examples

X <- matrix(seq(0, 1, length.out = 20), ncol = 1)
y <- sin(2 * pi * X[, 1])
fit <- fit.lps(
  X, y, support.grid = 8L, degree.grid = 1L,
  kernel.grid = "tricube", cv.folds = 2L, backend = "R"
)
head(fit$fitted.values)

Description

Fits the current model-averaged local polynomial smoother (MALPS) on a supplied coordinate matrix. MALPS fits local polynomial regressions around observed anchor points and averages all positive-weight local predictions at each training point. The current implementation supports supplied coordinates, observed or supplied coordinate anchors, fixed support parameters, deterministic cross-validation and exact dense GCV over support parameters, local PCA chart coordinates, ordinary and weighted new-response refits, coordinate-space new-point prediction, linear-smoother diagnostics, and fixed-selection bootstrap uncertainty summaries. Graph-geodesic supports are available for training and refit; graph-geodesic new-point prediction is deferred.

Usage

fit.malps(
  X,
  y,
  graph = NULL,
  adj.list = NULL,
  weight.list = NULL,
  graph.stage = "final",
  anchor.index = NULL,
  anchor.coordinates = NULL,
  degree = 2L,
  degree.grid = NULL,
  support.type = c("adaptive.radius", "knn", "fixed.radius"),
  support.size = NULL,
  support.grid = NULL,
  radius = NULL,
  radius.grid = NULL,
  min.support = NULL,
  min.support.grid = NULL,
  support.buffer = 3L,
  kernel = c("epanechnikov", "triangular", "gaussian", "tricube"),
  kernel.grid = NULL,
  model.weight.rule = c("none", "condition", "support", "boundary", "quality"),
  duplicate.action = c("keep", "error"),
  coordinate.method = c("coordinates", "local.pca"),
  chart.dim = NULL,
  support.metric = c("auto", "coordinates", "graph.geodesic"),
  auto.chart.support.metric = c("coordinates", "operator", "both"),
  auto.chart.selection.metric = c("coordinates", "operator"),
  support.selection = c("fixed", "cv", "gcv"),
  foldid = NULL,
  cv.folds = 5L,
  cv.loss = c("rmse", "mae", "mse"),
  cv.repeats = 1L,
  cv.seed = NULL,
  cv.one.se = FALSE,
  gcv.exact.max.n = 1000L,
  local.solver = c("auto", "normal.equations", "qr", "svd"),
  normal.equations.max.condition = 1e+08,
  robust.iterations = 0L,
  robust.tuning.constant = 6,
  verbose = FALSE,
  ...
)

Arguments

Details

Current development benchmarks suggest using the plain MALPS defaults for production-style fits unless a specific diagnostic motivates a non-default option: support.selection = "cv", robust.iterations = 0L, model.weight.rule = "none", and observed anchors (anchor.coordinates = NULL). Exact GCV selection (support.selection = "gcv") is a useful faster option for fixed-weight, non-robust MALPS fits, especially exploratory continuous smoothing runs, but it is not yet the universal default. Binary responses are treated as numeric conditional-expectation targets; fitted values are not clipped inside the smoother. Grid anchors, empirical-Bayes shrinkage, robust local residual reweighting, and conservative/smoothed GCV selection remain experimental controls rather than default recommendations.

Graph construction is supplied by dgraphs or by explicit adj.list/weight.list payloads. Coordinate MALPS paths and graph-geodesic payload validation are package-local; shortest-path distances are computed through dgraphs.

Value

A list of class "malps" with fitted values, local model coefficients, supports, averaging weights, diagnostics, and selection metadata.

Choosing controls

First fit: supply X, y, degree, and support.type with a suitable support.size or radius. Coordinate-distance fits support new predictions.

Selection: support and degree grids, foldid, and cv.folds govern held-out selection. Robust refits recompute residual weights; selection is not automatically repeated by refit.

Geometry: support.metric, graph inputs, anchors, and chart controls determine the local models. Graph-geodesic fits require aligned vertices and do not currently predict at new targets.

Numerics: local.solver, conditioning controls, and robust weighting have distinct meanings. Keep their defaults until diagnostics identify a need.

Diagnostics: inspect selected, model weights, and stored local fits. smoother.matrix allocates a dense matrix and applies only to fixed weighting; its max.n guard makes that cost explicit.

Examples

X <- matrix(seq(0, 1, length.out = 20), ncol = 1)
fit <- fit.malps(X, sin(2 * pi * X[, 1]), degree = 1L,
                 support.type = "knn", support.size = 8L)
head(fit$fitted.values)

Description

Estimate a conditional mean by weighted least squares with an unnormalized graph Laplacian penalty. Repeated observations at a vertex remain separate. Unobserved components return NA rather than an arbitrary constant.

Usage

fit.masked.graph.field(operator, vertex, y, lambda, weights = NULL)

Arguments

Value

A list with vertex fitted values, support flags, observation weights, objective, and relative linear-system residual.

Description

Fits graph-spectral low-pass regression on a supplied graph by transforming metric edge lengths into conductances and smoothing the response in the eigenbasis of the weighted graph Laplacian.

Usage

fit.metric.graph.lowpass(
  adj.list,
  weight.list,
  y,
  conductance.rule = c("inverse.length.power", "exp.length", "exp.length.squared",
    "self.tuned.gaussian"),
  conductance.epsilon = 1e-08,
  conductance.alpha = 1,
  conductance.sigma = NULL,
  conductance.sigma.rule = c("edge.quantile", "median", "local.k"),
  conductance.sigma.quantile = 0.75,
  conductance.local.k = 5L,
  laplacian.type = c("unnormalized", "symmetric.normalized"),
  n.eigenpairs = 50L,
  filter.type = c("heat_kernel", "tikhonov", "cubic_spline", "gaussian", "exponential",
    "butterworth"),
  eta.grid = NULL,
  n.candidates = 40L,
  eigen.solver = c("auto", "sparse", "dense"),
  dense.eigen.threshold = 200L,
  dense.fallback.threshold = 5000L,
  dense.fallback = c("auto", "never", "always"),
  verbose = FALSE,
  eta.search = c("fixed", "guarded.gcv"),
  eta.expansion.factor = 3,
  eta.max.expansions = 3L,
  eta.identity.departure = 0.01,
  eta.truncation.tol = 1e-04
)

Arguments

Value

A list of class "metric.graph.lowpass.fit".

References

Gajer, P. and Ravel, J. (2025). Adaptive geometric regression for high-dimensional structured data. arXiv:2511.03817. doi:10.48550/arXiv.2511.03817

Examples

adj <- list(2L, c(1L, 3L), c(2L, 4L), 3L)
lengths <- list(1, c(1, 1), c(1, 1), 1)
fit.metric.graph.lowpass(
  adj, lengths, y = c(0, 0.2, 1.8, 2), n.eigenpairs = 4L,
  eta.grid = c(0.1, 1), eigen.solver = "dense"
)

Description

Experimental R reference implementation of prediction-synchronized local polynomial smoothing (PS-LPS). For fixed local polynomial parameters, PS-LPS fits one local polynomial chart per anchor and adds a quadratic penalty that synchronizes chart predictions on overlap points.

Usage

fit.ps.lps(
  X,
  y,
  foldid = NULL,
  support.size = NULL,
  degree = 2L,
  kernel = "gaussian",
  chart.dim = NULL,
  support.grid = NULL,
  degree.grid = NULL,
  kernel.grid = NULL,
  auto.chart.support.metric = c("coordinates", "operator", "both"),
  auto.chart.selection.metric = c("coordinates", "operator"),
  chart.dim.grid = NULL,
  selection.strategy = c("grid", "sparse_kd", "plateau_kd"),
  chart.dim.max = NULL,
  design.margin = 2L,
  lambda.sync.grid = c(0, 0.001, 0.01, 0.1, 1, 10),
  lambda.sync.search = c("grid", "guarded"),
  lambda.sync.selection = c("cv", "fixed"),
  local.candidate.search = c("screened", "full", "subgrid"),
  local.candidate.search.control = list(),
  lambda.sync.search.control = list(),
  lambda.diagnostics = c("all", "selected"),
  lambda.ridge = 1e-08,
  design.basis = c("monomial", "weighted.qr", "weighted.qr.drop",
    "orthogonal.polynomial.drop"),
  design.drop.tol = sqrt(.Machine$double.eps),
  ridge.multiplier.grid = NULL,
  ridge.condition.max = Inf,
  sync.neighbor.size = NULL,
  overlap.weight = c("normalized.product", "product"),
  chart.activation = c("none", "subject.od"),
  chart.activation.response = NULL,
  chart.activation.control = list(),
  ps.lps.geometry.cache = NULL,
  ps.lps.local.pca.supports = NULL,
  cv.folds = 5L,
  cv.seed = 1L
)

Arguments

Value

A list with fitted values, selected lambda, CV table, diagnostics, and fitted chart coefficients.

Examples

X <- matrix(seq(0, 1, length.out = 20), ncol = 1)
fit.ps.lps(
  X, sin(2 * pi * X[, 1]), support.size = 7L, degree = 1L,
  chart.dim = 1L, lambda.sync.grid = 0,
  lambda.sync.selection = "fixed", cv.folds = 2L
)

Description

Fits experimental PTTF regression models from a transported difference operator produced by pttf.operator. The implementation focuses on one response vector and explicit operator-row policies, especially the 1-D boundary policy used to separate interior transported rows from endpoint rows.

Usage

fit.pttf.trend.filtering(
  geometry = NULL,
  operator = NULL,
  X = NULL,
  y,
  derivative.order = 3L,
  penalty = c("l1", "l2"),
  lambda.grid = NULL,
  lambda.selection = c("cv", "fixed"),
  n.lambda = 80L,
  nfolds = 5L,
  foldid = NULL,
  cv.loss = c("mse", "mae"),
  selection = c("min", "one.se"),
  lambda.extension = c("none", "auto"),
  solver = c("genlasso", "admm", "auto"),
  operator.row.policy = c("all", "drop.line.boundary", "diagnostic.only"),
  line.order = NULL,
  boundary.trim = NULL,
  row.mass.rule = c("none", "node.mass"),
  row.normalize = c("none", "l2"),
  weights = NULL,
  maxsteps = 2000L,
  minlam = 0,
  approx = FALSE,
  rtol = 1e-07,
  btol = 1e-07,
  eps = 1e-04,
  admm.rho = NULL,
  admm.maxiter = 2000L,
  admm.abstol = 1e-04,
  admm.reltol = 0.001,
  verbose = FALSE,
  admm.adaptive.rho = TRUE,
  ...
)

Arguments

Value

A list of class "pttf.trend.filtering.fit".

Examples

n <- 10L
X <- matrix(seq(0, 1, length.out = n), ncol = 1)
adj <- lapply(seq_len(n), function(i) intersect(c(i - 1L, i + 1L), 1:n))
lengths <- Map(function(i, j) abs(X[j, 1] - X[i, 1]), seq_len(n), adj)
geometry <- pttf.geometry(
  X, adj, lengths, graph = "supplied", tangent.dim = 1L
)
operator <- pttf.operator(geometry, derivative.order = 2L)
fit.pttf.trend.filtering(
  operator = operator, y = sin(2 * pi * X[, 1]), penalty = "l2",
  lambda.grid = 0.1, lambda.selection = "fixed"
)

Description

Fits the fixed-operator S-LPL-TF estimator

$\hat f = \arg\min_f \frac12\|y-f\|_2^2 + \lambda_1\|A_{\mathrm{LPL}}f\|_1 + \frac{\lambda_2}{2}\|C_{\mathrm{sync}}f\|_2^2.$

Usage

fit.slpl.tf(
  X = NULL,
  y,
  operator = NULL,
  lambda1 = NULL,
  lambda2 = 0,
  lambda1.grid = NULL,
  lambda2.grid = NULL,
  lambda.selection = c("fixed", "cv"),
  operator.grid = NULL,
  foldid = NULL,
  cv.folds = 5L,
  cv.loss = c("mse", "rmse", "mae"),
  cv.seed = NULL,
  cv.repeats = 1L,
  solver = c("genlasso"),
  selection = c("min", "one.se"),
  maxsteps = 2000L,
  minlam = 0,
  approx = FALSE,
  rtol = 1e-07,
  btol = 1e-07,
  eps = 1e-04,
  verbose = FALSE,
  ...
)

Arguments

Value

A list of class "slpl_tf".

Examples

if (requireNamespace("genlasso", quietly = TRUE)) {
  X <- matrix(seq(0, 1, length.out = 18), ncol = 1)
  fit.slpl.tf(X, sin(2 * pi * X[, 1]), degree = 1L,
              support.type = "knn", support.size = 7L,
              lambda1 = 0.1, lambda2 = 0, lambda.selection = "fixed")
}

Description

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

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.

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.

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"
)

Description

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

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$.

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.

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
)
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)
)
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)
)

Description

Convenience wrapper that converts subject visit indices into a density input. Density-native methods dispatch to fit.density. LPS/PS-LPS methods fit ordinary smoother objects first and then call normalize.density; they are subject-occupation workflows, not standalone density-native methods.

Usage

fit.subject.od(
  X,
  subject.index,
  method = c("empirical", "graph_random_walk", "lps_count", "ps_lps_count",
    "lps_logistic_binary", "chart_kernel", "local_likelihood_density",
    "local_likelihood_bernoulli"),
  graph = NULL,
  graph.control = list(),
  od.control = list(),
  return.details = TRUE,
  od.cv = c("none", "visit"),
  visit.foldid = NULL,
  visit.cv.folds = 5L,
  visit.cv.seed = 1L,
  visit.cv.epsilon = 1e-15,
  ...
)

Arguments

Details

The "lps_logistic_binary" method name is historical. In the current OD workflow it calls fit.lps(..., outcome.family = "bernoulli"), which fits a clipped probability field with the LPS least-squares core and then converts that field to a density with normalize.density. It is not the outcome.family = "binomial" local-logistic IRLS path.

When od.cv = "visit", graph random-walk candidates may pass OD-level grids inside graph.control: walk.step.grid, affinity.method.grid, affinity.scale.grid, affinity.epsilon.grid, and normalize.grid. Missing affinity.scale values mean the usual data-derived affinity scale.

When od.cv = "visit", chart-based and LPS-family methods may pass chart.dim.grid through .... For "chart_kernel", "local_likelihood_density", "local_likelihood_bernoulli", "lps_count", "lps_logistic_binary", and "ps_lps_count", this grid compares fixed integer chart dimensions with optional "auto" and "local.auto" chart-dimension policies under the same held-out-visit negative-log-occupation score. For LPS-family OD visit CV, each outer candidate is passed to the source smoother as a scalar local model configuration; this avoids nested row-level multi-candidate selection inside each held-out-visit fold.

Value

A list of class "density_fit". Its rho component is the estimated probability mass over the rows of X, and empirical.rho is the normalized visit-count mass. The subject component summarizes the number of visits, the number of distinct visited support points, the largest visit multiplicity, and the fraction of repeatedly visited support points. Other components describe the fitted method, normalization accounting, smoothing diagnostics, and warnings as documented in fit.density. With od.cv = "visit" and return.details = TRUE, the result also contains the visit-level cross-validation table, fold assignments, and held-out predicted masses.

Examples

X <- matrix(seq(0, 1, length.out = 6), ncol = 1)
fit.subject.od(X, subject.index = c(1L, 3L, 3L, 6L),
               method = "empirical")

Description

Identifies the boundary vertices of a specified region in a graph, defined as vertices in the region that have at least one neighbor outside the region.

Usage

get.region.boundary(adj.list, region)

Arguments

Details

This function determines which vertices in a specified region are positioned at the boundary - meaning they have at least one neighbor that is not part of the region. These boundary vertices are often treated differently in graph-based algorithms, such as harmonic smoothing, where their values are typically fixed as constraints.

The boundary definition used matches the one implemented in the C++ harmonic smoothing functions, focusing on vertices that have connections to the outside of the region.

Value

An integer vector containing the indices of the boundary vertices, which is a subset of the input region vector. Returns an empty integer vector if no boundary vertices exist.

See Also

perform.harmonic.smoothing, harmonic.smoother

Examples

# Create a simple grid graph adjacency list (5x5 grid)
create_grid_adj_list <- function(n_rows, n_cols) {
  n <- n_rows * n_cols
  adj_list <- vector("list", n)
  for (i in 1:n_rows) {
    for (j in 1:n_cols) {
      v <- (i-1) * n_cols + j
      neighbors <- numeric(0)

      # Add neighbors (up, down, left, right)
      if (i > 1) neighbors <- c(neighbors, (i-2) * n_cols + j)  # up
      if (i < n_rows) neighbors <- c(neighbors, i * n_cols + j)  # down
      if (j > 1) neighbors <- c(neighbors, (i-1) * n_cols + (j-1))  # left
      if (j < n_cols) neighbors <- c(neighbors, (i-1) * n_cols + (j+1))  # right

      adj_list[[v]] <- neighbors
    }
  }
  return(adj_list)
}

# Create a 5x5 grid graph
grid_adj_list <- create_grid_adj_list(5, 5)

# Define a region (central 3x3 subgrid)
central_region <- c(7:9, 12:14, 17:19)

# Find boundary vertices of the central region
boundary <- get.region.boundary(grid_adj_list, central_region)
print(boundary)
# Expected output: c(7, 8, 9, 12, 14, 17, 18, 19)
# (all except the center vertex 13)

Description

Builds the weighted graph difference operators used by graph trend filtering. The default operator.family = "graph.laplacian.recursive" uses recursive graph operators

$k=0:\quad \Delta_w^{(1)} = D_w,$

$k=1:\quad \Delta_w^{(2)} = L_w = D_w^\top D_w,$