tune_Lambda: Tune Smoothing and Ridge Penalties via Leave-One-Out or...

In lgspline: Lagrangian Multiplier Smoothing Splines for Smooth Function Estimation

Tune Smoothing and Ridge Penalties via Leave-One-Out or Generalized Cross-Validation

Description

Tunes the penalty matrix by minimizing either exact leave-one-out on the transformed problem or modified GCV. This is the top-level tuning entry point: it assembles \boldsymbol{\Lambda}, initializes from a small grid search, and then refines the penalty parameters by quasi-Newton optimization with either closed-form or finite-difference gradients.

Usage

tune_Lambda(
  y,
  X,
  X_gram,
  smoothing_spline_penalty,
  A,
  K,
  p_expansions,
  N_obs,
  opt,
  use_custom_bfgs,
  C,
  colnm_expansions,
  wiggle_penalty,
  flat_ridge_penalty,
  initial_wiggle,
  initial_flat,
  unique_penalty_per_predictor,
  unique_penalty_per_partition,
  penalty_vec,
  meta_penalty,
  family,
  unconstrained_fit_fxn,
  keep_weighted_Lambda,
  iterate,
  qp_score_function,
  quadprog,
  qp_Amat,
  qp_bvec,
  qp_meq,
  tol,
  sd_y,
  delta,
  tuning_criterion,
  gcv_gamma,
  constraint_value_vectors,
  parallel,
  parallel_eigen,
  parallel_trace,
  parallel_aga,
  parallel_matmult,
  parallel_qr,
  parallel_bfgs,
  parallel_grideval,
  qr_pivot_smoothing_constraints,
  parallel_unconstrained,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks,
  shared_env,
  custom_penalty_mat,
  order_list,
  glm_weight_function,
  schur_correction_function,
  need_dispersion_for_estimation,
  dispersion_function,
  observation_weights,
  homogenous_weights,
  blockfit,
  just_linear_without_interactions,
  Vhalf,
  VhalfInv,
  verbose,
  include_warnings,
  ...
)

Arguments

