A package for L0-regularized sparse approximation

Description

Trend filtering is a typical method for nonparametric regression. The commonly used trend filtering models is the L1 trend filtering model (a) based on the difference matrix \boldsymbol{D}^{(q+1)}, as illustrated below.

\min _{\boldsymbol{\alpha} \in \mathbb{R}^n} \frac{1}{2}\|\boldsymbol{y}-\boldsymbol{\alpha}\|_2^2 + \lambda\|\boldsymbol{D}^{(q+1)} \boldsymbol{\alpha}\|_{\ell_1}, \quad q=0,1,2, \ldots. \quad (a)

L0 trend filtering (b) has a advantage over other trend filtering methods, especially in the detection of change points. The expression for L0 trend filtering is as follows:

\min _{\boldsymbol{\alpha} \in \mathbb{R}^n} \frac{1}{2}\|\boldsymbol{y}-\boldsymbol{\alpha}\|_2^2 + \lambda\|\boldsymbol{D}^{(q+1)} \boldsymbol{\alpha}\|_{\ell_0}. \quad (b)

We explore transforming the problem (b) into a L0-regularized sparse format (c) by introducing an artificial design matrix \boldsymbol{X}^{(q+1)} that corresponds to the difference matrix, thereby reformulating the L0 trend filtering problem into the following format.

\min _{\boldsymbol{\beta} \in \mathbb{R}^n} \frac{1}{2}\|\boldsymbol{y}-\boldsymbol{X}^{(q+1)}\boldsymbol{\beta}\|_2^2 + \lambda \sum_{i=q+2}^n |\boldsymbol{\beta}_i|_{\ell_0}. \quad (c)

In our practical approach, we consider the maximum number of change points k_{\text{max}} as a constraint, transforming the aforementioned L0 penalty problem (c) into the following L0 constraint problem.

\text{ minimize }\frac{1}{2}\|\boldsymbol{y}-\boldsymbol{X}^{(q+1)}\boldsymbol{\beta}\|_2^2,\quad \text{ subject to } \sum_{i=q+2}^n |\boldsymbol{\beta}_i|_{\ell_0} \leq k_{\text{max}}. \quad (d)

For such L0 constraint problems (d), we employ a splicing-based approach to design algorithms for processing. This package has the following seven main methods:

matrix with special structure

\quadGenerate \boldsymbol{X}^{(q+1)} or \boldsymbol{D}^{(q+1)} matrix.

inverse of the crossprod matrix

\quadSimplify the calculation of the inverse matrix of (\boldsymbol{X}^{(q+1)}_A)^T \boldsymbol{X}^{(q+1)}_A for the cases where q=0 or q=1, which is frequently used in splicing algorithms.

inverse L0 trend filtering with fixed change points

\quadFit a piecewise constant or piecewise linear estimated trend with a given number of change points.

inverse L0 trend filtering with optimal change points

\quadFit a piecewise constant or piecewise linear estimated trend with a maximum number of change points, and select the optimal estimated trend using appropriate information criteria.

simulated data

\quadGenerate piecewise constant or piecewise linear data.

print/coef

\quadPrint a summary of the trend estimation results.

plot

\quadPlot a summary of the trend estimation results.

Details

_PACKAGE

• In previous studies, algorithms solving trend filtering problems (a) necessitate the computation of ((\boldsymbol{D}^{(q+1)})^T \boldsymbol{D}^{(q+1)})^{-1}. When n is large, just fitting the matrix into memory becomes an issue.

• In L0 trend filtering (b), the positions of non-zero elements in the L0 norm correspond with the locations of change points. We consider two subsets: the active set A for non-zero elements and the inactive set I for zero elements. Despite this, computing ((\boldsymbol{D}^{(q+1)}_I)^T \boldsymbol{D}^{(q+1)}_I)^{-1} remains a task involving a substantial matrix.

• Due to the connection between L0 constraint problems and L0 penalty problems, and considering that the sparsity of \boldsymbol{\beta} is is more meaningful in practical applications than the selection of the hyperparameter \lambda. We focus on the constraint that reflects our aim to achieve an estimated trend with a given number of change points. So we transform the L0 penalty problem (c) into the L0 constraint problem (d).

References

Kim SJ, Koh K, Boyd SP and Gorinevsky DM. L1 Trend Filtering. Society for Industrial and Applied Mathematics (2009).

Wen C, Wang X and Zhang A. L0 Trend Filtering. INFORMS Journal on Computing (2023).

