# General P-splines for non-uniform B-splines

Zheyuan<EMAIL_ADDRESS>[ Jiguo<EMAIL_ADDRESS>[
Zheyuan Li is Postdoctoral Fellow, Department of Statistics and Actuarial
Science, Simon Fraser University, Burnaby, BC, Canada .  Jiguo Cao is
Professor, Department of Statistics and Actuarial Science, Simon Fraser
University, Burnaby, BC, Canada .

###### Abstract

P-spline represents an unknown univariate function with uniform B-splines on
equidistant knots and penalizes their coefficients using a simple difference
matrix for smoothness. But for non-uniform B-splines on unevenly spaced knots,
such difference penalty fails, and the conventional derivative penalty is
hitherto the only choice. We proposed a general P-spline estimator to lift
this restriction by deriving a general difference penalty for non-uniform
B-splines. We also established a sandwich formula between derivative and
general difference penalties for a better understanding of their connections.
Simulations show that both P-spline variants have close MSE performance in
general. But in practice, one can yield a more satisfactory fit than the
other. For example, the bone mineral content (BMC) data favor general
P-spline, while the fossil shell data favor standard P-spline. We therefore
believe both variants to be useful tools for practical modeling. To implement
our general P-spline, we developed two new R packages: gps and gps.mgcv. The
latter creates a new “gps” smooth class for mgcv, so that a general P-spline
can be specified as s(x, bs = "gps") in a model formula and estimated in the
framework of generalized additive models.

P-splines,

B-splines,

difference penalty,

penalized regression,

generalized additive models,

###### keywords:

and

##  1 Introduction

P-spline (Eilers and Marx, 1996) proves very popular and successful, having
been applied in various statistical modeling frameworks (Eilers, Marx and
Durbán, 2015; Eilers and Marx, 2021). It represents an unknown univariate
function as a linear combination of basis called B-splines (de Boor, 2001),
and penalizes differences of nearby coefficients for its smoothness. The
appearance of B-splines depends on their knot spacing. Figure 1 shows cubic
B-splines on three types of knots: (a) uniform B-splines on equidistant knots;
(b) non-uniform B-splines on unevenly spaced knots; (c) clamped non-uniform
B-splines on unevenly spaced knots and repeated boundary knots. Note that (a)
uniform B-splines differ from each other by a horizontal shift so they have
identical shapes; (b) non-uniform B-splines have inhomogeneous shapes; (c) the
first/last clamped B-spline is open at the left/right boundary with a jump
discontinuity, while other B-splines are closed as usual, with each tail tied
to a knot. Clearly, B-splines only remain nonzero between a few adjacent
knots, thus they have local support.

Figure 1: (a) uniform B-splines with equidistant knots; (b) non-uniform
B-splines with unevenly spaced knots; (c) non-uniform and clamped B-splines
with unevenly spaced knots and in particular, duplicated boundary knots. In
each graph, there are 12 knots (black dots, stacked up if duplicated) and 8
cubic B-splines (in different colours).

Although Eilers and Marx have reiterated that uniform B-splines with
equidistant knots are all that is needed for P-spline, non-uniform B-splines
with unevenly spaced knots still turn out useful in reality. In this paper, we
propose a new P-spline variant that accommodates non-uniform B-splines,
hereafter called the general P-spline. By contrast, Eilers and Marx’s P-spline
based on uniform B-splines is referred to as the standard P-spline.

Figure 2: (A) BMC measurements (in grams) from 112 male subjects. (B) log-
transformed measurements. (C) standard P-spline estimates (with 10 uniform
cubic B-splines) of population BMC trajectory ($\exp(\hat{f})$, solid) and two
example subject trajectories ($\exp(\hat{f}+\hat{f}_{j})$, dashed). (D)
general P-spline estimates (with 10 non-uniform and clamped cubic B-splines)
for the same trajectories in (C). In (C) and (D), gray dots are data for the
two selected subjects and gray vertical dashed lines mark locations of domain
knots (see Figure 4 for more about knots).

Our general P-spline is motivated by the analysis of bone mineral content
(BMC) longitudinal data from the Pediatric Bone Mineral Accrual Study (PBMAS)
(Bailey, 1997), which aimed to investigate BMC accumulation in growing
children, with boys and girls aged from 8 to 15 attending two elementary
schools in the city of Saskatoon enrolled for annual scan since 1991. See
Figure 2(A) for illustration of the data from the male group. For either
gender, we build the following generalized additive mixed model (see Zuur,
Saveliev and Ieno (2014) for a thorough introduction):

$\log(\textrm{bmc})_{j,\textrm{age}}=f(\textrm{age})+f_{j}(\textrm{age})+e_{j,\textrm{age}},$

where $f(\textrm{age})$ (a fixed effect) is the population log(BMC)
trajectory, $f_{j}(\textrm{age})$ (a random effect) is the departure of
subject $j$’s log(BMC) from $f(\textrm{age})$ and $e_{j,\textrm{age}}$ is an
iid Gaussian error. Denote estimated functions by $\hat{f}$ and $\hat{f}_{j}$,
then estimated population and subject BMC trajectories are $\exp(\hat{f})$ and
$\exp(\hat{f}+\hat{f}_{j})$. As Figure 2(C) illustrates, standard P-spline
estimates of these trajectories exhibit suspicious wiggliness, whereas general
P-spline estimates in Figure 2(D) are very satisfactory.

Extending standard P-spline to general P-spline is not as naive as replacing
uniform B-splines by non-uniform B-splines in basis representation. The
difference penalty applied in the standard case must be adapted, too.
Likewise, we call the difference penalty associated with the standard/general
P-spline the standard/general difference penalty. It has not been said in
literature, but standard P-spline does not make sense for non-uniform knots.
As a result, derivative penalty as used in smoothing splines is hitherto the
only choice for penalizing non-uniform B-splines (Wood, 2017a). This
restriction is now lifted with the availability of our general difference
penalty.

Figure 3: Smoothing noisy observations (gray dots) from
$g(x)=\frac{1}{8}|x|^{3}$, using (a) smoothing spline equivalent, (b) standard
P-spline, (c) naive P-spline and (d) general P-spline. Gray vertical dashed
lines mark locations of knots.

To demonstrate the failure/success of standard/general difference penalty for
non-uniform B-splines, consider smoothing noisy observations with unevenly
spaced $x$-values on $[-3,3]$ from $g(x)=\frac{1}{8}|x|^{3}$ , using the
following combinations of basis and penalty: (a) non-uniform B-splines with
derivative penalty (equivalent to smoothing spline), (b) uniform B-splines
with standard difference penalty (standard P-spline), (c) non-uniform
B-splines with standard difference penalty (naive P-spline) and (d) non-
uniform B-splines with general difference penalty (general P-spline). We place
100 equidistant knots for uniform B-splines and 100 quantile knots (at equal
quantiles of $x$-values) for non-uniform ones. As Figure 3 shows, (c) fits
data poorly, while (a)(b)(d) appear equally good. Boxplot of mean squared
error (MSE) on 100 simulations confirms that such observation is not a random
chance, and in fact, our general P-spline has the best performance for this
particular example.

This paper has four major contributions. First of all, we proposed general
P-spline (section 4.3) to accommodate non-uniform B-splines on unevenly spaced
knots. It is a complement to Eilers and Marx’s standard P-spline that is only
sensible for uniform B-splines on equidistant knots. We implemented it via R
packages gps and gps.mgcv (supplementary materials), with satisfactory
application to the modeling of BMC data (section 6.1). Secondly, we discovered
a sandwich formula between the conventional derivative penalty and our general
difference penalty, so that the former can be calculated from the latter
(section 4.4). This not only expands our knowledge on these two types of
penalties, but also reveals a straightforward way to adapt derivative penalty
from $l^{2}$ to $l^{1}$ regularization while preserving sparsity (section
4.5). Thirdly, based on simulation studies with random examples (section 5.2),
we concluded that there is generally no MSE performance gap between using two
types of penalties or two types of B-splines, so that smoothing spline and
standard/general P-spline are expected to be equally good. However, for any
specific example, say the U-shaped curve in Figure 3, one may indeed be
slightly better than the others. Lastly, through two real data examples
(section 6), we brought more insight into the debate on the superiority of
uniform or non-uniform knots, hence that of standard or general P-spline.
Unlike Ruppert, Wand and Carroll (2009); Wand and Ormerod (2008); Eilers, Marx
and Durbán (2015) that argued for either choice, we claimed that the story has
two sides. On one hand, both P-spline variants tend to give comparable
performance when the number of knots is sufficiently large to mitigate
sensitivity in exact knot locations. On the other hand, when the number of
knots is small to moderate, one knot placement strategy may be more effective
than the other, hence one P-spline variant may produce a more satisfactory
fit. Therefore, both P-spline variants are useful from a practical
perspective.

The rest of this paper is planned as follows. Section 2 contains an in-depth
preliminary on spline and B-splines. Section 3 reviews spline estimators from
regression spline to standard P-spline. Section 4 first explains why standard
difference penalty goes wrong for non-uniform B-splines, then derives general
difference penalty and establishes its link with the conventional derivative
penalty. Section 5 conducts simulation studies to demonstrate the success of
general P-spline, with comparison to its competitors in terms of MSE
performance. Section 6 investigates BMC data and fossil shell data, showing
that standard or general P-spline can produce a more satisfactory fit than the
other on different datasets, thus they are both practically useful.
Supplementary materials (Li and Cao, 2022) include R packages gps and gps.mgcv
that implement general P-spline, as well as R code for all figures in this
paper.

##  2 A primer on polynomial spline and B-splines

An order-$d$ polynomial spline $f(x)$ comprises smoothly connected polynomial
segments of degree $d-1$, with $d-2$ continuous derivatives at their joints.
Given $k$ interior knots $a<s_{1}<s_{2}<\ldots<s_{k}<b$ in a domain $[a,b]$ of
interest, we may express this spline as:

$f(x)=p_{j}(x)=\sum_{q=0}^{d-1}c_{q,j}(x-s_{j})^{q},\kern
5.69054pts_{j}<x<s_{j+1},$ (1)

where $s_{0}=a$ and $s_{k+1}=b$ are domain boundary knots. For smoothness at
interior knots, it is required that

$p_{j-1}^{\scriptscriptstyle(m)}(s_{j}^{-})=p_{j}^{\scriptscriptstyle(m)}(s_{j}^{+}),\kern
5.69054pt\forall m=0,1,\ldots,d-2.$

In this piecewise polynomial parametrization, there are $d(k+1)$ coefficients
subject to $k(d-1)$ linear constraints, ending up with $p=k+d$ free
parameters.

Another way to represent $f(x)$ is to use B-splines. With the above interior
knots, $p$ order-$d$ B-splines $(\mathcal{B}_{j}(x))_{1}^{p}$ can be
constructed, giving a basis representation:

$f(x)=\sum_{j=1}^{p}\mathcal{B}_{j}(x)\beta_{j}$ (2)

with $p$ unconstrained coefficients $(\beta_{j})_{1}^{p}$.

The two representations are equivalent, although the B-spline approach (2) has
a computational advantage that all its coefficients are unconstrained, whereas
the piecewise polynomial (1) has more interpretable coefficients:
$c_{m,j}=\frac{1}{m!}f^{\scriptscriptstyle(m)}(s_{j}),\ \forall
m=0,1,\ldots,d-1.$ There is little point in arguing which is superior, but it
is common to use (1) for interpolation and (2) for regression.

B-splines may be confusing because their construction is not as transparent as
piecewise polynomial. One subtle issue is that while interior knots
$(s_{j})_{1}^{k}$ and domain boundaries $a$, $b$ are all that is needed for
constructing piecewise polynomial (1), $d-1$ extra auxiliary boundary knots
are required on each side of $[a,b]$ to construct
$(\mathcal{B}_{j}(x))_{1}^{p}$ in (2). The connection between two sets of
knots is illustrated in Figure 4. The upper deck shows the $k+2$ domain knots
required for piecewise polynomial parametrization, while the lower deck shows
the $K=k+2d$ knots $(t_{j})_{1}^{K}$ needed for the equivalent B-spline
representation. In principle, the $(d-1)$ auxiliary boundary knots on each
side can be arbitrarily positioned, as they do not affect the piecewise
polynomial representation on $[a,b]$. However, they do have impact on
B-spline’s shape. For example, duplicated boundary knots $t_{1}=\ldots=t_{d}$
and $t_{d+k+1}=\ldots=t_{2d+k}$ give rise to clamped B-splines that are open
at $a$ and $b$, as Figure 1(c) shows. (It is a good time to revisit this
figure for further understanding. In these B-spline demos, $k=4$, $d=4$,
$p=8$, $K=12$. Black dots mark $(t_{j})_{1}^{12}$ and the gray shaded area
highlights the domain $[a,b]$. Besides, (b) (c) have common interior knots but
different boundary ones.)

Figure 4: Knots for constructing piecewise polynomial (1) (upper deck) and
equivalent B-spline representation (2) (lower deck). B-spline knots may be
grouped into (i) interior knots and boundary knots; (ii) domain knots and
auxiliary knots. We prefer to (i) and denote the number of interior knots by
$k$, so the total number of knots is $K=k+2d$ and the number of B-splines is
$p=k+d$.

It is inconvenient to express B-splines in closed forms, because they are
constructed recursively. Given a full knot sequence $(t_{j})_{1}^{K}$, we
begin with $(K-1)$ order-1 B-splines:

$\mathcal{B}_{j,1}(x)=\begin{cases}1&x\in[t_{j},t_{j+1}]\\\
0&\textrm{otherwise},\end{cases}$

then sequentially construct higher-order B-splines from lower-order B-splines
using formula:

$\mathcal{B}_{j,d}(x)=(x-t_{j})\frac{\mathcal{B}_{j,d-1}(x)}{\gamma_{j,d-1}}+(t_{j+d}-x)\frac{\mathcal{B}_{j+1,d-1}(x)}{\gamma_{j+1,d-1}},$

where $\gamma_{j,q}=t_{j+q}-t_{j}$ is lag-$q$ difference of knots. Figure 5
(from bottom to top) illustrates construction of order-4 B-splines on 8 non-
uniform knots $(t_{j})_{1}^{8}$. An order-$d$ B-spline is nonzero on a span of
$d+1$ successive knots (or $d$ consecutive intervals). In consequence, only
$d$ order-$d$ B-splines are nonzero between two adjacent knots. In the figure,
B-splines that are zero on $[t_{4},t_{5}]$ are coloured gray. In this demo,
$[a,b]=[t_{4},t_{5}]$ and there is no interior knot, i.e., $k=0$.

Figure 5: Given 8 non-uniform knots (dots on the $x$-axis), 7 constant
B-splines $(\mathcal{B}_{j,1}(x))_{1}^{7}$ (not drawn), 6 linear B-splines
$(\mathcal{B}_{j,2}(x))_{1}^{6}$, 5 quadratic B-splines
$(\mathcal{B}_{j,3}(x))_{1}^{5}$ and 4 cubic B-splines
$(\mathcal{B}_{j,4}(x))_{1}^{4}$ are sequentially constructed. In this demo,
$[a,b]=[t_{4},t_{5}]$ (gray shaded) and $k=0$ (no interior knots).

We have to index a B-spline with two subscripts when working with B-splines of
different orders at the same time. We will do so in section 4.2, 4.3 when
deriving our general difference penalty. But when there is no such need in
other sections, we will drop subscript $d$. For example, when writing basis
representation (2), it is assumed that the B-spline order $d$ can be inferred
from the context.

##  3 From regression spline to standard P-spline

Suppose we have observed data $(x_{i},y_{i})_{1}^{n}$ with
$x_{1}<x_{2}<\ldots<x_{n}$ from a generative model:

$y_{i}=g(x_{i})+\varepsilon_{i},\kern
5.69054pt\varepsilon_{i}\overset{\scriptscriptstyle\textrm{iid}}{\sim}N(0,\sigma^{2}),$

where $g(x)$ is a signal function and $\varepsilon_{i}$ is a Gaussian white
noise. A scatterplot smoothing model hypothesizes that $(x_{i},y_{i})_{1}^{n}$
comes from:

$y_{i}=f(x_{i})+e_{i},\kern
5.69054pte_{i}\overset{\scriptscriptstyle\textrm{iid}}{\sim}N(0,\sigma_{e}^{2}),$

where $f(x)$ belongs to some class of smooth function and $e_{i}$ is an iid
Gaussian error. In general, $g(x)$ is true but unknown, and $f(x)$ is an
estimator of $g(x)$. For example, if $f(x)$ is a linear function of $x$, i.e.,
$f(x)=\beta_{0}+\beta_{1}x$, then the smoothing model reduces to linear
regression. For nonlinear regression, it may be assumed that $f(x)$ is a
polynomial spline.

In this section, we review several polynomial spline estimators of $g(x)$,
including the simplest regression B-splines (section 3.1), the most
complicated smoothing spline (section 3.2), their crossbreed, that is,
penalized regression B-splines with derivative penalty (section 3.3) and its
adaptation with standard difference penalty, i.e., the standard P-spline
(section 3.4). Finally in section 3.5, we describe automatic knot placement
strategies for both types of penalized B-splines.

###  3.1 Regression B-splines

After choosing an order $d$ and placing $k$ interior knots
$a<s_{1}<s_{2}<\ldots<s_{k}<b$ (as well as auxiliary boundary knots), we can
express $f(x)$ using basis form (2), with coefficients $(\beta_{j})_{1}^{p}$
to estimate. Least squares regression B-splines minimize:

$\sum_{i=1}^{n}(y_{i}-f(x_{i}))^{2}=\|\boldsymbol{y}-\boldsymbol{B\beta}\|^{2},$

where
$\boldsymbol{\beta}=(\beta_{1},\beta_{2},\ldots,\beta_{p})^{\scriptscriptstyle\textrm{T}}$,
$\boldsymbol{y}=(y_{1},y_{2},\ldots,y_{n})^{\scriptscriptstyle\textrm{T}}$ and
$\boldsymbol{B}$ is an $n\times p$ design matrix with element
$B_{ij}=\mathcal{B}_{j}(x_{i})$. This gives least squares solution
$\boldsymbol{\hat{\beta}}=(\boldsymbol{B^{\scriptscriptstyle\textrm{T}}B})^{-1}\boldsymbol{B^{\scriptscriptstyle\textrm{T}}y}$
and fitted values $\boldsymbol{\hat{y}}=\boldsymbol{B\hat{\beta}}$.

The estimator makes it easy to fit nonlinear smooth curves to data, but it is
difficult to decide how many knots are needed and where to place them in
practice. Figure 6 shows that when $k$ is small, $f(x)$ underfits data and its
shape is very sensitive to exact knot locations. Choosing a big $k$ mitigates
such sensitivity but overfits data. There are attempts to estimate or select
optimal knots, resulting in “free-knot” splines (Jupp, 1978; Spiriti et al.,
2013) or locally adaptive regression splines (Mammen and van de Geer, 1997),
but they are computationally intensive. The smoothing spline in the next
section takes an alternative route. It bypasses the knot placement problem by
placing knots at $(x_{i})_{1}^{n}$, then penalizes the roughness of $f(x)$ to
avoid getting an interpolation spline. At it turns out, this strategy is both
effective and computationally efficient.

Figure 6: Regression cubic B-splines fit (black) with different $k$ and knot
locations (gray dashed lines) fitted to 200 noisy observations (gray dots)
from $g(x)=(x-0.4)\sin(3\pi x)+0.2x$ (red) at equidistant $x$-values on
$[0,1]$.

###  3.2 Smoothing spline

For any $\lambda>0$, $a<x_{1}$, $b>x_{n}$ and positive integer $m$, a
smoothing spline (Wahba, 1990; Wang, 2011) minimizes the penalized least
squares objective:

$\sum_{i=1}^{n}(y_{i}-f(x_{i}))^{2}+\lambda\int_{a}^{b}f^{\scriptscriptstyle(m)}(x)^{2}\textrm{d}x.$
(3)

It is a nonparametric estimator, since there is no a priori assumption on
$f(x)$’s functional form, but Reinsch (1967, 1971) proved that the solution
$f(x)$ is actually an order-$2m$ polynomial spline with interior knots
$(x_{i})_{1}^{n}$. Referring to the piecewise polynomial representation (1),
$f(x)$ has $n+2m\geq n$ free parameters. However, minimizing (3) enforces
$\int_{a}^{x_{1}}f^{\scriptscriptstyle(m)}(x)^{2}\textrm{d}x=\int_{x_{n}}^{b}f^{\scriptscriptstyle(m)}(x)^{2}\textrm{d}x=0$,
implying that $f^{\scriptscriptstyle(m)}(x)=0$ on $[a,x_{1}]$ or $[x_{n},b]$,
i.e., $f(x)$ is an order-$m$ rather than an order-$2m$ polynomial segment on
$[a,x_{1}]$ or $[x_{n},b]$. This eliminates $2m$ coefficients so that $f(x)$
has as many free parameters as data. On one hand, it is an interpolation
spline when $\lambda=0$. On the other hand, as $\lambda\to+\infty$, (3) is
dominated by $\int_{a}^{b}f^{\scriptscriptstyle(m)}(x)^{2}\textrm{d}x$, which
is minimized when $f(x)$ is an order-$m$ polynomial on the entire $[a,b]$.
Intuitively, objective (3) seeks a trade-off between $f(x)$’s closeness to
data and its smoothness. A small $\lambda$ drags $f(x)$ toward the
interpolation spline (very wiggly), while a large $\lambda$ pulls $f(x)$
toward an order-$m$ polynomial (very smooth). For this reason, $\lambda$ is
called a smoothing parameter and
$\int_{a}^{b}f^{\scriptscriptstyle(m)}(x)^{2}\textrm{d}x$ a wiggliness
penalty.

Figure 7: Linear, cubic and quintic smoothing splines (black) fitted to 200
noisy observations (gray dots) from $g(x)=(x-0.4)\sin(3\pi x)+0.2x$ (red) at
equidistant $x$-values on $[0,1]$. In (A)(B)(C), $\hat{f}(x)$ with optimal
$\lambda$ is plotted. In (D), a smaller and a bigger $\lambda$ than its
optimal value are attempted. Edf of these curves are displayed for
interpretability.

Given $\lambda$, smoothing spline can be computed efficiently in $O(n)$ time.
The main problem is how to tune $\lambda$. Cross-validation may be used, and
in fact, an explicit formula for the leave-one-out cross-validation score
exists, known as the generalized cross-validation (GCV) (Craven and Wahba,
1979):

$\textrm{GCV}(\lambda)=\frac{n\|\boldsymbol{y}-\boldsymbol{\hat{y}}\|^{2}}{(n-\textrm{edf})^{2}}.$

It is calculable from the fitted spline $\hat{y}_{i}=\hat{f}(x_{i})$ using all
data.
$\boldsymbol{\hat{y}}=(\hat{y}_{1},\hat{y}_{2},\ldots,\hat{y}_{n})^{\scriptscriptstyle\textrm{T}}$
is the vector of fitted values and the effective degree of freedom in
$\hat{f}(x)$ is

$\textrm{edf}(\hat{f})=\textrm{tr}(\boldsymbol{A}),$

where $\boldsymbol{A}$ is a $\lambda$-dependent matrix such that
$\boldsymbol{\hat{y}}=\boldsymbol{Ay}$. Edf is an interpretable measure for
the complexity of $\hat{f}(x)$. It is monotone in $\lambda$: the higher
$\lambda$ is, the smaller edf is. It grows to $n$ as $\lambda\to 0$ and
diminishes to $m$ as $\lambda\to+\infty$. The best $\lambda$, hence the best
edf that minimizes GCV score can be selected. Other criteria for choosing
$\lambda$ exists. For example, the wiggliness penalty can be interpreted as a
Gaussian prior (Silverman, 1985), hence a smoothing spline can also be
estimated as a linear mixed effect model (Pinheiro and Bates, 2000), with
$\lambda$ selected to maximize restricted log-likelihood. See Ruppert, Wand
and Carroll (2003) for details.

Figure 7 illustrates linear ($m=1$), cubic ($m=2$) and quintic ($m=3$)
smoothing splines fitted to the same data. For prediction on $[0,1]$, linear
smoothing spline is more wiggly the others. For extrapolation outside $[0,1]$,
quintic smoothing spline performs quadratic extrapolation which quickly goes
to $\infty$. Since cubic smoothing spline provides both satisfactory fit and
plausible extrapolation, it is often the favorite choice in practice.

Smoothing spline is super flexible, in the sense that its edf ranges from $m$
to $n$ when $\lambda$ varies. However, there is a question of whether this is
really necessary, given that the optimal edf of a smoothing spline is often
much smaller than $n$. For example, in Figure 7(B), the optimal edf for cubic
smoothing spline is about 15 against $n=200$, thus 185 out of 200 free
parameters are completely suppressed by penalization, which is a big waste.
Observing that edf = 35 already overfits data, the excessive representation
power for smoothing spline with edf $>$ 35 is of little use. A spline
estimator whose edf only increases to 35 when $\lambda\to 0$ would have been
powerful enough. In general, it makes sense to bound the maximum edf of
$f(x)$, because complexity of $g(x)$ in a generative model is independent of
$n$. Clearly, smoothing spline whose maximum edf is always $n$ does not
satisfy this.

Another problem of smoothing spline is that the order of spline is $2m$, which
is implicitly bound with the order of penalized derivatives. For example, it
is not possible to penalize the first derivative while getting a cubic spline.
In addition, smoothing spline of even degree simply does not exist.
Fortunately, the O-spline estimator coming next overcomes both drawbacks of
smoothing spline.

###  3.3 O-spline

Consider penalized least squares objective (3) again. This time we
deliberately set up $f(x)$ as regression B-splines, by choosing an order $d$
and placing $k<n$ interior knots $a<x_{1}<s_{1}<s_{2}<\ldots<s_{k}<x_{n}<b$
(as well as auxiliary boundary knots). O’Sullivan (1986) finds that under
basis representation (2), the wiggliness penalty is a quadratic form:

$\int_{a}^{b}f^{\scriptscriptstyle(m)}(x)^{2}\textrm{d}x=\boldsymbol{\beta^{\scriptscriptstyle\textrm{T}}S\beta},$

where $\boldsymbol{S}$ is a $p\times p$ symmetric penalty matrix with entry
$S_{uv}=\int_{a}^{b}\mathcal{B}_{u}^{\scriptscriptstyle(m)}(x)\mathcal{B}_{v}^{\scriptscriptstyle(m)}(x)\textrm{d}x$.
Then (3) becomes:

$\|\boldsymbol{y}-\boldsymbol{B\beta}\|^{2}+\lambda\boldsymbol{\beta^{\scriptscriptstyle\textrm{T}}S\beta}.$
(4)

Minimizing this objective gives a $\lambda$-dependent solution
$\boldsymbol{\hat{\beta}}=(\boldsymbol{B^{\scriptscriptstyle\textrm{T}}B}+\lambda\boldsymbol{S})^{-1}\boldsymbol{B^{\scriptscriptstyle\textrm{T}}y}$.
Fitted values are therefore
$\boldsymbol{\hat{y}}=\boldsymbol{B\hat{\beta}}=\boldsymbol{Ay}$, and

$\boldsymbol{A}=\boldsymbol{B}(\boldsymbol{B^{\scriptscriptstyle\textrm{T}}B}+\lambda\boldsymbol{S})^{-1}\boldsymbol{B^{\scriptscriptstyle\textrm{T}}}$

is the matrix that defines the edf of the fitted spline. As with smoothing
spline, the optimal $\lambda$ can be selected to minimize GCV score.

Figure 8: O-spline (black) with 30 interior knots at different locations (gray
dashed lines) fitted to 200 noisy observations (gray dots) from
$g(x)=(x-0.4)\sin(3\pi x)+0.2x$ (red) at equidistant $x$-values. $\lambda$ is
selected to minimize GCV score.

Wand and Ormerod (2008) names this estimator the O-spline after O’Sullivan.
Combining basis representation with wiggliness penalty, it is a crossbreed of
regression B-splines and smoothing spline. It overcomes weaknesses of
regression B-splines. As Figure 8 shows, it does not overfit data when $k$ is
big, and exact knot locations are less relevant to the resulting fit. It also
has the merits that smoothing spline does not have. (1) Its maximum edf is $p$
instead of $n$. As $\lambda\to 0$, it approaches a regression spline instead
of an interpolation spline. (2) It allows us to mix and match the order of
spline and penalized derivatives (as long as the latter is smaller). Comparing
Figure 8 to 7, it is also clear that O-spline is as good as smoothing spline,
even if it positions much fewer interior knots.

Interior knots of O-spline are usually placed at equal quantiles of
$(x_{i})_{1}^{n}$. The reason is often not explained, but here is an
intuition. Suppose that $n$ is too large to fit a smoothing spline fast
enough, then we may draw $(\tilde{x}_{i},\tilde{y}_{i})_{1}^{k}$ from
$(x_{i},y_{i})_{1}^{n}$ and fit one to these subsamples instead. It makes
sense to pick $(\tilde{x}_{i})_{1}^{k}$ at equal quantiles of
$(x_{i})_{1}^{n}$, as they are closest to $(x_{i})_{1}^{n}$ in distribution.
Consequently, it is natural to use such
$(\tilde{x}_{i},\tilde{y}_{i})_{1}^{k}$ as knots for an O-spline. With this
connection, O-spline may be considered equivalent to smoothing spline. In
fact, R function smooth.spline fits an O-spline rather than a smoothing spline
when $n>50$, unless we set all.knots = TRUE.

O-spline is an appealing estimator, but it used to have a nontrivial setup
cost in computing $\boldsymbol{S}$ accurately for $m>4$, which is not lifted
only till recently (Wood, 2017a). For this historical reason, standard
P-spline was proposed in 1990s has been a popular alternative ever since.

###  3.4 Standard P-spline

Figure 9: Left column: standard P-spline (black) with 30 uniform interior
knots (gray dashed lines) fitted to 200 noisy observations (gray dots) from
$g(x)=(x-0.4)\sin(3\pi x)+0.2x$ (red) at equidistant $x$-values. Right column:
estimated B-spline coefficients. (A1)(B1): too big a $\lambda$; (A2)(B2):
optimal $\lambda$ minimizing GCV score; (A3)(B3): too small a $\lambda$.

Standard P-spline (Eilers and Marx, 1996) also features the combination of
B-splines and wiggliness penalty, and is a close relative of O-spline. But
unlike O-spline which penalizes
$\int_{a}^{b}f^{\scriptscriptstyle(m)}(x)^{2}\textrm{d}x$, it penalizes the
sum of squared order-$m$ differences of nearby B-spline coefficients. For
example, for $m=1,2,3$, it penalizes

$\displaystyle\sum_{j=1}^{p-1}(\beta_{j+1}-\beta_{j})^{2},\kern
14.22636pt\sum_{j=1}^{p-2}(\beta_{j+2}-2\beta_{j+1}+\beta_{j})^{2},$
$\displaystyle\sum_{j=1}^{p-3}(\beta_{j+3}-3\beta_{j+2}+3\beta_{j+1}-\beta_{j})^{3},$

respectively. They can be written as
$\|\boldsymbol{D}_{p}^{\scriptscriptstyle(1)}\boldsymbol{\beta}\|^{2}$,
$\|\boldsymbol{D}_{p}^{\scriptscriptstyle(2)}\boldsymbol{\beta}\|^{2}$, and
$\|\boldsymbol{D}_{p}^{\scriptscriptstyle(3)}\boldsymbol{\beta}\|^{2}$ for
briefness, where

$\boldsymbol{D}_{p}^{\scriptscriptstyle(m)}=\boldsymbol{\Delta}_{p-m}\cdots\boldsymbol{\Delta}_{p-2}\boldsymbol{\Delta}_{p-1}$
(5)

is the order-$m$ standard difference matrix of dimension $(p-m)\times p$, and

$\boldsymbol{\Delta}_{r}=\begin{bmatrix}-1&1\\\ &-1&1\\\ &&\ddots&\ddots\\\
&&&-1&1\end{bmatrix}_{r\times(r+1)}$ (6)

is a fundamental difference matrix. The following are some concrete examples
of $\boldsymbol{D}_{p}^{\scriptscriptstyle(m)}$.

$\displaystyle\boldsymbol{D}_{5}^{\scriptscriptstyle(1)}=\boldsymbol{\Delta}_{4}=\begin{bmatrix}-1&1\\\
&-1&1\\\ &&-1&1\\\ &&&-1&1\end{bmatrix},$
$\displaystyle\boldsymbol{D}_{6}^{\scriptscriptstyle(2)}=\boldsymbol{\Delta}_{4}\boldsymbol{\Delta}_{5}=\begin{bmatrix}1&-2&1\\\
&1&-2&1\\\ &&1&-2&1\\\ &&&1&-2&1\end{bmatrix},$
$\displaystyle\boldsymbol{D}_{7}^{\scriptscriptstyle(3)}=\boldsymbol{\Delta}_{4}\boldsymbol{\Delta}_{5}\boldsymbol{\Delta}_{6}=\begin{bmatrix}-1&3&-3&1\\\
&-1&3&-3&1\\\ &&-1&3&-3&1\\\ &&&-1&3&-3&1\end{bmatrix}.$

With the aid of $\boldsymbol{D}_{p}^{\scriptscriptstyle(m)}$, the penalized
least squares objective for standard P-spline is

$\|\boldsymbol{y}-\boldsymbol{B\beta}\|^{2}+\lambda\|\boldsymbol{D}_{p}^{\scriptscriptstyle(m)}\boldsymbol{\beta}\|^{2},$
(7)

or equivalently,

$\|\boldsymbol{y}-\boldsymbol{B\beta}\|^{2}+\lambda\boldsymbol{\beta^{\scriptscriptstyle\textrm{T}}}\boldsymbol{D}_{p}^{\scriptscriptstyle(m)\scriptscriptstyle\textrm{T}}\boldsymbol{D}_{p}^{\scriptscriptstyle(m)}\boldsymbol{\beta}.$

Comparing this with O-spline’s objective (4), we see that standard P-spline
uses
$\boldsymbol{D}_{p}^{\scriptscriptstyle(m)\scriptscriptstyle\textrm{T}}\boldsymbol{D}_{p}^{\scriptscriptstyle(m)}$
instead of $\boldsymbol{S}$ for the penalty matrix.

Equidistant knots are mandatory for standard P-splines. The reason is that
uniform B-splines have identical shapes and local support (Figure 1), so that
smoothness of their coefficients translates to that of their linear
combination. This is well illustrated in Figure 9, where we spot a direct
correspondence between the pattern in $(\hat{\beta}_{j})_{1}^{p}$ and
$\hat{f}(x)$. So if $(\beta_{j})_{1}^{p}$ is the “skeleton” of $f(x)$, then
$(\mathcal{B}_{j}(x))_{1}^{p}$ is its “flesh”. Since the order-$m$ difference
of $(\beta_{j})_{1}^{p}$ reasonably measures the roughness of
$(\beta_{j})_{1}^{p}$, the difference penalty is a sensible wiggliness penalty
for $f(x)$. But intuitively, this logic chain will be broken once B-splines
become non-uniform. Indeed, we have demonstrated in Figure 3 the failure of
standard difference penalty for non-uniform B-splines. We will provide a
mathematical reason for this later in section 4.1. Now we are to conclude our
review by describing automatic knot placement for uniform and non-uniform
B-splines.

###  3.5 Automatic knot placement

Knot placement starts with specification of $[a,b]$, where $f(x)$ is defined
and penalized. Domain boundaries $a$ and $b$ may be simply chosen as $x_{1}$
and $x_{n}$, but choosing $a<x_{1}$ and $b>x_{n}$ is also legitimate, giving a
wider domain than the range of $x$-values. This is necessary when we want to
make prediction outside $[x_{1},x_{n}]$. A wider domain is also preferred if
we wish to reduce estimated variance near $x_{1}$ and $x_{n}$.

For the $k$ interior knots, we have already seen quantile knots for O-spline
and equidistant knots for standard P-spline. We now highlight these two types
of knots as they facilitate automatic placement of uniform or non-uniform
knots given $k$.

  * •

uniform knots. $(k+2)$ equidistant domain knots are placed on $[a,b]$, with
step size $(b-a)/(k+1)$. Auxiliary knots are then automatically determined.
Such knot placement method spawns uniform B-splines, suitable for standard
P-spline.

  * •

non-uniform knots with clamped boundary. $(k+2)$ domain knots are placed at
equal quantiles of $\min(a,x_{1}),x_{2},\ldots,x_{n-1},\max(x_{n},b)$ and
auxiliary boundary knots are duplicated at $a$ and $b$. Such knot placement
method spawns non-uniform and clamped B-splines, suitable for O-spline and our
general P-spline as well.

Note that in the second case, if $(x_{i})_{1}^{n}$ are equally spaced, so are
quantile domain knots, yet the full knot sequence is not uniform due to
coalesced boundary knots.

##  4 The new general P-spline

In this section, we formally introduce general P-spline. At its heart, we
derive a general difference matrix
$\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}$ to replace the standard
one $\boldsymbol{D}_{p}^{\scriptscriptstyle(m)}$ tailored for standard
P-spline. We first explain why $\boldsymbol{D}_{p}^{\scriptscriptstyle(m)}$
fails for non-uniform B-splines in section 4.1. The mathematical justification
also lays down a principle to validate our
$\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}$. Derivation of
$\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}$ takes two steps. In
section 4.2 we first derive B-spline representation for
$f^{\scriptscriptstyle(m)}(x)$, then in section 4.3 we arrange the result to
work out $\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}$. Note that
although $\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}$ is called a
difference matrix, it actually has direct control over
$f^{\scriptscriptstyle(m)}(x)$, underlying its link to the derivative penalty
$\int_{a}^{b}f^{\scriptscriptstyle(m)}(x)^{2}\textrm{d}x$. In section 4.4 we
discuss relations between two types of penalties and establish a sandwich
formula for them. In section 4.5 we use such formula to obtain a sparse “root”
of the derivative penalty, making it straightforward to adapt the $l^{2}$
penalty to $l^{1}$ regularization. In section 4.6 we exemplify construction of
$\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}$ through a small
example.