y	List; response vectors by partition.
X	List; design matrices by partition.
X_gram	List; Gram matrices by partition.
smoothing_spline_penalty	Matrix; integrated squared second derivative penalty.
A	Matrix; smoothness constraints at knots.
K	Integer; number of interior knots in 1-D, number of partitions - 1 in higher dimensions.
p_expansions	Integer; columns per partition.
N_obs	Integer; total sample size.
opt	Logical; TRUE to optimize penalties, FALSE to use initial values.
use_custom_bfgs	Logical; TRUE for analytic gradient BFGS as natively implemented, FALSE for finite differences as implemented by stats::optim(). The analytic path is typically faster; however, for exact LOO tuning the leverage derivative can be numerically delicate in some problems, in which case use_custom_bfgs = FALSE remains a conservative fallback.
C	Matrix; polynomial expansion matrix (used for initialization). This is the monomial expansion design used to derive starting values and is not the inequality-constraint matrix sometimes denoted by \mathbf{C} elsewhere in the package documentation.
colnm_expansions	Character vector; column names of the expansion matrix.
wiggle_penalty, flat_ridge_penalty	Fixed penalty values if provided.
initial_wiggle, initial_flat	Numeric vectors; candidate values for grid search initialization on the raw (non-negative) scale. Converted to log scale internally for optimization.
unique_penalty_per_predictor, unique_penalty_per_partition	Logical; allow predictor/partition-specific penalties.
penalty_vec	Numeric vector; initial values for predictor/partition penalties on the raw (non-negative) scale. Converted to log scale internally for optimization. Use c() when no per-predictor or per-partition penalties are needed.
meta_penalty	The "meta" ridge penalty, a regularization for predictor/partition penalties to pull them towards 1 on the raw scale.
family	GLM family with optional custom tuning loss.
unconstrained_fit_fxn	Function for unconstrained fitting.
keep_weighted_Lambda, iterate	Logical controlling GLM fitting.
qp_score_function, quadprog, qp_Amat, qp_bvec, qp_meq	Quadratic programming parameters (see arguments of lgspline).
tol	Numeric; convergence tolerance. Used directly for strict absolute criterion and log-parameter checks, and as the floor for the adaptive plateau check in the custom BFGS path.
sd_y, delta	Response standardization parameters.
tuning_criterion	Character scalar selecting the tuning objective: "loo" for exact leave-one-out on the transformed tuning problem or "gcv" for generalized cross-validation. The LOO criterion itself is still evaluated exactly on that transformed problem; the caveat is that its fully analytic leverage derivative can be numerically sensitive in some datasets. For very large samples, "gcv" is often the more practical choice; as a rough guideline, it is recommended once the sample size is above about 250,000.
gcv_gamma	Numeric scalar, at least 1, used only when tuning_criterion = "gcv". Multiplies the effective degrees of freedom in the GCV denominator.
constraint_value_vectors	List; constraint values.
parallel	Logical; enable parallel computation.
parallel_eigen, parallel_trace, parallel_aga	Logical; specific parallel flags.
parallel_matmult, parallel_qr, parallel_bfgs, parallel_grideval	qr_pivot_smoothing_constraints,parallel_unconstrained Logical; specific parallel flags. When parallel_grideval = TRUE or parallel_bfgs = TRUE, inner parallel flags such as parallel_eigen are ignored for those tuning fits.
cl	Parallel cluster object.
chunk_size, num_chunks, rem_chunks	Integer; parallel computation parameters.
shared_env	Environment; shared variables exported to cluster workers.
custom_penalty_mat	Optional custom penalty matrix.
order_list	List; observation ordering by partition.
glm_weight_function, schur_correction_function	Functions for GLM weights and corrections.
need_dispersion_for_estimation, dispersion_function	Control dispersion estimation.
observation_weights	Optional observation weights.
homogenous_weights	Logical; TRUE if all weights equal.
blockfit	Logical; when TRUE, the backfitting block decomposition (blockfit_solve) is used in place of get_B for the coefficient estimation step at each tuning-objective evaluation, provided just_linear_without_interactions is non-empty and K > 0. The dispatch condition is identical to that in lgspline.fit so that penalties are tuned under the same estimator that will be used for the final fit. Falls back to get_B automatically on failure.
just_linear_without_interactions	Numeric; vector of columns for non-spline effects without interactions.
Vhalf, VhalfInv	Square root and inverse square root correlation structure matrices. These are passed through to the coefficient-estimation step used inside each tuning-objective evaluation, so any dense or Woodbury-accelerated correlated solve is the same one used for the final fitted model.
verbose	Logical; print progress.
include_warnings	Logical; print warnings/try-errors.
...	Additional arguments passed to fitting functions.

Details

The tuning procedure consists of the following steps:

1. Preprocessing: Convert raw-scale penalties to log scale, compute cross-products, determine pseudocount delta, ensure constraint matrix compatibility.

2. Blockfit dispatch: Pre-compute flat_cols and the use_blockfit flag so that every tuning-objective evaluation uses the same coefficient estimator as the final fit. flat_cols are identified by matching column names against paste0("_", just_linear_without_interactions, "_") in colnm_expansions.

3. Grid search: Evaluate the selected tuning criterion over a grid of (wiggle, ridge) penalty candidates to find a good starting point (see .tune_grid_search). When parallel_grideval = TRUE and a cluster is supplied, these candidate fits are spread across the workers, and if more than six workers are available the grid is enlarged with additional random raw-scale candidates drawn from [\min/10,\max\times 10] in each penalty direction.

4. BFGS optimization: Minimize the selected tuning criterion via either the custom damped BFGS with closed-form gradients (see .damped_bfgs, .compute_gcvu_gradient, .compute_loocv_gradient) or base R's stats::optim with finite-difference gradients. The analytic paths differentiate the fitted criterion once to obtain \mathbf{M}_k=\partial\ell/\partial\boldsymbol{\Lambda}_k; all log-penalty gradients then reduce to blockwise trace products with the corresponding current penalty direction. For LOO, the leverage derivative may be numerically delicate in some problems and can motivate the finite-difference fallback. After the first 10 iterations, the custom BFGS path also stops on a small scale-aware plateau in accepted criterion values. When parallel_bfgs = TRUE and a cluster is supplied, each damped step evaluates a batch of halved step sizes in parallel and takes the best improving candidate from that batch. If inequality constraints are present, the active set from the last accepted coefficient solve is used as the initial working set for the next tuning evaluation.

