fix(spline): correct penalty matrices and validate diff_order

docs/getting_started/quickstart.md

@@ -111,7 +111,7 @@ x = np.random.randn(200, 1)
 spline = NaturalCubicSplineTransformer(output_dim=8)
 spline.fit_transform(x)
 
-penalty = spline.get_penalty_matrix()   # second-difference penalty for GAM-style fitting
+penalty = spline.get_penalty_matrix()   # integrated-curvature penalty for GAM-style fitting
 ```
 
 The multivariate thin-plate spline models several columns jointly and is sized by

docs/representations/comparison_table.md

@@ -48,12 +48,15 @@ source of truth, and these tables mirror it.
 | Natural cubic spline | `naturalspline` | univariate | optional | yes | yes | yes |
 | Penalized spline (P-spline) | `pspline` | univariate | forbidden | no | yes | yes |
 | Tensor-product spline | `tensorspline` | multivariate | forbidden | no | yes | no |
-| Thin-plate spline | `tprs` | multivariate | forbidden | no | yes | no |
+| Thin-plate spline | `tprs` | multivariate | forbidden | no | experimental | no |
 
 ```{note}
 The multivariate splines (`tensorspline`, `tprs`) model several inputs jointly and are used
 standalone, not selected per column through `Preprocessor`. The alias `thinplate` resolves to
-`tprs`.
+`tprs`. `ThinPlateSplineTransformer.get_penalty_matrix()` is experimental: it is not guaranteed
+positive semi-definite (the retained eigenvalues of the projected landmark kernel can be
+negative) and emits a `ConfigWarning` on every call. `transform()` is unaffected; only the
+penalty is experimental.
 ```
 
 ## Functional expansions

docs/representations/spline_expansions.md

@@ -270,6 +270,14 @@ t.fit_transform(X2).shape   # (200, 10): output width is n_components, not input
 Constructor highlights: `n_components=10`, `landmark_strategy="kmeans"`, `rank_strategy="eigen"`,
 `include_bias=False`, `random_state`.
 
+```{warning}
+Unlike the other spline families on this page, `get_penalty_matrix()` on
+`ThinPlateSplineTransformer` is experimental: the retained eigenvalues are not guaranteed
+non-negative, so the returned penalty is not guaranteed positive semi-definite, and a
+`ConfigWarning` is raised on every call. `transform()` is unaffected; avoid this penalty for
+actual smoothing regularization until it is fully fixed.
+```
+
 ```{warning}
 The tensor-product and thin-plate splines are multivariate. They are standalone transformers
 and are not available as a per-column `numerical_method`. Fit them directly on the columns you

pretab/expansion/spline/cubic_regression.py

@@ -1,3 +1,5 @@
+import itertools
+
 import numpy as np
 from sklearn.base import BaseEstimator, TransformerMixin
 from sklearn.utils.validation import check_is_fitted
@@ -179,13 +181,18 @@ def fit(self, X, y=None):
 
         self.knots_ = []
         self.n_basis_ = []
+        self.x_min_ = []
+        self.x_max_ = []
         for i in range(X.shape[1]):
             xi = X[:, i]
             knots = self._place_interior_knots(
                 xi, y, n_interior, strategy, selector, self.task, min_interior, max_interior
             )
             self.knots_.append(knots)
             self.n_basis_.append(self._bspline_basis(xi, knots).shape[1])
+            finite_xi, _ = self._finite_column(xi, None)
+            self.x_min_.append(float(finite_xi.min()))
+            self.x_max_.append(float(finite_xi.max()))
 
         self.n_knots_ = [len(knots) for knots in self.knots_]
         return self
@@ -208,6 +215,11 @@ def fit_transform(self, X, y=None):
     def get_penalty_matrix(self, feature_index=0):
         """Return the curvature penalty matrix for a fitted feature.
 
+        Penalizes the integrated squared second derivative of every basis
+        column over the fitted feature range, including the polynomial
+        ``x**2`` / ``x**3`` columns (only ``x`` and the bias have an
+        identically-zero second derivative and so are unpenalized).
+
         Parameters
         ----------
         feature_index : int, default=0
@@ -221,18 +233,36 @@ def get_penalty_matrix(self, feature_index=0):
         """
         check_is_fitted(self, "n_basis_")
         n_basis = self.n_basis_[feature_index]
