Skip to content

Commit 3405f23

Browse files
jonathan-taylorjonathan-taylor
committed
passing tol as argument; both solvers agree
1 parent 1237c78 commit 3405f23

File tree

5 files changed

+90
-53
lines changed

5 files changed

+90
-53
lines changed

‎selectiveInference/R/funs.fixed.R

Lines changed: 13 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -338,7 +338,7 @@ InverseLinfty <- function(sigma, n, e, resol=1.2, mu=NULL, maxiter=50, threshold
338338
return(M)
339339
}
340340

341-
InverseLinftyOneRow <- function (Sigma, i, mu, maxiter=50, soln_result=NULL) {
341+
InverseLinftyOneRow <- function (Sigma, i, mu, maxiter=50, soln_result=NULL, kkt_tol=1.e-6, objective_tol=1.e-6) {
342342

343343
# If soln_result is not Null, it is used as a warm start.
344344
# It should be a list
@@ -365,10 +365,18 @@ InverseLinftyOneRow <- function (Sigma, i, mu, maxiter=50, soln_result=NULL) {
365365
linear_func = soln_result$linear_func
366366
}
367367

368-
result = find_one_row_debiasingM(Sigma, i, mu, maxiter, 0 * soln, gradient, ever_active, nactive) # C function uses 0-based indexing
369-
#result1 = solve_QP(Sigma, mu, maxiter, soln, linear_func, gradient, ever_active, nactive)
370-
#print("close?")
371-
#print(c(sqrt(sum((result1$soln-result$soln)^2)/sum(result$soln^2)), sum(result$soln^2)))
368+
soln1 = rep(0, p)
369+
gradient1 = rep(0, p)
370+
ever_active1 = rep(0, p)
371+
ever_active1[1] = i-1
372+
nactive1 = as.integer(1)
373+
result1 = find_one_row_debiasingM(Sigma, i-1, mu, maxiter, soln1, gradient1, ever_active1, nactive1, kkt_tol, objective_tol) # C function uses 0-based indexing
374+
result = solve_QP(Sigma, mu, maxiter, soln, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol)
375+
print("close?")
376+
print(c(sqrt(sum((result1$soln-result$soln)^2)/sum(result$soln^2)), sum(result$soln^2)))
377+
print(c(result1$iter, result$iter, sum(result1$soln^2)))
378+
379+
#result = find_one_row_debiasingM(Sigma, i-1, mu, maxiter, soln, gradient, ever_active, nactive, kkt_tol, objective_tol) # C function uses 0-based indexing
372380

373381
# Check feasibility
374382

‎selectiveInference/src/Rcpp-debias.cpp

Lines changed: 16 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,9 @@ Rcpp::List solve_QP(Rcpp::NumericMatrix Sigma,
99
Rcpp::NumericVector linear_func,
1010
Rcpp::NumericVector gradient,
1111
Rcpp::IntegerVector ever_active,
12-
Rcpp::IntegerVector nactive
12+
Rcpp::IntegerVector nactive,
13+
double kkt_tol,
14+
double objective_tol
1315
) {
1416

1517
int nrow = Sigma.nrow(); // number of features
@@ -37,14 +39,17 @@ Rcpp::List solve_QP(Rcpp::NumericMatrix Sigma,
3739
nrow,
3840
bound,
3941
(double *) theta.begin(),
40-
maxiter);
42+
maxiter,
43+
kkt_tol,
44+
objective_tol);
4145

4246
// Check whether feasible
4347

4448
int kkt_check = check_KKT_qp(theta.begin(),
4549
gradient.begin(),
4650
nrow,
47-
bound);
51+
bound,
52+
kkt_tol);
4853

4954
return(Rcpp::List::create(Rcpp::Named("soln") = theta,
5055
Rcpp::Named("gradient") = gradient,
@@ -64,7 +69,9 @@ Rcpp::List find_one_row_debiasingM(Rcpp::NumericMatrix Sigma,
6469
Rcpp::NumericVector theta,
6570
Rcpp::NumericVector gradient,
6671
Rcpp::IntegerVector ever_active,
67-
Rcpp::IntegerVector nactive
72+
Rcpp::IntegerVector nactive,
73+
double kkt_tol,
74+
double objective_tol
6875
) {
6976

7077
int nrow = Sigma.nrow(); // number of features
@@ -92,15 +99,18 @@ Rcpp::List find_one_row_debiasingM(Rcpp::NumericMatrix Sigma,
9299
bound,
93100
(double *) theta.begin(),
94101
maxiter,
95-
row);
102+
row,
103+
kkt_tol,
104+
objective_tol);
96105

97106
// Check whether feasible
98107

99108
int kkt_check = check_KKT(theta.begin(),
100109
gradient.begin(),
101110
nrow,
102111
row,
103-
bound);
112+
bound,
113+
kkt_tol);
104114

105115
return(Rcpp::List::create(Rcpp::Named("soln") = theta,
106116
Rcpp::Named("gradient") = gradient,

‎selectiveInference/src/debias.c

Lines changed: 37 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,4 @@
1+
#include <stdio.h>
12
#include <math.h> // for fabs
23

34
// Find an approximate row of \hat{Sigma}^{-1}
@@ -67,6 +68,7 @@ int update_ever_active(int coord,
6768
for (iactive=0; iactive<nactive; iactive++) {
6869
ever_active_ptr_tmp = ((int *) ever_active_ptr + iactive);
6970
active_var = (*ever_active_ptr_tmp);
71+
// fprintf(stderr, "%d %d\n", iactive, active_var);
7072
if (active_var == coord) {
7173
return(1);
7274
}
@@ -85,16 +87,16 @@ int update_ever_active(int coord,
8587
return(0);
8688
}
8789

88-
int check_KKT(double *theta, /* current theta */
90+
int check_KKT(double *theta, /* current theta */
8991
double *gradient_ptr, /* Sigma times theta */
90-
int nrow, /* how many rows in Sigma */
91-
int row, /* which row: 0-based */
92-
double bound) /* Lagrange multipler for \ell_1 */
92+
int nrow, /* how many rows in Sigma */
93+
int row, /* which row: 0-based */
94+
double bound, /* Lagrange multipler for \ell_1 */
95+
double tol) /* precision for checking KKT conditions */
9396
{
9497
// First check inactive
9598

9699
int irow;
97-
double tol = 1.e-6;
98100
double *theta_ptr, *gradient_ptr_tmp;
99101
double gradient;
100102

@@ -105,9 +107,12 @@ int check_KKT(double *theta, /* current theta */
105107
// Compute this coordinate of the gradient
106108

107109
gradient = *gradient_ptr_tmp;
108-
// if (row == irow) {
109-
// gradient -= 1;
110-
// }
110+
111+
// For the basis vector
112+
113+
if (row == irow) {
114+
gradient -= 1;
115+
}
111116

112117
if (*theta_ptr != 0) { // these coordinates of gradients should be equal to -bound
113118
if ((*theta_ptr > 0) && (fabs(gradient + bound) > tol * bound)) {
@@ -122,22 +127,23 @@ int check_KKT(double *theta, /* current theta */
122127
return(0);
123128
}
124129
}
130+
125131
}
126132

127133
return(1);
128134
}
129135

130136
double update_one_coord(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
131137
double *Sigma_diag_ptr, /* Diagonal entries of Sigma */
132-
double *gradient_ptr, /* Sigma times theta */
138+
double *gradient_ptr, /* Sigma times theta */
133139
int *ever_active_ptr, /* Ever active set: 0-based */
134-
int *nactive_ptr, /* Size of ever active set */
135-
int nrow, /* How many rows in Sigma */
136-
double bound, /* feasibility parameter */
137-
double *theta, /* current value */
138-
int row, /* which row: 0-based */
139-
int coord, /* which coordinate to update: 0-based */
140-
int is_active) /* Is this part of ever_active */
140+
int *nactive_ptr, /* Size of ever active set */
141+
int nrow, /* How many rows in Sigma */
142+
double bound, /* feasibility parameter */
143+
double *theta, /* current value */
144+
int row, /* which row: 0-based */
145+
int coord, /* which coordinate to update: 0-based */
146+
int is_active) /* Is this part of ever_active */
141147
{
142148

143149
double delta;
@@ -215,21 +221,23 @@ double update_one_coord(double *Sigma_ptr, /* A covariance matrix: X^T
215221

216222
int find_one_row_(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
217223
double *Sigma_diag_ptr, /* Diagonal entry of covariance matrix */
218-
double *gradient_ptr, /* Sigma times theta */
224+
double *gradient_ptr, /* Sigma times theta */
219225
int *ever_active_ptr, /* Ever active set: 0-based */
220-
int *nactive_ptr, /* Size of ever active set */
221-
int nrow, /* How many rows in Sigma */
222-
double bound, /* feasibility parameter */
223-
double *theta, /* current value */
224-
int maxiter, /* how many iterations */
225-
int row) /* which coordinate to update: 0-based */
226+
int *nactive_ptr, /* Size of ever active set */
227+
int nrow, /* How many rows in Sigma */
228+
double bound, /* feasibility parameter */
229+
double *theta, /* current value */
230+
int maxiter, /* how many iterations */
231+
int row, /* which coordinate to update: 0-based */
232+
double kkt_tol, /* precision for checking KKT conditions */
233+
double objective_tol) /* precision for checking relative decrease in objective value */
226234
{
227235

228236
int iter = 0;
229237
int icoord = 0;
230238
int iactive = 0;
231239
int *active_ptr;
232-
double old_value, new_value, tol=1.e-5;
240+
double old_value, new_value;
233241

234242
int check_objective = 1;
235243

@@ -272,7 +280,8 @@ int find_one_row_(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
272280
gradient_ptr,
273281
nrow,
274282
row,
275-
bound) == 1) {
283+
bound,
284+
kkt_tol) == 1) {
276285
break;
277286
}
278287

@@ -299,7 +308,8 @@ int find_one_row_(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
299308
gradient_ptr,
300309
nrow,
301310
row,
302-
bound) == 1) {
311+
bound,
312+
kkt_tol) == 1) {
303313
break;
304314
}
305315

@@ -312,7 +322,7 @@ int find_one_row_(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
312322
bound,
313323
theta);
314324

315-
if (((old_value - new_value) < tol * fabs(new_value)) && (iter > 0)) {
325+
if (((old_value - new_value) < objective_tol * fabs(new_value)) && (iter > 0)) {
316326
break;
317327
}
318328

‎selectiveInference/src/debias.h

Lines changed: 12 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -12,30 +12,36 @@ int solve_qp(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
1212
int nrow, /* How many rows in Sigma */
1313
double bound, /* feasibility parameter */
1414
double *theta, /* current value */
15-
int maxiter);
15+
int maxiter, /* how many iterations */
16+
double kkt_tol, /* precision for checking KKT conditions */
17+
double objective_tol); /* precision for checking relative decrease in objective value */
18+
1619

1720
int check_KKT_qp(double *theta, /* current theta */
1821
double *gradient_ptr, /* Current gradient of quadratic loss */
1922
int nrow, /* how many rows in Sigma */
20-
double bound); /* Lagrange multipler for \ell_1 */
23+
double bound, /* Lagrange multipler for \ell_1 */
24+
double tol); /* precision for checking KKT conditions */
2125

2226
int find_one_row_(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
23-
double *Sigma_ptr_ptr, /* Diagonal entry of covariance matrix */
27+
double *Sigma_diag_ptr, /* Diagonal entry of covariance matrix */
2428
double *gradient_ptr, /* Current gradient of quadratic loss */
2529
int *ever_active_ptr, /* Ever active set: 0-based */
2630
int *nactive_ptr, /* Size of ever active set */
2731
int nrow, /* How many rows in Sigma */
2832
double bound, /* feasibility parameter */
2933
double *theta, /* current value */
3034
int maxiter, /* how many iterations */
31-
int row); /* which coordinate to update: 0-based */
35+
int row, /* which coordinate to update: 0-based */
36+
double kkt_tol, /* precision for checking KKT conditions */
37+
double objective_tol); /* precision for checking relative decrease in objective value */
3238

3339
int check_KKT(double *theta, /* current theta */
3440
double *gradient_ptr, /* Current gradient of quadratic loss */
3541
int nrow, /* how many rows in Sigma */
3642
int row, /* which row: 0-based */
37-
double bound); /* Lagrange multipler for \ell_1 */
38-
43+
double bound, /* Lagrange multipler for \ell_1 */
44+
double kkt_tol); /* precision for checking KKT conditions */
3945