5. Criterion choice: "loo" uses exact per-observation leverages on the transformed tuning problem, while "gcv" uses the usual trace-based denominator with optional gcv_gamma inflation.

6. Sample-size adjustment: After optimization, divide the tuned penalties by (N+1)/(N-1) (equivalently multiply by (N-1)/(N+1)) for both GCV- and LOO-based tuning so that the final penalties are decreased at small sample sizes.

7. Final Lambda: Compute the final penalty matrix from optimized parameters via compute_Lambda.

Parameterization: initial penalty values are accepted on the raw (non-negative) scale and converted to natural log-scale internally, i.e. raw_penalty = exp(theta), so that raw penalties are always positive. The chain rule factor d(exp(theta))/d(theta) = exp(theta) = raw_penalty.

If \mathbf{L}_{j,k} is the current contribution of \theta_j = \log\lambda_j to partition k, then

\frac{\partial \ell}{\partial \theta_j} = \sum_{k=0}^{K} \mathrm{tr}\{\mathbf{M}_k\mathbf{L}_{j,k}\}.

This avoids repeating the full derivative through \mathbf{G}, \mathbf{U}, and \mathbf{H} for each penalty.

The resulting penalty follows the same notation used in the paper:

\boldsymbol{\Lambda}_k = \lambda_w\Bigl(\boldsymbol{\Lambda}_s + \lambda_r\boldsymbol{\Lambda}_r + \sum_{l=1}^{L}\lambda_{l,k}\boldsymbol{\Lambda}_{l,k}\Bigr),

with full penalty

\boldsymbol{\Lambda} = \mathrm{blockdiag}(\boldsymbol{\Lambda}_0, \ldots, \boldsymbol{\Lambda}_K).

The current package stores these pieces using the implementation labels L1, L2, L_predictor_list, and L_partition_list, corresponding to \boldsymbol{\Lambda}_s, \boldsymbol{\Lambda}_r, and the additional \boldsymbol{\Lambda}_{l,k} components.

If correlation structure inputs are supplied, each tuning-objective evaluation calls the same constrained correlated solver used by the final model fit. In the structured-correlation case, the Woodbury reduction described in lgspline-details is inherited automatically through get_B; there is no separate tuning-specific notation or solver.

Value

List containing:

Lambda

Baseline per-partition penalty matrix corresponding to \lambda_w(\boldsymbol{\Lambda}_s + \lambda_r\boldsymbol{\Lambda}_r + \sum_{l=1}^{L}\lambda_{l,k}\boldsymbol{\Lambda}_{l,k}).

flat_ridge_penalty

Optimized ridge penalty.

wiggle_penalty

Optimized smoothing penalty.

other_penalties

Optimized additional penalties \lambda_{l,k}.

L_predictor_list

Implementation lists storing additional penalty matrices \boldsymbol{\Lambda}_{l,k} for predictor-specific directions.

L_partition_list

Implementation lists storing additional penalty matrices \boldsymbol{\Lambda}_{l,k} for partition-specific directions.

See Also

• optim for Hessian-free optimization

• compute_Lambda for penalty matrix construction

• compute_G_eigen for eigendecomposition of penalized Gram matrices

• get_B for constrained coefficient estimation

• blockfit_solve for the backfitting block-decomposition estimator used when blockfit = TRUE

Examples

## Not run: 
## ## Example 1: Basic usage within lgspline ## ## ## ## ## ## ## ## ## ## ##
## tune_Lambda is called internally by lgspline; direct calls are
#  for advanced users. Here we verify that the refactored version
#  produces identical output to the original.

set.seed(1234)
t <- runif(200, -5, 5)
y <- sin(t) + rnorm(200, 0, 0.5)

## Fit with automatic penalty tuning (calls tune_Lambda internally)
fit1 <- lgspline(t, y, K = 3)
cat("Wiggle penalty:", fit1$penalties$wiggle_penalty, "\n")
cat("Ridge penalty:", fit1$penalties$flat_ridge_penalty, "\n")
cat("Trace (edf):", fit1$trace_XUGX, "\n")

