R: Linearized Bregman solver for linear, binomial, multinomial...
lb
R Documentation
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.
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.
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.
#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)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(Libra)
Loading required package: nnls
Loaded Libra 1.5
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/Libra/lb.Rd_%03d_medium.png", width=480, height=480)
> ### Name: lb
> ### Title: Linearized Bregman solver for linear, binomial, multinomial
> ### models with lasso, group lasso or column lasso penalty.
> ### Aliases: lb
> ### Keywords: regression
>
> ### ** Examples
>
> #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)
>
>
>
>
>
>
> dev.off()
null device
1
>