L1 Trend Filtering
Introduction
Kim et al. (2009) propose the \(l_1\) trend filtering method for trend estimation. The method solves an optimization problem of the form
\[ \begin{array}{ll} \underset{\beta}{\mbox{minimize}} & \frac{1}{2}\sum_{i=1}^m (y_i - \beta_i)^2 + \lambda ||D\beta||_1 \end{array} \] where the variable to be estimated is \(\beta\) and we are given the problem data \(y\) and \(\lambda\). The matrix \(D\) is the second-order difference matrix,
\[ D = \left[ \begin{matrix} 1 & -2 & 1 & & & & \\ & 1 & -2 & 1 & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & 1 & -2 & 1 & \\ & & & & 1 & -2 & 1\\ \end{matrix} \right]. \]
The implementation is in both C and Matlab. Hadley Wickham provides an R interface to the C code. This is on GitHub and can be installed via:
library(devtools)
install_github("hadley/l1tf")
Example
We will use the example in l1tf
to illustrate. The package provides
the function l1tf
which computes the trend estimate for a specified
\(\lambda\).
sp_data <- data.frame(x = sp500$date,
y = sp500$log,
l1_50 = l1tf(sp500$log, lambda = 50),
l1_100 = l1tf(sp500$log, lambda = 100))
The CVXR
version
CVXR
provides all the atoms and functions necessary to formulat the
problem in a few lines. For example, the \(D\) matrix above is provided
by the function diff(..., differences = 2)
. Notice how the
formulation tracks the mathematical construct above.
## lambda = 50
y <- sp500$log
lambda_1 <- 50
beta <- Variable(length(y))
objective_1 <- Minimize(0.5 * p_norm(y - beta) +
lambda_1 * p_norm(diff(x = beta, differences = 2), 1))
p1 <- Problem(objective_1)
betaHat_50 <- solve(p1)$getValue(beta)
## lambda = 100
lambda_2 <- 100
objective_2 <- Minimize(0.5 * p_norm(y - beta) +
lambda_2 * p_norm(diff(x = beta, differences = 2), 1))
p2 <- Problem(objective_2)
betaHat_100 <- solve(p2)$getValue(beta)
NOTE Of course, CVXR
is much slower since it is not optimized just
for one problem.
Comparison Plots
A plot of the estimates for two values of \(\lambda\) is shown below
using both approaches. First the l1tf
plot.
ggplot(data = sp_data) +
geom_line(mapping = aes(x = x, y = y), color = 'grey50') +
labs(x = "Date", y = "SP500 log-price") +
geom_line(mapping = aes(x = x, y = l1_50), color = 'red', size = 1) +
geom_line(mapping = aes(x = x, y = l1_100), color = 'blue', size = 1)

Figure 1: \(L_1\) trends for \(\lambda = 50\) (red) and \(\lambda = 100\) (blue).
Next the corresponding CVXR
plots.
cvxr_data <- data.frame(x = sp500$date,
y = sp500$log,
l1_50 = betaHat_50,
l1_100 = betaHat_100)
ggplot(data = cvxr_data) +
geom_line(mapping = aes(x = x, y = y), color = 'grey50') +
labs(x = "Date", y = "SP500 log-price") +
geom_line(mapping = aes(x = x, y = l1_50), color = 'red', size = 1) +
geom_line(mapping = aes(x = x, y = l1_100), color = 'blue', size = 1)

Figure 2: CVXR
estimated \(L_1\) trends for \(\lambda = 50\) (red) and \(\lambda = 100\) (blue).
Notes
The CVXR
solution is not quite exactly that of l1tf
: on the left it shows a larger difference for the two
\(\lambda\) values; in the middle, it is less flatter than l1tf
; and
on the right, it does not have as many knots as l1tf
.
Session Info
sessionInfo()
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin19.5.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS/LAPACK: /usr/local/Cellar/openblas/0.3.10_1/lib/libopenblasp-r0.3.10.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices datasets utils methods base
##
## other attached packages:
## [1] l1tf_0.0.0.9000 ggplot2_3.3.2 CVXR_1.0-9
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.1.0 xfun_0.15 slam_0.1-47 purrr_0.3.4
## [5] lattice_0.20-41 Rmosek_9.2.3 colorspace_1.4-1 vctrs_0.3.2
## [9] generics_0.0.2 htmltools_0.5.0 yaml_2.2.1 gmp_0.6-0
## [13] rlang_0.4.7 pillar_1.4.6 glue_1.4.1 Rmpfr_0.8-1
## [17] withr_2.2.0 Rcplex_0.3-3 bit64_0.9-7 lifecycle_0.2.0
## [21] stringr_1.4.0 munsell_0.5.0 blogdown_0.19 gtable_0.3.0
## [25] gurobi_9.0.3.1 codetools_0.2-16 evaluate_0.14 labeling_0.3
## [29] knitr_1.28 cccp_0.2-4 highr_0.8 Rcpp_1.0.5
## [33] scales_1.1.1 farver_2.0.3 bit_1.1-15.2 digest_0.6.25
## [37] stringi_1.4.6 bookdown_0.19 dplyr_1.0.0 grid_4.0.2
## [41] Rglpk_0.6-4 ECOSolveR_0.5.3 tools_4.0.2 magrittr_1.5
## [45] tibble_3.0.3 crayon_1.3.4 pkgconfig_2.0.3 ellipsis_0.3.1
## [49] rcbc_0.1.0.9001 Matrix_1.2-18 assertthat_0.2.1 rmarkdown_2.3
## [53] R6_2.4.1 compiler_4.0.2
Source
References
Kim, Seung-Jean, Kwangmoo Koh, Stephen Boyd, and Dimitry Gorinevsky. 2009. “\(l_1\) Trend Filtering.” SIAM Review 51 (2): 339–60. https://doi.org/doi:10.1137/070690274.