Smoothing

Splines

Spline basis constructors.

class aurora.smoothing.splines.BSplineBasis(knots, degree=3, cyclic=False)

Bases: object

B-spline basis functions via Cox-de Boor recursion.

B-splines are piecewise polynomials with compact support, defined recursively using the Cox-de Boor formula. They form a basis for the space of splines of a given degree with specified knots.

Parameters:

• knots (array-like) – Knot vector (should include boundary knots with multiplicity = degree + 1). For open B-splines, repeat boundary knots degree+1 times.

• degree (int, default=3) – Degree of the B-spline basis (3 for cubic B-splines).

• cyclic (bool, default=False) – Whether to use cyclic (periodic) boundary conditions.

Variables:

• knots (ndarray) – Full knot vector including repeated boundary knots.

• degree (int) – Degree of the splines.

• n_basis (int) – Number of basis functions.

• cyclic (bool) – Whether cyclic boundary conditions are used.

Notes

For degree p and n_basis basis functions, the knot vector must have length n_basis + p + 1.

References

de Boor, C. (1978). A Practical Guide to Splines. Springer-Verlag. Eilers, P.H.C. & Marx, B.D. (1996). Flexible smoothing with B-splines and

penalties. Statistical Science, 11(2), 89-121.

basis_matrix(x, sparse=False)

Compute B-spline basis matrix using Cox-de Boor recursion.

Parameters:

• x (array-like, shape (n_samples,)) – Points at which to evaluate basis functions.

• sparse (bool, default=False) – If True, return scipy.sparse.csr_matrix (only for NumPy backend). If False, return dense array compatible with input backend.

Returns:

B – Basis matrix where B[i,j] = B_j^p(x_i), the j-th B-spline of degree p evaluated at x_i.

Return type:

array or csr_matrix, shape (n_samples, n_basis)

Notes

Uses the Cox-de Boor recursion formula:

B_i^0(x) = 1 if t_i ≤ x < t_{i+1}, else 0 B_i^p(x) = w_i^p(x) B_i^{p-1}(x) + (1 - w_{i+1}^p(x)) B_{i+1}^{p-1}(x)

where w_i^p(x) = (x - t_i) / (t_{i+p} - t_i)

Sparse output exploits the compact support property: only (degree + 1) basis functions are non-zero at any point. For degree=3, each row has at most 4 non-zero entries, regardless of n_basis. This reduces: - Memory: O(n_samples × n_basis) → O(n_samples × degree) - Computation: O(n_samples × n_basis × degree²) → O(n_samples × degree²)

Sparse output is only available for NumPy backend. PyTorch and JAX backends always return dense arrays (sparse support is experimental).

penalty_matrix(order=2)

Compute difference penalty matrix.

For B-splines, we typically use a difference penalty that approximates the integrated squared derivative penalty.

Parameters:

order (int, default=2) – Order of differences (2 for approximating second derivative).

Returns:

S – Penalty matrix where β’Sβ approximates ∫ [f^(m)(x)]^2 dx.

Return type:

ndarray, shape (n_basis, n_basis)

Notes

The difference penalty is:

S = D’D

where D is the order-th difference matrix.

For order=2: D_ij penalizes (β_i - 2β_{i+1} + β_{i+2})² This approximates the integrated squared second derivative.

static create_knots(x, n_basis=10, degree=3, method='quantile')

Create knot vector for B-splines from data.

Parameters:

• x (array-like) – Data values.

• n_basis (int) – Number of basis functions desired.

• degree (int) – Degree of B-splines.

• method (str) – Knot placement method: - ‘quantile’: Place interior knots at quantiles - ‘uniform’: Place interior knots uniformly

Returns:

knots – Full knot vector with repeated boundary knots.

Return type:

ndarray

Notes

For open B-splines, boundary knots are repeated (degree + 1) times. This ensures the spline interpolates at boundaries.

derivative_basis_matrix(x, order=1)

Compute derivative of B-spline basis functions.

Parameters:

• x (array-like, shape (n_samples,)) – Evaluation points.

• order (int, default=1) – Order of derivative (1 for first derivative, etc.).

Returns:

dB – Matrix of derivative values.

Return type:

array, shape (n_samples, n_basis)

Notes

Uses the derivative formula:

dB_i^p(x)/dx = p * [B_i^{p-1}(x)/(t_{i+p}-t_i) -

B_{i+1}^{p-1}(x)/(t_{i+p+1}-t_{i+1})]

class aurora.smoothing.splines.CubicSplineBasis(knots, boundary_knots=None)

Bases: object

Natural cubic spline basis for smooth function representation.

Natural cubic splines provide a flexible basis for representing smooth functions. They are especially useful in GAMs for modeling non-linear relationships.

Parameters:

• knots (array-like) – Interior knot locations. Should be sorted in ascending order.

• boundary_knots (tuple of float, optional) – (lower, upper) boundary knots. If None, uses min and max of knots.

Variables:

• knots (ndarray) – Interior knots

• boundary_knots (tuple) – Boundary knots

• n_basis (int) – Number of basis functions (len(knots) + 2 for natural splines)

References

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

basis_matrix(x)

Compute the basis matrix B where B[i,j] = b_j(x_i).

Parameters:

x (array-like, shape (n_samples,)) – Points at which to evaluate the basis functions.

Returns:

