1
2 """
3 Created on Mon Jan 27 16:23:54 2014
4
5 Copyright (c) 2013-2014, CEA/DSV/I2BM/Neurospin. All rights reserved.
6
7 @author: Tommy Löfstedt
8 @email: tommy.loefstedt@cea.fr
9 @license: BSD 3-clause.
10 """
11 import numpy as np
12 from grad import grad_l1
13 from grad import grad_l2_squared
14 from grad import grad_gl
15 from utils import bisection_method
16
17 __all__ = ['load']
18
19
20 -def load(l, k, g, beta, M, e, A, snr=None, intercept=False):
21 """Returns data generated such that we know the exact solution.
22
23 The data generated by this function is fit to the Linear regression + L1 +
24 L2 + Group lasso function, i.e. to:
25
26 f(b) = (1 / 2).|Xb - y|² + l.|b|_1 + (k / 2).|b|² + g.GL(b),
27
28 where |.|_1 is the L1 norm, |.|² is the squared L2 norm and GL is the
29 group lasso penalty.
30
31 Parameters
32 ----------
33 l : Non-negative float. The L1 regularisation parameter.
34
35 k : Non-negative float. The L2 regularisation parameter.
36
37 g : Non-negative float. The group lasso regularisation parameter.
38
39 beta : Numpy array (p-by-1). The regression vector to generate data from.
40
41 M : Numpy array (n-by-p). The matrix to use when building data. This
42 matrix carries the desired correlation structure of the generated
43 data. The generated data will be a column-scaled version of this
44 matrix.
45
46 e : Numpy array (n-by-1). The error vector e = Xb - y. This vector carries
47 the desired distribution of the residual.
48
49 A : Numpy or (usually) scipy.sparse array (K-by-p). The linear operator
50 for the Nesterov function.
51
52 snr : Positive float. Signal-to-noise ratio between model and residual.
53
54 intercept : Boolean. Whether or not to include an intercept variable. This
55 variable is not penalised. Note that if intercept is True, then e
56 will be centred.
57
58 Returns
59 -------
60 X : Numpy array (n-by-p). The generated X matrix.
61
62 y : Numpy array (n-by-1). The generated y vector.
63
64 beta : Numpy array (p-by-1). The regression vector with the correct snr.
65 """
66 l = float(l)
67 k = float(k)
68 g = float(g)
69
70 if intercept:
71 e = e - np.mean(e)
72
73 if snr != None:
74 def f(x):
75 X, y = _generate(l, k, g, x * beta, M, e, A, intercept)
76
77
78
79
80
81
82 return (np.linalg.norm(np.dot(X, x * beta)) / np.linalg.norm(e)) \
83 - snr
84
85 snr = bisection_method(f, low=0.0, high=np.sqrt(snr), maxiter=30)
86
87 beta = beta * snr
88
89 X, y = _generate(l, k, g, beta, M, e, A, intercept)
90
91 return X, y, beta
92
93
94 -def _generate(l, k, g, beta, M, e, A, intercept):
95
96 p = beta.shape[0]
97
98 if intercept:
99 gradL1 = grad_l1(beta[1:, :])
100 gradL2 = grad_l2_squared(beta[1:, :])
101 gradGL = grad_gl(beta[1:, :], A)
102 else:
103 gradL1 = grad_l1(beta)
104 gradL2 = grad_l2_squared(beta)
105 gradGL = grad_gl(beta, A)
106
107 alpha = -(l * gradL1 + k * gradL2 + g * gradGL)
108 Mte = np.dot(M.T, e)
109 if intercept:
110 alpha = np.divide(alpha, Mte[1:, :])
111 else:
112 alpha = np.divide(alpha, Mte)
113
114 X = np.ones(M.shape)
115 if intercept:
116 for i in xrange(p - 1):
117 X[:, i + 1] = M[:, i + 1] * alpha[i, 0]
118 else:
119 for i in xrange(p):
120 X[:, i] = M[:, i] * alpha[i, 0]
121
122 y = np.dot(X, beta) - e
123
124 return X, y
125