###  4.1 Why standard difference penalty goes wrong with non-uniform B-splines

To understand why standard difference penalty (5) goes wrong with non-uniform
B-splines, let’s investigate the limiting behaviour of $f(x)$ as
$\lambda\to+\infty$. For O-spline (4), $\lambda=+\infty$ enforces
$f^{\scriptscriptstyle(m)}(x)=0$, i.e., $f(x)$ is an order-$m$ polynomial. For
standard P-spline (7), $\lambda=+\infty$ enforces
$\boldsymbol{D}_{p}^{\scriptscriptstyle(m)}\boldsymbol{\beta}=\boldsymbol{0}$,
i.e.,
$\boldsymbol{\beta}\in\textrm{null}(\boldsymbol{D}_{p}^{\scriptscriptstyle(m)})$.
Let $\boldsymbol{H}$ be a $p\times m$ null space basis matrix, then
$\boldsymbol{\beta}=\boldsymbol{H\alpha}\in\textrm{null}(\boldsymbol{D}_{p}^{\scriptscriptstyle(m)}),\
\forall\boldsymbol{\alpha}\in\mathbb{R}^{m}$. Therefore, when
$\lambda=+\infty$, (7) becomes an ordinary least squares objective:

$\|\boldsymbol{y}-\boldsymbol{BH\alpha}\|^{2}=\|\boldsymbol{y}-\boldsymbol{X\alpha}\|^{2}.$

