Skip to content

Commit 7fa5120

Browse files
jonathan-taylorjonathan-taylor
committed
NF: stop early if active set gets too big
1 parent 3383e31 commit 7fa5120

File tree

5 files changed

+61
-19
lines changed

5 files changed

+61
-19
lines changed

‎selectiveInference/R/funs.fixed.R

Lines changed: 21 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -271,7 +271,9 @@ fixedLasso.poly=
271271
## Approximates inverse covariance matrix theta
272272
InverseLinfty <- function(sigma, n, e, resol=1.2, mu=NULL, maxiter=50, threshold=1e-2, verbose = TRUE, max.try=10) {
273273

274-
isgiven <- 1;
274+
max_active = max(50, 2 * sqrt(n))
275+
276+
isgiven <- 1;
275277
if (is.null(mu)){
276278
isgiven <- 0;
277279
}
@@ -295,17 +297,19 @@ InverseLinfty <- function(sigma, n, e, resol=1.2, mu=NULL, maxiter=50, threshold
295297

296298
output = NULL
297299

300+
last.beta = NULL
301+
298302
while ((mu.stop != 1) && (try.no<max.try) ){
299-
last.beta <- beta
300303

301-
output <- InverseLinftyOneRow(sigma, i, mu, maxiter=maxiter, soln_result=output) # uses a warm start
304+
output <- InverseLinftyOneRow(sigma, i, mu, maxiter=maxiter, soln_result=output, max_active=max_active) # uses a warm start
302305
beta <- output$soln
303306

304307
iter <- output$iter
308+
305309
if (isgiven==1) {
306310
mu.stop <- 1
307311
}
308-
else{
312+
else {
309313
if (try.no==1){
310314
if (iter == (maxiter+1)){
311315
incr <- 1;
@@ -325,14 +329,19 @@ InverseLinfty <- function(sigma, n, e, resol=1.2, mu=NULL, maxiter=50, threshold
325329
if ((incr == 0)&&(iter < (maxiter+1))){
326330
mu <- mu/resol;
327331
}
328-
if ((incr == 0)&&(iter == (maxiter+1))){
332+
if ((incr == 0) && (iter == (maxiter+1))){
329333
mu <- mu*resol;
330334
beta <- last.beta;
331335
mu.stop <- 1;
332336
}
333337
}
338+
if (output$max_active_check) {
339+
mu.stop <- 1;
340+
beta <- last.beta;
341+
}
334342
}
335343
try.no <- try.no+1
344+
last.beta <- beta
336345
}
337346
M[i,] <- beta;
338347
}
@@ -341,14 +350,18 @@ InverseLinfty <- function(sigma, n, e, resol=1.2, mu=NULL, maxiter=50, threshold
341350

342351

343352
InverseLinftyOneRow <- function (Sigma, i, mu, maxiter=50, soln_result=NULL, kkt_tol=1.e-6, objective_tol=1.e-6,
344-
use_QP=TRUE) {
353+
use_QP=TRUE, max_active=NULL) {
345354

346355
# If soln_result is not Null, it is used as a warm start.
347356
# It should be a list
348357
# with entries "soln", "gradient", "ever_active", "nactive"
349358

350359
p = nrow(Sigma)
351360

361+
if (is.null(max_active)) {
362+
max_active = 50 # arbitrary?
363+
}
364+
352365
if (is.null(soln_result)) {
353366
soln = rep(0, p)
354367
ever_active = rep(0, p)
@@ -375,9 +388,9 @@ InverseLinftyOneRow <- function (Sigma, i, mu, maxiter=50, soln_result=NULL, kkt
375388
}
376389

377390
if (use_QP) {
378-
result = solve_QP(Sigma, mu, maxiter, soln, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol)
391+
result = solve_QP(Sigma, mu, maxiter, soln, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol, max_active)
379392
} else {
380-
result = find_one_row_debiasingM(Sigma, i, mu, maxiter, soln, gradient, ever_active, nactive, kkt_tol, objective_tol) # C function uses 1-based indexing for active set
393+
result = find_one_row_debiasingM(Sigma, i, mu, maxiter, soln, gradient, ever_active, nactive, kkt_tol, objective_tol, max_active) # C function uses 1-based indexing for active set
381394
}
382395

383396
# Check feasibility

‎selectiveInference/src/Rcpp-debias.cpp

Lines changed: 16 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,8 @@ Rcpp::List solve_QP(Rcpp::NumericMatrix Sigma,
1414
Rcpp::IntegerVector ever_active,
1515
Rcpp::IntegerVector nactive,
1616
double kkt_tol,
17-
double objective_tol
17+
double objective_tol,
18+
int max_active
1819
) {
1920

2021
int nrow = Sigma.nrow(); // number of features
@@ -44,7 +45,8 @@ Rcpp::List solve_QP(Rcpp::NumericMatrix Sigma,
4445
(double *) theta.begin(),
4546
maxiter,
4647
kkt_tol,
47-
objective_tol);
48+
objective_tol,
49+
max_active);
4850

4951
// Check whether feasible
5052

@@ -54,13 +56,16 @@ Rcpp::List solve_QP(Rcpp::NumericMatrix Sigma,
5456
bound,
5557
kkt_tol);
5658

59+
int max_active_check = (*(nactive.begin()) >= max_active);
60+
5761
return(Rcpp::List::create(Rcpp::Named("soln") = theta,
5862
Rcpp::Named("gradient") = gradient,
5963
Rcpp::Named("linear_func") = linear_func,
6064
Rcpp::Named("iter") = iter,
6165
Rcpp::Named("kkt_check") = kkt_check,
6266
Rcpp::Named("ever_active") = ever_active,
63-
Rcpp::Named("nactive") = nactive));
67+
Rcpp::Named("nactive") = nactive,
68+
Rcpp::Named("max_active_check") = max_active_check));
6469

6570
}
6671

@@ -74,7 +79,8 @@ Rcpp::List find_one_row_debiasingM(Rcpp::NumericMatrix Sigma,
7479
Rcpp::IntegerVector ever_active,
7580
Rcpp::IntegerVector nactive,
7681
double kkt_tol,
77-
double objective_tol
82+
double objective_tol,
83+
int max_active
7884
) {
7985

8086
int nrow = Sigma.nrow(); // number of features
@@ -104,7 +110,8 @@ Rcpp::List find_one_row_debiasingM(Rcpp::NumericMatrix Sigma,
104110
maxiter,
105111
row,
106112
kkt_tol,
107-
objective_tol);
113+
objective_tol,
114+
max_active);
108115

109116
// Check whether feasible
110117

@@ -115,11 +122,14 @@ Rcpp::List find_one_row_debiasingM(Rcpp::NumericMatrix Sigma,
115122
bound,
116123
kkt_tol);
117124

125+
int max_active_check = (*(nactive.begin()) >= max_active);
126+
118127
return(Rcpp::List::create(Rcpp::Named("soln") = theta,
119128
Rcpp::Named("gradient") = gradient,
120129
Rcpp::Named("iter") = iter,
121130
Rcpp::Named("kkt_check") = kkt_check,
122131
Rcpp::Named("ever_active") = ever_active,
123-
Rcpp::Named("nactive") = nactive));
132+
Rcpp::Named("nactive") = nactive,
133+
Rcpp::Named("max_active_check") = max_active_check));
124134

125135
}

‎selectiveInference/src/debias.c

Lines changed: 10 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -228,7 +228,8 @@ int find_one_row_(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
228228
int maxiter, /* how many iterations */
229229
int row, /* which coordinate to solve: 1-based */
230230
double kkt_tol, /* precision for checking KKT conditions */
231-
double objective_tol) /* precision for checking relative decrease in objective value */
231+
double objective_tol, /* precision for checking relative decrease in objective value */
232+
int max_active) /* Upper limit for size of active set -- otherwise break */
232233
{
233234

234235
int iter = 0;
@@ -311,6 +312,14 @@ int find_one_row_(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
311312
break;
312313
}
313314

315+
// Check size of active set
316+
317+
if (*nactive_ptr >= max_active) {
318+
break;
319+
}
320+
321+
// Check relative decrease of objective
322+
314323
if (check_objective) {
315324
new_value = objective(Sigma_ptr,
316325
ever_active_ptr,

‎selectiveInference/src/debias.h

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -14,8 +14,8 @@ int solve_qp(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
1414
double *theta, /* current value */
1515
int maxiter, /* how many iterations */
1616
double kkt_tol, /* precision for checking KKT conditions */
17-
double objective_tol); /* precision for checking relative decrease in objective value */
18-
17+
double objective_tol, /* precision for checking relative decrease in objective value */
18+
int max_active); /* Upper limit for size of active set -- otherwise break */
1919

2020
int check_KKT_qp(double *theta, /* current theta */
2121
double *gradient_ptr, /* Current gradient of quadratic loss */
@@ -34,7 +34,8 @@ int find_one_row_(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
3434
int maxiter, /* how many iterations */
3535
int row, /* which coordinate to update: 0-based */
3636
double kkt_tol, /* precision for checking KKT conditions */
37-
double objective_tol); /* precision for checking relative decrease in objective value */
37+
double objective_tol, /* precision for checking relative decrease in objective value */
38+
int max_active); /* Upper limit for size of active set -- otherwise break */
3839

3940
int check_KKT(double *theta, /* current theta */
4041
double *gradient_ptr, /* Current gradient of quadratic loss */

‎selectiveInference/src/quadratic_program.c

Lines changed: 10 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -220,7 +220,8 @@ int solve_qp(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
220220
double *theta, /* current value */
221221
int maxiter, /* max number of iterations */
222222
double kkt_tol, /* precision for checking KKT conditions */
223-
double objective_tol) /* precision for checking relative decrease in objective value */
223+
double objective_tol, /* precision for checking relative decrease in objective value */
224+
int max_active) /* Upper limit for size of active set -- otherwise break */
224225
{
225226

226227
int iter = 0;
@@ -303,6 +304,14 @@ int solve_qp(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
303304
break;
304305
}
305306

307+
// Check size of active set
308+
309+
if (*nactive_ptr >= max_active) {
310+
break;
311+
}
312+
313+
// Check relative decrease of objective
314+
306315
if (check_objective) {
307316
new_value = objective_qp(Sigma_ptr,
308317
linear_func_ptr,

0 commit comments

Comments
 (0)