Projection Pursuit with Order-2 Trend Filtering — fit.projection.trend.filtering • geosmooth

Fit a sum of piecewise-quadratic curves on fixed or learned linear projections. Gaussian models use squared loss; binomial models add the curves on the logit scale.

Usage

fit.projection.trend.filtering(
  X,
  y,
  lambda,
  family = c("gaussian", "binomial"),
  directions = NULL,
  direction.search = c("learned", "fixed", "gn", "pls"),
  nterms = NULL,
  weights = NULL,
  offset = NULL,
  curve.ridge = 0.001,
  preprocessing = c("scale", "isotropic", "none"),
  seed = 1L,
  control = list(),
  backend = c("auto", "dense", "banded"),
  warn = TRUE
)

Arguments

X: Finite numeric predictor matrix or data frame.
y: Finite numeric response. For family="binomial", values must lie in [0,1]; proportions can be accompanied by trial weights.
lambda: Nonnegative scalar or Inf, on the mean-loss scale.
family: "gaussian" or "binomial".
directions: Matrix with one direction per row and one column per original predictor. Directions refer to the preprocessed coordinates. Entries for dropped constant columns must be zero. In learned mode these are additional candidates; in fixed mode these are the fitted directions.
direction.search: "fixed" fits supplied directions (coordinate directions when omitted); "learned" selects and refines directions with pattern search. "gn" uses multistart damped Gauss–Newton proposals; "pls" uses one covariance start per component with the same Gauss–Newton refinement. The latter two modes construct their own starts and do not accept supplied directions.
nterms: Number of supplied directions in fixed mode; maximum number of components in learned mode. Defaults to the supplied row count, all nonconstant coordinates in fixed mode, or two in learned mode. Zero fits an intercept only.
weights: Positive finite observation weights, default one.
offset: Known additive predictor, default zero.
curve.ridge: Nonnegative component ridge coefficient on the mean-loss scale; default 0.001. Positive ridge stabilizes component decomposition.
preprocessing: "scale" centers and scales each nonconstant column; "isotropic" centers and uses a common scale, preserving rotations and distance ratios; "none" leaves retained columns unchanged. All transformations use training weights only. The default is "isotropic" for "gn" and "pls", and "scale" otherwise.
seed: Integer random seed. Local randomness leaves the caller's RNG state unchanged.
control: Named list. max.sweeps=200, tolerance=1e-7, gap.tolerance=1e-8, and objective.tolerance=1e-10 govern joint fitting. inner=list() passes controls to the univariate solver. joint.solver="backfit" retains alternating component updates; joint.solver="dense" attempts a simultaneous dense solve when there are exactly two components, with guarded fallback to backfitting. joint.max.workspace.bytes=1024^3 bounds the estimated dense workspace. The joint certificate and objective safeguard apply to both solvers. Search controls are starts=2, sparse.directions=8, dense.directions=8, refinement.passes=2, refinement.proposals=16, steps=c(0.25,0.1,0.04,0.01), and improvement.tolerance=1e-9. Zero refinement passes select from a fixed candidate bank. The "gn" and "pls" modes instead use gn.max.steps=8 accepted proposals and gn.max.fits=40 native fits per start, with gn.shared.fits=400 shared fits for previously selected components. stream.key=character() can identify a dataset or representation; local start streams also depend on the seed and component number.
backend: Univariate solver: "auto", "dense", or "banded".
warn: Warn when joint fitting or a terminal search failure makes the returned fit ineligible. Failed candidates are retained in search.history.

Value

A projection.trend.filtering.fit containing centered terms, intercept, canonical directions in retained predictor coordinates, component fits, saved preprocessing, checks, converged, status, objective.history, and a candidate-level search.history. checks$converged reports joint numerical optimality; converged also excludes terminal search failures. Search termination is recorded separately in search.status. Ordinary predictions are NA for failed fits. Gauss–Newton modes additionally return numerical.certified, cv.eligible, counters, curve.diagnostics, and search.banks. An interrupted search may preserve a numerical certificate while being ineligible for CV. Requested and effective component counts can differ. In counters, candidate.fits counts returned attempts, while completed.candidate.fits counts eligible, certified results.

Details

For training weight mass W, the objective is $$W^{-1}\sum_i w_i\ell(y_i,\eta_i) + \lambda\sum_j\|D_j^{(3)}h_j\|_1 + \frac{\rho}{2W}\sum_{ij}w_i h_j(u_j^T z_i)^2.$$ Each component has training weighted mean zero; directions have unit norm and a canonical sign. Projection ranges and the falling-factorial basis follow fit.gaussian.trend.filtering and fit.logistic.trend.filtering. Native lambda and curve ridge are each multiplied by W exactly once. Gaussian loss includes a factor of one half. lambda=Inf restricts components to quadratics; it does not select a constant model.

Fixed-direction fitting is convex. Backfitting updates a curve with the other components supplied as an offset and checks joint stationarity, dual feasibility, complementarity and the original objective after each sweep. The "learned" mode uses greedy additions and derivative-free tangent proposals, followed by joint refitting. The search is nonconvex; its recorded stopping status is not a global optimality certificate. Gauss–Newton modes use one-sided derivatives of the falling-factorial curve, solve in the tangent space of the unit direction, and accept only certified joint-objective improvements. Binomial curvature is floored at 1e-6 without changing the likelihood score. All curve fits use the original responses and offsets, including for binomial outcomes. Use preprocessing="isotropic" to preserve the geometry of supplied score coordinates. PCA and variance-threshold tuning are not performed automatically by this function. In "gn" mode the ordered starts are coordinate axes, weighted ridge, weighted covariance, sliced inverse regression, and twice as many random directions as retained coordinates. Sign-equivalent duplicates retain their first origin. SIR uses two outcome classes for the first Bernoulli component with unit trial weights; otherwise it slices the centered working response by weight mass, keeping exact response ties together. Both supervised seeds and fitted curves eliminate a weighted intercept. A start must obtain a completed joint certificate before it can be selected. A per-start limit can retain that certificate; exhaustion of the candidate or shared allowance makes the candidate ineligible. With q retained coordinates and requested K components, the default native-fit caps are 40*K*(3*q+3)+400 for "gn" and 40*K+400 for "pls". An intercept-only fit makes no native curve-fitting calls. Directions with fewer than four distinct projections are ineligible. Exact ties are retained. Distinct projected locations that collapse during floating-point normalization produce an explicit failed fit; locations are never rounded, merged or jittered to repair a candidate. One-class binomial responses have no finite unpenalized intercept and produce an error.

Examples

x <- seq(-1, 1, length.out = 24)
X <- cbind(x, cos(3*x))
fit <- fit.projection.trend.filtering(X, sin(2*x), lambda = 0.002,
    directions = matrix(c(1, 0), 1), direction.search = "fixed")
predict(fit, X[1:3, , drop = FALSE])
#> [1] -1.255071 -1.145935 -1.036798
#> attr(,"outside.range")
#>      term1
#> [1,] FALSE
#> [2,] FALSE
#> [3,] FALSE