+        knots = self.knots_[feature_index]
+        x_min = self.x_min_[feature_index]
+        x_max = self.x_max_[feature_index]
+
+        # Every basis column's second derivative is piecewise linear in x, with
+        # kinks only at the knots, so a single-panel Simpson's rule per
+        # knot-delimited segment integrates every pairwise product exactly.
+        breakpoints = np.unique(np.concatenate(([x_min, x_max], np.clip(knots, x_min, x_max))))
+        breakpoints.sort()
+
         P = np.zeros((n_basis, n_basis))
-        offset = 4 if self.include_bias else 3
-        for i in range(offset, n_basis):
-            for j in range(offset, n_basis):
-                ki = self.knots_[feature_index][i - offset]
-                kj = self.knots_[feature_index][j - offset]
-                P[i, j] = self._spline_penalty_entry(ki, kj, self.knots_[feature_index])
+        for lo, hi in itertools.pairwise(breakpoints):
+            if hi <= lo:
+                continue
+            xs = np.array([lo, 0.5 * (lo + hi), hi])
+            weights = (hi - lo) / 6.0 * np.array([1.0, 4.0, 1.0])
+            D = np.column_stack([self._second_derivative(xs, i, knots) for i in range(n_basis)])
+            P += (D * weights[:, None]).T @ D
         return P
 
-    def _spline_penalty_entry(self, ki, kj, knots):
-        kmax = max(ki, kj)
-        upper = knots[-1]
-        x_vals = np.linspace(kmax, upper, 100)
-        integrand = 36 * (x_vals - ki) * (x_vals - kj)
-        return np.trapezoid(integrand, x_vals)
+    def _second_derivative(self, x, basis_index, knots):
+        """Second derivative of one basis column, evaluated at ``x``."""
+        poly_offset = 1 if self.include_bias else 0
+        col = basis_index - poly_offset
+        if col <= 0:
+            return np.zeros_like(x, dtype=float)  # bias / x: identically zero
+        if col == 1:
+            return np.full_like(x, 2.0, dtype=float)  # x**2
+        if col == 2:
+            return 6.0 * x  # x**3
+        knot = knots[col - 3]
+        return 6.0 * np.maximum(x - knot, 0.0)
+

pretab/expansion/spline/multivariate/tensor_product.py

@@ -42,6 +42,8 @@ class TensorProductSplineTransformer(SplineBasisMixin, TransformerMixin, BaseEst
 
     diff_order : int, default=2
         Order of the finite difference penalty used to enforce smoothness along each input dimension.
+        Must be a positive integer and small enough that the marginal basis width supports it
+        (``diff_order <= output_dim - 1``, roughly); both are validated at ``fit``.
 
     include_bias : bool, default=False
         If True, prepend a constant column to each marginal basis before the
@@ -177,6 +179,8 @@ def fit(self, X, y=None):
                 f"output_dim must be >= degree + 1 = {self.degree + 1} for the tensor-product "
                 f"spline basis, got {output_dim}"
             )
+        if not isinstance(self.diff_order, (int, np.integer)) or self.diff_order < 1:
+            raise InvalidParamError(f"diff_order must be a positive integer (>= 1); got {self.diff_order!r}.")
 
         self.dim_ = X.shape[1]
         self.knots_ = []
@@ -193,7 +197,15 @@ def fit(self, X, y=None):
                 X[:, d], y, output_dim, self.degree, strategy, None, None, min_interior, max_interior
             )
             basis = self._basis_matrix(X[:, d], knots)
-            penalty = self._difference_penalty(len(knots) - self.degree - 1)
+            n_basis = len(knots) - self.degree - 1
+            if self.diff_order > n_basis - 1:
+                raise InvalidParamError(
+                    f"diff_order={self.diff_order} is too large for a marginal basis of width "
+                    f"{n_basis} (output_dim={output_dim}); the finite-difference penalty needs "
+                    f"diff_order <= {n_basis - 1} to stay non-trivial. Lower diff_order or raise "
+                    "output_dim."
+                )
+            penalty = self._difference_penalty(n_basis)
             if self.include_bias:
                 penalty = np.pad(penalty, ((1, 0), (1, 0)))
             self.knots_.append(knots)