Generate a difference matrix

Description

This function generates a matrix for computing differences of a certain order, useful in numerical methods and for creating specific matrix patterns.

Usage

DiffMat(n, q)

Arguments

Value

A matrix with dimensions n-q-1 by n, whose elements correspond to the combinatorial values of q.

See Also

XMat

Examples

Mat1 <- DiffMat(n = 10, q = 0)
print(Mat1)

Mat2 <- DiffMat(n = 15, q = 1)
print(Mat2)

Mat3 <- DiffMat(n = 15, q = 2)
print(Mat3)

The inverse L0 trend filtering with fixed change points

Description

Fit the input data points to a piecewise constant or piecewise linear trend with a given number of change points.

Usage

L0TFinv.fix(y = y, k = k, q = q, first = 0, last = 1)

Arguments

Value

An S3 object of type "L0TFinvfix". A list containing the fitted trend results:

See Also

L0TFinv.opt

Examples

tau = c(0.1, 0.3, 0.4, 0.7, 0.85)
h = c(-1, 5, 3, 0, -1, 2)
BlocksData <- SimuBlocksInv(n = 350, sigma = 0.2, seed = 50, tau = tau ,h = h)
res <- L0TFinv.fix(y=BlocksData$y, k=5, q=0, first=0.01, last=1)
print(res$Ak)
print(BlocksData$setA)
plot(BlocksData$x, BlocksData$y, xlab="", ylab="")
lines(BlocksData$x, BlocksData$y0, col = "red")
lines(BlocksData$x, res$yk, col = "lightgreen")

tau1 = c(0.1, 0.3, 0.4, 0.7, 0.85)
h1 = c(-1, 5, 3, 0, -1, 2)
a0 = -10
WaveData <- SimuWaveInv(n = 2000, sigma = 0.1, seed = 50, tau = tau1, h = h1, a0 = a0)
res1 <- L0TFinv.fix(y=WaveData$y, k=5, q=1, first=0, last=0.99)
print(res1$Ak)
print(WaveData$setA)
plot(WaveData$x, WaveData$y, xlab="", ylab="")
lines(WaveData$x, WaveData$y0, col = "red")
lines(WaveData$x, res1$yk, col = "lightgreen")

The inverse L0 trend filtering with optimal change points

Description

Fit the input data points to a piecewise constant or piecewise linear trend with optimal change points.

Usage

L0TFinv.opt(y = y, kmax = kmax, q = q, first = 0, last = 1, penalty = "bic")

Arguments

Details

Let the fitted trend be denoted as \hat{\boldsymbol{y}}, then

\text{sic} = n \times \log(\frac{1}{n}\|\boldsymbol{y}-\hat{\boldsymbol{y}}\|_2^2)+2\log(\log(n)) \times \log(n) \times \text{df}(\hat{\boldsymbol{y}})

and

\text{bic} = n \times \log(\frac{1}{n}\|\boldsymbol{y}-\hat{\boldsymbol{y}}\|_2^2)+2\times \log(n) \times \text{df}(\hat{\boldsymbol{y}}).

The term \text{df}(\hat{\boldsymbol{y}}) represents the degrees of freedom for the estimated trend, where \text{df}(\hat{\boldsymbol{y}})=k+q+1. Here, k refers to the number of change points in the estimated trend.

Value

An S3 object of type "L0TFinvopt". A list containing the fitted trend results:

Examples

tau = c(0.1, 0.3, 0.4, 0.7, 0.85)
h = c(-1, 5, 3, 0, -1, 2)
BlocksData <- SimuBlocksInv(n = 500, sigma = 0.2, seed = 50, tau = tau ,h = h)
res <- L0TFinv.opt(y=BlocksData$y, kmax=20, q=0, first=0.01, last=1, penalty="bic")
print(res$Aopt)
print(BlocksData$setA)
plot(BlocksData$x, BlocksData$y, xlab="", ylab="")
lines(BlocksData$x, BlocksData$y0, col = "red")
lines(BlocksData$x, res$yopt, col = "lightgreen")

tau1 = c(0.4, 0.6, 0.7)
h1 = c(-3, 5, -4, 6)
a0 = -10
WaveData <- SimuWaveInv(n = 500, sigma = 0.1, seed = 50, tau = tau1, h = h1, a0 = a0)
res1 <- L0TFinv.opt(y=WaveData$y, kmax=10, q=1, first=0, last=0.99, penalty="sic")
print(res1$Aopt)
print(WaveData$setA)
plot(WaveData$x, WaveData$y, xlab="", ylab="")
lines(WaveData$x, WaveData$y0, col = "red")
lines(WaveData$x, res1$yopt, col = "lightgreen")

