Skip to content

Commit 5a0884b

Browse files
jonathan-taylorjonathan-taylor
authored
Merge pull request #28 from jonathan-taylor/export_lasso_poly
exporting function to construct polyhedral constraints of LASSO
2 parents 0a3bc8d + f0010f3 commit 5a0884b

File tree

4 files changed

+122
-47
lines changed

4 files changed

+122
-47
lines changed

‎forLater/fixedLassoPoly.Rd

Lines changed: 81 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,81 @@
1+
\name{fixedLassoPoly}
2+
\alias{fixedLassoPoly}
3+
4+
\title{
5+
Compute polyhedral constraints for a LASSO problem with
6+
a fixed value of lambda.
7+
}
8+
\description{
9+
Compute polyhedral representation of the selection region of Lee et al. (2016).
10+
By construction, y should satisfy A %*% y elementwise less then or equal b.
11+
}
12+
\usage{
13+
fixedLassoPoly(X, y, lambda, beta, active, inactive = FALSE)
14+
}
15+
\arguments{
16+
\item{X}{
17+
Design matrix of LASSO problem.
18+
}
19+
\item{y}{
20+
Response of LASSO problem.
21+
}
22+
\item{lambda}{
23+
Value of regularization parameter.
24+
}
25+
\item{beta}{
26+
Solution of LASSO problem with regularization parameter set to lambda.
27+
}
28+
\item{active}{
29+
Active set of the LASSO problem as a boolean vector. Should correspond
30+
to the non-zeros of beta.
31+
}
32+
\item{inactive}{
33+
Form the inactive constraints as well?
34+
}
35+
}
36+
\details{
37+
This function computes
38+
the polyhedral representation of the selection region of Lee et al. (2016).
39+
}
40+
41+
\value{
42+
\item{A}{Linear part of the affine inequalities.}
43+
\item{b}{RHS offset the affine inequalities.}
44+
}
45+
46+
\references{
47+
Jason Lee, Dennis Sun, Yuekai Sun, and Jonathan Taylor (2016).
48+
Exact post-selection inference, with application to the lasso. Annals of Statistics, 44(3), 907-927.
49+
50+
Jonathan Taylor and Robert Tibshirani (2017) Post-selection inference for math L1-penalized likelihood models.
51+
Canadian Journal of Statistics, xx, 1-21. (Volume still not posted)
52+
}
53+
\author{Ryan Tibshirani, Rob Tibshirani, Jonathan Taylor, Joshua Loftus, Stephen Reid}
54+
55+
\examples{
56+
57+
set.seed(43)
58+
n = 50
59+
p = 10
60+
sigma = 1
61+
62+
x = matrix(rnorm(n*p),n,p)
63+
x = scale(x,TRUE,TRUE)
64+
65+
beta = c(3,2,rep(0,p-2))
66+
y = x\%*\%beta + sigma*rnorm(n)
67+
68+
# first run glmnet
69+
gfit = glmnet(x,y,standardize=FALSE)
70+
71+
# extract coef for a given lambda; note the 1/n factor!
72+
# (and we don't save the intercept term)
73+
lambda = .8
74+
beta = coef(gfit, s=lambda/n, exact=TRUE)[-1]
75+
active = (beta != 0)
76+
77+
fixedLassoPoly(x, y, lambda, beta, active)
78+
fixedLassoPoly(x, y, lambda, beta, active, inactive=TRUE)
79+
80+
}
81+

‎selectiveInference/R/funs.fixed.R

Lines changed: 29 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -93,14 +93,14 @@ fixedLassoInf <- function(x, y, beta, lambda, family=c("gaussian","binomial","co
9393
"'thresh' parameter, for a more accurate convergence."))
9494

9595
# Get lasso polyhedral region, of form Gy >= u
96-
if (type == 'full' & p > n) out = fixedLasso.poly(x,y,beta,lambda,vars,inactive=TRUE)
97-
else out = fixedLasso.poly(x,y,beta,lambda,vars)
98-
G = out$G
99-
u = out$u
96+
if (type == 'full' & p > n) out = fixedLassoPoly(x,y,lambda,beta,vars,inactive=TRUE)
97+
else out = fixedLassoPoly(x,y,lambda,beta,vars)
98+
A = out$A
99+
b = out$b
100100

