Type:	Package
Title:	Lagrangian Multiplier Smoothing Splines for Smooth Function Estimation
Version:	1.0.1
Description:	Implements Lagrangian multiplier smoothing splines for flexible nonparametric regression and function estimation. Provides tools for fitting, prediction, and inference using a constrained optimization approach to enforce smoothness. Supports generalized linear models, Weibull accelerated failure time (AFT) models, quadratic programming constraints, and customizable working-correlation structures, with options for parallel fitting. The core spline construction builds on Ezhov et al. (2018) <doi:10.1515/jag-2017-0029>. Quadratic-programming and SQP details follow Goldfarb & Idnani (1983) <doi:10.1007/BF02591962> and Nocedal & Wright (2006) <doi:10.1007/978-0-387-40065-5>. For smoothing spline and penalized spline background, see Wahba (1990) <doi:10.1137/1.9781611970128> and Wood (2017) <doi:10.1201/9781315370279>. For variance-component and correlation-parameter estimation, see Searle et al. (2006) <ISBN:978-0470009598>. The default multivariate partitioning step uses k-means clustering as in MacQueen (1967).
License:	MIT + file LICENSE
Language:	en-US
Depends:	R (≥ 3.5.0)
Imports:	Rcpp (≥ 1.0.7), RcppArmadillo, FNN, RColorBrewer, plotly, quadprog, methods, stats
LinkingTo:	Rcpp, RcppArmadillo
Suggests:	testthat (≥ 3.0.0), spelling, knitr, rmarkdown, parallel, survival, MASS, graphics
URL:	https://github.com/matthewlouisdavisBioStat/lgspline
BugReports:	https://github.com/matthewlouisdavisBioStat/lgspline/issues
Encoding:	UTF-8
RoxygenNote:	7.3.2
Config/testthat/edition:	3
NeedsCompilation:	yes
Packaged:	2026-03-15 06:45:34 UTC; matth
Author:	Matthew Davis [aut, cre]
Maintainer:	Matthew Davis <matthewlouisdavis@gmail.com>
Repository:	CRAN
Date/Publication:	2026-03-15 12:10:07 UTC

Lagrangian Multiplier Smoothing Splines

Description

Implements Lagrangian multiplier smoothing splines for flexible nonparametric regression and function estimation. Provides tools for fitting, prediction, and inference using a constrained optimization approach to enforce smoothness. Supports generalized linear models, Weibull accelerated failure time (AFT) models, quadratic programming constraints, and customizable working-correlation structures, with options for parallel fitting. The core spline construction builds on Ezhov et al. (2018) doi:10.1515/jag-2017-0029. Quadratic-programming and SQP details follow Goldfarb & Idnani (1983) doi:10.1007/BF02591962 and Nocedal & Wright (2006) doi:10.1007/978-0-387-40065-5. For smoothing spline and penalized spline background, see Wahba (1990) doi:10.1137/1.9781611970128 and Wood (2017) doi:10.1201/9781315370279. For variance-component and correlation-parameter estimation, see Searle et al. (2006) <ISBN:978-0470009598>. The default multivariate partitioning step uses k-means clustering as in MacQueen (1967).

Author(s)

Maintainer: Matthew Davis matthewlouisdavis@gmail.com (ORCID)

See Also

Useful links:

• https://github.com/matthewlouisdavisBioStat/lgspline

• Report bugs at https://github.com/matthewlouisdavisBioStat/lgspline/issues

Efficient Matrix Multiplication Operator

Description

Operator wrapper around C++ efficient_matrix_mult() for matrix multiplication syntax.

This is an internal function meant to provide improvement over base R's operator for certain large matrix operations, at a cost of potential slight slowdown for smaller problems.

Usage

x %**% y

Arguments

Value

Matrix product of x and y

Examples

M1 <- matrix(1:4, 2, 2)
M2 <- matrix(5:8, 2, 2)
M1 %**% M2

Partition-Wise Active-Set Refinement for Inequality Constraints

Description

Replaces the dense .qp_refine SQP loop with a partition-wise active-set method that reuses the existing Lagrangian projection machinery. The key insight is that an active inequality constraint can be treated as an additional equality constraint and absorbed into the augmented constraint matrix \mathbf{A}_{\mathrm{aug}}, after which the standard \mathbf{G}^{1/2}\mathbf{r}^* trick applies without ever forming the full P \times P information matrix.

Usage

.active_set_refine(
  result,
  X,
  y,
  K,
  p_expansions,
  A,
  R_constraints,
  constraint_value_vectors,
  Lambda,
  Ghalf,
  GhalfInv,
  family,
  qp_Amat,
  qp_bvec,
  qp_meq,
  Xy_or_uncon,
  is_path3,
  parallel_aga,
  parallel_matmult,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks,
  tol,
  max_as_iter = 50
)

Arguments

Details

The active-set loop:

1. Solve the equality-constrained subproblem (Lagrangian projection) with current active set.

2. Check primal feasibility: any violated inactive inequalities?

3. Check dual feasibility: any negative multipliers on active inequalities?

4. If both satisfied, return (KKT conditions met).

5. Otherwise, add most-violated constraint or drop most-negative multiplier, and repeat.

Falls back to .qp_refine if the active-set method does not converge within max_as_iter iterations.

Value

A list with components:

result

List of refined coefficient column vectors by partition.

qp_info

List with active constraint information, or NULL.

converged

Logical; TRUE if active-set method converged.

Assemble qp_info from a solve.QP Solution

Description

Thin wrapper around .solver_assemble_qp_info(). Keeping the local helper name avoids changing the public return shape or downstream references inside blockfit_solve().

Usage

.bf_assemble_qp_info(
  last_qp_sol,
  beta_block,
  qp_Amat_combined,
  qp_bvec_combined,
  qp_meq_combined,
  converged,
  final_deviance,
  info_matrix = NULL
)

Arguments

Value

A list suitable for the qp_info slot, or NULL if last_qp_sol is NULL.

Gaussian Identity + GEE: Closed-Form Solve (Case a)

Description

When the response is Gaussian with identity link and a working correlation structure is present, the whitened system \mathbf{V}^{-1/2}\mathbf{X} is not block-diagonal, so partition-wise backfitting does not apply. This function performs a single closed-form solve that replicates get_B Path 1a but in the backfitting context.

Usage

.bf_case_gauss_gee(
  X,
  y,
  K,
  p_expansions,
  order_list,
  VhalfInv,
  Lambda,
  L_partition_list,
  unique_penalty_per_partition,
  A,
  constraint_values,
  spline_cols,
  flat_cols,
  observation_weights
)

Arguments

Value

A named list:

beta_spline

List of K+1 spline coefficient vectors.

beta_flat

Numeric vector of flat coefficients.

result

List of K+1 full per-partition coefficient vectors.

Gaussian Identity Backfitting Without Correlation (Case b)

Description

Standard block-coordinate descent on the convex quadratic objective, alternating between the spline step (Lagrangian projection) and the flat step (pooled penalized regression).

Usage

.bf_case_gauss_no_corr(
  X_spline,
  X_flat,
  y,
  K,
  Ghalf_sp,
  GhalfA_sp,
  XfXf_pen_inv,
  constraint_values_spline,
  nc_spline,
  nc_flat,
  A_spline,
  tol,
  max_backfit_iter,
  verbose,
  split = NULL,
  constraint_values = list(),
  Lambda_flat = NULL,
  p_expansions = NULL
)

Arguments

Value

A named list with beta_spline and beta_flat.

GLM + GEE Two-Stage Estimation (Case c)

Description

Stage 1 runs damped Newton-Raphson backfitting on the unwhitened link-scale working response to produce a warm start. Stage 2 refines via damped SQP on the full whitened system, replicating get_B Path 1b.

Usage

.bf_case_glm_gee(
  X,
  y,
  K,
  p_expansions,
  flat_cols,
  split,
  family,
  order_list,
  glm_weight_function,
  schur_correction_function,
  qp_score_function,
  need_dispersion_for_estimation,
  dispersion_function,
  observation_weights,
  iterate,
  tol,
  max_backfit_iter,
  parallel_eigen,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks,
  schur_zero,
  quadprog,
  qp_Amat,
  qp_bvec,
  qp_meq,
  Lambda,
  L_partition_list,
  unique_penalty_per_partition,
  A,
  constraint_values,
  Vhalf,
  VhalfInv,
  verbose,
  ...
)

Arguments

Value

A named list with result, beta_spline, beta_flat, and qp_info.

GLM Without GEE: Damped Newton-Raphson + Backfitting (Case d)

Description

Damped Newton-Raphson outer loop wrapping the weighted backfitting inner loop. Each damped Newton-Raphson step computes weights from the current linear predictor, then calls .bf_inner_weighted.

Usage

.bf_case_glm_no_corr(
  X,
  y,
  K,
  p_expansions,
  split,
  family,
  order_list,
  glm_weight_function,
  need_dispersion_for_estimation,
  dispersion_function,
  observation_weights,
  iterate,
  tol,
  max_backfit_iter,
  parallel_eigen,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks,
  schur_zero,
  unique_penalty_per_partition,
  VhalfInv,
  verbose,
  constraint_values = list(),
  ...
)

Arguments

Value

A named list with beta_spline and beta_flat.

Constrained Flat Update Given Current Spline Solution

Description

Solves for flat coefficients subject to the residual equality constraint imposed by the full constraint system, conditional on the current spline block solution.

Usage

.bf_constrained_flat_update(
  XfWXf_pen,
  Xfr,
  A_full,
  constraint_values,
  beta_spline,
  flat_rows_all,
  spline_rows_all,
  spline_cols,
  flat_cols,
  p_expansions,
  K,
  nc_flat
)

Arguments

Details

When the full constraint system is \mathbf{A}^\top \boldsymbol{\beta} = \mathbf{c}, and we partition \boldsymbol{\beta} into spline and flat blocks, the flat update must satisfy:

\mathbf{A}_{\mathrm{flat}}^\top \boldsymbol{\beta}_{\mathrm{flat}} = \mathbf{c} - \mathbf{A}_{\mathrm{spline}}^\top \boldsymbol{\beta}_{\mathrm{spline}}

This function solves the penalized least-squares problem for flat coefficients subject to this residual equality constraint using a Lagrangian approach.

Value

Numeric vector of flat coefficients (nc_f \times 1).

Compute Deviance for Non-GEE Models

Description

Evaluates the mean deviance (or mean squared error as fallback) at the current fitted values. Used inside the damped Newton-Raphson outer loop of blockfit_solve for convergence monitoring.

Usage

.bf_deviance(y_vec, mu_vec, obs_wt, ord, fam, ...)

Arguments

Value

Scalar; the mean deviance contribution.

Compute Whitened Deviance for GEE Convergence Monitoring

Description

Evaluates deviance in the whitened (decorrelated) space when the family supplies fam$custom_dev.resids. In that case the raw deviance residuals are divided by \sqrt{W} and then pre-multiplied by \mathbf{V}^{-1/2} before squaring and averaging. If only fam$dev.resids is available, the function falls back to the usual mean deviance (and otherwise mean squared error). Used in the GEE refinement loop (Case c) of blockfit_solve.

Usage

.bf_gee_deviance(y_vec, mu_vec, W_vec, ord, obs_wt, VhInv, fam, ...)

Arguments

Value

Scalar; the mean whitened deviance.

Numerical Derivative of the Inverse Link Function

Description

Returns \partial\mu / \partial\eta for the supplied family object, using the family's own mu.eta method when available and falling back to analytic forms for common links or central differences otherwise.

Usage

.bf_get_mu_eta(eta, fam)

Arguments

Value

Numeric vector of the same length as eta.

GLM Backfitting Inner Loop

Description

Given damped Newton-Raphson weights and response partitioned into lists, runs the inner backfitting loop that alternates between a weighted spline step and a weighted flat step until convergence. Shared by both Case c (GEE warm start) and Case d (GLM without GEE).

Usage

.bf_inner_weighted(
  X_spline,
  X_flat,
  z_list,
  W_list,
  K,
  Lambda_spline,
  Lambda_flat,
  L_part_spline,
  unique_penalty_per_partition,
  A_spline,
  constraint_values_spline,
  nc_spline,
  nc_flat,
  beta_spline_init,
  beta_flat_init,
  tol,
  max_backfit_iter,
  parallel_eigen,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks,
  split = NULL,
  constraint_values = list(),
  p_expansions = NULL
)

Arguments

Value

A named list with beta_spline and beta_flat.

Lagrangian Projection for the Spline-Only Subproblem

Description

Projects the adjusted cross-product \mathbf{X}_s^{\top} \tilde{\mathbf{y}} through \mathbf{G}_s^{1/2} and removes the component along \mathbf{G}_s^{1/2}\mathbf{A}_s to enforce the smoothness constraints in the spline block. When cv_spline contains nonzero values (e.g.\ from no_intercept), the corresponding contribution is added back after projection.

Usage

.bf_lagrangian_project(
  Xy_adj,
  Ghalf_cur,
  GhalfA_cur,
  cv_spline,
  K,
  nc_spline,
  A_spline
)

Arguments

Value

List of K+1 coefficient vectors (nc_s \times 1).

Construct Full Block-Diagonal Penalty Matrix

Description

Thin wrapper around .solver_build_lambda_block(). The block/backfitting solver keeps its original helper name so the surrounding code stays unchanged, while the shared implementation now lives in solver_utils.R.

Usage

.bf_make_Lambda_block(
  Lambda,
  K,
  unique_penalty_per_partition,
  L_partition_list
)

Arguments

Value

A (p_expansions \cdot (K+1)) \times (p_expansions \cdot (K+1)) matrix.

Split Design and Penalty into Spline and Flat Components

Description

Partitions the per-partition design matrices, penalty matrices, and constraint matrix into "spline" (columns receiving partition-specific coefficients) and "flat" (columns receiving a single shared coefficient across all partitions) subsets.

Usage

.bf_split_components(
  X,
  flat_cols,
  p_expansions,
  K,
  Lambda,
  L_partition_list,
  A,
  constraint_values
)

Arguments

Value

A named list with elements:

X_spline

List of K+1 matrices, spline columns only.

X_flat

List of K+1 matrices, flat columns only.

Lambda_spline

nc_s \times nc_s penalty submatrix.

Lambda_flat

nc_f \times nc_f penalty submatrix.

L_part_spline

List of partition-specific penalty submatrices for the spline columns.

A_spline

Spline-only constraint matrix (columns pruned and rank-reduced via pivoted QR).

nca_spline

Integer; number of columns in A_spline.

constraint_values_spline

List of constraint RHS vectors restricted to spline rows.

spline_cols

Integer vector of spline column indices.

nc_spline

Integer; number of spline columns.

nc_flat

Integer; number of flat columns.

A_flat

Flat-only constraint matrix (rows for flat coefficients, columns pruned and rank-reduced). Used to enforce mixed constraints on the flat update step.

nca_flat

Integer; number of columns in A_flat.

A_full

The original full constraint matrix A, retained for mixed-constraint enforcement.

flat_rows_all

Integer vector of flat-coefficient row indices in the full P-space.

spline_rows_all

Integer vector of spline-coefficient row indices in the full P-space.

has_mixed_constraints

Logical; TRUE if any constraint column in A has nonzero entries on both spline and flat rows.

mixed_constraint_cols

Integer vector of column indices in the original A that are mixed (touch both spline and flat rows).

Damped SQP Refinement on the Full (Optionally Whitened) System

Description

Runs a damped sequential quadratic programming loop using quadprog::solve.QP at each iteration, enforcing all smoothness equalities, flat-equality constraints, and any user-supplied inequality constraints.

Usage

.bf_sqp_loop(
  X_design,
  y_design,
  X_block_raw,
  beta_init,
  Lambda_block,
  qp_Amat_combined,
  qp_bvec_combined,
  qp_meq_combined,
  K,
  p_expansions,
  family,
  order_list,
  glm_weight_function,
  schur_correction_function,
  qp_score_function,
  need_dispersion_for_estimation,
  dispersion_function,
  observation_weights,
  iterate,
  tol,
  VhalfInv,
  VhalfInv_perm,
  is_gee,
  deviance_fun,
  X_partitions,
  y_partitions,
  verbose,
  ...
)

Arguments

Details

This function serves two roles within blockfit_solve:

1. GEE Stage 2 (Case c): operates on the whitened system \tilde{\mathbf{X}} = \mathbf{V}^{-1/2} \mathbf{X}_{\mathrm{block}}, initialized from the damped Newton-Raphson backfitting warm start.

2. Non-GEE SQP refinement: operates on the unwhitened block-diagonal design, applying inequality constraints after backfitting convergence.

Value

A named list with elements beta_block (final coefficient vector), result (per-partition coefficient list), last_qp_sol (last successful solve.QP output or NULL), converged (logical), and final_deviance (scalar).

Build Derivative QP Constraints in Full P-Dimensional Space

Description

Internal helper that constructs the Amat / bvec pair enforcing derivative sign constraints at every row of a block-diagonal design matrix.

Given an N_{\mathrm{sub}} \times P block-diagonal design matrix X_block, where P = p \times (K+1), this function:

1. Recovers the per-partition p-column expansion matrix C_qp and records each row's partition assignment.

2. Calls make_derivative_matrix on C_qp to obtain first or second derivative matrices with respect to each predictor.

3. Optionally selects only derivatives for a subset of predictor variables (target_vars).

4. Maps each derivative row into the full P-dimensional coefficient space, yielding one constraint column per (observation, variable) pair.

The result is a constraint pair \mathbf{A}^{\top}\boldsymbol{\beta} \ge \mathbf{b} (with \mathbf{b} = \mathbf{0}) suitable for solve.QP.

Usage

.build_deriv_qp(
  X_block,
  sign_mult,
  just_first,
  p_expansions,
  K,
  colnm_expansions,
  power1_cols,
  power2_cols,
  nonspline_cols,
  interaction_single_cols,
  interaction_quad_cols,
  triplet_cols,
  include_2way_interactions,
  include_3way_interactions,
  include_quadratic_interactions,
  expansion_scales,
  target_vars = NULL,
  og_cols = NULL
)

Arguments

Value

A list with components:

Amat

P \times M constraint matrix, where M is the number of unique constraint columns (after deduplication).

bvec

Numeric vector of length M, all zeros.

meq

Integer, always 0 (inequality constraints).

See Also

process_qp, solve.QP

Build Block-Diagonal Penalty Matrix

Description

Thin wrapper around .solver_build_lambda_block(). Keeping the local helper name preserves existing call sites in get_B() while centralizing the shared implementation in solver_utils.R.

Usage

.build_lambda_block(Lambda, K, unique_penalty_per_partition, L_partition_list)

Arguments

Value

A P \times P block-diagonal matrix where P = p \times (K+1).

Build Tuning Environment

Description

Assembles a named list containing all pre-computed objects and configuration needed by the GCV_u evaluation and gradient functions during penalty tuning. This avoids deep nesting of closures and makes the dependencies explicit.

Usage

.build_tuning_env(
  y,
  X,
  X_gram,
  Xy,
  smoothing_spline_penalty,
  A,
  R_constraints,
  K,
  p_expansions,
  N_obs,
  custom_penalty_mat,
  colnm_expansions,
  unique_penalty_per_predictor,
  unique_penalty_per_partition,
  meta_penalty,
  family,
  delta,
  order_list,
  observation_weights,
  homogenous_weights,
  parallel,
  parallel_eigen,
  parallel_trace,
  parallel_aga,
  parallel_matmult,
  parallel_unconstrained,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks,
  unconstrained_fit_fxn,
  keep_weighted_Lambda,
  iterate,
  qp_score_function,
  quadprog,
  qp_Amat,
  qp_bvec,
  qp_meq,
  tol,
  sd_y,
  constraint_value_vectors,
  glm_weight_function,
  schur_correction_function,
  need_dispersion_for_estimation,
  dispersion_function,
  blockfit,
  just_linear_without_interactions,
  Vhalf,
  VhalfInv,
  verbose,
  include_warnings,
  flat_cols,
  use_blockfit
)

Arguments

Details

In addition to the standard tuning arguments, the environment stores two pre-computed blockfit dispatch items:

use_blockfit

Logical; TRUE when blockfit is enabled, flat_cols is non-empty, and K > 0. Mirrors the dispatch logic in lgspline.fit so that the same fitting path is used during tuning as during the final fit.

flat_cols

Integer vector; column indices of non-interactive linear terms derived from just_linear_without_interactions and colnm_expansions. Pre-computed once here rather than re-derived at every GCV evaluation.

Value

Named list (the "tuning environment").

Check KKT Conditions for Partition-Wise Active-Set Method

Description

Given current constrained coefficient estimates and a set of active inequality constraints (treated as equalities in the last Lagrangian projection), checks primal feasibility of inactive constraints and dual feasibility (non-negative multipliers) of active constraints.

Usage

.check_kkt_partitionwise(
  result,
  Ghalf,
  GhalfInv,
  Xy_or_uncon,
  is_path3,
  A_aug,
  n_eq_orig,
  qp_Amat,
  qp_bvec,
  active_ineq,
  K,
  p_expansions,
  family,
  parallel_matmult,
  parallel_aga,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks,
  tol
)

Arguments

Details

Multipliers for active inequality constraints are recovered from the OLS fit used in the Lagrangian projection: the fitted coefficients on \mathbf{X}^* = \mathbf{G}^{1/2}\mathbf{A}_{\mathrm{aug}} give the Lagrangian multipliers (up to sign and scaling).

Value

A list with components:

feasible

Logical; TRUE if all inactive inequality constraints are satisfied within tolerance.

dual_feasible

Logical; TRUE if all active inequality multipliers are non-negative within tolerance.

violated

Integer vector; indices of violated inactive constraints (into qp_Amat columns).

drop

Integer vector; indices of active constraints with negative multipliers that should be dropped.

multipliers

