Description Usage Arguments Details Value Author(s) References See Also Examples

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.

1 2 3 |

`y` |
Observed data, of length n. |

`D` |
Penalty matrix of size m x n (see Details). Default is |

`D_type` |
Types of |

`q` |
Nonnegative integer used to specify the order of penalty matrix. See |

`kmax` |
The maximum number of knots. The number of knots is thus |

`rho` |
The hyperparameter |

`tmax` |
The maximum number of iterations in the AMIAS algorithm. |

`eps` |
The tolerance |

`smooth` |
Whether to smooth the data, if |

`h` |
Bandwidth in smoothing data. See |

`adjust` |
Whether to adjust the indexes of the active set in the AMIAS algorithm. If |

`delta` |
The minimum gap etween the adjacent knots. Only used when |

`adjust.max` |
The number of iterations in the adjustment step. Only used when |

`...` |
Other arguments. |

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*.

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 |

`q` |
The order of penalty matrix. |

`k` |
The sequential list of number of knots, i.e., |

`alpha` |
The optimal coefficients |

`v` |
The optimal primal variable or split variable of the argumented lagrange form in |

`u` |
The optimal dual variable or lagrange operator of the argumented lagrange form in |

`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 |

`v.all` |
Matrix of primal variables with |

`u.all` |
Matrix of dual variables with |

`A.all` |
List of active sets, of length |

`bic` |
The BIC value for each value of |

`iter` |
A vector of length |

`smooth` |
Whether to smooth the data. |

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

Maintainer: Canhong Wen <wench@ustc.edu.cn>

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)
``` |

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.