# What are the exact factors used in smooth.spline by R?

What exactly is optimized mathematically when I use:

smooth.spline(x, y, lambda)


in terms of the integrated second derivative?

Is it $$\min_{f\in C^2} \sum_{i=1}^N (f(x_i)-y_i)^2 + \lambda\int_{-\infty}^{\infty}\left(f''(x)\right)^2dx$$ or is it $$\min_{f\in C^2}\frac{1}{{\color{"red"}N}}\sum_{i=1}^N (f(x_i)-y_i)^2 + \lambda\int_{-\infty}^{\infty}\left(f''(x)\right)^2dx$$ or is it something different?

Could it be that it is: $$\min_{f\in C^2} \sum_{i=1}^N (f(x_i)-y_i)^2 + \lambda {\color{"red"}c}\int_{-\infty}^{\infty}\left(f''(x)\right)^2dx,$$ where $$c\approx 8$$ (this would explain some of my experiments, but my experiments are quite complicated, so there could be a mistake in my experiment.)?

smooth.spline is a function in the programming language R to calculate natural cubic smoothing splines.

PS: I made more experiments that look like that there is a strange factor, but it might not be completely constant.

PPS(corrected mistakes): My experiment: I ran the R-code:

x=c(-1,-2**(-21),2**(-21),1)
y=c(1,0,0,1)

spline=smooth.spline(x,y,lambda = 0.1)
splineCorrected=smooth.spline(x,y,lambda = 0.1/8)
a=0.1875
b=1.875

x_smooth=seq(from = -1, to = 1, length.out = 1920)
plot(x,y)
lines(x_smooth,a+b*x_smooth^2/2-b*abs(x_smooth)^3/6,col='green',lwd=21)
lines(x_smooth,predict(spline,x_smooth)$$y,col='red',lwd=7) lines(x_smooth,predict(splineCorrected,x_smooth)$$y,col='blue',lwd=7)


which outputs:

and I calculated the solution for this simple 3 point example (which actually contains 4 points because the middle point exists kind of 2 times) by hand, which I plotted in green: Out of symetry-reasons the spline function $$f$$ has to have the form: $$f(x)=a+\frac{bt^2}{2}-\frac{b|t|^3}{6}$$ This results in a total loss function of: $$\sum_{i=1}^N (f(x_i)-y_i)^2 + \lambda\int_{-\infty}^{\infty}\left(f''(x)\right)^2dx=\\ 2\left(1-(a+\frac{b}{2}-\frac{b}{6})\right)^2 +2a^2 + \lambda 2 \left(b^2-b^2+\frac{b^2}{3}\right)$$

This optimization problem can be solved by Wolfram alpha: https://www.wolframalpha.com/input/?i=minimize+2*%281-%28a%2Bb%2F2-b%2F6%29%29%C2%B2%2B2a%C2%B2%2B0.1*2*%28b%C2%B2-b%C2%B2%2Bb%C2%B2%2F3%29+

https://www.wolframalpha.com/input/?i=minimize+2*%281-%28a%2Bb%2F2-b%2F6%29%29%C2%B2%2B2%28a%2Bb*2%5E%28-2*21%29%2F2-b*2%5E%28-3*21%29%2F6%29%C2%B2%2B0.1*2*%28b%C2%B2-b%C2%B2%2Bb%C2%B2%2F3%29+

The result of this optimization is used to plot the green line.

Conclusion: The spline calculated by R (red line) looks very different from the analytic result (green line). But when I change $$\lambda$$ by a factor of $$c=8$$, I get the blue line, which is (almost) identical to the analytic solution.

PPPS: Does somebody have an idea for a simpler example where one can calculate the analytic solution by hand and compare it to the result obtained by R? smooth.spline needs at least 4 data points to run (that's why I changed my example from 3 to 4 samples by taking two slightly shifted versions $$(-2^{-21},+2^{-21})$$ of the point $$(x_2,y_2)=(0,0)$$ in the middel).

short answer: $$\min_{f\in C^2} \sum_{i=1}^N (f(x_i)-y_i)^2 + \lambda {\color{"red"}c}\int_{-\infty}^{\infty}\left(f''(x)\right)^2dx,$$ where $$c=\left(\max_{i\in\{1,...,N\}}(x_i)-\min_{i\in\{1,...,N\}}(x_i)\right)^3$$ is calculated by smooth.spline(x, y, lambda).

I think the reason for this is that R internally scales the data to an interval of length one before it fits the spline and afterwards it scales it back. This result in the factor $$c=s^3$$, where s is the scaling factor $$s=\max x-\min x$$.

PS: Now, all my experiments are consistent with this theory.

l=2
lambda=0.3
x=c(-l,-2**(-20),2**(-20),l)
y=c(1,0,0,1)

spline=smooth.spline(x,y,lambda = lambda)
splineCorrected=smooth.spline(x,y,lambda = lambda/((2*l)**3))
a=3*lambda/(l^3+6*lambda)
b=3*l/(l^3+6*lambda)

x_smooth=seq(from = -l, to = l, length.out = 1920)
plot(x,y)
lines(x_smooth,a+b*x_smooth^2/2-b*abs(x_smooth)^3/(6*l),col='green',lwd=21)
lines(x_smooth,predict(spline,x_smooth)$$y,col='red',lwd=7) lines(x_smooth,predict(splineCorrected,x_smooth)$$y,col='blue',lwd=7)

smooth.spline(x = x, y = y, lambda=0.2)


works fine:

also for example l=2 works perfectly:

https://www.wolframalpha.com/input/?i=minimize+2*%281-%28a%2Bl%5E2*b%2F2-l%5E3*b%2F%286*l%29%29%29%C2%B2%2B2a%C2%B2%2Be*2*%28l*b%C2%B2-l%C2%B2b%C2%B2%2Fl%2Bl%C2%B3b%C2%B2%2F%283*l%C2%B2%29%29%2C+where+l%3Dpi

Did you read the output of ?smooth.spline? Everything is in there, and the calculated model object contains all this information. A very simple example:

testfun <- function(x) x+0.5*x^2 + 0.1 * x^3 + 0.001*x^4
x <- (-5):5
set.seed(7*11*13) # My public seed
Y <- testfun(x)+rnorm(x, sd=0.1)
test <- smooth.spline(x, Y, df=4)
test

Call:
smooth.spline(x = x, y = Y, df = 4)

Smoothing Parameter  spar= 0.4987503  lambda= 0.00262488 (14 iterations)
Equivalent Degrees of Freedom (Df): 3.999435

• Can you present some of your experiments and the reasoning leading to your $c$? Oct 5, 2019 at 21:16