Simulate Blocks Data

Description

This function generates data points of piecewise constant trends.

Usage

SimuBlocksInv(n, sigma, seed = NA, tau, h)

Arguments

Details

• To simplify the analysis, normalize the change point positions to a range between 0 and 1. Require that all elements of the input tau are within this range. Consequently, the change point positions in simulated data forms a subset of the set { \frac{1}{n}, \frac{2}{n}, \frac{3}{n},..., 1}.

• In fact, length(tau) change points can divide the interval into length(tau)+1 segments of constant function values. Therefore, ensure that the length of vector h is length(tau)+1.

Value

A list containing the piecewise constant simulated data and the underlying trend:

Examples

tau = c(0.1, 0.3, 0.4, 0.7, 0.85)
h = c(-1, 5, 3, 0, -1, 2)
BlocksData <- SimuBlocksInv(n = 350, sigma = 0.1, seed = 50, tau = tau ,h = h)
plot(BlocksData$x, BlocksData$y, xlab="", ylab="")
lines(BlocksData$x, BlocksData$y0, col = "red")
print(BlocksData$setA)
print(BlocksData$tau)

Simulate Wave Data

Description

This function generates data points of piecewise linear trends.

Usage

SimuWaveInv(n, sigma, seed = NA, tau, h, a0 = 0)

Arguments

Value

A list containing the piecewise linear simulated data and the underlying trend:

See Also

SimuBlocksInv

Examples

tau = c(0.1, 0.3, 0.4, 0.7, 0.85)
h = c(-1, 5, 3, 0, -1, 2)
a0 = -10
WaveData <- SimuWaveInv(n = 650, sigma = 0.1, seed = 50, tau = tau, h = h, a0 = a0)
plot(WaveData$x, WaveData$y, xlab="", ylab="")
lines(WaveData$x, WaveData$y0, col = "red")
print(WaveData$setA)
print(WaveData$tau)

Print four metrics about change point detection results

Description

Prints four metrics to compare the quality of change point detection results.

Usage

TFmetrics(y0, tau = NULL, yhat, cpts = NULL)

Arguments

Details

\hat{\boldsymbol{\tau}} represents the estimated change point positions, while \boldsymbol{\tau} denotes the locations of change points in the underlying trend.

d_H=\frac{1}{n} \max \{\max_k \min_j |\tau_j-\hat{\tau}_k|,\max_j \min_k |\tau_j-\hat{\tau}_k|\}.

Note that the number of \hat{\boldsymbol{\tau}} and \boldsymbol{\tau} does not need to be the same.

Value

Examples

tau = c(0.1, 0.3, 0.4, 0.7, 0.85)
h = c(-1, 5, 3, 0, -1, 2)
n = 500
BlocksData <- SimuBlocksInv(n = n, sigma = 0.2, seed = 50, tau = tau ,h = h)
res <- L0TFinv.opt(y=BlocksData$y, kmax=10, q=0, first=0.01, last=1, penalty="bic")
metrics <- TFmetrics(BlocksData$y0,BlocksData$tau,res$yopt,res$Aopt/n)
print(metrics)

tau1 = c(0.1, 0.3, 0.4, 0.7, 0.85)
h1 = c(-1, 5, 3, 0, -1, 2)
a0 = -10
n1 = 2000
WaveData <- SimuWaveInv(n = n1, sigma = 0.1, seed = 50, tau = tau1, h = h1, a0 = a0)
res1 <- L0TFinv.fix(y=WaveData$y, k=20, q=1, first=0, last=0.99)
metrics1 <- TFmetrics(WaveData$y0,WaveData$tau,res1$y.all[,5],res1$A.all[[5]]/n1)
print(metrics1)

Generate an artificial design matrix

Description

This matrix corresponds to the difference matrix, transforming the L0 trend filtering model into an inverse statistical problem.

Usage

XMat(n, q)

Arguments

Details