B – Basis matrix. Uses the same array backend as input x.

Return type:

array, shape (n_samples, n_basis)

Notes

The basis functions are constructed to satisfy: - Cubic polynomials between knots - Continuous second derivatives at knots - Zero second derivatives at boundaries (natural condition)

penalty_matrix()

Compute the penalty matrix S for integrated squared second derivative.

The penalty matrix satisfies: β’Sβ = ∫ [f’’(x)]² dx

Returns:

S – Penalty matrix (always NumPy array).

Return type:

ndarray, shape (n_basis, n_basis)

Notes

For natural cubic splines, the penalty can be computed analytically. The first two basis functions (constant and linear) receive zero penalty since their second derivatives are zero.

For the ESL (5.4-5.5) basis N_{j+2} = d_j - d_{K-1} with d_j(x) = [(x - xi_j)_+^3 - (x - xi_K)_+^3] / (xi_K - xi_j), the second derivative is d_j’’(x) = 6[(x - xi_j)_+ - (x - xi_K)_+]/(xi_K - xi_j). Since (x - xi_K)_+ = 0 on [xi_1, xi_K] and every N_j’’ vanishes outside the boundary knots, the integral over the real line reduces to

∫ d_i’’ d_j’’ dx = J(xi_i, xi_j) / [(xi_K - xi_i)(xi_K - xi_j)]

with J(t, s) = ∫ 6(x - t)_+ 6(x - s)_+ dx, and

S[i+2, j+2] = D(i,j) - D(i,K-1) - D(K-1,j) + D(K-1,K-1).

static create_knots(x, n_knots=10, method='quantile')

Create knot locations from data.

Parameters:

• x (array-like) – Data values used to determine knot positions.

• n_knots (int, default=10) – Number of interior knots to place.

• method (str, default='quantile') –

  Knot placement method:

  • 'quantile': Place knots at quantiles of x, adapting to the empirical data distribution.

  • 'uniform': Place knots uniformly between min(x) and max(x).

Returns:

knots – Sorted interior knot locations. May contain fewer than n_knots entries if duplicate quantiles are merged.

Return type:

ndarray, shape (n_knots,) or fewer

Raises:

ValueError – If n_knots is less than 1 or method is unknown.

Examples

>>> x = np.random.randn(500)
>>> knots = CubicSplineBasis.create_knots(x, n_knots=8)
>>> len(knots) <= 8
True

class aurora.smoothing.splines.PSplineBasis(n_basis=20, degree=3, penalty_order=2, domain=None)

Bases: object

P-Spline basis: B-splines with equally-spaced knots and difference penalty.

P-Splines use a large number of B-spline basis functions with a difference penalty on coefficients to achieve smooth fits. The amount of smoothing is controlled by a single parameter λ.

Parameters:

• n_basis (int, default=20) – Number of B-spline basis functions. More basis functions allow more flexibility but require larger penalty. Typically 15-40 is sufficient.

• degree (int, default=3) – B-spline degree. Cubic (degree=3) is standard and provides C² continuity.

• penalty_order (int, default=2) – Order of difference penalty: - 1: Penalizes first differences (roughness of coefficients) - 2: Penalizes second differences (curvature, most common) - 3: Penalizes third differences (rate of curvature change)

• domain (tuple of float, optional) – (min, max) domain for the basis. If None, determined from data.

Variables:

• n_basis (int) – Number of basis functions

• degree (int) – B-spline degree

• penalty_order (int) – Order of difference penalty

• knots (ndarray) – Knot locations (set after fitting)

Examples

>>> import numpy as np
>>> from aurora.smoothing.splines import PSplineBasis
>>>
>>> # Create basis
>>> basis = PSplineBasis(n_basis=25, degree=3, penalty_order=2)
>>>
>>> # Generate data
>>> x = np.linspace(0, 2*np.pi, 100)
>>> y = np.sin(x) + np.random.normal(0, 0.2, 100)
>>>
>>> # Fit with automatic smoothing selection
>>> result = basis.fit(x, y, lambda_='gcv')
>>> y_smooth = result.fitted_values_

Notes

The P-spline approach separates the basis complexity from smoothness control: - Use enough basis functions (n_basis) to capture the true function - Control smoothness entirely via the penalty parameter λ

For most applications, 20-30 basis functions with cubic degree and second-order penalty works well. The optimal λ is typically selected via GCV, AIC, or REML.

basis_matrix(x, sparse=False)

Evaluate B-spline basis at given points.

Parameters:

• x (array-like) – Points at which to evaluate basis functions

• sparse (bool, default=False) – Return sparse matrix (efficient for large n_basis)

Returns:

B – Basis matrix where B[i,j] = B_j(x[i])

Return type:

ndarray or sparse matrix, shape (len(x), n_basis)

penalty_matrix()

Return the difference penalty matrix S = D’D.

Returns:

S – Penalty matrix for penalized regression

Return type:

ndarray, shape (n_basis, n_basis)

Raises:

ValueError – If basis has not been set up yet (call basis_matrix first)

fit(x, y, lambda_='gcv', weights=None, lambda_range=(1e-6, 1e6))

Fit P-spline to data.

Parameters:

• x (array-like) – Predictor values

• y (array-like) – Response values

