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

Question

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).

Answer 1

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

Answer 2

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
Penalized Criterion (RSS): 16.03144
GCV: 3.598315

should answer your question.