picasso: Fit Sparse Regression Models with Convex and Non-Convex...

View source: R/picasso.R

picassoR Documentation

Fit Sparse Regression Models with Convex and Non-Convex Penalties

Description

picasso fits regularization paths for sparse Gaussian, logistic (family = "binomial"), Poisson, and sqrt-lasso models using lasso, MCP, or SCAD penalties.

Usage

picasso(X, Y, lambda = NULL, nlambda = 100, lambda.min.ratio =
                 0.05, family = "gaussian", method = "l1",
                 type.gaussian = "naive", gamma = 3, df = NULL,
                 dfmax = NULL, standardize = TRUE, intercept = TRUE,
                 prec = 1e-07, max.ite = 1000, verbose = FALSE)

Arguments

X

Numeric design matrix with n rows (samples) and d columns (features).

Y

Response vector of length n. Use numeric responses for family = "gaussian" and family = "sqrtlasso", a two-level response for family = "binomial", and non-negative integer counts for family = "poisson".

lambda

Optional sequence of decreasing regularization values. If NULL, the sequence is generated automatically from nlambda and lambda.min.ratio.

nlambda

Number of regularization values to generate when lambda = NULL. The default is 100.

lambda.min.ratio

Smallest generated lambda value as a fraction of \lambda_{\max} when lambda = NULL. The default is 0.05. For non-convex penalties in logistic and Poisson models, values much smaller than 0.05 may lead to ill-conditioning.

family

Model family. Supported values are "gaussian", "binomial", "poisson", and "sqrtlasso". The default is "gaussian".

method

Penalty type. Supported values are "l1" (lasso), "mcp", and "scad". The default is "l1".

type.gaussian

Gaussian solver mode. "naive" maintains residuals and updates them in O(n) per coordinate; "covariance" pre-computes the Gram matrix and updates gradients in O(d) per coordinate, which is faster when n \gg d. The default is "naive".

gamma

Concavity parameter for MCP and SCAD penalties. The default is 3.

df

Reserved for backward compatibility with older Gaussian covariance-update implementations. It is currently not used by the active solver backend.

dfmax

Maximum number of nonzero coefficients for early stopping on the lambda path. When the number of nonzero coefficients exceeds dfmax, the solver stops and returns only the solutions computed so far. If NULL (default), no limit is imposed. In addition to this user-controlled limit, the solver automatically stops when the deviance ratio exceeds 0.999 or the relative deviance change falls below 1e-5 (checked after at least 5 lambda values).

standardize

If TRUE, standardize columns of X before fitting. The default is TRUE.

intercept

If TRUE, include an intercept term. The default is TRUE.

prec

Convergence tolerance. The default is 1e-7.

max.ite

Maximum number of iterations. The default is 1000.

verbose

If TRUE, print tracing information during fitting. The default is FALSE.

Details

When lambda is not supplied, picasso constructs a logarithmically spaced path between \lambda_{\max} and \lambda_{\min} = \mathrm{lambda.min.ratio} \cdot \lambda_{\max}.

The method solves a penalized optimization problem:

\min_{\beta_0, \beta}\; \mathcal{L}(\beta_0,\beta) + \sum_{j=1}^d p_{\lambda,\gamma}(|\beta_j|),

where p_{\lambda,\gamma} is lasso, MCP, or SCAD.

Loss functions by family are:

  • "gaussian": \mathcal{L}(\beta_0,\beta)=\frac{1}{2n}\lVert Y-\beta_0\mathbf{1}-X\beta \rVert_2^2.

  • "binomial": \mathcal{L}(\beta_0,\beta)=\frac{1}{n}\sum_{i=1}^n[\log(1+\exp(\beta_0+x_i^\top\beta))-y_i(\beta_0+x_i^\top\beta)].

  • "poisson": \mathcal{L}(\beta_0,\beta)=\frac{1}{n}\sum_{i=1}^n[\exp(\beta_0+x_i^\top\beta)-y_i(\beta_0+x_i^\top\beta)].

  • "sqrtlasso": \mathcal{L}(\beta_0,\beta)=\frac{1}{\sqrt{n}}\lVert Y-\beta_0\mathbf{1}-X\beta \rVert_2.

Value

A fitted object of class "gaussian", "logit", "poisson", or "sqrtlasso" depending on family.

Common components include:

beta

Sparse matrix of fitted coefficients. Columns correspond to regularization parameters in the solution path.

intercept

Intercept values along the solution path.

lambda

Regularization sequence used for fitting.

nlambda

Number of regularization values.

df

Degrees of freedom (number of nonzero coefficients) along the path.

method

Penalty type used for fitting.

ite

Iteration information returned by the solver.

verbose

Input verbose value.

runtime

Elapsed fitting time.

family

Input family value.

Author(s)

Jason Ge, Xingguo Li, Haoming Jiang, Mengdi Wang, Tong Zhang, Han Liu and Tuo Zhao
Maintainer: Tuo Zhao <tourzhao@gatech.edu>

References