Numeric vector of Lagrangian multipliers for the active inequality constraints.

Compute Partitioned GLS Cross-Products

Description

Computes \mathbf{X}_k^{\top}\mathbf{V}^{-1}\mathbf{y} for each partition, accounting for cross-partition contributions from the correlation structure.

Usage

.compute_Xy_V(X, y, VhalfInv_perm, K, p_expansions, order_list)

Arguments

Value

A list of K+1 column vectors, each p\_expansions \times 1.

Compute Euclidean distance matrix for a cluster block

Description

Returns pairwise Euclidean distances within the cluster indexed by inds. When spacetime has multiple columns, squared distances are averaged across dimensions before taking the square root.

Usage

.compute_dist_block(spacetime, inds)

Evaluate GCV_u Criterion at a Given Penalty Configuration

Description

Computes the GCV_u (unbiased generalized cross-validation) criterion for a given set of penalty parameters. This is the objective function minimized during penalty tuning.

Usage

.compute_gcvu(par, log_penalty_vec, env, ...)

Arguments

Value

List containing:

GCV_u

Numeric; GCV_u criterion value including meta-penalty.

B

List; fitted coefficient vectors by partition.

GXX

List; \mathbf{G}_{k} \mathbf{X}_{k}^{\top}\mathbf{X}_{k} matrices.

G_list

List; eigendecomposition results from compute_G_eigen.

mean_W

Numeric; mean of hat matrix diagonal.

sum_W

Numeric; trace of hat matrix.

Lambda

Matrix; combined penalty matrix.

L1

Matrix; smoothing spline penalty component.

L2

Matrix; ridge penalty component.

L_predictor_list

List; predictor-specific penalty matrices.

L_partition_list

List; partition-specific penalty matrices.

numerator

Numeric; sum of squared residuals.

denominator

Numeric; GCV denominator N(1 - \bar{W})^{2}.

residuals

List; residual vectors by partition.

denom_sq

Numeric; squared denominator.

AGAInv

Matrix; (\mathbf{A}^{\top}\mathbf{G}\mathbf{A})^{-1}.

Compute Closed-Form Gradient of GCV_u Criterion

Description

Computes the gradient of the GCV_u criterion with respect to the log-scale penalty parameters using analytical derivatives of the hat matrix trace and residual sum of squares.

Usage

.compute_gcvu_gradient(par, log_penalty_vec, outlist = NULL, env, ...)

Arguments

Details

The gradient is computed via:

\frac{\partial \mathrm{GCV}_u}{\partial \theta} = \frac{1}{D^{2}} \left( \frac{\partial N}{\partial \theta} D - N \frac{\partial D}{\partial \theta} \right)

where N = \sum r_{i}^{2} (numerator), D = n(1 - \bar{W})^{2} (denominator), \theta is the log-scale penalty parameter, and the chain rule d\lambda / d\theta = \lambda (exp parameterization) is applied.

For predictor- and partition-specific penalties, a trace-ratio heuristic is used:

\frac{\partial \mathrm{GCV}_u}{\partial \lambda_{j}} \approx \frac{\mathrm{tr}(\mathbf{L}_{j})}{\mathrm{tr}(\boldsymbol{\Lambda})} \cdot \frac{\partial \mathrm{GCV}_u}{\partial \lambda_{w}}

Value

List containing:

GCV_u

Numeric; GCV_u criterion value including meta-penalty.

gradient

Numeric vector; gradient on the log penalty scale.

outlist

List; GCV_u components (for reuse to avoid recomputation).

Compute Regularization (Meta) Penalty on Penalty Parameters

Description

Computes the regularization term that pulls predictor- and partition-specific penalty parameters toward 1 on the raw (positive) scale. This acts as a "meta-penalty" on the penalty magnitudes themselves.

Usage

.compute_meta_penalty(
  wiggle_penalty,
  penalty_vec,
  meta_penalty_coef,
  unique_penalty_per_predictor,
  unique_penalty_per_partition
)

Arguments

Details

The penalty takes the form:

0.5 \times c_{\mathrm{meta}} \times \sum_{j} (\lambda_{j} - 1)^{2} + 0.5 \times 10^{-32} \times (\lambda_{w} - 1)^{2}

where \lambda_{j} are predictor/partition penalties and \lambda_{w} is the wiggle penalty.

Value

Numeric scalar; the regularization penalty value.

Compute Gradient of Regularization (Meta) Penalty

Description

Computes the gradient of the meta-penalty with respect to the log-scale penalty parameters, incorporating the exp parameterization chain rule.

Usage

.compute_meta_penalty_gradient(
  wiggle_penalty,
  penalty_vec,
  meta_penalty_coef,
  unique_penalty_per_predictor,
  unique_penalty_per_partition
)

Arguments

Details

Under exp parameterization \lambda = \exp(\theta):

\frac{\partial}{\partial \theta} \left[ 0.5 c (\exp(\theta) - 1)^{2} \right] = c (\lambda - 1) \lambda

Value

Numeric vector; gradient of the meta-penalty on the log scale. Length equals 2 + length(penalty_vec).

Compute Partitioned GEE Score Vector

Description

Computes the GEE score vector \mathbf{X}^{\top}\mathrm{diag} (\mathbf{W})\mathbf{V}^{-1}(\mathbf{y} - \boldsymbol{\mu}) split into per-partition pieces of dimension p \times 1 each. Uses the identity \mathbf{V}^{-1} = \mathbf{I} + \boldsymbol{\Delta}_V where \boldsymbol{\Delta}_V = \mathbf{V}^{-1} - \mathbf{I} is precomputed and fixed across iterations, avoiding explicit formation of \mathbf{V}^{-1}.

Usage

.compute_score_V_partitioned(
  X,
  X_block,
  y,
  result,
  K,
  p_expansions,
  family,
  W,
  Delta_V,
  observation_weights
)

Arguments

Value

A list of K+1 column vectors, each p\_expansions \times 1, representing the partition-wise components of the full GEE score.

Compute Pseudocount Delta for Link Function Stabilization

Description

Determines the pseudocount \delta used to stabilize link function transformations during GCV penalty tuning. For identity link or when the response is naturally in the domain of the link function, returns 0. Otherwise, finds the \delta that makes the transformed response distribution most closely approximate a t-distribution.

Usage

.compute_tuning_delta(family, unl_y, N_obs, observation_weights, opt)

Arguments

Value

Numeric scalar; the pseudocount \delta \geq 0.

Compute Predictions During Penalty Tuning

Description

Wrapper around matmult_block_diagonal for computing partition-wise predictions \mathbf{X}_{k} \boldsymbol{\beta}_{k} during GCV penalty tuning.

Usage

.compute_tuning_predictions(
  X,
  B,
  K,
  parallel,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks
)

Arguments

Value

List of prediction vectors, one per partition.

Compute Residuals for GCV Criterion During Penalty Tuning

Description

Computes residuals used in the numerator of the GCV criterion. Handles identity link, general GLM link functions with pseudocount delta, custom deviance residual functions, and observation weights.

Usage

.compute_tuning_residuals(
  y,
  preds,
  delta,
  family,
  observation_weights,
  K,
  order_list,
  ...
)

Arguments

Details

Three computation paths:

1. Identity link or no custom deviance: Standard g(y) - \hat{\eta} residuals on the link scale, optionally weighted by observation weights for non-Gaussian families. For Gaussian identity link with heterogeneous weights, the weights have already been absorbed into X and y prior to this call.

2. Custom deviance residuals: Delegates to family$custom_dev.resids(y, mu, order_indices, family, observation_weights, ...).

Value

List of residual vectors, one per partition.

Damped BFGS Optimizer for GCV Penalty Tuning

Description

Custom implementation of damped BFGS quasi-Newton optimization for minimizing the GCV_u criterion. Uses step-size damping with backtracking and Sherman-Morrison-Woodbury inverse Hessian updates.

Usage

.damped_bfgs(
  par,
  log_penalty_vec,
  gcvu_fxn,
  gr_fxn,
  env,
  tol,
  max_iter = 100,
  ...
)

Arguments

Details

The optimizer uses the following strategy:

1. Iterations 1-2: steepest descent with damping.

2. Iteration 3+: BFGS quasi-Newton with inverse Hessian approximation updated via the standard secant condition. Falls back to identity matrix when the update is numerically unstable.

3. Step acceptance: Armijo-like criterion (accept if \mathrm{GCV}_{u}^{(\mathrm{new})} \leq \mathrm{GCV}_{u}^{(\mathrm{old})}).

4. Backtracking: damping factor halved on rejection; terminates when damp < 2^{-10} (early iterations) or 2^{-12} (later iterations).

Value

List containing:

par

Numeric vector; best log-scale penalty parameters found.

gcv_u

Numeric; best GCV_u value achieved.

iterations

Integer; number of iterations performed.

Detect Whether Inequality Constraints Require Global (Dense) QP

Description

Thin wrapper around .solver_detect_qp_global(). This preserves the original helper name inside get_B() while using the shared detection logic defined in solver_utils.R.

Usage

.detect_qp_global(qp_Amat, p_expansions, K)

Arguments

Value

Logical; TRUE if dense QP is needed, FALSE if partition-wise active-set is valid.

Extract Per-Partition Diagonal Blocks from a Full Matrix

Description

Given a P \times P matrix (e.g., the full-system \mathbf{G} = (\mathbf{X}^{\top}\mathbf{V}^{-1}\mathbf{X} + \boldsymbol{\Lambda})^{-1}), extracts the p \times p diagonal block for each of the K+1 partitions.

Usage

.extract_G_diagonal(M, p_expansions, K)

Arguments

Value

A list of length K+1, each element a p\_expansions \times p\_expansions matrix corresponding to the diagonal block for partition k.

Fit Coefficients During GCV Tuning: blockfit_solve or get_B

Description

Dispatches to blockfit_solve when the blockfit conditions are met (i.e. env$use_blockfit is TRUE), otherwise calls get_B. On blockfit_solve failure, falls back to get_B automatically.

Usage

.fit_coefficients(G_list, Lambda, L_partition_list, env, return_G_getB, ...)

Arguments

Details

The blockfit condition mirrors lgspline.fit: blockfit && length(flat_cols) > 0 && K > 0, pre-computed in tune_Lambda and stored in env$use_blockfit.

return_G_getB is set to TRUE by the callers so that B_list$G_list contains the updated G matrices (after any GLM weight iteration inside get_B or blockfit_solve). These are needed immediately after this call for AGAmult_wrapper, GXX, and the trace computation.

Value

List; output of blockfit_solve or get_B, containing at minimum $B (coefficient list) and $G_list.

Call get_B During GCV Tuning

Description

Internal wrapper that calls get_B with all arguments drawn from the tuning environment env. Separated from .fit_coefficients so the fallback path in .fit_coefficients is clean and does not repeat the full argument list.

Usage

.fit_get_B(G_list, Lambda, L_partition_list, env, return_G_getB, ...)

Arguments

Path 2: Gaussian Identity Link, No Correlation

Description

For the canonical Gaussian case without correlation structures, the unconstrained estimate has the closed form \hat{\boldsymbol{\beta}}_k = \mathbf{G}_k \mathbf{X}_k^{\top} \mathbf{y}_k and the constrained estimate follows by a single Lagrangian projection. No iteration is needed.

Usage

.get_B_gaussian_nocorr(
  Xy,
  Ghalf,
  GhalfInv,
  A,
  K,
  p_expansions,
  R_constraints,
  constraint_value_vectors,
  family,
  return_G_getB,
  quadprog,
  qp_Amat,
  qp_bvec,
  qp_meq,
  qp_score_function,
  parallel_aga,
  parallel_matmult,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks,
  X,
  y,
  Lambda,
  order_list,
  observation_weights,
  iterate,
  tol,
  glm_weight_function,
  schur_correction_function,
  need_dispersion_for_estimation,
  dispersion_function,
  unique_penalty_per_partition,
  L_partition_list,
  VhalfInv,
  homogenous_weights,
  ...
)

Arguments

Details

When K = 0 and there are no inequality constraints or nonzero constraint values, the function returns \hat{\boldsymbol{\beta}} = \mathbf{G}\mathbf{X}^{\top}\mathbf{y} directly without forming the full P-dimensional OLS system.

Value

Same structure as get_B.

Path 1a: Gaussian Identity + GEE (Closed-Form Full-System Solve)

Description

Computes the constrained penalized GLS estimate for Gaussian response with identity link when a correlation structure is present. Because \mathbf{V}^{-1/2} couples all partitions, fitting must operate on the full P-dimensional whitened system rather than partition-wise.

Usage

.get_B_gee_gaussian(
  X_block,
  X_tilde,
  y_tilde,
  VhalfInv_perm,
  Lambda_block,
  A,
  K,
  p_expansions,
  constraint_value_vectors,
  family,
  return_G_getB,
  quadprog,
  qp_Amat,
  qp_bvec,
  qp_meq,
  qp_score_function,
  order_list,
  observation_weights,
  ...
)

Arguments

Details

The unconstrained GLS estimate is:

\hat{\boldsymbol{\beta}} = \mathbf{G}(\mathbf{X}^{*\top} \mathbf{y}^{*}), \quad \mathbf{G} = (\mathbf{X}^{*\top}\mathbf{X}^{*} + \boldsymbol{\Lambda})^{-1}

where \mathbf{X}^{*} = \mathbf{V}^{-1/2}\mathbf{X} and \mathbf{y}^{*} = \mathbf{V}^{-1/2}\mathbf{y}. The constrained estimate is then obtained via the \mathbf{G}^{1/2}\mathbf{r}^* trick in the full P-dimensional space.

Value

If return_G_getB = TRUE: list with B, G_list, and qp_info. Otherwise: list of B and qp_info.

Path 1b: Non-Gaussian GEE (Damped SQP with Full Whitened Design)

Description

Estimates constrained coefficients for non-Gaussian GLMs with correlation structures using damped Sequential Quadratic Programming (SQP) in the full P-dimensional whitened space. The whitened design \mathbf{X}^{*} = \mathbf{V}^{-1/2}\mathbf{X} is required because \mathbf{V}^{-1/2} couples all partitions.

Usage

.get_B_gee_glm(
  X_block,
  X_tilde,
  y_block,
  y_tilde,
  VhalfInv_perm,
  Lambda_block,
  A,
  K,
  p_expansions,
  constraint_value_vectors,
  family,
  return_G_getB,
  iterate,
  tol,
  qp_Amat,
  qp_bvec,
  qp_meq,
  qp_score_function,
  order_list,
  observation_weights,
  glm_weight_function,
  schur_correction_function,
  need_dispersion_for_estimation,
  dispersion_function,
  VhalfInv,
  ...
)

Arguments

Value

Same structure as .get_B_gee_gaussian.

Path 1b-Woodbury: Non-Gaussian GEE with Woodbury Acceleration

Description

Replacement for .get_B_gee_glm when the off-diagonal rank r is low. At each damped Newton iteration, the Woodbury decomposition is recomputed with updated GLM working weights \mathbf{W}(\boldsymbol{\beta}^{(s)}) (since the weighted perturbation \boldsymbol{\Delta}(\mathbf{W}) changes), and the constrained step is taken via .lagrangian_project_woodbury.

Usage

.get_B_gee_glm_woodbury(
  X,
  y,
  K,
  p_expansions,
  VhalfInv_perm,
  order_list,
  A,
  R_constraints,
  constraint_value_vectors,
  family,
  return_G_getB,
  iterate,
  tol,
  qp_Amat,
  qp_bvec,
  qp_meq,
  qp_score_function,
  observation_weights,
  Lambda,
  Lambda_block,
  unique_penalty_per_partition,
  L_partition_list,
  wb_decomp_init,
  wb_sqrt_init,
  glm_weight_function,
  schur_correction_function,
  need_dispersion_for_estimation,
  dispersion_function,
  VhalfInv,
  parallel_eigen,
  parallel_aga,
  parallel_matmult,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks,
  ...
)

Arguments

Details

The key precomputation exploited here is that \boldsymbol{\Delta}_V = \mathbf{V}^{-1} - \mathbf{I} and \boldsymbol{\Delta}_V \mathbf{X} are fixed across iterations. At each step, the weighted perturbation \boldsymbol{\Delta}(\mathbf{W}) = \mathbf{X}^{\top} \mathrm{diag}(\mathbf{W})(\mathbf{V}^{-1} - \mathbf{I})\mathbf{X} is computed in O(NP) using the precomputed \boldsymbol{\Delta}_V \mathbf{X}, then split into block-diagonal corrections (absorbed into the corrected Gram) and an off-diagonal low-rank remainder.

Falls back to the dense .get_B_gee_glm if the Woodbury decomposition becomes invalid (e.g.\ rank exceeds threshold or \mathbf{F} is not positive definite) at any iteration.

Value

Same structure as .get_B_gee_glm.

Path 1a-Woodbury: Gaussian GEE with Woodbury Acceleration

Description

Replacement for .get_B_gee_gaussian when the off-diagonal rank r of the \mathbf{V}^{-1} perturbation is low. Uses the conserved \mathbf{G}^{1/2}\mathbf{r}^* OLS trick with \mathbf{G}_V^{1/2} expressed as block-diagonal \mathbf{G}_c^{1/2} plus rank-r correction through \mathbf{U}_Q. The code structure closely mirrors .get_B_gaussian_nocorr.

Usage

.get_B_gee_woodbury(
  X,
  y,
  K,
  p_expansions,
  VhalfInv_perm,
  order_list,
  A,
  R_constraints,
  constraint_value_vectors,
  family,
  return_G_getB,
  quadprog,
  qp_Amat,
  qp_bvec,
  qp_meq,
  qp_score_function,
  observation_weights,
  wb_decomp,
  wb_sqrt,
  parallel_aga,
  parallel_matmult,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks,
  qp_global,
  ...
)

Arguments

Details

Cost is O(Kp^3 + Pr^2) compared to O(P^3) for the dense Path 1a.

Value

Same structure as .get_B_gee_gaussian.

Path 3: Non-Gaussian GLM, No Correlation

Description

For GLMs with non-identity links or non-Gaussian families (without correlation structures), unconstrained partition-wise estimates are first obtained via Newton-Raphson (or OLS for Gaussian identity), then the constrained estimate is computed by a single Lagrangian projection. For non-canonical links, \mathbf{G} depends on the current fitted values through the GLM working weights \mathbf{W}; the projection is therefore iterated, updating \mathbf{G} at the current constrained estimate until convergence.

Usage

.get_B_glm_nocorr(
  X,
  y,
  X_gram,
  Xy,
  Lambda,
  Ghalf,
  GhalfInv,
  A,
  K,
  p_expansions,
  R_constraints,
  constraint_value_vectors,
  family,
  return_G_getB,
  iterate,
  tol,
  quadprog,
  qp_Amat,
  qp_bvec,
  qp_meq,
  qp_score_function,
  unconstrained_fit_fxn,
  keep_weighted_Lambda,
  unique_penalty_per_partition,
  L_partition_list,
  parallel_eigen,
  parallel_aga,
  parallel_matmult,
  parallel_unconstrained,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks,
  order_list,
  observation_weights,
  glm_weight_function,
  schur_correction_function,
  need_dispersion_for_estimation,
  dispersion_function,
  VhalfInv,
  homogenous_weights,
  ...
)

Arguments

Value

Same structure as get_B.

Lagrangian Projection via OLS Reformulation

Description

Given unconstrained (or GEE-initialized) coefficient estimates \hat{\boldsymbol{\beta}}, computes the constrained estimate \tilde{\boldsymbol{\beta}} = \mathbf{U}\hat{\boldsymbol{\beta}} where \mathbf{U} = \mathbf{I} - \mathbf{G}\mathbf{A} (\mathbf{A}^{\top}\mathbf{G}\mathbf{A})^{-1}\mathbf{A}^{\top}.

Usage

.lagrangian_project(
  GhalfXy,
  Ghalf,
  A,
  K,
  p_expansions,
  R_constraints,
  constraint_value_vectors,
  family,
  parallel_aga,
  parallel_matmult,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks
)

Arguments

Details

Rather than forming \mathbf{U} directly, the projection is reformulated as a residual from an OLS problem (the \mathbf{G}^{1/2}\mathbf{r}^* trick):

\mathbf{y}^* = \mathbf{G}^{-1/2}\hat{\boldsymbol{\beta}}

\mathbf{X}^* = \mathbf{G}^{1/2}\mathbf{A}

\mathbf{r}^* = (\mathbf{I} - \mathbf{X}^*(\mathbf{X}^{*\top} \mathbf{X}^*)^{-1}\mathbf{X}^{*\top})\mathbf{y}^*

\tilde{\boldsymbol{\beta}} = \mathbf{G}^{1/2}\mathbf{r}^*

This avoids explicitly forming and inverting \mathbf{A}^{\top}\mathbf{G}\mathbf{A}; the most expensive step is the QR decomposition of the R \times R system inside .lm.fit, which is far cheaper than the full P \times P solve.

Value

A list of length K+1, each element a column vector of constrained coefficients \tilde{\boldsymbol{\beta}}_k.

Woodbury-Corrected Lagrangian Projection

Description

Performs the constrained Lagrangian projection using the conserved \mathbf{G}^{1/2}\mathbf{r}^* OLS trick, with \mathbf{G}_V^{1/2} expressed as block-diagonal \mathbf{G}_c^{1/2} plus a rank-r correction through \mathbf{U}_Q. The structure is identical to .lagrangian_project: form \mathbf{y}^*, form \mathbf{X}^*, run .lm.fit, back-transform. Steps 1, 2, and 4 each gain an additive rank-r correction at cost O(Pr), while step 3 is literally unchanged.

Usage

.lagrangian_project_woodbury(
  GhalfXy_V,
  Ghalf_corrected,
  A,
  K,
  p_expansions,
  R_constraints,
  constraint_value_vectors,
  family,
  wb_sqrt,
  parallel_aga,
  parallel_matmult,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks
)