• lambda (float or {'gcv', 'aic', 'reml'}) – Smoothing parameter. If string, automatically select via: - ‘gcv’: Generalized Cross-Validation (default, fast) - ‘aic’: Akaike Information Criterion - ‘reml’: Restricted Maximum Likelihood (better for small samples)

• weights (array-like, optional) – Observation weights. If None, uniform weights.

• lambda_range (tuple of float) – Search range for automatic lambda selection

• lambda_ (float | str)

Returns:

Fitted P-spline with coefficients, fitted values, etc.

Return type:

PSplineResult

Examples

>>> # Automatic smoothing selection
>>> result = basis.fit(x, y, lambda_='gcv')
>>>
>>> # Fixed smoothing parameter
>>> result = basis.fit(x, y, lambda_=100.0)
>>>
>>> # With observation weights
>>> result = basis.fit(x, y, weights=w, lambda_='gcv')

derivative_matrix(x, order=1)

Compute derivative basis matrix.

Parameters:

• x (array-like) – Points at which to evaluate derivatives

• order (int) – Order of derivative (1, 2, …)

Returns:

dB – Derivative basis matrix

Return type:

ndarray, shape (len(x), n_basis)

class aurora.smoothing.splines.PSplineResult(coef_, fitted_values_, lambda_, basis, edf_, gcv_, sigma2_, x_, y_, weights_)

Bases: object

Results from P-spline fitting.

Variables:

• coef (ndarray) – Fitted spline coefficients

• fitted_values (ndarray) – Fitted values at training points

• lambda (float) – Smoothing parameter used

• basis (PSplineBasis) – The basis object used for fitting

• edf (float) – Effective degrees of freedom (trace of hat matrix)

• gcv (float) – Generalized Cross-Validation score

• sigma2 (float) – Estimated residual variance

• x (ndarray) – Training predictor values

• y (ndarray) – Training response values

• weights (ndarray) – Observation weights used

Parameters:

• coef_ (NDArray)

• fitted_values_ (NDArray)

• lambda_ (float)

• basis (PSplineBasis)

• edf_ (float)

• gcv_ (float)

• sigma2_ (float)

• x_ (NDArray)

• y_ (NDArray)

• weights_ (NDArray)

coef_: NDArray

fitted_values_: NDArray

lambda_: float

basis: PSplineBasis

edf_: float

gcv_: float

sigma2_: float

x_: NDArray

y_: NDArray

weights_: NDArray

predict(x_new)

Predict at new points.

Parameters:

x_new (array-like) – Points at which to predict

Returns:

y_pred – Predicted values

Return type:

ndarray

derivative(x=None, order=1)

Evaluate derivative of fitted spline.

Parameters:

• x (array-like, optional) – Points at which to evaluate. If None, uses training points.

• order (int) – Order of derivative (1, 2, …)

Returns:

deriv – Derivative values

Return type:

ndarray

confidence_band(x=None, level=0.95)

Compute pointwise confidence band.

Parameters:

• x (array-like, optional) – Points at which to compute band. If None, uses training points.

• level (float) – Confidence level (e.g., 0.95 for 95%)

Returns:

lower, upper – Lower and upper confidence bounds

Return type:

ndarray

summary()

Return summary statistics.

Returns:

Summary including lambda, edf, GCV, R², etc.

Return type:

dict

aurora.smoothing.splines.fit_pspline(x, y, n_basis=20, degree=3, penalty_order=2, lambda_='gcv', weights=None, **kwargs)

Fit P-spline smoother to data.

Convenience function that creates a PSplineBasis and fits it to data.

Parameters:

• x (array-like) – Predictor values

• y (array-like) – Response values

• n_basis (int, default=20) – Number of B-spline basis functions

• degree (int, default=3) – B-spline degree (cubic by default)

• penalty_order (int, default=2) – Order of difference penalty

• lambda (float or str, default='gcv') – Smoothing parameter or selection method (‘gcv’, ‘aic’, ‘reml’)

• weights (array-like, optional) – Observation weights

• **kwargs – Additional arguments passed to PSplineBasis.fit()

• lambda_ (float | str)

Returns:

Fitted P-spline

Return type:

PSplineResult

Examples

>>> import numpy as np
>>> from aurora.smoothing.splines import fit_pspline
>>>
>>> # Generate noisy sine data
>>> x = np.linspace(0, 2*np.pi, 100)
>>> y = np.sin(x) + np.random.normal(0, 0.2, 100)
>>>
>>> # Fit P-spline with automatic smoothing
>>> result = fit_pspline(x, y, n_basis=25)
>>> print(f"EDF: {result.edf_:.2f}, λ: {result.lambda_:.2f}")
>>>
>>> # Predict at new points
>>> x_new = np.linspace(0, 2*np.pi, 200)
>>> y_smooth = result.predict(x_new)

class aurora.smoothing.splines.bspline.BSplineBasis(knots, degree=3, cyclic=False)

Bases: object

B-spline basis functions via Cox-de Boor recursion.

B-splines are piecewise polynomials with compact support, defined recursively using the Cox-de Boor formula. They form a basis for the space of splines of a given degree with specified knots.

Parameters:

• knots (array-like) – Knot vector (should include boundary knots with multiplicity = degree + 1). For open B-splines, repeat boundary knots degree+1 times.

• degree (int, default=3) – Degree of the B-spline basis (3 for cubic B-splines).

• cyclic (bool, default=False) – Whether to use cyclic (periodic) boundary conditions.