@@ -238,13 +250,16 @@ def get_penalty_matrices(self):
         penalties : list of ndarray
             One full penalty matrix per marginal direction, each formed as a
             Kronecker product of a marginal difference penalty with identity
-            matrices for the remaining dimensions.
+            matrices for the remaining dimensions, in dimension order so the
+            result lines up with the ``einsum`` + ``reshape`` flatten order used
+            by :meth:`transform` (dimension 0 slowest/outermost, the last
+            dimension fastest/innermost).
         """
         kron_penalties = []
         for i, Si in enumerate(self.penalties_):
-            mats = [np.eye(size) for j, size in enumerate(self.marginal_sizes_) if j != i]
-            P = Si
-            for M in mats:
+            mats = [Si if j == i else np.eye(size) for j, size in enumerate(self.marginal_sizes_)]
+            P = mats[0]
+            for M in mats[1:]:
                 P = np.kron(P, M)
             kron_penalties.append(P)
         return kron_penalties

pretab/expansion/spline/multivariate/thin_plate.py

@@ -1,3 +1,5 @@
+import warnings
+
 import numpy as np
 from scipy.linalg import eigh
 from scipy.spatial.distance import cdist
@@ -6,7 +8,7 @@
 from sklearn.utils import check_random_state
 from sklearn.utils.validation import check_is_fitted
 
-from ....exceptions import InsufficientSamplesError, InvalidParamError
+from ....exceptions import ConfigWarning, InsufficientSamplesError, InvalidParamError
 from ..mixins import SplineBasisMixin
 
 _LANDMARK_STRATEGIES = ("kmeans", "subsample")
@@ -56,7 +58,8 @@ class ThinPlateSplineTransformer(SplineBasisMixin, TransformerMixin, BaseEstimat
         The retained eigenvalues of the projected landmark kernel.
     penalty_ : ndarray
         Diagonal smoothing penalty of ``eigvals_`` (with an unpenalized leading
-        row/column when ``include_bias=True``).
+        row/column when ``include_bias=True``). **Experimental**: not guaranteed
+        positive semi-definite, see :meth:`get_penalty_matrix`.
     d_ : int
         Number of input features (also ``n_features_in_``).
     n_basis_ : list of int
@@ -202,6 +205,14 @@ def transform(self, X):
     def get_penalty_matrix(self, feature_index=0):
         """Return the smoothing penalty matrix for the fitted basis.
 
+        .. warning::
+           **Experimental.** The retained eigenvalues of the projected landmark
+           kernel are not guaranteed to be non-negative, so this penalty is not
+           guaranteed positive semi-definite. Using it for penalized-regression
+           smoothing can make an otherwise convex problem non-convex. This does
+           not affect :meth:`transform`, only this penalty. A warning is emitted
+           on every call until this is fully fixed.
+
         Parameters
         ----------
         feature_index : int, default=0
@@ -215,4 +226,11 @@ def get_penalty_matrix(self, feature_index=0):
             leading row/column when ``include_bias=True``).
         """
         check_is_fitted(self, "penalty_")
+        warnings.warn(
+            "ThinPlateSplineTransformer.get_penalty_matrix() is experimental: the "
+            "retained eigenvalues are not guaranteed non-negative, so the returned "
+            "penalty is not guaranteed positive semi-definite.",
+            ConfigWarning,
+            stacklevel=2,
+        )
         return self.penalty_

pretab/expansion/spline/natural_cubic.py

@@ -232,8 +232,10 @@ def get_penalty_matrix(self, feature_index=0):
         """
         check_is_fitted(self, "knots_")
         knots = self.knots_[feature_index]
