\(\newcommand{\ci}{\perp\!\!\!\perp}\) \(\newcommand{\cA}{\mathcal{A}}\) \(\newcommand{\cB}{\mathcal{B}}\) \(\newcommand{\cC}{\mathcal{C}}\) \(\newcommand{\cD}{\mathcal{D}}\) \(\newcommand{\cE}{\mathcal{E}}\) \(\newcommand{\cF}{\mathcal{F}}\) \(\newcommand{\cG}{\mathcal{G}}\) \(\newcommand{\cH}{\mathcal{H}}\) \(\newcommand{\cI}{\mathcal{I}}\) \(\newcommand{\cJ}{\mathcal{J}}\) \(\newcommand{\cK}{\mathcal{K}}\) \(\newcommand{\cL}{\mathcal{L}}\) \(\newcommand{\cM}{\mathcal{M}}\) \(\newcommand{\cN}{\mathcal{N}}\) \(\newcommand{\cO}{\mathcal{O}}\) \(\newcommand{\cP}{\mathcal{P}}\) \(\newcommand{\cQ}{\mathcal{Q}}\) \(\newcommand{\cR}{\mathcal{R}}\) \(\newcommand{\cS}{\mathcal{S}}\) \(\newcommand{\cT}{\mathcal{T}}\) \(\newcommand{\cU}{\mathcal{U}}\) \(\newcommand{\cV}{\mathcal{V}}\) \(\newcommand{\cW}{\mathcal{W}}\) \(\newcommand{\cX}{\mathcal{X}}\) \(\newcommand{\cY}{\mathcal{Y}}\) \(\newcommand{\cZ}{\mathcal{Z}}\) \(\newcommand{\bA}{\mathbf{A}}\) \(\newcommand{\bB}{\mathbf{B}}\) \(\newcommand{\bC}{\mathbf{C}}\) \(\newcommand{\bD}{\mathbf{D}}\) \(\newcommand{\bE}{\mathbf{E}}\) \(\newcommand{\bF}{\mathbf{F}}\) \(\newcommand{\bG}{\mathbf{G}}\) \(\newcommand{\bH}{\mathbf{H}}\) \(\newcommand{\bI}{\mathbf{I}}\) \(\newcommand{\bJ}{\mathbf{J}}\) \(\newcommand{\bK}{\mathbf{K}}\) \(\newcommand{\bL}{\mathbf{L}}\) \(\newcommand{\bM}{\mathbf{M}}\) \(\newcommand{\bN}{\mathbf{N}}\) \(\newcommand{\bO}{\mathbf{O}}\) \(\newcommand{\bP}{\mathbf{P}}\) \(\newcommand{\bQ}{\mathbf{Q}}\) \(\newcommand{\bR}{\mathbf{R}}\) \(\newcommand{\bS}{\mathbf{S}}\) \(\newcommand{\bT}{\mathbf{T}}\) \(\newcommand{\bU}{\mathbf{U}}\) \(\newcommand{\bV}{\mathbf{V}}\) \(\newcommand{\bW}{\mathbf{W}}\) \(\newcommand{\bX}{\mathbf{X}}\) \(\newcommand{\bY}{\mathbf{Y}}\) \(\newcommand{\bZ}{\mathbf{Z}}\) \(\newcommand{\ba}{\mathbf{a}}\) \(\newcommand{\bb}{\mathbf{b}}\) \(\newcommand{\bc}{\mathbf{c}}\) \(\newcommand{\bd}{\mathbf{d}}\) \(\newcommand{\be}{\mathbf{e}}\) \(\newcommand{\bg}{\mathbf{g}}\) \(\newcommand{\bh}{\mathbf{h}}\) \(\newcommand{\bi}{\mathbf{i}}\) \(\newcommand{\bj}{\mathbf{j}}\) \(\newcommand{\bk}{\mathbf{k}}\) \(\newcommand{\bl}{\mathbf{l}}\) \(\newcommand{\bm}{\mathbf{m}}\) \(\newcommand{\bn}{\mathbf{n}}\) \(\newcommand{\bo}{\mathbf{o}}\) \(\newcommand{\bp}{\mathbf{p}}\) \(\newcommand{\bq}{\mathbf{q}}\) \(\newcommand{\br}{\mathbf{r}}\) \(\newcommand{\bs}{\mathbf{s}}\) \(\newcommand{\bt}{\mathbf{t}}\) \(\newcommand{\bu}{\mathbf{u}}\) \(\newcommand{\bv}{\mathbf{v}}\) \(\newcommand{\bw}{\mathbf{w}}\) \(\newcommand{\bx}{\mathbf{x}}\) \(\newcommand{\by}{\mathbf{y}}\) \(\newcommand{\bz}{\mathbf{z}}\) \(\newcommand{\RR}{\mathbb{R}}\) \(\newcommand{\NN}{\mathbb{N}}\) \(\newcommand{\balpha}{\boldsymbol{\alpha}}\) \(\newcommand{\bbeta}{\boldsymbol{\beta}}\) \(\newcommand{\btheta}{\boldsymbol{\theta}}\) \(\newcommand{\hpi}{\widehat{\pi}}\) \(\newcommand{\bpi}{\boldsymbol{\pi}}\) \(\newcommand{\hbpi}{\widehat{\boldsymbol{\pi}}}\) \(\newcommand{\bxi}{\boldsymbol{\xi}}\) \(\newcommand{\bmu}{\boldsymbol{\mu}}\) \(\newcommand{\bepsilon}{\boldsymbol{\epsilon}}\) \(\newcommand{\bzero}{\mathbf{0}}\) \(\newcommand{\T}{\text{T}}\) \(\newcommand{\Trace}{\text{Trace}}\) \(\newcommand{\Cov}{\text{Cov}}\) \(\newcommand{\Var}{\text{Var}}\) \(\newcommand{\E}{\mathbb{E}}\) \(\newcommand{\Pr}{\text{Pr}}\) \(\newcommand{\pr}{\text{pr}}\) \(\newcommand{\pdf}{\text{pdf}}\) \(\newcommand{\P}{\text{P}}\) \(\newcommand{\p}{\text{p}}\) \(\newcommand{\One}{\mathbf{1}}\) \(\newcommand{\argmin}{\operatorname*{arg\,min}}\) \(\newcommand{\argmax}{\operatorname*{arg\,max}}\) \(\newcommand{\dtheta}{\frac{\partial}{\partial\theta} }\) \(\newcommand{\ptheta}{\nabla_\theta}\) \(\newcommand{\alert}[1]{\color{darkorange}{#1}}\) \(\newcommand{\alertr}[1]{\color{red}{#1}}\) \(\newcommand{\alertb}[1]{\color{blue}{#1}}\)

1 Function Approximation

1.1 Approximating a Function

In statistical modeling, we often observe data generated from

\[ y_i = f_0(\bx_i) + \varepsilon_i, \qquad \mathbb{E}(\varepsilon_i \mid \bx_i) = 0, \qquad i = 1, \ldots, n. \]

Here \(f_0(\bx) = \mathbb{E}(Y \mid \bX = \bx)\) is the unknown regression function, which is our target. Our goal is to use the observed data to construct an estimator \(\widehat f\) that closely resembles \(f_0\).

To be concrete, suppose we choose a class of candidate functions \(\cF\). A natural estimator is

\[ \widehat f = \underset{f \in \cF}{\arg\min} \frac{1}{n}\sum_{i=1}^n \big(y_i-f(\bx_i)\big)^2. \]

There are really two questions here. First, can the functions in \(\cF\) approximate \(f_0\)? Second, can we estimate a good function from a finite and noisy sample? These correspond to approximation error and estimation error. A small function class may not be flexible enough, while a very large function class may lead to unstable estimation or overfitting.

For example, if the true regression function is linear,

\[ f_0(\bx) = \bbeta_0^\top \bx, \]

where \(\bx\) includes an intercept, we have learned to use the ordinary least-squares estimator

\[ \widehat{\bbeta} = (\bX^\top\bX)^{-1}\bX^\top\by \]

to estimate \(\bbeta_0\). Under the usual Gauss-Markov conditions, including a full-rank design, \(\mathbb{E}(\bepsilon\mid\bX)=\bzero\), and \(\operatorname{Var}(\bepsilon\mid\bX)=\sigma^2\bI\), this estimator is unbiased and is the best linear unbiased estimator.

But oftentimes, we may be dealing with a more complex \(f_0\) that is not linear. There are several strategies that we can use:

  • Using higher-order polynomials or nonlinear transformations of the input variables, such as \(x^2\), \(\log(x)\), \(\sin(x)\), and \(I(x>3)\). Broadly speaking, this is called a basis expansion.
  • Using local averages, such as \(k\)-Nearest Neighbors (KNN) or kernel smoothing. This also includes adaptive weighting methods such as random forests.
  • Using neural networks, which can be viewed as a more complex form of basis expansion with multiple layers and nonlinear transformations.

When we consider what type of method to use, the following aspects are important:

  • Flexibility: How flexible is the method in capturing the true function? For example, a second-order polynomial may not be flexible enough to capture a complex function.
  • High-dimensionality: Can the method handle high-dimensional data? High-order polynomial expansions can introduce too many terms when the number of input variables increases.
  • Computational cost: Are the parameters easy to estimate? How much computation is required to fit the model?
  • Theoretical guarantees: Can we analyze the statistical properties of the estimator and understand what controls its complexity?

In this lecture, we will start with polynomial basis expansions, construct splines from local polynomial pieces, and use regularization to control smoothness. These ideas will lead naturally to the RKHS framework in the next lecture.

1.2 Example: Polynomial Basis

It could be helpful to first understand how a function-approximation method works in practice. We will use this example to demonstrate how to estimate a model and predict at new input values using our own code as much as possible. Since this is a simulation, we know that the true regression function is

\[ f_0(x)=0.5+\sin(x), \]

although the fitting procedure only observes the noisy data.

  # generate one realization from the model
  set.seed(123)

  n = 50
  f_true <- function(x) 0.5 + sin(x)
  x <- seq(-pi, pi, length.out = n)
  y <- f_true(x) + rnorm(length(x), sd = 0.2)

  # create a data frame for the observations
  train_data <- data.frame(x = x, y = y)

  # plot the observed data and the known truth
  plot(train_data$x, train_data$y, pch = 19, col = "blue",
       xlab = "x", ylab = "y", main = "Observed Data and Truth")
  lines(x, f_true(x), col = "darkorange", lwd = 3)
  legend("topleft", c("Truth", "Observed data"),
         col = c("darkorange", "blue"), lty = c(1, NA),
         pch = c(NA, 19), lwd = c(3, NA), bty = "n")

A typical statistical modeling procedure would follow these steps:

  • Specify a model form, such as a polynomial basis expansion of degree \(d\).
  • Specify an empirical loss function, such as the squared-error loss.
  • Optimize the loss function on the observed data to find the parameter estimates.
  • Evaluate the estimated function at new input values.

If we consider a polynomial basis expansion of degree \(d\), a candidate function can be written as

\[ f_{\bbeta}(x) = \sum_{j=0}^d \beta_j x^j, \]

where the \(\beta_j\)’s are coefficients to be estimated. The empirical loss based on the observations \(\{(x_i,y_i)\}_{i=1}^n\) is

\[\begin{align} L(\bbeta) &= \frac{1}{n}\sum_{i=1}^n\left(y_i-f_{\bbeta}(x_i)\right)^2 \\ &= \frac{1}{n}\sum_{i=1}^n \left(y_i-\sum_{j=0}^d\beta_jx_i^j\right)^2. \end{align}\]

Let’s consider \(d=3\). To calculate the loss, we only need to evaluate the polynomial basis functions at the observed \(x_i\)’s. The columns of the resulting design matrix are

  • \(x_i^0\), a column of 1 corresponding to the intercept;
  • \(x_i^1\), the original \(x_i\);
  • \(x_i^2\), the quadratic term;
  • \(x_i^3\), the cubic term.
  # construct the design matrix and estimate the coefficients
  X = cbind(1, train_data$x, train_data$x^2, train_data$x^3)
  beta_hat = solve(t(X) %*% X, t(X) %*% train_data$y)
  cat("Estimated coefficients:\n")
## Estimated coefficients:
  beta_hat
##              [,1]
## [1,]  0.485874303
## [2,]  0.865520516
## [3,]  0.006134783
## [4,] -0.094449751

  # fitted values at the observed inputs
  yfit <- X %*% beta_hat

  plot(train_data$x, train_data$y, pch = 19, col = "blue",
       xlab = "x", ylab = "y", main = "Cubic Polynomial Fit")
  lines(train_data$x, f_true(train_data$x), col = "darkorange", lwd = 3)
  lines(train_data$x, yfit, col = "red", lwd = 2)
  legend("topleft", c("Truth", "Fitted function"),
         col = c("darkorange", "red"), lty = 1,
         lwd = c(3, 2), bty = "n")

Within \([-\pi,\pi]\), a cubic polynomial gives a reasonable approximation. But of course, this does not mean that a low-degree polynomial can capture the same function over a wider interval. Let’s use the range \([-2\pi,2\pi]\) and compare degrees 3 and 5. The grid below is only used to display the fitted functions. It is not a held-out testing sample.

  # generate one realization on a wider input range
  set.seed(123)
  x <- seq(-2*pi, 2*pi, length.out = n)
  y <- f_true(x) + rnorm(length(x), sd = 0.2)
  train_data <- data.frame(x = x, y = y)

  # define a plotting grid
  x_grid = seq(-2*pi, 2*pi, length.out = 200)
  grid_data = data.frame(x = x_grid)

  # fit polynomial regression models
  poly3.fit <- lm(y ~ poly(x, degree = 3), data = train_data)
  y_hat3 <- predict(poly3.fit, newdata = grid_data)

  poly5.fit <- lm(y ~ poly(x, degree = 5), data = train_data)
  y_hat5 <- predict(poly5.fit, newdata = grid_data)

  # compare the fitted functions against the known truth
  par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))

  plot(train_data$x, train_data$y, pch = 19, col = "blue",
       xlab = "x", ylab = "y", main = "Polynomial Degree = 3")
  lines(x_grid, f_true(x_grid), col = "darkorange", lwd = 3)
  lines(x_grid, y_hat3, col = "red", lwd = 2)

  plot(train_data$x, train_data$y, pch = 19, col = "blue",
       xlab = "x", ylab = "y", main = "Polynomial Degree = 5")
  lines(x_grid, f_true(x_grid), col = "darkorange", lwd = 3)
  lines(x_grid, y_hat5, col = "red", lwd = 2)
  legend("topleft", c("Truth", "Fitted function"),
         col = c("darkorange", "red"), lty = 1,
         lwd = c(3, 2), bty = "n")

The polynomial basis is convenient because the resulting model is still linear in the coefficients. However, every polynomial basis function changes over the entire input range. Hence, changing one coefficient changes the fitted curve globally. A natural question is whether we can construct basis functions that make more local changes. Let’s use the U.S. birth rate data to make this limitation more visible.

2 Splines

A spline is still linear in its coefficients, but its basis functions are constructed from local polynomial pieces. This allows the fitted function to change its behavior in one region without forcing the same change over the entire input range.

2.1 A Motivating Example: U.S. Birth Rates

We use the U.S. birth rate data as an example. The data records birth rates from 1917 to 2003. The birth rate data show a clear nonlinear trend.

    load("../../dataset/birthrates.Rda")
    head(birthrates)
    par(mar = c(4,4,1,1))
    plot(birthrates, pch = 19, col = "darkorange")

It might be interesting to fit a linear regression with higher-order polynomial terms. The poly() function constructs an orthogonal polynomial basis spanning the polynomial space up to a specified degree. This is more numerically stable than directly computing terms such as I(Year^2) and I(Year^3), because the numerical values of Year are large.

    par(mfrow=c(1,2))

    par(mar = c(2,3,2,0))
    lmfit <- lm(Birthrate ~ poly(Year, 3), data = birthrates)
    plot(birthrates, pch = 19, col = "darkorange")
    lines(birthrates$Year, lmfit$fitted.values, lty = 1, col = "deepskyblue", lwd = 2)
    title("degree = 3")

    par(mar = c(2,3,2,0))
    lmfit <- lm(Birthrate ~ poly(Year, 5), data = birthrates)
    plot(birthrates, pch = 19, col = "darkorange")
    lines(birthrates$Year, lmfit$fitted.values, lty = 1, col = "deepskyblue", lwd = 2)
    title("degree = 5")

2.1.1 From Global to Local Approximation

These polynomial fits do not seem to perform very well. A natural question is whether we can model the curve locally. We know that \(k\)NN uses nearby observations for each prediction. Here, we will try something slightly different. Let’s divide the year range into non-overlapping intervals of 10 years and estimate a constant mean within each interval. For prediction, we only need to identify the interval containing the new input.

    par(mfrow = c(1, 1), mar = c(4, 4, 2, 1))

    breaks = seq(1917, 2007, by = 10)
    interval = cut(
        birthrates$Year,
        breaks = breaks,
        right = FALSE,
        include.lowest = TRUE,
        labels = FALSE
    )

    # one indicator basis function for each interval
    mybasis = sapply(
        seq_len(length(breaks) - 1),
        function(l) as.numeric(interval == l)
    )

    lmfit <- lm(birthrates$Birthrate ~ . - 1, data = data.frame(mybasis))

    plot(birthrates, pch = 19, col = "darkorange")
    abline(v = breaks[2:(length(breaks) - 1)], lty = 2, col = "gray70")
    lines(birthrates$Year, lmfit$fitted.values, type = "s",
          col = "deepskyblue", lwd = 2)
    title("Histogram Regression")

This method is called histogram regression. Suppose the interval containing a new input \(x\) is \(\phi(x)\). Then the fitted function is

\[ \widehat{f}(x) = \frac{\sum_{i=1}^n y_i I\{x_i \in \phi(x)\}}{\sum_{i=1}^n I\{x_i \in \phi(x)\}} \]

The fitted value is simply the average response within the corresponding interval. This is related to the usual histogram idea because both methods summarize observations within fixed bins. However, the fitted regression function has jumps at the interval boundaries. Hence, we need a more flexible framework.

2.2 Piecewise Polynomials

Instead of fitting constant functions within each interval (between two knots), we may consider fitting a line. Consider a simpler case, where we use 3 knots at 1936, 1960, and 1978, which gives 4 intervals.

    par(mfrow=c(1,2))

    myknots = c(1936, 1960, 1978)
    bounds = c(1917, myknots, 2003)

    # piecewise constant
    mybasis = cbind("x_1" = (birthrates$Year < myknots[1]),
                    "x_2" = (birthrates$Year >= myknots[1])*(birthrates$Year < myknots[2]),
                    "x_3" = (birthrates$Year >= myknots[2])*(birthrates$Year < myknots[3]),
                    "x_4" = (birthrates$Year >= myknots[3]))

    lmfit <- lm(birthrates$Birthrate ~ . -1, data = data.frame(mybasis))
    par(mar = c(2,3,2,0))
    plot(birthrates, pch = 19, col = "darkorange")
    abline(v = myknots, lty = 2)
    title("Piecewise constant")

    for (k in 1:4)
        points(c(bounds[k], bounds[k+1]), rep(lmfit$coefficients[k], 2), type = "l",
               lty = 1, col = "deepskyblue", lwd = 4)

    # piecewise linear
    mybasis = cbind("x_1" = (birthrates$Year < myknots[1]),
                    "x_2" = (birthrates$Year >= myknots[1])*(birthrates$Year < myknots[2]),
                    "x_3" = (birthrates$Year >= myknots[2])*(birthrates$Year < myknots[3]),
                    "x_4" = (birthrates$Year >= myknots[3]),
                    "x_11" = birthrates$Year*(birthrates$Year < myknots[1]),
                    "x_21" = birthrates$Year*(birthrates$Year >= myknots[1])*(birthrates$Year < myknots[2]),
                    "x_31" = birthrates$Year*(birthrates$Year >= myknots[2])*(birthrates$Year < myknots[3]),
                    "x_41" = birthrates$Year*(birthrates$Year >= myknots[3]))

    lmfit <- lm(birthrates$Birthrate ~ .-1, data = data.frame(mybasis))
    par(mar = c(2,3,2,0))
    plot(birthrates, pch = 19, col = "darkorange")
    abline(v = myknots, lty = 2)
    title("Piecewise linear")

    for (k in 1:4)
        points(c(bounds[k], bounds[k+1]),
               lmfit$coefficients[k] + c(bounds[k], bounds[k+1])*lmfit$coefficients[k+4],
               type = "l", lty = 1, col = "deepskyblue", lwd = 4)

2.3 Splines

However, these functions are not continuous. Hence we use a trick to construct continuous basis functions. We center Year at the first observed year only to make the basis functions easier to display. This changes the parameterization, but not the spline space.

    par(mfrow=c(1,1))
    pos <- function(x) x*(x>0)
    year_centered = birthrates$Year - min(birthrates$Year)
    knots_centered = myknots - min(birthrates$Year)

    mybasis = cbind("int" = 1, "x_1" = year_centered,
                    "x_2" = pos(year_centered - knots_centered[1]),
                    "x_3" = pos(year_centered - knots_centered[2]),
                    "x_4" = pos(year_centered - knots_centered[3]))

    par(mar = c(2,2,2,0))
    matplot(birthrates$Year, mybasis[, -1], type = "l", lty = 1,
            yaxt = 'n', ylim = range(mybasis[, -1]), lwd = 2)
    title("Spline Basis Functions")

With this definition, any fitted model will be

  • continuous everywhere;
  • linear between adjacent knots;
  • allowed to have a different slope in each region.

The resulting model is called a linear spline.

    lmfit <- lm(birthrates$Birthrate ~ .-1, data = data.frame(mybasis))
    par(mar = c(2,3,2,0))
    plot(birthrates, pch = 19, col = "darkorange")
    lines(birthrates$Year, lmfit$fitted.values, lty = 1, col = "deepskyblue", lwd = 4)
    abline(v = myknots, lty = 2)
    title("Linear Spline")

Of course, writing this out explicitly is very tedious. Hence, we use the bs() function in the splines package.

    par(mar = c(2,2,2,0))
    lmfit <- lm(Birthrate ~ splines::bs(Year, degree = 1, knots = myknots), data = birthrates)
    plot(birthrates, pch = 19, col = "darkorange")
    lines(birthrates$Year, lmfit$fitted.values, lty = 1, col = "deepskyblue", lwd = 4)
    title("Linear spline with the bs() function")

The next step is to increase the degree to account for more complicated functions. There are a few questions we need to consider here:

  • How many knots should be used?
  • Where should the knots be placed?
  • What should the polynomial degree be within each region?

For example, let’s consider the following setting:

    par(mar = c(2,2,2,0))
    lmfit <- lm(Birthrate ~ splines::bs(Year, degree = 3, knots = myknots), data = birthrates)
    plot(birthrates, pch = 19, col = "darkorange")
    lines(birthrates$Year, lmfit$fitted.values, lty = 1, col = "deepskyblue", lwd = 4)
    title("Cubic spline with 3 knots")

All of these choices affect the fitted function. In particular, the number of knots and the polynomial degree determine the model’s degrees of freedom. For simplicity, we can use the df argument to control the number of basis functions. In the following model, bs() returns five basis columns and lm() adds an intercept, giving six fitted coefficients.

    par(mar = c(2,2,2,0))
    lmfit <- lm(Birthrate ~ splines::bs(Year, degree = 3, df = 5), data = birthrates)
    plot(birthrates, pch = 19, col = "darkorange")
    lines(birthrates$Year, lmfit$fitted.values, lty = 1, col = "deepskyblue", lwd = 4)
    title("Cubic Spline with 6 Parameters")

2.4 Spline Basis

There are different ways to construct a basis for the same spline space. We previously used a truncated-power basis and the basis-spline construction, commonly called a B-spline basis. B-spline basis functions have local support and are often more convenient computationally. Here is a comparison of B-spline bases with different degrees.

    par(mfrow = c(4, 1), mar = c(0, 0, 2, 0))

    for (d in 0:3)
    {
        bs_d = splines2::bSpline(1:100, degree = d, knots = seq(10, 90, 10), intercept = TRUE)
        matplot(1:100, bs_d, type = ifelse(d == 0, "s", "l"), lty = 1, ylab = "spline",
                xaxt = 'n', yaxt = 'n', ylim = c(-0.05, 1.05), lwd = 2)
        title(paste("degree =", d))
    }

2.5 Natural Cubic Spline

Extrapolation is generally dangerous because a fitted function can behave extremely outside the range of the observed data. In linear models fit by bs(), extrapolation beyond the boundary knots triggers a warning.

    par(mfrow = c(1, 1), mar = c(4, 4, 2, 1))
    library(splines)

    fit.bs = lm(Birthrate ~ bs(Year, df = 6), data = birthrates)
    fit.ns = lm(Birthrate ~ ns(Year, df = 6), data = birthrates)

    year_grid = seq(1900, 2020)
    pred.bs = predict(fit.bs, data.frame(Year = year_grid))
    pred.ns = predict(fit.ns, data.frame(Year = year_grid))
    y_range = range(c(birthrates$Birthrate, pred.bs, pred.ns))

    plot(birthrates$Year, birthrates$Birthrate,
         ylim = y_range, xlim = range(year_grid),
         pch = 19, xlab = "Year", ylab = "Birth rate",
         col = "darkorange")
    lines(year_grid, pred.bs, col = "deepskyblue", lwd = 3)
    lines(year_grid, pred.ns, col = "darkgreen", lwd = 3)
    legend("topright", c("Cubic B-spline", "Natural Cubic Spline"),
           col = c("deepskyblue", "darkgreen"), lty = 1,
           lwd = 3, cex = 1.1)
    title("Birth Rate Extrapolation")

Hence, this motivates us to impose an additional constraint near the boundaries. A natural cubic spline on \([a,b]\) satisfies the boundary conditions

\[ g''(a)=g''(b)=0, \]

where \(a\) and \(b\) are the boundary knots. Under the usual natural continuation, the spline is linear beyond these two boundary knots. Thus, the fitted function has no curvature outside the observed range.

    par(mfrow = c(1, 1), mar = c(4, 4, 2, 1))
    ncs = ns(1:100, df = 6, intercept = TRUE)
    matplot(1:100, ncs, type = "l", lty = 1,
            xlab = "x", ylab = "Basis value", lwd = 3)
    title("Natural Cubic Spline")

2.6 Smoothing Spline

Once the knots are chosen, a regression spline is just a linear model. However, its performance can depend strongly on the number and placement of the knots. Can we use many knots and let regularization control the complexity instead?

Let’s start with a seemingly “horrible” idea by placing a knot at every observed input value. This gives us a very flexible natural cubic spline, but it may also lead to overfitting. Following the ridge-regression idea, we control the flexibility by adding a roughness penalty.

Let

\[ a=\min_i x_i, \qquad b=\max_i x_i, \]

and consider the second-order Sobolev space

\[ \cW_2^2[a,b] = \left\{ g: g' \text{ is absolutely continuous and } \int_a^b [g''(x)]^2\,dx < \infty \right\}. \]

For a tuning parameter \(\lambda>0\), we estimate the regression function by solving

\[\begin{equation} \widehat g_\lambda = \underset{g\in\cW_2^2[a,b]}{\arg\min} \left\{ \frac{1}{n}\sum_{i=1}^n\big(y_i-g(x_i)\big)^2 + \lambda\int_a^b[g''(x)]^2\,dx \right\}. \tag{2.1} \end{equation}\]

The first term measures how well \(g\) fits the observations. The second term measures curvature. A small \(\lambda\) allows more curvature, while a large \(\lambda\) pushes the fitted function toward a straight line.

2.6.1 Optional Theory: Why Is the Solution a Natural Cubic Spline?

The optimization in Equation (2.1) is over an infinite-dimensional space of functions. A natural question is whether we can reduce it to a finite-dimensional problem that only depends on the observed data.

For simplicity, suppose the input values are distinct and ordered as

\[ a=x_1<x_2<\cdots<x_n=b. \]

Take any candidate function \(g\in\cW_2^2[a,b]\), and let \(\widetilde g\) be the natural cubic spline that interpolates \(g\) at the observed input values:

\[ \widetilde g(x_i)=g(x_i), \qquad i=1,\ldots,n. \]

Hence, \(g\) and \(\widetilde g\) have exactly the same empirical loss. Define their difference as

\[ h(x)=g(x)-\widetilde g(x), \]

so that \(h(x_i)=0\) for every \(i\). The key question is whether replacing \(g\) by \(\widetilde g\) can increase the roughness penalty.

Because \(\widetilde g\) is a natural cubic spline,

\[ \widetilde g''(a)=\widetilde g''(b)=0, \]

and \(\widetilde g^{(3)}(x)=c_j\) is constant on each interval \((x_j,x_{j+1})\). Integration by parts gives

\[\begin{align} \int_a^b\widetilde g''(x)h''(x)\,dx &= \left[\widetilde g''(x)h'(x)\right]_a^b - \sum_{j=1}^{n-1} \int_{x_j}^{x_{j+1}}\widetilde g^{(3)}(x)h'(x)\,dx \\ &= -\sum_{j=1}^{n-1} c_j\left\{h(x_{j+1})-h(x_j)\right\} \\ &=0. \end{align}\]

The boundary term is zero because of the natural boundary conditions. The last line is zero because \(h(x_i)=0\) at every observed input. Therefore,

\[\begin{align} \int_a^b[g''(x)]^2\,dx &= \int_a^b[\widetilde g''(x)+h''(x)]^2\,dx \\ &= \int_a^b[\widetilde g''(x)]^2\,dx + \int_a^b[h''(x)]^2\,dx \\ &\geq \int_a^b[\widetilde g''(x)]^2\,dx. \end{align}\]

If \(J_\lambda(g)\) denotes the objective in Equation (2.1), we have

\[ J_\lambda(g) = J_\lambda(\widetilde g) + \lambda\int_a^b[h''(x)]^2\,dx \geq J_\lambda(\widetilde g). \]

This is the important point. Every candidate function can be replaced by a natural cubic spline with the same fitted values and no larger roughness. Therefore, it is enough to search for the minimizer within the finite-dimensional natural cubic spline space.

2.6.2 Optional Theory: The Finite-Dimensional Solution

The smoothing-spline estimator has the representation

\[ \widehat g_\lambda(x) = \sum_{j=1}^n\widehat\beta_jN_j(x), \]

where \(N_1,\ldots,N_n\) form a natural cubic spline basis with knots at the observed input values. Define

\[ F_{ij}=N_j(x_i) \]

and define the roughness penalty matrix by

\[ \Omega_{jk} = \int_a^bN_j''(x)N_k''(x)\,dx. \]

Then Equation (2.1) becomes

\[\begin{equation} J_\lambda(\bbeta) = \frac{1}{n}\lVert\by-\bF\bbeta\rVert^2 + \lambda\bbeta^\top\Omega\bbeta. \tag{2.2} \end{equation}\]

Notice that \(\Omega\) is a symmetric positive semidefinite penalty matrix, since

\[ \bbeta^\top\Omega\bbeta = \int_a^b \left\{ \sum_{j=1}^n\beta_jN_j''(x) \right\}^2dx \geq 0. \]

Taking the derivative with respect to \(\bbeta\), we obtain

\[\begin{align} \bzero &= -\frac{2}{n}\bF^\top(\by-\bF\widehat{\bbeta}) + 2\lambda\Omega\widehat{\bbeta}, \\ (\bF^\top\bF+n\lambda\Omega)\widehat{\bbeta} &= \bF^\top\by. \end{align}\]

Therefore, when the matrix is invertible,

\[ \widehat{\bbeta} = (\bF^\top\bF+n\lambda\Omega)^{-1}\bF^\top\by. \]

This is a generalized ridge-regression solution. Although we started by optimizing over a large space of functions, the solution has a finite representation determined by the observed \(x_i\)’s. This classical connection between smoothing penalties and finite representations is developed in Kimeldorf and Wahba (1970). This result is the spline prototype for the general result introduced in The Representer Theorem: an infinite-dimensional optimization problem can have a finite-dimensional solution determined by the observed data. Both arguments rely on directions that leave the fitted values unchanged and cannot improve the penalty.

For a fixed \(\lambda\), the fitted values are linear in the observed responses:

\[ \widehat{\by} = \bS_\lambda\by, \qquad \bS_\lambda = \bF(\bF^\top\bF+n\lambda\Omega)^{-1}\bF^\top. \]

The effective degrees of freedom are \(\operatorname{tr}(\bS_\lambda)\). A common choice of \(\lambda\) minimizes the generalized cross-validation criterion

\[ \operatorname{GCV}(\lambda) = \frac{n^{-1}\lVert\by-\widehat{\by}\rVert^2} {\left\{1-\operatorname{tr}(\bS_\lambda)/n\right\}^2}. \]

As \(\lambda\) increases, the effective degrees of freedom decrease and the fitted function becomes smoother.

2.7 Fitting Smoothing Splines

Fitting a smoothing spline can be done using the smooth.spline() function in R. By default, the function uses GCV to select the smoothing parameter. For the birth-rate data, GCV selects a relatively large number of effective degrees of freedom, so the curve follows much of the year-to-year variation. Let’s compare it with a deliberately smoother fit.

    par(mfrow = c(1, 1), mar = c(4, 4, 2, 1))

    fit.gcv = smooth.spline(birthrates$Year, birthrates$Birthrate)
    fit.df8 = smooth.spline(birthrates$Year, birthrates$Birthrate, df = 8)
    year_grid = seq(1917, 2003)

    plot(birthrates$Year, birthrates$Birthrate, pch = 19,
         xlab = "Year", ylab = "Birth rate", col = "darkorange")
    lines(year_grid, predict(fit.gcv, year_grid)$y,
          col = "deepskyblue", lwd = 3)
    lines(year_grid, predict(fit.df8, year_grid)$y,
          col = "darkgreen", lwd = 3)
    legend("topright", c("GCV fit", "Fit with df = 8"),
           col = c("deepskyblue", "darkgreen"),
           lty = 1, lwd = 3, bty = "n")
    title("The Effect of Smoothing")


    c("GCV degrees of freedom" = fit.gcv$df,
      "Fixed degrees of freedom" = fit.df8$df)
##   GCV degrees of freedom Fixed degrees of freedom 
##                60.769100                 7.998551

The green curve is smoother because its effective degrees of freedom are smaller. This comparison illustrates the main role of \(\lambda\), but one realization alone does not tell us which curve has the best predictive performance.

Let’s look at a simulation example where the true function is known.

    set.seed(1)
    n = 100
    x = seq(0, 1, length.out = n)
    g_true = function(x) sin(12*(x+0.2))/(x+0.2)
    y = g_true(x) + rnorm(n)

    fit = smooth.spline(x, y)

    par(mfrow = c(1, 1), mar = c(4, 4, 2, 1))
    plot(x, y, pch = 19, xlim = c(0, 1),
         xlab = "x", ylab = "y", col = "darkorange")
    lines(x, g_true(x), col = "red", lwd = 3)
    lines(x, predict(fit, x)$y, col = "deepskyblue", lwd = 3)
    legend("bottomright", c("Truth", "Smoothing spline"),
           col = c("red", "deepskyblue"), lty = 1,
           lwd = 3, cex = 1.1, bty = "n")
    title("Smoothing Spline with Known Truth")


    fit$df
## [1] 9.96443

The fitted curve recovers the main structure of the true function while smoothing over part of the observation noise. The figure illustrates the behavior for one simulated data set. It is not, by itself, a general statement about prediction error.

3 Kernel Basis and Regularization

3.1 Example: Kernel Basis

Natural cubic splines give us one particular family of local basis functions. A natural question is whether we can replace them with another family of local functions while keeping a similar computation. Let’s consider Gaussian basis functions and see what happens. Here, the word kernel refers to a localized basis function. The conditions under which a kernel defines an RKHS are introduced in RKHS and Kernel Functions.

Be aware that using a kernel as a basis function is different from using a kernel smoother, such as the Nadaraya-Watson estimator. To be concrete, with four knots \(\xi_1,\ldots,\xi_4\), we consider the candidate function

\[ f_{\bbeta}(x) = \beta_0+\sum_{j=1}^4\beta_jK_\sigma(x,\xi_j). \]

Hence, the empirical loss is

\[ L(\bbeta) = \frac{1}{n}\sum_{i=1}^n \left[ y_i-\beta_0-\sum_{j=1}^4\beta_jK_\sigma(x_i,\xi_j) \right]^2. \]

To match the dnorm() implementation below, we use

\[ K_\sigma(x,\xi) = \frac{1}{\sqrt{2\pi}\sigma} \exp\left\{-\frac{(x-\xi)^2}{2\sigma^2}\right\}, \]

with \(\sigma=0.5\). The normalizing constant could be absorbed into the coefficients, but the intercept remains a separate part of the fitted function. To implement this numerically, we will

  • define a set of knots \(\xi_j\), which are the centers of the Gaussian basis functions;
  • evaluate the basis functions at the observed inputs to construct a design matrix;
  • perform linear regression using these newly constructed covariates;
  • evaluate the same basis functions on a plotting grid to obtain the fitted curve.

For practice, let’s use four knots selected from the empirical quantiles.

  set.seed(123)
  n = 50
  f_true <- function(x) 0.5 + sin(x)
  x <- seq(-pi, pi, length.out = n)
  y <- f_true(x) + rnorm(length(x), sd = 0.2)
  train_data <- data.frame(x = x, y = y)

  nknots = 4
  knots <- quantile(train_data$x, probs = (1:nknots)/(nknots + 1))

  # Gaussian basis functions with sigma = 0.5
  basis <- function(x, knots, sigma = 0.5) {
    sapply(knots, function(knot) dnorm(x, mean = knot, sd = sigma))
  }

  design_matrix <- matrix(NA, nrow = n, ncol = nknots + 1)
  design_matrix[, 1] <- 1
  for (j in 1:nknots) {
    design_matrix[, j + 1] <- basis(train_data$x, knots[j])
  }

  par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
  plot(train_data$x, train_data$y, pch = 19, col = "blue",
       xlab = "x", ylab = "y", main = "Observed Data")

  matplot(train_data$x, design_matrix[, -1], type = "l",
          lty = 1, lwd = 2, xlab = "x",
          ylab = expression(K[sigma](x, xi[j])),
          main = "Gaussian Basis Functions")

The fitted function can now be written as a sum of the estimated basis contributions.

  beta_hat = solve(
      t(design_matrix) %*% design_matrix,
      t(design_matrix) %*% train_data$y
  )

  x_grid = seq(-pi, pi, length.out = 200)
  grid_design_matrix <- matrix(NA, nrow = length(x_grid), ncol = nknots + 1)
  grid_design_matrix[, 1] <- 1
  for (j in 1:nknots) {
    grid_design_matrix[, j + 1] <- basis(x_grid, knots[j])
  }

  y_hat <- grid_design_matrix %*% beta_hat
  components <- sweep(
      grid_design_matrix[, -1, drop = FALSE],
      2,
      as.numeric(beta_hat[-1]),
      FUN = "*"
  )

  y_range <- range(c(train_data$y, y_hat, components, beta_hat[1]))
  plot(train_data$x, train_data$y, pch = 19, col = "blue",
       ylim = y_range, xlab = "x", ylab = "y",
       main = "Basis Contributions and Fitted Function")
  matlines(x_grid, components, lty = 1, lwd = 2)
  abline(h = beta_hat[1], lty = 2, lwd = 2)
  lines(x_grid, y_hat, col = "black", lwd = 4)
  legend("topright", c("Intercept", "Fitted function"),
         col = c("black", "black"), lty = c(2, 1),
         lwd = c(2, 4), bty = "n")

Each colored curve is one fitted basis contribution \(\widehat\beta_jK_\sigma(x,\xi_j)\). The dashed line is the intercept, and the thick black curve adds the intercept and all four contributions. This makes the finite basis representation visible.

3.2 Increasing the Number of Basis Functions

As the functional form becomes more complicated, we can increase the number of basis functions. There are several questions we could ask while performing this task:

  • What basis or kernel function should we use?
  • How many knots should be used, and where should they be located?
  • How do we prevent overfitting when the number of knots increases?

3.2.1 From Many Basis Functions to Regularization

In this final step, let’s think about what happens when the number of knots is as large as, or even greater than, the number of observations. The design matrix may be rank deficient or nearly so, and ordinary least squares becomes unstable. A natural way to deal with this is ridge regularization.

Let \(\bB\) denote the design matrix, including an intercept in its first column. We solve

\[ \widehat{\bbeta} = \underset{\bbeta}{\arg\min} \left\{ \frac{1}{n}\lVert\by-\bB\bbeta\rVert^2 + \lambda\bbeta^\top\bP\bbeta \right\}, \qquad \bP=\operatorname{diag}(0,1,\ldots,1). \]

The leading zero means that the intercept is not penalized.

  set.seed(123)
  n = 100
  f_true <- function(x) 0.5 + sin(x)
  x <- seq(-2*pi, 2*pi, length.out = n)
  y <- f_true(x) + rnorm(length(x), sd = 0.2)
  train_data <- data.frame(x = x, y = y)

  # place a Gaussian basis function at every observed input
  nknots = n
  knots <- sort(train_data$x)

  basis <- function(x, knots, sigma = 0.5) {
    sapply(knots, function(knot) dnorm(x, mean = knot, sd = sigma))
  }

  design_matrix <- matrix(NA, nrow = n, ncol = nknots + 1)
  design_matrix[, 1] <- 1
  for (j in 1:nknots) {
    design_matrix[, j + 1] <- basis(train_data$x, knots[j])
  }

  # ridge fit with an unpenalized intercept
  lambda = 0.01
  penalty_matrix = diag(c(0, rep(1, nknots)))
  gram_matrix = t(design_matrix) %*% design_matrix / n
  rhs = t(design_matrix) %*% train_data$y / n
  beta_hat = solve(gram_matrix + lambda * penalty_matrix, rhs)

  x_grid = seq(-2*pi, 2*pi, length.out = 200)
  grid_design_matrix <- matrix(NA, nrow = length(x_grid), ncol = nknots + 1)
  grid_design_matrix[, 1] <- 1
  for (j in 1:nknots) {
    grid_design_matrix[, j + 1] <- basis(x_grid, knots[j])
  }

  y_hat <- grid_design_matrix %*% beta_hat

  plot(train_data$x, train_data$y, pch = 19, col = "blue",
       xlim = c(-2*pi, 2*pi), xlab = "x", ylab = "y",
       main = "Regularized Fitted Function")
  lines(x_grid, f_true(x_grid), col = "darkorange", lwd = 3)
  lines(x_grid, y_hat, col = "black", lwd = 4)
  legend("topleft", c("Truth", "Regularized fit"),
         col = c("darkorange", "black"), lty = 1,
         lwd = c(3, 4), bty = "n")

Here, \(\lambda=0.01\) is used only to illustrate stabilization. We have not tuned it for prediction. Also, this is ordinary ridge regression on the coefficients of a chosen basis. It is not yet kernel ridge regression, where the penalty is defined through an RKHS norm.

Basis expansion turned nonlinear function estimation into a linear coefficient problem. Splines made the construction local, while regularization allowed us to use many basis functions without fitting the noise directly.

3.2.2 Optional: A Special Sobolev Kernel

For an optional technical comparison, let’s see what can happen when we use many basis functions without regularization. Consider the localized function

\[ K(x,\xi) = \frac{\sqrt{3}}{3} \exp\left(-\frac{\sqrt{3}}{2}|x-\xi|\right) \sin\left(\frac{|x-\xi|}{2}+\frac{\pi}{6}\right). \]

With an appropriate inner product and boundary conditions, second-order Sobolev spaces can be formulated as RKHSs. The exact connection depends on the domain and the norm, so we will not use this special formula in the main development. See Berlinet and Thomas-Agnan (2011), Section 6.1.6, for more details.

The following code illustrates how 40 copies of this localized function can be combined.

  set.seed(123)
  n = 100
  f_true <- function(x) 0.5 + sin(x)
  x <- seq(-2*pi, 2*pi, length.out = n)
  y <- f_true(x) + rnorm(length(x), sd = 0.2)
  train_data <- data.frame(x = x, y = y)

  nknots = 40
  knots <- quantile(train_data$x, probs = (1:nknots)/(nknots + 1))

  sobolev_basis <- function(x, knots) {
    sapply(knots, function(knot) {
      (sqrt(3)/3) * exp(-sqrt(3) * abs(x - knot) / 2) *
      sin(abs(x - knot) / 2 + pi/6)
    })
  }

  sobolev_design <- matrix(NA, nrow = n, ncol = nknots + 1)
  sobolev_design[, 1] <- 1
  for (j in 1:nknots) {
    sobolev_design[, j + 1] <- sobolev_basis(train_data$x, knots[j])
  }

  beta_sobolev = solve(
      t(sobolev_design) %*% sobolev_design,
      t(sobolev_design) %*% train_data$y
  )

  x_grid = seq(-2*pi, 2*pi, length.out = 200)
  sobolev_grid <- matrix(NA, nrow = length(x_grid), ncol = nknots + 1)
  sobolev_grid[, 1] <- 1
  for (j in 1:nknots) {
    sobolev_grid[, j + 1] <- sobolev_basis(x_grid, knots[j])
  }

  y_hat <- sobolev_grid %*% beta_sobolev

  plot(train_data$x, train_data$y, pch = 19, col = "blue",
       xlim = c(-2*pi, 2*pi), xlab = "x", ylab = "y",
       main = "Fitted Function with 40 Basis Functions")
  lines(x_grid, f_true(x_grid), col = "darkorange", lwd = 3)
  lines(x_grid, y_hat, col = "black", lwd = 3)
  legend("topleft", c("Truth", "Fitted function"),
         col = c("darkorange", "black"), lty = 1,
         lwd = 3, bty = "n")

The fitted curve follows some of the observation noise. This is the instability that the ridge penalty is intended to control.

4 Extending Splines to Multiple Variables

So far, we have focused on a regression function with one input variable. In many applications, the response depends on several variables. A natural first step is to preserve the one-dimensional spline construction and combine several smooth functions.

4.1 Additive Models with Splines

Since spline models can be written as linear models in transformed covariates, an additive structure gives

\[ f(\bx) = \beta_0+ \sum_{j=1}^p h_j(x_j) = \beta_0+ \sum_{j=1}^p\sum_{k=1}^m N_{jk}(x_j)\beta_{jk}, \]

Here \(h_j(x_j)\) is a univariate function approximated using spline basis functions. To make the component functions identifiable, we impose the centering constraints

\[ \sum_{i=1}^n h_j(x_{ij})=0, \qquad j=1,\ldots,p. \]

In practice, additive-model software uses an equivalent constraint or basis parameterization. The additive model uses roughly \(pm\) basis functions and can be fitted as either a linear or generalized linear model.

This makes the additive model computationally convenient. However, the additive assumption also has an important restriction: the effect of \(x_1\) cannot change depending on the value of \(x_2\).

4.1.1 Optional Computation: A Logistic Additive Model

The following optional example uses gam() to fit a logistic additive model to the diabetes data in the recommended MASS package. We use natural-spline terms for glucose, blood pressure, body mass index, and age.

    library(gam)
    data("Pima.tr", package = "MASS")

    form = formula("type ~ splines::ns(glu, df=4) +
                           splines::ns(bp, df=4) +
                           splines::ns(bmi, df=4) +
                           splines::ns(age, df=4)")

    m = gam(form, data = Pima.tr, family = binomial)
    summary(m)
## 
## Call: gam(formula = form, family = binomial, data = Pima.tr)
## Deviance Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.0215 -0.6707 -0.3136  0.7039  2.6823 
## 
## (Dispersion Parameter for binomial family taken to be 1)
## 
##     Null Deviance: 256.4142 on 199 degrees of freedom
## Residual Deviance: 175.5097 on 183 degrees of freedom
## AIC: 209.5097 
## 
## Number of Local Scoring Iterations: 6 
## 
## Anova for Parametric Effects
##                           Df  Sum Sq Mean Sq F value    Pr(>F)    
## splines::ns(glu, df = 4)   4  21.360  5.3401  4.9518 0.0008192 ***
## splines::ns(bp, df = 4)    4   2.201  0.5502  0.5102 0.7283355    
## splines::ns(bmi, df = 4)   4   9.229  2.3072  2.1395 0.0776878 .  
## splines::ns(age, df = 4)   4  13.107  3.2769  3.0386 0.0186497 *  
## Residuals                183 197.351  1.0784                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
    plot(m, se = TRUE, residuals = TRUE,
         pch = 19, col = "darkorange")

Each panel shows one estimated component function while holding the other components fixed. This makes the additive structure easy to interpret. But it still does not allow two variables to interact.

4.2 Tensor-Product Basis Functions

How can we allow the effect of one variable to depend on another variable? Suppose \(\{B_k(x_1)\}_{k=1}^K\) and \(\{C_\ell(x_2)\}_{\ell=1}^L\) are two univariate bases. Their tensor-product basis functions are

\[ \phi_{k\ell}(x_1,x_2) = B_k(x_1)C_\ell(x_2), \qquad k=1,\ldots,K, \quad \ell=1,\ldots,L. \]

Therefore, \(K\) basis functions in the first direction and \(L\) basis functions in the second direction produce \(KL\) basis functions. At a point \((x_1,x_2)\), the complete set of basis values is the outer product

\[ \mathbf b(x_1)\mathbf c(x_2)^\top, \]

where \(\mathbf b(x_1)=(B_1(x_1),\ldots,B_K(x_1))^\top\) and \(\mathbf c(x_2)=(C_1(x_2),\ldots,C_L(x_2))^\top\).

Before using spline bases in both directions, it is useful to visualize the geometry with two simple localized basis functions. For the first variable, consider a Gaussian kernel basis centered at \(\xi\):

\[ G_{\xi,\sigma}(x_1) = \exp\left\{-\frac{(x_1-\xi)^2}{2\sigma^2}\right\}. \]

For the second variable, consider a triangular kernel basis centered at \(\zeta\):

\[ T_{\zeta,h}(x_2) = \left(1-\frac{|x_2-\zeta|}{h}\right)_+, \qquad (u)_+=\max(u,0). \]

Their product is

\[ \phi(x_1,x_2) = G_{\xi,\sigma}(x_1)T_{\zeta,h}(x_2). \]

This is a tensor-product basis function, but it is not yet a tensor-product spline because the Gaussian basis is not a spline. We use these two different shapes so that the product construction is easy to see.

    library(plotly)

    x1_grid = seq(0, 1, length.out = 61)
    x2_grid = seq(0, 1, length.out = 47)

    gaussian_basis = exp(-0.5 * ((x1_grid - 0.35) / 0.10)^2)
    triangular_basis = pmax(1 - abs(x2_grid - 0.65) / 0.20, 0)

    gaussian_surface = outer(
        rep(1, length(x2_grid)), gaussian_basis
    )
    triangular_surface = outer(
        triangular_basis, rep(1, length(x1_grid))
    )
    product_surface = outer(triangular_basis, gaussian_basis)

    scene_axes = list(
        xaxis = list(title = "x1"),
        yaxis = list(title = "x2"),
        zaxis = list(title = "Basis value", range = c(0, 1))
    )

    plot_ly(height = 420) %>%
        add_surface(
            x = x1_grid, y = x2_grid, z = gaussian_surface,
            scene = "scene", showscale = FALSE,
            colorscale = "Viridis"
        ) %>%
        add_surface(
            x = x1_grid, y = x2_grid, z = triangular_surface,
            scene = "scene2", showscale = FALSE,
            colorscale = "Viridis"
        ) %>%
        add_surface(
            x = x1_grid, y = x2_grid, z = product_surface,
            scene = "scene3", showscale = FALSE,
            colorscale = "Viridis"
        ) %>%
        layout(
            paper_bgcolor = "transparent",
            plot_bgcolor = "rgb(254, 247, 234)",
            scene = modifyList(
                scene_axes,
                list(domain = list(x = c(0, 0.31), y = c(0, 1)))
            ),
            scene2 = modifyList(
                scene_axes,
                list(domain = list(x = c(0.345, 0.655), y = c(0, 1)))
            ),
            scene3 = modifyList(
                scene_axes,
                list(domain = list(x = c(0.69, 1), y = c(0, 1)))
            ),
            annotations = list(
                list(text = "Gaussian basis", x = 0.155, y = 1.02,
                     xref = "paper", yref = "paper", showarrow = FALSE),
                list(text = "Triangular basis", x = 0.500, y = 1.02,
                     xref = "paper", yref = "paper", showarrow = FALSE),
                list(text = "Product basis", x = 0.845, y = 1.02,
                     xref = "paper", yref = "paper", showarrow = FALSE)
            ),
            margin = list(l = 0, r = 0, b = 0, t = 35)
        )

The first two surfaces extend a univariate basis function along the other coordinate. Their product is large only when both marginal basis functions are large. Along the \(x_1\) direction, the product keeps the smooth Gaussian shape. Along the \(x_2\) direction, it keeps the compact, piecewise-linear triangular shape. The resulting surface is therefore localized in both directions.

4.3 Tensor-Product Spline Regression

We now take both marginal bases to be spline bases. Let \(\{B_k(x_1)\}_{k=1}^K\) and \(\{C_\ell(x_2)\}_{\ell=1}^L\) be two B-spline bases. A tensor-product spline model is

\[ f_{\boldsymbol\theta}(x_1,x_2) = \sum_{k=1}^K\sum_{\ell=1}^L \theta_{k\ell}B_k(x_1)C_\ell(x_2) = \mathbf b(x_1)^\top \boldsymbol\Theta \mathbf c(x_2), \]

where \(\boldsymbol\Theta\) is a \(K\times L\) coefficient matrix. For observation \(i\), the design matrix contains the rowwise products

\[ \boldsymbol\Phi_{i,(k,\ell)} = B_k(x_{i1})C_\ell(x_{i2}). \]

Because both marginal B-spline bases contain a constant direction, their tensor product already contains an intercept. We do not add a separate intercept to the model.

Therefore, fitting a tensor-product spline is still a linear regression problem:

\[ \widehat{\boldsymbol\theta} = \underset{\boldsymbol\theta}{\arg\min} \frac{1}{n} \left\|\by-\boldsymbol\Phi\boldsymbol\theta\right\|_2^2. \]

4.3.1 Optional Computation: Fitting a Tensor-Product Spline

The following simulation uses the nonadditive regression function

\[ f_0(x_1,x_2) = \sin(2\pi x_1)\cos(2\pi x_2). \]

For a fixed \(x_2\), the shape in the \(x_1\) direction changes with \(x_2\). An additive model cannot represent this interaction. We construct the complete tensor-product design matrix by hand and fit its coefficients using least squares.

    library(splines)

    set.seed(546)
    n = 350
    x1 = runif(n)
    x2 = runif(n)

    f_true = function(x1, x2) {
        sin(2*pi*x1) * cos(2*pi*x2)
    }

    y = f_true(x1, x2) + rnorm(n, sd = 0.15)

    B1 = bs(
        x1, df = 7, degree = 3, intercept = TRUE,
        Boundary.knots = c(0, 1)
    )
    B2 = bs(
        x2, df = 7, degree = 3, intercept = TRUE,
        Boundary.knots = c(0, 1)
    )

    K = ncol(B1)
    L = ncol(B2)

    tensor_design = t(vapply(
        seq_len(n),
        function(i) as.vector(outer(B1[i, ], B2[i, ])),
        numeric(K*L)
    ))

    fit = lm.fit(x = tensor_design, y = y)
    theta_hat = matrix(fit$coefficients, nrow = K, ncol = L)

    x1_grid = seq(0, 1, length.out = 61)
    x2_grid = seq(0, 1, length.out = 47)
    B1_grid = predict(B1, x1_grid)
    B2_grid = predict(B2, x2_grid)

    fitted_surface = t(B1_grid %*% theta_hat %*% t(B2_grid))
    true_surface = outer(
        x2_grid, x1_grid,
        function(x2, x1) f_true(x1, x2)
    )

    grid_rmse = sqrt(mean((fitted_surface - true_surface)^2))
    round(grid_rmse, 3)
## [1] 0.092

    z_range = range(c(true_surface, fitted_surface))
    fit_scene_axes = list(
        xaxis = list(title = "x1"),
        yaxis = list(title = "x2"),
        zaxis = list(title = "f(x1, x2)", range = z_range)
    )

    plot_ly(height = 430) %>%
        add_surface(
            x = x1_grid, y = x2_grid, z = true_surface,
            scene = "scene", showscale = FALSE,
            colorscale = "Viridis", cmin = z_range[1], cmax = z_range[2]
        ) %>%
        add_surface(
            x = x1_grid, y = x2_grid, z = fitted_surface,
            scene = "scene2", showscale = FALSE,
            colorscale = "Viridis", cmin = z_range[1], cmax = z_range[2]
        ) %>%
        layout(
            paper_bgcolor = "transparent",
            plot_bgcolor = "rgb(254, 247, 234)",
            scene = modifyList(
                fit_scene_axes,
                list(domain = list(x = c(0, 0.47), y = c(0, 1)))
            ),
            scene2 = modifyList(
                fit_scene_axes,
                list(domain = list(x = c(0.53, 1), y = c(0, 1)))
            ),
            annotations = list(
                list(text = "True regression surface", x = 0.235, y = 1.02,
                     xref = "paper", yref = "paper", showarrow = FALSE),
                list(text = "Tensor-product spline fit", x = 0.765, y = 1.02,
                     xref = "paper", yref = "paper", showarrow = FALSE)
            ),
            margin = list(l = 0, r = 0, b = 0, t = 35)
        )

The fitted surface recovers the main interaction pattern from one noisy sample. This is only a simulation illustration, not a general guarantee. We have used unpenalized least squares so that the tensor-product construction remains visible. In practice, we may also add separate roughness penalties in the two input directions.

4.4 The High-Dimensional Limitation

The progression from polynomials to splines has shown how local basis functions and regularization can produce flexible estimates. However, extending the same construction to several input variables creates a new difficulty. The additive model uses roughly \(pm\) basis functions, but it only includes additive effects. A complete tensor-product construction includes interactions, but if we use \(m\) basis functions for each of \(p\) input variables, it may require

\[ \underbrace{m\times\cdots\times m}_{p\text{ input variables}} = m^p \]

basis functions. This quickly becomes computationally infeasible. An RKHS formulation can avoid explicitly constructing all \(m^p\) tensor-product basis functions, although it does not eliminate the statistical difficulty of high-dimensional estimation or the computational cost of working with an \(n\times n\) kernel matrix. Can we construct a more general space of functions whose solution still has a finite, data-dependent representation? This is the question that motivates RKHS.

Berlinet, Alain, and Christine Thomas-Agnan. 2011. Reproducing Kernel Hilbert Spaces in Probability and Statistics. Springer Science & Business Media.
Kimeldorf, George S, and Grace Wahba. 1970. “A Correspondence Between Bayesian Estimation on Stochastic Processes and Smoothing by Splines.” The Annals of Mathematical Statistics 41 (2): 495–502.