# PSpline


P-Spline: B-spline basis with m-th order difference penalty.


Usage

``` python
PSpline(
    k=10,
    degree=3,
    m=2,
)
```


A P-spline (Eilers & Marx 1996) combines a rich B-spline basis on equally-spaced knots with a discrete difference penalty directly on the coefficients, rather than an integrated-derivative penalty on the fitted curve. This decouples basis richness from smoothness -- you can use many more basis functions than you would dare with an unpenalized spline, because the difference penalty (together with an appropriately chosen `lambda`) does the work of controlling wiggliness. Equivalent to mgcv's `bs="ps"` basis. P-splines are a good default choice for a single continuous covariate: they are cheap to construct (no eigendecomposition or dense linear solve is needed to build the penalty), the penalty matrix is banded/sparse which keeps large-`k` fits fast, and B-splines extrapolate smoothly beyond the training range via their boundary basis functions. Unlike [CRS](CRS.md#whittaker.CRS), knots are equidistant rather than placed at data quantiles, so P-splines are somewhat less efficient when the data are very unevenly distributed but are simpler and faster to set up.


## Parameters


`k: int = ``10`  
Number of B-spline basis functions. Must satisfy `k >= degree + 1`. Larger `k` gives a richer basis (more, narrower B-splines) and lets the fit follow finer local structure; as with other penalized bases, wiggliness is ultimately governed by the penalty and `lambda`, not by `k` alone, so `k` can usually be set generously (e.g. `k=20-40`) without much downside. The default is `10`.

`degree: int = ``3`  
Polynomial degree of each B-spline piece. `degree=3` (cubic) is the conventional choice and matches most GAM software; `degree=1` gives a piecewise-linear basis useful for less smooth phenomena, `degree=0` a step-function basis. The default is `3` (cubic).

`m: int = ``2`  
Order of the difference penalty applied to adjacent B-spline coefficients. `m=2` (the default) penalizes second differences, which is the discrete analogue of penalizing curvature and is by far the most common choice; `m=1` penalizes changes in level (shrinks toward a constant), `m=3` penalizes changes in slope-of-slope for extra-smooth fits. Must satisfy `0 < m < k`. The default is `2` (second differences: penalizes curvature).


## Notes

The knot vector is built by `_bspline_knots()`: `degree + 1` repeated (clamped) knots at each of `x_min` and `x_max`, plus `k - degree - 1` interior knots equally spaced between them, giving `k` B-spline basis functions of the requested `degree` via `scipy`'s `BSpline` machinery (de Boor's algorithm). Because the knots are equally spaced rather than data-adaptive, the design matrix `basis_matrix(x)` can be evaluated in closed form for any `x`, including points beyond `[x_min, x_max]`, which the boundary B-splines extend smoothly.

The penalty acts directly on the coefficient vector `\boldsymbol{\beta}` through the `m`-th order finite-difference matrix `\mathbf{D}_m` (shape `(k - m, k)`, built by applying `numpy.diff` `m` times to the identity):

 \mathbf{S} = \mathbf{D}\_m^\top \mathbf{D}\_m, 

a `k \times k` positive semi-definite matrix of rank `k - m`. Its `m`-dimensional null space is spanned by the discrete polynomial sequences `[1, 1, \ldots, 1]`, `[0, 1, \ldots, k-1]`, …, up to degree `m - 1` in the coefficient index -- i.e. coefficient vectors that are themselves polynomial in index have zero penalty, mirroring how polynomials of degree `< m` have zero `m`-th derivative in the continuous case. Because `S` is banded (bandwidth `m`), it is sparse and cheap to factorize even for large `k`; the main numerical caveat is that `basis_matrix()` clips evaluation points to the B-spline's knot support before calling `scipy`'s `BSpline.design_matrix`, since points exactly at or beyond the padded boundary can otherwise trigger an out-of-support error.


## Examples


``` python
import numpy as np
from whittaker.smooths import PSpline

x = np.linspace(0, 1, 100)
basis = PSpline(k=10).fit(x)
B = basis.basis_matrix(x)
S = basis.penalty_matrix()
B.shape, S.shape
```


    ((100, 10), (10, 10))


## Attributes