-        B = self._basis(np.linspace(knots[0], knots[-1], 200), knots)
-        B_dd = np.gradient(np.gradient(B, axis=0), axis=0)
+        x_grid = np.linspace(knots[0], knots[-1], 200)
+        B = self._basis(x_grid, knots)
+        B_d = np.gradient(B, x_grid, axis=0)
+        B_dd = np.gradient(B_d, x_grid, axis=0)
 
         n_basis = B.shape[1]
         P = np.zeros((n_basis, n_basis))
@@ -242,6 +244,6 @@ def get_penalty_matrix(self, feature_index=0):
         for i in range(offset, n_basis):
             for j in range(offset, n_basis):
                 integrand = B_dd[:, i] * B_dd[:, j]
-                P[i, j] = np.trapezoid(integrand, np.linspace(knots[0], knots[-1], 200))
+                P[i, j] = np.trapezoid(integrand, x_grid)
 
         return P

pretab/expansion/spline/p_spline.py

@@ -40,6 +40,8 @@ class PSplineTransformer(SplineBasisMixin, TransformerMixin, BaseEstimator):
     diff_order : int, default=2
         The order of the difference penalty used to compute the smoothness penalty matrix.
         For example, 2 corresponds to a second-order difference penalty (encouraging smooth second derivatives).
+        Must be a positive integer and small enough that the fitted basis width supports it
+        (``diff_order <= output_dim - 1``, roughly); both are validated at ``fit``.
 
     include_bias : bool, default=False
         If True, prepend a constant intercept column per feature. The bias term is
@@ -148,6 +150,8 @@ def fit(self, X, y=None):
             raise InvalidParamError(
                 f"output_dim must be >= degree + 1 = {self.degree + 1} for the p-spline basis, got {output_dim}"
             )
+        if not isinstance(self.diff_order, (int, np.integer)) or self.diff_order < 1:
+            raise InvalidParamError(f"diff_order must be a positive integer (>= 1); got {self.diff_order!r}.")
 
         self.knots_ = []
         self.penalty_ = []
@@ -164,6 +168,13 @@ def fit(self, X, y=None):
                 x, y, output_dim, self.degree, strategy, None, None, min_interior, max_interior
             )
             n_basis = len(knots) - self.degree - 1
+            if self.diff_order > n_basis - 1:
+                raise InvalidParamError(
+                    f"diff_order={self.diff_order} is too large for a basis of width {n_basis} "
+                    f"(output_dim={output_dim}); the finite-difference penalty needs "
+                    f"diff_order <= {n_basis - 1} to stay non-trivial. Lower diff_order or "
+                    "raise output_dim."
+                )
             D = np.eye(n_basis)
             for _ in range(self.diff_order):
                 D = np.diff(D, n=1, axis=0)

tests/expansion/spline/test_cubic_transformer.py

@@ -49,6 +49,31 @@ def test_cubic_spline_penalty_matrix_shape():
     assert np.allclose(P, P.T, atol=1e-6)
 
 
+def test_cubic_spline_penalty_matrix_last_knot_value():
+    # Hand-derived closed form: 36 * integral of (x - knot) ** 2 dx from knot to
+    # x_max. With interior knots at [2.5, 5, 7.5] and x_max = 10, the last-knot
+    # entry is 187.5.
+    X = np.linspace(0, 10, 200).reshape(-1, 1)
+    transformer = CubicRegressionSplineTransformer(output_dim=6, placement_strategy="uniform")
+    transformer.fit(X)
+
+    assert transformer.knots_[0] == pytest.approx([2.5, 5.0, 7.5])
+    P = transformer.get_penalty_matrix()
+    assert P[-1, -1] == pytest.approx(187.5, rel=1e-6)
+
+
+def test_cubic_spline_penalty_matrix_penalizes_polynomial_columns():
+    # x**2 and x**3 have real curvature (f'' = 2 and f'' = 6x) and must not be
+    # silently unpenalized; only the linear x column has f'' = 0 everywhere.
+    X = np.linspace(0, 10, 200).reshape(-1, 1)
+    transformer = CubicRegressionSplineTransformer(output_dim=6).fit(X)
+    P = transformer.get_penalty_matrix()
+
+    assert not np.allclose(P[1, :], 0.0)  # x**2 row
+    assert not np.allclose(P[2, :], 0.0)  # x**3 row
+    assert np.allclose(P[0, :], 0.0)  # x row: identically zero second derivative
+
+
 def test_cubic_feature_names_out():
     X = np.random.rand(20, 2)
     transformer = CubicRegressionSplineTransformer(output_dim=8)