Arguments

Details

The correction matrix \mathbf{F} (identity plus rank-r modification) factors as \mathbf{F}^{1/2} = \mathbf{I}_P - \mathbf{U}_Q \mathbf{C} \mathbf{U}_Q^{\top} where \mathbf{C} is r \times r diagonal. The inverse square root is \mathbf{F}^{-1/2} = \mathbf{I}_P + \mathbf{U}_Q \mathbf{C}_{\mathrm{inv}} \mathbf{U}_Q^{\top}.

Value

A list of K+1 coefficient column vectors.

Quadratic Programming Refinement for Inequality Constraints

Description

After the Lagrangian projection handles smoothness equality constraints, this function refines the estimate to satisfy additional inequality constraints (monotonicity, derivative sign, range bounds, or user-supplied constraints) via quadprog::solve.QP.

Usage

.qp_refine(
  result,
  X,
  y,
  K,
  p_expansions,
  A,
  Lambda,
  Lambda_block,
  family,
  iterate,
  tol,
  qp_Amat,
  qp_bvec,
  qp_meq,
  qp_score_function,
  order_list,
  glm_weight_function,
  schur_correction_function,
  need_dispersion_for_estimation,
  dispersion_function,
  observation_weights,
  VhalfInv,
  ...
)

Arguments

Details

The subproblem at each iteration is a second-order Taylor approximation of the penalized log-likelihood around the current iterate \boldsymbol{\beta}^*. Collecting terms, this yields:

\tilde{\boldsymbol{\beta}} = \arg\min_{\boldsymbol{\beta}} \left\{-\mathbf{d}^{\top}\boldsymbol{\beta} + \frac{1}{2} \boldsymbol{\beta}^{\top}\mathbf{G}^{-1}\boldsymbol{\beta}\right\} \quad \text{s.t.} \quad \mathbf{A}^{\top}\boldsymbol{\beta} = \mathbf{0}, \quad \mathbf{C}^{\top}\boldsymbol{\beta} \succeq \mathbf{c}

where \mathbf{d} is the score adjusted by the current iterate and \mathbf{c} is the constraint value vector. Step acceptance uses damped updates with deviance monitoring; see Nocedal and Wright (2006) for the general SQP framework.

Value

A list with components:

result

List of refined coefficient column vectors by partition.

qp_info

List with QP solve metadata including Lagrangian multipliers and the active constraint matrix.

Compute rank-based distance matrix for AR(1) structures

Description

Converts the unique within-block distances returned by .compute_dist_block() to integer lags 0, 1, 2, ... in increasing order.

Usage

.rank_dists(spacetime, inds)

Recompute G at Current Coefficient Estimates

Description

Thin wrapper around .solver_recompute_G_at_estimate(). The numerical core is shared with blockfit_solve() so the final G, Ghalf, and GhalfInv objects are constructed through one code path.

Usage

.recompute_G_at_estimate(
  X,
  y,
  result,
  K,
  Lambda,
  family,
  order_list,
  glm_weight_function,
  schur_correction_function,
  need_dispersion_for_estimation,
  dispersion_function,
  observation_weights,
  VhalfInv,
  parallel_eigen,
  parallel_matmult,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks,
  unique_penalty_per_partition,
  L_partition_list,
  ...
)

Arguments

Value

A list with components G, Ghalf, and GhalfInv, each a list of K+1 matrices.

Assemble qp_info from a solve.QP Solution

Description

Packages the output of quadprog::solve.QP into the qp_info list expected by downstream code (inference, generate_posterior, and varcovmat construction).

Usage

.solver_assemble_qp_info(
  last_qp_sol,
  beta_block,
  qp_Amat_combined,
  qp_bvec_combined,
  qp_meq_combined,
  converged,
  final_deviance,
  info_matrix = NULL
)

Arguments

Details

Active constraint columns are identified as all equality columns (1:qp_meq_combined) plus any inequality column whose Lagrange multiplier exceeds sqrt(.Machine$double.eps). This matches the convention used in .qp_refine().

Value

A list with elements solution, lagrangian, active_constraints, iact, Amat_active, bvec_active, meq_active, converged, and final_deviance, plus info_matrix, Amat_combined, bvec_combined, and meq_combined when info_matrix is supplied. Returns NULL if last_qp_sol is NULL.

Build Block-Diagonal Penalty Matrix

Description

Assembles the full P \times P block-diagonal penalty matrix \boldsymbol{\Lambda} from a shared per-partition penalty Lambda and optional partition-specific additive terms. P = p \times (K+1) where p is the number of basis terms per partition.

Usage

.solver_build_lambda_block(
  Lambda,
  K,
  unique_penalty_per_partition,
  L_partition_list
)

Arguments

Details

When unique_penalty_per_partition = TRUE, the k-th diagonal block is Lambda + L_partition_list[[k]]; otherwise every block is Lambda.

Value

A P \times P block-diagonal matrix, P = p \times (K+1).

Detect Whether Inequality Constraints Require a Dense Global QP

Description

Inspects the columns of qp_Amat to determine whether every inequality constraint is confined to a single partition block (block-separable) or whether any constraint couples coefficients across partitions.

Usage

.solver_detect_qp_global(qp_Amat, p_expansions, K)

Arguments

Details

Returns FALSE (partition-wise active-set is valid) when all columns have nonzeros in at most one block. Returns TRUE (dense SQP required) when any column spans multiple blocks.

Value

Logical scalar.

Recompute G, Ghalf, and GhalfInv at a Supplied Coefficient Estimate

Description

Given current constrained coefficient estimates result, recomputes the penalized information matrix \mathbf{G}_k = (\mathbf{X}_k^{\top}\mathbf{W}_k \mathbf{D}_k\mathbf{X}_k + \boldsymbol{\Lambda}_k)^{-1} and its matrix square roots for each partition.

Usage

.solver_recompute_G_at_estimate(
  X,
  y,
  result,
  K,
  Lambda,
  family,
  order_list,
  glm_weight_function,
  schur_correction_function,
  need_dispersion_for_estimation,
  dispersion_function,
  observation_weights,
  VhalfInv,
  parallel_eigen,
  parallel_matmult,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks,
  unique_penalty_per_partition,
  L_partition_list,
  ...
)

Arguments

Details

This is used after Newton-Raphson convergence in Path 3 of get_B() and at the final return in blockfit_solve() when return_G_getB = TRUE. The implementation matches .recompute_G_at_estimate() exactly; it exists here so both solvers can call the same numerical core.

Value

A list with components G, Ghalf, and GhalfInv, each a list of K+1 matrices. G[[k]] is computed as tcrossprod(Ghalf[[k]]) to guarantee exact symmetry.

Grid Search Initialization for Penalty Tuning

Description

Evaluates the GCV_u criterion over a grid of initial wiggle and ridge penalty values to find a good starting point for BFGS optimization.

Usage

.tune_grid_search(
  log_initial_wiggle,
  log_initial_flat,
  log_penalty_vec,
  gcvu_fxn,
  env,
  include_warnings,
  ...
)

Arguments

Value

Numeric vector of length 2; the best (log_wiggle, log_flat) found.

Woodbury Decomposition of Correlation Structure

Description

Given \mathbf{V}^{-1/2}, the partition-specific design matrices, and the penalty, decomposes the GLS information matrix as \mathbf{G}_V^{-1} = \mathbf{G}_c^{-1} + \boldsymbol{\Delta}_{\mathrm{off}} where \mathbf{G}_c absorbs within-partition corrections from \mathbf{V}^{-1} - \mathbf{I} and \boldsymbol{\Delta}_{\mathrm{off}} captures cross-partition coupling with effective rank r.

Usage

.woodbury_decompose_V(
  VhalfInv_perm,
  X,
  K,
  p_expansions,
  Lambda,
  Lambda_block,
  unique_penalty_per_partition,
  L_partition_list,
  order_list,
  parallel_eigen,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks,
  rank_threshold_fraction = 1/3,
  family
)

Arguments

Details

When r = 0 (no cross-partition coupling) or r \ge P/3 (low-rank approximation not worthwhile), returns use_woodbury = FALSE and the caller should use the dense Path 1 approach instead.

Value

A list with components:

use_woodbury

Logical; FALSE triggers dense Path 1 fallback.

Ghalf_corrected

List of corrected \mathbf{G}_{c,k}^{1/2}.

GhalfInv_corrected

List of corrected \mathbf{G}_{c,k}^{-1/2}.

G_corrected

List of corrected \mathbf{G}_{c,k}.

L

Off-diagonal low-rank factor (P \times r).

S_signs

Length-r vector of \pm 1.

r

Integer effective rank.

M

Precomputed r \times r Woodbury inner inverse (\mathrm{diag}(1/S) + \mathbf{L}^{\top}\mathbf{G}_c \mathbf{L})^{-1}.

GL

List of K+1 matrices \mathbf{G}_{c,k} \mathbf{L}[rows_k,], each p \times r.

Compute Woodbury Half-Square-Root Components

Description

Given the low-rank Woodbury components from .woodbury_decompose_V, precomputes the \mathbf{U}_Q, \mathbf{C}, and \mathbf{C}_{\mathrm{inv}} matrices needed to express \mathbf{G}_V^{1/2} and \mathbf{G}_V^{-1/2} as block-diagonal plus rank-r correction.

Usage

.woodbury_halfsqrt_components(Ghalf_corrected, L, S_signs, r, K, p_expansions)

Arguments

Details

The correction matrix \mathbf{F} (identity plus rank-r modification) satisfies \mathbf{G}_V = \mathbf{G}_c^{1/2} \mathbf{F}\,\mathbf{G}_c^{1/2} and its square root factors as \mathbf{F}^{1/2} = \mathbf{I}_P - \mathbf{U}_Q \mathbf{C} \mathbf{U}_Q^{\top}.

Value

A list with components:

valid

Logical; FALSE if \mathbf{F} is not positive definite (signals fallback to dense).

U_Q

P \times r left singular vectors of \mathbf{Q} = \mathbf{G}_c^{1/2}\mathbf{L}.

C

r \times r diagonal matrix \mathbf{I}_r - \mathrm{diag}(\sqrt{1 - d_{F,i}}).

C_inv

r \times r diagonal matrix \mathrm{diag}(1/\sqrt{1 - d_{F,i}}) - \mathbf{I}_r.

GchalfUQ

List of K+1 matrices, each p \times r: \mathbf{G}_{c,k}^{1/2} \mathbf{U}_Q[rows_k,].

Per-Iteration Weighted Woodbury Redecomposition

Description

Given the FIXED \boldsymbol{\Delta}_V = \mathbf{V}^{-1} - \mathbf{I} and precomputed \boldsymbol{\Delta}_V \mathbf{X} (both invariant across Newton iterations), together with the current iteration's GLM working weights \mathbf{W} and weighted Gram matrices, computes all Woodbury components needed for the constrained step.

Usage

.woodbury_redecompose_weighted(
  Delta_V,
  X_block,
  DV_X,
  W,
  X,
  K,
  p_expansions,
  X_gram_weighted,
  Lambda,
  schur_corrections,
  unique_penalty_per_partition,
  L_partition_list,
  parallel_eigen,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks,
  rank_threshold_fraction = 1/3,
  family
)

Arguments

Details

The weighted perturbation is decomposed as \boldsymbol{\Delta}(\mathbf{W}) = \mathbf{X}^{\top} \mathrm{diag}(\mathbf{W})(\mathbf{V}^{-1} - \mathbf{I})\mathbf{X}, split into block-diagonal corrections (absorbed into the corrected Gram) and an off-diagonal remainder factored at low rank.

Value

A list with the same structure as .woodbury_decompose_V:

use_woodbury

Logical; FALSE if rank exceeds threshold.

Ghalf_corrected

List of corrected \mathbf{G}_{c,k}^{1/2}.

GhalfInv_corrected

List of corrected \mathbf{G}_{c,k}^{-1/2}.

G_corrected

List of corrected \mathbf{G}_{c,k}.

L

Off-diagonal low-rank factor (P \times r).

S_signs

Length-r sign vector.

r

Integer effective rank.

M

Precomputed r \times r Woodbury inner inverse.

GL

List of per-partition \mathbf{G}_{c,k}\mathbf{L}[rows_k,].

Efficient Matrix Multiplication for \textbf{A}^{T}\textbf{G}\textbf{A}

Description

Efficient Matrix Multiplication for \textbf{A}^{T}\textbf{G}\textbf{A}

Usage

AGAmult_wrapper(
  G,
  A,
  K,
  p_expansions,
  R_constraints,
  parallel,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks
)

Arguments

Details

Computes \textbf{A}^{T}\textbf{G}\textbf{A} efficiently in parallel chunks using AGAmult_chunk().

Value

Matrix product \textbf{A}^{T}\textbf{G}\textbf{A}

Lagrangian Multiplier Smoothing Splines: Mathematical Details

Description

This document provides the mathematical and implementation details for Lagrangian Multiplier Smoothing Splines as implemented in lgspline.

This package provides exhaustive resources for fitting multivariate smoothing smoothing splines with a monomial basis and analytical form cubic smoothing spline penalty.

The material is presented such that a programmer or statistician of reasonable experience and background can understand and implement the procedure from scratch, and also potentially critique some of the modelling choices that went into designing this package.

Informally, lgspline answers the following question: How can we best adapt a useful functionality of basis splines, under the alternative interpretation of smoothing as explicit external constraints instead?

The obvious benefit is a much more flexible and interpretable final model that for non-experienced users is simply easier to understand without post-hoc processing, and for experienced users can be used to customize models more easily.

The drawback is that the interpretation of constraints as external adds a new layer of complexity to each step of the model fitting process, whereas for implicit design matrix construction these complications are bypassed.

While it is true a B-spline can always be converted back into monomial form, tensor-product splines that generalize this to multiple dimensions often explodes the number and degree of interaction terms, the conversion may not be computationally stable, and it is not available in standard software.

Statistical Problem Formulation

Consider an N \times q matrix of predictors \mathbf{T} = (\mathbf{t}_1, \dots, \mathbf{t}_N)^{\top} and an N \times 1 response vector \mathbf{y} = (y_1, \dots, y_N)^{\top}. We assume the relationship follows a generalized linear model with unknown smooth function f:

y_i \sim \mathcal{D}(g^{-1}(f(\mathbf{t}_i)),\, \sigma^2)

where \mathcal{D} is a distribution (e.g. exponential family or related) with mean \mu_i = g^{-1}(f(\mathbf{t}_i)), link function g(\cdot), and dispersion parameter \sigma^2. For Gaussian response with identity link, observations are independently distributed as y_i \mid \mathbf{t}_i, \sigma^2 \sim \mathcal{N}(f(\mathbf{t}_i), \sigma^2).

The objective is to estimate f by:

1. Partitioning the predictor space into K+1 mutually exclusive regions.

2. Fitting local polynomial models within each partition.

3. Enforcing smoothness at partition boundaries via Lagrangian multipliers.

4. Penalizing the integrated squared second derivative to discourage roughness.

Unlike other smoothing spline formulations, no post-fitting algebraic rearrangement or disentanglement of a spline basis is needed to obtain interpretable models. The polynomial expansions are homogeneous across partitions, and the relationship between predictor and response is explicit at the coefficient level.

To anchor the notation, in the single-predictor cubic case one would write

\hat{f}(t_i) = \hat{\beta}_{(0)} + \hat{\beta}_{(1)}t_i + \hat{\beta}_{(2)}t_i^2 + \hat{\beta}_{(3)}t_i^3 = \mathbf{x}_i^{\top}\hat{\boldsymbol{\beta}},

where \mathbf{x}_i = (1, t_i, t_i^2, t_i^3)^{\top}. The LMSS formulation preserves exactly this kind of polynomial representation, but now does so within each partition and then forces neighboring pieces to agree in the smoothness conditions described below.

Core notation used throughout:

• \mathbf{y}_{(N \times 1)}: Response vector.

• \mathbf{T}_{(N \times q)}: Matrix of predictors.

• \mathbf{X}_{(N \times P)}: Block-diagonal matrix of polynomial expansions, with diagonal blocks \mathbf{X}_k of dimension n_k \times p.

• \boldsymbol{\Lambda}_{(P \times P)}: Block-diagonal penalty matrix, with blocks \boldsymbol{\Lambda}_k of dimension p \times p.

• \hat{\boldsymbol{\beta}}_{(P \times 1)}: Unconstrained penalized estimate.

• \tilde{\boldsymbol{\beta}}_{(P \times 1)}: Constrained coefficient estimates.

• \mathbf{G}_{(P \times P)}: Block-diagonal matrix with blocks \mathbf{G}_k = (\mathbf{X}_k^{\top}\mathbf{W}_k\mathbf{D}_k\mathbf{X}_k + \boldsymbol{\Lambda}_k)^{-1}, where \mathbf{W}_k and \mathbf{D}_k are defined below.

• \mathbf{A}_{(P \times r)}: Constraint matrix encoding smoothness conditions. Reduced to linearly independent columns via pivoted QR decomposition.

• \mathbf{U}_{(P \times P)}: \mathbf{I} - \mathbf{G}\mathbf{A}(\mathbf{A}^{\top}\mathbf{G}\mathbf{A})^{-1}\mathbf{A}^{\top}.

• \mathbf{D}_{(N \times N)}: Diagonal matrix of user-supplied observation weights (observation_weights or weights). Defaults to the identity. These play the role of prior precision on individual observations: a weight of 2 is equivalent to seeing that observation twice.

• \mathbf{W}_{(N \times N)}: Diagonal matrix of GLM working weights. In the implementation these diagonal entries are whatever is returned by glm_weight_function; by default this is family$variance(mu), optionally multiplied by user-supplied observation weights. For Gaussian response with identity link, \mathbf{W} = \mathbf{I}. For other families, \mathbf{W} depends on the current fitted values and is updated at each Newton-Raphson iteration. For the common canonical families used by default, this matches the familiar Fisher-scoring weighting role.

• \mathbf{V}_{(N \times N)}: Correlation matrix of errors. When no correlation structure is specified, \mathbf{V} = \mathbf{I}. Otherwise supplied via VhalfInv or estimated through VhalfInv_fxn.

In the Gaussian identity case with unit weights and no correlation, \mathbf{G}_k = (\mathbf{X}_k^{\top}\mathbf{X}_k + \boldsymbol{\Lambda}_k)^{-1} and most formulas simplify accordingly. When \mathbf{D} or \mathbf{W} appear in a formula, the product \mathbf{W}\mathbf{D} means “GLM working weights times observation weights”; whenever one of them is the identity it drops out.

Before these quantities reach the main fitting stage, the user-facing inputs are parsed, standardized, and organized by process_input. When the formula interface is used and auto_encode_factors = TRUE, that preprocessing step also relies on helpers such as create_onehot to encode factor levels before the design reaches lgspline.fit(). The notation in the remainder of this document therefore refers to the internal objects that actually enter lgspline.fit(), not necessarily the raw objects originally supplied by the user.

From the user side, many of the arguments that control these internal objects can be supplied either individually or through the grouped lists penalty_args, tuning_args, expansion_args, constraint_args, qp_args, parallel_args, covariance_args, return_args, and glm_args, as documented in lgspline. These grouped lists are unpacked before dispatch into the same fitting pipeline, so they are a convenience layer rather than a separate modeling abstraction. A closely related exploratory mode is dummy_fit = TRUE in lgspline or lgspline.fit, which runs the preprocessing, partition construction, expansion building, and penalty setup without solving for nonzero coefficients, making it a practical way to inspect objects such as X, A, the returned make_partition_list from make_partitions, and the assembled penalties from compute_Lambda before a full fit.

Model Formulation and Estimation

Piecewise Polynomial Structure

For K knots (one predictor) or K+1 partitions (multiple predictors) there are K+1 mutually exclusive partitions \mathcal{P}_0, \dots, \mathcal{P}_{K}. Each observation i belongs to exactly one partition. Within partition k, the function is represented as a polynomial of degree p-1 in each predictor:

\hat{f}_k(\mathbf{t}) = \mathbf{x}^{\top}\tilde{\boldsymbol{\beta}}_k

where \mathbf{x} collects the polynomial basis terms (intercept, linear, quadratic, cubic, and optionally quartic and interaction terms) and \tilde{\boldsymbol{\beta}}_k are the corresponding coefficients. In one predictor, the same idea can be written more explicitly as

\hat{f}(t_i) = \sum_{k=0}^{K} \mathbf{x}_{ik}^{\top}\hat{\boldsymbol{\beta}}_k \mathbf{1}(t_i \in \mathcal{P}_k),

which highlights that the unconstrained problem is just a collection of local polynomial regressions. The expansions are homogeneous across partitions, so coefficients are directly comparable. This is implemented via get_polynomial_expansions.

The exact contents of \mathbf{x} are controlled by the basis-expansion arguments documented in lgspline: include_quadratic_terms, include_cubic_terms, include_quartic_terms, include_2way_interactions, include_3way_interactions, include_quadratic_interactions, exclude_interactions_for, exclude_these_expansions, and custom_basis_fxn. Likewise, just_linear_with_interactions and just_linear_without_interactions determine which predictors remain structurally linear even though they still participate in the same partition-wise polynomial bookkeeping described here.

Letting p denote the number of basis terms per partition, P = p(K+1) is the total number of coefficients. The full design matrix \mathbf{X} and penalty matrix \boldsymbol{\Lambda} are block-diagonal with K+1 blocks, so unconstrained estimation reduces to K+1 independent penalized regressions, which appears as follows for the identity link case:

\hat{\boldsymbol{\beta}}_k = \mathbf{G}_k \mathbf{X}_k^{\top} \mathbf{W}_k\mathbf{D}_k\mathbf{y}_k, \quad \mathbf{G}_k = (\mathbf{X}_k^{\top}\mathbf{W}_k\mathbf{D}_k\mathbf{X}_k + \boldsymbol{\Lambda}_k)^{-1}.