| Name | Description |
|----|----|
| [interior_knots](#interior_knots) | Interior knot locations only. |
| [is_fitted](#is_fitted) | Whether the basis has been fitted. |
| [knots](#knots) | Full augmented B-spline knot vector. |
| [n_basis](#n_basis) | Total number of B-spline basis functions. |

------------------------------------------------------------------------


### interior_knots


Interior knot locations only.


`interior_knots: NDArray`


Excludes the repeated boundary knots at `x_min` and `x_max`, leaving just the `k - degree - 1` equally-spaced knots strictly between them.


------------------------------------------------------------------------


### is_fitted


Whether the basis has been fitted.


`is_fitted: bool`


------------------------------------------------------------------------


### knots


Full augmented B-spline knot vector.


`knots: NDArray`


Includes the `degree + 1` repeated boundary knots at each end (needed for the clamped B-spline construction) as well as the equally-spaced interior knots.


------------------------------------------------------------------------


### n_basis


Total number of B-spline basis functions.


`n_basis: int`


Equal to `k`, the requested basis dimension: every B-spline coefficient column is retained (no null-space or wrap-around columns are dropped, unlike some of the other smooth bases).


## Methods

| Name | Description |
|----|----|
| [basis_matrix()](#basis_matrix) | Evaluate the B-spline basis at `x`. |
| [fit()](#fit) | Fit the P-spline to training data `x`. |
| [identifiability_constraints()](#identifiability_constraints) | Return the sum-to-zero constraint row for the intercept. |
| [null_space_dimension()](#null_space_dimension) | Return the dimension of the penalty null space. |
| [penalty_matrix()](#penalty_matrix) | Return the `k x k` penalty matrix `S = D_m' D_m`. |

------------------------------------------------------------------------


### basis_matrix()


Evaluate the B-spline basis at `x`.


Usage

``` python
basis_matrix(x)
```


Values outside the training range `[x_min, x_max]` are extrapolated by the boundary B-spline functions (the basis smoothly extends beyond the knot range).


#### Parameters


`x: NDArray`  
1-D evaluation points, shape `(n,)` or `(n, 1)`.


#### Returns


`NDArray`  
Design matrix of shape `(n, k)`.


------------------------------------------------------------------------


### fit()


Fit the P-spline to training data `x`.


Usage

``` python
fit(x)
```


Determines the knot vector from the data range and pre-computes the penalty matrix. No per-observation computation is stored as the B-spline design matrix is always evaluated on demand.


#### Parameters


`x: NDArray`  
1-D training covariate, shape `(n,)` or `(n, 1)`.


#### Returns


`PSpline`  
Returns `self` for method chaining.


#### Raises


`ValueError`  
If `n < k`, or if `x` is not 1-D / `(n, 1)`.


------------------------------------------------------------------------


### identifiability_constraints()


Return the sum-to-zero constraint row for the intercept.


Usage

``` python
identifiability_constraints()
```


Evaluates the basis on a fine equally-spaced grid spanning the training range and averages each column, giving a `(1, k)` row `C` such that `C @ beta == 0` constrains the smooth to have mean zero over the training range.


#### Returns


`NDArray`  
Shape `(1, k)` matrix that enforces mean-zero contribution over the training knot range.


------------------------------------------------------------------------


### null_space_dimension()


Return the dimension of the penalty null space.


Usage

``` python
null_space_dimension()
```


For a P-spline with difference order `m`, the null space of `D_m^T D_m` is exactly `m`-dimensional: it is spanned by the coefficient sequences that are themselves polynomial in the coefficient index up to degree `m - 1`, since these have zero `m`-th finite difference.


#### Returns


`int`  
The difference order `m`.


------------------------------------------------------------------------


### penalty_matrix()


Return the `k x k` penalty matrix `S = D_m' D_m`.


Usage

``` python
penalty_matrix()
```


`S` is positive semi-definite with rank `k - m`. Its null space is spanned by the `m` polynomial sequences `[1, 1, ..., 1]`, `[0, 1, ..., k-1]`, …, `[0^(m-1), 1^(m-1), ..., (k-1)^(m-1)]`.


#### Returns


`NDArray`  
Shape `(k, k)`.