1
$\begingroup$

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: enter image description here

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

$\endgroup$
0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ So how much is the integral weighted? How large is the factor c, that I observe? What other parameters influence c? $\endgroup$ – Jakob Oct 5 '19 at 20:54
  • $\begingroup$ Can you present some of your experiments and the reasoning leading to your $c$? $\endgroup$ – kjetil b halvorsen Oct 5 '19 at 21:16
  • $\begingroup$ Yes, I added the experiment to my question (@kjetil b halvorsen). $\endgroup$ – Jakob Oct 5 '19 at 21:57
0
$\begingroup$

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: enter image description here

also for example l=2 works perfectly: enter image description here

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.