# DuchonSpline


Duchon spline basis.


Usage

``` python
DuchonSpline(
    k=10,
    m=2,
)
```


Duchon splines (Duchon, 1977) generalize thin plate splines by decoupling the radial basis exponent from the covariate dimension. Ordinary TPRS ties the exponent of its radial kernel to both the derivative order `m` being penalized and the covariate dimension [d](TPRS.md#whittaker.TPRS.d) (exponent `2m - d`), which means that for high-dimensional covariates the derivative order actually penalized can end up being uncomfortably high just to keep the kernel well-defined. [DuchonSpline](DuchonSpline.md#whittaker.DuchonSpline) introduces an independent exponent parameter `s`, so the radial kernel and the polynomial null-space order can be chosen separately. Like TPRS, it is built as a low-rank eigen-approximation to the full spline (Wood 2003's construction, generalized to the Duchon kernel), so it requires no knot placement and works for any covariate dimension. Choose [DuchonSpline](DuchonSpline.md#whittaker.DuchonSpline) over [TPRS](TPRS.md#whittaker.TPRS) when you want explicit control over the smoothness/exponent trade-off independent of dimension -- for example to match a specific derivative-penalty order in higher dimensions without also inflating the null-space order.


## Parameters


`k: int = ``10`  
Total number of basis functions, including the `M` polynomial null-space columns. Must satisfy `k > M` where `M = C(m_order - 1 + d, d)`. Larger `k` allows more wiggly fits at the cost of more computation; the roughness penalty (not `k`) ultimately controls smoothness once `lambda` is chosen. The default is `10`.

`m: int | list | tuple = ``2`  
Order specification, either:

- a single integer, interpreted as the polynomial null-space order `m_order`, with the radial exponent parameter `s` defaulting to `1.0`; or
- a two-element list/tuple `[s, m_order]`, where `s >= 0` is the real-valued radial basis exponent (the kernel behaves like `r^(2s)`, optionally with a log factor) and `m_order >= 1` is the polynomial null-space order (the null space consists of all monomials of total degree `<= m_order - 1`).

Setting `s = m_order - d / 2` for integer `m_order` recovers the ordinary TPRS basis for that order. The default is `2` (i.e. `s=1.0`, `m_order=2`).


## Notes

The Duchon radial kernel is

 \eta_s(r) = \begin{cases} r^{2s} & 2s \text{ is not an even integer} \\ r^{2s} \log(r) & 2s \text{ is an even integer} \end{cases}, 

evaluated at pairwise distances `r = ||x_i - x_j||`, with the convention `η(0) = 0`. Together with the polynomial null space of all monomials of total degree at most `m_order - 1` (dimension `M = C(m_order - 1 + d, d)`), the basis is constructed in the same two stages as [TPRS](TPRS.md#whittaker.TPRS):

1.  **Polynomial null space** (first `M` columns): unpenalized low-degree polynomials.
2.  **Truncated spline part** (remaining `k - M` columns): the full `n x n` kernel matrix is projected onto the orthogonal complement of the null space (via a QR decomposition of the null-space design matrix) and eigendecomposed; the `k - M` leading eigenvectors give the best rank-`(k - M)` approximation to the full Duchon spline for that basis dimension.

The resulting penalty matrix is block-diagonal,

 \mathbf{S} = \operatorname{diag}(0, \ldots, 0, \lambda_1, \ldots, \lambda\_{k-M}), 

with the first `M` rows/columns exactly zero (unpenalized null space) and the remainder equal to the retained eigenvalues of the projected kernel matrix. As with [TPRS](TPRS.md#whittaker.TPRS), columns of `x` with very different scales can cause numerical issues in the eigendecomposition, so centering and/or standardizing covariates before fitting is advisable. Non-integer or large `s` values can also make the kernel matrix increasingly ill-conditioned; if fitting becomes numerically unstable, try a smaller `s` or standardized covariates.


## Examples


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

rng = np.random.default_rng(0)
x = rng.uniform(0, 1, 100)

basis = DuchonSpline(k=10, m=[1.0, 2]).fit(x)
B = basis.basis_matrix(x)
S = basis.penalty_matrix()
B.shape, S.shape
```


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


## Attributes

| Name | Description |
|----|----|
| [is_fitted](#is_fitted) | Whether the basis has been fitted. |
| [n_basis](#n_basis) | Total number of basis functions. |

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


### is_fitted


Whether the basis has been fitted.


`is_fitted: bool`


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


### n_basis


Total number of basis functions.


`n_basis: int`


Equal to `k`: the `M` polynomial null-space columns plus the `k - M` truncated-spline columns retained from the eigendecomposition of the projected Duchon kernel matrix.


## Methods

| Name | Description |
|----|----|
| [basis_matrix()](#basis_matrix) | Evaluate the Duchon spline basis at `x`. |
| [fit()](#fit) | Fit the Duchon 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 `M`, the dimension of the polynomial null space. |
| [penalty_matrix()](#penalty_matrix) | Return the `k x k` penalty matrix `S`. |

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


### basis_matrix()


Evaluate the Duchon spline basis at `x`.


Usage

``` python
basis_matrix(x)
```


#### Parameters


`x: NDArray`  
Covariate values. Shape `(n,)` or `(n, d)` where [d](TPRS.md#whittaker.TPRS.d) must match the training dimension.


#### Returns


`NDArray`  
Design matrix of shape `(n, k)`. The first `M` columns are the polynomial null-space functions; the remaining `k - M` columns are the truncated spline functions.


#### Raises


`ValueError`  
If the covariate dimension of `x` does not match the training dimension.


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


### fit()


Fit the Duchon spline to training data `x`.


Usage

``` python
fit(x)
```


#### Parameters


`x: NDArray`  
Training covariates. Shape `(n,)` for univariate or `(n, d)` for multivariate.


#### Returns


`DuchonSpline`  
Returns `self` for method chaining.


#### Raises


`ValueError`  
If `k` is too small for the null-space dimension `M`, or too large for the number of observations.


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


### identifiability_constraints()


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


Usage

``` python
identifiability_constraints()
```


#### Returns


`NDArray`  
A `(1, k)` matrix whose product with the coefficient vector is zero when the smooth has mean zero over the training data.


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


### null_space_dimension()


Return `M`, the dimension of the polynomial null space.


Usage

``` python
null_space_dimension()
```


#### Returns


`int`  
`M = C(m_order - 1 + d, d)`, the number of monomials of total degree at most `m_order - 1` in [d](TPRS.md#whittaker.TPRS.d) variables.


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


### penalty_matrix()


Return the `k x k` penalty matrix `S`.


Usage

``` python
penalty_matrix()
```


`S` is block-diagonal: the first `M` rows/columns (the polynomial null space) are exactly zero, and the remaining `k - M` rows/columns hold the retained eigenvalues of the projected Duchon kernel matrix on the diagonal.


#### Returns


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