For Gaussian identity with unit weights this reduces to the familiar \mathbf{G}_k = (\mathbf{X}_k^{\top}\mathbf{X}_k + \boldsymbol{\Lambda}_k)^{-1}. The block-diagonal structure means these can be computed in parallel across partitions. In the user-facing interface this is realized by supplying a cluster through cl, optionally controlling work splitting with chunk_size, and enabling stages such as parallel_eigen for the eigendecompositions and, in non-Gaussian Path 3, parallel_unconstrained for the partition-wise unconstrained fits; nearby stages can likewise use parallel_penalty and parallel_make_constraint. The eigenvalue decomposition and matrix square roots of each \mathbf{G}_k are computed by compute_G_eigen, and can be returned in the fitted object as G and Ghalf when return_G = TRUE and return_Ghalf = TRUE.

Fitted values for the canonical Gaussian case appear as \tilde{\mathbf{y}} = \mathbf{X}\tilde{\boldsymbol{\beta}} = \mathbf{H}\mathbf{y} for \mathbf{H} = \mathbf{X}\mathbf{U}\mathbf{G}\mathbf{X}^{\top}.

Smoothing Constraints and the Constraint Matrix

Without further intervention the piecewise polynomial will be discontinuous. The central idea of LMSS is that smoothness is not hidden inside a special basis, but instead imposed directly where neighboring partitions meet. At each knot t_{k,k+1} between neighboring partitions k and k+1, up to three smoothing constraints are imposed:

1. Continuity: \mathbf{x}_{k,k+1}^{\top}\boldsymbol{\beta}_k = \mathbf{x}_{k,k+1}^{\top}\boldsymbol{\beta}_{k+1}.

2. First-derivative continuity: \mathbf{x}_{k,k+1}^{\prime\top}\boldsymbol{\beta}_k = \mathbf{x}_{k,k+1}^{\prime\top}\boldsymbol{\beta}_{k+1}.

3. Second-derivative continuity: \mathbf{x}_{k,k+1}^{\prime\prime\top}\boldsymbol{\beta}_k = \mathbf{x}_{k,k+1}^{\prime\prime\top}\boldsymbol{\beta}_{k+1}.

where \mathbf{x}^{\prime} and \mathbf{x}^{\prime\prime} are elementwise first and second derivatives of the basis with respect to \mathbf{t}. For the familiar cubic single-predictor basis \mathbf{x} = (1, t, t^2, t^3)^{\top}, these derivative vectors are

\mathbf{x}' = (0, 1, 2t, 3t^2)^{\top}, \qquad \mathbf{x}'' = (0, 0, 2, 6t)^{\top}.

With K knots this yields up to 3K scalar constraints (for a single predictor; more for multiple predictors with interactions), collected as linear equations \mathbf{A}^{\top}\boldsymbol{\beta} = \mathbf{0} in a P \times r matrix \mathbf{A}. The constraint matrix is built by make_constraint_matrix and returned in the fitted object as A.

In higher dimensions or with many partitions, the constraints can become over-specified and force the model toward a single global polynomial. In these cases it is recommended to drop second-derivative constraints or include quartic terms, allowing the model to fit a richer surface while maintaining perceived smoothness at knots. The appropriate constraint level can be controlled via include_constrain_fitted, include_constrain_first_deriv, and include_constrain_second_deriv. The companion flag include_constrain_interactions determines whether the analogous mixed-partial constraints are imposed for interaction terms, and no_intercept adds the special homogeneous equality constraint that fixes the intercept at zero (the same behavior triggered by using 0 + in the formula interface).

Before computing the projection \mathbf{U}, the constraint matrix is reduced to a linearly independent subset of columns via pivoted QR decomposition. This avoids numerical instability from redundant constraints and ensures \mathbf{A}^{\top}\mathbf{G}\mathbf{A} is invertible.

Lagrangian Projection

The constrained estimate is derived via Lagrangian multipliers. Define the P \times P projection matrix:

\mathbf{U} = \mathbf{I} - \mathbf{G}\mathbf{A}(\mathbf{A}^{\top}\mathbf{G}\mathbf{A})^{-1}\mathbf{A}^{\top}.

Then the constrained estimate is:

\tilde{\boldsymbol{\beta}} = \mathbf{U}\hat{\boldsymbol{\beta}}.

The matrix \mathbf{U} has the property that \mathbf{U}\mathbf{G}\mathbf{U}^{\top} = \mathbf{U}\mathbf{G}, which is used extensively in variance estimation and posterior draws. In words, the unconstrained penalized estimate is projected back into the coefficient space that satisfies the smoothness restrictions, and all subsequent uncertainty calculations inherit that same projected geometry. The projection is computed via get_U and, when requested, returned in the fitted object as U through return_U = TRUE.

When the constraints are inhomogeneous (\mathbf{A}^{\top}\boldsymbol{\beta} = \mathbf{c} with \mathbf{c} \neq \mathbf{0}), a particular solution \boldsymbol{\beta}_0 satisfying \mathbf{A}^{\top}\boldsymbol{\beta}_0 = \mathbf{c} is added back after projection, yielding the full Lagrangian solution \mathbf{U}\hat{\boldsymbol{\beta}} + (\mathbf{I} - \mathbf{U})\boldsymbol{\beta}_0. In lgspline and lgspline.fit, users realize this by supplying extra equality columns in constraint_vectors together with matching right-hand sides in constraint_values; null_constraint provides the alternate shorthand documented in lgspline when constraint_vectors is supplied and constraint_values is left empty.

In practice \mathbf{U} is never explicitly formed during fitting. The constrained estimate is obtained from a transformed OLS residual problem (the \mathbf{G}^{1/2}\mathbf{r}^{*} trick) in four steps:

1. Obtain the unconstrained partition-wise unconstrained estimate \hat{\boldsymbol{\beta}}.

2. Set \mathbf{y}^{*} = \mathbf{G}^{-1/2}\hat{\boldsymbol{\beta}} and \mathbf{X}^{*} = \mathbf{G}^{1/2}\mathbf{A}.

3. Fit the linear model \mathbb{E}[\mathbf{y}^{*}] = \mathbf{X}^{*}\boldsymbol{\gamma} by OLS using QR decomposition.

4. Compute the residuals \mathbf{r}^{*} = \mathbf{y}^{*} - \mathbf{X}^{*}(\mathbf{X}^{*\top}\mathbf{X}^{*})^{-1}\mathbf{X}^{*\top}\mathbf{y}^{*} from that transformed OLS fit and recover the constrained estimate by \tilde{\boldsymbol{\beta}} = \mathbf{G}^{1/2}\mathbf{r}^{*}.

A scaling factor 1/\sqrt{K+1} is applied to both \mathbf{X}^{*} and \mathbf{y}^{*} prior to the OLS call and divided out afterward, improving numerical conditioning when the constraint matrix has many rows.

The most expensive operation in this approach is the QR decomposition of the P \times r matrix \mathbf{X}^{*} = \mathbf{G}^{1/2}\mathbf{A}, which is far cheaper than working with the full P \times P system directly. Without correlation or SQP constraints, \mathbf{G} is stored and operated upon as a list of K+1 small p \times p matrices rather than the full P \times P block-diagonal, saving substantial memory when K is large and allowing for parallelism.

When correlation is present, \mathbf{X}^{\top}\mathbf{V}^{-1}\mathbf{X} is no longer block-diagonal, so the full-dimensional system must be handled directly unless the Woodbury acceleration (see .woodbury_decompose_V) applies. When additional inequality constraints are present, the code either augments the equality system with a partition-wise active-set refinement (block-separable case) or falls back to dense SQP via the Goldfarb-Idnani dual active-set method implemented in solve.QP.

GLM Extension and Iterative Updates

Working Quantities