1. J. Friedman, T. Hastie, H. Hofling, and R. Tibshirani. Pathwise coordinate optimization. The Annals of Applied Statistics, 2007.
2. C.-H. Zhang. Nearly unbiased variable selection under minimax concave penalty. The Annals of Statistics, 2010.
3. J. Fan and R. Li. Variable selection via nonconcave penalized likelihood and its oracle properties. Journal of the American Statistical Association, 2001.
4. R. Tibshirani, J. Bien, J. Friedman, T. Hastie, N. Simon, J. Taylor, and R. Tibshirani. Strong rules for discarding predictors in lasso-type problems. Journal of the Royal Statistical Society: Series B, 2012.
5. J. Ge, X. Li, H. Jiang, M. Wang, T. Zhang, H. Liu, and T. Zhao. PICASSO: A sparse learning library for high-dimensional data analysis in R and Python. arXiv preprint https://arxiv.org/abs/2006.15261.

See Also

picasso-package.

Examples

################################################################
## Sparse linear regression
## Generate the design matrix and regression coefficient vector
n = 100 # sample number 
d = 80 # sample dimension
c = 0.5 # correlation parameter
s = 20  # support size of coefficient
set.seed(2016)
X = scale(matrix(rnorm(n*d),n,d)+c*rnorm(n))/sqrt(n-1)*sqrt(n)
beta = c(runif(s), rep(0, d-s))

## Generate response using Gaussian noise, and fit sparse linear models
noise = rnorm(n)
Y = X%*%beta + noise

## l1 regularization solved with naive update
fitted.l1.naive = picasso(X, Y, nlambda=100, type.gaussian="naive")

## covariance update (faster when n >> d)
fitted.l1.covariance  = picasso(X, Y, nlambda=100, type.gaussian="covariance")

## early stopping: stop when more than 10 nonzero coefficients
fitted.l1.dfmax = picasso(X, Y, nlambda=100, dfmax=10)

## mcp regularization
fitted.mcp = picasso(X, Y, nlambda=100, method="mcp")

## scad regularization
fitted.scad = picasso(X, Y, nlambda=100, method="scad")

## lambdas used 
print(fitted.l1.naive$lambda)

## number of nonzero coefficients for each lambda
print(fitted.l1.naive$df)

## coefficients and intercept for the i-th lambda
i = 30
print(fitted.l1.naive$lambda[i])
print(fitted.l1.naive$beta[,i])
print(fitted.l1.naive$intercept[i])

## Visualize the solution path
plot(fitted.l1.naive)
plot(fitted.l1.covariance)
plot(fitted.mcp)
plot(fitted.scad)


################################################################
## Sparse logistic regression
## Generate the design matrix and regression coefficient vector
n <- 100  # sample number 
d <- 80   # sample dimension
c <- 0.5   # parameter controlling the correlation between columns of X
s <- 20    # support size of coefficient
set.seed(2016)
X <- scale(matrix(rnorm(n*d),n,d)+c*rnorm(n))/sqrt(n-1)*sqrt(n)
beta <- c(runif(s), rep(0, d-s))

## Generate response and fit sparse logistic models
p = 1/(1+exp(-X%*%beta))
Y = rbinom(n, rep(1,n), p)

## l1 regularization
fitted.l1 = picasso(X, Y, nlambda=100, family="binomial", method="l1")

## mcp regularization
fitted.mcp = picasso(X, Y, nlambda=100, family="binomial", method="mcp")

## scad regularization
fitted.scad = picasso(X, Y, nlambda=100, family="binomial", method="scad")

## lambdas used 
print(fitted.l1$lambda)

## number of nonzero coefficients for each lambda
print(fitted.l1$df)

## coefficients and intercept for the i-th lambda
i = 30
print(fitted.l1$lambda[i])
print(fitted.l1$beta[,i])
print(fitted.l1$intercept[i])

## Visualize the solution path
plot(fitted.l1)

## Estimate of Bernoulli parameters
param.l1 = fitted.l1$p


################################################################
## Sparse poisson regression
## Generate the design matrix and regression coefficient vector
n <- 100  # sample number 
d <- 80   # sample dimension
c <- 0.5   # parameter controlling the correlation between columns of X
s <- 20    # support size of coefficient
set.seed(2016)
X <- scale(matrix(rnorm(n*d),n,d)+c*rnorm(n))/sqrt(n-1)*sqrt(n)
beta <- c(runif(s), rep(0, d-s))/sqrt(s)

## Generate response and fit sparse poisson models
p = X%*%beta+rnorm(n)
Y = rpois(n, exp(p))

## l1 regularization
fitted.l1 = picasso(X, Y, nlambda=100, family="poisson", method="l1")

## mcp regularization
fitted.mcp = picasso(X, Y, nlambda=100, family="poisson", method="mcp")

## scad regularization
fitted.scad = picasso(X, Y, nlambda=100, family="poisson", method="scad")

## lambdas used 
print(fitted.l1$lambda)

## number of nonzero coefficients for each lambda
print(fitted.l1$df)

## coefficients and intercept for the i-th lambda
i = 30
print(fitted.l1$lambda[i])
print(fitted.l1$beta[,i])
print(fitted.l1$intercept[i])

## Visualize the solution path
plot(fitted.l1)

picasso documentation built on March 12, 2026, 5:06 p.m.