Skip to content

Commit bf65665

Browse files
jonathan-taylorjonathan-taylor
committed
WIP: C code for wide version working -- need to work on R code
1 parent c26a06b commit bf65665

File tree

5 files changed

+92
-118
lines changed

5 files changed

+92
-118
lines changed

‎selectiveInference/src/Makevars

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@ PKG_CFLAGS= -I.
22
PKG_CPPFLAGS= -I.
33
PKG_LIBS=-L.
44

5-
$(SHLIB): Rcpp Rcpp-matrixcomps.o Rcpp-debias.o RcppExports.o quadratic_program.o
5+
$(SHLIB): Rcpp Rcpp-matrixcomps.o Rcpp-debias.o RcppExports.o quadratic_program.o quadratic_program_wide.o
66

77
clean:
88
rm -f *o

‎selectiveInference/src/Rcpp-debias.cpp

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
#include <Rcpp.h> // need to include the main Rcpp header file
2-
#include <debias.h> // where find_one_row_void is defined
2+
#include <debias.h> // where solve_QP, solve_QP_wide are defined
33

44
// Below, the gradient should be equal to Sigma * theta + linear_func!!
55
// No check is done on this.
@@ -97,23 +97,23 @@ Rcpp::List solve_QP_wide(Rcpp::NumericMatrix X,
9797
Rcpp::IntegerVector need_update(nfeature);
9898

9999
// Extract the diagonal
100-
Rcpp::NumericVector X_diag(nfeature);
101-
double *X_diag_p = X_diag.begin();
100+
Rcpp::NumericVector nndef_diag(nfeature);
101+
double *nndef_diag_p = nndef_diag.begin();
102102

103-
for (icase=0; icase<ncase; icase++) {
104-
X_diag_p[icase] = 0;
105-
for (ifeature=0; ifeature<nfeature; ifeature++) {
106-
X_diag_p[icase] += X(icase, ifeature) * X(icase, ifeature);
103+
for (ifeature=0; ifeature<nfeature; ifeature++) {
104+
nndef_diag_p[ifeature] = 0;
105+
for (icase=0; icase<ncase; icase++) {
106+
nndef_diag_p[ifeature] += X(icase, ifeature) * X(icase, ifeature);
107107
}
108-
X_diag_p[icase] = X_diag_p[icase] / ncase;
108+
nndef_diag_p[ifeature] = nndef_diag_p[ifeature] / ncase;
109109
}
110110

111111
// Now call our C function
112112

113113
int iter = solve_wide((double *) X.begin(),
114114
(double *) X_theta.begin(),
115115
(double *) linear_func.begin(),
116-
(double *) X_diag.begin(),
116+
(double *) nndef_diag.begin(),
117117
(double *) gradient.begin(),
118118
(int *) need_update.begin(),
119119
(int *) ever_active.begin(),

‎selectiveInference/src/debias.h

Lines changed: 23 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -3,31 +3,31 @@ extern "C"
33
{
44
#endif /* __cplusplus */
55

6-
int solve_qp(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
6+
int solve_qp(double *nndef_ptr, /* A non-negative definite matrix */
77
double *linear_func_ptr, /* Linear term in objective */
8-
double *Sigma_diag_ptr, /* Diagonal entry of covariance matrix */
9-
double *gradient_ptr, /* Current gradient of quadratic loss */
10-
int *ever_active_ptr, /* Ever active set: 0-based */
8+
double *nndef_diag_ptr, /* Diagonal of nndef */
9+
double *gradient_ptr, /* nndef times theta */
10+
int *ever_active_ptr, /* Ever active set: 1-based */
1111
int *nactive_ptr, /* Size of ever active set */
12-
int nrow, /* How many rows in Sigma */
12+
int nrow, /* How many rows in nndef */
1313
double bound, /* feasibility parameter */
1414
double *theta, /* current value */
15-
int maxiter, /* how many iterations */
15+
int maxiter, /* max number of iterations */
1616
double kkt_tol, /* precision for checking KKT conditions */
1717
double objective_tol, /* precision for checking relative decrease in objective value */
18-
int max_active); /* Upper limit for size of active set -- otherwise break */
18+
int max_active); /* Upper limit for size of active set -- otherwise break */
1919

20-
int check_KKT_qp(double *theta, /* current theta */
21-
double *gradient_ptr, /* Current gradient of quadratic loss */
22-
int nrow, /* how many rows in Sigma */
23-
double bound, /* Lagrange multipler for \ell_1 */
24-
double tol); /* precision for checking KKT conditions */
20+
int check_KKT_qp(double *theta, /* current theta */
21+
double *gradient_ptr, /* nndef times theta + linear_func */
22+
int nrow, /* how many rows in nndef */
23+
double bound, /* Lagrange multipler for \ell_1 */
24+
double tol) /* precision for checking KKT conditions */
2525

26-
int solve_wide(double *X_ptr, /* A design matrix */
26+
int solve_wide(double *X_ptr, /* Sqrt of non-neg def matrix -- X^TX/ncase = nndef */
2727
double *X_theta_ptr, /* Fitted values */
2828
double *linear_func_ptr, /* Linear term in objective */
29-
double *X_diag_ptr, /* Diagonal entry of covariance matrix */
30-
double *gradient_ptr, /* X times theta */
29+
double *nndef_diag_ptr, /* Diagonal entries of non-neg def matrix */
30+
double *gradient_ptr, /* X^TX/ncase times theta + linear_func*/
3131
int *need_update_ptr, /* Keeps track of updated gradient coords */
3232
int *ever_active_ptr, /* Ever active set: 1-based */
3333
int *nactive_ptr, /* Size of ever active set */
@@ -38,23 +38,23 @@ int solve_wide(double *X_ptr, /* A design matrix */
3838
int maxiter, /* max number of iterations */
3939
double kkt_tol, /* precision for checking KKT conditions */
4040
double objective_tol, /* precision for checking relative decrease in objective value */
41-
int max_active); /* Upper limit for size of active set -- otherwise break */
41+
int max_active); /* Upper limit for size of active set -- otherwise break */
4242

4343
int check_KKT_wide(double *theta_ptr, /* current theta */
44-
double *gradient_ptr, /* X^TX/n times theta */
44+
double *gradient_ptr, /* X^TX/ncase times theta + linear_func*/
4545
double *X_theta_ptr, /* Current fitted values */
46-
double *X_ptr, /* A design matrix */
47-
double *linear_func_ptr, /* Linear term in objective */
46+
double *X_ptr, /* Sqrt of non-neg def matrix -- X^TX/ncase = nndef */
47+
double *linear_func_ptr, /* Linear term in objective */
4848
int *need_update_ptr, /* Which coordinates need to be updated? */
4949
int nfeature, /* how many columns in X */
5050
int ncase, /* how many rows in X */
5151
double bound, /* Lagrange multipler for \ell_1 */
52-
double tol); /* precision for checking KKT conditions */
52+
double tol) /* precision for checking KKT conditions */
5353

54-
void update_gradient_wide(double *gradient_ptr, /* X^TX/n times theta */
54+
void update_gradient_wide(double *gradient_ptr, /* X^TX/ncase times theta + linear_func */
5555
double *X_theta_ptr, /* Current fitted values */
56-
double *X_ptr, /* A design matrix */
57-
double *linear_func_ptr, /* Linear term in objective */
56+
double *X_ptr, /* Sqrt of non-neg def matrix -- X^TX/ncase = nndef */
57+
double *linear_func_ptr, /* Linear term in objective */
5858
int *need_update_ptr, /* Which coordinates need to be updated? */
5959
int nfeature, /* how many columns in X */
6060
int ncase); /* how many rows in X */

‎selectiveInference/src/quadratic_program.c

Lines changed: 28 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
#include <math.h> // for fabs
22

3-
// Find an approximate row of \hat{Sigma}^{-1}
3+
// Find an approximate row of \hat{nndef}^{-1}
44

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

@@ -11,17 +11,17 @@
1111
// Therefore we don't have to negate the answer to get theta.
1212
// Update one coordinate
1313

14-
double objective_qp(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
14+
double objective_qp(double *nndef_ptr, /* A non-negative definite matrix */
1515
double *linear_func_ptr, /* Linear term in objective */
1616
int *ever_active_ptr, /* Ever active set: 0-based */
1717
int *nactive_ptr, /* Size of ever active set */
18-
int nrow, /* how many rows in Sigma */
18+
int nrow, /* how many rows in nndef */
1919
double bound, /* Lagrange multipler for \ell_1 */
2020
double *theta) /* current value */
2121
{
2222
int irow, icol;
2323
double value = 0;
24-
double *Sigma_ptr_tmp = Sigma_ptr;
24+
double *nndef_ptr_tmp = nndef_ptr;
2525
double *linear_func_ptr_tmp = linear_func_ptr;
2626
double *theta_row_ptr, *theta_col_ptr;
2727
int *active_row_ptr, *active_col_ptr;
@@ -43,9 +43,9 @@ double objective_qp(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
4343
active_col = *active_col_ptr - 1; // Ever-active is 1-based
4444
theta_col_ptr = ((double *) theta + active_col);
4545

46-
Sigma_ptr_tmp = ((double *) Sigma_ptr + nrow * active_col + active_row); // Matrices are column-major order
46+
nndef_ptr_tmp = ((double *) nndef_ptr + nrow * active_col + active_row); // Matrices are column-major order
4747

48-
value += 0.5 * (*Sigma_ptr_tmp) * (*theta_row_ptr) * (*theta_col_ptr);
48+
value += 0.5 * (*nndef_ptr_tmp) * (*theta_row_ptr) * (*theta_col_ptr);
4949
}
5050
value += bound * fabs((*theta_row_ptr)); // the \ell_1 term
5151

@@ -91,8 +91,8 @@ int update_ever_active_qp(int coord,
9191
}
9292

9393
int check_KKT_qp(double *theta, /* current theta */
94-
double *gradient_ptr, /* Sigma times theta */
95-
int nrow, /* how many rows in Sigma */
94+
double *gradient_ptr, /* nndef times theta + linear_func */
95+
int nrow, /* how many rows in nndef */
9696
double bound, /* Lagrange multipler for \ell_1 */
9797
double tol) /* precision for checking KKT conditions */
9898
{
@@ -128,13 +128,13 @@ int check_KKT_qp(double *theta, /* current theta */
128128
return(1);
129129
}
130130

131-
double update_one_coord_qp(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
131+
double update_one_coord_qp(double *nndef_ptr, /* A non-negative definite matrix */
132132
double *linear_func_ptr, /* Linear term in objective */
133-
double *Sigma_diag_ptr, /* Diagonal entries of Sigma */
134-
double *gradient_ptr, /* Sigma times theta */
133+
double *nndef_diag_ptr, /* Diagonal of nndef */
134+
double *gradient_ptr, /* nndef times theta + linear_func */
135135
int *ever_active_ptr, /* Ever active set: 1-based */
136136
int *nactive_ptr, /* Size of ever active set */
137-
int nrow, /* How many rows in Sigma */
137+
int nrow, /* How many rows in nndef */
138138
double bound, /* feasibility parameter */
139139
double *theta, /* current value */
140140
int coord, /* which coordinate to update: 0-based */
@@ -145,12 +145,12 @@ double update_one_coord_qp(double *Sigma_ptr, /* A covariance matrix:
145145
double linear_term = 0;
146146
double value = 0;
147147
double old_value;
148-
double *Sigma_ptr_tmp;
148+
double *nndef_ptr_tmp;
149149
double *gradient_ptr_tmp;
150150
double *theta_ptr;
151151
int icol = 0;
152152

153-
double *quadratic_ptr = ((double *) Sigma_diag_ptr + coord);
153+
double *quadratic_ptr = ((double *) nndef_diag_ptr + coord);
154154
double quadratic_term = *quadratic_ptr;
155155

156156
gradient_ptr_tmp = ((double *) gradient_ptr + coord);
@@ -160,7 +160,7 @@ double update_one_coord_qp(double *Sigma_ptr, /* A covariance matrix:
160160
old_value = *theta_ptr;
161161

162162
// The coord entry of gradient_ptr term has a diagonal term in it:
163-
// Sigma[coord, coord] * theta[coord]
163+
// nndef[coord, coord] * theta[coord]
164164
// This removes it.
165165

166166
linear_term -= quadratic_term * old_value;
@@ -191,13 +191,13 @@ double update_one_coord_qp(double *Sigma_ptr, /* A covariance matrix:
191191
if (fabs(old_value - value) > 1.e-6 * (fabs(value) + fabs(old_value))) {
192192

193193
delta = value - old_value;
194-
Sigma_ptr_tmp = ((double *) Sigma_ptr + coord * nrow);
194+
nndef_ptr_tmp = ((double *) nndef_ptr + coord * nrow);
195195
gradient_ptr_tmp = ((double *) gradient_ptr);
196196

197197
for (icol=0; icol<nrow; icol++) {
198-
*gradient_ptr_tmp = *gradient_ptr_tmp + delta * (*Sigma_ptr_tmp);
198+
*gradient_ptr_tmp = *gradient_ptr_tmp + delta * (*nndef_ptr_tmp);
199199
gradient_ptr_tmp += 1;
200-
Sigma_ptr_tmp += 1;
200+
nndef_ptr_tmp += 1;
201201
}
202202

203203
theta_ptr = ((double *) theta + coord);
@@ -209,13 +209,13 @@ double update_one_coord_qp(double *Sigma_ptr, /* A covariance matrix:
209209

210210
}
211211

212-
int solve_qp(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
212+
int solve_qp(double *nndef_ptr, /* A non-negative definite matrix */
213213
double *linear_func_ptr, /* Linear term in objective */
214-
double *Sigma_diag_ptr, /* Diagonal entry of covariance matrix */
215-
double *gradient_ptr, /* Sigma times theta */
214+
double *nndef_diag_ptr, /* Diagonal of nndef */
215+
double *gradient_ptr, /* nndef times theta */
216216
int *ever_active_ptr, /* Ever active set: 1-based */
217217
int *nactive_ptr, /* Size of ever active set */
218-
int nrow, /* How many rows in Sigma */
218+
int nrow, /* How many rows in nndef */
219219
double bound, /* feasibility parameter */
220220
double *theta, /* current value */
221221
int maxiter, /* max number of iterations */
@@ -235,7 +235,7 @@ int solve_qp(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
235235

236236
if (check_objective) {
237237

238-
old_value = objective_qp(Sigma_ptr,
238+
old_value = objective_qp(nndef_ptr,
239239
linear_func_ptr,
240240
ever_active_ptr,
241241
nactive_ptr,
@@ -253,9 +253,9 @@ int solve_qp(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
253253
active_ptr = (int *) ever_active_ptr;
254254

255255
for (iactive=0; iactive < *nactive_ptr; iactive++) {
256-
update_one_coord_qp(Sigma_ptr,
256+
update_one_coord_qp(nndef_ptr,
257257
linear_func_ptr,
258-
Sigma_diag_ptr,
258+
nndef_diag_ptr,
259259
gradient_ptr,
260260
ever_active_ptr,
261261
nactive_ptr,
@@ -281,9 +281,9 @@ int solve_qp(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
281281

282282
for (icoord=0; icoord<nrow; icoord++) {
283283

284-
update_one_coord_qp(Sigma_ptr,
284+
update_one_coord_qp(nndef_ptr,
285285
linear_func_ptr,
286-
Sigma_diag_ptr,
286+
nndef_diag_ptr,
287287
gradient_ptr,
288288
ever_active_ptr,
289289
nactive_ptr,
@@ -313,7 +313,7 @@ int solve_qp(double *Sigma_ptr, /* A covariance matrix: X^TX/n */
313313
// Check relative decrease of objective
314314

315315
if (check_objective) {
316-
new_value = objective_qp(Sigma_ptr,
316+
new_value = objective_qp(nndef_ptr,
317317
linear_func_ptr,
318318
ever_active_ptr,
319319
nactive_ptr,

0 commit comments

Comments
 (0)