Noticing the correspondence between \boldsymbol{D}^{(q+1)} and \boldsymbol{X}^{(q+1)}, the result of their matrix multiplication is a combination of a zero matrix and an identity matrix. Expressed as \boldsymbol{D}^{(q+1)} \boldsymbol{X}^{(q+1)}=(\boldsymbol{O}_{(n-q-1)\times(q+1)},\quad \boldsymbol{I}_{(n-q-1)\times(n-q-1)}). The result is advantageous for the invertible processing of the original L0 trend filtering problem.

Value

A matrix with dimensions n by n, whose elements correspond to the difference matrix.

Examples

mat1 <- XMat(n = 10, q = 0)
print(mat1)

mat2 <- XMat(n = 15, q = 1)
print(mat2)

mat3 <- XMat(n = 15, q = 2)
print(mat3)

Mat1 <- DiffMat(n = 10, q = 0)
Mat2 <- DiffMat(n = 15, q = 1)
Mat3 <- DiffMat(n = 15, q = 2)
print(Mat1%*%mat1)
print(Mat2%*%mat2)
print(Mat3%*%mat3)

Extract estimated trends

Description

Extract the coefficients of the estimated trends under the constraint of a given number of change points.

Usage

## S3 method for class 'L0TFinvfix'
coef(object, k = NULL, ...)

## S3 method for class 'L0TFinvopt'
coef(object, k = NULL, ...)

Arguments

Value

Estimated parameters and trends for L0TFinvfix or L0TFinvopt models with specified numbers of change points

See Also

print.L0TFinvfix

Plot L0TFinvfix or L0TFinvopt object

Description

Plots a summary of L0TFinvfix or L0TFinvopt

Usage

## S3 method for class 'L0TFinvfix'
plot(x, type = NULL, k = NULL, ...)

## S3 method for class 'L0TFinvopt'
plot(x, type = "mse", k = NULL, ...)

Arguments

Value

Plot the information criteria for L0TFinvfix and L0TFinvopt across an increasing number of change points or display the estimated trend with a selected change point

See Also

print.L0TFinvfix

Print L0TFinvfix or L0TFinvopt object

Description

Prints a summary of L0TFinvfix or L0TFinvopt

Usage

## S3 method for class 'L0TFinvfix'
print(x, ...)

## S3 method for class 'L0TFinvopt'
print(x, ...)

Arguments

Value

Estimated parameters and information criteria for L0TFinvfix and L0TFinvopt across an increasing number of change points

Examples

library(ggplot2)

tau = c(0.1, 0.3, 0.4, 0.7, 0.85)
h = c(-1, 5, 3, 0, -1, 2)
BlocksData <- SimuBlocksInv(n = 500, sigma = 0.2, seed = 50, tau = tau ,h = h)
res <- L0TFinv.opt(y=BlocksData$y, kmax=10, q=0, first=0.01, last=1, penalty="bic")
print(res)
coef(res,k=res$kopt)
plot(res,type="yhat")
plot(res,type="bic")

tau1 = c(0.1, 0.3, 0.4, 0.7, 0.85)
h1 = c(-1, 5, 3, 0, -1, 2)
a0 = -10
WaveData <- SimuWaveInv(n = 2000, sigma = 0.1, seed = 50, tau = tau1, h = h1, a0 = a0)
res1 <- L0TFinv.fix(y=WaveData$y, k=20, q=1, first=0, last=0.99)
print(res1)
coef(res1,k=5)
plot(res1,type="yhat",k=5)
plot(res1,type="mse")

Generate the inverse of the crossprod matrix

Description

Generate the inverse matrix of (\boldsymbol{X}^{(q+1)}_A)^T \boldsymbol{X}^{(q+1)}_A for the cases where q=0 or q=1, commonly employed in splicing algorithms. Note that an explicit solution exists for the inverse when q=0, but not when q=1.

Usage

solMat(n, q, A)

Arguments

Value

The inverse matrix of (\boldsymbol{X}^{(q+1)}_A)^T \boldsymbol{X}^{(q+1)}_A for the cases where q= 0 or 1.

Examples

Mat1 <- XMat(n = 10, q = 0)
A1 = c(1,2,5,8)
mat1 = as.matrix(Mat1[,A1])
S1 <- solMat(n = 10, q = 0, A = A1)
print(S1)
print(round(S1%*%t(mat1)%*%mat1,10))

Mat2 <- XMat(n = 15, q = 1)
A2 = c(1,3,8,10,15)
mat2 = as.matrix(Mat2[,A2])
S2 <- solMat(n = 15, q = 1, A = A2)
print(S2)
print(round(S2%*%t(mat2)%*%mat2,10))