Variables:

• knots (ndarray) – Full knot vector including repeated boundary knots.

• degree (int) – Degree of the splines.

• n_basis (int) – Number of basis functions.

• cyclic (bool) – Whether cyclic boundary conditions are used.

Notes

For degree p and n_basis basis functions, the knot vector must have length n_basis + p + 1.

References

de Boor, C. (1978). A Practical Guide to Splines. Springer-Verlag. Eilers, P.H.C. & Marx, B.D. (1996). Flexible smoothing with B-splines and

penalties. Statistical Science, 11(2), 89-121.

basis_matrix(x, sparse=False)

Compute B-spline basis matrix using Cox-de Boor recursion.

Parameters:

• x (array-like, shape (n_samples,)) – Points at which to evaluate basis functions.

• sparse (bool, default=False) – If True, return scipy.sparse.csr_matrix (only for NumPy backend). If False, return dense array compatible with input backend.

Returns:

B – Basis matrix where B[i,j] = B_j^p(x_i), the j-th B-spline of degree p evaluated at x_i.

Return type:

array or csr_matrix, shape (n_samples, n_basis)

Notes

Uses the Cox-de Boor recursion formula:

B_i^0(x) = 1 if t_i ≤ x < t_{i+1}, else 0 B_i^p(x) = w_i^p(x) B_i^{p-1}(x) + (1 - w_{i+1}^p(x)) B_{i+1}^{p-1}(x)

where w_i^p(x) = (x - t_i) / (t_{i+p} - t_i)

Sparse output exploits the compact support property: only (degree + 1) basis functions are non-zero at any point. For degree=3, each row has at most 4 non-zero entries, regardless of n_basis. This reduces: - Memory: O(n_samples × n_basis) → O(n_samples × degree) - Computation: O(n_samples × n_basis × degree²) → O(n_samples × degree²)

Sparse output is only available for NumPy backend. PyTorch and JAX backends always return dense arrays (sparse support is experimental).

penalty_matrix(order=2)

Compute difference penalty matrix.

For B-splines, we typically use a difference penalty that approximates the integrated squared derivative penalty.

Parameters:

order (int, default=2) – Order of differences (2 for approximating second derivative).

Returns:

S – Penalty matrix where β’Sβ approximates ∫ [f^(m)(x)]^2 dx.

Return type:

ndarray, shape (n_basis, n_basis)

Notes

The difference penalty is:

S = D’D

where D is the order-th difference matrix.

For order=2: D_ij penalizes (β_i - 2β_{i+1} + β_{i+2})² This approximates the integrated squared second derivative.

static create_knots(x, n_basis=10, degree=3, method='quantile')

Create knot vector for B-splines from data.

Parameters:

• x (array-like) – Data values.

• n_basis (int) – Number of basis functions desired.

• degree (int) – Degree of B-splines.

• method (str) – Knot placement method: - ‘quantile’: Place interior knots at quantiles - ‘uniform’: Place interior knots uniformly

Returns:

knots – Full knot vector with repeated boundary knots.

Return type:

ndarray

Notes

For open B-splines, boundary knots are repeated (degree + 1) times. This ensures the spline interpolates at boundaries.

derivative_basis_matrix(x, order=1)

Compute derivative of B-spline basis functions.

Parameters:

• x (array-like, shape (n_samples,)) – Evaluation points.

• order (int, default=1) – Order of derivative (1 for first derivative, etc.).

Returns:

dB – Matrix of derivative values.

Return type:

array, shape (n_samples, n_basis)

Notes

Uses the derivative formula:

dB_i^p(x)/dx = p * [B_i^{p-1}(x)/(t_{i+p}-t_i) -

B_{i+1}^{p-1}(x)/(t_{i+p+1}-t_{i+1})]

class aurora.smoothing.splines.cubic.CubicSplineBasis(knots, boundary_knots=None)

Bases: object

Natural cubic spline basis for smooth function representation.

Natural cubic splines provide a flexible basis for representing smooth functions. They are especially useful in GAMs for modeling non-linear relationships.

Parameters:

• knots (array-like) – Interior knot locations. Should be sorted in ascending order.

• boundary_knots (tuple of float, optional) – (lower, upper) boundary knots. If None, uses min and max of knots.

Variables:

• knots (ndarray) – Interior knots

• boundary_knots (tuple) – Boundary knots

• n_basis (int) – Number of basis functions (len(knots) + 2 for natural splines)

References

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

basis_matrix(x)

Compute the basis matrix B where B[i,j] = b_j(x_i).

Parameters:

x (array-like, shape (n_samples,)) – Points at which to evaluate the basis functions.

Returns:

B – Basis matrix. Uses the same array backend as input x.

Return type:

array, shape (n_samples, n_basis)

Notes

The basis functions are constructed to satisfy: - Cubic polynomials between knots - Continuous second derivatives at knots - Zero second derivatives at boundaries (natural condition)

penalty_matrix()

Compute the penalty matrix S for integrated squared second derivative.

The penalty matrix satisfies: β’Sβ = ∫ [f’’(x)]² dx

Returns:

S – Penalty matrix (always NumPy array).

Return type:

ndarray, shape (n_basis, n_basis)

Notes

For natural cubic splines, the penalty can be computed analytically. The first two basis functions (constant and linear) receive zero penalty since their second derivatives are zero.

