# Second differences notion behind GAM penalties

I'm going back through Simon Wood's book on generalized additive models (GAMs) and came back across the definition of the penalty term employed, which is supposed to combat overfitting of smooths. The equation used is the following:

$$\Vert{y-X\beta}\Vert^2 + \lambda \sum_{j=2}^{k-1}\{ f(x^*_{j-1})-2f(x^*_j) + f(x^*_{j+1}) \}^2$$

The book goes on to explain this formula in the following sentence:

The summation term measures wiggliness as a sum of squared second differences of the function at the knots (which crudely approximates the integrated squared second derivative penalty used in cubic spline smoothing: see section 5.1.2, p. 198). When f is very wiggly the penalty will take high values and when f is ‘smooth’ the penalty will be low. If f is a straight line then the penalty is actually zero.

So here is what I get from this paragraph:

• The summation from the right side of the equation, if higher, indicates high wiggliness and will be penalized thereafter.
• The second differences are some kind of approximation of quadratic, cubic, etc. relations with the variable.
• There is no penalty assigned when the function is linear because there is no wiggliness to add to the right side of this equation.

My question is the following...how exactly does this equation achieve this? I'm not sure if the $$\lambda$$ component is part of the confusion, because from what I see it is simply labeled as the smoothing parameter on the next page. The three functions within the equation confuse me as to what they are doing exactly, and the exponent on the outside of the brackets also doesn't make it immediately clear what is going on.

The three bits of this equation

$$\Vert{y-X\beta}\Vert^2 + \lambda \sum_{j=2}^{k-1}\{ f(x^*_{j-1})-2f(x^*_j) + f(x^*_{j+1}) \}^2$$

are:

1. $$\Vert{y-X\beta}\Vert^2$$ is the sum of squared errors, which is measuring the lack of fit of the piecewise linear spline being fitted

2. $$\lambda$$ is the smoothing parameter, which we use to control how much penalty we pay for the wiggliness of the estimated spline

3. $$\sum_{j=2}^{k-1}\{ f(x^*_{j-1})-2f(x^*_j) + f(x^*_{j+1}) \}^2$$ is the sum of squared second differences

The last term is the same as a second order central finite difference of the values of the spline $$f()$$ at the knot locations. The second order central finite difference is

$$f^{\prime \prime}(x) \approx \frac{f(x + h) - 2f(x) + f(x - h)}{h^2}$$

which as shown is an approximation of the second order derivative of the function $$f()$$. As we are differencing evenly spaced knots, $$h = 1$$ and so the $$h$$s disappear yielding essentially the same thing as in the equation you quoted.

What we are doing here is getting an estimate of the "curvature" (second derivative) of the spline at each internal knot by differencing the values of the spline at the focal knot plus the previous and next knots.

If the spline were flat, the values of the function at the focal knot and the previous and next knots would be the same value so the bit in {} for a single focal knot would be equal to 0. If the spline was a straight line, then $$f(x^*_{j+1})$$ would be twice as large as $$f(x^*_{j-1})$$, but as we subtracted off $$2f(x^*_j)$$ so the sum also equals 0. The constant function or other straight line have 0 curvature and this is what we see when we think about those functions in terms of the second order central finite difference.

We square this second order difference because the spline could be curved down or up (concave or convex) and as such a spline that oscillated about some value would likely have zero second order finite difference when summed over the spline because of all the peaks and valleys cancelling each other out. But this would incorrectly measure the wiggliness of such a function; intuition tells us that an oscillating function (say a sine wave) has a lot of wiggliness. Squaring the finite difference estimate of the second derivative solves this problem. It is the same idea as summing squared errors rather than summing errors as a measure of lack of fit.

As for $$\lambda$$, this controls how much penalty we pay for wiggliness in the spline. If we let $$\lambda \rightarrow \infty$$, then any wiggliness in the spline will dominate the equation, so the smallest value of the entire equation will be when the penalty (the summation term) is 0, i.e. when we have a straight line. Then we recover a linear model (in this case). If we let $$\lambda \rightarrow 0$$ then we pay no penalty for the wiggliness and we recover an unpenalised piecewise linear regression fit. This latter fit model would likely be overfitted to the data, while the former would be underfitted if the true relationship was not linear. Hence there will be some optimal value of $$\lambda$$ that allows the spline/model to fit the data reasonably well without overfitting, and we can recover a non-linear estimate of the relationship between $$y$$ and $$x$$ if one exists (given sufficient data for the noise level).

You can think of $$\lambda$$ as controlling the trade off between the smoothness of the estimated function (the spline) and the closeness of the fit to the data.

If you flip the page to p. 168 Simon goes on to explain that the summation can be written more concisely if we realize that for the tent basis functions he is using in the example, the value of the function at the $$j$$th knot ($$f(x^*_j)$$) is the same as the coefficients of $$f$$, i.e. $$\beta_j = f(x^*_j)$$. Substituting these coefficients into the term in the sum we would get $${\beta_{j-1} - 2\beta_j + \beta_{j+1}}$$. If we write out all the terms in the sum (for $$j=2$$ through $$j=k-1$$), and stack them in a matrix, we get

$$\begin{bmatrix} \beta_1 - 2\beta_2 + \beta_3\\ \beta_2 - 2\beta_3 + \beta_4\\ \beta_3 - 2\beta_4 + \beta_5\\ \vdots \\ \end{bmatrix}$$

and we can decompose this matrix into a a vector of coefficients $$\boldsymbol{\beta}$$ and a difference matrix $$\mathbf{D}$$

$$\begin{bmatrix} \beta_1 - 2\beta_2 + \beta_3 \\ \beta_2 - 2\beta_3 + \beta_4 \\ \beta_3 - 2\beta_4 + \beta_5 \\ \vdots \end{bmatrix} = \begin{bmatrix} 1 & -2 & 1 & 0 & \cdots & \cdots & \vdots \\ 0 & 1 & -2 & 1 & 0 & \cdots & \vdots \\ 0 & 0 & 1 & -2 & 1 & 0 & \vdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \end{bmatrix} \begin{bmatrix} \beta_1 \\ \beta_2 \\ \beta_3 \\ \vdots \end{bmatrix} = \mathbf{D}\boldsymbol{\beta}$$

To make progress with other kinds of splines (and the more compact notation used by Simon in the rest of his book), we want to write the penalty in quadratic form, $$\boldsymbol{\beta}^{\mathsf{T}}\mathbf{S}\boldsymbol{\beta}$$ for a penalty matrix $$\mathbf{S}$$. We can turn $$\mathbf{D}$$ into $$\mathbf{S}$$ via

$$\mathbf{S} = \mathbf{D}^{\mathsf{T}}\mathbf{D}$$

and hence we can equate each of the main topics we have discussed into

$$\sum_{j=2}^{k-1}\{ f(x^*_{j-1})-2f(x^*_j) + f(x^*_{j+1}) \}^2 = \boldsymbol{\beta}^{\mathsf{T}}\mathbf{D}^{\mathsf{T}}\mathbf{D}\boldsymbol{\beta} = \boldsymbol{\beta}^{\mathsf{T}}\mathbf{S} \boldsymbol{\beta}$$

• Wow this is a wonderful answer. Thank you for spending the time writing out a comprehensive explanation. Mar 14 at 23:33