Setup
Simulation
We simulate a linear model with \(p = 95\) predictors and \(n = 50\) observations. The true coefficient vector \(\beta\) is piecewise constant: it is zero over three large segments (sizes 25, 25, 25) and takes values \(+1\) and \(-1\) over two smaller active blocks (sizes 10 each).
set.seed(42)
p <- 95
n <- 50
beta <- rep(c(0, 1, 0, -1, 0), c(25, 10, 25, 10, 25))
## Block-correlated predictors: high within-segment, zero across segments
rho <- 0.75
Soo <- toeplitz(rho^(0:24)) # Toeplitz correlation within zero segments
Sww <- matrix(rho, 10, 10) # constant correlation within active segments
Sigma <- bdiag(Soo, Sww, Soo, Sww, Soo)
diag(Sigma) <- 1
X <- as.matrix(matrix(rnorm(p * n), n, p) %*% chol(Sigma))
y <- X %*% beta + rnorm(n, 0, 10)
## Segment labels (for plot legends)
segments <- rep(1:5, c(25, 10, 25, 10, 25))
segment_names <- paste0("seg", 1:5)
seg_labels <- segment_names[segments]The five segments are clearly visible in the true signal:
data.frame(index = seq_along(beta), beta = beta, segment = seg_labels) |>
ggplot(aes(x = index, y = beta, colour = segment)) +
geom_step() + geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.4) +
labs(x = "Predictor index", y = expression(beta), title = "True coefficient vector") +
theme_bw()
Community graph and its Laplacian
The structural prior we wish to encode is: predictors within the same segment should have similar coefficients. This is captured by a community graph whose edges connect all pairs of predictors belonging to the same segment (a union of disjoint cliques).
## Adjacency matrix: clique within each segment, no edges across
Ioo <- matrix(1, 25, 25)
Iww <- matrix(1, 10, 10)
A <- as.matrix(bdiag(Ioo, Iww, Ioo, Iww, Ioo))
diag(A) <- 0
## Graph Laplacian: L = D - A (regularized for positive definiteness)
L2 <- -A
diag(L2) <- colSums(A) + 1e-2The Laplacian \(L_2\) defines the quadratic form
\[\beta^\top L_2\, \beta = \sum_{(i,j)\in E} (\beta_i - \beta_j)^2,\]
which penalizes differences between connected predictors.
Adding \(\lambda_2 \beta^\top L_2 \beta /
2\) to the loss therefore encourages within-segment coefficients
to be equal — a perfect match for the piecewise-constant signal. Note
that \(L_2\) is passed as the
struct argument to ridge() and
elastic_net().
For the Fused Lasso, the graph structure is encoded differently:
struct takes the adjacency matrix \(A\) directly, since the fusion penalty
is
\[\lambda_2 \sum_{(i,j)\in E} |\beta_i - \beta_j|,\]
an \(\ell_1\) (not \(\ell_2\)) penalty on pairwise differences. The graph topology is the same, but the penalty form differs.
if (requireNamespace("igraph", quietly = TRUE)) {
g <- igraph::graph_from_adjacency_matrix(A, mode = "undirected")
cols <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#CC79A7")
vertex_colors <- cols[segments]
igraph::plot.igraph(
g, vertex.color = vertex_colors, vertex.size = 3,
vertex.label = NA, edge.width = 0.1,
main = "Community graph (colour = segment)"
)
legend("topleft", legend = segment_names, fill = cols, bty = "n", cex = 0.8)
}