Denote its solution by $\boldsymbol{\hat{\alpha}}$, then
$\boldsymbol{\hat{y}}=\boldsymbol{X\hat{\alpha}}$ is the vector of fitted
values. The question is: does $\boldsymbol{\hat{y}}$ also represent an
order-$m$ polynomial? The answer is yes for uniform B-splines but no for non-
uniform B-splines. To recover the correct limiting behaviour of P-spline in
case of unevenly spaced knots, we need to replace
$\boldsymbol{D}_{p}^{\scriptscriptstyle(m)}$ by a proper matrix with the right
null space. This motivates the derivation of our general difference penalty
matrix $\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}$.

Figure 10: Limiting behaviour of several cubic P-spline variants of 2nd order
difference penalty. 400 noisy data with unevenly spaced $x$-values on $[-3,3]$
are generated from $g(x)=x$ so that $\lambda=+\infty$ is the optimal smoothing
parameter. Columns (A) and (B) experiment $k$ = 10 and 50 interior knots,
respectively. Gray vertical dashed lines mark locations of all knots.

To visualize the null space and the P-spline fit at $\lambda=+\infty$, let’s
take $m=2$ and do a simple experiment. We simulate 400 noisy observations with
unevenly spaced $x$-values on $[-3,3]$ from $g(x)=x$ then fit the following
cubic P-splines: (1) standard P-spline (uniform B-splines with standard
difference penalty), (2) naive P-spline (non-uniform B-splines with standard
difference penalty) and (3) general P-spline (non-uniform B-splines with
general difference penalty). The results are shown in Figure 10, and to see if
the number of interior knots matters, $k=10$ (column A) and $k=50$ (column B),
are attempted. In each graph, the two black dashed lines are columns of
$\boldsymbol{X}$ and the black solid line is their linear combination, i.e.,
the fitted spline $\boldsymbol{\hat{y}}=\boldsymbol{X\hat{\alpha}}$. (A1) and
(B1) show that for standard P-spline, $\boldsymbol{\hat{y}}$ is a straight
line regardless of $k$. But (A2) and (B2) show that this is not true for non-
uniform B-splines with standard difference penalty. Fortunately, (A3) and (B3)
verify that our general P-spline corrects this fault. We will now derive the
general penalty matrix $\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}$
that characterizes our general P-spline.

###  4.2 B-spline representation for $f^{\scriptscriptstyle(m)}(x)$

