# Find good smoothing spline factor

I'd like an automatic way to find the "best" smoothing factor s for a spline fit to a given set of data points. Here's a sample visualization of some data and the fit splines for various s values:

In this case, clearly s=2 (and at a lesser degree s=1) is not a good fit. On the other hand s=0.5 fits the data almost as good as s=0.1 but with less than half the number of knots and thus is less susceptible to overfitting. So my question is, what's a robust method to determine the "optimal", or at least a good enough s to fit the data?

• You might consider cross-validation approaches. – Glen_b Apr 12 '18 at 5:59
• Are you working in any particular programming language? Also, what type of spline are you working with? – Jim Apr 12 '18 at 6:16
• I'm using scipy's UnivariateSpline – gsakkis Apr 13 '18 at 6:49
• If you iterate over a range of smoothing factors, each time calculating R-squared (R2) as "R2 = 1.0 - (absolute_error_variance / dependent_data_variance)" (this would run quickly) and then plot smoothing factor vs. R-squared, you should be able to see where the approximate optimal smoothing factor is on that plot. This is easy to program and might be worth a try as it is simple to do. For a numpy array named X, the variance is X.var() so the R-squared calculation becomes trivial with numpy. – James Phillips Apr 14 '18 at 15:06
• Wow, the content in the scipy doc is incorrect. This is certainly not a smoothing spline. Could you post the data? I would like to provide at least a clean R solution. – Jim Apr 14 '18 at 19:04

What is a smoothing spline?

The Wikipedia article on smoothing splines does a good job in explaining that. To recap, given a set of data points, $\{ (x_i, y_i)_{i=1}^n \}$, a smoothing spline is a solution to the interpolation problem:

$$\underset{f}{\arg\min} \sum_{i=1}^n (y_i - f(x_i))^2 + \lambda \int_{x_{(1)}}^{x_{(n)}} f''(x)^2 dx,$$

with $f$ constrained to be piecewise cubic between different $x_i$. The first part measures the goodness of fit of such an $f$ to the observed data. The second part is a penalty term for the wiggliness (non-smoothness) of $f$.

Leaving it to us to find a good trade-off between fit and smoothness by means of $\lambda$.

Smoothing splines in R

Luckily R has the splines package that does the heavy lifting for us.

library(splines)

myspline <- smooth.spline(x = mydata$x, y = mydata$y
, lambda = 8e-9 # optim 8.332658e-11
, cv = TRUE)

xgrid <- sort(union(mydata$x , seq(from = min(mydata$x), to = max(mydata$x), by = 1)) , decreasing = FALSE) yhat_xgrid <- predict(myspline, x = xgrid)$y

plot(x = mydata$x, y = mydata$y, log = "x", ylim = c(0,1)
, xlab = "x (log-scale)", ylab = "y"
, col  = "lightblue", pch = 19)
lines(x = xgrid, y = yhat_xgrid, type = "l", col = "darkorange")
grid()
legend(...)


And we obtain this lovely plot.

The optimal values for $\lambda$ are $\hat{\lambda}^*_{\text{LOO}} = 8.33 \cdot 10^{-11}$ and $\hat{\lambda}^*_{\text{GCV}} = 5.81 \cdot 10^{-13}$. I like the one plotted: $\hat{\lambda}^*_{\text{Jim}} = 8 \cdot 10^{-9}$.

• In stackoverflow.com/q/55363639/10698244 this post is referred to and it is asked for a smoothing spline implementation in Python. As it seems, there are different thoughts whether UnivariateSpline is a smoothing spline or not. Perhaps you could participate in the discussion there? – user7468395 Mar 28 '19 at 20:17