Skip to content

Commit 2bd33e8

Browse files
tibshiranitibshirani
authored
Merge pull request #18 from jonathan-taylor/master
comparing test
2 parents 1672cee + f6ae738 commit 2bd33e8

File tree

2 files changed

+228
-8
lines changed

2 files changed

+228
-8
lines changed

‎selectiveInference/R/funs.fixed.R

Lines changed: 42 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -159,7 +159,8 @@ sigma=NULL, alpha=0.1,
159159
hsigmaSinv <- solve(hsigmaS) # pinv(hsigmaS)
160160

161161
# Approximate inverse covariance matrix for when (n < p) from lasso_Inference.R
162-
htheta <- InverseLinfty(hsigma, n, length(S), verbose=FALSE)
162+
useC = TRUE
163+
htheta <- InverseLinfty(hsigma, n, length(S), verbose=FALSE, useC=useC)
163164
# htheta <- InverseLinfty(hsigma, n, verbose=FALSE)
164165

165166
FS = rbind(diag(length(S)),matrix(0,pp-length(S),length(S)))
@@ -269,7 +270,7 @@ fixedLasso.poly=
269270
### Functions borrowed and slightly modified from lasso_inference.R
270271

271272
## Approximates inverse covariance matrix theta
272-
InverseLinfty <- function(sigma, n, e, resol=1.5, mu=NULL, maxiter=50, threshold=1e-2, verbose = TRUE) {
273+
InverseLinfty <- function(sigma, n, e, resol=1.5, mu=NULL, maxiter=50, threshold=1e-2, verbose = TRUE, useC = FALSE) {
273274
# InverseLinfty <- function(sigma, n, resol=1.5, mu=NULL, maxiter=50, threshold=1e-2, verbose = TRUE) {
274275
isgiven <- 1;
275276
if (is.null(mu)){
@@ -294,7 +295,11 @@ InverseLinfty <- function(sigma, n, e, resol=1.5, mu=NULL, maxiter=50, threshold
294295
incr <- 0;
295296
while ((mu.stop != 1)&&(try.no<10)){
296297
last.beta <- beta
297-
output <- InverseLinftyOneRow(sigma, i, mu, maxiter=maxiter, threshold=threshold)
298+
if (useC == FALSE) {
299+
output <- InverseLinftyOneRow(sigma, i, mu, maxiter=maxiter, threshold=threshold)
300+
} else {
301+
output <- InverseLinftyOneRowC(sigma, i, mu, maxiter=maxiter)
302+
}
298303
beta <- output$optsol
299304
iter <- output$iter
300305
if (isgiven==1){
@@ -334,17 +339,46 @@ InverseLinfty <- function(sigma, n, e, resol=1.5, mu=NULL, maxiter=50, threshold
334339
return(M)
335340
}
336341

342+
InverseLinftyOneRowC <- function (Sigma, i, mu, maxiter=50) {
343+
344+
p = nrow(Sigma)
345+
basis_vector = rep(0, p)
346+
basis_vector[i] = 1.
347+
theta = rep(0, p)
348+
349+
val = .C("find_one_row",
350+
Sigma=as.double(Sigma),
351+
Sigma_diag=as.double(diag(Sigma)),
352+
Sigma_theta=as.double(rep(0, p)),
353+
nrow=as.integer(p),
354+
bound=as.double(mu),
355+
theta=as.double(theta),
356+
maxiter=as.integer(50),
357+
row=as.integer(i-1),
358+
coord=as.integer(i-1),
359+
dup=FALSE,
360+
package="selectiveInference")
361+
362+
# Check feasibility
363+
364+
if (max(abs(Sigma %*% val$theta - basis_vector)) > 1.01 * mu) {
365+
warning("Solution for row of M does not seem to be feasible")
366+
}
367+
368+
return(val$theta)
369+
}
370+
337371
InverseLinftyOneRow <- function ( sigma, i, mu, maxiter=50, threshold=1e-2 ) {
338372
p <- nrow(sigma);
339373
rho <- max(abs(sigma[i,-i])) / sigma[i,i];
340374
mu0 <- rho/(1+rho);
341375
beta <- rep(0,p);
342376

343-
if (mu >= mu0){
344-
beta[i] <- (1-mu0)/sigma[i,i];
345-
returnlist <- list("optsol" = beta, "iter" = 0);
346-
return(returnlist);
347-
}
377+
#if (mu >= mu0){
378+
# beta[i] <- (1-mu0)/sigma[i,i];
379+
# returnlist <- list("optsol" = beta, "iter" = 0);
380+
# return(returnlist);
381+
#}
348382

349383
diff.norm2 <- 1;
350384
last.norm2 <- 1;

‎selectiveInference/src/debiasing_matrix.c

Lines changed: 186 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,186 @@
1+
#include <stdio.h>
2+
#include <math.h> // for fabs
3+
4+
// Find an approximate row of \hat{Sigma}^{-1}
5+
6+
// Solves a dual version of problem (4) of https://arxiv.org/pdf/1306.3171.pdf
7+
8+
// Dual problem: \text{min}_{\theta} 1/2 \theta^T \Sigma \theta - e_i^T\theta + \mu \|\theta\|_1
9+
10+
// This is the "negative" of the problem as in https://gist.github.com/jonathan-taylor/07774d209173f8bc4e42aa37712339bf
11+
// Therefore we don't have to negate the answer to get theta.
12+
// Update one coordinate
13+
14+
double objective(double *Sigma, /* A covariance matrix: X^TX/n */
15+
int nrow, /* how many rows in Sigma */
16+
int row, /* which row: 0-based */
17+
double bound, /* Lagrange multipler for \ell_1 */
18+
double *theta) /* current value */
19+
{
20+
int irow, icol;
21+
double value = 0;
22+
double *Sigma_ptr = Sigma;
23+
double *theta_row_ptr, *theta_col_ptr;
24+
25+
theta_row_ptr = theta;
26+
theta_col_ptr = theta;
27+
28+
for (irow=0; irow<nrow; irow++) {
29+
double *theta_col_ptr = theta;
30+
if (*theta_row_ptr != 0) {
31+
for (icol=0; icol<nrow; icol++) {
32+
value += 0.5 * (*Sigma_ptr) * (*theta_row_ptr) * (*theta_col_ptr);
33+
Sigma_ptr++;
34+
theta_col_ptr++;
35+
}
36+
}
37+
if (irow == row) {
38+
value -= (*theta_row_ptr); // the elementary basis vector term
39+
}
40+
value = value + bound * fabs((*theta_row_ptr)); // the \ell_1 term
41+
theta_row_ptr++;
42+
}
43+
44+
return(value);
45+
}
46+
47+
48+
double update_one_coord(double *Sigma, /* A covariance matrix: X^TX/n */
49+
double *Sigma_diag, /* Diagonal entries of Sigma */
50+
double *Sigma_theta, /* Sigma times theta */
51+
int nrow, /* How many rows in Sigma */
52+
double bound, /* feasibility parameter */
53+
double *theta, /* current value */
54+
int row, /* which row: 0-based */
55+
int coord) /* which coordinate to update: 0-based */
56+
{
57+
58+
double delta;
59+
double linear_term = 0;
60+
double value = 0;
61+
double old_value;
62+
double *Sigma_ptr;
63+
double *Sigma_theta_ptr;
64+
double *theta_ptr;
65+
int icol = 0;
66+
67+
double *quadratic_ptr = ((double *) Sigma_diag + coord);
68+
double quadratic_term = *quadratic_ptr;
69+
70+
Sigma_theta_ptr = ((double *) Sigma_theta + coord);
71+
linear_term = *Sigma_theta_ptr;
72+
73+
theta_ptr = ((double *) theta + coord);
74+
old_value = *theta_ptr;
75+
76+
// The coord entry of Sigma_theta term has a diagonal term in it:
77+
// Sigma[coord, coord] * theta[coord]
78+
// This removes it.
79+
linear_term -= quadratic_term * old_value;
80+
81+
if (row == coord) {
82+
linear_term -= 1;
83+
}
84+
85+
// Now soft-threshold the coord entry of theta
86+
87+
// Objective is t \mapsto q/2 * t^2 + l * t + bound |t|
88+
// with q=quadratic_term and l=linear_term
89+
90+
// With a negative linear term, solution should be
91+
// positive
92+
93+
if (linear_term < -bound) {
94+
value = (-linear_term - bound) / quadratic_term;
95+
}
96+
else if (linear_term > bound) {
97+
value = -(linear_term - bound) / quadratic_term;
98+
}
99+
100+
if (fabs(old_value - value) > 1.e-6 * (fabs(value) + fabs(old_value))) { // Update the linear term
101+
102+
delta = value - old_value;
103+
Sigma_ptr = ((double *) Sigma + coord * nrow);
104+
Sigma_theta_ptr = ((double *) Sigma_theta);
105+
106+
for (icol=0; icol<nrow; icol++) {
107+
*Sigma_theta_ptr = *Sigma_theta_ptr + delta * (*Sigma_ptr);
108+
Sigma_theta_ptr += 1;
109+
Sigma_ptr += 1;
110+
}
111+
112+
theta_ptr = ((double *) theta + coord);
113+
*theta_ptr = value;
114+
115+
}
116+
117+
return(value);
118+
119+
}
120+
121+
void find_one_row(double *Sigma, /* A covariance matrix: X^TX/n */
122+
double *Sigma_diag, /* Diagonal entry of covariance matrix */
123+
double *Sigma_theta, /* Sigma times theta */
124+
int *nrow_ptr, /* How many rows in Sigma */
125+
double *bound_ptr, /* feasibility parameter */
126+
double *theta, /* current value */
127+
int *maxiter_ptr, /* how many iterations */
128+
int *row_ptr) /* which coordinate to update: 0-based */
129+
{
130+
131+
int maxiter = *maxiter_ptr;
132+
int iter = 0;
133+
int icoord = 0;
134+
int row = *row_ptr;
135+
double bound = *bound_ptr;
136+
int nrow = *nrow_ptr;
137+
138+
double old_value = objective(Sigma,
139+
nrow,
140+
row,
141+
bound,
142+
theta);
143+
double new_value;
144+
double tol=1.e-10;
145+
146+
for (iter=0; iter<maxiter; iter++) {
147+
148+
// Update the diagonal first
149+
150+
update_one_coord(Sigma,
151+
Sigma_diag,
152+
Sigma_theta,
153+
nrow,
154+
bound,
155+
theta,
156+
row,
157+
row);
158+
159+
for (icoord=0; icoord<nrow; icoord++) {
160+
161+
update_one_coord(Sigma,
162+
Sigma_diag,
163+
Sigma_theta,
164+
nrow,
165+
bound,
166+
theta,
167+
row,
168+
icoord);
169+
}
170+
171+
new_value = objective(Sigma,
172+
nrow,
173+
row,
174+
bound,
175+
theta);
176+
177+
if (((old_value - new_value) < tol * fabs(new_value)) && (iter > 0)) {
178+
break;
179+
}
180+
181+
old_value = new_value;
182+
}
183+
184+
*nrow_ptr = iter-1;
185+
}
186+

0 commit comments

Comments
 (0)