In this and the next sections, we need to index a B-spline and its coefficient
by its order as well. To begin with, we rewrite basis representation (2) as

$f(x)=\sum_{j=1}^{p}\mathcal{B}_{j,d}(x)\beta_{j,d}.$

It follows that

$f^{\prime}(x)=\sum_{j=1}^{p}\mathcal{B}_{j,d}^{\prime}(x)\beta_{j,d}.$

de Boor (2001) shows that

$\mathcal{B}_{j,d}^{\prime}(x)=(d-1)\left(\frac{\mathcal{B}_{j,d-1}(x)}{\gamma_{j,d-1}}-\frac{\mathcal{B}_{j+1,d-1}(x)}{\gamma_{j+1,d-1}}\right).$
(8)

Therefore,

$\frac{f^{\prime}(x)}{d-1}=\sum_{i=1}^{p}\frac{\mathcal{B}_{j,d-1}(x)}{\gamma_{j,d-1}}\beta_{j,d}-\sum_{i=1}^{p}\frac{\mathcal{B}_{j+1,d-1}(x)}{\gamma_{j+1,d-1}}\beta_{j,d}.$

In fact, the first term $\mathcal{B}_{1,d-1}(x)$ in the first summation and
the last term $\mathcal{B}_{p+1,d-1}(x)$ in the second summation are zero,
because their supports are outside the domain. We may revisit Figure 5 for a
hint, paying attention to those low-order B-splines in gray. As a result,

$\begin{split}&\frac{f^{\prime}(x)}{d-1}=\sum_{j=2}^{p}\frac{\mathcal{B}_{j,d-1}(x)}{\gamma_{j,d-1}}\beta_{j,d}-\sum_{j=1}^{p-1}\frac{\mathcal{B}_{j+1,d-1}(x)}{\gamma_{j+1,d-1}}\beta_{j,d}\\\
=&\sum_{j^{\star}=1}^{p-1}\frac{\mathcal{B}_{j^{\star}+1,d-1}(x)}{\gamma_{j^{\star}+1,d-1}}\beta_{j^{\star}+1,d}-\sum_{j=1}^{p-1}\frac{\mathcal{B}_{j+1,d-1}(x)}{\gamma_{j+1,d-1}}\beta_{j,d}\\\
=&\sum_{j=1}^{p-1}\mathcal{B}_{j+1,d-1}(x)\frac{\beta_{j+1,d}-\beta_{j,d}}{\gamma_{j+1,d-1}},\end{split}$

where a change of variable $j=j^{\star}+1$ is applied to the first summation
from line 1 to 2. The result may alternatively be expressed as

$\displaystyle\beta_{j+1,d-1}=\frac{\beta_{j+1,d}-\beta_{j,d}}{\gamma_{j+1,d-1}/(d-1)},$
(9a) $\displaystyle
f^{\prime}(x)=\sum_{j=1}^{p-1}\mathcal{B}_{j+1,d-1}(x)\beta_{j+1,d-1}.$ (9b)

Further differentiating (9b) gives

$f^{\prime\prime}(x)=\sum_{j=1}^{p-1}\mathcal{B}_{j+1,d-1}^{\prime}(x)\beta_{j+1,d-1}.$

Substituting $j$ with $j+1$ and $d$ with $d-1$ in (8) gives

$\mathcal{B}_{j+1,d-1}^{\prime}(x)=(d-2)\left(\frac{\mathcal{B}_{j+1,d-2}(x)}{\gamma_{j+1,d-2}}-\frac{\mathcal{B}_{j+2,d-2}(x)}{\gamma_{j+2,d-2}}\right).$

By the similar argument above, we can show

$\frac{f^{\prime\prime}(x)}{d-2}=\sum_{j=1}^{p-2}\mathcal{B}_{j+2,d-2}(x)\frac{\beta_{j+2,d-1}-\beta_{j+1,d-1}}{\gamma_{j+2,d-2}},$

which can be rearranged as:

$\displaystyle\beta_{j+2,d-2}=\frac{\beta_{j+2,d-1}-\beta_{j+1,d-1}}{\gamma_{j+2,d-2}/(d-2)},$
(10a) $\displaystyle
f^{\prime\prime}(x)=\sum_{j=1}^{p-2}\mathcal{B}_{j+2,d-2}(x)\beta_{j+2,d-2}.$
(10b)

Applying induction to (9a)(9b)(10a)(10b), we obtain the result for
$f^{\scriptscriptstyle(m)}(x)$ ($1\leq m\leq d-1$):

$\displaystyle\beta_{j+m,d-m}=\frac{\beta_{j+m,d-m+1}-\beta_{j+m-1,d-m+1}}{\gamma_{j+m,d-m}/(d-m)},$
(11a) $\displaystyle
f^{\scriptscriptstyle(m)}(x)=\sum_{j=1}^{p-m}\mathcal{B}_{j+m,d-m}(x)\beta_{j+m,d-m}.$
(11b)

These results are not new. de Boor (2001) gives similar formulae on page 117
without explicit summation limits. However, de Boor never thought of deriving
a difference penalty from these results. His interest was the evaluation of
$f^{\scriptscriptstyle(m)}(x)$ from its B-spline form (FORTRAN program
BVALUE), or the evaluation of
$(\mathcal{B}_{j,d}^{\scriptscriptstyle(m)}(x))_{1}^{p}$ (BSPLVB and BSPLVD).
We also comment that results for arbitrary knot sequence are rarely seen in
statistical literature. It is more common to see special results for uniform
knots.

###  4.3 General difference matrix and general P-spline

We now work out our general different matrix
$\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}$ from (11a)(11b). To get
rid of distracting subscripts in these equations, let’s define

$\displaystyle\boldsymbol{\bar{\mathcal{B}}}_{d-m}(x)=\begin{bmatrix}\mathcal{B}_{1+m,d-m}(x)\\\
\mathcal{B}_{2+m,d-m}(x)\\\ \vdots\\\ \mathcal{B}_{p,d-m}(x)\end{bmatrix},\
\boldsymbol{\bar{\beta}}_{d-m}=\begin{bmatrix}\beta_{1+m,d-m}\\\
\beta_{2+m,d-m}\\\ \vdots\\\ \beta_{p,d-m}\end{bmatrix},$
$\displaystyle\boldsymbol{R}_{d-m}=\begin{bmatrix}\frac{\gamma_{1+m,d-m}}{d-m}\\\
&\frac{\gamma_{2+m,d-m}}{d-m}\\\ &&\ddots\\\
&&&\frac{\gamma_{p,d-m}}{d-m}\end{bmatrix}.$

We may revisit Figure 5 to understand these quantities.
$\boldsymbol{\bar{\mathcal{B}}}_{d-m}(x)$ collects all order-$(d-m)$ B-splines
that are nonzero on the domain, i.e., all B-splines that are not gray in the
figure. Coefficients of these B-splines are in
$\boldsymbol{\bar{\beta}}_{d-m}$. The support of these B-splines is indicated
by $(t_{j})_{1+m}^{K-m}$, i.e., knots that are not gray in the figure. The
lag-$(d-m)$ difference of this knot sequence is $(\gamma_{j,d-m})_{m+1}^{p}$.
Dividing it by $(d-m)$ gives the average distance between adjacent knots,
which makes the diagonal elements of the weight matrix $\boldsymbol{R}_{d-m}$.

Using these quantities, (11a)(11b) can be written in a neat matrix-vector
form:

$\displaystyle\boldsymbol{\bar{\beta}}_{d-m}=\boldsymbol{\bar{\Delta}}_{p-m}\boldsymbol{\bar{\beta}}_{d-m+1},$
(12a) $\displaystyle
f^{\scriptscriptstyle(m)}(x)=\boldsymbol{\bar{\mathcal{B}}}_{d-m}(x)^{\scriptscriptstyle\textrm{T}}\boldsymbol{\bar{\beta}}_{d-m},$
(12b)

where
$\boldsymbol{\bar{\Delta}}_{p-m}=\boldsymbol{R}_{d-m}^{-1}\boldsymbol{\Delta}_{p-m}$
is a weighted version of the fundamental difference matrix
$\boldsymbol{\Delta}_{p-m}$ as defined in (6). Expanding (12a) gives

$\begin{split}\boldsymbol{\bar{\beta}}_{d-m}&=\boldsymbol{\bar{\Delta}}_{p-m}\boldsymbol{\bar{\beta}}_{d-m+1}\\\
&=\boldsymbol{\bar{\Delta}}_{p-m}\boldsymbol{\bar{\Delta}}_{p-m+1}\boldsymbol{\bar{\beta}}_{d-m+2}\\\
&=\cdots=\boldsymbol{\bar{\Delta}}_{p-m}\cdots\boldsymbol{\bar{\Delta}}_{p-2}\boldsymbol{\bar{\Delta}}_{p-1}\boldsymbol{\bar{\beta}}_{d},\end{split}$

from which we observe our general difference matrix:

$\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}=\boldsymbol{\bar{\Delta}}_{p-m}\cdots\boldsymbol{\bar{\Delta}}_{p-2}\boldsymbol{\bar{\Delta}}_{p-1}.$
(13)

Compared with the standard difference matrix
$\boldsymbol{D}_{p}^{\scriptscriptstyle(m)}$ as defined in (5),
$\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}$ takes knot spacing into
consideration by the weighting process. For uniform knots,
$\boldsymbol{R}_{d-m}$ have identical diagonal entries so that
$\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}$ is proportional to
$\boldsymbol{D}_{p}^{\scriptscriptstyle(m)}$. Therefore, our general
difference matrix contains the standard one as a special case.

Using this new difference matrix, we define our general P-spline as the
minimizer of the following penalized least squares objective:

$\sum_{i=1}^{n}\big{(}y_{i}-\boldsymbol{\bar{\mathcal{B}}}_{d}(x_{i})^{\scriptscriptstyle\textrm{T}}\boldsymbol{\bar{\beta}}_{d}\big{)}^{2}+\lambda\|\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}\boldsymbol{\bar{\beta}}_{d}\|^{2}.$

Going back to the simpler notations as used in section 3,
$\boldsymbol{\bar{\beta}}_{d}$ coincides with $\boldsymbol{\beta}$ and
evaluating $\boldsymbol{\bar{\mathcal{B}}}_{d}(x)$ at $(x_{i})_{1}^{n}$ gives
design matrix $\boldsymbol{B}$. Therefore, the objective can be expressed as

$\|\boldsymbol{y}-\boldsymbol{B\beta}\|^{2}+\lambda\|\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}\boldsymbol{\beta}\|^{2}.$
(14)

Comparing with objective (7) of the standard P-spline, we have only made a
substitution of difference matrix. But the restriction that $\boldsymbol{B}$
must be set up with uniform B-splines is now removed, as
$\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}$ can properly deal with
unevenly spaced knots.

Our general P-spline has the correct limiting behaviour at $\lambda=+\infty$,
as is demonstrated in Figure 10. A formal proof is also straightforward. We
first rewrite (12a) as (note: the simpler notation $\boldsymbol{\beta}$ is
used in place of $\boldsymbol{\bar{\beta}}_{d}$):

$\boldsymbol{\bar{\beta}}_{d-m}=\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}\boldsymbol{\beta},$
(15)

Together with (12b), it is easy to establish the following equivalence:

$\boldsymbol{\beta}\in\textrm{null}(\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)})\kern
5.69054pt\Leftrightarrow\kern
5.69054pt\boldsymbol{\bar{\beta}}_{d-m}=\boldsymbol{0}\kern
5.69054pt\Leftrightarrow\kern 5.69054ptf^{\scriptscriptstyle(m)}(x)=0.$

As a result,

$\lambda=+\infty\kern 5.69054pt\Rightarrow\kern
5.69054pt\boldsymbol{\beta}\in\textrm{null}(\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)})\kern
5.69054pt\Rightarrow\kern 5.69054ptf^{\scriptscriptstyle(m)}(x)=0.$

###  4.4 The link to derivative penalty and O-spline

In this section, we will compare our general difference penalty
$\|\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}\boldsymbol{\beta}\|^{2}$
with derivative penalty
$\int_{a}^{b}f^{\scriptscriptstyle(m)}(x)^{2}\textrm{d}x$ from two aspects.

From an interpretation point of view, they are squared $l^{2}$ norms defined
on different vector spaces. Our general difference penalty
$\|\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}\boldsymbol{\beta}\|^{2}=\|\boldsymbol{\bar{\beta}}_{d-m}\|^{2}$
is the squared $l^{2}$ norm of B-spline coefficients in
$f^{\scriptscriptstyle(m)}(x)$, defined on a Euclidean space. By contrast, the
derivative penalty $\int_{a}^{b}f^{\scriptscriptstyle(m)}(x)^{2}\textrm{d}x$
is the squared $l^{2}$ norm of $f^{\scriptscriptstyle(m)}(x)$, defined on a
function space.

From a computational point of view, they are quadratic forms of different
penalty matrices. Plugging (15) into (12b), we get

$f^{\scriptscriptstyle(m)}(x)=\boldsymbol{\bar{\mathcal{B}}}_{d-m}(x)^{\scriptscriptstyle\textrm{T}}\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}\boldsymbol{\beta},$

so that we can calculate the derivative penalty as

$\begin{split}\int_{a}^{b}f^{\scriptscriptstyle(m)}(x)^{2}\textrm{d}x&=\int_{a}^{b}f^{\scriptscriptstyle(m)}(x)^{\scriptscriptstyle\textrm{T}}f^{\scriptscriptstyle(m)}(x)\textrm{d}x\\\
&=\boldsymbol{\beta}^{\scriptscriptstyle\textrm{T}}\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)\scriptscriptstyle\textrm{T}}\boldsymbol{\bar{S}}_{d-m}\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}\boldsymbol{\beta}.\end{split}$

$\boldsymbol{\bar{S}}_{d-m}$ is a $(p-m)\times(p-m)$ symmetric matrix with
element

$(\boldsymbol{\bar{S}}_{d-m})_{uv}=\int_{a}^{b}\mathcal{B}_{u,d-m}(x)\mathcal{B}_{v,d-m}(x)\textrm{d}x$
(16)

