get_2ndDerivPenalty: Compute Integrated Squared Second Derivative Penalty Matrix
In lgspline: Lagrangian Multiplier Smoothing Splines for Smooth Function Estimation
View source: R/HelperFunctions.R
get_2ndDerivPenalty R Documentation
Compute Integrated Squared Second Derivative Penalty Matrix
Description
Computes the p \times p integrated squared second-derivative penalty matrix
\boldsymbol{\Lambda}_s for one partition of a monomial spline,
such that
\boldsymbol{\beta}_k^\top \boldsymbol{\Lambda}_s
\boldsymbol{\beta}_k
= \int_{\mathbf{a}}^{\mathbf{b}}
\|\tilde{f}_k''(\mathbf{t})\|^2 \, d\mathbf{t},
where \mathbf{a} and \mathbf{b} are the observed
predictor minimums and maximums, and
\tilde{f}_k(\mathbf{t}) = \mathbf{x}_k^\top
\boldsymbol{\beta}_k is the fitted function for partition
\mathcal{P}_k.
The implementation uses a general monomial derivative rule that
handles all term types (marginal powers, two-way interactions,
quadratic interactions, three-way interactions) in a single unified
loop.
Usage
get_2ndDerivPenalty(
colnm_expansions,
C,
power1_cols,
power2_cols,
power3_cols,
power4_cols,
interaction_single_cols,
interaction_quad_cols,
triplet_cols,
p_expansions,
select_cols = NULL
)
Arguments
colnm_expansions
Character vector of length p giving
column names for the basis expansion matrix \mathbf{X}_k.
Each name encodes the monomial structure: predictor names appear
literally for linear terms, with ^d suffixes for degree
d, and an x separator between factors in interaction
terms (e.g.\ "_1_alphax_2_beta^2" for
t_1 t_2^2). Predictor names must not be substrings of one
another.
C
Numeric N \times p matrix of basis expansions
\mathbf{X}_k. Used only to read the observed range
[a_j, b_j] of each predictor via min/max of
the linear columns.
power1_cols
Integer vector. Column indices of linear terms
t_j.
power2_cols
Integer vector. Column indices of quadratic terms
t_j^2.
power3_cols
Integer vector. Column indices of cubic terms
t_j^3.
power4_cols
Integer vector. Column indices of quartic terms
t_j^4.
interaction_single_cols
Integer vector. Column indices of
linear-by-linear interaction terms t_j t_l.
interaction_quad_cols
Integer vector. Column indices of
linear-by-quadratic interaction terms t_j t_l^2 and
t_j^2 t_l.
triplet_cols
Integer vector. Column indices of three-way
interaction terms t_j t_l t_m.
select_cols
Optional integer vector of predictor indices
(positions within power1_cols) whose curvature operators
D_v are summed. Defaults to all q predictors.
Details
Mathematical framework
Let the basis expansion for partition \mathcal{P}_k be
\mathbf{x}_k = (\phi_1(\mathbf{t}), \ldots,
\phi_p(\mathbf{t}))^\top where each \phi_i is a
multivariate monomial
\phi_i(\mathbf{t})
= \prod_{j=1}^{q} t_j^{\alpha_{ij}}.
The second derivative of \tilde{f}_k decomposes into
q total curvature operators, one per predictor. For
predictor v:
D_v = \frac{\partial^2}{\partial t_v^2}
+ \sum_{s \neq v} \frac{\partial^2}{\partial t_v \, \partial t_s}.
The monomial derivative rule gives each second partial derivative
in closed form. For the pure second derivative (r = s = v):
\frac{\partial^2}{\partial t_v^2}
\prod_{j} t_j^{\alpha_j}
= \alpha_v(\alpha_v - 1) \;
t_v^{\alpha_v - 2} \prod_{j \neq v} t_j^{\alpha_j}.
For a mixed second derivative (s \neq v):
\frac{\partial^2}{\partial t_v \, \partial t_s}
\prod_{j} t_j^{\alpha_j}
= \alpha_v \alpha_s \;
t_v^{\alpha_v - 1} t_s^{\alpha_s - 1}
\prod_{j \neq v,s} t_j^{\alpha_j}.
In both cases a term is zero when the required exponent would be
negative (e.g.\ \alpha_v < 2 for the pure case). Applying
D_v to \phi_i produces a sum of monomials with known
coefficients and exponent vectors.
Integration
Because every D_v(\phi_i) is polynomial, the product
D_v(\phi_i) \, D_v(\phi_j) is also polynomial and the
multivariate integral factorises over predictors:
\int_{\mathbf{a}}^{\mathbf{b}}
\prod_{j=1}^{q} t_j^{e_j} \, d\mathbf{t}
= \prod_{j=1}^{q}
\frac{b_j^{e_j+1} - a_j^{e_j+1}}{e_j + 1}.
Notably, this integral runs over all q predictor
ranges, including predictors that do not appear in the integrand
(for which e_j = 0 and the factor reduces to
b_j - a_j).
Single-predictor verification
For q = 1 with expansion
\mathbf{x} = (1, t, t^2, t^3)^\top on [a, b], the
penalty matrix reduces to
\boldsymbol{\Lambda}_s
= \int_a^b \mathbf{x}'' \mathbf{x}''^\top \, dt
= \begin{pmatrix}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 4(b - a) & 6(b^2 - a^2) \\
0 & 0 & 6(b^2 - a^2) & 12(b^3 - a^3)
\end{pmatrix},
matching the formula in Section 2.3.
Value
A symmetric positive semi-definite p \times p matrix
\boldsymbol{\Lambda}_s with
[\boldsymbol{\Lambda}_s]_{ij}
= \sum_{v=1}^{q} \int_{\mathbf{a}}^{\mathbf{b}}
D_v(\phi_i) \, D_v(\phi_j) \, d\mathbf{t}.
Naming convention
Column names in colnm_expansions must encode each monomial's
predictor content. The parser checks each predictor name (taken
from the linear columns) against the column name:
-
predname^d (literal caret) signals exponent d.
-
predname without ^ signals exponent 1.
Absence of predname signals exponent 0.
Multiple factors in an interaction are separated by x.
Predictor names must be unique and must not be substrings of one
another.
References
Reinsch, C. H. (1967). Smoothing by spline functions.
Numerische Mathematik, 10(3), 177–183.
See Also
lgspline for the full model fitting
interface.
Examples
## Verification example: 3-predictor model with all term types:
# This example constructs the penalty matrix analytically and then
# verifies selected entries against closed-form hand calculations.
# Users can extend or adapt this example to audit new basis
# expansions.
set.seed(1234)
n <- 2000
t1 <- runif(n, 1, 4) # predictor 1, support [1, 4]
t2 <- runif(n, -2, 3) # predictor 2, support [-2, 3]
t3 <- runif(n, 0.5, 2) # predictor 3, support [0.5, 2]
## Predictor names (must not be substrings of one another)
pn <- c("_1_aa", "_2_bb", "_3_cc")
## Build column names encoding the monomial structure
col_names <- c(
pn, # linear
paste0(pn, "^2"), # quadratic
paste0(pn, "^3"), # cubic
paste0(pn[1], "x", pn[2]), # t1*t2
paste0(pn[1], "x", pn[3]), # t1*t3
paste0(pn[2], "x", pn[3]), # t2*t3
paste0(pn[2], "x", pn[1], "^2"), # t1^2*t2
paste0(pn[1], "x", pn[2], "^2"), # t1*t2^2
paste0(pn[1], "x", pn[2], "x", pn[3]) # t1*t2*t3
)
p_expansions <- length(col_names) # 14 basis functions
## Build the expansion matrix C
C <- cbind(
t1, t2, t3,
t1^2, t2^2, t3^2,
t1^3, t2^3, t3^3,
t1*t2, t1*t3, t2*t3,
t2*t1^2, t1*t2^2,
t1*t2*t3
)
colnames(C) <- col_names
## Compute the penalty matrix
Ls <- get_2ndDerivPenalty(
colnm_expansions = col_names,
C = C,
power1_cols = 1:3,
power2_cols = 4:6,
power3_cols = 7:9,
power4_cols = integer(0),
interaction_single_cols = 10:12,
interaction_quad_cols = 13:14,
triplet_cols = 15,
p_expansions = p_expansions,
select_cols = 1:3
)
## Hand-computed reference values (exact, using true bounds) ---
# Notation: dt1 = 4-1 = 3, dt2 = 3-(-2) = 5, dt3 = 2-0.5 = 1.5
#
# Entry [4,4]: t1^2 diagonal.
# D_1(t1^2) = 2. No other D_v contributes.
# integral (2)^2 dt1 dt2 dt3 = 4 * 3 * 5 * 1.5 = 90
#
# Entry [10,10]: t1*t2 diagonal.
# D_1(t1*t2) = 1 (mixed d^2/dt1 dt2).
# D_2(t1*t2) = 1 (mixed d^2/dt2 dt1).
# integral 1 dt1 dt2 dt3 + integral 1 dt1 dt2 dt3 = 2*3*5*1.5 = 45
#
# Entry [15,15]: t1*t2*t3 diagonal.
# D_1(t1*t2*t3) = t3 + t2 (mixed partials d^2/dt1 dt2 and d^2/dt1 dt3)
# D_2(t1*t2*t3) = t1 + t3 (similarly)
# D_3(t1*t2*t3) = t1 + t2 (similarly)
# Full integral = sum of 3 terms = 120 + 337.5 + 266.25 = 723.75
## Compare (allowing ~0.5% tolerance for data-derived bounds)
cat("Entry [4,4]: analytical =", round(Ls[4,4], 2),
" exact = 90.00\n")
cat("Entry [10,10]: analytical =", round(Ls[10,10], 2),
" exact = 45.00\n")
cat("Entry [15,15]: analytical =", round(Ls[15,15], 2),
" exact = 723.75\n")
cat("Symmetric:", isSymmetric(Ls), "\n")
lgspline documentation built on Aug. 5, 2026, 1:10 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
View source: R/HelperFunctions.R
| get_2ndDerivPenalty | R Documentation |
Compute Integrated Squared Second Derivative Penalty Matrix
Description
Computes the p \times p integrated squared second-derivative penalty matrix
\boldsymbol{\Lambda}_s for one partition of a monomial spline,
such that
\boldsymbol{\beta}_k^\top \boldsymbol{\Lambda}_s
\boldsymbol{\beta}_k
= \int_{\mathbf{a}}^{\mathbf{b}}
\|\tilde{f}_k''(\mathbf{t})\|^2 \, d\mathbf{t},
where \mathbf{a} and \mathbf{b} are the observed
predictor minimums and maximums, and
\tilde{f}_k(\mathbf{t}) = \mathbf{x}_k^\top
\boldsymbol{\beta}_k is the fitted function for partition
\mathcal{P}_k.
The implementation uses a general monomial derivative rule that handles all term types (marginal powers, two-way interactions, quadratic interactions, three-way interactions) in a single unified loop.
Usage
get_2ndDerivPenalty(
colnm_expansions,
C,
power1_cols,
power2_cols,
power3_cols,
power4_cols,
interaction_single_cols,
interaction_quad_cols,
triplet_cols,
p_expansions,
select_cols = NULL
)
Arguments
colnm_expansions |
Character vector of length |
C |
Numeric |
power1_cols |
Integer vector. Column indices of linear terms
|
power2_cols |
Integer vector. Column indices of quadratic terms
|
power3_cols |
Integer vector. Column indices of cubic terms
|
power4_cols |
Integer vector. Column indices of quartic terms
|
interaction_single_cols |
Integer vector. Column indices of
linear-by-linear interaction terms |
interaction_quad_cols |
Integer vector. Column indices of
linear-by-quadratic interaction terms |
triplet_cols |
Integer vector. Column indices of three-way
interaction terms |
select_cols |
Optional integer vector of predictor indices
(positions within |
Details
Mathematical framework
Let the basis expansion for partition \mathcal{P}_k be
\mathbf{x}_k = (\phi_1(\mathbf{t}), \ldots,
\phi_p(\mathbf{t}))^\top where each \phi_i is a
multivariate monomial
\phi_i(\mathbf{t})
= \prod_{j=1}^{q} t_j^{\alpha_{ij}}.
The second derivative of \tilde{f}_k decomposes into
q total curvature operators, one per predictor. For
predictor v:
D_v = \frac{\partial^2}{\partial t_v^2}
+ \sum_{s \neq v} \frac{\partial^2}{\partial t_v \, \partial t_s}.
The monomial derivative rule gives each second partial derivative
in closed form. For the pure second derivative (r = s = v):
\frac{\partial^2}{\partial t_v^2}
\prod_{j} t_j^{\alpha_j}
= \alpha_v(\alpha_v - 1) \;
t_v^{\alpha_v - 2} \prod_{j \neq v} t_j^{\alpha_j}.
For a mixed second derivative (s \neq v):
\frac{\partial^2}{\partial t_v \, \partial t_s}
\prod_{j} t_j^{\alpha_j}
= \alpha_v \alpha_s \;
t_v^{\alpha_v - 1} t_s^{\alpha_s - 1}
\prod_{j \neq v,s} t_j^{\alpha_j}.
In both cases a term is zero when the required exponent would be
negative (e.g.\ \alpha_v < 2 for the pure case). Applying
D_v to \phi_i produces a sum of monomials with known
coefficients and exponent vectors.
Integration
Because every D_v(\phi_i) is polynomial, the product
D_v(\phi_i) \, D_v(\phi_j) is also polynomial and the
multivariate integral factorises over predictors:
\int_{\mathbf{a}}^{\mathbf{b}}
\prod_{j=1}^{q} t_j^{e_j} \, d\mathbf{t}
= \prod_{j=1}^{q}
\frac{b_j^{e_j+1} - a_j^{e_j+1}}{e_j + 1}.
Notably, this integral runs over all q predictor
ranges, including predictors that do not appear in the integrand
(for which e_j = 0 and the factor reduces to
b_j - a_j).
Single-predictor verification
For q = 1 with expansion
\mathbf{x} = (1, t, t^2, t^3)^\top on [a, b], the
penalty matrix reduces to
\boldsymbol{\Lambda}_s
= \int_a^b \mathbf{x}'' \mathbf{x}''^\top \, dt
= \begin{pmatrix}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 4(b - a) & 6(b^2 - a^2) \\
0 & 0 & 6(b^2 - a^2) & 12(b^3 - a^3)
\end{pmatrix},
matching the formula in Section 2.3.
Value
A symmetric positive semi-definite p \times p matrix
\boldsymbol{\Lambda}_s with
[\boldsymbol{\Lambda}_s]_{ij}
= \sum_{v=1}^{q} \int_{\mathbf{a}}^{\mathbf{b}}
D_v(\phi_i) \, D_v(\phi_j) \, d\mathbf{t}.
Naming convention
Column names in colnm_expansions must encode each monomial's
predictor content. The parser checks each predictor name (taken
from the linear columns) against the column name:
-
predname^d(literal caret) signals exponentd. -
prednamewithout^signals exponent 1. Absence of
prednamesignals exponent 0.
Multiple factors in an interaction are separated by x.
Predictor names must be unique and must not be substrings of one
another.
References
Reinsch, C. H. (1967). Smoothing by spline functions. Numerische Mathematik, 10(3), 177–183.
See Also
lgspline for the full model fitting
interface.
Examples
## Verification example: 3-predictor model with all term types:
# This example constructs the penalty matrix analytically and then
# verifies selected entries against closed-form hand calculations.
# Users can extend or adapt this example to audit new basis
# expansions.
set.seed(1234)
n <- 2000
t1 <- runif(n, 1, 4) # predictor 1, support [1, 4]
t2 <- runif(n, -2, 3) # predictor 2, support [-2, 3]
t3 <- runif(n, 0.5, 2) # predictor 3, support [0.5, 2]
## Predictor names (must not be substrings of one another)
pn <- c("_1_aa", "_2_bb", "_3_cc")
## Build column names encoding the monomial structure
col_names <- c(
pn, # linear
paste0(pn, "^2"), # quadratic
paste0(pn, "^3"), # cubic
paste0(pn[1], "x", pn[2]), # t1*t2
paste0(pn[1], "x", pn[3]), # t1*t3
paste0(pn[2], "x", pn[3]), # t2*t3
paste0(pn[2], "x", pn[1], "^2"), # t1^2*t2
paste0(pn[1], "x", pn[2], "^2"), # t1*t2^2
paste0(pn[1], "x", pn[2], "x", pn[3]) # t1*t2*t3
)
p_expansions <- length(col_names) # 14 basis functions
## Build the expansion matrix C
C <- cbind(
t1, t2, t3,
t1^2, t2^2, t3^2,
t1^3, t2^3, t3^3,
t1*t2, t1*t3, t2*t3,
t2*t1^2, t1*t2^2,
t1*t2*t3
)
colnames(C) <- col_names
## Compute the penalty matrix
Ls <- get_2ndDerivPenalty(
colnm_expansions = col_names,
C = C,
power1_cols = 1:3,
power2_cols = 4:6,
power3_cols = 7:9,
power4_cols = integer(0),
interaction_single_cols = 10:12,
interaction_quad_cols = 13:14,
triplet_cols = 15,
p_expansions = p_expansions,
select_cols = 1:3
)
## Hand-computed reference values (exact, using true bounds) ---
# Notation: dt1 = 4-1 = 3, dt2 = 3-(-2) = 5, dt3 = 2-0.5 = 1.5
#
# Entry [4,4]: t1^2 diagonal.
# D_1(t1^2) = 2. No other D_v contributes.
# integral (2)^2 dt1 dt2 dt3 = 4 * 3 * 5 * 1.5 = 90
#
# Entry [10,10]: t1*t2 diagonal.
# D_1(t1*t2) = 1 (mixed d^2/dt1 dt2).
# D_2(t1*t2) = 1 (mixed d^2/dt2 dt1).
# integral 1 dt1 dt2 dt3 + integral 1 dt1 dt2 dt3 = 2*3*5*1.5 = 45
#
# Entry [15,15]: t1*t2*t3 diagonal.
# D_1(t1*t2*t3) = t3 + t2 (mixed partials d^2/dt1 dt2 and d^2/dt1 dt3)
# D_2(t1*t2*t3) = t1 + t3 (similarly)
# D_3(t1*t2*t3) = t1 + t2 (similarly)
# Full integral = sum of 3 terms = 120 + 337.5 + 266.25 = 723.75
## Compare (allowing ~0.5% tolerance for data-derived bounds)
cat("Entry [4,4]: analytical =", round(Ls[4,4], 2),
" exact = 90.00\n")
cat("Entry [10,10]: analytical =", round(Ls[10,10], 2),
" exact = 45.00\n")
cat("Entry [15,15]: analytical =", round(Ls[15,15], 2),
" exact = 723.75\n")
cat("Symmetric:", isSymmetric(Ls), "\n")
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.