4046
#ifdef __cplusplus
4147
} /* extern "C" */

‎selectiveInference/src/quadratic_program.c

Lines changed: 12 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -88,15 +88,15 @@ int update_ever_active_qp(int coord,
8888
return(0);
8989
}
9090

91-
int check_KKT_qp(double *theta, /* current theta */
91+
int check_KKT_qp(double *theta, /* current theta */
9292
double *gradient_ptr, /* Sigma times theta */
93-
int nrow, /* how many rows in Sigma */
94-
double bound) /* Lagrange multipler for \ell_1 */
93+
int nrow, /* how many rows in Sigma */
94+
double bound, /* Lagrange multipler for \ell_1 */
95+
double tol) /* precision for checking KKT conditions */
9596
{
9697
// First check inactive
9798

9899
int irow;
99-
double tol = 1.e-6;
100100
double *theta_ptr, *gradient_ptr_tmp;
101101
double gradient;
102102

@@ -216,7 +216,9 @@ int solve_qp(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
216216
int nrow, /* How many rows in Sigma */
217217
double bound, /* feasibility parameter */
218218
double *theta, /* current value */
219-
int maxiter)
219+
int maxiter, /* max number of iterations */
220+
double kkt_tol, /* precision for checking KKT conditions */
221+
double objective_tol) /* precision for checking relative decrease in objective value */
220222
{
221223

222224
int iter = 0;
@@ -227,7 +229,6 @@ int solve_qp(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
227229
int check_objective = 1;
228230

229231
double old_value, new_value;
230-
double tol=1.e-8;
231232

232233
if (check_objective) {
233234

@@ -268,7 +269,8 @@ int solve_qp(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
268269
if (check_KKT_qp(theta,
269270
gradient_ptr,
270271
nrow,
271-
bound) == 1) {
272+
bound,
273+
kkt_tol) == 1) {
272274
break;
273275
}
274276

@@ -294,7 +296,8 @@ int solve_qp(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
294296
if (check_KKT_qp(theta,
295297
gradient_ptr,
296298
nrow,
297-
bound) == 1) {
299+
bound,
300+
kkt_tol) == 1) {
298301
break;
299302
}
300303

@@ -307,7 +310,7 @@ int solve_qp(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
307310
bound,
308311
theta);
309312

310-
if (((old_value - new_value) < tol * fabs(new_value)) && (iter > 0)) {
313+
if ((fabs(old_value - new_value) < objective_tol * fabs(new_value)) && (iter > 0)) {
311314
break;
312315
}
313316
old_value = new_value;

0 commit comments

Comments
 (0)