The penalty matrix in this quadratic form is

$\boldsymbol{S}_{d}=\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)\scriptscriptstyle\textrm{T}}\boldsymbol{\bar{S}}_{d-m}\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}.$
(17)

By contrast, our general difference penalty is

$\|\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}\boldsymbol{\beta}\|^{2}=\boldsymbol{\beta}^{\scriptscriptstyle\textrm{T}}\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)\scriptscriptstyle\textrm{T}}\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}\boldsymbol{\beta}.$

The penalty matrix in this quadratic form is
$\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)\scriptscriptstyle\textrm{T}}\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}$.
The discrepancy between two penalty matrices lies in the matrix that is
sandwiched by
$\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)\scriptscriptstyle\textrm{T}}$
and $\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}$, which is
$\boldsymbol{\bar{S}}_{d-m}$ for the derivative penalty and an identity matrix
for our general difference penalty.

So how far is $\boldsymbol{\bar{S}}_{d-m}$ from an identity matrix? Recall
that B-splines have local support, so $\boldsymbol{\bar{S}}_{d-m}$ is a band
matrix with $(d-m)$ sub-diagonals and supper-diagonals (due to symmetry). Here
are some examples to illustrate the structure of such band matrices, where
elements on each diagonal are denoted by the same symbol (which does not imply
the same value):

$\boldsymbol{\bar{S}}_{1}=\begin{bmatrix}\times&\circ\\\ \circ&\times&\circ\\\
&\ddots&\ddots&\ddots\\\ &&\circ&\times&\circ\\\
&&&\circ&\times\end{bmatrix},\
\boldsymbol{\bar{S}}_{2}=\begin{bmatrix}\times&\circ&\square\\\
\circ&\times&\circ&\square\\\ \square&\circ&\times&\circ&\square\\\
&\ddots&\ddots&\ddots&\ddots&\ddots\\\ &&\square&\circ&\times&\circ&\square\\\
&&&\square&\circ&\times&\circ\\\ &&&&\square&\circ&\times\end{bmatrix}.$

The closer $m$ is to $d$, the closer $\boldsymbol{\bar{S}}_{d-m}$ is to
diagonal. But it is not possible to have $m=d$ because $\boldsymbol{R}_{0}$
(i.e., lag-0 difference) and
$\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(d)}$ do not exist.
Therefore, our general difference penalty is always a simplification of
derivative penalty. Nevertheless, both penalties have the same null space, so
as $\lambda\to+\infty$, both O-spline (4) and our general P-spline (14) tend
to an order-$m$ polynomial ($f^{\scriptscriptstyle(m)}(x)=0$). To the opposite
limit, as $\lambda\to 0$, both estimators tend to the same regression
B-splines. Thus, both wiggliness penalties have the same role and the two
estimators can be seen equivalent.

###  4.5 Sparse “root” of derivative penalty

To our knowledge, we are the first to get the sandwich formula (17) that gives
an explicit link between derivative and difference penalties. It offers a new
way to compute the derivative penalty. Previously, $\boldsymbol{S}_{d}$ is
computed by its definition, i.e.,
$(\boldsymbol{S}_{d})_{uv}=\int_{a}^{b}\mathcal{B}_{u,d}^{\scriptscriptstyle(m)}(x)\mathcal{B}_{v,d}^{\scriptscriptstyle(m)}(x)\textrm{d}x$,
using the algorithm of Wood (2017a). Now we can first compute the general
difference matrix $\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}$ by
its definition (13), then compute $\boldsymbol{\bar{S}}_{d-m}$ from (16), and
finally assemble them for $\boldsymbol{S}_{d}$ using (17). Computation of
$\boldsymbol{\bar{S}}_{d-m}$ can also make use of Wood’s algorithm, but now
there is no need to evaluate derivatives of B-splines.

The benefit to break down the computation of $\boldsymbol{S}_{d}$ is that it
yields a sparse “root” of this matrix, that is, a sparse matrix
$\boldsymbol{K}$ such that $\boldsymbol{K^{\prime}K}=\boldsymbol{S}_{d}$. It
appears that $\boldsymbol{K}$ is the Cholesky factor of $\boldsymbol{S_{d}}$,
but $\boldsymbol{S}_{d}$ is only positive semi-definite (with $m$ zero
eigenvalues related to the order-$m$ polynomial in its null space) so Cholesky
factorization is not doable. A modified algorithm with pivoting works, but
pivoting destroys sparsity in $\boldsymbol{S_{d}}$ (a symmetric band matrix
with $d$ sub-diagonals and super-diagonals), making $\boldsymbol{K}$ a dense
matrix. Alternatively, we may perform a symmetric eigendecomposition
$\boldsymbol{S}_{d}=\boldsymbol{E\Lambda E}^{\scriptscriptstyle\textrm{T}}$.
When eigenvalues are ordered in decreasing size, the last $m$ eigenvalues in
$\boldsymbol{\Lambda}$ are exactly zeros and can be discarded. Let
$\boldsymbol{\bar{E}}$ be the first $(p-m)$ columns of $\boldsymbol{E}$ and
$\boldsymbol{\bar{\Lambda}}$ be the corresponding block in
$\boldsymbol{\Lambda}$, there is
$\boldsymbol{S}_{d}=\boldsymbol{\bar{E}\bar{\Lambda}\bar{E}}^{\scriptscriptstyle\textrm{T}}$
with
$\boldsymbol{K}=\boldsymbol{\bar{\Lambda}}^{1/2}\boldsymbol{\bar{E}}^{\scriptscriptstyle\textrm{T}}$.
However, eigendecomposition also destroys sparsity, making $\boldsymbol{K}$ a
dense matrix.

Thanks to the sandwich formula (17), a sparse $\boldsymbol{K}$ can now be
found. Note that $\boldsymbol{\bar{S}}$ is positive definite so that its
Cholesky factorization
$\boldsymbol{\bar{S}}_{d-m}=\boldsymbol{U}_{d-m}^{\scriptscriptstyle\textrm{T}}\boldsymbol{U}_{d-m}$
is doable. Sparsity in $\boldsymbol{\bar{S}}$ is preserved, and
$\boldsymbol{U}_{d-m}$ is a $(p-m)\times(p-m)$ band matrix with $(d-m)$ super-
diagonals. For example, $\boldsymbol{U}_{1}$ and $\boldsymbol{U}_{2}$ have the
following structure:

$\boldsymbol{U}_{1}=\begin{bmatrix}\times&\circ\\\ &\times&\circ\\\
&&\ddots&\ddots\\\ &&&\times&\circ\\\ &&&&\times\end{bmatrix},\
\boldsymbol{U}_{2}=\begin{bmatrix}\times&\circ&\square\\\
&\times&\circ&\square\\\ &&\times&\circ&\square\\\ &&&\ddots&\ddots&\ddots\\\
&&&&\times&\circ&\square\\\ &&&&&\times&\circ\\\ &&&&&&\times\end{bmatrix}.$

Furthermore, $\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}$ is a
$(p-m)\times p$ band matrix with $m$ super-diagonals, for example:

$\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(1)}=\begin{bmatrix}\times&\circ\\\
&\times&\circ\\\ &&\ddots&\ddots\\\ &&&\times&\circ\end{bmatrix},\
\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(2)}=\begin{bmatrix}\times&\circ&\square\\\
&\times&\circ&\square\\\ &&\ddots&\ddots&\ddots\\\
&&&\times&\circ&\square\end{bmatrix}.$

As a result,
$\boldsymbol{K}=\boldsymbol{U}_{d-m}\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(m)}$
is a $(p-m)\times p$ band matrix with $(d-m)+m=d$ super-diagonals. To
illustrate the structure of $\boldsymbol{K}$, it helps to partition it into
two row chunks: one is $(p-d)\times p$, the other is $(d-m)\times p$. See
below for an example with $d=3$ and $m=1$:

$\left[\begin{array}[]{ccccccc}\times&\circ&\square&\triangle\\\
&\times&\circ&\square&\triangle\\\ &&\ddots&\ddots&\ddots&\ddots\\\
&&&\times&\circ&\square&\triangle\\\ \hline\cr&&&&\times&\circ&\square\\\
&&&&&\times&\circ\end{array}\right].$

Note that $\boldsymbol{K}$ has full row rank, equal to the rank of
$\boldsymbol{S}_{d}$.

So why is $\boldsymbol{K}$ interesting? Because with such a matrix, we can
write the derivative penalty as a squared $l^{2}$ norm:

$\int_{a}^{b}f^{\scriptscriptstyle(m)}(x)^{2}\textrm{d}x=\boldsymbol{\beta}^{\scriptscriptstyle\textrm{T}}\boldsymbol{S}_{d}\boldsymbol{\beta}=\|\boldsymbol{K\beta}\|^{2}.$

Then the O-spline objective (4) is evidently a generalized ridge problem:

$\|\boldsymbol{y}-\boldsymbol{B\beta}\|^{2}+\lambda\|\boldsymbol{K\beta}\|^{2}.$

It is well known that a ridge penalty (or $l^{2}$ regularization) does not
encourage sparsity as a LASSO penalty (or $l^{1}$ regularization) will do
(Tibshirani, 1996). Therefore, it may be worth attempting the $l^{1}$
adaptation of O-spline:

$\|\boldsymbol{y}-\boldsymbol{B\beta}\|^{2}+\lambda\|\boldsymbol{K\beta}\|_{1},$

which is a generalized LASSO problem (Tibshirani and Taylor, 2011; Arnold and
Tibshirani, 2016). In fact, it is a trend filtering problem (Tibshirani, 2014;
Ramdas and Tibshirani, 2016) as the variable $y$ is one-dimensional. The key
to its computational efficiency is that both $\boldsymbol{B}$ and
$\boldsymbol{K}$ should be sparse. And ideally, $\boldsymbol{K}$ should have
full row rank. Although sparsity in $\boldsymbol{B}$ is guaranteed by
B-splines’ local support, sparsity in $\boldsymbol{K}$ is not made possible
until we provide the sandwich formula (17). Further, our proposed
$\boldsymbol{K}$ also meets the requirement of full row rank.

We do not pursue $l^{1}$ regularization in this paper, but it is a promising
research direction. Trend filtering is closely related to locally adaptive
regression splines (Mammen and van de Geer, 1997). Therefore, the $l^{1}$
adaptation of O-spline will likely lead to adaptive smoothing.

###  4.6 Computing general difference matrix: an example

In general, for $p$ order-$d$ B-splines on $K=p+d$ knots $(t_{j})_{1}^{K}$,
general difference matrices from order 1 to $(d-1)$ can be calculated
iteratively:
$\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(1)}=\boldsymbol{\bar{\Delta}}_{p-1}$,
$\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(2)}=\boldsymbol{\bar{\Delta}}_{p-2}\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(1)}$,
… ,
$\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(d-1)}=\boldsymbol{\bar{\Delta}}_{p-(d-1)}\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(d-2)}$.
The key is to calculate differences of the correct set of knots at the right
lag, summarized below.

| $\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(d-1)}$ | … | $\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(2)}$ | $\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(1)}$
---|---|---|---|---
knots | $(t_{j})_{1+(d-1)}^{K-(d-1)}$ | … | $(t_{j})_{1+2}^{K-2}$ | $(t_{j})_{1+1}^{K-1}$
lag | 1 | … | $d-2$ | $d-1$

As a demonstrative example, consider clamped non-uniform cubic B-splines
($d=4$) on domain $[a,b]=[0,4]$ with $k=2$ interior knots $s_{1}=1$, $s_{2}=3$
and duplicated boundary knots. In this case, $p=k+d=6$ and $K=p+d=10$ (revisit
Figure 4 if necessary). The complete knot sequence $(t_{j})_{1}^{K}$ is:

$t_{1}$ | $t_{2}$ | $t_{3}$ | $t_{4}$ | $t_{5}$ | $t_{6}$ | $t_{7}$ | $t_{8}$ | $t_{9}$ | $t_{10}$
---|---|---|---|---|---|---|---|---|---
0 | 0 | 0 | 0 | 1 | 3 | 4 | 4 | 4 | 4

Step 1: Construct $5\times 6$ matrix
$\boldsymbol{\mathcal{D}}_{6}^{\scriptscriptstyle(1)}$. Lag-$3$ difference of
$(t_{j})_{2}^{9}$ is $t_{5}-t_{2}=1$, $t_{6}-t_{3}=3$, $t_{7}-t_{4}=4$,
$t_{8}-t_{5}=3$, $t_{9}-t_{6}=1$. We construct
$\boldsymbol{\mathcal{D}}_{6}^{\scriptscriptstyle(1)}$ by weighting standard
difference matrix $\boldsymbol{\Delta}_{5}$:

$\begin{split}\boldsymbol{\mathcal{D}}_{6}^{\scriptscriptstyle(1)}=\boldsymbol{\bar{\Delta}}_{5}&=\begin{bmatrix}\frac{1}{3}\\\
&\frac{3}{3}\\\ &&\frac{4}{3}\\\ &&&\frac{3}{3}\\\
&&&&\frac{1}{3}\end{bmatrix}^{-1}\begin{bmatrix}-1&1\\\ &-1&1\\\ &&-1&1\\\
&&&-1&1\\\ &&&&&-1&1\end{bmatrix}\\\ &=\begin{bmatrix}-3&3\\\ &-1&1\\\
&&-\frac{3}{4}&\frac{3}{4}\\\ &&&-1&1\\\ &&&&&-3&3\end{bmatrix}.\end{split}$

Step 2: Construct $4\times 6$ matrix
$\boldsymbol{\mathcal{D}}_{6}^{\scriptscriptstyle(2)}$. Lag-$2$ difference of
$(t_{j})_{3}^{8}$ is $t_{5}-t_{3}=1$, $t_{6}-t_{4}=3$, $t_{7}-t_{5}=3$,
$t_{8}-t_{6}=1$. We first weight standard difference matrix
$\boldsymbol{\Delta}_{4}$:

