Skip to content

Commit c7c41e4

Browse files
jonathan-taylorjonathan-taylor
committed
optimized the C code a bit -- still has debug statements
1 parent 210d0fe commit c7c41e4

File tree

2 files changed

+119
-55
lines changed

2 files changed

+119
-55
lines changed

‎selectiveInference/R/funs.fixed.R

Lines changed: 14 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -348,17 +348,25 @@ InverseLinftyOneRowC <- function (Sigma, i, mu, maxiter=50) {
348348

349349
val = .C("find_one_row",
350350
Sigma=as.double(Sigma),
351+
Sigma_diag=as.double(diag(Sigma)),
352+
Sigma_theta=as.double(rep(0, p)),
351353
nrow=as.integer(p),
352354
bound=as.double(mu),
353355
theta=as.double(theta),
354-
maxiter=as.integer(maxiter),
356+
maxiter=as.integer(50),
355357
row=as.integer(i-1),
356358
coord=as.integer(i-1),
357359
dup=FALSE,
358360
package="selectiveInference")
359361

360362
# Check feasibility
361363

364+
# DEBUG statements
365+
#print(diag(Sigma))
366+
#print(0.5 * sum(val$theta * (Sigma %*% val$theta)) - val$theta[i] + mu * sum(abs(val$theta)))
367+
#print(Sigma %*% val$theta - val$Sigma_theta)
368+
#print(val$nrow) # number of iterations
369+
362370
if (max(abs(Sigma %*% val$theta - basis_vector)) > 1.01 * mu) {
363371
warning("Solution for row of M does not seem to be feasible")
364372
}
@@ -372,11 +380,11 @@ InverseLinftyOneRow <- function ( sigma, i, mu, maxiter=50, threshold=1e-2 ) {
372380
mu0 <- rho/(1+rho);
373381
beta <- rep(0,p);
374382

375-
if (mu >= mu0){
376-
beta[i] <- (1-mu0)/sigma[i,i];
377-
returnlist <- list("optsol" = beta, "iter" = 0);
378-
return(returnlist);
379-
}
383+
#if (mu >= mu0){
384+
# beta[i] <- (1-mu0)/sigma[i,i];
385+
# returnlist <- list("optsol" = beta, "iter" = 0);
386+
# return(returnlist);
387+
#}
380388

381389
diff.norm2 <- 1;
382390
last.norm2 <- 1;

‎selectiveInference/src/debiasing_matrix.c

Lines changed: 105 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -11,34 +11,72 @@
1111
// Therefore we don't have to negate the answer to get theta.
1212
// Update one coordinate
1313

14-
double update_one_coord(double *Sigma, /* A covariance matrix: X^TX/n */
15-
int nrow, /* How many rows in Sigma */
16-
double bound, /* feasibility parameter */
17-
double *theta, /* current value */
18-
int row, /* which row: 0-based */
19-
int coord) /* which coordinate to update: 0-based */
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 */
2019
{
20+
int irow, icol;
21+
double value = 0;
22+
double *Sigma_ptr = Sigma;
23+
double *theta_row_ptr, *theta_col_ptr;
2124

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+
for (icol=0; icol<nrow; icol++) {
31+
value += 0.5 * (*Sigma_ptr) * (*theta_row_ptr) * (*theta_col_ptr);
32+
Sigma_ptr++;
33+
theta_col_ptr++;
34+
}
35+
if (irow == row) {
36+
value -= (*theta_row_ptr); // the elementary basis vector term
37+
}
38+
39+
value = value + bound * fabs((*theta_row_ptr)); // the \ell_1 term
40+
theta_row_ptr++;
41+
}
42+
43+
return(value);
44+
}
45+
46+
47+
double update_one_coord(double *Sigma, /* A covariance matrix: X^TX/n */
48+
double *Sigma_diag, /* Diagonal entries of Sigma */
49+
double *Sigma_theta, /* Sigma times theta */
50+
int nrow, /* How many rows in Sigma */
51+
double bound, /* feasibility parameter */
52+
double *theta, /* current value */
53+
int row, /* which row: 0-based */
54+
int coord) /* which coordinate to update: 0-based */
55+
{
56+
57+
double delta;
2258
double linear_term = 0;
23-
double quadratic_term = 0;
2459
double value = 0;
60+
double old_value;
2561
double *Sigma_ptr;
26-
double *theta_ptr = theta;
62+
double *Sigma_theta_ptr;
63+
double *theta_ptr;
2764
int icol = 0;
2865

29-
Sigma_ptr = ((double *) Sigma + nrow * coord);
66+
double *quadratic_ptr = ((double *) Sigma_diag + coord);
67+
double quadratic_term = *quadratic_ptr;
68+
69+
Sigma_theta_ptr = ((double *) Sigma_theta + coord);
70+
linear_term = *Sigma_theta_ptr;
71+
72+
theta_ptr = ((double *) theta + coord);
73+
old_value = *theta_ptr;
74+
75+
// The coord entry of Sigma_theta term has a diagonal term in it:
76+
// Sigma[coord, coord] * theta[coord]
77+
// This removes it.
78+
linear_term -= quadratic_term * old_value;
3079

31-
for (icol=0; icol < nrow; icol++) {
32-
if (icol != coord) {
33-
linear_term += (*Sigma_ptr) * (*theta_ptr);
34-
}
35-
else {
36-
quadratic_term = (*Sigma_ptr);
37-
}
38-
Sigma_ptr += 1;
39-
theta_ptr += 1;
40-
}
41-
4280
if (row == coord) {
4381
linear_term -= 1;
4482
}
@@ -58,60 +96,72 @@ double update_one_coord(double *Sigma, /* A covariance matrix: X^TX/n */
5896
value = -(linear_term - bound) / quadratic_term;
5997
}
6098

99+
if (fabs(old_value - value) > 1.e-6 * (fabs(value) + fabs(old_value))) { // Update the linear term
100+
101+
delta = value - old_value;
102+
Sigma_ptr = ((double *) Sigma + coord * nrow);
103+
Sigma_theta_ptr = ((double *) Sigma_theta);
104+
105+
for (icol=0; icol<nrow; icol++) {
106+
*Sigma_theta_ptr = *Sigma_theta_ptr + delta * (*Sigma_ptr);
107+
Sigma_theta_ptr += 1;
108+
Sigma_ptr += 1;
109+
}
110+
111+
double before = objective(Sigma,
112+
nrow,
113+
row,
114+
bound,
115+
theta);
116+
fprintf(stderr, "before %f\n", before);
117+
61118
theta_ptr = ((double *) theta + coord);
62119
*theta_ptr = value;
63-
return(value);
64120

65-
}
121+
double after = objective(Sigma,
122+
nrow,
123+
row,
124+
bound,
125+
theta);
66126

67-
double objective(double *Sigma, /* A covariance matrix: X^TX/n */
68-
int nrow, /* how many rows in Sigma */
69-
int row, /* which row: 0-based */
70-
double bound, /* Lagrange multipler for \ell_1 */
71-
double *theta) /* current value */
72-
{
73-
int irow, icol;
74-
double value = 0;
75-
double *Sigma_ptr = Sigma;
76-
double *theta_row_ptr, *theta_col_ptr;
127+
fprintf(stderr, "after %f\n", after);
128+
if (after > before) {
129+
fprintf(stderr, "not a descent step!!!!!!!!!!!!!!!!!!!!!\n");
130+
}
77131

78-
theta_row_ptr = theta;
79-
theta_col_ptr = theta;
80132

81-
for (irow=0; irow<nrow; irow++) {
82-
double *theta_col_ptr = theta;
133+
}
134+
135+
Sigma_ptr = ((double *) Sigma + coord * nrow);
136+
Sigma_theta_ptr = ((double *) Sigma_theta);
83137
for (icol=0; icol<nrow; icol++) {
84-
value += (*Sigma_ptr) * (*theta_row_ptr) * (*theta_col_ptr);
85-
Sigma_ptr++;
86-
theta_col_ptr++;
87-
}
88-
if (irow == row) {
89-
value += (*theta_row_ptr); // the elementary basis vector term
138+
139+
Sigma_theta_ptr += 1;
140+
Sigma_ptr += 1;
90141
}
91142

92-
value = value + bound * fabs((*theta_row_ptr)); // the \ell_1 term
93-
theta_row_ptr++;
94-
}
95143

96144
return(value);
145+
97146
}
98147

99148
void find_one_row(double *Sigma, /* A covariance matrix: X^TX/n */
149+
double *Sigma_diag, /* Diagonal entry of covariance matrix */
150+
double *Sigma_theta, /* Sigma times theta */
100151
int *nrow_ptr, /* How many rows in Sigma */
101152
double *bound_ptr, /* feasibility parameter */
102153
double *theta, /* current value */
103154
int *maxiter_ptr, /* how many iterations */
104-
int *row_ptr, /* which row: 0-based */
105-
int *coord_ptr) /* which coordinate to update: 0-based */
155+
int *row_ptr) /* which coordinate to update: 0-based */
106156
{
107157

108158
int maxiter = *maxiter_ptr;
109159
int iter = 0;
110160
int icoord = 0;
111-
int coord = *coord_ptr;
112161
int row = *row_ptr;
113162
double bound = *bound_ptr;
114163
int nrow = *nrow_ptr;
164+
115165
double old_value = objective(Sigma,
116166
nrow,
117167
row,
@@ -125,6 +175,8 @@ void find_one_row(double *Sigma, /* A covariance matrix: X^TX/n */
125175
// Update the diagonal first
126176

127177
update_one_coord(Sigma,
178+
Sigma_diag,
179+
Sigma_theta,
128180
nrow,
129181
bound,
130182
theta,
@@ -134,6 +186,8 @@ void find_one_row(double *Sigma, /* A covariance matrix: X^TX/n */
134186
for (icoord=0; icoord<nrow; icoord++) {
135187

136188
update_one_coord(Sigma,
189+
Sigma_diag,
190+
Sigma_theta,
137191
nrow,
138192
bound,
139193
theta,
@@ -147,9 +201,11 @@ void find_one_row(double *Sigma, /* A covariance matrix: X^TX/n */
147201
bound,
148202
theta);
149203

150-
if (((old_value - new_value) < tol * fabs(new_value)) && (iter > 3)) {
204+
if (((old_value - new_value) < tol * fabs(new_value)) && (iter > 5)) {
151205
break;
152206
}
207+
208+
fprintf(stderr, "%f %f value\n", old_value, new_value);
153209
old_value = new_value;
154210
}
155211

0 commit comments

Comments
 (0)