# 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

 1 2 3 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 <wench@ustc.edu.cn>

## References

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

print.samias, coef.samias and plot.samias method, and the amias function.
  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50  ##----- 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)