For the ESL (5.4-5.5) basis N_{j+2} = d_j - d_{K-1} with d_j(x) = [(x - xi_j)_+^3 - (x - xi_K)_+^3] / (xi_K - xi_j), the second derivative is d_j’’(x) = 6[(x - xi_j)_+ - (x - xi_K)_+]/(xi_K - xi_j). Since (x - xi_K)_+ = 0 on [xi_1, xi_K] and every N_j’’ vanishes outside the boundary knots, the integral over the real line reduces to

∫ d_i’’ d_j’’ dx = J(xi_i, xi_j) / [(xi_K - xi_i)(xi_K - xi_j)]

with J(t, s) = ∫ 6(x - t)_+ 6(x - s)_+ dx, and

S[i+2, j+2] = D(i,j) - D(i,K-1) - D(K-1,j) + D(K-1,K-1).

static create_knots(x, n_knots=10, method='quantile')

Create knot locations from data.

Parameters:

• x (array-like) – Data values used to determine knot positions.

• n_knots (int, default=10) – Number of interior knots to place.

• method (str, default='quantile') –

  Knot placement method:

  • 'quantile': Place knots at quantiles of x, adapting to the empirical data distribution.

  • 'uniform': Place knots uniformly between min(x) and max(x).

Returns:

knots – Sorted interior knot locations. May contain fewer than n_knots entries if duplicate quantiles are merged.

Return type:

ndarray, shape (n_knots,) or fewer

Raises:

ValueError – If n_knots is less than 1 or method is unknown.

Examples

>>> x = np.random.randn(500)
>>> knots = CubicSplineBasis.create_knots(x, n_knots=8)
>>> len(knots) <= 8
True

class aurora.smoothing.splines.pspline.PSplineBasis(n_basis=20, degree=3, penalty_order=2, domain=None)

Bases: object

P-Spline basis: B-splines with equally-spaced knots and difference penalty.

P-Splines use a large number of B-spline basis functions with a difference penalty on coefficients to achieve smooth fits. The amount of smoothing is controlled by a single parameter λ.

Parameters:

• n_basis (int, default=20) – Number of B-spline basis functions. More basis functions allow more flexibility but require larger penalty. Typically 15-40 is sufficient.

• degree (int, default=3) – B-spline degree. Cubic (degree=3) is standard and provides C² continuity.

• penalty_order (int, default=2) – Order of difference penalty: - 1: Penalizes first differences (roughness of coefficients) - 2: Penalizes second differences (curvature, most common) - 3: Penalizes third differences (rate of curvature change)

• domain (tuple of float, optional) – (min, max) domain for the basis. If None, determined from data.

Variables:

• n_basis (int) – Number of basis functions

• degree (int) – B-spline degree

• penalty_order (int) – Order of difference penalty

• knots (ndarray) – Knot locations (set after fitting)

Examples

>>> import numpy as np
>>> from aurora.smoothing.splines import PSplineBasis
>>>
>>> # Create basis
>>> basis = PSplineBasis(n_basis=25, degree=3, penalty_order=2)
>>>
>>> # Generate data
>>> x = np.linspace(0, 2*np.pi, 100)
>>> y = np.sin(x) + np.random.normal(0, 0.2, 100)
>>>
>>> # Fit with automatic smoothing selection
>>> result = basis.fit(x, y, lambda_='gcv')
>>> y_smooth = result.fitted_values_

Notes

The P-spline approach separates the basis complexity from smoothness control: - Use enough basis functions (n_basis) to capture the true function - Control smoothness entirely via the penalty parameter λ

For most applications, 20-30 basis functions with cubic degree and second-order penalty works well. The optimal λ is typically selected via GCV, AIC, or REML.

knots_: NDArray | None

basis_matrix(x, sparse=False)

Evaluate B-spline basis at given points.

Parameters:

• x (array-like) – Points at which to evaluate basis functions

• sparse (bool, default=False) – Return sparse matrix (efficient for large n_basis)

Returns:

B – Basis matrix where B[i,j] = B_j(x[i])

Return type:

ndarray or sparse matrix, shape (len(x), n_basis)

penalty_matrix()

Return the difference penalty matrix S = D’D.

Returns:

S – Penalty matrix for penalized regression

Return type:

ndarray, shape (n_basis, n_basis)

Raises:

ValueError – If basis has not been set up yet (call basis_matrix first)

fit(x, y, lambda_='gcv', weights=None, lambda_range=(1e-6, 1e6))

Fit P-spline to data.

Parameters:

• x (array-like) – Predictor values

• y (array-like) – Response values

• lambda (float or {'gcv', 'aic', 'reml'}) – Smoothing parameter. If string, automatically select via: - ‘gcv’: Generalized Cross-Validation (default, fast) - ‘aic’: Akaike Information Criterion - ‘reml’: Restricted Maximum Likelihood (better for small samples)

• weights (array-like, optional) – Observation weights. If None, uniform weights.

• lambda_range (tuple of float) – Search range for automatic lambda selection

• lambda_ (float | str)

Returns:

Fitted P-spline with coefficients, fitted values, etc.

Return type:

PSplineResult

Examples

>>> # Automatic smoothing selection
>>> result = basis.fit(x, y, lambda_='gcv')
>>>
>>> # Fixed smoothing parameter
>>> result = basis.fit(x, y, lambda_=100.0)
>>>
>>> # With observation weights
>>> result = basis.fit(x, y, weights=w, lambda_='gcv')

