diff --git a/.Rbuildignore b/.Rbuildignore index 3b3ecc8..9f6bb39 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -12,3 +12,4 @@ README.qmd ^memcheck ^\.positai$ ^\.claude$ +^vignettes/*_files$ diff --git a/.gitignore b/.gitignore index 3a34ba2..2461146 100644 --- a/.gitignore +++ b/.gitignore @@ -62,3 +62,4 @@ docs .vscode/ tools/ .positai +inst/doc diff --git a/DESCRIPTION b/DESCRIPTION index d2cddef..ee69836 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -5,10 +5,10 @@ Version: 1.0-0 Date: 2026-05-xx Authors@R: person("Julien", "Chiquet", , "julien.chiquet@inrae.fr", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-3629-3429")) -Description: This package fits the solution paths of classical sparse regression - models with efficient active set algorithms by solving small sub-problems. - Include LASSO, SCAD, MCP, (Sparse) Group-LASSO, Cooperative-LASSO, LAVA and - Elastic-Net variants. Also provides a few methods for model selection purpose +Description: This package fits the solution paths of classical sparse regression + models with efficient active set algorithms by solving small sub-problems. + Include LASSO, SCAD, MCP, (Sparse) Group-LASSO, Cooperative-LASSO, LAVA and + Elastic-Net variants. Also provides a few methods for model selection purpose (cross-validation, stability selection). URL: https://github.com/jchiquet/quadrupen, https://jchiquet.github.io/quadrupen/ BugReports: https://github.com/jchiquet/quadrupen/issues @@ -36,11 +36,14 @@ Suggests: elasticnet, glmnet, igraph, + knitr, lars, MASS, ncvreg, + grpreg, + rmarkdown, testthat -LinkingTo: +LinkingTo: Rcpp, RcppArmadillo Collate: @@ -67,3 +70,4 @@ Collate: Packaged: 2013-03-25 19:41:00 UTC; jchiquet Roxygen: list(markdown = TRUE) Config/roxygen2/version: 8.0.0 +VignetteBuilder: knitr diff --git a/NEWS.md b/NEWS.md index 54950bb..4b06e73 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,6 @@ NEWS/Changelog -# quadrupen dev version (2026-05-xx) +# quadrupen dev version (2026-06-xx) - major updates - complete rewriting of R code using R6 classes @@ -13,6 +13,7 @@ NEWS/Changelog - added mcp and scad (concave penalties, Chernozukov et al, 2017) - added refit version of Lasso/Elastic-Net ("relaxed" Lasso/Enet) - changing many parameters (badly) named, do not expect backward compatibility + - added vignettes - minor updates - Integration of changes from CRAN versions from 0.2-4 to 0.2-13 - set up github workflow for pkgdown page diff --git a/R/CrossValidation-R6Class.R b/R/CrossValidation-R6Class.R index fa0c94a..fa63a1d 100644 --- a/R/CrossValidation-R6Class.R +++ b/R/CrossValidation-R6Class.R @@ -77,9 +77,10 @@ CrossValidation <- R6::R6Class( print = function() { self$show() }, #' @description Plot 1-dimensional cross-validation #' @param log_scale logical, should a log-scale be used for the x-axis + #' @param se logical, should confidence band be displayed (TRUE by default) #' @param title graph title #' @return a ggplot object - plotCV_1D = function(log_scale=TRUE, title = "Cross-validation error") { + plotCV_1D = function(log_scale=TRUE, title = "Cross-validation error", se = TRUE) { ## SIMPLE CROSS-VALIDATION GRAPH if (length(self$lambda1) > length(self$lambda2)) { d <- ggplot(self$data, aes(x=.data$lambda1,y=.data$mean)) + theme_bw() @@ -97,7 +98,7 @@ CrossValidation <- R6::R6Class( d <- d + xlab(xlab) + ylab("Mean square error") + geom_point(alpha=0.3) + geom_smooth(aes(ymin=.data$mean-.data$serr, ymax=.data$mean+.data$serr, group=as.factor(.data$lambda2), - colour=as.factor(.data$lambda2)), data=self$data, alpha=0.2, stat="identity") + colour=as.factor(.data$lambda2)), se=se, data=self$data, alpha=0.2, stat="identity") if (log_scale) { d <- d + scale_x_log10() + annotation_logticks(sides="b") } diff --git a/R/QuadrupenFit-R6Class.R b/R/QuadrupenFit-R6Class.R index 545b446..0aa96db 100644 --- a/R/QuadrupenFit-R6Class.R +++ b/R/QuadrupenFit-R6Class.R @@ -540,7 +540,7 @@ QuadrupenFit <- R6::R6Class( "path" = self$plot_path(log_scale = log_scale, labels = labels), "criteria" = self$information_criteria$plot(log_scale = log_scale), "crossval" = self$cross_validation$plot(log_scale = log_scale), - "stability" = self$stability_path$plot(log_scale = log_scale, labels = labels) + "stability" = self$stability_path$plot(labels = labels) ) d diff --git a/R/StabilityPath-R6Class.R b/R/StabilityPath-R6Class.R index 55a6b4e..a9fff26 100644 --- a/R/StabilityPath-R6Class.R +++ b/R/StabilityPath-R6Class.R @@ -91,8 +91,6 @@ StabilityPath <- R6::R6Class( #' (first penalty level) or `"fraction"` (fraction of the #' penalty level applied tune by \eqn{\lambda_1}{lambda1}). Default #' is `"lambda"`. - #' @param log_scale logical; indicates if a log-scale should be used - #' when \code{xvar="lambda"}. Default is `TRUE`. #' @param plot logical; indicates if the graph should be #' plotted. Default is `TRUE`. If `FALSE`, only the #' \pkg{ggplot2} object is sent back. @@ -135,7 +133,7 @@ StabilityPath <- R6::R6Class( #' plot = function( xvar = "lambda", title = "Stability path", labels = rep("unknown status",self$nvar), - sel_mode = c("rank", "PFER"), cutoff=0.75, PFER=2, nvarsel=floor(self$nobs/log(self$nvar))) { + sel_mode = c("rank", "PFER"), cutoff=0.75, PFER=2, nvarsel=min(self$nvar, floor(self$nobs/log(self$nvar)))) { sel_mode <- match.arg(sel_mode) selection <- rep("unselected",length(self$nonzero)) diff --git a/_pkgdown.yml b/_pkgdown.yml index efb0193..10c33a0 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -4,6 +4,18 @@ template: mathjax: true math-rendering: mathjax +articles: +- title: "Getting started" + navbar: ~ + contents: + - sparse-regression + - group-sparse-regression +- title: "Advanced topics" + navbar: ~ + contents: + - structured-signal-recovery + - lava + reference: - title: "Main functions for user" desc: > diff --git a/inst/examples/ridge.R b/inst/examples/ridge.R index 7606b49..b87c690 100644 --- a/inst/examples/ridge.R +++ b/inst/examples/ridge.R @@ -1,4 +1,5 @@ library(quadrupen) +require(Matrix) beta <- rep(c(0,1,0,-1,0), c(25,10,25,10,25)) cor <- 0.75 @@ -10,13 +11,11 @@ n <- 50 x <- as.matrix(matrix(rnorm(95*n),n,95) %*% chol(Sigma)) y <- 10 + x %*% beta + rnorm(n,0,10) p <- ncol(x) +plot(beta, type = "s") - -## neighborhood prior -C <- bandSparse(p,k=0:1,diagonals=list(rep(1,p),rep(-1,p-1))) -L1 <- t(C) %*% C - -## clique prior +## Structuring matrix : +# Laplacian of a community graph +# (clique prior) Ioo <- matrix(1,25,25) Iww <- matrix(1,10,10) A <- bdiag(Ioo,Iww,Ioo,Iww,Ioo) @@ -24,25 +23,39 @@ diag(A) <- 0 L2 <- -A diag(L2) <- colSums(A) + 1e-2 -ridge.classical <- ridge(x,y,intercept=F) -plot(ridge.classical) -ridge.neighborhood <- ridge(x,y, struct=L1, lambda_max=1000) -plot(ridge.neighborhood) -ridge.block <- ridge(x,y, struct=L2, lambda_max=1000) -plot(ridge.block) - - +labels <- paste("segment", rep(1:5, c(25,10,25,10,25)), sep="-" ) +lasso_fit <- lasso(x,y,intercept=F) +plot(lasso_fit, labels=labels) +fused_lasso_fit <- fused_lasso(x,y,lambda2=5,intercept=F) +plot(fused_lasso_fit, labels=labels) +ridge_fit <- ridge(x,y,intercept=F) +plot(ridge_fit, labels=labels) +ridge_struct <- ridge(x,y, struct=L2, lambda_max=1000) +plot(ridge_struct, labels=labels) +enet_fit <- elastic_net(x,y, lambda2 = 5, intercept=F) +plot(enet_fit, labels=labels) +enet_struct <- elastic_net(x,y, lambda2 = 5, struct=L2) +plot(enet_struct, labels=labels) +enet_struct$debias <- TRUE +plot(enet_struct, labels=labels) + +beta_lasso <- as.numeric(lasso(x,y,lambda=10)$coefficients) +beta_fused_lasso <- as.numeric(fused_lasso(x,y,lambda1=20, lambda2 = 5)$coefficients) beta_ridge <- as.numeric(ridge(x,y,lambda=1)$coefficients) -beta_ridge_neighbor <- as.numeric(ridge(x,y,lambda=10,struct=L1)$coefficients) -beta_ridge_block <- as.numeric(ridge(x,y,lambda=10,struct=L2)$coefficients) -beta_enet_block <- as.numeric(elastic_net(x,y,lambda2=5,lambda1=20,struct=L2,refit=TRUE)$coefficients) +beta_ridge_struct <- as.numeric(ridge(x,y,lambda=10,struct=L2)$coefficients) +beta_enet_struct <- as.numeric(elastic_net(x,y,lambda2=5,lambda1=20,struct=L2)$coefficients) +beta_enet_struct_refit <- as.numeric(elastic_net(x,y,lambda2=5,lambda1=20,struct=L2,refit=TRUE)$coefficients) -par(mfrow=c(2,2)) +par(mfrow=c(3,2)) +plot(beta_lasso, type="b", ylim=c(-1,1)) +points(beta, pch="+",col="red") +plot(beta_fused_lasso, type="b", ylim=c(-1,1)) +points(beta, pch="+",col="red") plot(beta_ridge, type="b", ylim=c(-1,1)) points(beta, pch="+",col="red") -plot(beta_ridge_neighbor, type="b", ylim=c(-1,1)) +plot(beta_ridge_struct, type="b", ylim=c(-1,1)) points(beta, pch="+",col="red") -plot(beta_ridge_block, type="b", ylim=c(-1,1)) +plot(beta_enet_struct, type="b", ylim=c(-1,1)) points(beta, pch="+",col="red") -plot(beta_enet_block, type="b", ylim=c(-1,1)) +plot(beta_enet_struct_refit, type="b", ylim=c(-1,1)) points(beta, pch="+",col="red") diff --git a/inst/experimental/group_lasso.R b/inst/experimental/group_lasso.R deleted file mode 100644 index 7d491c4..0000000 --- a/inst/experimental/group_lasso.R +++ /dev/null @@ -1,338 +0,0 @@ -#' #' Fit a linear model with group-lasso (either l1/l2 or l1/l-inf) regularization -#' #' -#' #' Adjust a linear model with group-lasso regularization, that is a -#' #' mixture of either a (possibly weighted) -#' #' \eqn{\ell_1/\ell_2}{l1/l2}- or -#' #' \eqn{\ell_1/\ell_\infty}{l1/linf}-norm, and a (possibly -#' #' structured) \eqn{\ell_2}{l2}-norm (ridge-like). The solution path -#' #' is computed at a grid of values for the -#' #' \eqn{\ell_1/\ell_q}{l1/lq}-penalty. See details for the criterion -#' #' optimized. -#' #' -#' #' @param x matrix of features, possibly sparsely encoded -#' #' (experimental). Do NOT include intercept. When normalized os -#' #' \code{TRUE}, coefficients will then be rescaled to the original -#' #' scale. -#' #' -#' #' @param y response vector. -#' #' -#' #' @param group vector of interers indicating group belonging. Must -#' #' match the number fo column in \code{x}. Must be SORTED integers -#' #' starting from 1. -#' #' -#' #' @param type string indicating wether the -#' #' \eqn{\ell_1/\ell_2}{l1/l2} or the -#' #' \eqn{\ell_1/\ell_\infty}{l1/linf} group-Lasso must be fitted. -#' #' -#' #' @param lambda1 sequence of decreasing \eqn{\ell_1}{l1}-penalty -#' #' levels. If \code{NULL} (the default), a vector is generated with -#' #' \code{nlambda1} entries, starting from a guessed level -#' #' \code{lambda1.max} where only the intercept is included, then -#' #' shrunken to \code{min.ratio*lambda1.max}. -#' #' -#' #' @param lambda2 real scalar; tunes the \eqn{\ell_2}{l2} -#' #' penalty. Default is 0.01. Set to 0 to recover the usual -#' #' group-Lasso. -#' #' -#' #' @param penscale vector with real positive values that weight the -#' #' \eqn{\ell_1}{l1}-penalty of each feature. -#' #' -#' #' @param struct matrix structuring the coefficients, possibly -#' #' sparsely encoded. Must be at least positive semidefinite (this is -#' #' checked internally if the \code{checkarg} argument is -#' #' \code{TRUE}). If \code{NULL} (the default), the identity matrix is -#' #' used. See details below. -#' #' -#' #' @param intercept logical; indicates if an intercept should be -#' #' included in the model. Default is \code{TRUE}. -#' #' -#' #' @param normalize logical; indicates if variables should be -#' #' normalized to have unit L2 norm before fitting. Default is -#' #' \code{TRUE}. -#' #' -#' #' @param naive logical; Compute either 'naive' of classic version: -#' #' as defined for the Elastic-net in Zou and Hastie (2006), the -#' #' vector of parameters is rescaled by a coefficient -#' #' \code{(1+lambda2)} when \code{naive} equals \code{FALSE}. No -#' #' rescaling otherwise. Default is \code{FALSE}. -#' #' -#' #' @param nlambda1 integer that indicates the number of values to put -#' #' in the \code{lambda1} vector. Ignored if \code{lambda1} is -#' #' provided. -#' #' -#' #' @param min.ratio minimal value of \eqn{\ell_1}{l1}-part of the -#' #' penalty that will be tried, as a fraction of the maximal -#' #' \code{lambda1} value. A too small value might lead to unstability -#' #' at the end of the solution path corresponding to small -#' #' \code{lambda1}. The default value tries to avoid this, adapting -#' #' to the '\eqn{n1'}, the -#' #' Fenchel duality gap is computed along the algorithm.} -#' #' } -#' #' -#' #' @param checkargs logical; should arguments be checked to -#' #' (hopefully) avoid internal crashes? Default is -#' #' \code{TRUE}. Automatically set to \code{FALSE} when calls are made -#' #' from cross-validation or stability selection procedures. -#' #' -#' #' @return an object with class \code{quadrupen}, see the -#' #' documentation page \code{\linkS4class{quadrupen}} for details. -#' #' -#' #' @note The optimized criterion is the following: \if{latex}{\deqn{% -#' #' \hat{\beta}_{\lambda_1} = \arg \min_{\beta} \frac{1}{2} (y - X -#' #' \beta)^T (y - X \beta) + \lambda_1 \|D \beta \|_{1}, }} -#' #' \if{html}{\out{
βhat -#' #' λ1 = -#' #' argminβ 1/2 RSS(&beta) + λ1 -#' #' | D β |1,
}} -#' #' \if{text}{\deqn{beta.hat(lambda1) = argmin_beta 1/2 -#' #' RSS(beta) + lambda1 |D beta|1,}} where -#' #' \eqn{D}{D} is a diagonal matrix, whose diagonal terms are provided -#' #' as a vector by the \code{penscale} argument. -#' #' -#' #' @seealso See also \code{\linkS4class{quadrupen}}, -#' #' \code{\link{plot.quadrupen}} and \code{\link{crossval}}. -#' #' @name group.lasso -#' #' @rdname group.lasso -#' #' @keywords models, regression -#' #' -#' #' -#' #' @export -#' group.lasso <- function(x, -#' y, -#' group, -#' type = 2, ## could be "inf" or 2, default is 2 -#' lambda1 = NULL, -#' lambda2 = 0.01, -#' penscale = switch(as.character(type), "2" = rep(sqrt(pk), pk), "inf" = rep(1,p)), -#' struct = NULL, -#' intercept = TRUE, -#' normalize = TRUE, -#' naive = FALSE, -#' nlambda1 = ifelse(is.null(lambda1),100,length(lambda1)), -#' min.ratio = ifelse(n<=p,1e-2,1e-4), -#' max.feat = ifelse(lambda2<1e-2,min(n,p),min(4*n,p)), -#' beta0 = rep(0,p), -#' control = list(), -#' checkargs = TRUE) { -#' -#' p <- ncol(x) # problem size -#' n <- nrow(x) # sample size -#' pk <- tabulate(group) # group size -#' -#' ## =================================================== -#' ## CHECKS TO (PARTIALLY) AVOID CRASHES OF THE C++ CODE -#' if (checkargs) { -#' if (is.data.frame(x)) -#' x <- as.matrix(x) -#' if(!inherits(x, c("matrix", "dgCMatrix"))) -#' stop("x has to be of class 'matrix' or 'dgCMatrix'.") -#' if(any(is.na(x))) -#' stop("NA value in x not allowed.") -#' if(!is.numeric(y)) -#' stop("y has to be of type 'numeric'") -#' if(n != length(y)) -#' stop("x and y have not correct dimensions") -#' if(length(penscale) != p) -#' stop("penscale must have ncol(x) entries") -#' if (any(penscale <= 0)) -#' stop("weights in penscale must be positive") -#' if(!inherits(lambda2, "numeric") | length(lambda2) > 1) -#' stop("lambda2 must be a scalar.") -#' if(lambda2 < 0) -#' stop("lambda2 must be a non negative scalar.") -#' if (!is.null(lambda1)) { -#' if(any(lambda1 <= 0)) -#' stop("entries inlambda1 must all be postive.") -#' if(is.unsorted(rev(lambda1))) -#' stop("lambda1 values must be sorted in decreasing order.") -#' if (length(lambda1)>1 & !is.null(beta0)) -#' warning("providing beta0 for a serie of l1 penalties mught be inefficient.") -#' } -#' if(min.ratio < 0) -#' stop("min.ratio must be non negative.") -#' if (!is.null(struct)) { -#' if (ncol(struct) != p | ncol(struct) != p) -#' stop("struct must be a (square) positive semidefinite matrix.") -#' if (any(eigen(struct,only.values=TRUE)$values<0)) -#' stop("struct must be a (square) positive semidefinite matrix.") -#' if(!inherits(struct, "dgCMatrix")) -#' struct <- as(struct, "dgCMatrix") -#' } -#' ## Check the group vector -#' if(length(group) != p) -#' stop("group length does not match ncol(x)") -#' if(is.unsorted(group)) -#' stop("the group vector must be sorted (do the manual permutation of the columns of x please)") -#' if(any(diff(unique(group)) != 1)) -#' stop("the groups must be labeled from 1 to nb groups with increment of 1.") -#' -#' if (!is.null(beta0)) { -#' beta0 <- as.numeric(beta0) -#' if (length(beta0) != p) -#' stop("beta0 must be a vector with p entries.") -#' } -#' if (length(max.feat)>1) -#' stop("max.feat must be an integer.") -#' if(is.numeric(max.feat) & !is.integer(max.feat)) -#' max.feat <- as.integer(max.feat) -#' } -#' -#' ## ============================================ -#' ## RECOVERING LOW LEVEL OPTIONS -#' quadra <- TRUE -#' if (!is.null(control$method)) { -#' if (control$method != "quadra") { -#' quadra <- FALSE -#' } -#' } -#' ctrl <- list(verbose = FALSE, # default control options -#' timer = FALSE, -#' zero = .Machine$double.eps, -#' max.iter = 500, -#' method = ifelse(type=='inf', "quadra","fista"), -#' threshold = ifelse(quadra, 1e-7, 1e-2)) -#' -#' ctrl[names(control)] <- control # overwritten by user specifications -#' if (ctrl$timer) {r.start <- proc.time()} -#' -#' ## ====================================================== -#' ## STARTING C++ CALL TO group_l1linf -#' if (ctrl$timer) {cpp.start <- proc.time()} -#' out <- .Call("group_lasso", -#' beta0 , -#' x , -#' y , -#' switch(type, "inf" = 1, 2), -#' pk , -#' struct , -#' lambda1 , -#' nlambda1 , -#' min.ratio , -#' penscale , -#' lambda2 , -#' intercept , -#' normalize , -#' rep(1,n) , -#' naive , -#' ctrl$thresh , -#' ctrl$max.iter, -#' max.feat , -#' switch(ctrl$method, -#' quadra = 0, -#' pathwise = 1, -#' fista = 2, 0), -#' ctrl$verbose, -#' inherits(x, "sparseMatrix"), -#' package = "quadrupen") -#' coefficients <- sparseMatrix(i = out$iA+1, -#' j = out$jA+1, -#' x = c(out$nzeros), -#' dims=c(length(out$lambda1),p)) -#' active.set <- sparseMatrix(i = out$iA+1, -#' j = out$jA+1, -#' dims=c(length(out$lambda1),p)) -#' ## END OF CALL -#' if (ctrl$timer) { -#' internal.timer <- (proc.time() - cpp.start)[3] -#' external.timer <- (proc.time() - r.start)[3] -#' } else { -#' internal.timer <- NULL -#' external.timer <- NULL -#' } -#' -#' ## ====================================================== -#' -#' ## BUILDING THE PENFIT OBJECT -#' out$converge[out$converge == 0] <- "converged" -#' out$converge[out$converge == 1] <- "max # of iterate reached" -#' out$converge[out$converge == 2] <- "max # of feature reached" -#' monitoring <- list(it.active = c(out$it.active ), -#' it.optim = c(out$it.optim ), -#' max.grad = c(out$max.grd ), -#' status = c(out$converge ), -#' pensteps.timer = c(out$timing ), -#' external.timer = external.timer , -#' internal.timer = internal.timer ) -#' -#' dimnames(coefficients)[[1]] <- round(c(out$lambda1),3) -#' if (is.null(colnames(x))) { -#' dimnames(coefficients)[[2]] <- 1:p -#' } else { -#' dimnames(coefficients)[[2]] <- colnames(x) -#' } -#' mu <- drop(out$mu) -#' df <- drop(out$df) -#' -#' ## FITTED VALUES AND RESIDUALS... -#' if (intercept) { -#' fitted <- sweep(tcrossprod(x,coefficients),2L,-mu,check.margin=FALSE) -#' df <- df + 1 -#' } else { -#' mu <- 0 -#' fitted <- tcrossprod(x,coefficients) -#' } -#' residuals <- apply(fitted,2,function(y.hat) y-y.hat) -#' r.squared <- 1-colSums(residuals^2)/ifelse(intercept,sum((y-mean(y))^2),sum(y^2)) -#' -#' return(new("quadrupen", -#' coefficients = coefficients , -#' active.set = active.set , -#' intercept = intercept , -#' mu = mu , -#' normx = drop(out$normx), -#' fitted = fitted , -#' residuals = residuals , -#' df = df , -#' r.squared = r.squared , -#' penscale = penscale , -#' penalty = switch(type, "inf" = "l1/linf group Lasso", "l1/l2 group Lasso"), -#' naive = naive , -#' lambda1 = c(out$lambda1) , -#' lambda2 = lambda2 , -#' monitoring = monitoring , -#' control = ctrl)) -#' -#' } diff --git a/inst/experimental/tmp_src/active_set.cpp b/inst/experimental/tmp_src/active_set.cpp deleted file mode 100644 index 5c6e177..0000000 --- a/inst/experimental/tmp_src/active_set.cpp +++ /dev/null @@ -1,167 +0,0 @@ -#include "active_set.h" - -using namespace Rcpp; -using namespace arma; - -// ______________________________________________________ -// ACTIVE SET CLASS -// -// handle vector of active variable, with method for adding/removing -// elements - -ACTIVE_SET::ACTIVE_SET(SEXP BETA0) { - vec beta0 = as(BETA0) ; - activeVar = find(beta0 != 0) ; - betaActive = beta0.elem(activeVar) ; - nbrIn = activeVar.n_elem ; - isActiveVar.zeros(beta0.n_elem) ; - for (int i=0; i=first_indVarOut; k--) { - del_var(k); - if (k == 0) {break;} - } -} - -// ______________________________________________________ -// ACTIVE DATA CLASS -// -// Inherits from the ACTIVE_SET class -// -// Contains "activated" part of the data required for optimization - -ACTIVE_DATA::ACTIVE_DATA(SEXP BETA0, REGRESSION_DATA data) : - ACTIVE_SET(BETA0) { - if (!betaActive.is_empty()) { - xtxA = mat(data.xt * data.x.cols(activeVar)) ; - smoothGrad = -data.xty + xtxA * betaActive ; - xAtxA = xtxA.rows(activeVar) ; - } -} - -void ACTIVE_DATA::add_var(uword varIn, REGRESSION_DATA data) { - - ACTIVE_SET::add_var(varIn) ; // update the active set - - vec new_col = data.xt * data.x.col(varIn); - - // if (lambda2 > 0) { - // Adding the column corresponding to the structurating matrix - // new_col += spS.col(var_in); - // } - - // UPDATE THE xtxA AND xAtxA MATRICES - if (!xAtxA.is_empty()) { - xAtxA = join_cols(xAtxA, xtxA.row(varIn)) ; - } - xtxA = join_rows(xtxA, new_col) ; - xAtxA = join_rows(xAtxA, trans(xtxA.row(varIn))) ; -} - -void ACTIVE_DATA::add_vars(uword first_varIn, uword last_varIn, REGRESSION_DATA data) { - - ACTIVE_SET::add_vars(first_varIn, last_varIn) ; // update the active set - - mat new_cols = data.xt * data.x.cols(first_varIn, last_varIn); - - // UPDATE THE xtxA AND xAtxA MATRICES - if (!xAtxA.is_empty()) { - xAtxA = join_cols(xAtxA, xtxA.rows(first_varIn, last_varIn)) ; - } - xtxA = join_rows(xtxA, new_cols) ; - xAtxA = join_rows(xAtxA, trans(xtxA.rows(first_varIn, last_varIn))) ; -} - - -void ACTIVE_DATA::del_var(uword indVarOut) { - ACTIVE_SET::del_var(indVarOut) ; - xtxA.shed_col(indVarOut) ; - xAtxA.shed_col(indVarOut) ; - xAtxA.shed_row(indVarOut) ; -} - -// indVarsOut MUST be sorted in the decreasing order -void ACTIVE_DATA::del_vars(uword first_indVarOut, uword last_indVarOut) { - ACTIVE_SET::del_vars(first_indVarOut, last_indVarOut) ; - xtxA.shed_cols(first_indVarOut, last_indVarOut) ; - xAtxA.shed_cols(first_indVarOut, last_indVarOut) ; - xAtxA.shed_rows(first_indVarOut, last_indVarOut) ; -} - -// handle vector of active variable, with method for adding/removing -// elements by group of variables - - -// GROUP_STRUCTURE:: -// GROUP_STRUCTURE(SEXP PK): -// pk (as(PK)) // a vector with the group numbers -// { -// K = pk.n_elem ; // number of groups -// gk_start.zeros(K) ; // starting index for each group -// gk_end.zeros (K) ; // starting index for each group -// group.zeros(sum(pk)) ; // indicator of group belonging -// is_active_grp.zeros(K) ; -// nbrGrpIn = 0 ; - -// gk_start[0] = 0 ; -// gk_end[0] = pk[0]; -// for (int k = 1; k < K; k++) { -// gk_start[k] = gk_start[k-1] + pk[k-1]; -// gk_end[k] = gk_start[k] + pk[k]; -// } -// for (int k = 0; k < K; k++) { -// group.subvec(gk_start[k], gk_end[k]-1) = k*arma::ones(pk[k]); -// } -// } - -// void GROUP_STRUCTURE::add_var(uword grpIn) { - -// active_grp.resize(nbrGrpIn+1) ; // update the active set -// active_grp[nbrGrpIn] = grpIn ; -// nbrGrpIn++ ; // update the number of active variable -// is_active_grp[grpIn] = 1 ; - -// for (int k=gk_start[grpIn]; k -#include "data_reg.h" - -using namespace Rcpp; -using namespace arma; - -class ACTIVE_SET { - -protected: - // VARIABLES FOR HANDLING THE ACTIVE SET - vec betaActive ; // vector of currently activated features - uvec activeVar ; // set of currently activated variables - uvec isActiveVar ; // a vector to check if a variable is already in the active set - uword nbrIn ; // number of active variables - -public: - ACTIVE_SET (SEXP BETA0); - - // ACTIVE SET HANDLING - void add_var(uword varIn) ; // add a variable in the active set - void del_var(uword indVarOut) ; // remove the variable activated in postition ind_var_out - void add_vars(uword first_varIn , uword last_varIn ) ; // the same for a set - void del_vars(uword first_indVarOut, uword last_indVarOut) ; // of contiguous variables - - const vec & get_betaActive() const { return betaActive ; } - const uvec & get_activeVar() const { return activeVar ; } - bool is_activeVar(uword i) { return (isActiveVar(i) == 1) ; } -}; - - -class ACTIVE_DATA: public ACTIVE_SET { -private: - bool l2Reg ; // TODO - mat xAtxA ; - mat xtxA ; - vec smoothGrad ; - -public: - ACTIVE_DATA(SEXP BETA0, REGRESSION_DATA data) ; - void add_var (uword varIn, REGRESSION_DATA data) ; // add a variable - void del_var(uword indVarOut) ; // remove a variable - void add_vars(uword first_varIn, uword last_varIn, REGRESSION_DATA data) ; // the same for a set - void del_vars(uword first_indVarOut, uword last_indVarOut) ; // of contiguous variables - - const mat& GramMat() const { return xAtxA ; } -} ; - -#endif - - -// class ACTIVE_SET_GROUPWISE: public ACTIVE_SET { - -// protected: -// // VARIABLES FOR HANDLING THE ACTIVE SET -// uvec active_grp ; // set of currently activated groups -// uvec is_active_grp ; // a vector to check if a variable is already in the active set -// uword nbrGrpIn ; // number of active variables -// uvec pk ; // vector of group numbers -// uword K ; // number of groups -// uvec gk_start ; // starting index for each group -// uvec gk_end ; // starting index for each group -// uvec group ; // indicator of group belonging for each variable - -// public: -// ACTIVE_SET_GROUPWISE (SEXP BETA0, SEXP PK); - -// // ACTIVE SET HANDLING -// void add_var(uword grpIn) ; // add a group of variables in the active set -// void del_var(uword indGrpOut) ; // remove the group of variable activated in postition ind_var_class - -// const uvec & get_active_grp() const { return active_grp ; } - -// }; - -// // ACTIVE_ALGORITHM -// // vec GROUP_LASSO::get_dist2opt(double lambda) { - -// // vec normGrad = max(zeros(K), grp_norm(smooth_grad) - lambda) ; -// // // normGrad.elem(active_group) = grp_norm(smooth_grad + lambda * betaA); - -// // return(normGrad) ; -// // } - diff --git a/inst/experimental/tmp_src/data_reg.cpp b/inst/experimental/tmp_src/data_reg.cpp deleted file mode 100644 index 772ca41..0000000 --- a/inst/experimental/tmp_src/data_reg.cpp +++ /dev/null @@ -1,62 +0,0 @@ -#include "data_reg.h" - -using namespace Rcpp; -using namespace arma; - -REGRESSION_DATA:: -REGRESSION_DATA(SEXP X , // matrix of features - SEXP Y , // vector of response - SEXP INTERCEPT , // should intercept be considered? - SEXP NORMALIZE , // should predictors be standardized? - SEXP PREWEIGHTS , // prediction weights - SEXP OBSWEIGHTS): // maximal number of features allowed to enter the path - - // BASICALLY READING THE INPUT ARGUMENTS FROM R - x (as (X)) , - y (as (Y)) , - intercept (as (INTERCEPT)) , - normalize (as (NORMALIZE)) , - preweights (as (PREWEIGHTS)), - obsweights (as (OBSWEIGHTS)) { - n = x.n_rows ; - p = x.n_cols ; -} - -void REGRESSION_DATA::standardize() { - - // TODO: OBSERVATION WEIGHTS - - if (intercept == 1) { - xbar = mean(x).t(); - ybar = mean(y) ; - for (int i=0; i - -using namespace Rcpp; -using namespace arma; - -// class OPTIM_DATA ; - -// TODO: use template to handle dense or sparse encoding (mat/sp_mat in armadillo) -class REGRESSION_DATA { - friend class ACTIVE_DATA; // NEED ACCES TO x, xt, y, xty, n, p - -private: - // DATA VARIABLES FOR REGRESSION PURPOSE - uword n ; // sample size - uword p ; // # of features - vec xbar ; // mean of the predictors - vec normx ; // norm of the predictors - double ybar ; // mean of the response - double normy ; // norm of the response - bool intercept ; // should intercept be considered? - bool normalize ; // should predictors be standardized? - vec obsweights ; // observation weights - vec preweights ; // predictor weights - mat x ; // matrix of predictors - mat xt ; // transpose matrix of predictor once and keep it to save time - vec xty ; // responses to predictors vector - vec y ; // vector of response - -public: - REGRESSION_DATA(SEXP X , // matrix of features - SEXP Y , // vector of response - SEXP INTERCEPT , // should intercept be considered? - SEXP NORMALIZE , // should predictors be standardized? - SEXP PREWEIGHTS , // prediction weights - SEXP OBSWEIGHTS); // observation weights - - // DATA NORMALIZATION - // intercept treatment, predictor standardization, predictor weihgts and observation weights - void standardize(); - - // various function to acces private members - const mat& X() const { return x ; } - const vec& Y() const { return y ; } - const vec& XTY() const { return xty ; } - -}; - -#endif diff --git a/inst/experimental/tmp_src/group_lasso_main.cpp b/inst/experimental/tmp_src/group_lasso_main.cpp deleted file mode 100644 index 81ed790..0000000 --- a/inst/experimental/tmp_src/group_lasso_main.cpp +++ /dev/null @@ -1,156 +0,0 @@ -/* - * Author: Julien CHIQUET - * Statistique et Génome - */ - -#include "group_lasso_main.h" - -using namespace Rcpp; -using namespace arma; - -SEXP group_lasso(SEXP BETA0 , - SEXP X , // matrix of features - SEXP Y , // vector of response - SEXP TYPE , // -1="l1/linf", 2="l1/l2" group-Lasso - SEXP PK , // successive groups length - SEXP OMEGA , // structuring matrix - SEXP LAMBDA1 , - SEXP NLAMBDA1 , - SEXP MINRATIO , - SEXP PENSCALE , - SEXP LAMBDA2 , - SEXP INTERCEPT, - SEXP NORMALIZE, - SEXP WEIGHTS , - SEXP NAIVE , - SEXP EPS , - SEXP MAXITER , - SEXP MAXFEAT , - SEXP FUN_OPTIM, - SEXP VERBOSE , - SEXP SPARSE ) { - - // ============================================================== - // INSTANTIATION OF THE REQUIRED OBJECTS - - // data instantiation, normalization and basics pretreatments - REGRESSION_DATA data(X, Y, INTERCEPT, NORMALIZE, PENSCALE, WEIGHTS) ; - data.standardize(); - - PENALTY lasso, bounded, grpLasso1(PK), grpLasso2(PK), coopLasso(PK) ; - grpLasso1.setPenalty("l1l2") ; - - // path instanciation and generation of the grid of lambda values - PATH path(MAXFEAT) ; - path.grid_penLevels(LAMBDA1, NLAMBDA1, MINRATIO, grpLasso1.dual_norm(data.XTY())); - - // ACTIVE SET + current Gram matrix - ACTIVE_DATA activeSet(BETA0, data); - - // ============================================================== - // COMPUTE THE PATH OF SOLUTIONS - // _____________________________________________________________ - // - // START THE LOOP OVER THE LEADING PENALTY - // timer.tic(); - uword npen = path.get_penLevels().n_elem; - bool verbose = as(VERBOSE) ; - double thres = as(EPS) ; - - // vec dualGap = zeros(npen) ; - // vec itActive = zeros(npen) ; - - // uword maxIter = as(MAXITER) ; - uword grpIn ; - - // vec dualGap.zeros(get_penlevels().n_elem) ; - for (int m=0; m eps & itActive[m] < maxIter) { - // _____________________________________________________________ - // - // (1) GROUP ACTIVATION IF APPLICABLE - // ____________________________________________________________ - - // if (is_active_group[grpIn] == 0) { - // add_group(grpIn) ; - // if (verbose == 2) {Rprintf("Group %i newly added.\n",grpIn);} - // } else { - // if (verbose == 2) {Rprintf("Group %i already in.\n",grpIn);} - // } - // // _____________________________________________________________ - // // - // // (2) OPTIMIZATION OVER THE CURRENTLY ACTIVATED VARIABLES - // // _____________________________________________________________ - // // - - // // _____________________________________________________________ - // // - // // (3) GROUP DELETION IF APPLICABLE - // // _____________________________________________________________ - // // - - - // // _____________________________________________________________ - // // - // // (4) OPTIMALITY TESTING - // // _____________________________________________________________ - - // itActive[m]++; - // R_CheckUserInterrupt(); // let the possibility to interrupt from R - // } // ending the active set algorihm - - } // ending the loop over the penalty - - // return List::create(Named("coefficients") = grp_lasso.get_coef(), - // Named("mu") = grp_lasso.get_mu() , - // Named("normx") = grp_lasso.get_normx() , - // Named("lambda1") = grp_lasso.get_penalty(), - // Named("df") = grp_lasso.get_df() ); - -} - - -// // _____________________________________________________________ -// // -// // NEW DECLARATION -// uword grpIn ; // currently added group -// uword varIn ; // currently added variable - -// // _____________________________________________________________ -// // -// // SOME INITIALIZATION -// itActive.zeros(penalty.n_elem) ; -// smooth_grad = -xty ; -// timing.zeros(penalty.n_elem) ; // successive timings -// df.zeros(penalty.n_elem) ; // degrees of freedom -// dualGap.zeros(penalty.n_elem) ; // distance to optimality -// is_active_group.zeros(K) ; - - // this works: - // lasso.setPenalty("l1") ; - // bounded.setPenalty("linf") ; - // grpLasso2.setPenalty("l1linf") ; - // coopLasso.setPenalty("coop") ; - - // lasso.proximal(data.get_xty(),1).print(); - // std::cout << "\n"; - // bounded.proximal(data.get_xty(),1).print(); - // std::cout << "\n"; - // grpLasso1.proximal(data.get_xty(),1).print(); - // std::cout << "\n"; - // grpLasso2.proximal(data.get_xty(),1).print(); - // std::cout << "\n"; - // coopLasso.proximal(data.get_xty(),1).print(); - // std::cout << "\n"; diff --git a/inst/experimental/tmp_src/group_lasso_main.h b/inst/experimental/tmp_src/group_lasso_main.h deleted file mode 100644 index 30dbd89..0000000 --- a/inst/experimental/tmp_src/group_lasso_main.h +++ /dev/null @@ -1,43 +0,0 @@ -/* - * Author: Julien CHIQUET - * Statistique et Génome - */ -#ifndef _quadrupen_GROUP_LASSO_MAIN_H -#define _quadrupen_GROUP_LASSO_MAIN__H - -#define ARMA_NO_DEBUG -#define ARMA_USE_LAPACK -#define ARMA_USE_BLAS - -#ifndef ARMA_HAVE_GETTIMEOFDAY -#define ARMA_HAVE_GETTIMEOFDAY -#endif - -#include "quadrupen_headers.h" - -using namespace Rcpp; -using namespace arma; - -RcppExport SEXP group_lasso(SEXP BETA0 , - SEXP X , // matrix of features - SEXP Y , // vector of response - SEXP TYPE , // -1="l1/linf", 2="l1/l2" group-Lasso - SEXP PK , // successive groups length - SEXP OMEGA , // structuring matrix - SEXP LAMBDA1 , - SEXP NLAMBDA1 , - SEXP MIN_RATIO, - SEXP PENSCALE , - SEXP LAMBDA2 , - SEXP INTERCEPT, - SEXP NORMALIZE, - SEXP WEIGHTS , - SEXP NAIVE , - SEXP EPS , - SEXP MAXITER , - SEXP MAXFEAT , - SEXP FUN_OPTIM, - SEXP VERBOSE , - SEXP SPARSE ) ; - -#endif diff --git a/inst/experimental/tmp_src/path.cpp b/inst/experimental/tmp_src/path.cpp deleted file mode 100644 index fbdef73..0000000 --- a/inst/experimental/tmp_src/path.cpp +++ /dev/null @@ -1,20 +0,0 @@ -#include "path.h" - -using namespace Rcpp; -using namespace arma; - -PATH::PATH (SEXP MAXFEAT) : - maxFeat (as(MAXFEAT)) -{} - -void PATH::grid_penLevels(SEXP PENLEVELS, - SEXP PENLEN , - SEXP MINRATIO , - double LMAX ) { - double lmax; - if (PENLEVELS != R_NilValue) { // when the iser provide one, just import it - penLevels = as(PENLEVELS) ; - } else { // otherwise compute a grid of penalties ona log scale. - penLevels = exp10(linspace(log10(LMAX), log10(as(MINRATIO)*LMAX), as(PENLEN))) ; - } -} diff --git a/inst/experimental/tmp_src/path.h b/inst/experimental/tmp_src/path.h deleted file mode 100644 index a2ef524..0000000 --- a/inst/experimental/tmp_src/path.h +++ /dev/null @@ -1,42 +0,0 @@ -/* - * Author: Julien CHIQUET - * Statistique et Génome - */ -#ifndef _quadrupen_PATH_H -#define _quadrupen_PATH_H - -#include - -using namespace Rcpp; -using namespace arma; - -class PATH { - -protected: - vec penLevels ; // the successive levels of penalties along the path - uword maxFeat ; // maximal number of features to enter the path - - // VARIABLES INDEXED BY PENALTY LEVELS - sp_mat beta ; // final matrix of fitted coefficients - vec mu ; // vector of the successive intercept term - vec df ; // vector of degrees of freedom along the path - mat nonzeros ; // contains non-zero value of beta - mat iA ; // contains row indices of the non-zero values - mat jA ; // contains column indices of the non-zero values - -public: - PATH(SEXP MAXFEAT) ; - - void grid_penLevels(SEXP PENLEVELS, - SEXP PENLEN , - SEXP MINRATIO , - double LMAX ) ; - - const vec & get_penLevels() const { return penLevels ; } - const uword & get_maxFeat() const { return maxFeat ; } - const vec & get_mu() const { return mu ; } - const vec & get_df() const { return df ; } - -}; - -#endif diff --git a/inst/experimental/tmp_src/penalties.cpp b/inst/experimental/tmp_src/penalties.cpp deleted file mode 100644 index 288f15e..0000000 --- a/inst/experimental/tmp_src/penalties.cpp +++ /dev/null @@ -1,272 +0,0 @@ -/* - * Author: Julien CHIQUET - * Statistique et Génome - */ - -#include "penalties.h" - -using namespace Rcpp; -using namespace arma; - -// Empty constructor for basic penalties -PENALTY::PENALTY() {} - -// Constructor requiring group information -PENALTY::PENALTY(SEXP PK) : pk (as(PK)){} - -void PENALTY::setPenalty(std::string penName){ - - // link public member functions to l1 norm functions - if (penName == "l1") { - _current_elt_norm_ptr = &PENALTY::elt_norm_L1; - _current_pen_norm_ptr = &PENALTY::pen_norm_L1; - _current_dual_norm_ptr = &PENALTY::dual_norm_L1; - _current_proximal_ptr = &PENALTY::proximal_L1; - - // link public member functions to linf norm functions - } else if (penName == "linf") { - _current_elt_norm_ptr = &PENALTY::elt_norm_LINF; - _current_pen_norm_ptr = &PENALTY::pen_norm_LINF; - _current_dual_norm_ptr = &PENALTY::dual_norm_LINF; - _current_proximal_ptr = &PENALTY::proximal_LINF; - - // link public member functions to l1/l2 norm functions - } else if (penName == "l1l2") { - _current_elt_norm_ptr = &PENALTY::elt_norm_L1L2; - _current_pen_norm_ptr = &PENALTY::pen_norm_L1L2; - _current_dual_norm_ptr = &PENALTY::dual_norm_L1L2; - _current_proximal_ptr = &PENALTY::proximal_L1L2; - - // link public member functions to l1/linf norm functions - } else if (penName == "l1linf") { - _current_elt_norm_ptr = &PENALTY::elt_norm_L1LINF; - _current_pen_norm_ptr = &PENALTY::pen_norm_L1LINF; - _current_dual_norm_ptr = &PENALTY::dual_norm_L1LINF; - _current_proximal_ptr = &PENALTY::proximal_L1LINF; - - // link public member functions to coop norm functions - } else if (penName == "coop") { - _current_elt_norm_ptr = &PENALTY::elt_norm_COOP; - _current_pen_norm_ptr = &PENALTY::pen_norm_COOP; - _current_dual_norm_ptr = &PENALTY::dual_norm_COOP; - _current_proximal_ptr = &PENALTY::proximal_COOP; - } -} - -vec PENALTY::elt_norm (vec x) {return((this->*_current_elt_norm_ptr)(x));} -double PENALTY::pen_norm (vec x) {return((this->*_current_pen_norm_ptr)(x));} -double PENALTY::dual_norm (vec x) {return((this->*_current_dual_norm_ptr)(x));} -vec PENALTY::proximal (vec x, double lambda) {return((this->*_current_proximal_ptr)(x,lambda));} - -// ______________________________________________________ -// L1 NORM A.K.A LASSO -vec PENALTY::elt_norm_L1(vec x) { - return(arma::abs(x)); -} - -double PENALTY::pen_norm_L1(vec x) { - return(sum(elt_norm(x))); -} - -double PENALTY::dual_norm_L1(vec x) { - return(max(elt_norm(x))) ; -} - -vec PENALTY::proximal_L1(vec x, double lambda) { - return(max(zeros(x.n_elem), 1-lambda/elt_norm(x)) % x); -} - -// ______________________________________________________ -// LINF NORM A.K.A BOUNDED REGRESSION -vec PENALTY::elt_norm_LINF(vec x) { - return(arma::abs(x)); -} - -double PENALTY::pen_norm_LINF(vec x) { - return(max(elt_norm(x))); -} - -double PENALTY::dual_norm_LINF(vec x) { - return(sum(elt_norm(x))) ; -} - -vec PENALTY::proximal_LINF(vec x, double lambda) { - int p = x.n_elem; - vec u, proj; - vec res = zeros(p); - - if ( as_scalar(sum(abs(x) / lambda)) >= 1) { - - // Reordering absolute values - u = sort(abs(x), "descend"); - - // values of the projected coordinate if non zero (dual problem) - proj = (cumsum(u) - lambda)/linspace(1,p,p); - - // selecting nonnull entries (dual) - uvec maxs = sort(find(u-proj>ZERO), "descend") ; - double thresh = proj[maxs[0]]; - - // solving primal problem - // We keep the smallest values and threshold the common values to +- thresh - for (int k=0; k < p;k++) { - if (fabs(x(k)) > ZERO) { - if (x(k) > 0) { - res(k) = fmin(fabs(x(k)),thresh); - } else { - res(k) = -fmin(fabs(x(k)),thresh); - } - } - } - } - return(res); -} - -// ______________________________________________________ -// L1/L2 NORM A.K.A GROUP-LASSO -vec PENALTY::elt_norm_L1L2(vec x) { - - vec res = zeros (pk.n_elem) ; // output with group norms - int ind = 0 ; // index to go through the groups - - for (int k=0; k(x.n_elem); - int k,ind = 0 ; - - vec tmp = max(zeros(pk.n_elem), 1-lambda/elt_norm(x)) ; - - for (k=0; k (pk.n_elem) ; // output with group norms - int ind = 0 ; // index to go through the groups - - for (int k=0; k (pk.n_elem) ; // output with group norms - int ind = 0 ; // index to go through the groups - - for (int k=0; k(sum(pk)); - - for (int k=0; k= 1) { - // Reordering absolute values - u = sort(abs(v), "descend"); - - // values of the projected coordinate if non zero (dual problem) - proj = (cumsum(u) - lambda)/linspace(1,p,p); - - // selecting nonnull entries (dual) - uvec maxs = sort(find(u-proj>ZERO), "descend") ; - double thresh = proj[maxs[0]]; - - // solving primal problem - // We keep the smallest values and threshold the common values to +- thresh - for (int j=0; j