$\begin{split}\boldsymbol{\bar{\Delta}}_{4}&=\begin{bmatrix}\frac{1}{2}\\\
&\frac{3}{2}\\\ &&\frac{3}{2}\\\
&&&\frac{1}{2}\end{bmatrix}^{-1}\begin{bmatrix}-1&1\\\ &-1&1\\\ &&-1&1\\\
&&&-1&1\end{bmatrix}\\\ &=\begin{bmatrix}-2&2\\\ &-\frac{2}{3}&\frac{2}{3}\\\
&&-\frac{2}{3}&\frac{2}{3}\\\ &&&-2&2\end{bmatrix},\end{split}$

then calculate

$\boldsymbol{\mathcal{D}}_{6}^{\scriptscriptstyle(2)}=\boldsymbol{\bar{\Delta}}_{4}\boldsymbol{\mathcal{D}}_{6}^{\scriptscriptstyle(1)}=\begin{bmatrix}6&-8&2\\\
&\frac{2}{3}&-\frac{7}{6}&\frac{1}{2}\\\
&&\frac{1}{2}&-\frac{7}{6}&\frac{2}{3}\\\ &&&2&-8&6\end{bmatrix}.$

Step 3: Construct $3\times 6$ matrix
$\boldsymbol{\mathcal{D}}_{6}^{\scriptscriptstyle(3)}$. Lag-$1$ difference of
$(t_{j})_{4}^{7}$ is $t_{5}-t_{4}=1$, $t_{6}-t_{5}=2$, $t_{7}-t_{6}=1$. We
first weight standard difference matrix $\boldsymbol{\Delta}_{3}$:

$\boldsymbol{\bar{\Delta}}_{3}=\begin{bmatrix}1\\\ &2\\\
&&1\end{bmatrix}^{-1}\begin{bmatrix}-1&1\\\ &-1&1\\\
&&-1&1\end{bmatrix}=\begin{bmatrix}-1&1\\\ &-\frac{1}{2}&\frac{1}{2}\\\
&&-1&1\end{bmatrix},$

then calculate

$\boldsymbol{\mathcal{D}}_{6}^{\scriptscriptstyle(3)}=\boldsymbol{\bar{\Delta}}_{3}\boldsymbol{\mathcal{D}}_{6}^{\scriptscriptstyle(2)}=\begin{bmatrix}-6&\frac{26}{3}&-\frac{19}{6}&\frac{1}{2}\\\
&-\frac{1}{3}&\frac{5}{6}&-\frac{5}{6}&\frac{1}{3}\\\
&&-\frac{1}{2}&\frac{19}{6}&-\frac{26}{3}&6\end{bmatrix}.$

##  5 Simulation studies

The U-shaped curve example in Figure 3 has already demonstrated the success of
our general P-spline, but is probably too simple an example as it is unimodal
and in fact, a cubic spline with a single knot at $x=0$. In this section, we
compose more challenging $(x,y)$ scatterplot smoothing examples, fit 4 types
of penalized regression B-splines, namely (a) O-spline (non-uniform B-splines
with derivative penalty), (b) standard P-spline (uniform B-splines with
standard difference penalty), (c) naive P-spline (non-uniform B-splines with
standard difference penalty) and (d) general P-spline (non-uniform B-splines
with general difference penalty), then compare their MSE performance through
simulations. For all candidates, the number of interior knots $k$, B-spline
order $d$ and penalty order $m$ are controlled to be identical; only knot
location and penalty type differ. Knots are automatically placed using
strategies described in section 3.5, that is, we place equidistant knots for
uniform B-splines and quantile knots for non-uniform B-splines. To make two
knot sequences (and their B-splines) substantially different, all examples are
designed with unevenly spaced $x$-values.

###  5.1 Normal mixture example

Figure 11: Fitting 4 types of penalized regression B-splines to normal mixture
example under sampling scheme 1. In all graphs, red curve is the true function
$g(x)$, solid black curve is the estimated spline curve $\hat{f}(x)$ and gray
vertical dashed lines mark positions of knots. The black dashed curve in that
last graph is the probability density of $x$-values. Figure 12: Fitting 4
types of penalized regression B-splines to normal mixture example under
sampling scheme 2. This Figure has the same interpretation with Figure 11.

Consider a normal mixture: $g(x)=\frac{1}{2}f_{1}(x)+\frac{1}{2}f_{2}(x)$,
where $f_{1}(x)$ and $f_{2}(x)$ are probability density functions of
$N(-1,0.5^{2})$ and $N(1,0.8^{2})$, respectively. On $[-2,2]$, $g(x)$ has a
local minimum at $x=0$ and two local maxima at $x=\pm 1$. We apply two
sampling schemes for unevenly spaced $x$-values. One is to draw $x_{i}$ from a
normal mixture
$\frac{1}{3}N(-1,0.35^{2})+\frac{1}{3}N(1,0.35^{2})+\frac{1}{3}N(0,0.2^{2})$,
the other is to draw $x_{i}$ from
$\frac{1}{2}N(-1,0.35^{2})+\frac{1}{2}N(1,0.35^{2})$. Under scheme 1,
$x$-values are clustered around all local extrema, while under scheme 2, there
are few $x$-values near the local minimum. See Figure 11 and 12 for 1000
observations under two schemes. We now smooth those data using 4 types of
penalized regression B-splines with $k=100$, $d=4$ and $m=2$. Estimated spline
curves with their MSE performance under two sampling schemes are shown in
Figure 11 and 12, respectively. In both cases, (a)(b)(d) fit data equally well
and they all outperform (c). The result is maybe surprising under scheme 2,
since the distribution of quantile knots has a big gap near $x=0$, with only a
single knot in approximately $(-0.4,0.4)$, yet (d) general P-spline captures
the local minimum better than its competitors.

###  5.2 Random curve example

To better compare all candidates, it is a good idea to design a simulation
study that is not based on a specific signal, but a family of random curves.
To be precise, we repeat the following procedure for $l=1,2,\ldots,N$:

  1. 1.

generate a random curve $g^{\scriptscriptstyle[l]}(x)$ on $[0,1]$;

  2. 2.

evaluate
$g_{i}^{\scriptscriptstyle[l]}=g^{\scriptscriptstyle[l]}(x_{i}^{\scriptscriptstyle[l]})$
at unevenly spaced $(x_{i}^{\scriptscriptstyle[l]})_{1}^{n}$;

  3. 3.

simulate noisy observations
$y_{i}^{\scriptscriptstyle[l]}=g_{i}^{\scriptscriptstyle[l]}+\varepsilon_{i}^{\scriptscriptstyle[l]}$,
with
$\varepsilon_{i}^{\scriptscriptstyle[l]}\overset{\scriptscriptstyle\textrm{iid}}{\sim}N(g_{i}^{\scriptscriptstyle[l]},\sigma_{\scriptscriptstyle
l}^{2})$;

  4. 4.

estimate
$\hat{f}_{i}^{\scriptscriptstyle[l]}=\hat{f}^{\scriptscriptstyle[l]}(x_{i}^{\scriptscriptstyle[l]})$
from $(x_{i}^{\scriptscriptstyle[l]},y_{i}^{\scriptscriptstyle[l]})$ using 4
types of penalized regression B-splines;

  5. 5.

compute relative MSE $\delta_{\scriptscriptstyle
l}=\frac{1}{n}\sum_{i=1}^{n}(\hat{f}_{i}^{\scriptscriptstyle[l]}-g_{i}^{\scriptscriptstyle[l]})^{2}/\sigma_{\scriptscriptstyle
l}^{2}$ for each of the 4 fitted splines.

At step 1, we construct $g^{\scriptscriptstyle[l]}(x)$ as an order-$d$ spline:

$g^{\scriptscriptstyle[l]}(x)=\sum_{j=1}^{3d}\mathcal{B}_{j,d}^{\star}(x)\beta_{j,d}^{\scriptscriptstyle[l]},\kern
5.69054ptx\in[0,1],$

with uniform B-splines $\mathcal{B}_{j,d}^{\star}(x)$ on $K^{\star}=4d$
equidistant knots $(t_{j}^{\star})_{1}^{K}$ and random coefficients
$\beta_{j,d}^{\scriptscriptstyle[l]}$ sampled from uniform distribution on
$[-1,1]$. This setup gives $k^{\star}=2d$ interior knots in $(0,1)$,
$k^{\star}+2$ domain knots on $[0,1]$ and $p^{\star}=k^{\star}+d=3d$ B-splines
and coefficients. Given $d$, the knot sequence can be uniquely determined. It
is legitimate to specify a higher $k^{\star}$ for the experiment, but
unnecessary. Placing more knots inevitably increases the complexity in
$g^{\scriptscriptstyle[l]}(x)$ (like producing more extrema), hence requires
more observations (i.e., a bigger $n$) for reasonable curve estimation. Figure
13 sketches some example cubic ($d=4$) and quintic ($d=6$) splines,
complicated enough for our simulation study. It is worth mentioning that
placing non-uniform knots is not recommended here, as they produce
$g^{\scriptscriptstyle[l]}(x)$ with spatially inhomogeneous smoothness which
requires adaptive smoothing (Mammen and van de Geer, 1997; Luo and Wahba,
1997) for proper estimation. Note that “knots” in this step are not placed for
spline estimator $f^{\scriptscriptstyle[l]}(x)$, but placed for the data-
generating signal $g^{\scriptscriptstyle[l]}(x)$.

At step 2, we generate more $x$-values around $g^{\scriptscriptstyle[l]}(x)$’s
local extrema. This is reasonable from the design point of view, because local
extrema provide key information of $g^{\scriptscriptstyle[l]}(x)$. We sample
$(x_{i}^{\scriptscriptstyle[l]})_{1}^{n}$ from “tent” density function, as the
dashed line in Figure 13 illustrates. Note that local extrema and
$(x_{i}^{\scriptscriptstyle[l]})_{1}^{n}$ vary from signal to signal.

Figure 13: Three realizations (solid) of random cubic ($d=4$) and quintic
($d=6$) splines on domain [0, 1]. Locations of $2d+2$ domain knots are marked
by gray vertical dashed lines. Unevenly spaced $x$-values are sampled from
“tent” probability density function (black dashed), which peaks at local
extrema (black dots) of the spline.

At step 3, the error standard deviation is specified as
$\sigma_{\scriptscriptstyle
l}=\gamma\times\textrm{s.e.}\\{(g_{i}^{\scriptscriptstyle[l]})_{1}^{n}\\}$,
where $\textrm{s.e.}\\{(g_{i}^{\scriptscriptstyle[l]})_{1}^{n}\\}$ is the
standard error of signal values and $\gamma$ is a noise-to-signal ratio.
$\gamma$ is held fixed through $l=1,2,\ldots,N$, but
$\sigma_{\scriptscriptstyle l}$ will vary from signal to signal. This explains
why in step 5 we scale raw MSE value by $\sigma_{\scriptscriptstyle l}^{2}$.

At step 4, we choose a reasonable $k$ for estimator
$f^{\scriptscriptstyle[l]}(x)$ as follows. It is known from step 1 that
$g^{\scriptscriptstyle[l]}(x)$ comprises $p^{\star}=3d$ B-splines, i.e., its
true degree of freedom is $3d$. For estimation purpose it makes sense to use
more, say $p=10p^{\star}$ regression B-splines to reduce bias (because we
won’t be able to place knots exactly at its true locations), then let the
penalty to suppress excess amount of degree of freedom. It is then implied
that we should place $k=p-d=29d$ interior knots. Attempting a even larger $p$
is legitimate, but unnecessary as more coefficients will be suppressed by the
penalty without improving the quality of fit.

Figure 14: MSE performance of 4 types of penalized regression B-splines,
namely (a) O-spline, (b) standard P-spline, (c) naive P-spline and (d) general
P-spline, in our simulation study with $N=100$ (number of simulations),
$d=4,6$ (B-spline order), $m=1,2,\ldots,d-1$ (number of derivatives or
differences to penalize), and various choices of $n$ (sample size) and
$\gamma$ (noise-to-signal ratio). (A1)(A2)(A3)(A4) are the results for cubic
splines ($d=4$); (B1)(B2)(B3)(B4) are the results for quintic splines ($d=6$).

In summary, to start the simulation we need to specify the number of
simulations $N$, B-spline order $d$, number of observations $n$ and noise-to-
signal ratio $\gamma$. There is no need to choose the penalty order $m$ (i.e.,
the number of derivatives or differences to penalize); instead, we attempt all
$m=1,2,\ldots,d-1$. In this way, we can assess whether using derivative
penalty and difference penalty results in different MSE performance as $m$
varies. This is indeed an interesting question, because as the sandwich
formula (17) suggests, the two penalties get increasingly different in
structure as $d-m$ grows. Thus, it is natural to ask if a MSE performance gap
exists and would increase as $m$ decreases. To our knowledge, no simulation
study like this has been done in literature

We experiment with cubic ($d=4$) and quintic ($d=6$) splines. The quintic case
is of particular interest as there is a lot of room for $d-m$ to vary. We
specify $N=100$ and four different $(n,\gamma)$ pairs: (1000, 0.1), (1000,
0.5), (5000, 0.1) and (5000, 0.5). Simulation results are shown in Figure 14.
It is clear that apart from $m=1$, there is generally no concerning
performance difference between (a)(b)(d) and they all outperform (c) without
surprise. The conclusion is that our general P-spline (d) successfully
integrates non-uniform B-splines and difference penalty, inferior to neither
standard P-spline (b) nor O-spline (a). Moreover, there is no suspected
performance gap between using difference penalty and derivative penalty.

All candidates have an unpredictable performance for $m=1$, because none of
the random curves exemplified in Figure 13 is close to constant or piecewise
constant. There is therefore no justification for penalizing the first
derivative or difference. It should also be pointed out that our simulation
results do not favor any particular $m$. It is false to either claim that
$m=3$ performs better than $m=2$, or that $m=1$ is never a good idea. For
example, in section 6.1 we will use $m=1$ for the random effect spline in our
additive mixed model for the BMC data. The most appropriate choice of $m$ is
arguably data dependent, and sometimes multiple penalties of several orders
may be used to achieve desired effects (see section 8.5 of Eilers and Marx
(2021)). So to recap, the key message of our simulation study is that, for any
decent choice of $d$ and $m$, our general P-spline is a successful smoothing
method, whose theoretical property is as good as standard P-spline and
O-spline.