derivative_matrix(x, order=1)

Compute derivative basis matrix.

Parameters:

• x (array-like) – Points at which to evaluate derivatives

• order (int) – Order of derivative (1, 2, …)

Returns:

dB – Derivative basis matrix

Return type:

ndarray, shape (len(x), n_basis)

Tensor Products

Tensor product smooths for modeling interactions between variables.

Tensor product smooths allow for flexible modeling of interactions between two or more variables. They construct a smooth surface by taking the tensor product of marginal basis functions.

References

Wood, S.N. (2017). Generalized Additive Models: An Introduction with R.

Chapman and Hall/CRC, 2nd edition. Chapter 5.

aurora.smoothing.tensor.tensor_product_basis(X1, X2, basis1, basis2)

Compute tensor product basis matrix.

Creates a tensor product basis by multiplying basis functions from two marginal bases: B(x₁, x₂) = B₁(x₁) ⊗ B₂(x₂)

Parameters:

• X1 (ndarray, shape (n,)) – First variable values.

• X2 (ndarray, shape (n,)) – Second variable values.

• basis1 (object) – First marginal basis (must have basis_matrix method).

• basis2 (object) – Second marginal basis (must have basis_matrix method).

Returns:

B – Tensor product basis matrix.

Return type:

ndarray, shape (n, n_basis1 * n_basis2)

Notes

For each observation i, the tensor product basis is:

B[i, :] = vec(B₁[i, :] ⊗ B₂[i, :])

where ⊗ denotes the outer product.

Examples

>>> import numpy as np
>>> from aurora.smoothing.splines.bspline import BSplineBasis
>>> from aurora.smoothing.tensor import tensor_product_basis
>>> # Create marginal bases
>>> x1 = np.linspace(0, 1, 100)
>>> x2 = np.linspace(0, 1, 100)
>>> knots1 = BSplineBasis.create_knots(x1, n_basis=8, degree=3)
>>> knots2 = BSplineBasis.create_knots(x2, n_basis=8, degree=3)
>>> basis1 = BSplineBasis(knots1, degree=3)
>>> basis2 = BSplineBasis(knots2, degree=3)
>>> # Compute tensor product
>>> B_tensor = tensor_product_basis(x1, x2, basis1, basis2)
>>> B_tensor.shape
(100, 64)

aurora.smoothing.tensor.tensor_product_penalty(S1, S2, p1, p2)

Compute tensor product penalty matrices.

For tensor products, we need two penalty matrices: - S_x1: penalizes variation in x₁ direction - S_x2: penalizes variation in x₂ direction

Parameters:

• S1 (ndarray, shape (p1, p1)) – Penalty matrix for first marginal basis.

• S2 (ndarray, shape (p2, p2)) – Penalty matrix for second marginal basis.

• p1 (int) – Number of basis functions in first marginal basis.

• p2 (int) – Number of basis functions in second marginal basis.

Returns:

• S_x1 (ndarray, shape (p1*p2, p1*p2)) – Penalty matrix for variation in x₁ direction.

• S_x2 (ndarray, shape (p1*p2, p1*p2)) – Penalty matrix for variation in x₂ direction.

Return type:

tuple[ndarray, ndarray]

Notes

The tensor product penalty is:

S_x1 = S1 ⊗ I2 (penalizes wiggliness in x₁) S_x2 = I1 ⊗ S2 (penalizes wiggliness in x₂)

where ⊗ is the Kronecker product.

The full penalty is: λ₁ β’ S_x1 β + λ₂ β’ S_x2 β

Examples

>>> import numpy as np
>>> from aurora.smoothing.tensor import tensor_product_penalty
>>> # Create marginal penalties
>>> S1 = np.eye(5)
>>> S2 = np.eye(6)
>>> # Compute tensor penalties
>>> S_x1, S_x2 = tensor_product_penalty(S1, S2, 5, 6)
>>> S_x1.shape
(30, 30)
>>> S_x2.shape
(30, 30)

aurora.smoothing.tensor.fit_tensor_product(X1, X2, y, basis1, basis2, S1, S2, lambda1=1.0, lambda2=1.0, weights=None)

Fit a tensor product smooth.

Parameters:

• X1 (ndarray, shape (n,)) – First variable values.

• X2 (ndarray, shape (n,)) – Second variable values.

• y (ndarray, shape (n,)) – Response variable.

• basis1 (object) – First marginal basis.

• basis2 (object) – Second marginal basis.

• S1 (ndarray, shape (p1, p1)) – Penalty matrix for first marginal basis.

• S2 (ndarray, shape (p2, p2)) – Penalty matrix for second marginal basis.

• lambda1 (float, default=1.0) – Smoothing parameter for x₁ direction.

• lambda2 (float, default=1.0) – Smoothing parameter for x₂ direction.

• weights (ndarray, shape (n,), optional) – Observation weights.

Returns:

result – Dictionary containing: - ‘coefficients’: Coefficient estimates - ‘fitted_values’: Fitted values - ‘basis_matrix’: Tensor product basis matrix - ‘edf’: Effective degrees of freedom

Return type:

dict

Notes

Minimizes: ||W^(1/2)(y - Bβ)||² + λ₁ β’ S_x1 β + λ₂ β’ S_x2 β

