f2sSP: Overlap Group Least Absolute Shrinkage and Selection Operator for scalar-on-function regression model

Description

Overlap Group-LASSO for scalar-on-function regression model solves the following optimization problem $$\textrm{min}_{\psi,\gamma} ~ \frac{1}{2} \sum_{i=1}^n \left( y_i - \int x_i(t) \psi(t) dt-z_i^\intercal\gamma \right)^2 + \lambda \sum_{g=1}^{G} \Vert S_{g}T\psi \Vert_2$$ to obtain a sparse coefficient vector $\psi\in\mathbb{R}^{M}$ for the functional penalized predictor $x(t)$ and a coefficient vector $\gamma\in\mathbb{R}^q$ for the unpenalized scalar predictors $z_1,\dots,z_q$. The regression function is $\psi(t)=\varphi(t)^\intercal\psi$ where $\varphi(t)$ is a B-spline basis of order $d$ and dimension $M$. For each group $g$, each row of the matrix $S_g\in\mathbb{R}^{d\times M}$ has non-zero entries only for those bases belonging to that group. These values are provided by the arguments groups and group_weights (see below). Each basis function belongs to more than one group. The diagonal matrix $T\in\mathbb{R}^{M\times M}$ contains the basis-specific weights. These values are provided by the argument var_weights (see below). The regularization path is computed for the overlap group-LASSO penalty at a grid of values for the regularization parameter $\lambda$ using the alternating direction method of multipliers (ADMM). See Boyd et al. (2011) and Lin et al. (2022) for details on the ADMM method.

Usage

f2sSP(
  vY,
  mX,
  mZ = NULL,
  M,
  group_weights = NULL,
  var_weights = NULL,
  standardize.data = TRUE,
  splOrd = 4,
  lambda = NULL,
  nlambda = 30,
  lambda.min.ratio = NULL,
  intercept = FALSE,
  overall.group = FALSE,
  control = list()
)

Value

A named list containing

sp.coefficients: a length-$M$ solution vector for the parameters $\psi$, which corresponds to the minimum in-sample MSE.
sp.coef.path: an $(n_\lambda\times M)$ matrix of estimated $\psi$ coefficients for each lambda.
sp.fun: a length-$r$ vector providing the estimated functional coefficient for $\psi(t)$.
sp.fun.path: an $(n_\lambda\times r)$ matrix providing the estimated functional coefficients for $\psi(t)$ for each lambda.
coefficients: a length-$q$ solution vector for the parameters $\gamma$, which corresponds to the minimum in-sample MSE. It is provided only when either the matrix $Z$ in input is not NULL or the intercept is set to TRUE.
coef.path: an $(n_\lambda\times q)$ matrix of estimated $\gamma$ coefficients for each lambda. It is provided only when either the matrix $Z$ in input is not NULL or the intercept is set to TRUE.
lambda: sequence of lambda.
lambda.min: value of lambda that attains the minimum in-sample MSE.
mse: in-sample mean squared error.
min.mse: minimum value of the in-sample MSE for the sequence of lambda.
convergence: logical. 1 denotes achieved convergence.
elapsedTime: elapsed time in seconds.
iternum: number of iterations.

When you run the algorithm, output returns not only the solution, but also the iteration history recording following fields over iterates,

objval: objective function value.
r_norm: norm of primal residual.
s_norm: norm of dual residual.
eps_pri: feasibility tolerance for primal feasibility condition.
eps_dual: feasibility tolerance for dual feasibility condition.

Iteration stops when both r_norm and s_norm values become smaller than eps_pri and eps_dual, respectively.

Arguments