For GLM responses with mean \mu_i = g^{-1}(\eta_i) and working weights w_i = [V(\mu_i)\{g'(\mu_i)\}^2]^{-1}, the penalized information becomes \mathbf{G}_k^{-1} = \mathbf{X}_k^{\top}\mathbf{W}_k\mathbf{X}_k + \boldsymbol{\Lambda}_k with \mathbf{W}_k = \mathrm{diag}(w_i : i \in \mathcal{P}_k). Because \mathbf{W} is diagonal, \mathbf{X}^{\top}\mathbf{W}\mathbf{X} remains block-diagonal and the four-step procedure carries through with \mathbf{G}_k = (\mathbf{X}_k^{\top}\mathbf{W}_k\mathbf{X}_k + \boldsymbol{\Lambda}_k)^{-1} in place of (\mathbf{X}_k^{\top}\mathbf{X}_k + \boldsymbol{\Lambda}_k)^{-1}. Under the default glm_weight_function, the diagonal entries reduce to family$variance(mu) (optionally scaled by observation weights), which for canonical families is the usual Fisher scoring weight.

The partition-wise unconstrained estimates are obtained by unconstrained_fit_fxn, by default unconstrained_fit_default, which initializes via an augmented ridge trick (appending \boldsymbol{\Lambda}^{1/2} as pseudo-observations to glm.fit) then refines using damped_newton_r with nr_iterate. For non-canonical log-link gamma regression, the keep_weighted_Lambda = TRUE option correctly returns the augmented ridge estimate directly.

Three Fitting Paths

Path 1: Correlation structure present. When a working correlation \mathbf{V} is present, such as for marginal means models or generalized estimating equation (GEE)-like models, \mathbf{V}^{-1} couples observations across partitions, so partition-wise fitting is not directly available. Two sub-paths handle this.

Path 1a (Gaussian identity + GEE) solves the whitened system directly via \tilde{\mathbf{G}} = (\tilde{\mathbf{X}}^{\top}\tilde{\mathbf{X}} + \boldsymbol{\Lambda}_{\mathrm{block}})^{-1} where \tilde{\mathbf{X}} = \mathbf{V}^{-1/2}\mathbf{X}, then applies the Lagrangian projection in the full P-space. See .get_B_gee_gaussian.

Path 1b (non-Gaussian GEE) uses a damped SQP iteration on the whitened system. The first iterate is a constrained Newton step from the projection matrix \mathbf{U}; subsequent iterates solve the quadratic subproblem via solve.QP. See .get_B_gee_glm.

Both sub-paths have Woodbury-accelerated variants described below.

Path 2: Gaussian identity, no correlation. The constrained estimate is obtained by a single Lagrangian projection from the per-partition penalized least-squares cross-products (the four steps in closed form). No outer iteration is needed. When K = 0 and there are no additional constraints, this reduces to the ordinary penalized closed form \hat{\boldsymbol{\beta}} = \mathbf{G}\mathbf{X}^{\top}\mathbf{y}. Implemented in .get_B_gaussian_nocorr.

Path 3: Non-Gaussian GLM, no correlation. Unconstrained estimates are obtained separately within each partition, then projected onto the constraint space. When the link is non-identity, the information \mathbf{G} depends on the current fitted values through the working weights, so the projection must be iterated. The unconstrained anchor \hat{\boldsymbol{\beta}} is held fixed while \mathbf{G} is recomputed at each iterate's constrained estimate:

\tilde{\boldsymbol{\beta}}^{(s+1)} = \mathbf{U}^{(s)}\hat{\boldsymbol{\beta}}, \quad \mathbf{U}^{(s)} = \mathbf{I} - \mathbf{G}^{(s)}\mathbf{A} (\mathbf{A}^{\top}\mathbf{G}^{(s)}\mathbf{A})^{-1}\mathbf{A}^{\top}.

Iteration stops when the mean absolute coefficient change falls below tol, or when the update begins increasing (in which case the previous iterate is restored). The recomputation of weighted Gram matrices, Schur corrections, and square-root information factors at each step is handled by .solver_recompute_G_at_estimate. Implemented in .get_B_glm_nocorr.

Woodbury Acceleration for Structured Correlation

For structured correlation matrices (AR(1), exchangeable, banded), the perturbation \mathbf{V}^{-1} - \mathbf{I} is sparse. Writing the GLS information matrix as

\mathbf{G}_V^{-1} = \underbrace{(\mathbf{X}^{\top}\mathbf{X} + \boldsymbol{\Lambda})}_{\text{block-diagonal}} + \underbrace{\mathbf{X}^{\top}(\mathbf{V}^{-1} - \mathbf{I})\mathbf{X}}_{\boldsymbol{\Delta}},

the block-diagonal part of \boldsymbol{\Delta} (within-partition corrections) is absorbed into per-partition eigendecompositions, yielding a corrected block-diagonal inverse \mathbf{G}_c. The off-diagonal remainder \boldsymbol{\Delta}_{\mathrm{off}} captures cross-partition coupling and has effective rank r determined numerically by .woodbury_decompose_V. For example, AR(1) correlation partitioned by time gives r \approx 2K since only observation pairs straddling knot boundaries contribute.

Factoring \boldsymbol{\Delta}_{\mathrm{off}} \approx \mathbf{L}\mathbf{S}\mathbf{L}^{\top} where \mathbf{L} is P \times r and \mathbf{S} is r \times r diagonal with entries \pm 1, the Woodbury identity gives

\mathbf{G}_V = \mathbf{G}_c - \mathbf{G}_c\mathbf{L} (\mathbf{S}^{-1} + \mathbf{L}^{\top}\mathbf{G}_c\mathbf{L})^{-1} \mathbf{L}^{\top}\mathbf{G}_c,

where the inner matrix is only r \times r.

The four-step projection is preserved by expressing \mathbf{G}_V^{1/2} through \mathbf{G}_c^{1/2}. Define \mathbf{Q} = \mathbf{G}_c^{1/2}\mathbf{L} (computed block-diagonally) with thin SVD \mathbf{Q} = \mathbf{O}_Q\boldsymbol{\Sigma}_Q\mathbf{R}_Q^{\top}. Then

\mathbf{G}_V^{1/2} = \mathbf{G}_c^{1/2}\mathbf{F}^{1/2}, \quad \mathbf{F}^{1/2} = \mathbf{I}_P - \mathbf{O}_Q\mathbf{C}\mathbf{O}_Q^{\top},

where \mathbf{C} is an r \times r diagonal matrix computed in .woodbury_halfsqrt_components. Every step of the projection decomposes as a block-diagonal operation through \mathbf{G}_c^{1/2} plus an additive rank-r correction through \mathbf{O}_Q:

\mathbf{y}^* = \mathbf{G}_c^{1/2}(\mathbf{X}^{\top} \mathbf{V}^{-1}\mathbf{y}) - \mathbf{G}_c^{1/2}\mathbf{O}_Q\mathbf{C}\mathbf{O}_Q^{\top} \mathbf{G}_c^{1/2}(\mathbf{X}^{\top}\mathbf{V}^{-1}\mathbf{y}),

and similarly for \mathbf{X}^* and the back-transformation. The OLS residual step itself is unchanged.

For the non-Gaussian Woodbury path (.get_B_gee_glm_woodbury), the perturbation \boldsymbol{\Delta}_V = \mathbf{V}^{-1} - \mathbf{I} and the product \boldsymbol{\Delta}_V\mathbf{X} are precomputed once and held fixed across Newton iterations. At each step, the weighted perturbation \boldsymbol{\Delta}(\mathbf{W}) = \mathbf{X}^{\top}\mathbf{W} (\mathbf{V}^{-1} - \mathbf{I})\mathbf{X} is formed using the precomputed product, then split and re-factored via .woodbury_redecompose_weighted.

The Woodbury path is used when r < P/3; otherwise the code falls back to the dense whitened approach. If the correction matrix \mathbf{F} is not positive definite at any point, the dense path is also used as fallback. See .get_B_gee_woodbury (Gaussian) and .get_B_gee_glm_woodbury (non-Gaussian).

Step Control

Different paths use different step-control strategies.

In the GEE paths (Paths 1a and 1b), the coefficient update is damped:

\boldsymbol{\beta}^{(s+1)} = (1 - \alpha_s)\boldsymbol{\beta}^{(s)} + \alpha_s\boldsymbol{\beta}_{\mathrm{cand}}^{(s)},

with \alpha_s = 2^{-d_s} and d_s increased whenever the candidate gives non-finite or larger deviance. The loop terminates after 10 consecutive rejections or 100 total iterations, and exits early when both coefficient change and deviance improvement fall below tol after a burn-in period.

In Path 3 (non-Gaussian, no correlation), the outer projection loop has no line search. The algorithm recomputes \mathbf{G} at the current constrained estimate and applies a fresh projection. If the mean absolute coefficient change begins increasing, the previous iterate is restored and the loop stops. The partition-wise unconstrained estimates themselves are obtained by damped Newton-Raphson inside unconstrained_fit_default, where damped_newton_r computes a Newton direction once per iteration and halves the step size until the penalized log-likelihood improves.

Accommodating Correlation Structures

Parametric Correlation Structures

Suppose \mathrm{Cov}(\mathbf{y}) = \sigma^2\,\mathbf{V}(\boldsymbol{\theta}) for a known parametric family indexed by \boldsymbol{\theta} (e.g., AR(1) with \theta = \rho, Matern with \boldsymbol{\theta} = (\ell, \nu), exchangeable with \theta = \rho). The penalized generalized least-squares problem becomes

\min_{\boldsymbol{\beta}}\; (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^{\top}\mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) + \boldsymbol{\beta}^{\top}\boldsymbol{\Lambda}\boldsymbol{\beta} \quad \text{s.t.}\; \mathbf{A}^{\top}\boldsymbol{\beta} = \mathbf{0}.

The correlation matrix \mathbf{V} is supplied through the fitted-object components Vhalf and VhalfInv, either directly or via user functions Vhalf_fxn and VhalfInv_fxn. When both are non-NULL, get_B dispatches to the GEE paths (Path 1a or 1b). For built-in correlation structures, the required square-root matrices are assembled numerically using matsqrt, matinvsqrt, and invert.

Whitening and Permutation

Because the data are stored in partition ordering (all observations from partition 0, then partition 1, etc.) while \mathbf{V} is in the original observation ordering, a permutation is applied internally: \mathbf{V}_{\mathrm{perm}}^{-1/2} = \mathbf{V}^{-1/2}[\boldsymbol{\pi}, \boldsymbol{\pi}], where \boldsymbol{\pi} = \texttt{unlist(order\_list)} maps original indices to partition-ordered indices. The whitened design and response are \tilde{\mathbf{X}} = \mathbf{V}_{\mathrm{perm}}^{-1/2}\mathbf{X}_{\mathrm{block}} and \tilde{\mathbf{y}} = \mathbf{V}_{\mathrm{perm}}^{-1/2}\mathbf{y}.

The \mathbf{X} and \mathbf{y} inputs to lgspline.fit are preserved in their unwhitened form. Whitening is applied inside get_B and blockfit_solve where the full N \times P block-diagonal design is available, since applying \mathbf{V}^{-1/2} to only the diagonal blocks of the partitioned design would silently discard cross-partition contributions and corrupt the Gram matrix.

Loss of Block-Diagonal Structure

Unlike the independent-errors case, \mathbf{X}^{\top}\mathbf{V}^{-1}\mathbf{X} is not block-diagonal because \mathbf{V}^{-1} introduces cross-partition coupling. The unconstrained estimator

\hat{\boldsymbol{\beta}} = (\mathbf{X}^{\top}\mathbf{V}^{-1}\mathbf{X} + \boldsymbol{\Lambda})^{-1} \mathbf{X}^{\top}\mathbf{V}^{-1}\mathbf{y}

requires a full P \times P solve, and the constraint projection proceeds with \mathbf{G} = (\mathbf{X}^{\top}\mathbf{V}^{-1}\mathbf{X} + \boldsymbol{\Lambda})^{-1}.

For structured correlation matrices (AR(1), exchangeable, banded), the perturbation \mathbf{V}^{-1} - \mathbf{I} is sparse and the cross-partition coupling has low effective rank. In these cases the Woodbury-accelerated paths (.get_B_gee_woodbury, .get_B_gee_glm_woodbury) recover the partition-wise computational structure by decomposing the coupling into a block-diagonal correction plus a low-rank remainder, as described in the GLM Extension section. For dense or high-rank \mathbf{V}^{-1}, the code falls back to the full whitened system (.get_B_gee_gaussian, .get_B_gee_glm).

GEE Deviance Monitoring

For non-Gaussian models with correlation, the deviance used for convergence monitoring is computed in the whitened space by .bf_gee_deviance. When the family supplies custom_dev.resids, the raw deviance residuals r_i = \mathrm{sign}(d_i)\sqrt{|d_i|} are divided by \sqrt{w_i} and pre=dmultiplied by \mathbf{V}_{\mathrm{perm}}^{-1/2} before squaring and averaging:

D_{\mathrm{GEE}} = \frac{1}{N}\left\| \mathbf{V}_{\mathrm{perm}}^{-1/2}\, \mathbf{w}^{-1/2}\mathbf{r}\right\|^{2},

where \mathbf{w} is the vector of working weights at the current iterate, clamped below at \sqrt{\varepsilon_{\mathrm{mach}}}. When only dev.resids is available, the function falls back to the standard mean deviance; otherwise it uses mean squared error.

REML Estimation of Correlation Parameters

Correlation parameters \boldsymbol{\theta} are estimated by minimizing a negative restricted log-likelihood (REML) objective. The criterion implemented in lgspline is a central-limit-theorem-based working approximation to a Laplace-style marginal likelihood criterion, applied here solely to correlation structure estimation rather than penalty parameter selection. Let \mathbf{D} = \mathrm{diag}(d_i) be the observation weight matrix, \tilde{\mathbf{W}} = \mathrm{diag}(\tilde{w}_i) the GLM working weight matrix at the current fitted values, \mathbf{V} the correlation matrix parameterized by \boldsymbol{\rho} (a vector on the unconstrained real line), and \tilde{\sigma}^{2} the dispersion profiled at its restricted maximum likelihood estimate. The negative REML objective implemented in lgspline, scaled by 1/N, is

-\ell_{R}(\boldsymbol{\rho}) = \frac{1}{N}\!\left[ -\log|\mathbf{V}^{-1/2}| + \frac{N}{2}\log\tilde{\sigma}^{2} + \frac{1}{2\tilde{\sigma}^{2}} (\mathbf{y} - \boldsymbol{\mu})^{\top} \mathbf{D}\tilde{\mathbf{W}}^{-1}\mathbf{V}^{-1} (\mathbf{y} - \boldsymbol{\mu}) + \frac{1}{2}\log|(\tilde{\sigma}^{2}\mathbf{U}\mathbf{G})^{-1}|^{+} \right],

where |\cdot|^{+} denotes the generalized determinant (product of nonzero eigenvalues), and \boldsymbol{\mu} = g^{-1}(\mathbf{X}\tilde{\boldsymbol{\beta}}) are the fitted values on the response scale.

Gradients with respect to correlation parameters are available in closed form for all built-in structures except Matern, which uses finite-difference approximation due to the complexity of differentiating the modified Bessel function K_{\nu} with respect to \nu. See reml_grad_from_dV for the full gradient derivation and notation. Custom analytic gradients can be supplied through REML_grad, and a fully custom criterion can replace REML through custom_VhalfInv_loss. The Toeplitz example in lgspline demonstrates how to supply custom correlation structures with user-defined gradient functions. Optimization over these working correlation parameters is then carried out by the same quasi-Newton engine used elsewhere in the package, namely efficient_bfgs with fallback to approx_grad when needed. In the user-facing interface, this machinery is activated through correlation_structure, correlation_id, spacetime, VhalfInv, Vhalf, VhalfInv_fxn, Vhalf_fxn, VhalfInv_par_init, REML_grad, custom_VhalfInv_loss, and VhalfInv_logdet.

The gradient of the negative REML has three terms per parameter:

1. \frac{1}{2}\mathrm{tr}(\mathbf{V}^{-1}\partial\mathbf{V}/\partial\theta_j): the log-determinant contribution.

2. -\frac{1}{2\tilde{\sigma}^{2}}\mathbf{r}^{\top} (\partial\mathbf{V}/\partial\theta_j)\mathbf{r}: the residual quadratic form contribution, where \mathbf{r} = \mathrm{diag}(\sqrt{d_i}/\sqrt{\tilde{w}_i}) \mathbf{V}^{-1/2}(\mathbf{y} - \boldsymbol{\mu}).

3. -\frac{1}{2}\mathrm{tr}\!\left(\mathbf{M}^{+}\mathbf{X}_{*}^{\top} \mathbf{V}^{-1}(\partial\mathbf{V}/\partial\theta_j)\mathbf{V}^{-1} \tilde{\mathbf{W}}\mathbf{D}\mathbf{X}_{*}\right): the REML correction, where \mathbf{X}_{*} = \mathbf{X}\mathbf{U} is the constrained design and \mathbf{M} = \mathbf{X}_{*}^{\top}\mathbf{V}^{-1} \tilde{\mathbf{W}}\mathbf{D}\mathbf{X}_{*} + \mathbf{U}^{\top}\boldsymbol{\Lambda}\mathbf{U} is the projected penalized information.

For each supported correlation family, the derivatives \partial\mathbf{V}/\partial\theta_j are available in closed form, enabling analytic gradient computation for use with the quasi-Newton optimizer.

Connection to Standard Mixed Model REML

The quadratic penalty \boldsymbol{\Lambda} acts as the inverse prior covariance of a Gaussian random effect on the spline coefficients, with the smoothing parameter satisfying \lambda = \tau^{2}/\sigma^{2}; this is the same mixed model representation used by mgcv, and Monte Carlo draws under the resulting Laplace-style approximation are available via generate_posterior. For Gaussian responses with identity link and no penalization, the implemented criterion coincides exactly with classical REML. For non-Gaussian responses, the criterion substitutes the Fisher information for the full penalized log-likelihood Hessian, exploiting the CLT approximation that \tilde{\mathbf{W}}^{-1/2}(\mathbf{y} - \boldsymbol{\mu}) is approximately Gaussian; this yields a method-of-moments style estimator for the correlation and variance parameters that depends only on mean-variance relationships through \mathbf{W}, and therefore generalizes naturally to quasi-likelihood and other settings where a fully specified log-likelihood is unavailable. The REML correction term \log|(\tilde{\sigma}^{2}\mathbf{U}\mathbf{G})^{-1}|^{+} uses a generalized log-determinant so that only nonzero eigenvalues contribute when rank deficiency arises from smoothness constraints or identifiability conditions in \mathbf{A}. When \boldsymbol{\Lambda} is full rank, the criterion coincides with the exact marginal likelihood from integrating out \boldsymbol{\beta} under its Gaussian prior; when \boldsymbol{\Lambda} is rank-deficient, unpenalized coefficients are projected out in the REML sense while penalized coefficients are integrated through their prior. During penalty tuning, the block-diagonal approximation is retained for GCV criteria and gradients; since GCV is rotation-invariant the practical effect on automatic selection is expected to be negligible, though this is not confirmed and the tuned penalties can always be overridden.

Built-In Correlation Structures

The package provides several built-in correlation structures for modeling spatial and temporal dependence. These are specified via correlation_structure with group membership in correlation_id and spatial or temporal coordinates in spacetime (an N-row matrix). Exchangeable correlation does not require spacetime.

All positive scale parameters are estimated on the log scale, with back-transform \exp(\cdot). Parameters constrained to (0, 1) use a double-exponential back-transform of the form \exp(-\exp(\eta)), so optimization still occurs on the unconstrained real line while the correlation remains bounded.

Exchangeable

Aliases: 'exchangeable', 'cs', 'CS', 'compoundsymmetric', 'compound-symmetric'.

A constant correlation \nu between any two observations within the same cluster. Parameterization: \nu = \exp(-\exp(\rho)), so \nu \in (0, 1). Only positive within-cluster correlation is supported under this parameterization.

Spatial Exponential

Aliases: 'spatial-exponential', 'spatialexponential', 'exp', 'exponential'.

Correlation decays exponentially with distance: \exp(-\omega d) where d is Euclidean distance and \omega > 0. Parameterization: \omega = \exp(\rho). Mathematically equivalent to the power correlation \theta^{d} with \theta = e^{-\omega}, but with better numerical properties during optimization.

AR(1)

Aliases: 'ar1', 'ar(1)', 'AR(1)', 'AR1'.

Correlation depends on rank difference between observations: \nu^{r} where r is the rank difference within cluster. Parameterization: \nu = \exp(-\exp(\rho)), so \nu \in (0, 1). Only positive autocorrelation is supported.

Gaussian / Squared Exponential

Aliases: 'gaussian', 'rbf', 'squared-exponential'.

Smooth decay with squared distance: \exp(-d^{2}/(2\ell^{2})) where \ell is the length scale. Parameterization: \ell = \exp(\rho).

Spherical

Aliases: 'spherical', 'Spherical', 'cubic', 'sphere'.

Polynomial decay with a hard cutoff at range r: 1 - 1.5(d/r) + 0.5(d/r)^{3} for d \le r, and 0 otherwise. Parameterization: r = \exp(\rho).

Matern

Aliases: 'matern', 'Matern'.

Flexible correlation with adjustable smoothness: (2^{1-\nu}/\Gamma(\nu))(\sqrt{2\nu}\,d/\ell)^{\nu}K_{\nu}(\sqrt{2\nu}\,d/\ell). Two parameters: length scale \ell = \exp(\rho_1) and smoothness \nu = \exp(\rho_2). No analytical gradient is available for \nu due to the difficulty of differentiating the modified Bessel function K_{\nu} with respect to its order, so finite differences are used; this makes Matern slower and potentially less stable than other structures.

Gamma-Cosine

Aliases: 'gamma-cosine', 'gammacosine', 'GammaCosine'.

Oscillatory dependence: (d^{\alpha-1}e^{-\gamma d})/(\Gamma(\alpha)/\gamma^{\alpha})\cdot\cos(\omega d). Three parameters: shape \alpha = \exp(\rho_1), rate \gamma = \exp(\rho_2), frequency \omega = \exp(\rho_3). Reduces to exponential when \alpha = 1 and \omega \approx 0.

Gaussian-Cosine

Aliases: 'gaussian-cosine', 'gaussiancosine', 'GaussianCosine'.

Smooth oscillatory correlation: \exp(-d^{2}/(2\ell^{2}))\cdot\cos(\omega d). Two parameters: length scale \ell = \exp(\rho_1) and frequency \omega = \exp(\rho_2). Reduces to Gaussian when \omega \approx 0.

Interpreting Estimated Correlation Parameters

Correlation parameters are estimated on transformed scales; they must be back-transformed for interpretation. When confint.lgspline is called and the inverse Hessian from BFGS is available, confidence intervals are returned on the untransformed (working) scale and should be back-transformed as described in the examples for lgspline.

Custom Correlation Structures

Custom correlation structures can be specified through:

• VhalfInv_fxn: Creates \mathbf{V}^{-1/2}.

• Vhalf_fxn: Creates \mathbf{V}^{1/2}. When omitted, the code computes it by explicit inversion of VhalfInv.

• REML_grad: Provides the analytical gradient of the REML objective.

• VhalfInv_logdet: Efficient log-determinant computation.

• custom_VhalfInv_loss: Replaces the REML objective entirely.

These functions enter lgspline through correlation_structure, VhalfInv_fxn, Vhalf_fxn, REML_grad, and custom_VhalfInv_loss, and the fitted object retains the resulting correlation machinery in components such as VhalfInv_fxn, Vhalf_fxn, and VhalfInv_params_estimates. When VhalfInv is supplied but Vhalf is not, Vhalf is computed unconditionally as the inverse of VhalfInv for all family/link combinations, since both get_B and blockfit_solve require it for GEE estimation.

Variance and Dispersion Estimation

Once the constrained estimate has been obtained, the next questions are how much flexibility the fitted model effectively used and how uncertainty should be propagated through the same constrained geometry. The quantities in this section are therefore all built on the projected information matrices from the previous sections.

Effective Degrees of Freedom and Dispersion

The effective degrees of freedom is the trace of the hat matrix. In the Gaussian identity case with observation weights and correlation, the fitted linear operator is built from the dense GLS analogue \mathbf{G}_{\mathrm{correct}} and can be written schematically as

\mathbf{H} = \mathbf{V}^{-1/2}(\mathbf{W}\mathbf{D})^{1/2} \mathbf{X}\mathbf{U}\mathbf{G}_{\mathrm{correct}}\mathbf{X}^{\top} (\mathbf{W}\mathbf{D})^{1/2}\mathbf{V}^{-1/2},

where for Gaussian identity \mathbf{W} = \mathbf{I}. In the no-correlation Gaussian case this reduces to the familiar \mathbf{H} = \mathbf{X}\mathbf{U}\mathbf{G}\mathbf{X}^{\top}\mathbf{D}.

For Gaussian identity fits, the dispersion estimate is computed as a weighted mean squared residual, optionally scaled by N/(N - \mathrm{tr}(\mathbf{H})) when unbias_dispersion = TRUE:

\tilde{\sigma}^{2} = \frac{1}{N - \mathrm{tr}(\mathbf{H})}\|\mathbf{y} - \mathbf{\tilde{y}}_i \|^{2}.

More generally with weights, a correlation structure and non-linear link function:

\tilde{\sigma}^{2} = \frac{1}{N - \mathrm{tr}(\mathbf{H})}\| \mathbf{V}^{-1/2}\mathbf{W}^{-1/2}\mathbf{D}^{1/2}(\mathbf{y} - \tilde{\mathbf{y}})\|^{2}.

This estimated dispersion is returned as sigmasq_tilde, and the corresponding effective degrees of freedom trace is returned as trace_XUGX. For non-Gaussian families, the fitting code delegates dispersion estimation to dispersion_function; thus the package does not assume a single closed-form Pearson-style formula outside the Gaussian identity setting. The hat-matrix trace itself is assembled by compute_trace_H in the dense correlation-aware case and by the same blockwise products summarized by trace_XUGX in the simpler no-correlation paths. A concrete built-in non-Gaussian example is the Weibull AFT path, which pairs weibull_family with weibull_dispersion_function, weibull_glm_weight_function, and weibull_schur_correction. Users who want these quantities available for downstream inference should keep estimate_dispersion = TRUE and return_varcovmat = TRUE (the defaults), since wald_univariate, confint.lgspline, and the prediction-standard-error path in predict.lgspline all rely on the post-fit dispersion and covariance components documented here.

Variance-Covariance Matrix

The variance-covariance matrix of \tilde{\boldsymbol{\beta}} is estimated as:

\mathrm{Var}(\tilde{\boldsymbol{\beta}}) = \tilde{\sigma}^{2}(\mathbf{U}\mathbf{G}^{1/2})(\mathbf{U}\mathbf{G}^{1/2})^{\top}

using the outer-product form for numerical stability. The result is returned as varcovmat when return_varcovmat = TRUE. The algebraically equivalent expression \tilde{\sigma}^{2}\mathbf{U}\mathbf{G} is not used because \mathbf{G} is only positive semi-definite when the penalty matrix \boldsymbol{\Lambda} has zero eigenvalues (e.g., the intercept and linear terms under the smoothing spline penalty when flat_ridge_penalty = 0), which can introduce negative diagonal entries in finite precision arithmetic. The outer-product form also guarantees symmetry.

This is the Bayesian posterior covariance, treating the penalty as a Gaussian prior on the coefficients. When exact_varcovmat = TRUE, a frequentist correction is additionally computed:

\mathrm{Var}_{\mathrm{exact}}(\tilde{\boldsymbol{\beta}}) = \tilde{\sigma^2}\mathbf{U}\mathbf{G}^{1/2}(\mathbf{X}^{\top}\mathbf{W}\mathbf{D}\mathbf{V}^{-1}\mathbf{X})\mathbf{G}^{1/2}\mathbf{U}^{\top} = \tilde{\sigma}^{2}\mathbf{U}\mathbf{G}\mathbf{U}^{\top} - \tilde{\sigma}^{2}\mathbf{U}\mathbf{G}\boldsymbol{\Lambda}\mathbf{G}\mathbf{U}^{\top}.

The first term is the Bayesian posterior covariance; the second is a frequentist correction such that for Gaussian identity link (with or without correlation), this is the exact variance-covariance matrix of the constrained estimator.

When a correlation structure is present (VhalfInv non-NULL), the block-diagonal \mathbf{G} is replaced by the full weighted GLS analogue

\mathbf{G}_{\mathrm{correct}} = \left(\mathbf{X}^{\top}\mathbf{W}\mathbf{D}\mathbf{V}^{-1}\mathbf{X} + \boldsymbol{\Lambda}\right)^{-1},

where \mathbf{W} = \mathbf{I} in the Gaussian identity case. This dense matrix is what enters the correlation-aware \mathbf{U}, \mathrm{Var}(\tilde{\boldsymbol{\beta}}), and \mathrm{Var}_{\mathrm{exact}}(\tilde{\boldsymbol{\beta}}) computations.

In user-facing terms, return_varcovmat controls whether this matrix is stored at all, while exact_varcovmat switches between the default posterior/Laplace approximation and the exact frequentist correction in the Gaussian-identity setting. The stored covariance is what powers wald_univariate, confint.lgspline, and se.fit = TRUE in predict.lgspline; the critical_value argument supplied at fit time is carried forward as the default cutoff for those interval-producing helpers.

Recomputation of G at Convergence

At the final iterate, \mathbf{G} is recomputed to reflect the converged working weights and Schur corrections. The implementation computes the weighted design \mathbf{X}_{w}^{(k)} = \mathbf{X}_k \cdot \mathrm{diag}(\sqrt{\mathbf{w}_k}), forms the weighted Gram matrix \mathbf{X}_{w}^{(k)\top}\mathbf{X}_{w}^{(k)}, adds the Schur correction, and performs eigendecomposition via compute_G_eigen to obtain \mathbf{G}_k, \mathbf{G}_k^{1/2}, and \mathbf{G}_k^{-1/2}. The relationship \mathbf{G}_k = \mathbf{G}_k^{1/2}(\mathbf{G}_k^{1/2})^{\top} is enforced exactly by construction, and the fitted object can retain these as G and Ghalf when return_G = TRUE and return_Ghalf = TRUE. The square-root factors are numerically stabilized through the helper routines matsqrt and matinvsqrt, which are also used elsewhere in the package when dense analogues of \mathbf{G}^{1/2} or \mathbf{G}^{-1/2} are required.

Bayesian Interpretation

The penalty has a natural Gaussian-prior interpretation, so once the constrained estimator and its covariance are available, Bayesian-style posterior simulation follows almost immediately. This section records the interpretation that is already implicit in the fitted object and in the package's posterior simulation helpers.

A Bayesian interpretation follows from viewing the penalty as a Gaussian prior on the coefficients. Conditional on the fitted smoothing parameters, the code samples on the coefficient scale from

\boldsymbol{\beta}^{(m)} = \tilde{\boldsymbol{\beta}} + \sqrt{\tilde{\sigma}^{2}},\mathbf{U}\mathbf{G}^{1/2}\mathbf{z}^{(m)}, \qquad \mathbf{z}^{(m)} \sim \mathcal{N}(\mathbf{0}, \mathbf{I}).

The coefficients are then back-transformed to the original response and predictor scales.

When inequality constraints are absent, these draws are i.i.d. Gaussian posterior draws around the fitted mode. The underlying coefficient-draw closure also contains an elliptical slice sampling route for active inequality constraints, using the same covariance factor to keep retained draws in the feasible region.

At the implementation level, standard Gaussian posterior draws may place positive mass on the infeasible region, so the constrained-draw closure instead targets the truncated posterior

\pi(\boldsymbol{\beta} \mid \mathbf{y}) \propto \exp\!\left(-\frac{1}{2}(\boldsymbol{\beta} - \tilde{\boldsymbol{\beta}})^{\top} \mathbf{G}^{-1}(\boldsymbol{\beta} - \tilde{\boldsymbol{\beta}})\right) \mathbf{1}(\mathbf{C}^{\top}\boldsymbol{\beta} \succeq \mathbf{c}),

yielding credible intervals that respect the constraint boundaries. The public generate_posterior wrapper forwards enforce_qp_constraints to the stored constrained-draw closure, so constrained draws can be requested directly from the user-facing interface. When a working correlation structure is present, the companion helper generate_posterior_correlation extends this idea by propagating uncertainty in the fitted correlation parameters through the same VhalfInv_fxn/Vhalf_fxn machinery described in the correlation section, rather than conditioning only on fixed covariance parameters. Correlation parameters are drawn from a multivariate normal distribution centered about their estimates with the inverse approxiamte BFGS Hessian of the REML optimization problem VhalfInv_params_vcov used by default (or a custom alternative as supplied to the argument correlation_param_vcov_sc).

Inequality Constraints via Sequential Quadratic Programming

Overview

Inequality constraints of the form \mathbf{C}^{\top}\boldsymbol{\beta} \succeq \mathbf{c} handle shape restrictions such as monotonicity, convexity, or boundedness. In the monomial basis, these are linear in \boldsymbol{\beta}. Monotonicity at a grid of points requires the first derivative polynomial to be non-negative there; convexity requires the second derivative to be non-negative; range constraints bound the fitted values directly. All of these translate to linear inequality constraints on the coefficient vector.

The inequality pieces are assembled by process_qp, which returns the qp_Amat, qp_bvec, and qp_meq objects passed to solve.QP, together with a quadprog flag indicating whether any inequality constraints are active. For derivative-sign constraints, process_qp calls .build_deriv_qp, which uses make_derivative_matrix on the expansion-standardized design and maps the derivative rows into the full P-dimensional coefficient space partition by partition.

Partition-Wise Active-Set Method

The sparsity pattern of \mathbf{C} is inspected automatically by .detect_qp_global (equivalently .solver_detect_qp_global). When every column of \mathbf{C} has nonzero entries in only a single partition block, the constraint system is block-separable and an active-set method can replace the dense QP.

At each iteration, the active set \mathcal{A} (constraints satisfied at equality) is appended to \mathbf{A} as additional equality constraints:

\mathbf{A}_{\mathrm{aug}} = [\mathbf{A} \mid \mathbf{C}_{\mathcal{A}}].

The constrained estimate is obtained by the same OLS projection with \mathbf{A}_{\mathrm{aug}} in place of \mathbf{A}, and since \mathbf{A}_{\mathrm{aug}} retains block-diagonal compatibility, all operations remain partition-wise.

The active set is updated by checking primal feasibility (adding the most-violated inactive constraint) and dual feasibility (dropping the active constraint with the most negative Lagrange multiplier) until the KKT conditions are satisfied. Lagrange multipliers for active inequalities are recovered from the OLS fit used in the Lagrangian projection: the fitted coefficients on \mathbf{X}^* = \mathbf{G}^{1/2}\mathbf{A}_{\mathrm{aug}} give the multipliers up to sign. Implemented in .active_set_refine and .check_kkt_partitionwise.

When any column of \mathbf{C} spans multiple partition blocks (for example, cross-knot monotonicity constraints), the block-diagonal structure is broken and the dense SQP approach is required. The selection is made automatically. If the active-set method does not converge within its iteration limit (default 50), the code falls back to dense SQP as well.

Dense SQP Iteration

The dense SQP approach, implemented in .qp_refine, solves a sequence of quadratic subproblems approximating the penalized log-likelihood. At each iteration s:

1. Compute the information matrix \mathbf{M}^{(s)} = \mathbf{X}^{\top}\mathbf{W}^{(s)}\mathbf{X} + \boldsymbol{\Lambda}_{\mathrm{block}} + \mathbf{S}^{(s)}, where \mathbf{S}^{(s)} is the Schur complement correction.

2. Compute the score vector via qp_score_function.

3. Solve the QP with solve.QP, where the combined constraint matrix is [\mathbf{A} \mid \mathbf{C}] with the first R columns treated as equalities.

4. Apply a damped update \boldsymbol{\beta}^{(s+1)} = (1 - \alpha)\boldsymbol{\beta}^{(s)} + \alpha\boldsymbol{\beta}_{\mathrm{QP}}, with \alpha = 2^{-d} and d incremented upon deviance increase.

A rescaling factor \mathrm{sc} = \sqrt{\mathrm{mean}(|\mathbf{M}^{(s)}|)} is applied to the Hessian and linear term before calling the QP solver. The equality-constrained estimate from the Lagrangian projection serves as a warm start. The qp_score_function defaults to the canonical GLM score \mathbf{X}^{\top}(\mathbf{y} - \boldsymbol{\mu}); for custom models a different score can be supplied.

Active Set and Lagrange Multipliers

The active set at the solution identifies binding inequality constraints. The Lagrange multipliers quantify the cost of each: a multiplier of zero means the constraint is not binding. The implementation stores active constraint indices, the corresponding submatrix of the constraint matrix, and the multiplier vector in the qp_info list returned alongside coefficient estimates, with components lagrangian, iact, and Amat_active. When return_lagrange_multipliers = TRUE, the fitted object also stores the final multiplier vector directly as lagrange_multipliers.

The original assembled inequality data are retained in the fitted object's quadprog_list component (containing qp_Amat, qp_bvec, and qp_meq from process_qp) so the final active set can be interpreted relative to the full specification. The final \mathbf{U} used in constructing the posterior variance-covariance matrix is built from both the equality constraints and the active inequality constraints at the solution.

Built-In Constraints

Built-in inequality constraints include:

• Monotonicity: qp_monotonic_increase, qp_monotonic_decrease. Enforced by requiring consecutive fitted values to be non-decreasing (or non-increasing): (\mathbf{x}_i - \mathbf{x}_{i-1})^{\top}\boldsymbol{\beta} \geq 0. These are constructed by process_qp from the partition-stacked block design reordered to observation order.

• Derivative sign: qp_positive_derivative, qp_negative_derivative. Enforced through the first-derivative design matrix from make_derivative_matrix. May be TRUE/FALSE or a character/integer vector selecting specific predictors.

• Second-derivative sign: qp_positive_2ndderivative, qp_negative_2ndderivative. Same construction using the second-derivative design matrix.

• Response range: qp_range_lower, qp_range_upper. Bounds on the linear predictor; for non-identity links, the bounds are transformed to the link scale inside process_qp.

• Custom constraints via qp_Amat_fxn, qp_bvec_fxn, and qp_meq_fxn, which receive the design matrix structure and return the constraint matrix, bound vector, and number of equalities. These are commonly paired with a custom qp_score_function when the quadratic approximation uses a non-default likelihood. The low-level objects qp_Amat, qp_bvec, and qp_meq remain documented in lgspline but in the current implementation serve as activation markers rather than being merged into the constraint set assembled by process_qp.

All constraints can be thinned to a user-specified subset of rows via qp_observations, which process_qp applies before assembly.

Blockfit Backfitting for Linear Non-Interactive Effects

Motivation

When a model contains both spline terms (receiving K+1 partition-specific coefficient vectors constrained to smoothness) and non-interactive linear terms (“flat” terms, specified via just_linear_without_interactions, receiving a single shared coefficient vector \mathbf{v} across all partitions), the standard solver carries K+1 copies of \mathbf{v} linked by equality constraints. Backfitting avoids this inflation by solving a lower-dimensional problem at each step. Write the partition-k design as \mathbf{X}_k = [\mathbf{Z}_k \mid \mathbf{X}_{\mathrm{flat}}^{(k)}], where \mathbf{Z}_k contains the spline columns and \mathbf{X}_{\mathrm{flat}}^{(k)} the flat columns. This is invoked when blockfit = TRUE, flat columns are non-empty, and K > 0.

The design, penalty, and constraint matrices are split into spline and flat components by .bf_split_components, which extracts spline rows from \mathbf{A}, drops null columns, rank-reduces via QR, and detects mixed constraints (columns of \mathbf{A} with nonzero entries on both spline and flat rows).

Block-Coordinate Descent

Spline step. Holding \mathbf{v} fixed, the code forms the adjusted response \mathbf{y}_k - \mathbf{X}_{\mathrm{flat}}^{(k)}\mathbf{v} and applies the spline-only Lagrangian projection via .bf_lagrangian_project:

\tilde{\boldsymbol{\beta}}_{\mathrm{spline}}^{(k)} = \mathbf{U}_{\mathrm{spline}}\mathbf{G}_{\mathrm{spline}} \mathbf{Z}_k^{\top}(\mathbf{y}_k - \mathbf{X}_{\mathrm{flat}}^{(k)}\mathbf{v}).

Flat step. Holding spline coefficients fixed, the shared flat vector is updated by pooled penalized regression on residuals:

\mathbf{v} = \left(\sum_{k=0}^{K} \mathbf{X}_{\mathrm{flat}}^{(k)\top} \mathbf{X}_{\mathrm{flat}}^{(k)} + \boldsymbol{\Lambda}_{\mathrm{flat}}\right)^{-1} \sum_{k=0}^{K} \mathbf{X}_{\mathrm{flat}}^{(k)\top} (\mathbf{y}_k - \mathbf{Z}_k \tilde{\boldsymbol{\beta}}_{\mathrm{spline}}^{(k)}).

When the constraint matrix \mathbf{A} has mixed columns (nonzero entries on both spline and flat rows), the flat update instead solves a KKT system enforcing the residual equality constraint \mathbf{A}_{\mathrm{flat}}^{\top}\mathbf{v} = \mathbf{c} - \mathbf{A}_{\mathrm{spline}}^{\top} \tilde{\boldsymbol{\beta}}_{\mathrm{spline}} via .bf_constrained_flat_update.

Convergence is checked using the maximum absolute change across both spline and flat coefficients. In the weighted inner loop used by the GLM solvers, the flat-block change alone determines stopping.

Four Estimation Cases

blockfit_solve dispatches to one of four paths.

Case (a): Gaussian identity + GEE (.bf_case_gauss_gee). Whitening destroys block-diagonal structure, so the code skips backfitting and performs the same full-system Gaussian GEE projection used by get_B Path 1a. The result is split back into spline and flat components for downstream assembly.

Case (b): Gaussian identity, no correlation (.bf_case_gauss_no_corr). Standard block-coordinate descent as described above. The spline-only \mathbf{G}_{\mathrm{spline}} factors are precomputed once, and the pooled flat penalized inverse is reused across iterations.

Case (c): GLM + GEE (.bf_case_glm_gee). Two stages. Stage 1 forms a warm start by running damped Newton-Raphson on the unwhitened working response: each outer iteration computes working responses and weights at the current linear predictor, then calls .bf_inner_weighted for the inner backfitting loop. Stage 2 refines the warm start on the full whitened system using the damped SQP loop (.bf_sqp_loop), replicating the approach used by .get_B_gee_glm.

Case (d): GLM without GEE (.bf_case_glm_no_corr). A damped Newton-Raphson outer loop updates working responses and weights, while each inner iteration alternates between weighted spline and weighted flat updates via .bf_inner_weighted. Deviance is monitored across outer iterations for convergence and damping.

Inequality Constraints and Reassembly

Because flat coefficients are shared by construction, the corresponding equality constraints are satisfied exactly. Smoothness constraints on the spline block are enforced by the spline-only Lagrangian projection. After backfitting convergence, inequality constraints are enforced via the same partition-wise active-set or dense SQP refinement used by get_B. The method is selected automatically by .solver_detect_qp_global: block-separable constraints use the active-set method through .active_set_refine, while cross-partition constraints trigger the dense SQP loop through .bf_sqp_loop. For GEE (Case c), inequality handling occurs inside Stage 2 on the whitened system.

After convergence, the shared flat vector \mathbf{v} is copied into each partition's coefficient vector, yielding \boldsymbol{\beta}_k = [\tilde{\boldsymbol{\beta}}_{\mathrm{spline}}^{(k)\top}, \mathbf{v}^{\top}]^{\top} for compatibility with downstream inference. If blockfit_solve throws an error, a warning is issued and the code falls back to get_B.

Knot Selection and Partitioning

The topic of knot selection is not the main focus of the package, but the partition structure is central because every later design matrix, penalty, and smoothness constraint depends on it. The defaults in lgspline are therefore meant to be practical and transparent rather than theoretically final.

Univariate Case

For a single predictor, the default partitioning is now handled by make_partitions in the same k-means framework used more generally: K+1 centers are fit on an internally standardized copy of the predictor, controlled by standardize_predictors_for_knots, and then returned on the raw scale. Custom knots can still be supplied via custom_knots, in which case partition assignment is built directly from those raw-scale breakpoints. The default number of knots K is chosen adaptively based on N, p, q, and the GLM family. For multivariate fits, the resulting partition metadata are returned as make_partition_list and can be re-used in later calls to lgspline. This is particularly useful when one wants to hold the partition geometry fixed across repeated fits, for example while varying penalties, families, or correlation structures.

Multivariate Case

For multiple predictors, K+1 cluster centers are identified by k-means on an internally standardized predictor matrix via make_partitions. This is the partitioning mechanism used to determine the multivariate spline regions; see MacQueen (1967) for the classical clustering formulation and Kisi et al. (2025) for a recent applied example of k-means-driven partitioning in a nonlinear prediction setting. Midpoints between neighboring centers (those whose midpoint does not fall into a third cluster) serve as knot locations. Observations are assigned to the nearest cluster center using get.knnx, and the returned centers and knots are on the original predictor scale. The resulting partition structure is a type of Voronoi diagram and is stored in the fitted object as make_partition_list. The do_not_cluster_on_these argument can exclude certain predictors from clustering (e.g., a treatment indicator that should not drive partitioning). The lower-level clustering behavior can be further controlled by cluster_args and neighbor_tolerance, while cluster_on_indicators determines whether binary predictors are allowed to influence the partition geometry at all.

Standardizing Predictors

Higher-order polynomial terms can dramatically inflate or deflate the magnitude of basis expansions, introducing numerical instability. All polynomial basis expansions are scaled by q_{0.69} - q_{0.31}, where q_{\zeta} is the \zeta-th quantile of the expansion. For a standard normal distribution this quantity is approximately 1, so the scaling is close to one standard deviation for symmetric distributions. This fitting-stage rescaling is controlled by standardize_expansions_for_fitting, while knot construction is controlled separately by standardize_predictors_for_knots. The same scaling is applied to the constraint matrix to maintain smoothness, and coefficients are back-transformed to the original scale after fitting.

Smoothing Spline Penalty

Penalty Construction

The penalty matrix \boldsymbol{\Lambda}_s penalizes the integrated squared total curvature of the fitted function over the observed predictor ranges. This is the step that makes the piecewise polynomial fit genuinely behave like a smoothing spline rather than merely a constrained regression spline. The package computes this penalty directly from the monomial structure of the basis rather than by appealing to a pre-tabulated spline basis. For a single partition k with basis expansion \mathbf{x}_k = (\phi_1(\mathbf{t}), \ldots, \phi_p(\mathbf{t}))^{\top} where each \phi_i(\mathbf{t}) = \prod_{j=1}^{q} t_j^{\alpha_{ij}} is a multivariate monomial:

\boldsymbol{\beta}_k^{\top}\boldsymbol{\Lambda}_s\boldsymbol{\beta}_k = \int_{\mathbf{a}}^{\mathbf{b}} \|\tilde{f}_k''(\mathbf{t})\|^{2}\,d\mathbf{t},

where \mathbf{a} and \mathbf{b} are the observed predictor minimums and maximums (computed globally from the data, not partition-specific), and \tilde{f}_k(\mathbf{t}) = \mathbf{x}_k^{\top}\boldsymbol{\beta}_k is the fitted function for partition \mathcal{P}_k.

Total curvature operator. The integrated squared second derivative decomposes into q curvature operators, one per predictor. For predictor v, the curvature operator D_v is defined as

D_v = \frac{\partial^{2}}{\partial t_v^{2}} + \sum_{s \neq v}\frac{\partial^{2}}{\partial t_v\,\partial t_s}.

That is, D_v captures both the pure second derivative with respect to t_v and all mixed second partial derivatives involving t_v. The penalty matrix entries are then

[\boldsymbol{\Lambda}_s]_{ij} = \sum_{v=1}^{q}\int_{\mathbf{a}}^{\mathbf{b}} D_v(\phi_i)\,D_v(\phi_j)\,d\mathbf{t}.

Monomial derivative rule. For a monomial \phi(\mathbf{t}) = \prod_j t_j^{\alpha_j}, the derivatives entering D_v have closed forms. The pure second derivative is

\frac{\partial^{2}}{\partial t_v^{2}}\prod_j t_j^{\alpha_j} = \alpha_v(\alpha_v - 1)\,t_v^{\alpha_v - 2}\prod_{j \neq v} t_j^{\alpha_j},

which is zero when \alpha_v < 2. The mixed second derivative is

\frac{\partial^{2}}{\partial t_v\,\partial t_s}\prod_j t_j^{\alpha_j} = \alpha_v\alpha_s\,t_v^{\alpha_v - 1}t_s^{\alpha_s - 1} \prod_{j \neq v,s} t_j^{\alpha_j},

which is zero when \alpha_v < 1 or \alpha_s < 1. Applying D_v to a monomial \phi_i produces a sum of monomials with known coefficients and exponent vectors.

Factorized integration. Because every D_v(\phi_i) is polynomial, the product D_v(\phi_i)\,D_v(\phi_j) is also polynomial and the multivariate integral factorizes over predictors:

\int_{\mathbf{a}}^{\mathbf{b}}\prod_{j=1}^{q} t_j^{e_j}\,d\mathbf{t} = \prod_{j=1}^{q}\frac{b_j^{e_j+1} - a_j^{e_j+1}}{e_j + 1}.

Crucially, this integral runs over all q predictor ranges, including predictors that do not appear in the integrand (for which e_j = 0 and the factor reduces to b_j - a_j). This ensures that the penalty is properly scaled relative to the volume of the predictor space.

Single-predictor verification. For q = 1 with expansion \mathbf{x} = (1, t, t^{2}, t^{3})^{\top} on [a, b], the curvature operator reduces to D_1 = \partial^{2}/\partial t^{2} (no mixed partials exist), and the penalty matrix reduces to

\boldsymbol{\Lambda}_s = \int_a^b \mathbf{x}''\mathbf{x}''^{\top}\,dt = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 4(b - a) & 6(b^{2} - a^{2}) \\ 0 & 0 & 6(b^{2} - a^{2}) & 12(b^{3} - a^{3}) \end{pmatrix},

equivalent to the classical cubic smoothing spline penalty originally proposed by Reinsch.

Handling of non-spline predictors. Predictors specified via just_linear_without_interactions or just_linear_with_interactions do not receive higher-order polynomial expansions in the design matrix. To ensure their curvature contributions are still correctly computed (particularly through interaction terms), the implementation temporarily appends phantom higher-order columns (with zero data) for these predictors, computes the full curvature penalty on the augmented basis, and then subsets the result back to the original p \times p dimensions. This ensures that interaction terms involving non-spline predictors receive appropriate penalty contributions without affecting the rest of the estimation pipeline.

Parallel computation. Because the total penalty is an additive sum over predictors (\boldsymbol{\Lambda}_s = \sum_{v=1}^{q}\boldsymbol{\Lambda}_{s,v}), the computation can be parallelized by distributing the per-predictor curvature matrices across workers via parallel::parLapply and summing the results. This is controlled by the parallel_penalty argument and is beneficial when q is large.

The penalty is computed by get_2ndDerivPenalty (single predictor or subset) and get_2ndDerivPenalty_wrapper (full assembly with optional parallelism and non-spline handling).

Because the smoothing penalty has zero eigenvalues for the intercept and linear terms (whose second derivatives vanish), an optional ridge penalty on lower-order terms is added for computational stability. The full penalty block for partition k is:

\boldsymbol{\Lambda}_k = \lambda_w\bigl(\boldsymbol{\Lambda}_s + \lambda_r\boldsymbol{\Lambda}_r + \sum_{m=1}^{M}\xi_{mk}\mathbf{L}_{mk}\bigr)

where \lambda_w is the global wiggle penalty (wiggle_penalty), \lambda_r is ridge penalty on linear and intercept terms (flat_ridge_penalty) multiplied by the wiggle penalty, and \xi_{mk} and \mathbf{L}_{mk} denote optional additional penalty multipliers and matrices, including the predictor- and partition-specific components activated through unique_penalty_per_predictor, unique_penalty_per_partition, predictor_penalties, and partition_penalties. This assembly is handled by compute_Lambda.

The penalty matrix \boldsymbol{\Lambda} is stored as a list of K+1 p \times p square, symmetric, positive semi-definite matrices.

Penalty Optimization via Generalized Cross-Validation

After the structural pieces of the model are fixed, the main remaining question is how much smoothing to apply. In lgspline, that tuning is performed with generalized cross-validation, but it is carried out using the same constrained estimator that will be used in the final model fit. Penalty parameters are estimated on the log scale via exponential parameterization (\lambda = \exp(\theta), \theta \in \mathbb{R}), ensuring positivity. The chain rule factor \partial\exp(\theta)/\partial\theta = \exp(\theta) = \lambda is applied throughout. User-facing arguments (initial_wiggle, initial_flat, predictor_penalties, partition_penalties) accept values on the raw, natural scale; conversion to log scale is handled internally. The final tuned values and assembled penalty pieces are returned in the fitted object's penalties component.

The total penalty matrix \boldsymbol{\Lambda} is constructed as:

\boldsymbol{\Lambda} = \lambda_w\mathbf{L}_{w} + \lambda_r\mathbf{L}_{r} + \sum_{j}\nu_j\mathbf{L}_j^{(\mathrm{pred})} + \sum_{k}\tau_k\mathbf{L}_k^{(\mathrm{part})},

where \mathbf{L}_{w} is the integrated squared second-derivative penalty (i.e., \boldsymbol{\Lambda}_s above), \mathbf{L}_{r} is a ridge penalty on intercept and linear coefficients, \mathbf{L}_j^{(\mathrm{pred})} are predictor-specific penalties, and \mathbf{L}_k^{(\mathrm{part})} are partition-specific penalties. The scalars \lambda_w (wiggle_penalty), \lambda_r (flat_ridge_penalty), \{\nu_j\} (predictor_penalties), and \{\tau_k\} (partition_penalties) are tuned. The unbiased generalized cross-validation criterion is

\mathrm{GCV}_{u} = \frac{\sum_{i=1}^{N}D_{ii}\,r_i^{2}}{N(1 - \bar{W})^{2}},

where r_i are residuals on the link scale and \bar{W} = \mathrm{tr}(\mathbf{H})/N is the mean of the hat-matrix diagonal. For identity link, r_i = y_i - \hat{\eta}_i. For non-identity links, the residuals are r_i = g((y_i + \delta)/(1+2\delta)) - (\hat{\eta}_i + \delta)/(1+2\delta), where \delta \geq 0 is a pseudocount that stabilizes the link transformation, automatically tuned within tune_Lambda if not supplied to delta. Non-Gaussian families with observation weights \omega_i have their residuals scaled by \omega_i. When the family provides a custom deviance residual function, that function is used in place of the link-scale residuals.

Pseudocount selection. The pseudocount \delta is chosen to make the transformed response distribution most closely approximate a t-distribution with N-1 degrees of freedom, in the sense of minimizing the (optionally weighted) mean absolute deviation between the sorted standardized transformed responses and the corresponding t-quantiles. This is solved via Brent's method over [10^{-64}, 1]. When the link is identity, or when the response naturally lies in the domain of the link function, \delta = 0. This behavior is exposed through the delta argument in lgspline: supplying a fixed numeric value bypasses the internal search, while leaving it NULL allows the tuning code to choose the stabilizing pseudocount automatically when needed.

Meta-penalty regularization. A regularization term pulls the predictor- and partition-specific penalty parameters toward 1 on the raw scale:

P_{\mathrm{meta}}(\lambda_w, \nu_j, \tau_k) = \frac{1}{2}c_{\mathrm{meta}}\sum_j(\nu_j - 1)^{2} + \frac{1}{2}\cdot 10^{-32}(\lambda_w - 1)^{2},

where c_{\mathrm{meta}} is a user-specified coefficient (meta_penalty). The gradient of P_{\mathrm{meta}} on the log scale, incorporating the exp chain rule, is \partial P_{\mathrm{meta}}/\partial\theta_j = c_{\mathrm{meta}}(\nu_j - 1)\nu_j and \partial P_{\mathrm{meta}}/\partial\theta_1 = 10^{-32}(\lambda_w - 1)\lambda_w. The total objective is \mathrm{GCV}_{u} + P_{\mathrm{meta}}.

Closed-Form Gradient of GCV

The gradient of \mathrm{GCV}_{u} with respect to \theta_1 = \log\lambda_w is computed analytically via the quotient rule:

\frac{\partial\mathrm{GCV}_{u}}{\partial\theta_1} = \frac{1}{D^{2}}\left(\frac{\partial\mathcal{N}}{\partial\theta_1}D - \mathcal{N}\frac{\partial D}{\partial\theta_1}\right),

where \mathcal{N} = \sum r_i^{2} (numerator) and D = N(1 - \bar{W})^{2} (denominator). The key intermediates are:

• \partial\mathbf{G}/\partial\lambda_w, computed from the matrix identity \partial(\mathbf{X}^{\top}\mathbf{X} + \boldsymbol{\Lambda})^{-1}/\partial\lambda = -\mathbf{G}(\partial\boldsymbol{\Lambda}/\partial\lambda)\mathbf{G}.

• \partial\mathbf{G}^{1/2}/\partial\lambda_w, derived from \partial\mathbf{G}/\partial\lambda_w via the eigendecomposition chain rule.

• \partial\bar{W}/\partial\lambda_w, the derivative of the trace of the hat matrix \mathbf{H} = \mathbf{X}\mathbf{U}\mathbf{G}\mathbf{X}^{\top}, which depends on both \partial\mathbf{G}/\partial\lambda_w and \partial\mathbf{G}^{1/2}/\partial\lambda_w.

• \partial\mathcal{N}/\partial\theta_1 = -2\mathbf{r}^{\top}\mathbf{X}(\partial(\mathbf{U}\mathbf{G})/\partial\lambda_w)\mathbf{X}^{\top}\mathbf{y}\cdot\lambda_w, via the chain rule applied to the residual vector.

• \partial D/\partial\theta_1 = 2(1 - \bar{W})(-\partial\bar{W}/\partial\lambda_w)\cdot\lambda_w.

In the implementation, these quantities are assembled by a small set of helper routines: compute_dG_dlambda for \partial\mathbf{G}/\partial\lambda, compute_dGhalf for \partial\mathbf{G}^{1/2}/\partial\lambda, compute_dW_dlambda_wrapper for derivatives of the effective degrees-of-freedom term, compute_trace_UGXX_wrapper for the trace pieces entering GCV, and compute_dG_u_dlambda_xy for the derivative of the fitted-value quadratic form.

The full gradient is scaled by N before adding the meta-penalty gradient.

For the ridge penalty and predictor-/partition-specific penalties, a trace-ratio heuristic is used:

\frac{\partial\mathrm{GCV}_{u}}{\partial\lambda_l} \approx \frac{\mathrm{mean}(\mathrm{diag}(\mathbf{L}_l))}{\mathrm{mean}(\mathrm{diag}(\boldsymbol{\Lambda}))} \frac{\partial\mathrm{GCV}_{u}}{\partial\lambda_w},

and analogously for predictor- and partition-specific penalties, where \mathbf{L}_j^{(\mathrm{pred})} or \mathbf{L}_k^{(\mathrm{part})} replaces \mathbf{L}_l in the numerator. This follows from a chain-rule argument: by the Leibniz rule and the inverse derivative, \partial\lambda_w/\partial\lambda_l = (\partial\boldsymbol{\Lambda}/\partial\lambda_l)(\partial\boldsymbol{\Lambda}/\partial\lambda_w)^{-1}. Since the derivative appears as a matrix rather than a scalar, the mean-diagonal ratio provides a scalar summary. Once the derivative for \lambda_w is in hand, the derivatives of other penalties are cheap to compute. The exp chain rule is then applied: \partial/\partial\theta = (\partial/\partial\lambda)\cdot\lambda.

Optimization Procedure

Grid search initialization. The \mathrm{GCV}_{u} criterion is evaluated over a grid of candidate values for (\lambda_w, \lambda_r) on the log scale. All combinations of user-supplied candidate vectors (initial_wiggle and initial_flat) are formed, and the combination yielding the smallest finite \mathrm{GCV}_{u} is selected as the starting point for BFGS optimization. Grid points producing non-finite \mathrm{GCV}_{u} are discarded. If all grid points fail, an error is raised advising the user to check the data or adjust the grid.

Damped BFGS optimizer. A custom damped BFGS quasi-Newton optimizer, implemented in efficient_bfgs, minimizes \mathrm{GCV}_{u} + P_{\mathrm{meta}}. When analytic gradients are not usable, the fallback finite-difference helper is approx_grad.

Iterations 1-2: steepest descent. The first two iterations use steepest descent with a damping factor \alpha: \boldsymbol{\phi}^{(t+1)} = \boldsymbol{\phi}^{(t)} - \alpha\nabla_{\boldsymbol{\phi}}.

Iterations 3+: BFGS. From iteration 3, an inverse Hessian approximation \mathbf{J}^{(t)} is maintained via the standard secant update. Let \mathbf{s}^{(t)} = \boldsymbol{\phi}^{(t)} - \boldsymbol{\phi}^{(t-1)} and \mathbf{v}^{(t)} = \nabla^{(t)} - \nabla^{(t-1)}. The BFGS update is:

\mathbf{J}^{(t+1)} = (\mathbf{I} - \mathbf{u}\mathbf{s}\mathbf{v}^{\top}) \mathbf{J}^{(t)} (\mathbf{I} - \mathbf{u}\mathbf{v}\mathbf{s}^{\top}) + \mathbf{u}\mathbf{s}\mathbf{s}^{\top}, \qquad \mathbf{u} = (\mathbf{v}^{\top}\mathbf{s})^{-1}.

When |\mathbf{v}^{\top}\mathbf{s}| < 10^{-64}, the approximation is reset to \mathbf{I} and the iteration is flagged for restart. The search direction is \mathbf{d}^{(t)} = -\mathbf{J}^{(t)}\nabla^{(t)}.

Step acceptance. A step is accepted if \mathrm{GCV}_{u}^{(\mathrm{new})} \leq \mathrm{GCV}_{u}^{(\mathrm{old})}. On rejection, \alpha is halved. If \alpha < 2^{-10} (early iterations) or \alpha < 2^{-12} (later iterations), the optimizer terminates with the best solution found.

Convergence. The optimizer terminates when |\mathrm{GCV}_{u}^{(t)} - \mathrm{GCV}_{u}^{(t-1)}| < \epsilon or \|\boldsymbol{\phi}^{(t)} - \boldsymbol{\phi}^{(t-1)}\|_{\infty} < \epsilon, provided at least 10 iterations have elapsed, for penalties \boldsymbol{\phi} = \lambda_w, \lambda_r, \nu_j, \tau_j, ....

Alternative. A base-R optim call with method "BFGS" and finite-difference gradients can be used used instead via use_custom_bfgs=FALSE, which uses stats::optim.

Post-optimization inflation. After optimization, the penalty parameters are inflated by a factor ((N+2)/(N-2))^{2} to counteract the in-sample bias toward underpenalization inherent in GCV-type criteria.

The tuning loop is implemented in tune_Lambda.

Incorporating Non-Spline Effects

Multiple fixed effects are accommodated naturally in the LMSS framework because spline effects, linear effects, and many interaction terms all live in the same partition-wise polynomial expansion. The distinction is therefore not whether a term is "allowed" by the solver, but whether it receives full spline treatment or remains structurally linear across partitions.

The constrained framework naturally accommodates non-spline terms. If only linear terms are included for a predictor (via just_linear_without_interactions or just_linear_with_interactions), the first-derivative smoothing constraint forces the linear coefficient to be identical across all partitions, since the derivative of a linear function is its slope. This is not an algorithmic modification but a natural consequence of the constraint structure.

For example, a model with one spline effect and a linear treatment indicator interaction will naturally keep the treatment-time interaction coefficient constant across partitions while allowing the time effect to vary nonlinearly. This conveniently extends to arbitrary combinations of spline and linear terms without requiring special handling.

When blockfit = TRUE is specified alongside just_linear_without_interactions, the flat-block path provides an alternative enforcement mechanism. Rather than relying on constraint projection, flat coefficients are pooled structurally across partitions during backfitting. The two approaches agree at the point estimate but differ in their uncertainty quantification; see the Blockfit section above.

Integration

Because the fitted object retains an explicit polynomial representation in each partition, numerical integration can be carried out in a fairly direct way. The package wraps that calculation in a user-facing S3 method so the user does not need to manage knot boundaries or partition membership by hand.

In the user-facing interface, numerical integration is applied through integrate.lgspline, which applies Gauss-Legendre quadrature to predictions from the fitted model produced by predict.lgspline.

Implementation

For a user-supplied rectangular domain, integrate.lgspline constructs a tensor-product grid of Gauss-Legendre nodes, evaluates the fitted model at those points, and forms the weighted sum. This works for both univariate and multivariate models, respects the fitted partition structure automatically, and avoids requiring the user to keep track of knot boundaries by hand.

The vars argument selects which predictors are integrated over. Predictors not listed in vars are held fixed at initial_values when supplied, or otherwise at the midpoint of their observed training range. The optional B_predict argument makes it possible to integrate posterior draws or other alternate coefficient sets, and n_quad controls the number of Gauss-Legendre nodes used per integrated dimension.

Integration is performed on the response scale by default. Setting link_scale = TRUE instead integrates the linear predictor \eta = f(\mathbf{t}). For identity-link Gaussian models the two scales coincide.

Lagrange Multipliers

When return_lagrange_multipliers = TRUE, the multiplier vector

\boldsymbol{\lambda} = (\mathbf{A}^{\top}\mathbf{G}\mathbf{A})^{-1}\mathbf{A}^{\top}\hat{\boldsymbol{\beta}}

is returned. These quantify the sensitivity of the penalized objective to relaxing each smoothness or user-supplied equality constraint. When constraint target values are nonzero (\mathbf{A}^{\top}\boldsymbol{\beta}_0 \neq \mathbf{0}), the modified formulation is used:

\boldsymbol{\lambda} = (\mathbf{A}^{\top}\mathbf{G}\mathbf{A})^{-1}\mathbf{A}^{\top}(\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}_0)

where \mathbf{A}^{\top}\boldsymbol{\beta}_0 is the vector of constraint target values. Multipliers are NULL when no constraints are active (\mathbf{A} is NULL or K = 0).

For inequality constraints, multipliers are returned as computed by solve.QP. The Lagrange multipliers for active inequality constraints can be used diagnostically to identify which shape constraints are most costly in terms of goodness of fit.

S3 Methods

Standard S3 methods are provided for objects of class lgspline:

• print.lgspline and summary.lgspline: Provide concise model summaries, with print.summary.lgspline formatting coefficient tables in a familiar regression-style layout.

• logLik.lgspline: Returns a standard logLik object. For Gaussian responses with identity link, the exact log-likelihood is computed. When a correlation structure is present via VhalfInv, the log-likelihood includes the \log|\mathbf{V}^{-1/2}| adjustment and the corresponding whitened quadratic form. For other families, the method falls back to family$aic or a deviance-based approximation. An include_prior argument (default TRUE) optionally adds the Gaussian prior penalty interpretation of the smoothing spline penalty -\frac{1}{2\sigma^{2}}\tilde{\boldsymbol{\beta}}^{\top}\boldsymbol{\Lambda}\tilde{\boldsymbol{\beta}} to obtain a penalized MAP log-likelihood.

• predict.lgspline: Produces fitted values and related quantities (e.g., derivatives and standard errors through se.fit = TRUE), lets new_predictors override newdata, accepts alternate coefficient lists through B_predict, and supports prediction on new predictor matrices consistent with the original spline expansions.

• coef.lgspline: Extracts partition-specific coefficient vectors.

• confint.lgspline: Extracts confidence intervals. When the inverse Hessian from BFGS optimization is available for correlation parameters, intervals for those correlation parameters are returned on the working (transformed) scale and should be back-transformed as described in the correlation section.

• plot.lgspline: For one-dimensional fits, produces base R graphics showing the fitted function (with optional partition-wise formulas) and supports overlay via add = TRUE. For two or more predictors, an interactive plotly-based visualization is returned. Specific predictors may be selected via vars.

• integrate.lgspline: Computes definite integrals of the fitted surface over rectangular domains by Gauss-Legendre quadrature.

Additional user-facing helpers include wald_univariate for coefficient-wise Wald inference, generate_posterior for posterior and posterior-predictive sampling, generate_posterior_correlation for correlation-aware posterior simulation, equation for closed-form display of the fitted partition formulas, and find_extremum for optimizing the fitted surface or a custom acquisition function built from it.

References

Buse, A. and Lim, L. (1977). Cubic Splines as a Special Case of Restricted Least Squares. Journal of the American Statistical Association, 72, 64-68.

Eilers, P. H. and Marx, B. D. (1996). Flexible Smoothing with B-splines and Penalties. Statistical Science, 11(2), 89-121.

Ezhov, N., Neitzel, F. and Petrovic, S. (2018). Spline Approximation, Part 1: Basic Methodology. Journal of Applied Geodesy, 12(2), 139-155.

Goldfarb, D. and Idnani, A. (1983). A Numerically Stable Dual Method for Solving Strictly Convex Quadratic Programs. Mathematical Programming, 27(1), 1-33.

Harville, D. A. (1977). Maximum Likelihood Approaches to Variance Component Estimation and to Related Problems. Journal of the American Statistical Association, 72(358), 320-338.

Hastie, T. J. and Tibshirani, R. J. (1990). Generalized Additive Models. Chapman & Hall/CRC.

Kisi, O., Heddam, S., Parmar, K. S., Petroselli, A., Kulls, C. and Zounemat-Kermani, M. (2025). Integration of Gaussian Process Regression and K Means Clustering for Enhanced Short Term Rainfall Runoff Modeling. Scientific Reports, 15, 7444.

MacQueen, J. B. (1967). Some Methods for Classification and Analysis of Multivariate Observations. In Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1, 281-297. University of California Press.

McCullagh, P. and Nelder, J. A. (1989). Generalized Linear Models. Chapman & Hall, 2nd edition.

Murray, I., Adams, R. P. and MacKay, D. J. C. (2010). Elliptical Slice Sampling. Proceedings of the 13th International Conference on Artificial Intelligence and Statistics (AISTATS), 9, 541-548.

Nocedal, J. and Wright, S. J. (2006). Numerical Optimization (2nd ed.). Springer.

Patterson, H. D. and Thompson, R. (1971). Recovery of Inter-Block Information When Block Sizes Are Unequal. Biometrika, 58, 545-554.

Pya, N. and Wood, S. N. (2015). Shape Constrained Additive Models. Statistics and Computing, 25(3), 543-559.

Reinsch, C. H. (1967). Smoothing by Spline Functions. Numerische Mathematik, 10, 177-183.

Ruppert, D., Wand, M. P. and Carroll, R. J. (2003). Semiparametric Regression. Cambridge University Press.

Searle, S. R., Casella, G. and McCulloch, C. E. (2006). Variance Components. Wiley.

Wahba, G. (1990). Spline Models for Observational Data. SIAM.

Wood, S. N. (2006). On Confidence Intervals for Generalized Additive Models Based on Penalized Regression Splines. Australian & New Zealand Journal of Statistics, 48(4), 445-464.

Wood, S. N. (2011). Fast Stable Restricted Maximum Likelihood and Marginal Likelihood Estimation of Semiparametric Generalized Linear Models. Journal of the Royal Statistical Society: Series B, 73(1), 3-36.

Wood, S. N. (2017). Generalized Additive Models: An Introduction with R. CRC Press, 2nd edition.

Efficient Matrix Multiplication of G and A Matrices

Description

Efficient Matrix Multiplication of G and A Matrices

Usage

GAmult_wrapper(
  G,
  A,
  K,
  p_expansions,
  R_constraints,
  parallel,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks
)

Arguments

Details

Computes G Processes in parallel chunks if enabled.

Value

List of matrix products, one per partition

Finite-difference Gradient Computer

Description

Computes a finite-difference approximation derived from the objective in fn at x, returned in the sign convention currently used by efficient_bfgs.

Usage

approx_grad(x, fn, eps = sqrt(.Machine$double.eps))

Arguments

Details

Used within efficient_bfgs when fn does not supply a usable gradient. In the main lgspline() correlation-optimization path, stats::optim() is used instead when no analytic REML gradient is available. Internally this helper returns the negated central-difference approximation rather than the raw derivative.

Value

Numeric vector of finite-difference approximated gradient values in the current optimizer sign convention

Matrix Inversion using Armadillo

Description

Computes the inverse of a matrix using Armadillo's inversion method

Usage

armaInv(x)

Arguments

Value

Inverted matrix

Backfitting Solver for Blockfit Models

Description

Fits models with mixed spline and non-interactive linear ("flat") terms using an iterative backfitting approach. Flat terms receive a single shared coefficient across partitions (rather than K+1 partition-specific coefficients constrained to equality), reducing the effective parameter count and improving efficiency when the number of flat terms is large relative to spline terms.

The backfitting loop alternates between:

1. Spline step

   Fit spline terms on the response adjusted for the current flat contribution using a constrained penalized least squares solve.

   If \mathbf{y}_k is the response vector for partition k, \mathbf{Z}_k the spline design matrix, \mathbf{X}_{\mathrm{flat}}^{(k)} the flat design matrix, and \mathbf{v} the flat coefficient vector, then

   \boldsymbol{\beta}_{\mathrm{spline}}^{(k)} = \arg\min_{\boldsymbol{\beta}} \left\| \mathbf{y}_k - \mathbf{Z}_k \boldsymbol{\beta} - \mathbf{X}_{\mathrm{flat}}^{(k)} \mathbf{v} \right\|^2 + \boldsymbol{\beta}^{\top} \mathbf{\Lambda}_{\mathrm{spline}} \boldsymbol{\beta}

   subject to

   \mathbf{A}_{\mathrm{spline}}^{\top} \boldsymbol{\beta} = \mathbf{c}_{\mathrm{spline}}.

2. Flat step

   Update flat coefficients via pooled penalized regression on residuals, subject to the residual equality constraint from any mixed constraints:

   \mathbf{v} = \arg\min_{\mathbf{v}} \sum_{k=0}^{K} \left\| \mathbf{y}_k - \mathbf{Z}_k \boldsymbol{\beta}_{\mathrm{spline}}^{(k)} - \mathbf{X}_{\mathrm{flat}}^{(k)} \mathbf{v} \right\|^2 + \mathbf{v}^{\top} \mathbf{\Lambda}_{\mathrm{flat}} \mathbf{v}

   subject to

   \tilde{\mathbf{A}}_{\mathrm{flat}}^{\top} \mathbf{v} = \mathbf{c} - \mathbf{A}_{\mathrm{spline}}^{\top} \boldsymbol{\beta}_{\mathrm{spline}}.

Four estimation paths are selected automatically based on the model configuration:

Case (a)

Gaussian identity + GEE: closed-form solve via .bf_case_gauss_gee.

Case (b)

Gaussian identity, no correlation: standard backfitting via .bf_case_gauss_no_corr.

Case (c)

GLM + GEE: two-stage (damped Newton-Raphson warm start then damped SQP) via .bf_case_glm_gee.

Case (d)

GLM without GEE: damped Newton-Raphson + backfitting via .bf_case_glm_no_corr.

When quadprog = TRUE without GEE, inequality constraints are enforced after backfitting convergence. The constraint handling method is selected automatically by inspecting the sparsity pattern of qp_Amat: if every column has nonzeros in only one partition block (e.g.\ derivative sign or range constraints), a partition-wise active-set method is used via .active_set_refine, avoiding the dense P \times P system. If any column spans multiple partition blocks (e.g.\ cross-knot monotonicity), or if the active-set method does not converge, the dense SQP fallback via .bf_sqp_loop is used. GEE paths always use the dense system.

After convergence, coefficients are reassembled into the standard per-partition format (flat coefficients replicated across partitions) for compatibility with downstream inference.

Usage

blockfit_solve(
  X,
  y,
  flat_cols,
  K,
  p_expansions,
  Lambda,
  L_partition_list,
  unique_penalty_per_partition,
  A,
  R_constraints,
  constraint_values,
  X_gram,
  Ghalf_full,
  GhalfInv_full,
  family,
  order_list,
  glm_weight_function,
  schur_correction_function,
  need_dispersion_for_estimation,
  dispersion_function,
  observation_weights,
  homogenous_weights = TRUE,
  iterate,
  tol,
  parallel_eigen,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks,
  return_G_getB,
  quadprog = FALSE,
  qp_Amat = NULL,
  qp_bvec = NULL,
  qp_meq = NULL,
  qp_score_function = NULL,
  keep_weighted_Lambda = FALSE,
  max_backfit_iter = 100,
  Vhalf = NULL,
  VhalfInv = NULL,
  include_warnings = TRUE,
  verbose = FALSE,
  ...
)

Arguments

Value

A list with elements:

B

List of K+1 coefficient vectors (p_expansions \times 1), flat coefficients replicated across partitions.

G_list

List containing \mathbf{G}, \mathbf{G}^{1/2}, and \mathbf{G}^{-1/2}, each a list of K+1 p_expansions \times p_expansions matrices, or NULL if return_G_getB is FALSE.

\mathbf{G} satisfies \mathbf{G} = \mathbf{G}^{1/2} (\mathbf{G}^{1/2})^{\top} exactly.

When a correlation structure is present, these matrices are obtained from the dense whitened information matrix and then split back into per-partition blocks. When flat columns are present, the returned matrices also reflect the pooled flat-column information across partitions.

qp_info

NULL when no inequality-refinement metadata are produced. Otherwise a list of available QP or active-set diagnostics. Dense SQP solves include solution, lagrangian, active_constraints, iact, Amat_active, bvec_active, meq_active, converged, and final_deviance; non-GEE dense solves also carry info_matrix, Amat_combined, bvec_combined, and meq_combined. Partition-wise active-set refinement returns the corresponding active-constraint summary only.

See Also

lgspline, lgspline.fit, get_B

Examples

## Not run: 
## Minimal verification example: Gaussian identity, no correlation,
## 2 partitions, 3 spline columns + 1 flat column per partition.
##
## This confirms that blockfit_solve returns compatible output and
## that flat coefficients are replicated across partitions.

set.seed(1234)
n1 <- 50; n2 <- 50
p_expansions <- 4  # intercept + x + x^2 + z (flat)
K  <- 1  # 2 partitions

X1 <- cbind(1, rnorm(n1), rnorm(n1)^2, rnorm(n1))
X2 <- cbind(1, rnorm(n2), rnorm(n2)^2, rnorm(n2))
beta_true <- c(1, 0.5, -0.3, 0.8)
y1 <- X1 %*% beta_true + rnorm(n1, 0, 0.5)
y2 <- X2 %*% beta_true + rnorm(n2, 0, 0.5)

## Constraint: spline coefficients equal across partitions
## (columns 1:3 constrained, column 4 is flat)
A <- matrix(0, nrow = p_expansions * (K + 1), ncol = 3)
for(j in 1:3){
  A[j, j]      <-  1
  A[p_expansions + j, j]  <- -1
}
qr_A <- qr(A)
A <- qr.Q(qr_A)[, 1:qr_A$rank, drop = FALSE]

Lambda <- diag(p_expansions) * 1e-4
L_partition_list <- list(0, 0)
X_gram <- list(crossprod(X1), crossprod(X2))

G_list <- compute_G_eigen(
  X_gram, Lambda, K,
  parallel = FALSE, cl = NULL,
  chunk_size = NULL, num_chunks = NULL, rem_chunks = NULL,
  family = gaussian(),
  unique_penalty_per_partition = FALSE,
  L_partition_list = L_partition_list,
  keep_G = TRUE,
  schur_corrections = list(0, 0))

result <- blockfit_solve(
  X = list(X1, X2),
  y = list(y1, y2),
  flat_cols = 4L,
  K = K, p_expansions = p_expansions,
  Lambda = Lambda,
  L_partition_list = L_partition_list,
  unique_penalty_per_partition = FALSE,
  A = A, R_constraints = ncol(A),
  constraint_values = list(),
  X_gram = X_gram,
  Ghalf_full = G_list$Ghalf,
  GhalfInv_full = G_list$GhalfInv,
  family = gaussian(),
  order_list = list(1:n1, (n1+1):(n1+n2)),
  glm_weight_function = function(mu, y, oi, fam, d, ow, ...) rep(1, length(y)),
  schur_correction_function = function(X, y, B, d, ol, K, fam, ow, ...) {
    lapply(1:(K + 1), function(k) 0)
  },
  need_dispersion_for_estimation = FALSE,
  dispersion_function = function(...) 1,
  observation_weights = list(rep(1, n1), rep(1, n2)),
  homogenous_weights = TRUE,
  iterate = TRUE, tol = 1e-8,
  parallel = FALSE, cl = NULL,
  chunk_size = NULL, num_chunks = NULL, rem_chunks = NULL,
  return_G_getB = TRUE,
  verbose = FALSE)

## Verify: flat coefficient (col 4) is identical across partitions
stopifnot(abs(result$B[[1]][4] - result$B[[2]][4]) < 1e-10)

## Verify: output structure
stopifnot(length(result$B) == K + 1)
stopifnot(length(result$B[[1]]) == p_expansions)
stopifnot(!is.null(result$G_list))

## End(Not run)

Extract Coefficients from a Fitted lgspline

Description

Returns the per-partition polynomial coefficient lists from a fitted lgspline model.

Usage

## S3 method for class 'lgspline'
coef(object, ...)

Arguments

Details

Coefficient names reflect the polynomial expansion terms, e.g.:

• intercept

• v: linear term for predictor v

• v_^2: quadratic term

• v^3: cubic term

• _v_x_w_: two-way interaction

Column/variable names replace numeric indices when available.

To get all coefficients as a single matrix: Reduce('cbind', coef(model_fit)).

Value

A list of per-partition coefficient vectors. Returns NULL with a warning if object$B is not found.

See Also

lgspline

Examples

set.seed(1234)
t <- runif(1000, -10, 10)
y <- 2*sin(t) + -0.06*t^2 + rnorm(length(t))
model_fit <- lgspline(t, y)

coefficients <- coef(model_fit)
print(coefficients[[1]])
print(Reduce('cbind', coefficients))

Extract Coefficients from a wald_lgspline Object

Description

Extract Coefficients from a wald_lgspline Object

Usage

## S3 method for class 'wald_lgspline'
coef(object, ...)

Arguments

Value

Named numeric vector of coefficient estimates, or NULL if no estimate column is available.

Collapse Matrix List into a Single Dense Block-Layout Matrix

Description

Transforms a list of matrices into a single dense block-layout matrix. This is useful especially for quadratic programming problems, where operating on lists of blocks is not convenient.

Usage

collapse_block_diagonal(matlist)

Arguments

Value

Dense matrix with each input block placed in its own column range

Compute Eigenvalues and Related Matrices for G

Description

Computes partition-wise inverse-information matrices and their matrix square roots from the penalized information matrix via eigendecomposition.

Computes partition-wise inverse-information matrices and their matrix square roots from the penalized information matrix via eigendecomposition.

Usage

compute_G_eigen(
  X_gram,
  Lambda,
  K,
  parallel,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks,
  family,
  unique_penalty_per_partition,
  L_partition_list,
  keep_G = TRUE,
  schur_corrections
)

compute_G_eigen(
  X_gram,
  Lambda,
  K,
  parallel,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks,
  family,
  unique_penalty_per_partition,
  L_partition_list,
  keep_G = TRUE,
  schur_corrections
)

Arguments

Value

A list with components G (or NULL blocks when keep_G = FALSE), Ghalf, and optionally GhalfInv.

A list with components G (or NULL blocks when keep_G = FALSE), Ghalf, and optionally GhalfInv.

Compute Component \textbf{G}^{1/2}\textbf{A}(\textbf{A}^{T}\textbf{G}\textbf{A})^{-1}\textbf{A}^{T}\textbf{G}\textbf{X}^{T}\textbf{y}

Description

Compute Component \textbf{G}^{1/2}\textbf{A}(\textbf{A}^{T}\textbf{G}\textbf{A})^{-1}\textbf{A}^{T}\textbf{G}\textbf{X}^{T}\textbf{y}

Usage

compute_GhalfXy_temp_wrapper(
  G,
  Ghalf,
  A,
  AGAInv,
  Xy,
  nc,
  K,
  parallel,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks
)

Arguments

Details

Computes the least-squares projection component \mathbf{G}^{1/2}\mathbf{A}(\mathbf{A}^{T}\mathbf{G}\mathbf{A})^{-1} \mathbf{A}^{T}\mathbf{G}\mathbf{X}^{T}\mathbf{y} together with the intermediate product (\mathbf{A}^{T}\mathbf{G}\mathbf{A})^{-1}\mathbf{A}^{T}\mathbf{G}\mathbf{X}^{T}\mathbf{y} for reuse downstream.

Value

Unnamed two-element list containing the projected result vector and the intermediate (\mathbf{A}^{T}\mathbf{G}\mathbf{A})^{-1}\mathbf{A}^{T}\mathbf{G}\mathbf{X}^{T}\mathbf{y} product.

Construct Smoothing Spline Penalty Matrix

Description

Builds penalty matrix combining smoothing spline and ridge penalties with optional predictor/partition-specific components.

Usage

compute_Lambda(
  custom_penalty_mat,
  L1,
  wiggle_penalty,
  flat_ridge_penalty,
  K,
  p_expansions,
  unique_penalty_per_predictor,
  unique_penalty_per_partition,
  penalty_vec,
  colnm_expansions,
  just_Lambda = TRUE
)

Arguments

Value

List containing Lambda, L1, L2, L_predictor_list, L_partition_list; or just Lambda if just_Lambda=TRUE and no partition penalties.

Compute Derivative of Inverse-Information Matrix G with Respect to Lambda

Description

Calculates the derivative of the inverse-information matrix \textbf{G} with respect to the smoothing parameter \lambda, supporting both shared and partition-specific penalties.

Usage

compute_dG_dlambda(
  G,
  L,
  K,
  lambda,
  unique_penalty_per_partition,
  L_partition_list,
  parallel,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks
)

Arguments

Value

A list of derivative matrices d\textbf{G}/d\lambda for each partition

Derivative of Constrained Penalized Coefficients with Respect to Lambda

Description

Computes d(\mathbf{U}\mathbf{G}\mathbf{X}^{T}\mathbf{y})/d\lambda, the sensitivity of the constrained coefficient estimates \tilde{\boldsymbol{\beta}} = \mathbf{U}\mathbf{G}\mathbf{X}^T\mathbf{y} to changes in the smoothing parameter \lambda.

Usage

compute_dG_u_dlambda_xy(
  AGAInv_AGXy,
  AGAInv,
  G,
  A,
  dG_dlambda,
  nc,
  nca,
  K,
  Xy,
  Ghalf,
  dGhalf,
  GhalfXy_temp,
  parallel,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks
)

Arguments

Details

This function is called during GCV/penalty optimization and computes how the constrained coefficient vector changes as the penalty weight varies. Two implementations are provided depending on problem size.

Derivation

The constrained estimate is \tilde{\boldsymbol{\beta}} = \mathbf{U}\hat{\boldsymbol{\beta}} where \hat{\boldsymbol{\beta}} = \mathbf{G}\mathbf{X}^T\mathbf{y} and \mathbf{U} = \mathbf{I} - \mathbf{G}\mathbf{A}(\mathbf{A}^T\mathbf{G}\mathbf{A})^{-1}\mathbf{A}^T. Both \mathbf{G} and \mathbf{U} depend on \lambda through \mathbf{G} = (\mathbf{X}^T\mathbf{X} + \lambda\boldsymbol{\Lambda})^{-1}.

Differentiating the product \mathbf{U}\mathbf{G}\mathbf{X}^T\mathbf{y} requires the chain rule applied to three \lambda-dependent components:

\frac{d}{d\lambda}(\mathbf{U}\mathbf{G}\mathbf{X}^T\mathbf{y})

= \frac{d\mathbf{U}}{d\lambda}\hat{\boldsymbol{\beta}} + \mathbf{U}\frac{d\hat{\boldsymbol{\beta}}}{d\lambda}

The unconstrained part is straightforward: d\hat{\boldsymbol{\beta}}/d\lambda = (d\mathbf{G}/d\lambda)\mathbf{X}^T\mathbf{y}, computed partition-wise since \mathbf{G} is block-diagonal.

The constraint projection derivative is more involved. Writing \mathbf{C} = (\mathbf{A}^T\mathbf{G}\mathbf{A})^{-1} and expanding d\mathbf{U}/d\lambda gives three terms (after applying the product rule and the identity d\mathbf{C}/d\lambda = -\mathbf{C}(d(\mathbf{A}^T\mathbf{G}\mathbf{A})/d\lambda)\mathbf{C}):

term2a

(d\mathbf{G}/d\lambda)\mathbf{A}\mathbf{C}\mathbf{A}^T\hat{\boldsymbol{\beta}} direct effect of d\mathbf{G}/d\lambda on the projection

term2b

\mathbf{G}\mathbf{A}\mathbf{C}(d(\mathbf{A}^T\mathbf{G}\mathbf{A})/d\lambda)\mathbf{C}\mathbf{A}^T\hat{\boldsymbol{\beta}} effect through the change in (\mathbf{A}^T\mathbf{G}\mathbf{A})^{-1}

term2c

\mathbf{G}\mathbf{A}\mathbf{C}\mathbf{A}^T(d\hat{\boldsymbol{\beta}}/d\lambda) projection of the unconstrained derivative back through the constraint

The full derivative is d\hat{\boldsymbol{\beta}}/d\lambda minus the sum of these three correction terms.

Large problem path (K >= 10, nc > 4)

Computes the three correction terms explicitly using the block structure of \mathbf{G} and d\mathbf{G}/d\lambda. The intermediate quantity d(\mathbf{A}^T\mathbf{G}\mathbf{A})/d\lambda = \mathbf{A}^T(d\mathbf{G}/d\lambda)\mathbf{A} is accumulated partition-wise. Shared vectors \mathbf{A}\mathbf{C}\mathbf{A}^T\hat{\boldsymbol{\beta}} and related products are precomputed once and reused across partitions. Parallelism is over chunks of partitions.

Small problem path

Reformulates the derivative using the matrix square root \mathbf{G}^{1/2} and its derivative d\mathbf{G}^{1/2}/d\lambda. The constraint \mathbf{A}^T\boldsymbol{\beta} = 0 is imposed via least squares projection: the residuals from regressing \mathbf{G}^{1/2}\mathbf{X}^T\mathbf{y} onto \mathbf{G}^{1/2}\mathbf{A} give the constrained component, and differentiating this projection with respect to \lambda yields the derivative. Uses .lm.fit for speed with a stabilizing rescaling factor to prevent numerical issues when the constraint matrix is poorly scaled.

Value

P \times 1 vector of derivatives d\tilde{\boldsymbol{\beta}}/d\lambda

Compute Eigen-Based Square-Root Factors for d\textbf{G}/d\lambda

Description

Applies an eigendecomposition-based matrix square root to each partition-wise d\textbf{G}/d\lambda matrix after replacing NA entries with 0 and treating all-0 inputs with an identity fallback.

Usage

compute_dGhalf(
  dG_dlambda,
  p_expansions,
  K,
  parallel,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks
)

Arguments

Value

List of p \times p eigen-based square-root factors derived from the partition-wise d\textbf{G}/d\lambda matrices

Compute Derivative of Inverse-Information Matrix G with Respect to Lambda (Wrapper)

Description

Wrapper for the derivative of the trace term used in GCV / effective- degrees-of-freedom calculations with respect to \lambda.

Usage

compute_dW_dlambda_wrapper(
  G,
  A,
  GXX,
  Ghalf,
  dG_dlambda,
  dGhalf_dlambda,
  AGAInv,
  nc,
  K,
  parallel,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks
)

Arguments

Value

Scalar value representing the trace derivative component.

Compute Gram Matrix for Block Diagonal Structure

Description

Compute Gram Matrix for Block Diagonal Structure

Usage

compute_gram_block_diagonal(
  list_in,
  parallel,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks
)

Arguments

Details

For a list of matrices, will compute the gram matrix of each element of the list.

Value

List of Gram matrices (\textbf{X}^{T}\textbf{X}) for each block

Effective degrees of freedom via trace of the hat matrix

Description

Computes \mathrm{tr}(\mathbf{X}\mathbf{U}\mathbf{G}\mathbf{X}^\top \mathbf{W}^{1/2}\mathbf{V}^{-1}\mathbf{W}^{1/2}) without forming the N \times N hat matrix, by reducing the problem to P \times P and r \times r operations.

Usage

compute_trace_H(
  G,
  Lambda,
  A,
  AGAInv,
  nc,
  nca,
  K,
  parallel,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks,
  unique_penalty_per_partition,
  L_partition_list
)

Arguments

Details

Derivation

Let \mathbf{G} = (\mathbf{X}^\top\mathbf{W}^{1/2}\mathbf{V}^{-1}\mathbf{W}^{1/2}\mathbf{X} + \boldsymbol{\Lambda})^{-1} and \mathbf{U} = \mathbf{I} - \mathbf{G}\mathbf{A}(\mathbf{A}^\top\mathbf{G}\mathbf{A})^{-1}\mathbf{A}^\top. By the cyclic property of the trace:

\mathrm{tr}(\mathbf{X}\mathbf{U}\mathbf{G}\mathbf{X}^\top\mathbf{W}^{1/2}\mathbf{V}^{-1}\mathbf{W}^{1/2}) = \mathrm{tr}(\mathbf{U}\mathbf{G}\,\mathbf{X}^\top\mathbf{W}^{1/2}\mathbf{V}^{-1}\mathbf{W}^{1/2}\mathbf{X})

Substituting \mathbf{X}^\top\mathbf{W}^{1/2}\mathbf{V}^{-1}\mathbf{W}^{1/2}\mathbf{X} = \mathbf{G}^{-1} - \boldsymbol{\Lambda}:

= \mathrm{tr}(\mathbf{U}\mathbf{G}(\mathbf{G}^{-1} - \boldsymbol{\Lambda})) = \mathrm{tr}(\mathbf{U}) - \mathrm{tr}(\mathbf{U}\mathbf{G}\boldsymbol{\Lambda})

Since \mathbf{U} is idempotent, \mathrm{tr}(\mathbf{U}) = \mathrm{rank}(\mathbf{U}) = P - r where r = \mathrm{rank}(\mathbf{A}). For the second term, expand \mathbf{U} = \mathbf{I} - \mathbf{G}\mathbf{A}(\mathbf{A}^\top\mathbf{G}\mathbf{A})^{-1}\mathbf{A}^\top:

\mathrm{tr}(\mathbf{U}\mathbf{G}\boldsymbol{\Lambda}) = \mathrm{tr}(\mathbf{G}\boldsymbol{\Lambda}) - \mathrm{tr}\!\left((\mathbf{A}^\top\mathbf{G}\mathbf{A})^{-1}\mathbf{A}^\top\mathbf{G}\boldsymbol{\Lambda}\mathbf{G}\mathbf{A}\right)

using the cyclic property on the second term. Both \mathbf{G} and \boldsymbol{\Lambda} are block-diagonal with K+1 blocks of dimension p \times p, so:

\mathrm{tr}(\mathbf{G}\boldsymbol{\Lambda}) = \sum_{k=1}^{K+1}\mathrm{tr}(\mathbf{G}_k\boldsymbol{\Lambda}_k)

and \mathbf{G}\boldsymbol{\Lambda}\mathbf{G} is also block-diagonal with blocks \mathbf{G}_k\boldsymbol{\Lambda}_k\mathbf{G}_k.

Combining:

\mathrm{tr}(\mathbf{X}\mathbf{U}\mathbf{G}\mathbf{X}^\top\mathbf{W}^{1/2}\mathbf{V}^{-1}\mathbf{W}^{1/2}) = (P - r) - \sum_{k=1}^{K+1}\mathrm{tr}(\mathbf{G}_k\boldsymbol{\Lambda}_k) + \mathrm{tr}\!\left((\mathbf{A}^\top\mathbf{G}\mathbf{A})^{-1}\mathbf{A}^\top\mathbf{G}\boldsymbol{\Lambda}\mathbf{G}\mathbf{A}\right)

The correction term \mathrm{tr}((\mathbf{A}^\top\mathbf{G}\mathbf{A})^{-1}\mathbf{A}^\top\mathbf{G}\boldsymbol{\Lambda}\mathbf{G}\mathbf{A}) is an r \times r trace and captures the degrees of freedom recovered because the constraints pin certain linear combinations of penalized coefficients, removing them from estimation. When \boldsymbol{\Lambda} = \mathbf{0} (no penalty), both \mathrm{tr}(\mathbf{G}\boldsymbol{\Lambda}) and the correction vanish, giving \mathrm{edf} = P - r. When \mathbf{U} = \mathbf{I} (no constraints), r = 0 and the correction vanishes, giving \mathrm{edf} = P - \mathrm{tr}(\mathbf{G}\boldsymbol{\Lambda}).

The result is invariant to the correlation structure \mathbf{V} and GLM weights \mathbf{W}, which enter only through \mathbf{G}.

Partition-specific penalties

When unique_penalty_per_partition = TRUE, each partition k has an additional penalty matrix \mathbf{L}_k from L_partition_list, so the effective per-partition penalty becomes \boldsymbol{\Lambda}_k = \boldsymbol{\Lambda} + \mathbf{L}_k. All block-wise traces and products use \boldsymbol{\Lambda}_k in place of the shared \boldsymbol{\Lambda}.

Computational cost

The partition-wise traces \mathrm{tr}(\mathbf{G}_k\boldsymbol{\Lambda}_k) cost O(p^2) each and are parallelizable. The correction term requires forming \mathbf{G}\boldsymbol{\Lambda}\mathbf{G}\mathbf{A} (block-diagonal times sparse, O((K+1)p^2 r)) and a single r \times r trace (O(r^2)). The total cost is O((K+1)p^2 + r^2), compared to O(N^2) or O(NP) for forming and tracing the full hat matrix.

Value

Scalar effective degrees of freedom, clamped to [0, \, (K+1)p].

Calculate Trace of Matrix Product \text{trace}(\textbf{X}\textbf{U}\textbf{G}\textbf{X}^{T})

Description

Calculate Trace of Matrix Product \text{trace}(\textbf{X}\textbf{U}\textbf{G}\textbf{X}^{T})

Usage

compute_trace_UGXX_wrapper(
  G,
  A,
  GXX,
  AGAInv,
  nc,
  nca,
  K,
  parallel,
  cl,
  chunk_size,
  num_chunks,
  rem_chunks
)

Arguments

Details

Computes \text{trace}(\textbf{X}\textbf{U}\textbf{G}\textbf{X}^{T}) where \textbf{U} = \textbf{I} - \textbf{G}\textbf{A}(\textbf{A}^{T}\textbf{G}\textbf{A})^{-1}\textbf{A}^{T}. Handles parallel computation by splitting into chunks.

Value

Trace value

Confidence Intervals for lgspline Coefficients

Description

Wald-based confidence intervals for regression coefficients and, when available, correlation parameters (on the working scale).

Usage

## S3 method for class 'lgspline'
confint(object, parm, level = 0.95, ...)

Arguments

Details

For Gaussian identity-link models, t-distribution quantiles are used with effective degrees of freedom N - \mathrm{trace}(\mathbf{XUGX}^\top). All other families use normal quantiles.

Correlation parameter intervals (if VhalfInv_params_estimates and VhalfInv_params_vcov are present) are computed on the unbounded working scale via a Wald interval.

Value

A matrix with columns giving lower and upper confidence limits, named e.g. 2.5 % and 97.5 % for 95% intervals. When available, rows for working-scale correlation parameters are appended after the regression coefficients.

Extract Confidence Intervals from a wald_lgspline Object

Description

Extract Confidence Intervals from a wald_lgspline Object

Usage

## S3 method for class 'wald_lgspline'
confint(object, parm = NULL, level = NULL, ...)

Arguments

Value

Matrix with columns lower and upper, or NULL if confidence limits are not available.

Cox PH Dispersion Function

Description

Returns 1 unconditionally. Cox PH has no dispersion parameter; this function exists solely for interface compatibility with lgspline's dispersion_function argument.

Usage

cox_dispersion_function(
  mu,
  y,
  order_indices,
  family,
  observation_weights,
  VhalfInv,
  ...
)

Arguments