Skip to content

Commit eedf0fd

Browse files
jonathan-taylorjonathan-taylor
committed
created the matrices for the affine transform, wrapper for calling Gaussian density
1 parent e01d8f1 commit eedf0fd

File tree

7 files changed

+308
-3
lines changed

7 files changed

+308
-3
lines changed

‎selectiveInference/R/funs.fixed.R

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -435,8 +435,7 @@ debiasingRow = function (Xinfo, # could be X or t(X) %*% X / n dep
435435
last_output = NULL
436436

437437
if (is_wide) {
438-
n = nrow(Xinfo)
439-
Xsoln = as.numeric(rep(0, n))
438+
Xsoln = as.numeric(rep(0, nrow(Xinfo)))
440439
}
441440

442441
while (counter_idx < max_try) {

‎selectiveInference/R/funs.randomized.R

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,7 @@ randomizedLASSO = function(X,
3737
if (length(lam) == 1) {
3838
lam = rep(lam, p)
3939
}
40+
4041
if (length(lam) != p) {
4142
stop("Lagrange parameter should be single float or of length ncol(X)")
4243
}
@@ -65,5 +66,38 @@ randomizedLASSO = function(X,
6566
objective_stop, # objective_stop
6667
kkt_stop, # kkt_stop
6768
param_stop) # param_stop
69+
70+
71+
sign_soln = sign(result$soln)
72+
73+
unpenalized = lam == 0
74+
active = !unpenalized * (sign_soln != 0)
75+
inactive = !unpenzlied * (sign_soln == 0)
76+
77+
unpenalized_set = which(unpenalized)
78+
active_set = which(active)
79+
inactive_set = which(inactive)
80+
81+
coef_term = t(X) %*% X[,c(unpenalized_set, # the coefficients
82+
active_set)]
83+
coef_term = coef_term %*% diag(c(rep(1, sum(unpenalized)), sign_soln[active])) # coefficients are non-negative
84+
coef_term[active,] = coef_term[active,] + ridge_term * diag(rep(1, sum(active))) # ridge term
85+
86+
subgrad_term = cbind(matrix(0, sum(inactive), sum(active) + sum(unpenalized)),
87+
diag(rep(1, sum(inactive))))
88+
linear_term = rbind(coef_term,
89+
subgrad_term)
90+
91+
offset_term = rep(0, p)
92+
offset_term[active] = lam[active] * sign_soln[active]
93+
94+
95+
96+
list(active_set = active_set,
97+
inactive_set = inactive_set,
98+
unpenalized_set = unpenalized_set,
99+
sign_soln = sign_soln)
100+
101+
68102
return(result)
69103
}

‎selectiveInference/src/Rcpp-randomized.cpp

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,40 @@
1+
#include <Rcpp.h> // need to include the main Rcpp header file
2+
#include <randomized_lasso.h> // where densities are defined
3+
4+
// [[Rcpp::export]]
5+
Rcpp::NumericVector log_density_gaussian_(double noise_scale, // Scale of randomization
6+
Rcpp::NumericMatrix internal_linear, // A_D -- linear part for data
7+
Rcpp::NumericMatrix internal_state, // D -- data state -- matrix of shape (nopt, npts)
8+
Rcpp::NumericMatrix optimization_linear, // A_O -- linear part for optimization variables
9+
Rcpp::NumericMatrix optimization_state, // O -- optimization state -- matrix of shape (ninternal, npts)
10+
Rcpp::NumericMatrix offset) { // h -- offset in affine transform -- "p" dimensional
11+
12+
int npt = internal_state.ncol(); // Function is vectorized
13+
if (optimization_state.ncol() != npt) { // Assuming each column is an internal or opt state because arrays are column major
14+
Rcpp::stop("Number of optimization samples should equal the number of (internally represented) data.");
15+
}
16+
17+
int ndim = optimization_linear.nrow();
18+
if (internal_linear.nrow() != ndim) {
19+
Rcpp::stop("Dimension of optimization range should be the same as the dimension of the data range.");
20+
}
21+
int ninternal = internal_linear.ncol();
22+
int noptimization = optimization_linear.ncol();
23+
24+
Rcpp::NumericVector result(npt);
25+
26+
int ipt;
27+
for (ipt=0; ipt<npt; ipt++) {
28+
result[ipt] = log_density_gaussian(noise_scale,
29+
ndim,
30+
ninternal,
31+
noptimization,
32+
(double *) internal_linear.begin(),
33+
((double *) internal_state.begin() + ipt * ninternal),
34+
(double *) optimization_linear.begin(),
35+
((double *) optimization_state.begin() + ipt * noptimization),
36+
(double *) offset.begin());
37+
}
38+
39+
return(result);
40+
}

‎selectiveInference/src/randomized_lasso.c

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,8 @@
1111
// Laplace is product of IID Laplace with scale noise_scale
1212
// Also evaluated at A_D D + A_O O + h
1313

14+
// Matrices are assumed in column major order!
15+
1416
double log_density_gaussian(double noise_scale, // Scale of randomization
1517
int ndim, // Number of features -- "p"
1618
int ninternal, // Dimension of internal data representation often 1

‎selectiveInference/src/randomized_lasso.c~

Lines changed: 188 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,188 @@
1+
#include <math.h> // for fabs
2+
3+
// Augmented density for randomized LASSO after
4+
// Gaussian randomization
5+
6+
// Described in https://arxiv.org/abs/1609.05609
7+
8+
// Gaussian is product of IID N(0, noise_scale^2) density
9+
// Evaluated at A_D D + A_O O + h
10+
11+
// Laplace is product of IID Laplace with scale noise_scale
12+
// Also evaluated at A_D D + A_O O + h
13+
14+
double log_density_gaussian(double noise_scale, // Scale of randomization
15+
int ndim, // Number of features -- "p"
16+
int ninternal, // Dimension of internal data representation often 1
17+
int noptimization, // Dimension of optimization variables -- "p"
18+
double *internal_linear, // A_D -- linear part for data
19+
double *internal_state, // D -- data state
20+
double *optimization_linear, // A_O -- linear part for optimization variables
21+
double *optimization_state, // O -- optimization state
22+
double *offset) // h -- offset in affine transform -- "p" dimensional
23+
{
24+
int irow, icol;
25+
double denom = 2 * noise_scale * noise_scale;
26+
double value = 0;
27+
double reconstruction = 0;
28+
double *offset_ptr;
29+
double *internal_linear_ptr;
30+
double *internal_state_ptr;
31+
double *optimization_linear_ptr;
32+
double *optimization_state_ptr;
33+
34+
for (irow=0; irow<ndim; irow++) {
35+
36+
// Compute the irow-th entry of the ndim vector
37+
38+
offset_ptr = ((double *) offset + irow);
39+
reconstruction = *offset_ptr;
40+
41+
// Internal (i.e. data) contribution
42+
for (icol=0; icol<ninternal; icol++) {
43+
44+
internal_linear_ptr = ((double *) internal_linear + icol * ndim + irow);
45+
internal_state_ptr = ((double *) internal_state + icol);
46+
47+
reconstruction += (*internal_linear_ptr) * (*internal_state_ptr);
48+
}
49+
50+
// Optimization contribution
51+
for (icol=0; icol<noptimization; icol++) {
52+
53+
optimization_linear_ptr = ((double *) optimization_linear + icol * ndim + irow);
54+
optimization_state_ptr = ((double *) optimization_state + icol);
55+
56+
reconstruction += (*optimization_linear_ptr) * (*optimization_state_ptr);
57+
}
58+
59+
value -= (reconstruction * reconstruction) / denom;
60+
}
61+
62+
return(value);
63+
}
64+
65+
double log_density_laplace(double noise_scale, // Scale of randomization
66+
int ndim, // Number of features -- "p"
67+
int ninternal, // Dimension of internal data representation often 1
68+
int noptimization, // Dimension of optimization variables -- "p"
69+
double *internal_linear, // A_D -- linear part for data
70+
double *internal_state, // D -- data state
71+
double *optimization_linear, // A_O -- linear part for optimization variables
72+
double *optimization_state, // O -- optimization state
73+
double *offset) // h -- offset in affine transform -- "p" dimensional
74+
{
75+
int irow, icol;
76+
double value = 0;
77+
double reconstruction = 0;
78+
double *offset_ptr;
79+
double *internal_linear_ptr;
80+
double *internal_state_ptr;
81+
double *optimization_linear_ptr;
82+
double *optimization_state_ptr;
83+
84+
for (irow=0; irow<ndim; irow++) {
85+
86+
// Compute the irow-th entry of the ndim vector
87+
88+
offset_ptr = ((double *) offset + irow);
89+
reconstruction = *offset_ptr;
90+
91+
// Internal (i.e. data) contribution
92+
for (icol=0; icol<ninternal; icol++) {
93+
94+
internal_linear_ptr = ((double *) internal_linear + icol * ndim + irow);
95+
internal_state_ptr = ((double *) internal_state + icol);
96+
97+
reconstruction += (*internal_linear_ptr) * (*internal_state_ptr);
98+
}
99+
100+
// Optimization contribution
101+
for (icol=0; icol<noptimization; icol++) {
102+
103+
optimization_linear_ptr = ((double *) optimization_linear + icol * ndim + irow);
104+
optimization_state_ptr = ((double *) optimization_state + icol);
105+
106+
reconstruction += (*optimization_linear_ptr) * (*optimization_state_ptr);
107+
}
108+
109+
value -= fabs(reconstruction) / noise_scale;
110+
}
111+
112+
return(value);
113+
}
114+
115+
// Keeping internal (data) state fixed
116+
117+
double log_density_gaussian_conditional(double noise_scale, // Scale of randomization
118+
int ndim, // Number of features -- "p"
119+
int noptimization, // Dimension of optimization variables -- "p"
120+
double *optimization_linear, // A_O -- linear part for optimization variables
121+
double *optimization_state, // O -- optimization state
122+
double *offset) // h -- offset in affine transform -- "p" dimensional
123+
{
124+
int irow, icol;
125+
double value = 0;
126+
double denom = 2 * noise_scale * noise_scale;
127+
double reconstruction = 0;
128+
double *offset_ptr;
129+
double *optimization_linear_ptr;
130+
double *optimization_state_ptr;
131+
132+
for (irow=0; irow<ndim; irow++) {
133+
134+
// Compute the irow-th entry of the ndim vector
135+
136+
offset_ptr = ((double *) offset + irow);
137+
reconstruction = *offset_ptr;
138+
139+
// Optimization contribution
140+
for (icol=0; icol<noptimization; icol++) {
141+
142+
optimization_linear_ptr = ((double *) optimization_linear + icol * ndim + irow);
143+
optimization_state_ptr = ((double *) optimization_state + icol);
144+
145+
reconstruction += (*optimization_linear_ptr) * (*optimization_state_ptr);
146+
}
147+
148+
value -= reconstruction * reconstruction / denom;
149+
}
150+
151+
return(value);
152+
}
153+
154+
double log_density_laplace_conditional(double noise_scale, // Scale of randomization
155+
int ndim, // Number of features -- "p"
156+
int noptimization, // Dimension of optimization variables -- "p"
157+
double *optimization_linear, // A_O -- linear part for optimization variables
158+
double *optimization_state, // O -- optimization state
159+
double *offset) // h -- offset in affine transform -- "p" dimensional
160+
{
161+
int irow, icol;
162+
double value = 0;
163+
double reconstruction = 0;
164+
double *offset_ptr;
165+
double *optimization_linear_ptr;
166+
double *optimization_state_ptr;
167+
168+
for (irow=0; irow<ndim; irow++) {
169+
170+
// Compute the irow-th entry of the ndim vector
171+
172+
offset_ptr = ((double *) offset + irow);
173+
reconstruction = *offset_ptr;
174+
175+
// Optimization contribution
176+
for (icol=0; icol<noptimization; icol++) {
177+
178+
optimization_linear_ptr = ((double *) optimization_linear + icol * ndim + irow);
179+
optimization_state_ptr = ((double *) optimization_state + icol);
180+
181+
reconstruction += (*optimization_linear_ptr) * (*optimization_state_ptr);
182+
}
183+
184+
value -= fabs(reconstruction) / noise_scale;
185+
}
186+
187+
return(value);
188+
}

‎selectiveInference/src/randomized_lasso.h

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -35,7 +35,7 @@ double log_density_laplace_conditional(double noise_scale, // Scale
3535
int noptimization, // Dimension of optimization variables -- "p"
3636
double *optimization_linear, // A_O -- linear part for optimization variables
3737
double *optimization_state, // O -- optimization state
38-
double *offset). // h -- offset in affine transform -- "p" dimensional
38+
double *offset); // h -- offset in affine transform -- "p" dimensional
3939

4040
#ifdef __cplusplus
4141
} /* extern "C" */

‎selectiveInference/src/randomized_lasso.h~

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,42 @@
1+
#ifdef __cplusplus
2+
extern "C"
3+
{
4+
#endif /* __cplusplus */
5+
6+
double log_density_gaussian(double noise_scale, // Scale of randomization
7+
int ndim, // Number of features -- "p"
8+
int ninternal, // Dimension of internal data representation often 1
9+
int noptimization, // Dimension of optimization variables -- "p"
10+
double *internal_linear, // A_D -- linear part for data
11+
double *internal_state, // D -- data state
12+
double *optimization_linear, // A_O -- linear part for optimization variables
13+
double *optimization_state, // O -- optimization state
14+
double *offset); // h -- offset in affine transform -- "p" dimensional
15+
16+
double log_density_laplace(double noise_scale, // Scale of randomization
17+
int ndim, // Number of features -- "p"
18+
int ninternal, // Dimension of internal data representation often 1
19+
int noptimization, // Dimension of optimization variables -- "p"
20+
double *internal_linear, // A_D -- linear part for data
21+
double *internal_state, // D -- data state
22+
double *optimization_linear, // A_O -- linear part for optimization variables
23+
double *optimization_state, // O -- optimization state
24+
double *offset); // h -- offset in affine transform -- "p" dimensional
25+
26+
double log_density_gaussian_conditional(double noise_scale, // Scale of randomization
27+
int ndim, // Number of features -- "p"
28+
int noptimization, // Dimension of optimization variables -- "p"
29+
double *optimization_linear, // A_O -- linear part for optimization variables
30+
double *optimization_state, // O -- optimization state
31+
double *offset); // h -- offset in affine transform -- "p" dimensional
32+
33+
double log_density_laplace_conditional(double noise_scale, // Scale of randomization
34+
int ndim, // Number of features -- "p"
35+
int noptimization, // Dimension of optimization variables -- "p"
36+
double *optimization_linear, // A_O -- linear part for optimization variables
37+
double *optimization_state, // O -- optimization state
38+
double *offset). // h -- offset in affine transform -- "p" dimensional
39+
40+
#ifdef __cplusplus
41+
} /* extern "C" */
42+
#endif /* __cplusplus */

0 commit comments

Comments
 (0)