vY: a length-$n$ vector of observations of the scalar response variable.
mX: a $(n\times r)$ matrix of observations of the functional covariate.
mZ: an $(n\times q)$ full column rank matrix of scalar predictors that are not penalized.
M: number of elements of the B-spline basis vector $\varphi(t)$.
group_weights: a vector of length $G$ containing group-specific weights. The default is square root of the group cardinality, see Bernardi et al. (2022).
var_weights: a vector of length $M$ containing basis-specific weights. The default is a vector where each entry is the reciprocal of the number of groups including that basis. See Bernardi et al. (2022) for details.
standardize.data: logical. Should data be standardized?
splOrd: the order $d$ of the spline basis.
lambda: either a regularization parameter or a vector of regularization parameters. In this latter case the routine computes the whole path. If it is NULL values for lambda are provided by the routine.
nlambda: the number of lambda values - default is 30.
lambda.min.ratio: smallest value for lambda, as a fraction of the maximum lambda value. If $n>M$, the default is 0.0001, and if $n<M$, the default is 0.01.
intercept: logical. If it is TRUE, a column of ones is added to the design matrix.
overall.group: logical. If it is TRUE, an overall group including all penalized covariates is added.
control: a list of control parameters for the ADMM algorithm. See ‘Details’.

Details

The control argument is a list that can supply any of the following components:

adaptation: logical. If it is TRUE, ADMM with adaptation is performed. The default value is TRUE. See Boyd et al. (2011) for details.
rho: an augmented Lagrangian parameter. The default value is 1.
tau.ada: an adaptation parameter greater than one. Only needed if adaptation = TRUE. The default value is 2. See Boyd et al. (2011) and Lin et al. (2022) for details.
mu.ada: an adaptation parameter greater than one. Only needed if adaptation = TRUE. The default value is 10. See Boyd et al. (2011) and Lin et al. (2022) for details.
abstol: absolute tolerance stopping criterion. The default value is sqrt(sqrt(.Machine$double.eps)).
reltol: relative tolerance stopping criterion. The default value is sqrt(.Machine$double.eps).
maxit: maximum number of iterations. The default value is 100.
print.out: logical. If it is TRUE, a message about the procedure is printed. The default value is TRUE.

References

bernardi_etal.2022fdaSP

boyd_etal.2011fdaSP

jenatton_etal.2011fdaSP

lin_etal.2022fdaSP

Examples

Run this code


## generate sample data
set.seed(1)
n     <- 40
p     <- 18                                  # number of basis to GENERATE beta
r     <- 100
s     <- seq(0, 1, length.out = r)

beta_basis <- splines::bs(s, df = p, intercept = TRUE)    # basis
coef_data  <- matrix(rnorm(n*floor(p/2)), n, floor(p/2))        
fun_data   <- coef_data %*% t(splines::bs(s, df = floor(p/2), intercept = TRUE))     

x_0   <- apply(matrix(rnorm(p, sd=1),p,1), 1, fdaSP::softhresh, 1)  # regression coefficients 
x_fun <- beta_basis %*% x_0                

b     <- fun_data %*% x_fun + rnorm(n, sd = sqrt(crossprod(fun_data %*% x_fun ))/10)
l     <- 10^seq(2, -4, length.out = 30)
maxit <- 1000


## set the hyper-parameters
maxit          <- 1000
rho_adaptation <- TRUE
rho            <- 1
reltol         <- 1e-5
abstol         <- 1e-5

mod <- f2sSP(vY = b, mX = fun_data, M = p,
             group_weights = NULL, var_weights = NULL, standardize.data = FALSE, splOrd = 4,
             lambda = NULL, nlambda = 30, lambda.min = NULL, overall.group = FALSE, 
             control = list("abstol" = abstol, 
                            "reltol" = reltol, 
                            "adaptation" = rho_adaptation, 
                            "rho" = rho, 
                            "print.out" = FALSE)) 

# plot coefficiente path
matplot(log(mod$lambda), mod$sp.coef.path, type = "l", 
        xlab = latex2exp::TeX("$\\log(\\lambda)$"), ylab = "", bty = "n", lwd = 1.2)

Run the code above in your browser using DataLab