##  6 Real data examples

In this section, we apply standard P-spline (candidate (b) in our simulation
studies) and our general P-spline (candidate (d)) to real datasets and compare
their quality of fit. Although they have comparable MSE performance according
to our simulation studies, it turns out that one can gain an edge over the
other in practical smoothing. That is, placing uniform or non-uniform knots
can yield a more satisfactory fit for some dataset, thus appears to be more
advantageous. This is most likely to happen when $k$ (number of interior
knots) is not big enough to eliminate sensitivity in exact knot locations.
Thus, such observation should not be interpreted as one spline estimator being
theoretically better than the other. Rather, we should state cautiously that
for a particular dataset, one type of knots is more effective than the other
when $k$ is restricted in size. For the two datasets exemplified in this
section, we will see that for small to moderate $k$, the BMC longitudinal data
favor non-uniform knots and general P-spline, while the fossil shell data
favor uniform knots and standard P-spline. Nevertheless, as $k$ becomes
adequately big, exact choice of knot placement strategy or P-spline estimator
no longer matters.

Although not relevant to our examples here, we would like to point out that
circumstances where we want to or have to use a small to moderate $k$ include
(1) the iterative algorithm used for model estimation does not converge for
big $k$, (2) computation with a big $k$ takes too much time, etc.

###  6.1 BMC longitudinal data

The Pediatric Bone Mineral Accrual Study (PBMAS) (Bailey, 1997) was launched
in 1991 by the University of Saskatchewan, aiming to investigate the
accumulation of bone mineral content (BMC) in growing children. At first, 109
boys and 113 girls aged 8 $\sim$ 15 from two elementary schools in Saskatoon
(the largest city in Saskatchewan, Canada) were involved for annual scan. 31
new subjects were recruited in 1992 $\sim$ 1993\. Over the years, the study
was interrupted several times, with data collection in 1991 $\sim$ 1997, 2003
$\sim$ 2005, 2007 $\sim$ 2011 and 2016 $\sim$ 2017\. Many subjects dropped out
as they left school, but some managed to return for follow-up scans as young
adults. As of 2017, the dataset has 982 and 1175 measurements from males and
females, respectively. It allows researchers to track the gain of BMC in the
early life and its subsequent loss.

Baxter-Jones et al. (2011) conducted some empirical statistical analysis of
this dataset, finding that the most substantial increase in BMC happens during
6.8 $\sim$ 16.8 years in females and 8.5 $\sim$ 18.5 years in males. After
that, significantly smaller quantity of BMC is gained until the late twenties
or very early thirties. Therefore, optimizing gains during growth allows for a
higher starting point at which BMC loss begins, lowering the risk for
developing osteoporosis.

Elhakeem et al. (2021) studied white ethnic subjects in the dataset
(comprising 112 males and 127 females) using a linear mixed effect model
(Pinheiro and Bates, 2000). For either group, the model is:

$\textrm{bmc}_{j,\textrm{age}}=\sum_{u=1}^{p_{1}}\mathcal{B}_{u}^{\scriptscriptstyle[1]}(\textrm{age})\beta_{u}+\sum_{v=1}^{p_{2}}\mathcal{B}_{v}^{\scriptscriptstyle[2]}(\textrm{age})\alpha_{j,v}+e_{j,\textrm{age}}.$

where $(\mathcal{B}_{u}^{\scriptscriptstyle[1]}(\textrm{age}))_{1}^{p_{1}}$
are cubic B-splines with fixed-effect coefficients $(\beta_{u})_{1}^{p1}$,
$(\mathcal{B}_{v}^{\scriptscriptstyle[2]}\textrm{age})_{1}^{p_{2}}$ are cubic
B-splines with random-effect coefficients
$(\alpha_{j,1},\ldots,\alpha_{j,p_{2}})^{\scriptscriptstyle\textrm{T}}\overset{\scriptscriptstyle\textrm{iid}}{\sim}N(\boldsymbol{0},\boldsymbol{\Sigma_{\alpha}})$
with an unknown covariance $\boldsymbol{\Sigma_{\alpha}}$ and
$e_{j,\textrm{age}}\overset{\scriptscriptstyle\textrm{iid}}{\sim}N(0,\sigma_{e}^{2})$
is a measurement error. To clarify, the random effects are uncorrelated
between subjects (i.e., $\textrm{cov}(\alpha_{i,s},\alpha_{j,t})=0$, $\forall
i\neq j$ and $\forall s,t=1,2,\ldots,p_{2}$) but correlated between components
(i.e., $\textrm{cov}(\alpha_{j,s},\alpha_{j,t})=\sigma_{s,t}$, $\forall s\neq
t$ and $\forall j$). They also differ in marginal variance (i.e.,
$\textrm{var}(\alpha_{j,s})=\sigma_{s}^{2}$, $\forall j$). The first summation
in the model is the population BMC trajectory, and the second is a random
deviation of subject $j$’s BMC trajectory from the former. The two sets of
B-splines may be constructed on different knots. Often, there is
$p_{2}<p_{1}$, so that a random trajectory is expected to be simpler than the
population trajectory.

This model has several limitations. Firstly, smoothness of those trajectories
is not regulated. This is most true for the population trajectory which is
estimated by regression B-splines. As is discussed in section 3.1, such
estimator is sensitive to exact knot location when $p_{1}$ is small and would
overfit data when $p_{1}$ is big. The random trajectory may also end up
wiggly. Although the covariance matrix $\boldsymbol{\Sigma_{\alpha}}$ for the
random-effect B-spline coefficients is kind of a penalty, it has no assumed
structure, so it will not act as a proper wiggliness penalty.

Secondly, estimation of random effects is problematic when using standard R
packages for linear mixed effect models, namely nlme (Pinheiro et al., 2021)
and lme4 (Bates et al., 2015, 2021). Either the estimation algorithm fails to
converge, or the estimated covariance matrix
$\boldsymbol{\hat{\Sigma}_{\alpha}}$ is singular. This is not resolved even if
$p_{2}$ is as small as 5 or 6, and $p_{2}=4$ is the only feasible choice. Note
that in this case a cubic spline has no interior knots so that it is actually
a cubic polynomial. As a result, complexity of the random trajectory allowed
by this model is very restricted in practice.

One way to tackle the above computational difficulty is to assume independent
random effects instead, i.e., $\boldsymbol{\Sigma_{\alpha}}$ is diagonal. Then
it is possible to attempt $p_{2}=10$. But this does not make sense because
$(\mathcal{B}_{v}^{\scriptscriptstyle[2]})_{1}^{p_{2}}$ are not orthogonal to
each other. A more reasonable workaround is to first orthogonalize
$(\mathcal{B}_{v}^{\scriptscriptstyle[2]})_{1}^{p_{2}}$ by QR decomposition,
then assume independent random-effect coefficients for the new basis vectors.
Nevertheless, the concern that trajectories might end up wiggly is still not
addressed.

A better strategy is to consider the following additive mixed model (Zuur,
Saveliev and Ieno, 2014; Pedersen et al., 2019):

$\textrm{bmc}_{j,\textrm{age}}=f(\textrm{age})+f_{j}(\textrm{age})+e_{j,\textrm{age}}.$

Actually, as Figure 2(A) and 2(B) show, BMC has non-constant variance while
log(BMC) does not, thus, a model for log(BMC) is more appropriate (as given in
section 1):

$\log(\textrm{bmc})_{j,\textrm{age}}=f(\textrm{age})+f_{j}(\textrm{age})+e_{j,\textrm{age}},$
(18)

Here, $f$ (a smooth fixed-effect function) is the population log(BMC)
trajectory, $f_{j}$ (smooth random-effect function) is the departure of
subject $j$’s log(BMC) trajectory from $f$ and
$e_{j,\textrm{age}}\overset{\scriptscriptstyle\textrm{iid}}{\sim}N(0,\sigma_{e}^{2})$
is a measurement error. When $f$ and $f_{j}$ are spline estimators regularized
by some wiggliness penalties, model (18) can be estimated in the framework of
generalized additive models (GAMs) (Wood, 2017b) using R package mgcv (Wood,
2021). Advantages of this modeling approach include (1) neither $f$ nor
$f_{j}$ will end up wiggly and (2) estimation of $f$ and $f_{j}$ is
numerically stable.

Figure 15: Estimated BMC trajectories for all subjects, using two P-spline
variants and different $p$ (number of cubic B-splines). Gray vertical dashed
lines mark positions of domain knots ($k=p-4$ interior knots and 2 domain
boundaries). Figure 16: Estimated BMC trajectories with general P-spline
(solid) and standard P-spline (dashed) for subject 13, 42, 68 and 102, when
$p=$ 10 (top row) and 20 (bottom row). Gray dots are BMC measurements. Figure
17: Cross-validation score (solid) and its 2 standard error bounds (dashed)
for general P-spline (black) and standard P-spline (red) at $p=5,6,\ldots,20$.
For BMC dataset, general P-spline with non-uniform knots is more effective
than standard P-spline with uniform knots.

Relevant to the methodology development in this paper, we are to construct $f$
and $f_{j}$ using our proposed general P-spline (section 4.3). For comparison,
we also construct them using Eilers and Marx’s standard P-spline (section
3.4). Let’s still denote the number of B-splines in $f$ and $f_{j}$ by $p_{1}$
and $p_{2}$. We consider the simplest setup where they are represented with
the same set of B-splines so that $p_{1}=p_{2}=p$. For B-spline order, we
choose the usual $d=4$ for cubic B-splines. For the order of difference
penalty, we choose the usual $m=2$ (penalizing the 2nd difference) for $f$,
but $m=1$ for $f_{j}$. This subtlety comes from the fact that fewer subjects
came back for measurements as they grew older. The number of subjects (out of
112 in total) with no observations after the age of 20, 25 and 30 are 36, 50
and 60, respectively. Using $m=2$ for $f_{j}$ includes a random slope in the
null space of $\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(2)}$ and
$\boldsymbol{D}_{p}^{\scriptscriptstyle(2)}$ that would lead to linear
extrapolation for these subjects. But as Figure 2(B) shows, log(BMC)
trajectories are relatively flat after the age of 20, it is thus more
plausible to exclude this random slope. By choosing $m=1$, there will only be
a random intercept in the null space of
$\boldsymbol{\mathcal{D}}_{p}^{\scriptscriptstyle(1)}$ and
$\boldsymbol{D}_{p}^{\scriptscriptstyle(1)}$.

We use automatic strategies (section 3.5) to place knots. We specify the
domain as $[a,b]=[\min(\textrm{age}),\max(\textrm{age})]$. To construct $p$
cubic B-splines we need to place $k=p-4$ interior knots in $(a,b)$ or $k+2$
domain knots on $[a,b]$. Since age data are right-skewed (i.e., more data at
early ages), quantile knots will be very different from uniform knots, so are
non-uniform B-splines from uniform ones. Therefore, a comparison between our
general P-spline and Eilers and Marx’s standard P-spline for varying $p$ (or
equivalently, $k$) is of great interest.

We attempt $p$ = 10 and 20. Estimated BMC trajectories (i.e.,
$\exp(\hat{f}+\hat{f}_{j})$) for all subjects are sketched in Figure 15. When
$p=10$, standard P-spline is more wiggly than general P-spline. It is more
evident in Figure 16, where trajectories for 4 arbitrary subjects are plotted.
Its top row ($p=10$) shows that standard P-spline fit (dashed) exhibits
suspicious wiggliness, whereas general P-spline fit (solid) is satisfactory.
However, Both figures show that when $p=20$, two P-spline variants give
similar fit.

The message is that for BMC dataset, general P-spline on non-uniform knots is
more effective than standard P-spline on uniform knots. Let’s demonstrate this
in a more quantitative way via cross-validation. Note that splitting an
unbalanced longitudinal dataset into training set and testing set needs care,
because if there are no training data for a subject, there is no way to make
prediction for this subject in the testing set. Suppose we want to use 10% of
the data for testing in a random partition,our strategy is to exclude subjects
with very few (like 1 or 2) observations, then sample approximately equal
number of test data from the remaining subjects. This ensures that every
subject is included in a training set. Now for $p=5,6,\ldots,20$, we compute a
cross-validation score for each P-spline variant based on 100 random
partitions. The result is shown in Figure 17. We see that general P-spline is
indeed more effective than standard P-spline for BMC dataset, as its
prediction squared error decreases with $p$ faster. It already achieves a
satisfactory fit at $p=9$, but standard P-spline fit requires $p=17$ to
attains about the same quality. For higher $p$ however, performance of two
P-spline variants is indistinguishable.

###  6.2 Fossil shell data

The fossil dataset (available in R package SemiPar (Wand, 2018)) has raised a
debate on which of uniform and non-uniform knots are superior (Ruppert, Wand
and Carroll, 2009; Wand and Ormerod, 2008; Eilers, Marx and Durbán, 2015).
Here we don’t describe this dataset further and simply treat it as a set of
$(x,y)$ observations. The true data-generating function $g(x)$ is unknown, but
aligned with these discussions, an estimate by R function smooth.spline is
taken as $g(x)$. Note that previous demonstration with non-uniform knots was
via O-spline. Here we replace it with our general P-spline.

Figure 18: General P-spline (black) and standard P-spline (black) v.s.
smooth.spline (red) for fossil data (grey dots). 20 interior knots (grey
vertical dashed lines) are placed.

Eilers, Marx and Durbán (2015) claimed that uniform knots are better than
quantile knots, because as Figure 18 shows, with 20 interior knots, general
P-spline with non-uniform knots fails to capture the local minimum near $x=98$
correctly, while standard P-spline with uniform knots has no trouble with
this. They criticized that placing quantile knots is “the root of all evil”,
as the distribution of knots has a gap near that local minimum, with just a
single knot in approximately $(94,103)$. However, we do not consider it a
sufficient reason to disapprove quantile knots in general. Recall that
quantile knots in our normal mixture example under sampling scheme 2 (Figure
12) are exposed to the same issue, but our general P-spline seems to capture
that local minimum better than standard P-spline. In addition, smooth.spline
actually fits an O-spline rather than an authentic smoothing spline and it
places quantile knots, too! Therefore, if this estimate is deemed truth, it is
self-contradicting to oppose quantile knots.

