Libra (version 1.7)

lb: Linearized Bregman solver for linear, binomial, multinomial models with lasso, group lasso or column lasso penalty.

Description

Solver for the entire solution path of coefficients for Linear Bregman iteration.

Usage

lb(
  X,
  y,
  kappa,
  alpha,
  c = 1,
  tlist,
  nt = 100,
  trate = 100,
  family = c("gaussian", "binomial", "multinomial"),
  group = FALSE,
  index = NA,
  intercept = TRUE,
  normalize = TRUE,
  print = FALSE
)

Arguments

X

An n-by-p matrix of predictors

y

Response Variable

kappa

The damping factor of the Linearized Bregman Algorithm that is defined in the reference paper. See details.

alpha

Parameter in Linearized Bregman algorithm which controls the step-length of the discretized solver for the Bregman Inverse Scale Space. See details.

c

Normalized step-length. If alpha is missing, alpha is automatically generated by alpha=n*c/(kappa*||X^T*X||_2). It should be in (0,2) for family = "gaussian"(Default is 1), (0,8) for family = "binomial"(Default is 4), (0,4) for family = "multinomial"(Default is 2). If beyond these range the path may be oscillated at large t values.

tlist

Parameters t along the path.

nt

Number of t. Used only if tlist is missing. Default is 100.

trate

tmax/tmin. Used only if tlist is missing. Default is 100.

family

Response type

group

Whether to use a group penalty, Default is FALSE.

index

For group models, the index is a vector that determines the group of the parameters. Parameters of the same group should have equal value in index. Be careful that multinomial group model default assumes the variables in same column are in the same group, and a empty value of index means each variable is a group.

intercept

if TRUE, an intercept is included in the model (and not penalized), otherwise no intercept is included. Default is TRUE.

normalize

if TRUE, each variable is scaled to have L2 norm square-root n. Default is TRUE.

print

If TRUE, the percentage of finished computation is printed.

Value

A "lb" class object is returned. The list contains the call, the type, the path, the intercept term a0 and value for alpha, kappa, iter, and meanvalue, scale factor of X, meanx and normx. For gaussian and bonomial, path is a p-by-nt matrix, and for multinomial, path is a k-by-p-by-nt array, each dimension represents class, predictor and parameter t.

Details

The Linearized Bregman solver computes the whole regularization path for different types of lasso-penalty for gaussian, binomial and multinomial models through iterations. It is the Euler forward discretized form of the continuous Bregman Inverse Scale Space Differential Inclusion. For binomial models, the response variable y is assumed to be a vector of two classes which is transformed in to {1,-1}. For the multinomial models, the response variable y can be a vector of k classes or a n-by-k matrix that each entry is in {0,1} with 1 indicates the class. Under all circumstances, two parameters, kappa and alpha need to be specified beforehand. The definitions of kappa and alpha are the same as that defined in the reference paper. Parameter alpha is defined as stepsize and kappa is the damping factor of the Linearized Bregman Algorithm that is defined in the reference paper.

References

Ohser, Ruan, Xiong, Yao and Yin, Sparse Recovery via Differential Inclusions, https://arxiv.org/abs/1406.7728

Examples

# NOT RUN {
#Examples in the reference paper
library(MASS)
n = 80;p = 100;k = 30;sigma = 1
Sigma = 1/(3*p)*matrix(rep(1,p^2),p,p)
diag(Sigma) = 1
A = mvrnorm(n, rep(0, p), Sigma)
u_ref = rep(0,p)
supp_ref = 1:k
u_ref[supp_ref] = rnorm(k)
u_ref[supp_ref] = u_ref[supp_ref]+sign(u_ref[supp_ref])
b = as.vector(A%*%u_ref + sigma*rnorm(n))
kappa = 16
alpha = 1/160
object <- lb(A,b,kappa,alpha,family="gaussian",group=FALSE,
             trate=20,intercept=FALSE,normalize=FALSE)
plot(object,xlim=c(0,3),main=bquote(paste("LB ",kappa,"=",.(kappa))))


#Diabetes, linear case
library(Libra)
data(diabetes)
attach(diabetes)
object <- lb(x,y,100,1e-3,family="gaussian",group=FALSE)
plot(object)
detach(diabetes)

#Simulated data, binomial case
data('west10')
y<-2*west10[,1]-1;
X<-as.matrix(2*west10[,2:10]-1);
path <- lb(X,y,kappa = 1,family="binomial",trate=100,normalize = FALSE)
plot(path,xtype="norm",omit.zeros=FALSE)

#Simulated data, multinomial case
X <- matrix(rnorm(500*100), nrow=500, ncol=100)
alpha <- matrix(c(rnorm(30*3), rep(0,70*3)),nrow=3)
P <- exp(alpha%*%t(X))
P <- scale(P,FALSE,apply(P,2,sum))
y <- rep(0,500)
rd <- runif(500)
y[rd<P[1,]] <- 1
y[rd>1-P[3,]] <- -1
result <- lb(X,y,kappa=5,alpha=0.1,family="multinomial",
 group=TRUE,intercept=FALSE,normalize = FALSE)
plot(result)

# }