where B is the tensor product basis matrix.

Examples

>>> import numpy as np
>>> from aurora.smoothing.splines.bspline import BSplineBasis
>>> from aurora.smoothing.tensor import fit_tensor_product
>>> # Generate data
>>> n = 200
>>> x1 = np.random.uniform(0, 1, n)
>>> x2 = np.random.uniform(0, 1, n)
>>> y = np.sin(2*np.pi*x1) * np.cos(2*np.pi*x2) + 0.1*np.random.randn(n)
>>> # Create bases
>>> knots1 = BSplineBasis.create_knots(x1, n_basis=10, degree=3)
>>> knots2 = BSplineBasis.create_knots(x2, n_basis=10, degree=3)
>>> basis1 = BSplineBasis(knots1, degree=3)
>>> basis2 = BSplineBasis(knots2, degree=3)
>>> S1 = basis1.penalty_matrix(order=2)
>>> S2 = basis2.penalty_matrix(order=2)
>>> # Fit tensor product
>>> result = fit_tensor_product(x1, x2, y, basis1, basis2, S1, S2)
>>> result['fitted_values'].shape
(200,)

Penalties

Penalties for smoothing parameter estimation.

Smoothing Parameter Selection

Smoothing parameter selection strategies.

aurora.smoothing.selection.gcv_score(y, X, S, lambda_, weights=None)

Compute GCV score for given smoothing parameter.

The GCV score is:

GCV(λ) = (n_eff * RSS) / (n_eff - tr(H))²

where RSS is residual sum of squares and H is the hat matrix.

Parameters:

• y (ndarray, shape (n,)) – Response variable.

• X (ndarray, shape (n, p)) – Design matrix (basis matrix).

• S (ndarray, shape (p, p)) – Penalty matrix.

• lambda (float) – Smoothing parameter.

• weights (ndarray, shape (n,), optional) – Observation weights (for weighted least squares).

• lambda_ (float)

Returns:

score – GCV score (lower is better).

Return type:

float

Notes

The penalized least squares problem is:

min_β ||y - Xβ||² + λ β’Sβ

Solution: β̂ = (X’WX + λS)⁻¹ X’Wy Hat matrix: H = X(X’WX + λS)⁻¹X’W

GCV approximates leave-one-out CV by:

GCV = RSS / (1 - tr(H)/n)²

With observation weights, n is replaced by the effective sample size n_eff = #{i : w_i > 0} (Wood 2017, §4.5): zero-weight observations contribute neither to RSS nor to tr(H), so counting them in n would bias the score towards undersmoothing.

References

Craven, P. & Wahba, G. (1978). Smoothing noisy data with spline functions.

Numerische Mathematik, 31, 377-403. https://doi.org/10.1007/BF01404567

Wood, S.N. (2017). Generalized Additive Models: An Introduction with R.

aurora.smoothing.selection.select_smoothing_parameter(y, X, S, weights=None, lambda_min=1e-6, lambda_max=1e6, method='golden')

Select optimal smoothing parameter using GCV.

Parameters:

• y (ndarray, shape (n,)) – Response variable.

• X (ndarray, shape (n, p)) – Design matrix (basis matrix).

• S (ndarray, shape (p, p)) – Penalty matrix.

• weights (ndarray, shape (n,), optional) – Observation weights.

• lambda_min (float, default=1e-6) – Minimum smoothing parameter to consider.

• lambda_max (float, default=1e6) – Maximum smoothing parameter to consider.

• method (str, default='golden') – Optimization method: - ‘golden’: Golden section search (robust, no derivatives) - ‘brent’: Brent’s method (faster, still no derivatives)

Returns:

result – Dictionary containing: - ‘lambda_opt’: Optimal smoothing parameter - ‘gcv_score’: GCV score at optimum - ‘edf’: Effective degrees of freedom at optimum - ‘coefficients’: Fitted coefficients at optimum - ‘fitted_values’: Fitted values at optimum

Return type:

dict

Notes

GCV tends to undersmooth slightly compared to REML/ML, but is computationally simpler and doesn’t require distributional assumptions.

For multiple smoothing parameters (multi-dimensional λ), use REML instead.

Examples

>>> from aurora.smoothing.splines.bspline import BSplineBasis
>>> x = np.linspace(0, 1, 100)
>>> y = np.sin(2 * np.pi * x) + 0.1 * np.random.randn(100)
>>> knots = BSplineBasis.create_knots(x, n_basis=10, degree=3)
>>> basis = BSplineBasis(knots, degree=3)
>>> X = basis.basis_matrix(x)
>>> S = basis.penalty_matrix(order=2)
>>> result = select_smoothing_parameter(y, X, S)
>>> print(f"Optimal λ: {result['lambda_opt']:.4f}")
>>> print(f"Effective DoF: {result['edf']:.2f}")

aurora.smoothing.selection.reml_score(y, X, S, lambda_, weights=None)

Compute REML score for given smoothing parameter.

The REML criterion (Wood 2011, eq. 4; Wood 2017, §4.6), profiled over the scale parameter and up to additive constants in lambda, is:

-2 log(REML) = (n - M_p) log(D) + log|X’WX + λS| - r log(λ) - log|S|_+

