Skip to content

Commit d65e437

Browse files
jonathan-taylorjonathan-taylor
committed
BF: everything was thrown into active set
2 parents bfad505 + 190c0a0 commit d65e437

File tree

4 files changed

+97
-69
lines changed

4 files changed

+97
-69
lines changed

‎Makefile

Lines changed: 11 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,14 @@
11
Rcpp:
22
- rm -f selectiveInference/src/RcppExports.cpp
33
- rm -f selectiveInference/R/RcppExports.R
4-
Rscript -e "library(Rcpp); Rcpp::compileAttributes('selectiveInference')"
4+
Rscript -e "library(Rcpp); Rcpp::compileAttributes('selectiveInference')"
5+
6+
install: Rcpp
7+
R CMD install selectiveInference
8+
9+
build:
10+
R CMD build selectiveInference
11+
12+
check: Rcpp build
13+
R CMD build selectiveInference
14+
R CMD check selectiveInference_1.2.2.tar.gz # fix this to be a script variable

‎selectiveInference/R/funs.fixed.R

Lines changed: 10 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -297,7 +297,7 @@ InverseLinfty <- function(sigma, n, e, resol=1.2, mu=NULL, maxiter=50, threshold
297297

298298
while ((mu.stop != 1)&&(try.no<10)){
299299
last.beta <- beta
300-
#print(c("trying ", try.no))
300+
print(c("#######################trying ", try.no))
301301
output <- InverseLinftyOneRow(sigma, i, mu, maxiter=maxiter, soln_result=output) # uses a warm start
302302
beta <- output$soln
303303
iter <- output$iter
@@ -344,9 +344,11 @@ InverseLinftyOneRow <- function (Sigma, i, mu, maxiter=50, soln_result=NULL) {
344344
# It should be a list
345345
# with entries "soln", "gradient", "ever_active", "nactive"
346346

347+
p = nrow(Sigma)
348+
347349
if (is.null(soln_result)) {
348-
soln = rep(0, nrow(Sigma))
349-
ever_active = rep(0, nrow(Sigma))
350+
soln = rep(0, p)
351+
ever_active = rep(0, p)
350352
ever_active[1] = i-1 # 0-based
351353
ever_active = as.integer(ever_active)
352354
nactive = as.integer(1)
@@ -363,7 +365,11 @@ InverseLinftyOneRow <- function (Sigma, i, mu, maxiter=50, soln_result=NULL) {
363365
linear_func = soln_result$linear_func
364366
}
365367

366-
result = find_one_row_debiasingM(Sigma, mu, maxiter, soln, linear_func, gradient, ever_active, nactive) # C function uses 0-based indexing
368+
result = solve_QP(Sigma, mu, maxiter, soln, linear_func, gradient, ever_active, nactive) # C function uses 0-based indexing
369+
result2 = find_one_row_debiasingM(Sigma, i, mu, maxiter, soln, gradient, ever_active, nactive) # C function uses 0-based indexing
370+
371+
print('close?')
372+
print(c(sqrt(sum((result$soln-result2$soln)^2)/sum(result$soln^2)), sqrt(sum(result$soln^2)), result2$nactive))
367373

368374
# Check feasibility
369375

‎selectiveInference/src/debias.c

Lines changed: 69 additions & 63 deletions
Original file line numberDiff line numberDiff line change
@@ -4,25 +4,23 @@
44

55
// Solves a dual version of problem (4) of https://arxiv.org/pdf/1306.3171.pdf
66

7-
// Dual problem: \text{min}_{\theta} 1/2 \theta^T \Sigma \theta - l^T\theta + \mu \|\theta\|_1
8-
// where l is `linear_func` below
7+
// Dual problem: \text{min}_{\theta} 1/2 \theta^T \Sigma \theta - e_i^T\theta + \mu \|\theta\|_1
98

109
// This is the "negative" of the problem as in https://gist.github.com/jonathan-taylor/07774d209173f8bc4e42aa37712339bf
1110
// Therefore we don't have to negate the answer to get theta.
1211
// Update one coordinate
1312

1413
double objective(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
15-
double *linear_func_ptr, /* Linear term in objective */
1614
int *ever_active_ptr, /* Ever active set: 0-based */
17-
int *nactive_ptr, /* Size of ever active set */
18-
int nrow, /* how many rows in Sigma */
19-
double bound, /* Lagrange multipler for \ell_1 */
20-
double *theta) /* current value */
15+
int *nactive_ptr, /* Size of ever active set */
16+
int nrow, /* how many rows in Sigma */
17+
int row, /* which row: 0-based */
18+
double bound, /* Lagrange multipler for \ell_1 */
19+
double *theta) /* current value */
2120
{
2221
int irow, icol;
2322
double value = 0;
2423
double *Sigma_ptr_tmp = Sigma_ptr;
25-
double *linear_func_ptr_tmp = linear_func_ptr;
2624
double *theta_row_ptr, *theta_col_ptr;
2725
int *active_row_ptr, *active_col_ptr;
2826
int active_row, active_col;
@@ -47,15 +45,12 @@ double objective(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
4745

4846
value += 0.5 * (*Sigma_ptr_tmp) * (*theta_row_ptr) * (*theta_col_ptr);
4947
}
50-
value += bound * fabs((*theta_row_ptr)); // the \ell_1 term
51-
52-
// The linear term in the objective
48+
value = value + bound * fabs((*theta_row_ptr)); // the \ell_1 term
49+
}
5350

54-
linear_func_ptr_tmp = ((double *) linear_func_ptr + active_row);
55-
value += (*linear_func_ptr_tmp) * (*theta_row_ptr);
51+
theta_row_ptr = ((double *) theta + row);
52+
value -= (*theta_row_ptr); // the elementary basis vector term
5653

57-
}
58-
5954
return(value);
6055
}
6156

@@ -71,14 +66,14 @@ int update_ever_active(int coord,
7166

7267
for (iactive=0; iactive<nactive; iactive++) {
7368
ever_active_ptr_tmp = ((int *) ever_active_ptr + iactive);
74-
active_var = *ever_active_ptr_tmp;
69+
active_var = (*ever_active_ptr_tmp);
7570
if (active_var == coord) {
7671
return(1);
7772
}
7873
}
79-
80-
// If we haven't returned yet, this means the coord was not in
81-
// ever_active.
74+
75+
// If we have not returned yet, this variable
76+
// was not in ever_active
8277

8378
// Add it to the active set and increment the
8479
// number of active variables
@@ -93,6 +88,7 @@ int update_ever_active(int coord,
9388
int check_KKT(double *theta, /* current theta */
9489
double *gradient_ptr, /* Sigma times theta */
9590
int nrow, /* how many rows in Sigma */
91+
int row, /* which row: 0-based */
9692
double bound) /* Lagrange multipler for \ell_1 */
9793
{
9894
// First check inactive
@@ -110,36 +106,39 @@ int check_KKT(double *theta, /* current theta */
110106
// Compute this coordinate of the gradient
111107

112108
gradient = *gradient_ptr_tmp;
109+
if (row == irow) {
110+
gradient -= 1;
111+
}
113112

114113
if (*theta_ptr != 0) { // these coordinates of gradients should be equal to -bound
115114
if ((*theta_ptr > 0) && (fabs(gradient + bound) > tol * bound)) {
116-
return(0);
115+
fail += 1;
117116
}
118117
else if ((*theta_ptr < 0) && (fabs(gradient - bound) > tol * bound)) {
119-
return(0);
118+
fail += 1;
120119
}
121120
}
122121
else {
123122
if (fabs(gradient) > (1. + tol) * bound) {
124-
return(0);
123+
fail += 1;
125124
}
126125
}
127126
}
128127

129-
return(1);
128+
return(fail == 0);
130129
}
131130

132131
double update_one_coord(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
133-
double *linear_func_ptr, /* Linear term in objective */
134132
double *Sigma_diag_ptr, /* Diagonal entries of Sigma */
135-
double *gradient_ptr, /* Sigma times theta */
133+
double *gradient_ptr, /* Sigma times theta */
136134
int *ever_active_ptr, /* Ever active set: 0-based */
137-
int *nactive_ptr, /* Size of ever active set */
138-
int nrow, /* How many rows in Sigma */
139-
double bound, /* feasibility parameter */
140-
double *theta, /* current value */
141-
int coord, /* which coordinate to update: 0-based */
142-
int is_active) /* Is this part of ever_active */
135+
int *nactive_ptr, /* Size of ever active set */
136+
int nrow, /* How many rows in Sigma */
137+
double bound, /* feasibility parameter */
138+
double *theta, /* current value */
139+
int row, /* which row: 0-based */
140+
int coord, /* which coordinate to update: 0-based */
141+
int is_active) /* Is this part of ever_active */
143142
{
144143

145144
double delta;
@@ -154,6 +153,8 @@ double update_one_coord(double *Sigma_ptr, /* A covariance matrix: X^T
154153
double *quadratic_ptr = ((double *) Sigma_diag_ptr + coord);
155154
double quadratic_term = *quadratic_ptr;
156155

156+
// int *ever_active_ptr_tmp;
157+
157158
gradient_ptr_tmp = ((double *) gradient_ptr + coord);
158159
linear_term = *gradient_ptr_tmp;
159160

@@ -163,9 +164,12 @@ double update_one_coord(double *Sigma_ptr, /* A covariance matrix: X^T
163164
// The coord entry of gradient_ptr term has a diagonal term in it:
164165
// Sigma[coord, coord] * theta[coord]
165166
// This removes it.
166-
167167
linear_term -= quadratic_term * old_value;
168168

169+
if (row == coord) {
170+
linear_term -= 1;
171+
}
172+
169173
// Now soft-threshold the coord entry of theta
170174

171175
// Objective is t \mapsto q/2 * t^2 + l * t + bound |t|
@@ -183,7 +187,7 @@ double update_one_coord(double *Sigma_ptr, /* A covariance matrix: X^T
183187

184188
// Add to active set if necessary
185189

186-
if (is_active == 0) {
190+
if ((is_active == 0) && (value != 0)) {
187191
update_ever_active(coord, ever_active_ptr, nactive_ptr);
188192
}
189193

@@ -211,31 +215,31 @@ double update_one_coord(double *Sigma_ptr, /* A covariance matrix: X^T
211215
}
212216

213217
int find_one_row_(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
214-
double *linear_func_ptr, /* Linear term in objective */
215218
double *Sigma_diag_ptr, /* Diagonal entry of covariance matrix */
216-
double *gradient_ptr, /* Sigma times theta */
219+
double *gradient_ptr, /* Sigma times theta */
217220
int *ever_active_ptr, /* Ever active set: 0-based */
218-
int *nactive_ptr, /* Size of ever active set */
219-
int nrow, /* How many rows in Sigma */
220-
double bound, /* feasibility parameter */
221-
double *theta, /* current value */
222-
int maxiter)
221+
int *nactive_ptr, /* Size of ever active set */
222+
int nrow, /* How many rows in Sigma */
223+
double bound, /* feasibility parameter */
224+
double *theta, /* current value */
225+
int maxiter, /* how many iterations */
226+
int row) /* which coordinate to update: 0-based */
223227
{
224228

225229
int iter = 0;
226230
int icoord = 0;
227231
int iactive = 0;
228232
int *active_ptr;
229233

230-
/* double old_value = objective(Sigma_ptr, */
231-
/* linear_func_ptr, */
232-
/* ever_active_ptr, */
233-
/* nactive_ptr, */
234-
/* nrow, */
235-
/* bound, */
236-
/* theta); */
234+
double old_value = objective(Sigma_ptr,
235+
ever_active_ptr,
236+
nactive_ptr,
237+
nrow,
238+
row,
239+
bound,
240+
theta);
237241
double new_value;
238-
double tol=1.e-8;
242+
double tol=1.e-5;
239243

240244
for (iter=0; iter<maxiter; iter++) {
241245

@@ -245,14 +249,14 @@ int find_one_row_(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
245249

246250
for (iactive=0; iactive < *nactive_ptr; iactive++) {
247251
update_one_coord(Sigma_ptr,
248-
linear_func_ptr,
249252
Sigma_diag_ptr,
250253
gradient_ptr,
251254
ever_active_ptr,
252255
nactive_ptr,
253256
nrow,
254257
bound,
255258
theta,
259+
row,
256260
*active_ptr,
257261
1);
258262
active_ptr++;
@@ -263,6 +267,7 @@ int find_one_row_(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
263267
if (check_KKT(theta,
264268
gradient_ptr,
265269
nrow,
270+
row,
266271
bound) == 1) {
267272
break;
268273
}
@@ -272,14 +277,14 @@ int find_one_row_(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
272277
for (icoord=0; icoord<nrow; icoord++) {
273278

274279
update_one_coord(Sigma_ptr,
275-
linear_func_ptr,
276280
Sigma_diag_ptr,
277281
gradient_ptr,
278282
ever_active_ptr,
279283
nactive_ptr,
280284
nrow,
281285
bound,
282286
theta,
287+
row,
283288
icoord,
284289
0);
285290
}
@@ -289,23 +294,24 @@ int find_one_row_(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
289294
if (check_KKT(theta,
290295
gradient_ptr,
291296
nrow,
297+
row,
292298
bound) == 1) {
293299
break;
294300
}
295301

296-
/* new_value = objective(Sigma_ptr, */
297-
/* linear_func_ptr, */
298-
/* ever_active_ptr, */
299-
/* nactive_ptr, */
300-
/* nrow, */
301-
/* bound, */
302-
/* theta); */
303-
304-
/* if (((old_value - new_value) < tol * fabs(new_value)) && (iter > 0)) { */
305-
/* break; */
306-
/* } */
307-
308-
// old_value = new_value;
302+
new_value = objective(Sigma_ptr,
303+
ever_active_ptr,
304+
nactive_ptr,
305+
nrow,
306+
row,
307+
bound,
308+
theta);
309+
310+
if (((old_value - new_value) < tol * fabs(new_value)) && (iter > 0)) {
311+
break;
312+
}
313+
314+
old_value = new_value;
309315
}
310316
return(iter);
311317
}

‎selectiveInference/src/quadratic_program.c

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
#include <math.h> // for fabs
2+
#include <stdio.h>
23

34
// Find an approximate row of \hat{Sigma}^{-1}
45

@@ -81,6 +82,8 @@ int update_ever_active_qp(int coord,
8182
// Add it to the active set and increment the
8283
// number of active variables
8384

85+
fprintf(stderr, "adding %d\n", coord);
86+
8487
ever_active_ptr_tmp = ((int *) ever_active_ptr + *nactive_ptr);
8588
*ever_active_ptr_tmp = coord;
8689
*nactive_ptr += 1;
@@ -181,7 +184,7 @@ double update_one_coord_qp(double *Sigma_ptr, /* A covariance matrix:
181184

182185
// Add to active set if necessary
183186

184-
if (is_active == 0) {
187+
if ((is_active == 0) && (value != 0)) {
185188
update_ever_active_qp(coord, ever_active_ptr, nactive_ptr);
186189
}
187190

@@ -230,6 +233,8 @@ int solve_qp(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
230233
double old_value, new_value;
231234
double tol=1.e-8;
232235

236+
fprintf(stderr, "%d nactive start\n", *nactive_ptr);
237+
233238
if (check_objective) {
234239

235240
old_value = objective_qp(Sigma_ptr,
@@ -245,6 +250,7 @@ int solve_qp(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
245250

246251
for (iter=0; iter<maxiter; iter++) {
247252

253+
fprintf(stderr, "%d nactive loop \n", *nactive_ptr);
248254
// Update the active variables first
249255

250256
active_ptr = (int *) ever_active_ptr;

0 commit comments

Comments
 (0)