samias: Solve Generalized l_0 Problem with SAMIAS Method

In AMIAS: Alternating Minimization Induced Active Set Algorithms

Description

This function solves the generalized l_0 problem with a sequential list of numbers of knots. The generalized l_0 coefficient is computed via the sequential alternating minimization induced active set (SAMIAS) algorithm. Can deal with any polynomial order or even a general penalty matrix for structural filtering.

Usage

samias(y, D = NULL, D_type = c("tf0", "tfq", "user"), q = 0, kmax = min(20,m-1), 
     rho = n^(q+1), tmax = 10, eps = 0, smooth = FALSE, h = 5, 
     adjust = FALSE, delta = 10, adjust.max = 10, ...)

Arguments

y	Observed data, of length n.
D	Penalty matrix of size m x n (see Details). Default is NULL.
D_type	Types of D. Either "tf0", "tfq" or "user", depending on what types of contraint that user wants to impose on the data. For D_type = "tfq" and D = NULL, we solve the l_0 trend filtering of order q+1, where q is determined by the argument q. For D_type = "tf0" and D = NULL, we solve the l_0 trend filtering of order 0, i.e., with piecewise constant constraint on the data. For D_type = "user", the penalty matrix D must be specified.
q	Nonnegative integer used to specify the order of penalty matrix. See genDtf1d for Details.
kmax	The maximum number of knots. The number of knots is thus {1,2,…,kmax}. Default is min(20,m-1).
rho	The hyperparameter ρ in the augmented Lagrangian of the l_0 trend filtering problem. Default is n^{q+1}.
tmax	The maximum number of iterations in the AMIAS algorithm.
eps	The tolerance ε for an early stoping rule in the SAMIAS algorithm. The early stopping rule is defined as \|y- α\|_2^2/n ≤ ε, where α is the fitting coefficient.
smooth	Whether to smooth the data, if TRUE, it smoothes the input data. Default is FALSE.
h	Bandwidth in smoothing data. See my.rollmean for details.
adjust	Whether to adjust the indexes of the active set in the AMIAS algorithm. If TRUE, it implements the adjustment when the indexes in the active set are not well separated. Default is FALSE.
delta	The minimum gap etween the adjacent knots. Only used when adjust = TRUE.
adjust.max	The number of iterations in the adjustment step. Only used when adjust = TRUE.
...	Other arguments.

Details

The generalized l_0 problem with a maximal number of possible knots kmax is

\min_α 1/2 \|y-α\|_2^2 \;\;{\rm s.t.}\;\; \|Dα\|_0 ≤ kmax.

The penalty matrix D is either the discrete difference operator of order q + 1 or a general matrix specified by users to enforce some special structure. We solve this problem by sequentially considering the following sub-problem.

\min_α 1/2 \|y-α\|_2^2 \;\;{\rm s.t.}\;\; \|Dα\|_0 = k

for k ranging from 1 to kmax. This sub-problem can be solved via the AMIAS method, see amias. Thus a sequential AMIAS algorithm named SAMIAS is proposed, which is a simple extension of AMIAS by adopting the warm start strategy.

Since it outputs the whole sequence of results for k = 1, . . . , kmax, we can naturally draw the solution path for increasing k and apply the Bayesian information criterion (BIC) for hyperparameter determination. The BIC is defined by

n log(mse) + 2*log(n)*df,

where mse is the mean squared error, i.e., \|y-α\|_2^2/n, and df = k+q+1.

Value

A list with class attribute 'samias' and named components:

call	The call that produces this object.
y	Observe sequence, if smooth, the smooth y will be returned.
D_type	Types of D.
q	The order of penalty matrix.
k	The sequential list of number of knots, i.e., {1,2,…,kmax}
alpha	The optimal coefficients α determined by minimizing BIC.
v	The optimal primal variable or split variable of the argumented lagrange form in Dα.
u	The optimal dual variable or lagrange operator of the argumented lagrange form in Dα for linear item.
A	The optimal estimate of active sets, i.e., set of detected knots.
df	Degree of freedom for all the candidate models.
alpha.all	Coefficients Matrix with k ranging from 1 to kmax, of dimension n x kmax.
v.all	Matrix of primal variables with k ranging from 1 to kmax, of dimension (n-q-1) x kmax.
u.all	Matrix of dual variables with k ranging from 1 to kmax, of dimension (n-q-1) x kmax.
A.all	List of active sets, of length kmax.
bic	The BIC value for each value of k.
iter	A vector of length kmax, record the iterations in the AMIAS algorithm.
smooth	Whether to smooth the data.

Author(s)

Canhong Wen, Xueqin Wang, Shijie Quan, Zelin Hong and Aijun Zhang.

Maintainer: Canhong Wen <[e-mail]>

References

Wen, C., Zhu, J., Wang, X., and Zhang, A. (2019) L0 trend filtering, technique report.

See Also

print.samias, coef.samias and plot.samias method, and the amias function.

Examples

  ##----- A toy example of piecewise constant signal -------
  set.seed(0)
  n <- 100
  x = seq(1/n, 1,length.out = n)
  y0 = 0*x; y0[x>0.5] = 1
  y = y0 + rnorm(n, sd = 0.1)
  fit <- samias(y, kmax = 5) 
  op <- par(mfrow=c(1,2))
  plot(fit, type= "coef", add.knots = FALSE, main = "Piecewise Constant")
  plot(fit, type = "vpath", main = "Piecewise Constant")
  par(op)
  
  
  ##----- A toy example of piecewise linear trend -------
  set.seed(0)
  y0 = 2*(0.5-x); y0[x>0.5] = 2*(x[x>0.5]-0.5)
  y = y0 + rnorm(n, sd = 0.1)
  fit <- samias(y, D_type = "tfq", q = 1, kmax = 5) 
  print(fit)
  

  ##------ Piecewise constant trend filtering example in Wen et al.(2018).-----
  set.seed(1)
  data <- SimuBlocks(2048)
  fit <- samias(data$y, kmax = 15)           # With default input argument
  plot(fit, type="coef", main = "Blocks")    # Plot the optimal estimate
  lines(data$x, data$y0, type="s")           # Add the true signal for reference
  plot(fit, type= "vpath", main = "Blocks")  # Plot the solution path

  
  ##------ Piecewise linear trend filtering example in Wen et al.(2018).-----
  set.seed(1)
  data <- SimuWave(512)
  
  # samias With adjustment
  fit <- samias(data$y, q = 1, D_type = "tfq", kmax = 15, adjust = TRUE, delta = 20) 
  plot(fit, main = "Wave", k = 10)          # Plot the estimate with user-specified k
  lines(data$x, data$y0, type="l")          # Add the true signal for reference
  plot(fit, type= "vpath", main = "Wave")   # Plot the solution path
  
  
  ##------ Doppler example in Wen et al.(2018).-----
  set.seed(1)
  data <- SimuDoppler(1024)
  op <- par(mfrow=c(1,2))
  fit1 <- samias(data$y, q = 2, D_type = "tfq", kmax = 30) # piecewise quadratic
  plot(fit1, main = "Doppler: q = 2")
  fit2 <- samias(data$y, q = 3, D_type = "tfq", kmax = 30) # piecewise Cubic polynomial
  plot(fit2, main = "Doppler: q = 3")
  par(op)

AMIAS documentation built on May 2, 2019, 2:10 a.m.