where D is the penalised deviance ||W^{1/2}(y - Xβ̂)||² + λβ̂'Sβ̂, r = rank(S), M_p = dim null(S), and log|S|_+ is the log pseudo-determinant of S.

Parameters:

• y (ndarray, shape (n,)) – Response variable.

• X (ndarray, shape (n, p)) – Design matrix (basis functions).

• S (ndarray, shape (p, p)) – Penalty matrix.

• lambda (float) – Smoothing parameter (must be positive).

• weights (ndarray, shape (n,), optional) – Observation weights.

• lambda_ (float)

Returns:

reml – REML score (lower is better).

Return type:

float

Notes

The REML criterion penalizes both overfitting (via the deviance) and model complexity (via determinant terms). It is invariant to fixed effects and provides unbiased estimates of variance components.

The -r log(λ) term acts as a barrier preventing λ → 0 (interpolation); omitting it collapses the optimum onto the lower search bound.

REML is generally preferred over GCV for: - Multiple smoothing parameters - Small sample sizes - Models with fixed effects

aurora.smoothing.selection.select_smoothing_parameter_reml(y, X, S, weights=None, lambda_min=1e-6, lambda_max=1e6, method='brent')

Select smoothing parameter via REML optimization.

Minimizes the REML criterion over the smoothing parameter λ.

Parameters:

• y (ndarray, shape (n,)) – Response variable.

• X (ndarray, shape (n, p)) – Design matrix (basis functions).

• S (ndarray, shape (p, p)) – Penalty matrix.

• weights (ndarray, shape (n,), optional) – Observation weights.

• lambda_min (float, default=1e-6) – Minimum smoothing parameter to search.

• lambda_max (float, default=1e6) – Maximum smoothing parameter to search.

• method ({'brent', 'golden'}, default='brent') –

  Optimization method: - ‘brent’: Bounded Brent search on the log(λ) interval

  (scipy.optimize.minimize_scalar(method='bounded'); recommended)

  • ’golden’: Golden section search over a bracket spanning the bounds

Returns:

result – Dictionary containing: - ‘lambda_opt’: Optimal smoothing parameter - ‘reml_score’: REML score at optimum - ‘coefficients’: Coefficient estimates - ‘fitted_values’: Fitted values - ‘edf’: Effective degrees of freedom

Return type:

dict

Examples

>>> import numpy as np
>>> from aurora.smoothing.selection.reml import select_smoothing_parameter_reml
>>> from aurora.smoothing.splines.bspline import BSplineBasis
>>> # Generate data
>>> x = np.linspace(0, 1, 100)
>>> y = np.sin(2 * np.pi * x) + 0.1 * np.random.randn(100)
>>> # Create basis
>>> knots = BSplineBasis.create_knots(x, n_basis=10, degree=3)
>>> basis = BSplineBasis(knots, degree=3)
>>> X = basis.basis_matrix(x)
>>> S = basis.penalty_matrix(order=2)
>>> # Select lambda via REML
>>> result = select_smoothing_parameter_reml(y, X, S)
>>> print(f"Optimal lambda: {result['lambda_opt']:.4f}")
>>> print(f"EDF: {result['edf']:.2f}")

Notes

REML is preferred over GCV when: - You have multiple smoothing parameters - Sample size is small - You want theoretical justification (ML-based)

The optimization is performed on log(λ) scale for numerical stability.

References

Wood, S.N. (2017). Generalized Additive Models: An Introduction with R.

Chapman and Hall/CRC, 2nd edition.

aurora.smoothing.selection.select_multiple_smoothing_parameters_reml(y, X_list, S_list, weights=None, lambda_init=None, lambda_min=1e-6, lambda_max=1e6, max_iter=20, tol=1e-4)

Select multiple smoothing parameters via REML optimization.

Uses alternating optimization: fix all λⱼ except one, optimize that one, repeat until convergence.

Parameters:

• y (ndarray, shape (n,)) – Response variable.

• X_list (list of ndarray) – List of design matrices, one per smooth term.

• S_list (list of ndarray) – List of penalty matrices, one per smooth term.

• weights (ndarray, shape (n,), optional) – Observation weights.

• lambda_init (list of float, optional) – Initial smoothing parameters. If None, starts with all 1.0.

• lambda_min (float, default=1e-6) – Minimum smoothing parameter.

• lambda_max (float, default=1e6) – Maximum smoothing parameter.

• max_iter (int, default=20) – Maximum number of alternating iterations.

• tol (float, default=1e-4) – Convergence tolerance for relative change in λ.

Returns:

result – Dictionary containing: - ‘lambda_opt’: List of optimal smoothing parameters - ‘reml_score’: REML score at optimum - ‘coefficients’: Full coefficient vector - ‘fitted_values’: Fitted values - ‘edf_values’: List of EDF per smooth term - ‘converged’: Whether optimization converged - ‘n_iter’: Number of iterations

Return type:

dict

Notes

This implements a simple alternating optimization scheme. More sophisticated methods exist (e.g., Newton-Raphson on all λ simultaneously) but are more complex to implement.

The alternating scheme is: 1. Fix λ₁, …, λⱼ₋₁, λⱼ₊₁, …, λₘ 2. Optimize λⱼ via REML 3. Cycle through j = 1, …, m 4. Repeat until convergence

References

Wood, S.N. (2011). Fast stable restricted maximum likelihood and marginal

likelihood estimation of semiparametric generalized linear models.