tests/expansion/spline/test_naturalcubic_transformer.py

@@ -48,6 +48,22 @@ def test_natural_spline_penalty_matrix_symmetry():
     assert np.allclose(P, P.T, atol=1e-6)
 
 
+def test_natural_spline_penalty_matrix_is_grid_density_invariant():
+    # Regression test: the old np.gradient(..., axis=0) call (missing the x_grid
+    # spacing argument) differentiated with respect to grid *index*, so the
+    # returned penalty shrunk by roughly dx**4 whenever the evaluation range
+    # changed. The penalty is computed on a fixed internal 200-point grid
+    # regardless of range, so as a proxy we assert it is not vanishingly small
+    # relative to a hand-checkable order of magnitude on a modest range.
+    X = np.linspace(0, 10, 200).reshape(-1, 1)
+    transformer = NaturalCubicSplineTransformer(output_dim=5, placement_strategy="uniform").fit(X)
+    P = transformer.get_penalty_matrix()
+
+    # The buggy implementation produced a diagonal on the order of 1e-3 for this
+    # range; the correctly-scaled penalty is on the order of 1e2-1e3.
+    assert np.diag(P).max() > 10.0
+
+
 def test_natural_spline_feature_names_out():
     X = np.random.rand(20, 2)
     transformer = NaturalCubicSplineTransformer(output_dim=5)

tests/expansion/spline/test_pspline_transformer.py

@@ -4,6 +4,7 @@
 import pytest
 from sklearn.exceptions import NotFittedError
 
+from pretab.exceptions import InvalidParamError
 from pretab.transformers import PSplineTransformer
 
 
@@ -50,6 +51,19 @@ def test_pspline_penalty_matrix_shape_and_symmetry():
     assert np.allclose(P, P.T, atol=1e-6)
 
 
+@pytest.mark.parametrize("diff_order", [-1, 0])
+def test_pspline_rejects_nonpositive_diff_order(diff_order):
+    X = np.linspace(0, 1, 30).reshape(-1, 1)
+    with pytest.raises(InvalidParamError, match="diff_order must be a positive integer"):
+        PSplineTransformer(output_dim=8, diff_order=diff_order).fit(X)
+
+
+def test_pspline_rejects_diff_order_too_large_for_output_dim():
+    X = np.linspace(0, 1, 30).reshape(-1, 1)
+    with pytest.raises(InvalidParamError, match="diff_order=50 is too large"):
+        PSplineTransformer(output_dim=8, diff_order=50).fit(X)
+
+
 def test_pspline_feature_names_out():
     X = np.random.rand(20, 2)
     transformer = PSplineTransformer(output_dim=5)

tests/expansion/spline/test_spline_api_parity.py

@@ -7,7 +7,7 @@
 import numpy as np
 import pytest
 
-from pretab.exceptions import IncompatibleParamsError
+from pretab.exceptions import ConfigWarning, IncompatibleParamsError
 from pretab.transformers import (
     CubicRegressionSplineTransformer,
     NaturalCubicSplineTransformer,
@@ -142,6 +142,7 @@ def test_tensor_penalty_matrix_signature_parity():
 def test_thinplate_penalty_matrix_accepts_feature_index():
     X = np.linspace(0, 1, 40).reshape(-1, 1)
     transformer = ThinPlateSplineTransformer(n_components=6, include_bias=True, random_state=0).fit(X)
-    P = transformer.get_penalty_matrix(feature_index=0)
+    with pytest.warns(ConfigWarning, match="experimental"):
+        P = transformer.get_penalty_matrix(feature_index=0)
     assert P.shape == (7, 7)
     assert np.allclose(P[0, :], 0.0) and np.allclose(P[:, 0], 0.0)