Skip to content

Commit 3da6c50

Browse files
jonathan-taylorjonathan-taylor
committed
NF: density functions for randomized LASSO
1 parent 781609d commit 3da6c50

File tree

2 files changed

+243
-0
lines changed

2 files changed

+243
-0
lines changed

‎selectiveInference/src/randomized_lasso.c

Lines changed: 201 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,201 @@
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 *internal_linear_ptr;
130+
double *internal_state_ptr;
131+
double *optimization_linear_ptr;
132+
double *optimization_state_ptr;
133+
134+
for (irow=0; irow<ndim; irow++) {
135+
136+
// Compute the irow-th entry of the ndim vector
137+
138+
offset_ptr = ((double *) offset + irow);
139+
reconstruction = *offset_ptr;
140+
141+
// Optimization contribution
142+
for (icol=0; icol<noptimization; icol++) {
143+
144+
optimization_linear_ptr = ((double *) optimization_linear + icol * ndim + irow);
145+
optimization_state_ptr = ((double *) optimization_state + icol);
146+
147+
reconstruction += (*optimization_linear_ptr) * (*optimization_state_ptr);
148+
}
149+
150+
value -= reconstruction * reconstruction / denom;
151+
}
152+
153+
return(value);
154+
}
155+
156+
double log_density_laplace_conditional(double noise_scale, // Scale of randomization
157+
int ndim, // Number of features -- "p"
158+
int noptimization, // Dimension of optimization variables -- "p"
159+
double *optimization_linear, // A_O -- linear part for optimization variables
160+
double *optimization_state, // O -- optimization state
161+
double *offset) // h -- offset in affine transform -- "p" dimensional
162+
{
163+
int irow, icol;
164+
double value = 0;
165+
double reconstruction = 0;
166+
double *offset_ptr;
167+
double *internal_linear_ptr;
168+
double *internal_state_ptr;
169+
double *optimization_linear_ptr;
170+
double *optimization_state_ptr;
171+
172+
for (irow=0; irow<ndim; irow++) {
173+
174+
// Compute the irow-th entry of the ndim vector
175+
176+
offset_ptr = ((double *) offset + irow);
177+
reconstruction = *offset_ptr;
178+
179+
// Internal (i.e. data) contribution
180+
for (icol=0; icol<ninternal; icol++) {
181+
182+
internal_linear_ptr = ((double *) internal_linear + icol * ndim + irow);
183+
internal_state_ptr = ((double *) internal_state + icol);
184+
185+
reconstruction += (*internal_linear_ptr) * (*internal_state_ptr);
186+
}
187+
188+
// Optimization contribution
189+
for (icol=0; icol<noptimization; icol++) {
190+
191+
optimization_linear_ptr = ((double *) optimization_linear + icol * ndim + irow);
192+
optimization_state_ptr = ((double *) optimization_state + icol);
193+
194+
reconstruction += (*optimization_linear_ptr) * (*optimization_state_ptr);
195+
}
196+
197+
value -= fabs(reconstruction) / noise_scale;
198+
}
199+
200+
return(value);
201+
}

‎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)