ZERO) { - v(j) = fmin(fabs(v(j)),thresh); - } else { - v(j) = -fmin(fabs(v(j)),thresh); - } - } - } - - res.subvec(ind,ind+pk(k)-1) = v ; - ind += pk(k); - } - return(res); -} - -// ______________________________________________________ -// COOP(ERATIVE) NORM A.K.A COOPERATIVE-LASSO -vec PENALTY::elt_norm_COOP(vec x) { - - vec res = zeros (pk.n_elem) ; // output with group norms - int ind = 0 ; // index to go through the groups - - for (int k=0; k(x.n_elem); - int k,ind = 0 ; - - vec tmp = max(zeros(pk.n_elem), 1-lambda/elt_norm(x)) ; - - for (k=0; k - -#define ZERO 2e-16 // practical zero - -using namespace Rcpp; -using namespace arma; - -class PENALTY { -public: - PENALTY() ; - PENALTY(SEXP PK) ; - - void setPenalty(std::string penName) ; - - vec elt_norm (vec x) ; - double pen_norm (vec x) ; - double dual_norm (vec x) ; - vec proximal(vec x, double lambda) ; - -private: - uvec pk ; - - typedef vec (PENALTY::*elt_norm_ptr)(vec x) ; - elt_norm_ptr _current_elt_norm_ptr ; - - typedef double (PENALTY::*pen_norm_ptr)(vec x) ; - pen_norm_ptr _current_pen_norm_ptr ; - - typedef double (PENALTY::*dual_norm_ptr)(vec x) ; - dual_norm_ptr _current_dual_norm_ptr ; - - typedef vec (PENALTY::*proximal_ptr)(vec x, double lambda) ; - proximal_ptr _current_proximal_ptr ; - - vec elt_norm_L1 (vec x) ; - double pen_norm_L1 (vec x) ; - double dual_norm_L1 (vec x) ; - vec proximal_L1 (vec x, double lambda) ; - - vec elt_norm_LINF (vec x) ; - double pen_norm_LINF (vec x) ; - double dual_norm_LINF (vec x) ; - vec proximal_LINF (vec x, double lambda) ; - - vec elt_norm_L1L2 (vec x) ; - double pen_norm_L1L2 (vec x) ; - double dual_norm_L1L2 (vec x) ; - vec proximal_L1L2 (vec x, double lambda) ; - - vec elt_norm_L1LINF (vec x) ; - double pen_norm_L1LINF (vec x) ; - double dual_norm_L1LINF (vec x) ; - vec proximal_L1LINF (vec x, double lambda) ; - - vec elt_norm_COOP (vec x) ; - double pen_norm_COOP (vec x) ; - double dual_norm_COOP (vec x) ; - vec proximal_COOP (vec x, double lambda) ; -}; - -#endif diff --git a/inst/experimental/tmp_src/ridge.h b/inst/experimental/tmp_src/ridge.h deleted file mode 100644 index b3d187c..0000000 --- a/inst/experimental/tmp_src/ridge.h +++ /dev/null @@ -1,187 +0,0 @@ -/* - * Author: Julien CHIQUET - * UMR MIA Paris-Saclay - * julien.chiquet@inrae.fr - */ -#ifndef _quadrupen_RIDGE_H -#define _quadrupen_RIDGE_H - -#define ARMA_NO_DEBUG -#define ARMA_USE_LAPACK -#define ARMA_USE_BLAS - -#ifndef ARMA_HAVE_GETTIMEOFDAY -#define ARMA_HAVE_GETTIMEOFDAY -#endif - -#include "quadrupen_headers.h" - -using namespace Rcpp; -using namespace arma; - -class RIDGE { -public: - RIDGE(const arma::mat X, - const arma::vec Y, - const arma::mat C, - const arma::vec WEIGHTS, - const uword VERBOSE); - - void standardize(const bool INTERCEPT, const bool NORMALIZE); - - void get_lambda(SEXP LAMBDA, - const arma::uword NLAMBDA, - const double LAMBDAMIN, - const double LAMBDAMAX); - - void LetsRoll(); - - // various function to access private members - const arma::mat & get_coef() const { return beta; } - const arma::vec & get_lambda() const { return lambda; } - const arma::vec & get_df() const { return df; } - const arma::rowvec & get_normx() const { return normx; } - const arma::vec & get_mu() const { return mu; } - -private: - arma::mat x ; // matrix of predictors - arma::vec y ; // vector of reponses - arma::mat c ; // Cholesky decomposition of the structuring matrix - arma::mat w ; // vector of obervation weights - int verbose ; // integer for verbose mode - - arma::rowvec xty ; // responses to predictors vector - arma::rowvec xbar ; // mean of the predictors - arma::rowvec normx ; // norm of the predictors - double normy ; // norm of the response - double ybar ; // mean of the response - - arma::vec lambda ; // vector of ridge tuning parameters - arma::mat beta ; // matrix of coefficients - arma::vec df ; // degrees of freedom along the path - arma::vec mu ; // vector of intercept term - -}; - -#endif - -// RIDGE::RIDGE(const arma::mat X, -// const arma::vec Y, -// const arma::mat C, -// const arma::vec WEIGHTS, -// const arma::uword VERBOSE) : -// x (X), // predictors -// y (Y), // responses -// c (trimatu(C)), // Cholesky of L -// w (WEIGHTS), // observation weights -// verbose (VERBOSE) // verbose mode (unused) -// {} -// -// -// // [[Rcpp::export]] -// Rcpp::List ridge_old_cpp( -// const arma::mat& X , // matrix of features -// const arma::vec& Y , // vector of response -// const arma::mat& C , // Cholesky decompostion of the structuring matrix -// SEXP LAMBDA, -// const arma::uword& NLAMBDA , -// const double& LAMBDAMIN, -// const double& LAMBDAMAX, -// const bool& INTERCEPT, -// const bool& NORMALIZE, -// const arma::vec& WEIGHTS , -// const arma::uword& VERBOSE ) { -// -// // ============================================================== -// // INSTANTIATE THE RIDGE PATH -// RIDGE ridge(X, Y, C, WEIGHTS, VERBOSE); -// -// // ============================================================== -// // DATA NORMALIZATION -// ridge.standardize(INTERCEPT, NORMALIZE) ; -// -// // ============================================================== -// // GET LAMBDA VALUES -// ridge.get_lambda(LAMBDA, NLAMBDA, LAMBDAMIN, LAMBDAMAX) ; -// -// // ============================================================== -// // COMPUTE THE PATH OF SOLUTIONS -// ridge.LetsRoll(); -// -// return List::create(Named("coefficients") = ridge.get_coef() , -// Named("mu") = ridge.get_mu() , -// Named("normx") = ridge.get_normx() , -// Named("lambda2") = ridge.get_lambda(), -// Named("df") = ridge.get_df() ); -// -// } -// -// void RIDGE::LetsRoll() { -// -// arma::mat cinv = inv(trimatu(c)) ; // inverting the Cholesky decomp. of the structuring matrix -// -// // SVD DECOMPOSITION OF ( X * C^-1) -// arma::vec eta ; // eigen value of X cinv -// arma::mat U ; // left singular vectors of X -// arma::mat V ; // right singular vectors of X -// svd_econ(U, eta, V, x*cinv) ; -// -// arma::mat cinvV = cinv * V ; -// arma::mat Uty = trans(U) * y ; -// -// beta.resize(lambda.n_elem, x.n_cols); -// df.resize(lambda.n_elem, 1); -// -// for (uword i; i(p) ; -// ybar = 0; -// } -// -// if (normalize == 1) { -// normx = sqrt(sum(square(x),0)); -// x.each_row() /= normx ; -// } else { -// normx = ones(p); -// } -// normy = sqrt(sum(square(y))) ; -// -// if (intercept == 1) { -// xty = trans(y-ybar) * x ; -// } else { -// xty = y.t()*x ; -// } -// } -// -// void RIDGE::get_lambda(SEXP LAMBDA, -// const arma::uword NLAMBDA, -// const double LAMBDAMIN, -// const double LAMBDAMAX) { -// -// if (LAMBDA != R_NilValue) { -// lambda = as(LAMBDA) ; -// } else { -// lambda = exp10(linspace(log10(LAMBDAMIN), log10(LAMBDAMAX), NLAMBDA)) ; -// } -// -// } diff --git a/inst/oldies/elastic_net.cpp b/inst/oldies/elastic_net.cpp deleted file mode 100644 index 66bfad7..0000000 --- a/inst/oldies/elastic_net.cpp +++ /dev/null @@ -1,626 +0,0 @@ -/* - * Author: Julien CHIQUET, INRAE - * julien.chiquet@inrae.fr - * Statistique et Génome - * MIA Paris-Saclay - */ - -#include "quadrupen_headers.h" - -using namespace Rcpp; -using namespace arma; - -// [[Rcpp::export]] -Rcpp::List elastic_net_cpp( - const Environment &dataModel , // data structure - const bool &intercept , // boolean for intercept - const List ®Param , // config of the optimisation - const List &controlFit // config of the optimisation -) { - - const uword n = dataModel["n"] ; // sample size - const uword p = dataModel["d"] ; // problem size - const SEXP &R_X = dataModel["X"] ; // design matrix - const arma::vec &y = dataModel["y"] ; // response vector - const arma::sp_mat& S = dataModel["S"] ; // Structuring matrix - const arma::vec &wobs = dataModel["wy"] ; // observation weights (not use at the moment) - const bool sparse = dataModel["sparse_encoding"] ; // boolean for sparse mode - - const SEXP R_LAMBDA1 = regParam["l1"] ; // vector of L1 penalties ; if NULL, automatically set - const arma::vec &wlambda1 = regParam["l1_weights"] ; // l1-penalty weights - const double lambda_l2 = regParam["l2"] ; // scalar for the amount L2 penalty - uword n_lambda = regParam["n_lambda1"] ; // # of l1-penalty levels - const double min_ratio = regParam["min_ratio"] ; // minimum penlaty value as a ratio of lambda1 max - - const bool normalize = controlFit["normalize"] ; // boolean for standardizing the predictor - const vec beta0 = controlFit["beta0"] ; // initial vector of coefficients - const double eps = controlFit["threshold"] ; // precision required - const double eps2 = pow(eps, 2) ; - const arma::uword maxiter = controlFit["maxiter"] ; // max # of iterates of the active set - const arma::uword maxfeat = controlFit["maxfeat"] ; // max # of variables activated - const arma::uword fun = controlFit["method"] ; // solver (0=quadra, 1=pathwise, 2=fista) - const arma::uword verbose = controlFit["verbose"] ; // int for verbose mode (0/1/2) - const bool usechol = controlFit["usechol"] ; // use Cholesky decomposition or not - const arma::uword monitor = controlFit["monitor"] ; // convergence monitoring (1 == Grandvalet's bound ;-) 2 == Fenchel duality gap) - - vec xty ; // responses to predictors vector - vec xbar ; // mean of the predictors - vec normx ; // norm of the predictors - double normy ; // norm of the response - double ybar ; // mean of the response - - // STRUCTURATING MATRIX (embed lambda_l2) - sp_mat diag_S = spdiags(sqrt(lambda_l2)*pow(wlambda1,-1/2), ivec({0}), p, p) ; - const sp_mat S_lambda_l2 = diag_S * S * diag_S ; // sparsely encoded structuring matrix - - // DESIGN MATRIX (either sparsely encoded or not) - mat x ; - mat xt ; - sp_mat sp_x ; - sp_mat sp_xt ; - if (sparse) { // Check how x is encoded for reading - sp_x = as(R_X) ; - standardize(sp_x, y, intercept, normalize, wlambda1, xty, normx, normy, xbar, ybar) ; - sp_xt = sp_x.t() ; - } else { - x = as(R_X) ; - standardize(x, y, intercept, normalize, wlambda1, xty, normx, normy, xbar, ybar) ; - xt = x.t(); - } - - // Sparsifying penalty - vec lambda_l1 = get_lambda1(R_LAMBDA1, n_lambda, min_ratio, max(abs(xty))) ; - n_lambda = lambda_l1.n_elem ; // # of penalty levels - - // Initializing "first level" variables (outside of the lambda_l1 loop) - mat R ; // Cholesky decomposition of XAtXA - uvec A ; // set of currently activated variables - vec betaA ; // vector of currently activated parameters - mat xtxA ; // t(x) * x_A covariance matrix - mat xAtxA ; // t(x_A) * x_A + S covariance matrix of the activated variable plus SAA matrix - vec xtxw ; // t(x_A) * x_A * beta(A) - vec grd = -xty ; // smooth part of the gradient - vec mu = zeros(n_lambda) ; // the intercept term - vec max_grd = zeros(n_lambda) ; // a vector with the successively reach duality gap - vec converge = zeros(n_lambda) ; // a vector indicating if convergence occurred (0/1/2) - uvec it_active = zeros(n_lambda) ; // # of loop in the active set for each lambda_l1 - uvec it_optim ; // # of loop in the optimization process for each loop of the active se - double L0 = 1.0 + lambda_l2 ; // Lipschitz constant for proximal methods - vec timing (n_lambda) ; // successive timing in - vec df (n_lambda) ; // degrees of freedom - wall_clock timer ; // clock - - // Initializing "second level" variables (within the active set - for a fixed value of lambda_l1) - uword var_in ; // currently added variable - uword nbr_in = 0 ; // # of currently added variables - uword nbr_opt = 0 ; // # of current calls to the optimization routine - uvec are_in = zeros(p) ; // a vector to check if a variable is already in the active set - List out_optim ; // the list of output of the optimization function - bool success_optim = true ; // was the internal system resolution successful? - uvec null ; // stores the variables which go to zero during optimization - vec grd_norm (p) ; // current value of the grd_norm for each variable - vec nonzeros ; // contains non-zero value of beta - urowvec iA ; // contains row indices of the non-zero values - urowvec jA ; // contains column indices of the non-zero values - - // WARM START - A = find(beta0 != 0) ; - betaA = beta0.elem(A) ; - if (A.n_elem > 0) { - if (sparse) { - xtxA = mat(sp_xt * sp_x.col(0)) - n * xbar * (xbar.elem(A)).t(); - } else { - xtxA = mat(xt * x.cols(A)) - n * xbar * (xbar.elem(A)).t(); - } - if (lambda_l2 > 0) { - for (uword i=0; i eps) && (it_active[m] < maxiter)) { - // _____________________________________________________________ - // - // (1) VARIABLE ACTIVATION IF APPLICABLE - // _____________________________________________________________ - - if (verbose == 2) {Rprintf("max grad %f\n",max_grd[m]);} - - // Check if the variable is already in the active set - if (are_in[var_in] == 0) { - if (sparse) { - add_var_enet(n, nbr_in, var_in, betaA, A, sp_x, sp_xt, xtxA, xAtxA, xtxw, R, lambda_l2, xbar, S_lambda_l2, usechol, fun) ; - } else { - add_var_enet(n, nbr_in, var_in, betaA, A, x, xt, xtxA, xAtxA, xtxw, R, lambda_l2, xbar, S_lambda_l2, usechol, fun) ; - } - if (verbose == 2) {Rprintf("newly added variable %i\n",var_in);} - are_in[var_in] = 1; - nbr_in++; - } else { - if (verbose == 2) {Rprintf("already in %i\n",var_in);} - } - - // _____________________________________________________________ - // - // (2) OPTIMIZATION OVER THE CURRENTLY ACTIVATED VARIABLES - // _____________________________________________________________ - // - - it_optim.reshape(nbr_opt + 1,1) ; - switch (fun) { - case 1 : - it_optim[nbr_opt] = pathwise_enet(betaA, xAtxA, xty.elem(A), xtxw, lambda_l1[m], null, lambda_l2, eps2); - break; - case 2 : - it_optim[nbr_opt] = fista_lasso(betaA, xAtxA, xty.elem(A), lambda_l1[m], null, L0, eps2); - break; - default: - try { - it_optim[nbr_opt] = quadra_enet(betaA, R, xAtxA, xty.elem(A), sign(grd.elem(A)), lambda_l1[m], null, usechol, eps); - } catch (std::runtime_error &error) { - if (verbose > 0) { - Rprintf("\nWarning: singular system at this stage of the solution path, cutting here.\n"); - } - success_optim = false ; - } - } - // update the smooth part of the gradient - grd = -xty + xtxA * betaA; - nbr_opt++; - - // _____________________________________________________________ - // - // (3) VARIABLE DELETION IF APPLICABLE - // _____________________________________________________________ - // - // removing variables zeroed during optimization - if (!null.is_empty()) { - if (verbose == 2) { - for (uword j=0; j 0) { - // _____________________________________________________________ - // - // (OPTIONAL) FOLLOWING CONVERGENCE BY COMPLETE MONITORING - // _____________________________________________________________ - bound_to_optimal(betaA, xAtxA, xty, grd, lambda_l1[m], lambda_l2, normy, A, monitor, J_hat, D_hat) ; - } - - // Moving to the next iterate - it_active[m]++; - - R_CheckUserInterrupt(); - } - - // degrees of freedom - df[m] = get_df_enet(lambda_l2, R, xAtxA, S_lambda_l2, A, usechol, fun); - - // the reference parameter (obtained once optimum is met) - if (monitor > 0) { - if (it_active[m] > 0) { - J_star = join_cols(J_star, ones(it_active[m],1) * J_hat[nbr_opt-1]) ; - } - } - - // Record the time ellapsed - timing[m] = timer.toc() ; - - // Checking convergence status - if (it_active[m] >= maxiter) converge[m] = 1 ; - if (nbr_in > maxfeat) converge[m] = 2 ; - if (!success_optim) converge[m] = 3 ; - - // Stop now if relevant - if (converge[m] == 2 || converge[m] == 3) { - lambda_l1 = lambda_l1.subvec(0,m-1) ; - mu = mu.subvec(0,m-1) ; - converge = converge.subvec(0,m) ; - max_grd = max_grd.subvec(0,m-1) ; - it_active = it_active.subvec(0,m) ; - timing = timing.subvec(0,m) ; - df = df.subvec(0,m-1) ; - break; - } else { - nonzeros = join_cols(nonzeros, betaA/(normx.elem(A) % wlambda1.elem(A))); - mu[m] = dot(betaA, xbar.elem(A)) ; - iA = join_rows(iA, m*ones(betaA.n_elem) ); - jA = join_rows(jA, A.t()) ; - } - } - - // Updating monitored quantities - if (monitor > 0) D_star = J_hat - J_star; - - return List::create( - Named("tuning_param") = List::create( - Named("l1") = lambda_l1, - Named("l2") = lambda_l2 - ), - Named("beta") = sp_mat(join_cols(iA, jA), nonzeros, lambda_l1.n_elem, p), - Named("active") = sp_mat(join_cols(iA, jA), vec(iA.n_elem, fill::ones), lambda_l1.n_elem, p), - Named("mu") = ybar - mu, - Named("df") = df , - Named("monitoring") = - List::create( - Named("it_active") = it_active, - Named("it_optim") = it_optim , - Named("max_grd") = max_grd , - Named("convergence") = converge , - Named("pensteps_timer") = timing , - Named("delta_hat") = D_hat , - Named("delta_star") = D_star - ) - ); - -} - -// // ============================================================================ -// // OLD VERSION -// // -// -// // [[Rcpp::export]] -// Rcpp::List elastic_net_old_cpp( -// SEXP BETA0 , // -// SEXP X , // regressor matrix -// const arma::vec y , // response vector -// const arma::sp_mat Struct , // Structuring matrix -// SEXP LAMBDA1 , // -// double n_lambda , // -// const double min_ratio , // -// const arma::vec penscale , // penalty weights -// const double lambda2 , // the smooth (ridge) penalty -// const bool intercept , // boolean for intercept mode -// const bool normalize , // boolean for standardizing the predictor -// const arma::vec weights , // observation weights (not use at the moment) -// const bool naive , // naive elastic-net or not -// const double eps , // precision required -// const arma::uword maxiter , // max # of iterates of the active set -// const arma::uword maxfeat , // max # of variables activated -// const arma::uword fun , // solver (0=quadra, 1=pathwise, 2=fista) -// const arma::uword verbose , // int for verbose mode (0/1/2) -// const bool sparse , // boolean for sparse mode -// const bool usechol , // use Cholesky decomposition or not -// const arma::uword monitor // convergence monitoring (1 == Grandvalet's bound ;-) 2 == Fenchel duality gap) -// ) { -// -// vec xty ; // responses to predictors vector -// vec xbar ; // mean of the predictors -// vec normx ; // norm of the predictors -// double normy ; // norm of the response -// double ybar ; // mean of the response -// const double eps2 = pow(eps, 2) ; -// const uword n(y.n_elem) ; // sample size -// const uword p(Struct.n_cols) ; // problem size -// -// mat x ; -// mat xt ; -// sp_mat sp_x ; -// sp_mat sp_xt ; -// if (sparse == 1) { // Check how x is encoded for reading -// sp_x = as(X) ; -// standardize(sp_x, y, intercept, normalize, penscale, xty, normx, normy, xbar, ybar) ; -// sp_xt = sp_x.t() ; -// } else { -// x = as(X) ; -// standardize(x, y, intercept, normalize, penscale, xty, normx, normy, xbar, ybar) ; -// xt = x.t(); -// } -// -// // STRUCTURATING MATRIX -// sp_mat diag_S = spdiags(sqrt(lambda2)*pow(penscale,-1/2), ivec({0}), p, p) ; -// const sp_mat S = diag_S * Struct * diag_S ; // sparsely encoded structuring matrix -// mat SAA ; // densely encoded active counterpart -// -// // VECTOR OF TUNING PARAMETER FOR THE L1-PENALTY -// vec lambda1 = get_lambda1(LAMBDA1, n_lambda, min_ratio, max(abs(xty))); -// n_lambda = lambda1.n_elem ; // # of penalty levels -// -// // Initializing "first level" variables (outside of the lambda1 loop) -// mat R ; // Cholesky decomposition of XAtXA -// uvec A ; // set of currently activated variables -// vec betaA ; // vector of currently activated parameters -// mat xtxA ; // t(x) * x_A covariance matrix -// mat xAtxA ; // t(x_A) * x_A + S covariance matrix of the activated variable plus SAA matrix -// vec xtxw ; // t(x_A) * x_A * beta(A) -// vec grd = -xty ; // smooth part of the gradient -// vec mu = zeros(n_lambda) ; // the intercept term -// vec max_grd = zeros(n_lambda) ; // a vector with the successively reach duality gap -// vec converge = zeros(n_lambda) ; // a vector indicating if convergence occured (0/1/2) -// uvec it_active = zeros(n_lambda) ; // # of loop in the active set for each lambda1 -// uvec it_optim ; // # of loop in the optimization process for each loop of the active se -// double L0 = 1.0 + lambda2 ; // Lipschitz constant for proximal methods -// vec timing (n_lambda) ; // succesive timing in -// vec df (n_lambda) ; // degrees of freedom -// wall_clock timer ; // clock -// -// // Initializing "second level" variables (within the active set - for a fixed value of lamdba) -// uword var_in ; // currently added variable -// uword nbr_in = 0 ; // # of currently added variables -// uword nbr_opt = 0 ; // # of current calls to the optimization routine -// uvec are_in = zeros(p) ; // a vector to check if a variable is already in the active set -// List out_optim ; // the list of output of the optimization function -// bool success_optim = true ; // was the internal system resolution successful? -// uvec null ; // stores the variables which go to zero during optimization -// vec grd_norm (p) ; // current value of the grd_norm for each variable -// mat nonzeros ; // contains non-zero value of beta -// mat iA ; // contains row indices of the non-zero values -// mat jA ; // contains column indices of the non-zero values -// -// vec beta0 ; -// // WARM START -// if (BETA0 != R_NilValue) { -// beta0 = as(BETA0) ; -// A = find(beta0 != 0) ; -// betaA = beta0.elem(A) ; -// if (sparse == 1) { -// // WRONG - DO IT THE RIGHT WAY -// xtxA = mat(sp_xt * sp_x.col(0)) ; -// } else { -// xtxA = mat(xt * x.cols(A)) ; -// } -// if (lambda2 > 0) { -// for (uword i=0; i eps) && (it_active[m] < maxiter)) { -// // _____________________________________________________________ -// // -// // (1) VARIABLE ACTIVATION IF APPLICABLE -// // _____________________________________________________________ -// -// // Check if the variable is already in the active set -// if (are_in[var_in] == 0) { -// if (sparse == 1) { -// add_var_enet(n, nbr_in, var_in, betaA, A, sp_x, sp_xt, xtxA, xAtxA, xtxw, R, lambda2, xbar, S, usechol, fun) ; -// } else { -// add_var_enet(n, nbr_in, var_in, betaA, A, x, xt, xtxA, xAtxA, xtxw, R, lambda2, xbar, S, usechol, fun) ; -// } -// if (verbose == 2) {Rprintf("newly added variable %i\n",var_in);} -// are_in[var_in] = 1; -// nbr_in++; -// } else { -// if (verbose == 2) {Rprintf("already in %i\n",var_in);} -// } -// -// // _____________________________________________________________ -// // -// // (2) OPTIMIZATION OVER THE CURRENTLY ACTIVATED VARIABLES -// // _____________________________________________________________ -// // -// -// it_optim.reshape(nbr_opt + 1,1) ; -// switch (fun) { -// case 1 : -// it_optim[nbr_opt] = pathwise_enet(betaA, xAtxA, xty.elem(A), xtxw, lambda1[m], null, lambda2, eps2); -// break; -// case 2 : -// it_optim[nbr_opt] = fista_lasso(betaA, xAtxA, xty.elem(A), lambda1[m], null, L0, eps2); -// break; -// default: -// try { -// it_optim[nbr_opt] = quadra_enet(betaA, R, xAtxA, xty.elem(A), sign(grd.elem(A)), lambda1[m], null, usechol, eps); -// } catch (std::runtime_error &error) { -// if (verbose > 0) { -// Rprintf("\nWarning: singular system at this stage of the solution path, cutting here.\n"); -// } -// success_optim = false ; -// } -// } -// // update the smooth part of the gradient -// grd = -xty + xtxA * betaA; -// nbr_opt++; -// -// // _____________________________________________________________ -// // -// // (3) VARIABLE DELETION IF APPLICABLE -// // _____________________________________________________________ -// // -// // removing variables zeroed during optimization -// if (!null.is_empty()) { -// if (verbose == 2) { -// for (uword j=0; j 0) { -// // _____________________________________________________________ -// // -// // (OPTIONAL) FOLLOWING CONVERGENCE BY COMPLETE MONITORING -// // _____________________________________________________________ -// bound_to_optimal(betaA, xAtxA, xty, grd, lambda1[m], lambda2, normy, A, monitor, J_hat, D_hat) ; -// } -// -// // Moving to the next iterate -// it_active[m]++; -// -// R_CheckUserInterrupt(); -// } -// -// // degrees of freedom -// df[m] = get_df_enet(lambda2, R, xAtxA, S, A, fun); -// -// // the reference parameter (obtained once optimum is met) -// if (monitor > 0) { -// if (it_active[m] > 0) { -// J_star = join_cols(J_star, ones(it_active[m],1) * J_hat[nbr_opt-1]) ; -// } -// } -// -// // Record the time ellapsed -// timing[m] = timer.toc() ; -// -// // Checking convergence status -// if (it_active[m] >= maxiter) { -// converge[m] = 1; -// } -// if (nbr_in > maxfeat) { -// converge[m] = 2 ; -// } -// if (!success_optim) { -// converge[m] = 3; -// } -// -// // Stop now if relevant -// if (converge[m] == 2 || converge[m] == 3) { -// lambda1 = lambda1.subvec(0,m-1) ; -// converge = converge.subvec(0,m) ; -// max_grd = max_grd.subvec(0,m-1) ; -// it_active = it_active.subvec(0,m) ; -// timing = timing.subvec(0,m) ; -// df = df.subvec(0,m) ; -// break; -// } else { -// nonzeros = join_cols(nonzeros, betaA/(normx.elem(A) % penscale.elem(A))); -// iA = join_cols(iA, m*ones(betaA.n_elem,1) ); -// jA = join_cols(jA, conv_to::from(A) ) ; -// if (intercept == 1) { -// mu[m] = dot(betaA, xbar.elem(A)) ; -// } -// } -// } -// // END OF THE LOOP OVER LAMBDA -// if (!naive) { -// nonzeros *= 1+lambda2; -// mu = ybar - (1+lambda2) * mu; -// } else { -// mu = ybar - mu; -// } -// -// // Updating monitored quantities -// if (monitor > 0) { -// D_star = J_hat - J_star; -// } -// -// return List::create(Named("nzeros") = nonzeros , -// Named("iA") = iA , -// Named("jA") = jA , -// Named("mu") = mu , -// Named("normx") = normx , -// Named("lambda1") = lambda1 , -// Named("df") = df , -// Named("nbr.in") = nbr_in , -// Named("it.active") = it_active, -// Named("it.optim") = it_optim , -// Named("max.grd") = max_grd , -// Named("timing") = timing , -// Named("delta.hat") = D_hat , -// Named("delta.star") = D_star , -// Named("converge") = converge ); -// -// } -// diff --git a/inst/oldies/first_order.cpp b/inst/oldies/first_order.cpp deleted file mode 100644 index 04d941d..0000000 --- a/inst/oldies/first_order.cpp +++ /dev/null @@ -1,317 +0,0 @@ -/* - * Author: Julien CHIQUET - * Statistique et Génome - */ -#include "first_order.h" - -int fista_lasso(vec &x0 , - mat &xtx , - vec xty , - double &pen, - uvec &null, - double &L0 , - double eps) { - - colvec xk = x0 ; // output vector - colvec s = x0 ; - int iter = 0 ; // current iterate - double delta = 2*eps ; // change in beta - int max_iter = 10000 ; // max. number of iteration - double L = L0 ; // Lipschitz's constant - - double t0 = 1.0, tk; - bool found=false; - double f0, fk ; - colvec df ; - - double l_num, l_den ; - - while ((delta > eps/x0.n_elem ) && (iter < max_iter)) { - - f0 = as_scalar(.5 * strans(s) * xtx * s - strans(xty) * s) ; - df = - xty + xtx * s ; - - // Line search over L - while(!found) { - // apply proximal operator of the Lasso - xk = s - df/L ; - for (uword j=0; j= l_num) || (sqrt(l_den) < eps)) { - found = true; - } else { - L = fmax(2*L, l_num/l_den); - } - - R_CheckUserInterrupt(); - } - - // updating t - tk = 0.5 * (1+sqrt(1+4*t0*t0)); - - // updating s - s = xk + (t0-1)/tk * ( xk - x0 ); - - // preparing next iterate - delta = sqrt(l_num); - x0 = xk; - t0 = tk; - found = false; - iter++; - - R_CheckUserInterrupt(); - } - - null = sort(find(abs(xk) + (abs(df) - pen) < ZERO), "descend") ; - - return(iter); -} - -int fista_grp(vec &x0, - uvec &pk, - mat &xtx, - vec xty, - vec &pen, - uword &grpnorm, - double &L0 , - double eps ) { - - colvec xk = x0 ; // output vector - colvec s = x0 ; - int iter = 0 ; // current iterate - double delta = 2*eps ; // change in beta - int max_iter = 10000 ; // max. number of iteration - double L = L0 ; // Lipschitz's constant - - double t0 = 1.0, tk; - bool found=false; - double f0, fk ; - colvec df ; - - double l_num, l_den ; - - while ( (delta > eps/x0.n_elem) && (iter < max_iter)) { - - f0 = as_scalar(.5 * strans(s) * xtx * s - strans(xty) * s) ; - df = (xtx * s - xty) ; - - // Line search over L - while(!found) { - // apply proximal operator of the group-Lasso - if (grpnorm == 2) { - xk = proximal_grp2(s - df/L, pk, pen, L); - } else { - xk = proximal_grpinf(s - df/L, pk, pen, L); - } - - fk = as_scalar(.5 * strans(xk) * xtx * xk - strans(xty) * xk) ; - l_num = as_scalar(2 * (fk - f0 - dot(df, xk-s) )); - l_den = as_scalar(pow(norm(xk-s,2),2)); - - if ( (L * l_den >= l_num) | (sqrt(l_den) < eps)) { - found = true; - } else { - L = fmax(2*L, l_num/l_den); - } - - R_CheckUserInterrupt(); - } - - // updating t - tk = 0.5 * (1+sqrt(1+4*t0*t0)); - - // updating s - s = xk + (t0-1)/tk * ( xk - x0 ); - - // preparing next iterate - delta = sqrt(l_num); - x0 = xk; - t0 = tk; - found = false; - iter++; - - R_CheckUserInterrupt(); - } - - return(iter); -} - -vec proximal_grp2(vec u, - uvec pk, - vec lambda, - double L) { - - uword ind = 0 ; - - vec grp_norm2 = grp_norm(u, pk, 2, 1); - for (uword k=0; k(p); - - if ( as_scalar(sum(abs(v) / lambda)) >= 1) { - - // Reordonnons les valeurs absolues - u = sort(abs(v), "descend"); - - // valeurs des coordonnées projetees si non nulles (problème dual) - proj = (cumsum(u) - lambda)/linspace(1,p,p); - - // selection des coordonnees non nulles (problème dual) - uvec maxs = sort(find(u-proj>ZERO), "descend") ; - double thresh = proj[maxs[0]]; - - // solution du programme primal - // On garde les valeurs les plus petites, et on seuille le reste à une valeur commune +- thresh - for (uword k=0; k < p;k++) { - if (fabs(v(k)) > ZERO) { - if (v(k) > 0) { - out(k) = fmin(fabs(v(k)),thresh); - } else { - out(k) = -fmin(fabs(v(k)),thresh); - } - } - } - } - return(out); -} - -int pathwise_enet(vec& x0, - mat& xtx, - vec xty, - vec& xtxw, - double& pen, - uvec &null, - const double& gam , - const double eps ) { - - colvec xk = x0 ; // output vector - int iter = 0 ; // current iterate - int max_iter = 10000 ; // max. number of iteration - double delta = 2*eps ; // change in beta - double u, d ; // temporary scalar - - while ((delta > eps/x0.n_elem ) && (iter < max_iter)) { - - delta = 0; - for (uword j=0; j eps/x0.n_elem ) && (iter < max_iter)) { - -// delta = 0; -// ind = 0; -// // blockwise coordinate descent -// for (k=0; k pen(ind)) { - -// // valeurs des coordonnees projetees si non nulles -// vec proj = (cumsum(us) - pen(ind)) /linspace(1,x0.n_elem,x0.n_elem); - -// // selection des coordonnees non nulles -// uword jmax; -// proj.max(jmax); - -// for (j=ind; j<(ind+pk.at(k)); j++) { - -// if (abs(u(iter)) > eps/100) { -// if (u(iter) > 0) { -// xk(j) = fmin(abs(u(i)),proj(jmax)); -// } else { -// xk(j) = -fmin(abs(u(i)),proj(jmax)); -// } -// } -// } -// } else { -// xk.subvec(ind,ind+pk.at(k)-1) = zeros(pk.at(k)); -// } -// d = xk.subvec(ind,ind+pk.at(k)-1)-x0.subvec(ind,ind+pk.at(k)-1); -// delta += norm(d,2); -// xtxw += xtx.cols(ind,ind+pk.at(k)-1) * d; - -// ind += pk.at(k); - -// } - -// // preparing next iterate -// delta = sqrt(delta); -// x0 = xk; -// iter++; - -// R_CheckUserInterrupt(); -// } - -// return(iter); -// } diff --git a/inst/oldies/first_order.h b/inst/oldies/first_order.h deleted file mode 100644 index f9f8007..0000000 --- a/inst/oldies/first_order.h +++ /dev/null @@ -1,28 +0,0 @@ -/* - * Author: Julien CHIQUET - * Statistique et Génome - */ -#ifndef _quadrupen_FIRSTORDER_H -#define _quadrupen_FIRSTORDER_H - -#include -#include "utils.h" - -using namespace Rcpp; -using namespace arma; - -#define ZERO 2e-16 // practical zero - -int fista_lasso(vec &x0, mat &xtx, vec xty, double &pen, uvec &null, double &L0, const double eps) ; -int fista_grp(vec &x0, uvec &pk, mat &xtx, vec xty, vec &pen, uword &grpnorm, double &L0, const double eps) ; - -vec proximal_grp2(vec u, uvec pk, vec lambda, double L) ; -vec proximal_grpinf(vec u, uvec pk, vec lambda, double L) ; -vec proximal_inf(vec v, vec lambda) ; - -int pathwise_enet(vec& x0, mat& xtx, vec xty, vec& xtxw, double& pen, uvec &null, const double& gam, const double eps) ; - -// int pathwise_grpl1linf(vec &x0, uvec &pk, mat &xtx, vec xty, vec &xtxw, vec &pen, double &gam, double eps) ; - - -#endif diff --git a/inst/oldies/quadratic.cpp b/inst/oldies/quadratic.cpp deleted file mode 100644 index 5746a0d..0000000 --- a/inst/oldies/quadratic.cpp +++ /dev/null @@ -1,75 +0,0 @@ -/* - * Author: Julien CHIQUET - */ -#include "quadratic.h" - -int quadra_enet(vec &x0, - mat &R, - mat &xAtxA, - vec xty, - vec sgn_grd, - double &pen , - uvec &null, - bool usechol, - double tol) { - - uword iter = 1; // current iterate - - uvec A = find(abs(x0) > ZERO) ; // vector of active variables - vec theta = -sgn_grd ; // vector of sign of the solution - theta.elem(A) = sign(x0.elem(A)); - - // Solving the quadratic problem - vec x1 ; - if (usechol) { - x1 = solve(trimatu(R), solve(trimatl(strans(R)), xty - pen * theta)); - } else { - x1 = cg(xAtxA, xty - pen*theta, x0, tol) ; - } - - // Check for swapping variables - uvec swap = find(abs(sign(x1.elem(A)) - theta.elem(A)) > ZERO); - if (swap.is_empty()) { - null = swap ; // this is empty - x0 = x1; - } else { - swap = A.elem(swap); - colvec x0_swap = x0.elem(swap); - colvec x1_swap = x1.elem(swap); - // first, go to zero for the swapped variable which cost the minimum - vec gamma = -x0_swap / (x1_swap-x0_swap); - uword i_swap = gamma.index_min(); - double scale = gamma(i_swap) ; - null = swap[i_swap]; - x1 = x0 + (x1-x0) * scale ; - // second, solve the problem after swaping the signs of the - // incriminated variable - x0 = x1; - x0(null[0]) = -x1_swap[i_swap]; - - A = find(abs(x0) > ZERO) ; // vector of active variables - theta = -sgn_grd ; // vector of sign of the solution - theta.elem(A) = sign(x0.elem(A)); - - vec x2 ; - if (usechol) { - x2 = solve(trimatu(R), solve( trimatl(strans(R)), xty - pen * theta)); - } else { - x2 = cg(xAtxA, xty - pen*theta, x1, tol) ; - } - iter++; - - // This is the gradient on the active part of the parameters - vec grd = -xty + xAtxA * x2; - // if the sign is coherent, keep that one... - if (fabs(grd(null[0]) + pen * as_scalar(sign(x2(null)))) <= ZERO) { - null = swap; // this is empty - x0 = x2 ; - } else { - // otherwise, backtrack to x1 - x0 = x1 ; - } - } - - return(iter); -} diff --git a/inst/oldies/quadratic.h b/inst/oldies/quadratic.h deleted file mode 100644 index 4cc633d..0000000 --- a/inst/oldies/quadratic.h +++ /dev/null @@ -1,19 +0,0 @@ -/* - * Author: Julien CHIQUET - * Statistique et Génome - */ - -#ifndef _quadrupen_QUADRATIC_H -#define _quadrupen_QUADRATIC_H - -#define ARMA_WARN_LEVEL 1 - -#include -#include "utils.h" - -using namespace Rcpp; -using namespace arma; - -int quadra_enet(vec& x0, mat& R, mat& xAtxA, vec xty, vec sgn_grd, double &pen, uvec& null, bool usechol, double tol) ; - -#endif diff --git a/inst/oldies/quadrupen_headers.h b/inst/oldies/quadrupen_headers.h deleted file mode 100644 index d00f716..0000000 --- a/inst/oldies/quadrupen_headers.h +++ /dev/null @@ -1,25 +0,0 @@ -/*** - * @file interface.h - * - */ -#ifndef __QUADRUPEN_HEADERS_HPP -#define __QUADRUPEN_HEADERS_HPP - -#define RCPP_ARMADILLO_RETURN_ANYVEC_AS_VECTOR - -// [[Rcpp::depends(RcppArmadillo)]] - -// Include Armadillo / Rcpp / R to C/C++ basics -#include "RcppArmadillo.h" -#include -#include - -#define ZERO 2e-16 // practical zero - -// Include utils and optimization routines -#include "utils.h" -#include "quadratic.h" -#include "first_order.h" - -#endif - diff --git a/inst/oldies/utils.cpp b/inst/oldies/utils.cpp deleted file mode 100644 index 4684d4f..0000000 --- a/inst/oldies/utils.cpp +++ /dev/null @@ -1,300 +0,0 @@ -#include "utils.h" - -using namespace Rcpp; -using namespace arma; - -// arma::uvec setdiff(arma::uvec& x, arma::uvec& y) { -// -// std::vector a = arma::conv_to< std::vector >::from(arma::sort(x)); -// std::vector b = arma::conv_to< std::vector >::from(arma::sort(y)); -// std::vector out; -// -// std::set_difference(a.begin(), a.end(), b.begin(), b.end(), -// std::inserter(out, out.end())); -// -// return arma::conv_to::from(out); -// } - -vec grp_norm(vec x, uvec pk, int norm, int rep) { - - vec res ; - if (rep > 0) { - res = zeros (x.n_elem) ; - } else { - res = zeros (pk.n_elem) ; - } - - int ind = 0; - double current_norm ; - - switch (norm) { - // -1 is the convention for the sup norm - case -1 : - for (uword k=0; k 0) { - res.subvec(ind, ind + pk(k) - 1) = current_norm*ones(pk(k)); - } else { - res(k) = current_norm; - } - ind += pk(k); - } - break; - case 1 : - for (uword k=0; k 0) { - res.subvec(ind, ind + pk(k) - 1) = current_norm*ones(pk(k)); - } else { - res(k) = current_norm; - } - ind += pk(k); - } - break; - case 2: - for (uword k=0; k 0) { - res.subvec(ind, ind + pk(k) - 1) = current_norm*ones(pk(k)); - } else { - res(k) = current_norm; - } - ind += pk(k); - } - break; - default: - for (uword k=0; k 0) { - res.subvec(ind, ind + pk(k) - 1) = current_norm*ones(pk(k)); - } else { - res(k) = current_norm; - } - ind += pk(k); - } - break; - } - - return(res); -} - -vec grp_sign(vec x, uvec pk) { - - vec signs = zeros (x.n_elem); - int ind = 0; - double eps = 1e-12; - - vec norm_grp = grp_norm(x, pk, -1, 0) ; - - for (uword k=0; k eps) { - signs(j) = 1; - } else { - signs(j) = -1; - } - } - } - ind += pk(k); - } - return(signs); -} - -void cholupdate(mat &R , mat &XtX) { - int p = XtX.n_cols; - - if (p == 1) { - R = sqrt(XtX); - } else { - colvec rp = zeros(p,1); - rp.subvec(0,p-2) = solve (trimatl(strans(R)), XtX.submat(0,p-1,p-2,p-1)); - rp(p-1) = sqrt(XtX(p-1,p-1) - dot(rp,rp)); - R = join_rows( join_cols(R, zeros(1,p-1)) , rp); - } -} - -void choldowndate(mat &R, int j) { - - vec x = zeros(2,1); - mat G = zeros(2,2); - mat H = zeros(2,2); - - R.shed_col(j); - int p = R.n_cols; - double r; - for (int k=j; k 0) { - if ((fun == 0) & (usechol)) { - B = solve(trimatu(R), eye(R.n_cols, R.n_cols)); - B = B * B.t(); - } else { - B = inv_sympd(xAtxA); - } - // have to do this due to sparse encoding - // either wait for Armadillo's guy to develop non contiguous - // subview for sparse matrix or iterate over the n_zeros only... - for (uword i=0;i nu); // set of adversarial variables outside the boundary - // Grandvalet's bound - D_hat(dim) = J_hat(dim) - (lambda1/nu) * J_hat(dim) - (pow(lambda1,2)/(2*lambda2))*((lambda1*(p-Ac.n_elem))/nu + pow(norm(gamma.elem(Ac),2)/nu,2)-p); - } else { - // Fenchel's bound - if (nu < lambda1) { - nu = lambda1; - } - D_hat(dim) = 0.5 * quad_loss * (1+pow(lambda1/nu,2)) + sum(abs(lambda1*betaA)) + (lambda1/nu)*(dot(betaA,xty.elem(A))-pow(normy,2)); - } - - // keep the smallest bound reached so far for a given lambda value - if (dim>0) { - if (J_hat[dim-1] < J_hat[dim]) { - if (D_hat[dim] > D_hat[dim-1] - (J_hat[dim-1] - J_hat[dim])) { - D_hat[dim] = D_hat[dim-1]; - } - } - } - -} - -vec cg(mat A, vec b, vec x, double tol) { - - vec r = b - A * x; - vec p = r ; - double rs_old = sum(square(r)) ; - - double rs_new = rs_old ; - int i = 0; - double alpha ; - mat Ap ; - - while ((sqrt(rs_new) > tol) & (i < 1e3)) { - Ap = A * p; - alpha = rs_old/dot(p,Ap) ; - x += alpha * p ; - r -= alpha * Ap ; - // Polak-Ribière update - rs_new = dot(r,-alpha * Ap); - p = r + rs_new/rs_old*p; - rs_old = rs_new; - i++; - } - - // Rprintf("\n nb of iterate %d",i); - return(x); -} - -// Can't find a reasonable Preconditioner that does not -// require a computational burden equivalent to a Cholesky decomposition -vec pcg(mat A, mat P, vec b, vec x, double tol) { - - vec r = b - A * x; - vec z = P * r; - vec p = z ; - //double rs_old = sum(square(r)) ; - double rs_old = dot(r,z) ; - - double rs_new = rs_old ; - int i = 0; - double alpha ; - mat Ap ; - - while ((sqrt(rs_new) > tol) & (i < 1e3)) { - Ap = A * p; - alpha = rs_old/dot(p,Ap) ; - x += alpha * p ; - r -= alpha * Ap ; - // Polak-Ribière update - z = P * r ; - rs_new = dot(z,-alpha * Ap); - p = z + rs_new/rs_old*p; - rs_old = rs_new; - i++; - } - - // Rprintf("\n nb of iterate %d",i); - return(x); -} diff --git a/inst/oldies/utils.h b/inst/oldies/utils.h deleted file mode 100644 index 3d25a62..0000000 --- a/inst/oldies/utils.h +++ /dev/null @@ -1,118 +0,0 @@ -/* - * Author: Julien CHIQUET - * Statistique et Génome - */ -#ifndef _quadrupen_UTILS_H -#define _quadrupen_UTILS_H - -#include - -using namespace Rcpp; -using namespace arma; - -#define ZERO 2e-16 // practical zero - -inline vec get_lambda1(SEXP LAMBDA1, uword n_lambda, const double min_ratio, double lmax) { - vec lambda1 ; - if (LAMBDA1 != R_NilValue) { - lambda1 = as(LAMBDA1) ; - } else { - lambda1 = exp10(linspace(log10(lmax), log10(min_ratio*lmax), n_lambda)) ; - } - return(lambda1); -} - -vec grp_norm(vec x, uvec pk, int norm, int rep) ; -vec grp_sign(vec x, uvec pk) ; - -double get_df_enet(const double &lambda2, mat &R, mat &xAtxA, const sp_mat &S, uvec &A, const bool &usechol, const uword &fun) ; - -vec cg(mat A, vec b, vec x, double tol) ; -vec pcg(mat A, mat P, vec b, vec x, double tol) ; - -void cholupdate(mat &R, mat& XtX) ; - -void choldowndate(mat &R, int j) ; - -void bound_to_optimal(vec &betaA, mat &xAtxA, const vec &xty, vec &grd, double &lambda1, const double &lambda2, const double &normy, uvec &A, const uword &monitor, vec &J_hat, vec &D_hat) ; - -template -void add_var_enet(const uword &n, uword &nbr_in, uword &var_in, vec &betaA, uvec &A, any_mat &x, any_mat &xt, mat &xtxA, mat &xAtxA, mat &xtxw, mat &R, const double &lambda2, const vec &xbar, const sp_mat &spS, const bool &usechol, const uword &fun) { - - vec new_col ; // column currently added to xtxA - - A.resize(nbr_in+1) ; // update the active set - A[nbr_in] = var_in ; - betaA.resize(nbr_in+1) ; // update the vector of active parameters - betaA[nbr_in] = 0.0 ; - - new_col = xt * x.col(var_in) - n * xbar * as_scalar(xbar[var_in]); - if (lambda2 > 0) { - // Adding the column corresponding to the structurating matrix - new_col += spS.col(var_in); - } - - // UPDATE THE xtxA AND xAtxA MATRICES - if (nbr_in > 0) { - xAtxA = join_cols(xAtxA, xtxA.row(var_in)) ; - } - xtxA = join_rows(xtxA, new_col) ; - xAtxA = join_rows(xAtxA, trans(xtxA.row(var_in))) ; - - if ((fun == 0) & (usechol == 1)) { - cholupdate(R, xAtxA) ; - } - - if (fun == 1) { - xtxw.resize(nbr_in+1) ; - xtxw(nbr_in) = dot(xAtxA.col(nbr_in),betaA); - } -} - -void remove_var_enet(uword &nbr_in, uvec &are_in, vec &betaA, uvec &A, mat &xtxAS, mat &xAtxA, mat &xtxw, mat &R, uvec &null, const bool &usechol, const uword &fun) ; - -template -void standardize(any_mat &x, const vec &y, const bool &intercept, const bool &normalize, const vec &wlambda1, - vec &xty, vec &normx, double &normy, vec &xbar, double &ybar) { - - uword n = x.n_rows; - uword p = x.n_cols; - - if (intercept) { - xbar = trans(rowvec(mean(x, 0))); - ybar = mean(y) ; - } else { - xbar = zeros(p) ; - ybar = 0; - } - - if (normalize) { - normx = sqrt(trans(sum(square(x),0)) - n * square(xbar)); - for (uword i=0; i}} - \preformatted{CrossValidation$plotCV_1D(log_scale = TRUE, title = "Cross-validation error")} + \preformatted{CrossValidation$plotCV_1D( + log_scale = TRUE, + title = "Cross-validation error", + se = TRUE +)} \if{html}{\out{}} } \subsection{Arguments}{ @@ -112,6 +116,7 @@ Should be called internally by an object \code{\link[=QuadrupenFit]{QuadrupenFit \describe{ \item{\code{log_scale}}{logical, should a log-scale be used for the x-axis} \item{\code{title}}{graph title} + \item{\code{se}}{logical, should confidence band be displayed (TRUE by default)} } \if{html}{\out{}} } diff --git a/man/StabilityPath.Rd b/man/StabilityPath.Rd index f700d23..b1565b2 100644 --- a/man/StabilityPath.Rd +++ b/man/StabilityPath.Rd @@ -174,7 +174,7 @@ selection. sel_mode = c("rank", "PFER"), cutoff = 0.75, PFER = 2, - nvarsel = floor(self$nobs/log(self$nvar)) + nvarsel = min(self$nvar, floor(self$nobs/log(self$nvar))) )} \if{html}{\out{}} } @@ -201,8 +201,6 @@ Buhlmannn's paper. Default is \code{"rank"}.} relevant when \code{sel_mode} equals \code{"PFER"}).} \item{\code{nvarsel}}{number of variables selected (only relevant when \code{sel_mode} equals \code{"rank"}. Default is \code{floor(n/log(p))}.} - \item{\code{log_scale}}{logical; indicates if a log-scale should be used -when \code{xvar="lambda"}. Default is \code{TRUE}.} \item{\code{plot}}{logical; indicates if the graph should be plotted. Default is \code{TRUE}. If \code{FALSE}, only the \pkg{ggplot2} object is sent back.} diff --git a/memcheck/tests.R b/memcheck/tests.R deleted file mode 100644 index 4f81ee3..0000000 --- a/memcheck/tests.R +++ /dev/null @@ -1,4 +0,0 @@ -# Code to be run with -# R -d "valgrind --tool=memcheck --leak-check=full --error-exitcode=1" --vanilla < memcheck/tests.R -library("quadrupen") -devtools::test() diff --git a/src/Quadrupen/BaseLava.h b/src/Quadrupen/BaseLava.h index 36b9ab3..5aa6275 100644 --- a/src/Quadrupen/BaseLava.h +++ b/src/Quadrupen/BaseLava.h @@ -95,7 +95,7 @@ class BaseLava : public ConcreteBase { this->set_.inverse_Gram() ; mat K = diagmat(ones(this->data_.n_)) - this->data_.X_.cols(this->set_.A_) * this->set_.XATXAinv_ * this->data_.X_.cols(this->set_.A_).t() ; - df -= trace(K * Proj_) ; + df += trace(K * Proj_) ; return df ; } diff --git a/vignettes/.gitignore b/vignettes/.gitignore new file mode 100644 index 0000000..ba7abf3 --- /dev/null +++ b/vignettes/.gitignore @@ -0,0 +1,3 @@ +*.html +*.R +*_files diff --git a/vignettes/group-sparse-regression.Rmd b/vignettes/group-sparse-regression.Rmd new file mode 100644 index 0000000..c946dea --- /dev/null +++ b/vignettes/group-sparse-regression.Rmd @@ -0,0 +1,233 @@ +--- +title: "Sparse Group Regression with Quadrupen" +subtitle: "A tour of standard analysis (fit, cross-validation, information criteria, stability)" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{group-sparse-regression} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.width = 7, + fig.height = 4 +) +``` + +## Setup + +```{r setup} +library(quadrupen) +data("Birthwt", package = "grpreg") +y <- Birthwt$bwt[-130] ## outlier +X <- Birthwt$X[-130, ] +group <- as.integer(Birthwt$group)[-130] +``` + +The `Birthwt` dataset contains `r nrow(Birthwt$X)` observations and `r ncol(Birthwt$X)` predictors organized into `r nlevels(Birthwt$group)` clinically meaningful groups; the response is birth weight (in kg). Observation 130 is excluded as an outlier throughout this vignette. The `group` argument expected by `group_sparse_lm()` is a sorted integer vector with one entry per column of `X`. + +```{r data-prep} +group_names <- levels(Birthwt$group) +var_labels <- group_names[group] +cat("Groups (", length(group_names), "):", paste(group_names, collapse = ", "), "\n") +cat("Group sizes:", tabulate(group), "\n") +``` + +## Fitting group-sparse linear models + +The main entry point is `group_sparse_lm()`, which fits a regularization path for several group-sparse penalties controlled by the `type` and `alpha` arguments. Convenience wrappers are available for the most common variants: + +| Wrapper | `type` | `alpha` | Penalty | +|---|---|---|---| +| `group_lasso()` | `"l2"` | 0 | Group Lasso ($\ell_1/\ell_2$): group-level sparsity | +| `coop_lasso()` | `"coop"` | 0 | Cooperative Lasso: group sparsity + within-group sign coherence | +| `sparse_group_lasso()` | `"l2"` | > 0 | Sparse Group Lasso: group + individual sparsity | + +### Group Lasso + +```{r gl-fit} +fit_gl <- group_lasso(X, y, group) +fit_gl +``` + +### Cooperative Lasso + +The Cooperative Lasso promotes coherent group selection by penalizing non-zero within-group coefficients that differ in sign. When a group enters the model, all its active coefficients are forced to share the same sign. + +```{r cl-fit} +fit_cl <- coop_lasso(X, y, group) +fit_cl +``` + +### Sparse Group Lasso + +`alpha` controls the mixture between the group-level penalty ($\alpha = 0$, pure Group Lasso) and an element-wise $\ell_1$ penalty ($\alpha = 1$, pure Lasso). Intermediate values such as `alpha = 0.5` enforce group-level sparsity while also allowing individual predictors within an active group to be zeroed out. + +```{r sgl-fit} +fit_sgl <- sparse_group_lasso(X, y, group, alpha = 0.5) +fit_sgl +``` + +### Fit objects + +All objects returned by `group_sparse_lm()` and its wrappers are R6 instances of the class `SparseGroupFit`, inheriting from `QuadrupenFit`. These objects (see [`QuadrupenFit`]) store all data related to the fit, accessible via named fields (e.g., `fit_gl$coefficients`, `fit_gl$deviance`, `fit_gl$degrees_freedom`; see `str(fit_gl)` and the documentation). + +They also provide methods for visualizing and analyzing the fit, and for pursuing complementary analyses such as model selection, cross-validation, and stability selection. For users unfamiliar with R6 classes, S3 methods are exported for the most common operations (e.g., `plot(fit_gl)` is equivalent to `fit_gl$plot()`), but the R6 methods expose more options. + +## Regularization paths + +`$plot_path()` (or `plot(fit, type = "path")`) displays how coefficients evolve along the penalty path. The `labels` argument adds a legend coloured by group name. + +```{r path-plots, fig.show='hold', out.width='50%'} +fit_gl$plot_path(xvar = "fraction", log_scale = FALSE, labels = var_labels) +fit_cl$plot_path(labels = var_labels) +fit_sgl$plot_path(labels = var_labels) +fit_sgl$plot_path(standardize = FALSE) +``` + +## Information criteria + +`$criteria()` computes AIC, BIC, mBIC, eBIC and GCV from the estimated degrees of freedom, storing the result as an [`InformationCriteria`] R6 object in the field `$information_criteria` of the current fit. This method is called automatically during fitting at no extra cost, so users rarely need to invoke it directly — the main exception being when a custom penalty term is desired. + +```{r criteria} +fit_gl$criteria() +``` + +```{r criteria-plots, fig.show='hold', out.width='50%'} +fit_gl$plot(type = "criteria") +``` + +For greater flexibility, use the `plot` method of the `InformationCriteria` object directly: + +```{r criteria-plots-2, fig.show='hold', out.width='50%'} +fit_gl$information_criteria$plot(c("AIC", "BIC", "mBIC")) +fit_sgl$information_criteria$plot("GCV") +``` + +## Model extraction + +`$get_model()` returns the coefficient vector selected by a given criterion. With group penalties it is particularly informative to examine which groups are entirely selected, which are partially active (Sparse Group Lasso only), and which are excluded. + +```{r get-model} +coef_gl <- fit_gl$get_model("BIC") +coef_sgl <- fit_sgl$get_model("BIC") + +active_gl <- unique(group[coef_gl[-1] != 0]) +active_sgl <- unique(group[coef_sgl[-1] != 0]) +cat("Group Lasso — active groups (BIC):", group_names[active_gl], "\n") +cat("Sparse GL — active groups (BIC):", group_names[active_sgl], "\n") +``` + +The Sparse Group Lasso can additionally zero out individual predictors within an active group: + +```{r sgl-within} +active_vars_sgl <- which(coef_sgl[-1] != 0) +cat("Sparse GL — active predictors (BIC):", colnames(X)[active_vars_sgl], "\n") +``` + +An existing lambda value can also be used directly: + +```{r get-model-lambda} +lambda_gl <- fit_gl$major_tuning +coef_gl_lambda <- fit_gl$get_model(lambda_gl[20]) +cat("Non-zero Group Lasso coefficients for lambda =", round(lambda_gl[20], 3), "\n") +print(coef_gl_lambda[coef_gl_lambda != 0]) +``` + +## Cross-validation + +`$cross_validate()` performs K-fold cross-validation over the penalty grid, storing the result as a [`CrossValidation`] R6 object in the field `$cross_validation` of the current fit. + +```{r cv, message=FALSE} +set.seed(42) +fit_gl$cross_validate(K = 10, verbose = FALSE) +fit_sgl$cross_validate(K = 10, verbose = FALSE) +``` + +```{r cv-plots, fig.show='hold', out.width='50%'} +fit_gl$plot(type = "crossval") +fit_sgl$plot(type = "crossval") +``` + +Model selection using the CV-minimizing penalty: + +```{r cv-model} +coef_gl_cv <- fit_gl$get_model("CV_min") +cat("Group Lasso — active groups (CV_min):", + group_names[unique(group[coef_gl_cv[-1] != 0])], "\n") +``` + +## Cross-validation on $\lambda_1 \times \lambda_2$ + +For regularizers with two penalties, cross-validation can be run on a two-dimensional grid. Here we explore a range of `lambda2` values jointly with the Group Lasso path: + +```{r cv2D, message=FALSE} +set.seed(42) +fit_gl$cross_validate(K = 5, lambda2 = 10^seq(1, -3, len = 20), verbose = FALSE) +``` + +```{r cv-plots-2D, fig.show='hold', out.width='50%'} +fit_gl$plot(type = "crossval") +fit_gl$cross_validation$plotCV_1D(se = FALSE) + ggplot2::ylim(c(0.04, 0.08)) +``` + +## Prediction + +`$predict()` returns fitted values for the whole path, or for a specific model selected by a criterion or a lambda value. + +```{r predict} +y_hat_gl <- fit_gl$predict(selection = "CV_min") +y_hat_cl <- fit_cl$predict(selection = "BIC") +y_hat_sgl <- fit_sgl$predict(selection = "CV_min") + +r2_gl <- fit_gl$r_squared[fit_gl$get_model("CV_min", type = "index")] +r2_cl <- fit_cl$r_squared[fit_cl$get_model("BIC", type = "index")] +r2_sgl <- fit_sgl$r_squared[fit_sgl$get_model("CV_min", type = "index")] +cat(sprintf( + "R² Group Lasso (CV_min): %.3f\nR² Coop Lasso (BIC) : %.3f\nR² Sparse GL (CV_min): %.3f\n", + r2_gl, r2_cl, r2_sgl +)) +``` + +## Stability selection + +`$stability()` estimates selection probabilities via repeated sub-sampling (Meinshausen & Bühlmann, 2010). Variables that appear consistently across sub-samples receive high probabilities and are considered robustly selected. The stability path is stored as a [`StabilityPath`] R6 object in the field `$stability_path` of the current fit. + +```{r stability, message=FALSE} +set.seed(42) +fit_gl$stability(n_subsamples = 200, verbose = FALSE) +``` + +```{r stability-plot} +fit_gl$plot(type = "stability", labels = var_labels) +fit_gl$stability_path$plot(nvarsel = 5, labels = var_labels) +colnames(X)[fit_gl$stability_path$selection(nvarsel = 5)] +``` + +## Debiased Estimator + +Group penalties induce shrinkage bias on the magnitude of estimated coefficients. A standard remedy is to refit the model on the selected support without penalization. This is implemented in **quadrupen** at no extra computational cost — the refitted coefficients are computed alongside the regularization path — and can be activated by setting the `$debias` field to `TRUE`. + +```{r debias} +fit_gl$debias <- TRUE +fit_gl$plot_path(labels = var_labels) +``` + +Since both versions are stored in the same object, it is straightforward to compare the regularized and debiased solutions — for the path as well as for CV error — by toggling `$debias`. + +```{r debias-cv} +fit_gl$cross_validate(verbose = FALSE) +fit_gl$plot(type = "crossval") +fit_gl$debias <- FALSE +fit_gl$plot_path(labels = var_labels) +``` + +Setting `$debias` back to `FALSE` restores the original regularized coefficients. + +```{r no-debias} +fit_gl$debias <- FALSE +fit_gl$plot_path(labels = var_labels) +``` diff --git a/vignettes/lava.Rmd b/vignettes/lava.Rmd new file mode 100644 index 0000000..5d703aa --- /dev/null +++ b/vignettes/lava.Rmd @@ -0,0 +1,298 @@ +--- +title: "LAVA: Recovering Sums of Sparse and Dense Signals" +subtitle: "The lava() and group_lava() estimators in quadrupen" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{lava} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.width = 7, + fig.height = 4 +) +``` + +## Motivation + +Classical sparse estimators such as the Lasso assume that the true coefficient vector $\theta$ is itself sparse: only a few entries are nonzero. This assumption fails when the signal has a **mixed structure** — a small number of large, identifiable effects superimposed on a broad background of many small effects. In this case neither the Lasso (too sparse) nor ridge regression (too dense) is well-suited. + +The **LAVA** estimator (Chernozhukov, Hansen & Liao, *Annals of Statistics*, 2017) explicitly decomposes the coefficient vector as + +$$\theta = \beta + \delta,$$ + +where $\beta$ is the **sparse** component (penalized by $\lambda_1 \|\beta\|_1$) and $\delta$ is the **dense** component (penalized by $\frac{\lambda_2}{2} \delta^\top S\, \delta$). The joint criterion is + +$$\min_{\theta = \beta + \delta} \frac{1}{2n} \|y - X(\beta + \delta)\|^2 + \lambda_1 \|\beta\|_1 + \frac{\lambda_2}{2}\, \delta^\top S\, \delta,$$ + +where $S$ is a positive semidefinite structuring matrix (identity by default). By varying $\lambda_1$ along a regularization path (with $\lambda_2$ fixed), **quadrupen** fits this model via `lava()` and returns an object of class `LavaFit` that gives access to both components separately. + +## Setup + +```{r setup, message = FALSE} +library(quadrupen) +library(ggplot2) +``` + +### Simulation: a mixed sparse + dense signal + +We simulate a linear model where the true signal is the sum of a sparse part $\beta$ (five isolated large effects) and a dense part $\delta$ (small Gaussian perturbations on every predictor). Neither component alone characterizes the signal. + +```{r simulation} +set.seed(42) +n <- 100 +p <- 50 + +## Sparse component: 5 large effects, the rest zero +beta <- numeric(p) +beta[c(5, 15, 25, 35, 45)] <- c(3, -3, 3, -3, 3) + +## Dense component: small effects on all predictors +delta <- rnorm(p, mean = 0, sd = 0.2) + +## Total signal +theta <- beta + delta + +## Predictors with moderate AR(1) correlation +rho <- 0.5 +Sigma <- toeplitz(rho^(0:(p - 1))) +X <- matrix(rnorm(n * p), n, p) %*% chol(Sigma) +y <- X %*% theta + rnorm(n, sd = 2) +``` + +```{r plot-truth, fig.height = 5} +df_truth <- data.frame( + index = rep(1:p, 3), + value = c(beta, delta, theta), + component = factor( + rep(c("Sparse (β)", "Dense (δ)", "Total (θ = β + δ)"), each = p), + levels = c("Sparse (β)", "Dense (δ)", "Total (θ = β + δ)") + ) +) + +ggplot(df_truth, aes(x = index, y = value, fill = component)) + + geom_col() + + facet_wrap(~ component, ncol = 1, scales = "free_y") + + scale_fill_manual(values = c("#56B4E9", "#E69F00", "#009E73"), guide = "none") + + labs(x = "Predictor index", y = "Value", + title = "True signal decomposition") + + theme_bw() +``` + +## Fitting with `lava()` + +`lava()` takes the same interface as `sparse_lm()`. The `lambda2` argument controls the ridge-like penalty on the dense component. + +```{r lava-fit} +fit_lava <- lava(X, y, lambda2 = 1) +fit_lava +``` + +The returned object is an R6 instance of class `LavaFit`, inheriting from `QuadrupenFit`. In addition to all the standard fields and methods (`$coefficients`, `$criteria()`, `$cross_validate()`, etc.), it exposes two component-specific fields: + +- `$sparse_coef`: the sparse part $\hat\beta$ — a $p \times |\lambda_1|$ matrix along the path. +- `$dense_coef`: the dense part $\hat\delta = \hat\theta - \hat\beta$ — same dimensions. + +## Regularization paths by component + +`$plot_path()` accepts a `component` argument (`"both"`, `"sparse"`, `"dense"`) to display each part of the decomposition independently. + +```{r paths, fig.show='hold', out.width='33%'} +fit_lava$plot_path(component = "both", labels = NULL) +fit_lava$plot_path(component = "sparse", labels = NULL) +fit_lava$plot_path(component = "dense", labels = NULL) +``` + +As $\lambda_1$ increases (right to left on the path): + +- The **sparse** path shows the usual Lasso-like variable selection: large effects enter early and small ones remain zero. +- The **dense** path is smooth and non-zero everywhere: the ridge penalty ensures that $\hat\delta$ is distributed across all predictors. +- The **total** path combines both: it never reaches exactly zero because the dense component picks up residual signal. + +## Model Selection + +Let us plot the information criteria along the path. The BIC and eBIC select a model with a couple of large effects in the sparse part, while AIC selects a more complex model with more non-zero entries in the sparse component. mBIC is more conservative, as expected. + +```{r lava-criteria} +fit_lava$information_criteria$plot(c("AIC", "BIC", "eBIC", "mBIC")) +``` + +## Model extraction and decomposition + +`$get_model()` extracts the total estimated coefficient $\hat\theta$ at a selected $\lambda_1$. The corresponding sparse and dense components can be recovered from `$sparse_coef` and `$dense_coef` at the same index. + +```{r extraction} +idx <- fit_lava$get_model("BIC", type = "index") + +theta_hat <- fit_lava$get_model("BIC")[-1] # drop intercept +beta_hat <- as.numeric(fit_lava$sparse_coef[, idx]) +delta_hat <- as.numeric(fit_lava$dense_coef[, idx]) + +cat("Non-zero entries in sparse component (BIC):", + sum(beta_hat != 0), "\n") +cat("Non-zero entries in dense component (BIC): ", + sum(delta_hat != 0), "\n") +cat("Positions with |beta_hat| > 1:", + which(abs(beta_hat) > 1), "\n") +``` + +The sparse component correctly locates the five large effects — those with `|beta_hat| > 1` match the true positions exactly. The dense component is non-zero everywhere, absorbing the background of small signals. + +## Comparing estimators + +We compare LAVA against the Lasso and Elastic-net at their BIC-selected models. In the LAVA panel, bars are split into two layers drawn with `geom_col`: the **total** estimate $\hat\theta = \hat\beta + \hat\delta$ is drawn first (pink, dense component), and the **sparse** component $\hat\beta$ is overlaid on top (green). The remaining pink portion at the tips of active bars represents the dense contribution $\hat\delta$; the small pink bars at inactive positions are the dense signal alone. + +```{r comparison-fit} +fit_lasso <- lasso(X, y) +fit_enet <- elastic_net(X, y, lambda2 = 1) +``` + +```{r comparison-plot, fig.height = 6} +get_coef <- function(fit) { + b <- fit$get_model("BIC") + if ("Intercept" %in% names(b)) b <- b[-1] + as.numeric(b) +} + +method_levels <- c("True signal", "Lasso (BIC)", "Elastic-net (BIC)", "LAVA (BIC)") + +## Non-LAVA panels: one bar per predictor +df_base <- rbind( + data.frame(method = "True signal", index = 1:p, + value = theta, fill_grp = "True signal"), + data.frame(method = "Lasso (BIC)", index = 1:p, + value = get_coef(fit_lasso), fill_grp = "Lasso"), + data.frame(method = "Elastic-net (BIC)", index = 1:p, + value = get_coef(fit_enet), fill_grp = "Elastic-net") +) +df_base$method <- factor(df_base$method, levels = method_levels) + +## LAVA panel: dense total drawn first, sparse component overlaid on top +df_lava_dense <- data.frame( + method = factor("LAVA (BIC)", levels = method_levels), + index = 1:p, + value = as.numeric(theta_hat), + fill_grp = "LAVA dense (δ̂)" +) +df_lava_sparse <- data.frame( + method = factor("LAVA (BIC)", levels = method_levels), + index = 1:p, + value = beta_hat, + fill_grp = "LAVA sparse (β̂)" +) + +fill_palette <- c( + "True signal" = "grey20", + "Lasso" = "#56B4E9", + "Elastic-net" = "#E69F00", + "LAVA dense (δ̂)" = "#CC79A7", + "LAVA sparse (β̂)" = "#009E73" +) + +ggplot(mapping = aes(x = index, y = value, fill = fill_grp)) + + geom_col(data = df_base) + + geom_col(data = df_lava_dense) + + geom_col(data = df_lava_sparse) + + facet_wrap(~ method, ncol = 2) + + scale_fill_manual( + values = fill_palette, + breaks = c("LAVA sparse (β̂)", "LAVA dense (δ̂)"), + name = "LAVA decomposition" + ) + + labs(x = "Predictor index", y = "Estimated coefficient", + title = "Coefficient recovery: Lasso vs. Elastic-net vs. LAVA") + + theme_bw() + + theme(legend.position = "bottom") +``` + +- The **Lasso** selects a handful of predictors but cannot represent the dense background. +- The **Elastic-net** spreads a small weight across more predictors, but conflates the two structural components. +- **LAVA** explicitly decomposes the signal: the green bars identify the large sparse effects precisely; the pink portions (small rims at active positions, full bars at inactive positions) capture the diffuse dense background. + +## Cross-validation + +`$cross_validate()` works identically to the other `QuadrupenFit` subclasses. + +```{r cv, message = FALSE} +set.seed(42) +fit_lava$cross_validate(K = 10, verbose = FALSE) +fit_lava$plot(type = "crossval") +``` + +```{r cv-model} +set.seed(42) +fit_lasso$cross_validate(K = 10, verbose = FALSE) +fit_enet$cross_validate(K = 10, verbose = FALSE) + +idx_cv <- fit_lava$get_model("CV_min", type = "index") +beta_cv <- as.numeric(fit_lava$sparse_coef[, idx_cv]) +cat("Non-zero entries in sparse component (CV_min):", + sum(beta_cv != 0), "\n") +cat("Positions with |beta_cv| > 1:", + which(abs(beta_cv) > 1), "\n") +``` + +--- + +## Group LAVA: group-sparse + dense decomposition + +`group_lava()` extends LAVA by replacing the element-wise $\ell_1$ penalty on $\beta$ with a group-wise mixed norm $\Omega_g(\beta)$: + +$$\min_{\theta = \beta + \delta} \frac{1}{2n} \|y - X(\beta + \delta)\|^2 + \lambda_1\, \Omega_g(\beta) + \frac{\lambda_2}{2}\, \delta^\top S\, \delta.$$ + +This is appropriate when the sparse component has **group structure**: entire groups of predictors are either active or silent. The `type` argument selects the group norm ($\ell_1/\ell_2$ by default, or $\ell_1/\ell_\infty$). + +### Simulation with a group-sparse signal + +```{r group-simulation} +## 5 groups of 10 predictors; groups 2 and 4 are active in the sparse part +group <- rep(1:5, each = 10) +beta_g <- rep(c(0, 2, 0, -2, 0), each = 10) +delta_g <- rnorm(p, mean = 0, sd = 0.2) +theta_g <- beta_g + delta_g + +y2 <- X %*% theta_g + rnorm(n, sd = 2) + +group_names <- paste0("grp", 1:5) +var_labels <- group_names[group] +``` + +### Fitting + +```{r group-lava-fit} +fit_gl <- group_lava(X, y2, group = group, lambda2 = 1) +fit_gl +``` + +### Paths of the sparse and dense components + +```{r group-paths, fig.show='hold', out.width='50%'} +fit_gl$plot_path(component = "sparse", labels = var_labels) +fit_gl$plot_path(component = "dense", labels = var_labels) +``` + +The sparse path shows entire groups entering the model simultaneously (Group Lasso behaviour), while the dense path remains diffuse and non-zero everywhere. + +### Model selection + +For group selection problems, eBIC (which carries a stronger penalty of $\log(n) + 2\log(p)$) is generally more conservative than BIC and better adapted to identifying which groups are truly active. + +### Group identification + +```{r group-extraction} +idx_ebic <- fit_gl$get_model("eBIC", type = "index") +beta_g_hat <- as.numeric(fit_gl$sparse_coef[, idx_ebic]) + +active_groups <- unique(group[beta_g_hat != 0]) +cat("Active groups in sparse component (eBIC):", group_names[active_groups], "\n") +``` + +--- + +## Reference + +Chernozhukov, V., Hansen, C., & Liao, Y. (2017). A lava attack on the recovery of sums of dense and sparse signals. *The Annals of Statistics*, 45(1), 39–76. diff --git a/vignettes/sparse-regression.Rmd b/vignettes/sparse-regression.Rmd new file mode 100644 index 0000000..1ff0084 --- /dev/null +++ b/vignettes/sparse-regression.Rmd @@ -0,0 +1,219 @@ +--- +title: "Sparse Linear Regression with Quadrupen" +subtitle: "A tour of standard analysis (fit, cross-validation, information criteria, stability) " +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{quadrupen-basics} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.width = 7, + fig.height = 4 +) +``` + +## Setup + +```{r setup} +library(quadrupen) +data("Birthwt", package = "grpreg") +y <- Birthwt$bwt[-130] ## outlier +X <- Birthwt$X[-130, ] +``` + +The `Birthwt` dataset contains birth weight data with `r nrow(Birthwt$X)` observations and `r ncol(Birthwt$X)` predictors. We use it throughout this vignette to illustrate sparse linear regression. + +## Fitting sparse linear models + +The main entry point is `sparse_lm()`, which fits a regularization path for several sparse penalties: `"l1"` (Lasso/Elastic-net), `"mcp"` (Minimax Concave Penalty) and `"scad"` (Smoothly Clipped Absolute Deviation). Convenience wrappers `lasso()`, `mcp()`, `scad()` and `elastic_net()` are also available. + +### Lasso + +```{r lasso-fit} +fit_lasso <- lasso(X, y, minratio = 1e-2) +fit_lasso +``` + +### SCAD + +SCAD is a non-convex penalty that matches the Lasso near the origin and tapers to a constant for large coefficients, thereby eliminating shrinkage bias for strong signals. The shape parameter `eta` must be greater than 2. + +```{r scad-fit} +fit_scad <- scad(X, y, eta = 4, minratio = 1e-2) +fit_scad +``` + +### MCP + +MCP applies a linearly diminishing marginal penalty that vanishes once a coefficient exceeds `eta × lambda1`, yielding near-unbiased estimates for large signals. The shape parameter `eta` must be greater than 1. + +```{r mcp-fit} +fit_mcp <- mcp(X, y, eta = 3, minratio = 1e-2) +fit_mcp +``` + +### Elastic-Net + +The Elastic-net combines an $\ell_1$ and an $\ell_2$ penalty, trading off sparsity against grouping of correlated predictors. The `lambda2` argument controls the ridge component. + +```{r enet-fit} +fit_enet <- elastic_net(X, y, lambda2 = 1, minratio = 1e-2) +fit_enet +``` + +### Fit objects + +All objects returned by `sparse_lm()` and its wrappers are R6 instances of the class `SparseFit`, inheriting from `QuadrupenFit`. These objects (see [`QuadrupenFit`]) store all data related to the fit, accessible directly via named fields (e.g., `fit_lasso$coefficients`, `fit_lasso$deviance`, `fit_lasso$degrees_freedom`; see `str(fit_lasso)` and the documentation). + +They also provide methods for visualizing and analyzing the fit, and for pursuing complementary analyses such as model selection, cross-validation, and stability selection. For users unfamiliar with R6 classes, S3 methods are exported for the most common operations (e.g., `plot(fit_lasso)` is equivalent to `fit_lasso$plot()`), but the R6 methods expose more options. + +## Regularization paths + +`$plot_path()` (or `plot(fit, type = "path")`) displays how coefficients evolve along the penalty path. + +```{r path-plots, fig.show='hold', out.width='50%'} +fit_lasso$plot_path(xvar="fraction", log_scale = FALSE) +fit_mcp$plot_path(labels = colnames(X)) +fit_scad$plot_path() +fit_enet$plot_path(standardize=FALSE) +``` + +## Information criteria + +`$criteria()` computes AIC, BIC, mBIC, eBIC and GCV from the estimated degrees of freedom, storing the result as an [`InformationCriteria`] R6 object in the field `$information_criteria` of the current fit. This method is called automatically during fitting (`sparse_lm()` or any wrapper) at no additional computational cost, so users rarely need to invoke it manually — the main exception being when a custom penalty term is desired. + +```{r criteria} +fit_lasso$criteria() +``` + +```{r criteria-plots, fig.show='hold', out.width='50%'} +fit_lasso$plot(type = "criteria") +``` + +For greater flexibility, use the `plot` method of the `InformationCriteria` object directly to select which criteria to display: + +```{r criteria-plots-2, fig.show='hold', out.width='50%'} +fit_lasso$information_criteria$plot(c("AIC", "BIC", "mBIC")) +fit_mcp$information_criteria$plot("GCV") +``` + +## Model extraction + +`$get_model()` extracts the coefficient vector selected by a given criterion. + +```{r get-model} +coef_lasso_bic <- fit_lasso$get_model("BIC") +coef_mcp_bic <- fit_mcp$get_model("BIC") + +cat("Non-zero Lasso coefficients (BIC):\n") +print(coef_lasso_bic[coef_lasso_bic != 0]) + +cat("\nNon-zero MCP coefficients (BIC):\n") +print(coef_mcp_bic[coef_mcp_bic != 0]) +``` + +An existing lambda value can also be passed directly: + +```{r get-model-lambda} +lambda_lasso <- fit_lasso$major_tuning +coef_lasso_lambda <- fit_lasso$get_model(lambda_lasso[18]) + +cat("Non-zero Lasso coefficients for lambda = ", round(lambda_lasso[18], 2), "\n") +print(coef_lasso_lambda[coef_lasso_lambda != 0]) +``` + +## Cross-validation + +`$cross_validate()` performs K-fold cross-validation over the penalty grid, storing the result as a [`CrossValidation`] R6 object in the field `$cross_validation` of the current fit. + +```{r cv, message=FALSE} +set.seed(42) +fit_lasso$cross_validate(K = 10, verbose = FALSE) +fit_mcp$cross_validate(K = 10, verbose = FALSE) +``` + +```{r cv-plots, fig.show='hold', out.width='50%'} +fit_lasso$plot(type = "crossval") +fit_mcp$plot(type = "crossval") +``` + +Model selection using the CV-minimizing penalty: + +```{r cv-model} +coef_lasso_cv <- fit_lasso$get_model("CV_min") +cat("Non-zero Lasso coefficients (CV min):\n") +print(coef_lasso_cv[coef_lasso_cv != 0]) +``` + +## Cross-validation on $\lambda_1 \times \lambda_2$ + +For regularizers with two penalties (typically the Elastic-net), cross-validation can be run on a two-dimensional grid: + +```{r cv2D, message=FALSE} +set.seed(42) +fit_enet$cross_validate(K = 5, lambda2 = 10^seq(1, -3, len=20), verbose = FALSE) +``` + +```{r cv-plots-2D, fig.show='hold', out.width='50%'} +fit_enet$plot(type = "crossval") +fit_enet$cross_validation$plotCV_1D(se = FALSE) + ggplot2::ylim(c(0.4, 0.55)) +``` + +## Prediction + +`$predict()` returns fitted values for the whole path, or for a specific model selected by a criterion or a lambda value. + +```{r predict} +## Predictions at the CV_min-selected model +y_hat_lasso <- fit_lasso$predict(selection = "CV_min") +y_hat_mcp <- fit_mcp$predict(selection = "CV_min") +y_hat_enet <- fit_enet$predict(selection = "CV_min") + +## R² at the selected model +r2_lasso <- fit_lasso$r_squared[fit_lasso$get_model("CV_min", type = "index")] +r2_mcp <- fit_mcp$r_squared[fit_mcp$get_model("CV_min", type = "index")] +r2_enet <- fit_enet$r_squared[fit_enet$get_model("CV_min", type = "index")] +cat(sprintf("R² Lasso (CV_min): %.3f\nR² MCP (CV_min): %.3f\nR² Enet (CV_min): %.3f\n", r2_lasso, r2_mcp, r2_enet)) +``` + +## Stability selection + +`$stability()` estimates selection probabilities via repeated sub-sampling (Meinshausen & Bühlmann, 2010). Variables that appear consistently across sub-samples receive high probabilities and are considered robustly selected. The stability path is stored as a [`StabilityPath`] R6 object in the field `$stability_path` of the current fit. + +```{r stability, message=FALSE} +set.seed(42) +fit_lasso$stability(n_subsamples = 200, verbose = FALSE) +``` + +```{r stability-plot} +fit_lasso$plot(type = "stability") +fit_lasso$stability_path$plot(nvarsel = 7) +colnames(X)[fit_lasso$stability_path$selection(nvarsel = 7) ] +``` + +## Debiased Estimators + +Every $\ell_1$-based regularizer is known to be biased due to the shrinkage effect. A standard remedy is to refit the model on the selected support without penalization. This is implemented in **quadrupen** at no extra computational cost — the refitted coefficients are computed alongside the regularization path — and can be activated by setting the `$debias` field to `TRUE`. + + +```{r debias} +fit_lasso$debias <- TRUE +fit_lasso$plot_path() +``` + +Since both versions are stored in the same object, it is straightforward to compare the regularized and debiased solutions — for the path as well as for CV error — by toggling `$debias`. + +```{r debias-cv} +fit_lasso$cross_validate(verbose = FALSE) +fit_lasso$plot(type = "crossval") +fit_lasso$debias <- FALSE +fit_lasso$plot_path() +``` + +Setting `$debias` back to `FALSE` restores the original regularized coefficients. + diff --git a/vignettes/structured-signal-recovery.Rmd b/vignettes/structured-signal-recovery.Rmd new file mode 100644 index 0000000..1b38fd5 --- /dev/null +++ b/vignettes/structured-signal-recovery.Rmd @@ -0,0 +1,292 @@ +--- +title: "Recovering a Structured Signal with Quadrupen" +subtitle: "From the Lasso to structured regularization via graph Laplacians" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{structured-signal-recovery} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.width = 7, + fig.height = 4 +) +``` + +## Motivation + +Sparse penalties such as the Lasso treat all predictors symmetrically: coefficients are shrunk independently without any notion of which predictors are related. In many applications, however, the signal has a known structure — predictors are grouped, ordered, or connected through a graph — and penalties that exploit this structure can substantially improve estimation. + +This vignette builds a progressively richer family of penalties on a single simulation, showing how incorporating structural knowledge step by step leads to better coefficient recovery. + +## Setup + +```{r setup, message = FALSE} +library(quadrupen) +library(Matrix) +library(ggplot2) +``` + +### Simulation + +We simulate a linear model with $p = 95$ predictors and $n = 50$ observations. The true coefficient vector $\beta$ is **piecewise constant**: it is zero over three large segments (sizes 25, 25, 25) and takes values $+1$ and $-1$ over two smaller active blocks (sizes 10 each). + +```{r simulation} +set.seed(42) +p <- 95 +n <- 50 +beta <- rep(c(0, 1, 0, -1, 0), c(25, 10, 25, 10, 25)) + +## Block-correlated predictors: high within-segment, zero across segments +rho <- 0.75 +Soo <- toeplitz(rho^(0:24)) # Toeplitz correlation within zero segments +Sww <- matrix(rho, 10, 10) # constant correlation within active segments +Sigma <- bdiag(Soo, Sww, Soo, Sww, Soo) +diag(Sigma) <- 1 + +X <- as.matrix(matrix(rnorm(p * n), n, p) %*% chol(Sigma)) +y <- X %*% beta + rnorm(n, 0, 10) + +## Segment labels (for plot legends) +segments <- rep(1:5, c(25, 10, 25, 10, 25)) +segment_names <- paste0("seg", 1:5) +seg_labels <- segment_names[segments] +``` + +The five segments are clearly visible in the true signal: + +```{r plot-beta, fig.height = 3} +data.frame(index = seq_along(beta), beta = beta, segment = seg_labels) |> + ggplot(aes(x = index, y = beta, colour = segment)) + + geom_step() + geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.4) + + labs(x = "Predictor index", y = expression(beta), title = "True coefficient vector") + + theme_bw() +``` + +### Community graph and its Laplacian + +The structural prior we wish to encode is: *predictors within the same segment should have similar coefficients*. This is captured by a **community graph** whose edges connect all pairs of predictors belonging to the same segment (a union of disjoint cliques). + +```{r graph} +## Adjacency matrix: clique within each segment, no edges across +Ioo <- matrix(1, 25, 25) +Iww <- matrix(1, 10, 10) +A <- as.matrix(bdiag(Ioo, Iww, Ioo, Iww, Ioo)) +diag(A) <- 0 + +## Graph Laplacian: L = D - A (regularized for positive definiteness) +L2 <- -A +diag(L2) <- colSums(A) + 1e-2 +``` + +The Laplacian $L_2$ defines the quadratic form + +$$\beta^\top L_2\, \beta = \sum_{(i,j)\in E} (\beta_i - \beta_j)^2,$$ + +which penalizes *differences between connected predictors*. Adding $\lambda_2 \beta^\top L_2 \beta / 2$ to the loss therefore encourages within-segment coefficients to be equal — a perfect match for the piecewise-constant signal. Note that $L_2$ is passed as the `struct` argument to `ridge()` and `elastic_net()`. + +For the Fused Lasso, the graph structure is encoded differently: `struct` takes the **adjacency matrix** $A$ directly, since the fusion penalty is + +$$\lambda_2 \sum_{(i,j)\in E} |\beta_i - \beta_j|,$$ + +an $\ell_1$ (not $\ell_2$) penalty on pairwise differences. The graph topology is the same, but the penalty form differs. + +```{r graph-viz, fig.height = 4, fig.width = 5, fig.align = "center"} +if (requireNamespace("igraph", quietly = TRUE)) { + g <- igraph::graph_from_adjacency_matrix(A, mode = "undirected") + cols <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#CC79A7") + vertex_colors <- cols[segments] + igraph::plot.igraph( + g, vertex.color = vertex_colors, vertex.size = 3, + vertex.label = NA, edge.width = 0.1, + main = "Community graph (colour = segment)" + ) + legend("topleft", legend = segment_names, fill = cols, bty = "n", cex = 0.8) +} +``` + +--- + +## Lasso: a baseline without structural knowledge + +The Lasso applies an element-wise $\ell_1$ penalty and treats all predictors independently. With highly correlated predictors it tends to select one representative per block arbitrarily, producing an erratic path. + +```{r lasso} +fit_lasso <- lasso(X, y, intercept = FALSE) +fit_lasso +``` + +```{r lasso-path} +fit_lasso$plot_path(labels = seg_labels) +``` + +--- + +## Fused Lasso: penalizing differences between consecutive predictors + +The Fused Lasso adds an $\ell_1$ penalty on *differences between adjacent coefficients* on top of the standard $\ell_1$ penalty: + +$$\min_\beta \frac{1}{2}\|y - X\beta\|^2 + \lambda_1 \|\beta\|_1 + \lambda_2 \sum_{i=1}^{p-1} |\beta_{i+1} - \beta_i|.$$ + +By default, `fused_lasso()` uses a **chain graph** (consecutive pairs), which is appropriate for ordered or spatial signals. + +```{r fused-chain} +fit_fused_chain <- fused_lasso(X, y, lambda2 = 5, intercept = FALSE) +fit_fused_chain +``` + +```{r fused-chain-path} +fit_fused_chain$plot_path(labels = seg_labels) +``` + +### Fused Lasso with the community graph + +Replacing the chain by the community adjacency matrix $A$ fuses predictors *within the same segment* rather than consecutive ones. This directly targets the true block structure of $\beta$. + +```{r fused-community} +fit_fused_comm <- fused_lasso(X, y, lambda2 = 5, struct = A, intercept = FALSE) +fit_fused_comm +``` + +```{r fused-community-path} +fit_fused_comm$plot_path(labels = seg_labels) +``` + +--- + +## Structured Ridge: shrinking toward block-constant solutions + +The standard ridge regression shrinks all coefficients toward **zero**. The structured ridge replaces the identity by a positive semidefinite matrix $S$: + +$$\min_\beta \frac{1}{2}\|y - X\beta\|^2 + \frac{\lambda_2}{2}\, \beta^\top S\, \beta.$$ + +With $S = L_2$ (the community Laplacian), the penalty becomes $\frac{\lambda_2}{2}\beta^\top L_2 \beta$, which shrinks *within-segment differences* rather than the coefficients themselves. + +### Standard ridge + +```{r ridge-std} +fit_ridge_std <- ridge(X, y, intercept = FALSE) +fit_ridge_std$plot_path(labels = seg_labels) +``` + +The path is a blur: all predictors are shrunk toward zero with no differentiation between segments. + +### Structured ridge ($S = L_2$) + +```{r ridge-struct} +fit_ridge_struct <- ridge(X, y, struct = L2, lambda_max = 1000, intercept = FALSE) +fit_ridge_struct$plot_path(labels = seg_labels) +``` + +The Laplacian prior forces coefficients within the same segment to track each other. The two active segments (2 and 4) now emerge as coherent bundles. + +--- + +## Structured Elastic-net: sparsity + graph-Laplacian regularization + +The Elastic-net combines an $\ell_1$ penalty with a quadratic penalty: + +$$\min_\beta \frac{1}{2}\|y - X\beta\|^2 + \lambda_1\|\beta\|_1 + \frac{\lambda_2}{2}\, \beta^\top S\, \beta.$$ + +Setting $S = L_2$ yields the **structured Elastic-net**: the $\ell_1$ term drives entire segments to zero (sparsity at the segment level), while the Laplacian term forces the non-zero segments to be internally smooth. + +```{r enet-struct} +fit_enet_struct <- elastic_net(X, y, lambda2 = 5, struct = L2, intercept = FALSE) +fit_enet_struct +``` + +```{r enet-struct-path} +fit_enet_struct$plot_path(labels = seg_labels) +``` + +The two active segments now enter the path as well-separated, near-constant bundles while the three zero segments stay at zero. + +### Debiased structured Elastic-net + +The $\ell_1$ penalty introduces shrinkage bias. Refitting the model on the selected support without penalization recovers near-unbiased estimates. + +```{r enet-debias} +fit_enet_struct$debias <- TRUE +fit_enet_struct$plot_path(labels = seg_labels) +fit_enet_struct$debias <- FALSE +``` + +--- + +## Coefficient recovery: a side-by-side comparison + +We extract the BIC-selected coefficient vector for each method and compare visually against the true $\beta$. + +```{r recovery} +methods <- list( + "Lasso" = fit_lasso, + "Fused Lasso (chain)" = fit_fused_chain, + "Fused Lasso (community)"= fit_fused_comm, + "Ridge (standard)" = fit_ridge_std, + "Ridge (structured)" = fit_ridge_struct, + "Elastic-net (structured)"= fit_enet_struct +) + +extract_coef <- function(fit) { + b <- fit$get_model("BIC") + if ("Intercept" %in% names(b)) b <- b[-1] # drop intercept + b +} + +df_coef <- do.call(rbind, lapply(names(methods), function(nm) { + b <- extract_coef(methods[[nm]]) + data.frame( + method = factor(nm, levels = names(methods)), + index = seq_along(b), + value = as.numeric(b), + segment = seg_labels + ) +})) + +df_true <- data.frame( + index = seq_along(beta), + beta = beta, + segment = seg_labels +) +``` + +```{r recovery-plot, fig.height = 9, fig.width = 7} +## Background shading for each segment +seg_bounds <- data.frame( + xmin = c(1, 26, 36, 61, 71), + xmax = c(25, 35, 60, 70, 95), + fill = factor(1:5) +) + +ggplot(df_coef, aes(x = index, y = value)) + + geom_rect(data = seg_bounds, + aes(xmin = xmin, xmax = xmax, ymin = -Inf, ymax = Inf, fill = fill), + alpha = 0.08, inherit.aes = FALSE) + + scale_fill_manual(values = c("#E69F00","#56B4E9","#009E73","#F0E442","#CC79A7"), + guide = "none") + + geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.3) + + geom_step(data = df_true, aes(y = beta), colour = "black", linewidth = 0.8, + linetype = "dashed") + + geom_point(aes(colour = segment), size = 0.8, alpha = 0.7) + + scale_colour_manual(values = c("#E69F00","#56B4E9","#009E73","#F0E442","#CC79A7")) + + facet_wrap(~ method, ncol = 2) + + labs(x = "Predictor index", y = "Estimated coefficient", + title = "Coefficient recovery by method", + subtitle = "Dashed line: true β", + colour = "Segment") + + theme_bw() + + theme(legend.position = "bottom") +``` + +The progression illustrates the benefit of structural priors: + +- **Lasso**: selects a few predictors erratically within the active segments. +- **Fused Lasso (chain)**: smoother path, but the chain graph blurs transitions between segments rather than reinforcing block structure. +- **Fused Lasso (community)**: within-segment fusion correctly pools the active blocks, but the $\ell_1$ fusion penalty may still miss some predictors. +- **Ridge (standard)**: all coefficients shrunk toward zero; the block pattern is not recovered. +- **Ridge (structured)**: within-segment coefficients cluster together; the active segments are identifiable but not zeroed out. +- **Elastic-net (structured)**: sparsity zeros out the inactive segments while the Laplacian smooths the active ones — the closest recovery to the true signal.