101101
# Check polyhedral region
102102
tol.poly = 0.01
103-
if (min(G %*% y - u) < -tol.poly * sqrt(sum(y^2)))
103+
if (max(A %*% y - b) > tol.poly * sqrt(sum(y^2)))
104104
stop(paste("Polyhedral constraints not satisfied; you must recompute beta",
105105
"more accurately. With glmnet, make sure to use exact=TRUE in coef(),",
106106
"and check whether the specified value of lambda is too small",
@@ -191,7 +191,7 @@ fixedLassoInf <- function(x, y, beta, lambda, family=c("gaussian","binomial","co
191191
sign[j] = sign(sum(vj*y))
192192
vj = sign[j] * vj
193193

194-
limits.info = TG.limits(y, -G, -u, vj, Sigma=diag(rep(sigma^2, n)))
194+
limits.info = TG.limits(y, A, b, vj, Sigma=diag(rep(sigma^2, n)))
195195
a = TG.pvalue.base(limits.info, null_value=null_value[j], bits=bits)
196196
pv[j] = a$pv
197197
vlo[j] = a$vlo * mj # Unstandardize (mult by norm of vj)
@@ -221,45 +221,39 @@ fixedLassoInf <- function(x, y, beta, lambda, family=c("gaussian","binomial","co
221221
#############################
222222

223223

224-
fixedLasso.poly=
225-
function(x, y, beta, lambda, a, inactive = FALSE) {
226-
xa = x[,a,drop=F]
227-
xac = x[,!a,drop=F]
228-
xai = pinv(crossprod(xa))
229-
xap = xai %*% t(xa)
230-
za = sign(beta[a])
224+
fixedLassoPoly =
225+
function(X, y, lambda, beta, active, inactive = FALSE) {
226+
Xa = X[,active,drop=F]
227+
Xac = X[,!active,drop=F]
228+
Xai = pinv(crossprod(Xa))
229+
Xap = Xai %*% t(Xa)
230+
231+
za = sign(beta[active])
231232
if (length(za)>1) dz = diag(za)
232233
if (length(za)==1) dz = matrix(za,1,1)
233234

234-
if (inactive) {
235-
P = diag(1,nrow(xa)) - xa %*% xap
235+
if (inactive) { # should we include the inactive constraints?
236+
R = diag(1,nrow(Xa)) - Xa %*% Xap # R is residual forming matrix of selected model
236237

237-
G = -rbind(
238-
1/lambda * t(xac) %*% P,
239-
-1/lambda * t(xac) %*% P,
240-
-dz %*% xap
238+
A = rbind(
239+
1/lambda * t(Xac) %*% R,
240+
-1/lambda * t(Xac) %*% R,
241+
-dz %*% Xap
241242
)
242243
lambda2=lambda
243-
if(length(lambda)>1) lambda2=lambda[a]
244-
u = -c(
245-
1 - t(xac) %*% t(xap) %*% za,
246-
1 + t(xac) %*% t(xap) %*% za,
247-
-lambda2 * dz %*% xai %*% za)
244+
if(length(lambda)>1) lambda2=lambda[active]
245+
b = c(
246+
1 - t(Xac) %*% t(Xap) %*% za,
247+
1 + t(Xac) %*% t(Xap) %*% za,
248+
-lambda2 * dz %*% Xai %*% za)
248249
} else {
249-
G = -rbind(
250-
# 1/lambda * t(xac) %*% P,
251-
# -1/lambda * t(xac) %*% P,
252-
-dz %*% xap
253-
)
250+
A = -dz %*% Xap
254251
lambda2=lambda
255-
if(length(lambda)>1) lambda2=lambda[a]
256-
u = -c(
257-
# 1 - t(xac) %*% t(xap) %*% za,
258-
# 1 + t(xac) %*% t(xap) %*% za,
259-
-lambda2 * dz %*% xai %*% za)
252+
if(length(lambda)>1) lambda2=lambda[active]
253+
b = -lambda2 * dz %*% Xai %*% za
260254
}
261255

262-
return(list(G=G,u=u))
256+
return(list(A=A, b=b))
263257
}
264258

265259
##############################

‎selectiveInference/R/funs.fs.R

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -295,13 +295,13 @@ fsInf <- function(obj, sigma=NULL, alpha=0.1, k=NULL, type=c("active","all","aic
295295
for (j in 1:k) {
296296
if (verbose) cat(sprintf("Inference for variable %i ...\n",vars[j]))
297297

298-
Gj = G[1:nconstraint[j],]
299-
uj = rep(0,nconstraint[j])
298+
Aj = -G[1:nconstraint[j],]
299+
bj = -rep(0,nconstraint[j])
300300
vj = vreg[j,]
301301
mj = sqrt(sum(vj^2))
302302
vj = vj / mj # Standardize (divide by norm of vj)
303303

304-
limits.info = TG.limits(y, -Gj, -uj, vj, Sigma=diag(rep(sigma^2, n)))
304+
limits.info = TG.limits(y, Aj, bj, vj, Sigma=diag(rep(sigma^2, n)))
305305
a = TG.pvalue.base(limits.info, bits=bits)
306306

307307
pv[j] = a$pv
@@ -353,10 +353,10 @@ fsInf <- function(obj, sigma=NULL, alpha=0.1, k=NULL, type=c("active","all","aic
353353
vj = vj / mj # Standardize (divide by norm of vj)
354354
sign[j] = sign(sum(vj*y))
355355
vj = sign[j] * vj
356-
Gj = rbind(G,vj)
357-
uj = c(u,0)
356+
Aj = -rbind(G,vj)
357+
bj = -c(u,0)
358358

359-
limits.info = TG.limits(y, -Gj, -uj, vj, Sigma=diag(rep(sigma^2, n)))
359+
limits.info = TG.limits(y, Aj, bj, vj, Sigma=diag(rep(sigma^2, n)))
360360
a = TG.pvalue.base(limits.info, bits=bits)
361361
pv[j] = a$pv
362362
sxj = sx[vars[j]]

‎selectiveInference/R/funs.lar.R

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -367,13 +367,13 @@ larInf <- function(obj, sigma=NULL, alpha=0.1, k=NULL, type=c("active","all","ai
367367
for (j in 1:k) {
368368
if (verbose) cat(sprintf("Inference for variable %i ...\n",vars[j]))
369369

370-
Gj = G[1:nk[j],]
371-
uj = rep(0,nk[j])
370+
Aj = -G[1:nk[j],]
371+
bj = -rep(0,nk[j])
372372
vj = vreg[j,]
373373
mj = sqrt(sum(vj^2))
374374
vj = vj / mj # Standardize (divide by norm of vj)
375375

376-
limits.info = TG.limits(y, -Gj, -uj, vj, Sigma=diag(rep(sigma^2, n)))
376+
limits.info = TG.limits(y, Aj, bj, vj, Sigma=diag(rep(sigma^2, n)))
377377
a = TG.pvalue.base(limits.info, bits=bits)
378378
pv[j] = a$pv
379379
sxj = sx[vars[j]]
@@ -428,10 +428,10 @@ larInf <- function(obj, sigma=NULL, alpha=0.1, k=NULL, type=c("active","all","ai
428428
vj = vj / mj # Standardize (divide by norm of vj)
429429
sign[j] = sign(sum(vj*y))
430430
vj = sign[j] * vj
431-
Gj = rbind(G,vj)
432-
uj = c(u,0)
431+
Aj = -rbind(G,vj)
432+
bj = -c(u,0)
433433

434-
limits.info = TG.limits(y, -Gj, -uj, vj, Sigma=diag(rep(sigma^2, n)))
434+
limits.info = TG.limits(y, Aj, bj, vj, Sigma=diag(rep(sigma^2, n)))
435435
a = TG.pvalue.base(limits.info, bits=bits)
436436

437437
pv[j] = a$pv

0 commit comments

Comments
 (0)