To complete the picture, there is a need to conduct two experiments with
general P-spline: either (1) insert some knots in the gap or (2) extract the
quantile knots placed by smooth.spline, then estimate general P-spline on the
new knot sequence to see if it improves the fit. Figure 19 (top row) shows the
result of experiment 1. It now turns out that both P-spline variants are
indistinguishable from the smooth.spline fit. Figure 19 (bottom row) shows the
result of experiment 2. It seems that general P-spline is even better than
smoothing spline in approximating data!

Figure 19: Two experiments on fossil dataset. (1) Two knots at 96 and 98
(black squares) are inserted to those 20 interior knots in Figure 18 and this
improves general P-spline fit. Standard P-spline is re-estimated using 22
uniform interior knots. (2) General P-spline is estimated with the same 62
quantile knots used by smooth.spline. Standard P-spline is re-estimated using
62 uniform interior knots.

Once again, we argue that which type of knots appears more advantageous is
conditional $k$ (number of interior knots). To illuminate this, let’s compare
the predictive performance of two P-spline variants at $k=10,11,\ldots,60$ via
cross-validation. In this example, CV score is readily available as $\lambda$
is selected to minimize the GCV score (by contrast, $\lambda$ is chosen to
maximize restricted likelihood in the BMC example). The result is shown in
Figure 20. It is astounding to spot a zigzag in general P-spline’s score. A
clue may be found in Figure 21. For fossil data, it turns out that general
P-spline is very sensitivity to exact knot locations, especially the existence
of a knot at 98. Such sensitivity gets reduced when $k$ increases. After
$k=40$, general P-spline outperforms standard P-spline by a small margin. But
before that, standard P-spline almost always gains the upper hand. Hence, for
small to moderate $k$, standard P-spline is more effective on fossil dataset.

Figure 20: GCV score for general P-spline (black) and standard P-spline (red)
at $k=10,11,\ldots,60$. The score is deterministic hence no error bounds
exist. For fossil dataset, standard P-spline with uniform knots is more
effective than general P-spline with non-uniform knots. Figure 21: A zoomed-in
display on $[92,106]$. General P-spline (black) v.s. smooth.spline (red) for
fossil data. 10, 13, and 16 interior knots are attempted for the former.
Whether a knot exists at 98 matters.

##  7 Conclusions

We have proposed a new penalized regression spline estimator called the
general P-spline, that accommodates B-splines on arbitrary knot sequences,
like non-uniform B-splines on unevenly spaced knots. It is a complement to
Eilers and Marx’s standard P-spline that only makes sense for uniform
B-splines on equidistant knots. At its core, we derived a novel general
difference matrix that accounts for knot spacing, to replace the standard
difference matrix tailored for equidistant knots. Both P-spline variants are
practically useful, because while they show comparable MSE performance in
simulation studies, for a particular real-world dataset one often produces a
more satisfactory fit than the other, especially when the number of knots is
not sufficiently large to mitigate sensitivity in exact knot locations.

Our general P-spline is also closely related to O-spline and smoothing spline.
To our knowledge, we are the first to precisely interpret P-spline’s
difference penalty as the squared $l^{2}$ norm of B-spline coefficients in
$f^{\scriptscriptstyle(m)}(x)$. By contrast, the classic derivative penalty
$\int_{a}^{b}f^{\scriptscriptstyle(m)}(x)^{2}\textrm{d}x$ is the squared
$l^{2}$ norm of $f^{\scriptscriptstyle(m)}(x)$. The former is always a
simplification of the latter, but they have the same null space and equivalent
control over the wiggliness of $f(x)$. We are also the first to link two types
of penalties via the explicit sandwich formula (17). This not only provides a
new method to compute the derivative penalty, but also yields its sparse
“root”, making it straightforward to adapt it to $l^{1}$ regularization which
is worth future research.

Derivation of our general difference matrix has utilized (11a)(11b) for
$f^{\scriptscriptstyle(m)}(x)$ with B-splines on arbitrary knot sequences. We
realized that these results are rarely seen in statistical literature. Even in
methodology development of monotone and shape constrained smoothing (Pya and
Wood, 2015; Wang and Small, 2015), where $f^{\scriptscriptstyle(m)}(x)$ is
subject to constraints so that these results are relevant, it turns out that
equidistant knots were chosen to simplify matters. In fact, standard P-spline
is the foundation of these more advanced smoothing techniques. Clearly, with
the availability of our general P-spline, these smoothing methods can also be
extended for non-uniform B-splines.

[Acknowledgments] Zheyuan Li is grateful to Ahmed Elhakeem (Bristol Medical
School, University of Bristol, UK) and Adam Baxter-Jones (University of
Saskatchewan, Canada), who shared their experience with modeling BMC data in
our previous collaboration. BMC data are not open data. For this paper, a
high-quality synthetic version is used and provided in gps.mgcv. Researchers
who are interested in accessing the authentic data from PBMAS should contact
Adam Baxter-Jones (baxter.jones@usask.ca).

Zheyuan Li is also Assistant Professor, School of Mathematics and Statistics,
Henan University, Kaifeng, Henan, China (zheyuan.li@vip.henu.edu.cn).

Zheyuan Li was supported by National Natural Science Foundation of China under
Young Scientists Fund NSFC-12001166.

Jiguo Cao was supported by Natural Sciences and Engineering Research Council
of Canada (NSERC) under Discovery Grant 2018-06008.

Supplement to “General P-splines for non-uniform B-splines” The supplement
contains R packages gps and gps.mgcv, as well as R code associated with the
paper.

## References

  * Arnold and Tibshirani (2016) [author] Arnold, Taylor B.T. B. and Tibshirani, Ryan J.R. J. (2016). Efficient Implementations of the Generalized Lasso Dual Path Algorithm. Journal of Computational and Graphical Statistics 25 1-27.
  * Bailey (1997) [author] Bailey, Donald AD. A. (1997). The Saskatchewan Pediatric Bone Mineral Accrual Study: Bone Mineral Acquisition During the Growing Years. International Journal of Sports Medicine 18 s191–s194.
  * Bates et al. (2015) [author] Bates, DouglasD., Mächler, MartinM., Bolker, BenB. and Walker, SteveS. (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software 67 1–48.
  * Bates et al. (2021) [author] Bates, DouglasD., Mächler, MartinM., Bolker, BenB. and Walker, SteveS. (2021). lme4: Linear Mixed-Effects Models using ’Eigen’ and S4 R package version 1.1-27.1.
  * Baxter-Jones et al. (2011) [author] Baxter-Jones, Adam DGA. D., Faulkner, Robert AR. A., Forwood, Mark RM. R., Mirwald, Robert LR. L. and Bailey, Donald AD. A. (2011). Bone mineral accrual from 8 to 30 years of age: An estimation of peak bone mass. Journal of Bone and Mineral Research 26 1729–1739.
  * Craven and Wahba (1979) [author] Craven, P.P. and Wahba, G.G. (1979). Smoothing noisy data with spline functions: Estimating the correct degree of smoothing by the method of generalized cross-validation. Numerische Mathematik 31 377–403.
  * de Boor (2001) [author] de Boor, C.C. (2001). A Practical Guide to Splines (Revised Edition). Applied Mathematical Sciences 27\. Springer New York.
  * Eilers and Marx (1996) [author] Eilers, P. H. C.P. H. C. and Marx, B. D.B. D. (1996). Flexible smoothing with b-splines and penalties. Statistical Science 11 89–102.
  * Eilers, Marx and Durbán (2015) [author] Eilers, P. H. C.P. H. C., Marx, B. D.B. D. and Durbán, M.M. (2015). Twenty years of p-splines. Statistics and Operations Research Transactions 39 149–186.
  * Eilers and Marx (2021) [author] Eilers, P. H.P. H. and Marx, B. D.B. D. (2021). Practical smoothing: the joys of P-splines. Cambridge University Press.
  * Elhakeem et al. (2021) [author] Elhakeem, AhmedA., Hughes, RachaelR., Tilling, KateK., Cousminer, DianaD., Jackowski, StefanS., Cole, TimT., Kwong, AlexA., Li, ZheyuanZ., Grant, StruanS., Baxter-Jones, AdamA., Zemel, BabetteB. and Lawlor, DeborahD. (2021). Using linear and natural cubic splines, SITAR, and latent trajectory models to characterise nonlinear longitudinal growth trajectories in cohort studies. medRxiv. 10.1101/2021.05.26.21257519
  * Jupp (1978) [author] Jupp, David L. BD. L. B. (1978). Approximation to Data by Splines with Free Knots. SIAM journal on numerical analysis 15 328–343.
  * Li and Cao (2022) [author] Li, ZheyuanZ. and Cao, JiguoJ. (2022). Supplement to ”General P-splines for non-uniform B-splines”. 10.1214/
  * Luo and Wahba (1997) [author] Luo, ZhenZ. and Wahba, GraceG. (1997). Hybrid Adaptive Splines. Journal of the American Statistical Association 92 107–116.
  * Mammen and van de Geer (1997) [author] Mammen, EnnoE. and van de Geer, SaraS. (1997). Locally Adaptive Regression Splines. The Annals of statistics 25 387–413.
  * O’Sullivan (1986) [author] O’Sullivan, F.F. (1986). A statistical perspective on ill-posed inverse problems. Statistical Science 1 502–518.
  * Pedersen et al. (2019) [author] Pedersen, Eric JE. J., Miller, David LD. L., Simpson, Gavin LG. L. and Ross, NoamN. (2019). Hierarchical generalized additive models in ecology: an introduction with mgcv. PeerJ 7 e6876–e6876.
  * Pinheiro and Bates (2000) [author] Pinheiro, José CJ. C. and Bates, Douglas MD. M. (2000). Mixed-Effects Models in S and S-PLUS. Springer.
  * Pinheiro et al. (2021) [author] Pinheiro, JoseJ., Bates, DouglasD., DebRoy, SaikatS., Sarkar, DeepayanD. and R Core Team (2021). nlme: Linear and Nonlinear Mixed Effects Models R package version 3.1-153.
  * Pya and Wood (2015) [author] Pya, NatalyaN. and Wood, Simon NS. N. (2015). Shape constrained additive models. Statistics and Computing 25 543–559.
  * Ramdas and Tibshirani (2016) [author] Ramdas, AadityaA. and Tibshirani, Ryan JR. J. (2016). Fast and Flexible ADMM Algorithms for Trend Filtering. Journal of computational and graphical statistics 25 839–858.
  * Reinsch (1967) [author] Reinsch, C. H.C. H. (1967). Smoothing by spline functions. Numerische Mathematik 10 177–183.
  * Reinsch (1971) [author] Reinsch, C. H.C. H. (1971). Smoothing by spline functions. II. Numerische Mathematik 16 451–454.
  * Ruppert, Wand and Carroll (2003) [author] Ruppert, D.D., Wand, M. P.M. P. and Carroll, R. J.R. J. (2003). Semiparametric Regression. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press.
  * Ruppert, Wand and Carroll (2009) [author] Ruppert, D.D., Wand, M. P.M. P. and Carroll, R. J.R. J. (2009). Semiparametric regression during 2003-2007. Electronic journal of statistics 3 1193–1256.
  * Silverman (1985) [author] Silverman, B. WB. W. (1985). Some Aspects of the Spline Smoothing Approach to Non-Parametric Regression Curve Fitting. Journal of the Royal Statistical Society: Series B (Statistical Methodological) 47 1–52.
  * Spiriti et al. (2013) [author] Spiriti, S.S., Eubank, R.R., Smith, P. W.P. W. and Young, D.D. (2013). Knot selection for least-squares and penalized splines. Journal of statistical computation and simulation 83 1020–1036.
  * Tibshirani (1996) [author] Tibshirani, RobertR. (1996). Regression Shrinkage and Selection Via the Lasso. Journal of the Royal Statistical Society: Series B (Methodological) 58 267-288.
  * Tibshirani (2014) [author] Tibshirani, Ryan J.R. J. (2014). Adaptive Piecewise Polynomial Estimation via Trend Filtering. The Annals of statistics 42 285–323.
  * Tibshirani and Taylor (2011) [author] Tibshirani, Ryan J.R. J. and Taylor, JonathanJ. (2011). The Solution Path of the Generalized LASSO. The Annals of statistics 39 1335–1371.
  * Wahba (1990) [author] Wahba, G.G. (1990). Spline models for observational data. Society for Industrial and Applied Mathematics.
  * Wand (2018) [author] Wand, MattM. (2018). SemiPar: Semiparametic Regression R package version 1.0-4.2.
  * Wand and Ormerod (2008) [author] Wand, M. P.M. P. and Ormerod, J. T.J. T. (2008). On semiparametric regression with O’Sullivan penalized splines. Australian & New Zealand journal of statistics 50 179–198.
  * Wang (2011) [author] Wang, Y.Y. (2011). Smoothing Splines: Methods and Applications. Monographs on statistics and applied probability 121\. Chapman and Hall.
  * Wang and Small (2015) [author] Wang, WeiW. and Small, Dylan SD. S. (2015). Monotone B-Spline Smoothing for a Generalized Linear Model Response. The American Statistician 69 28–33.
  * Wood (2017a) [author] Wood, S. N.S. N. (2017a). P-splines with derivative based penalties and tensor product smoothing of unevenly distributed data. Statistics and Computing 27 985–989.
  * Wood (2017b) [author] Wood, S. N.S. N. (2017b). Generalized Additive Models: An Introduction with R, 2nd ed. Chapman and Hall.
  * Wood (2021) [author] Wood, SimonS. (2021). mgcv: Mixed GAM Computation Vehicle with Automatic Smoothness Estimation R package version 1.8-38.
  * Zuur, Saveliev and Ieno (2014) [author] Zuur, Alain F.A. F., Saveliev, Anatoly A.A. A. and Ieno, Elena N.E. N. (2014). A Beginner’s Guide to Generalized Additive Mixed Models with R. Highland Statistics.