## ## Example 2: Fixed penalties (no tuning) ## ## ## ## ## ## ## ## ## ## ##
fit2 <- lgspline(t, y, K = 3, opt = FALSE,
                 wiggle_penalty = 1e-4,
                 flat_ridge_penalty = 0.1)
cat("Fixed wiggle:", fit2$penalties$wiggle_penalty, "\n")

## ## Example 3: blockfit path #### ## ## ## ## ## ## ## ## ## ## ## ## ## ##
## When blockfit = TRUE and just_linear_without_interactions is non-empty,
#  tune_Lambda dispatches to blockfit_solve at each GCV evaluation,
#  ensuring penalties are tuned under the same estimator used in the final
#  fit. Verify that tuned penalties are consistent across both paths.

set.seed(7)
n  <- 300
t1 <- runif(n, 0, 5)
t2 <- rnorm(n)
y2 <- sin(t1) + 0.5 * t2 + rnorm(n, 0, 0.3)
df <- data.frame(t1 = t1, t2 = t2)

## blockfit = TRUE uses blockfit_solve during tuning
fit_bf  <- lgspline(df, y2, K = 2, blockfit = TRUE,
                    just_linear_without_interactions = 2)
## blockfit = FALSE uses get_B during tuning (original path)
fit_std <- lgspline(df, y2, K = 2, blockfit = FALSE,
                    just_linear_without_interactions = 2)

cat("blockfit wiggle  :", fit_bf$penalties$wiggle_penalty, "\n")
cat("standard wiggle  :", fit_std$penalties$wiggle_penalty, "\n")
## Penalties may differ slightly; predictions should be close.
cat("Max pred diff:", max(abs(fit_bf$ytilde - fit_std$ytilde)), "\n")

## ## Example 4: Verify refactored subfunctions ## ## ## ## ## ## ## ## ## ##
## The internal .compute_meta_penalty should match hand calculation
mp <- lgspline:::.compute_meta_penalty(
  wiggle_penalty = 0.5,
  penalty_vec = c(predictor1 = 1.2, partition1 = 0.8),
  meta_penalty_coef = 1e-8,
  unique_penalty_per_predictor = TRUE,
  unique_penalty_per_partition = TRUE
)
expected <- 0.5 * 1e-8 * ((1.2 - 1)^2 + (0.8 - 1)^2) +
            0.5 * 1e-32 * (0.5 - 1)^2
stopifnot(abs(mp - expected) < 1e-20)
cat("Meta-penalty check passed.\n")

## ## Example 5: Verify gradient of meta-penalty ## ## ## ## ## ## ## ## ## #
gr <- lgspline:::.compute_meta_penalty_gradient(
  wiggle_penalty = 2.0,
  penalty_vec = c(predictor1 = 1.5),
  meta_penalty_coef = 1e-8,
  unique_penalty_per_predictor = TRUE,
  unique_penalty_per_partition = FALSE
)
# gr[1] should be 1e-32 * (2 - 1) * 2 = 2e-32
# gr[2] should be 0
# gr[3] should be 1e-8 * (1.5 - 1) * 1.5 = 7.5e-9
stopifnot(abs(gr[1] - 2e-32) < 1e-40)
stopifnot(gr[2] == 0)
stopifnot(abs(gr[3] - 7.5e-9) < 1e-17)
cat("Meta-penalty gradient check passed.\n")

## ## Example 6: Residual computation paths ## ##
## Identity link ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ##
r1 <- lgspline:::.compute_tuning_residuals(
  y = list(c(1, 2, 3)),
  preds = list(c(1.1, 1.9, 3.2)),
  delta = 0,
  family = gaussian(),
  observation_weights = list(NULL),
  K = 0,
  order_list = list(1:3)
)
stopifnot(max(abs(r1[[1]] - c(-0.1, 0.1, -0.2))) < 1e-10)
cat("Residual check passed.\n")

## End(Not run)

lgspline documentation built on May 8, 2026, 5:07 p.m.