Source Code for Module parsimony.datasets.simulate.grad

  1  # -*- coding: utf-8 -*- 
  2  """ 
  3  Created on Thu Sep 26 12:06:07 2013 
  4   
  5  Copyright (c) 2013-2014, CEA/DSV/I2BM/Neurospin. All rights reserved. 
  6   
  7  @author:  Tommy Löfstedt, Edouard Duchesnay 
  8  @email:   tommy.loefstedt@cea.fr, edouard.duchesnay@cea.fr 
  9  @license: BSD 3-clause. 
 10  """ 
 11  import abc 
 12   
 13  import numpy as np 
 14   
 15  from utils import TOLERANCE 
 16  from utils import RandomUniform 
 17  from utils import norm2 
 18   
 19  __all__ = ["grad_l1", "grad_l1mu", "grad_l2", "grad_l2", "grad_l2_squared", 
 20             "grad_tv", "grad_tvmu", "grad_grouptvmu"] 
 21   
 22   
 23 -class Function(object): 
 24   
 25      __metaclass__ = abc.ABCMeta 
 26   
 27 -    def __init__(self, l, **kwargs): 
 28   
 29          self.l = float(l) 
 30   
 31          for k in kwargs: 
 32              setattr(self, k, kwargs[k]) 
 33   
 34      @abc.abstractmethod 
 35 -    def grad(self, x): 
 36          raise NotImplementedError("Abstract method 'grad' must be " 
 37                                    "specialised!") 
 38   
 39   
 40 -class L1(Function): 
 41   
 42 -    def __init__(self, l, rng=RandomUniform(-1, 1)): 
 43   
 44          super(L1, self).__init__(l, rng=rng) 
 45   
 46 -    def grad(self, x): 
 47          """Sub-gradient of the function 
 48   
 49              f(x) = |x|_1, 
 50   
 51          where |x|_1 is the L1-norm. 
 52          """ 
 53          grad = np.zeros((x.shape[0], 1)) 
 54          grad[x >= TOLERANCE] = 1.0 
 55          grad[x <= -TOLERANCE] = -1.0 
 56          between = (x > -TOLERANCE) & (x < TOLERANCE) 
 57          grad[between] = self.rng(between.sum()) 
 58   
 59          return self.l * grad 
 60   
 61   
 62 -def grad_l1(beta, rng=RandomUniform(-1, 1)): 
 63      """Sub-gradient of the function 
 64   
 65          f(x) = |x|_1, 
 66   
 67      where |x|_1 is the L1-norm. 
 68      """ 
 69      grad = np.zeros((beta.shape[0], 1)) 
 70      grad[beta >= TOLERANCE] = 1.0 
 71      grad[beta <= -TOLERANCE] = -1.0 
 72      between = (beta > -TOLERANCE) & (beta < TOLERANCE) 
 73      grad[between] = rng(between.sum()) 
 74   
 75      return grad 
 76   
 77   
 78 -class SmoothedL1(Function): 
 79   
 80 -    def __init__(self, l, mu=TOLERANCE): 
 81   
 82          super(SmoothedL1, self).__init__(l, mu=mu) 
 83   
 84 -    def grad(self, x): 
 85          """Gradient of the function 
 86   
 87              f(x) = L1(mu, x), 
 88   
 89          where L1(mu, x) is the Nesterov smoothed L1-norm. 
 90          """ 
 91          alpha = (1.0 / self.mu) * x 
 92          asnorm = np.abs(alpha) 
 93          i = asnorm > 1.0 
 94          alpha[i] = np.divide(alpha[i], asnorm[i]) 
 95   
 96          return self.l * alpha 
 97   
 98   
 99 -def grad_l1mu(beta, mu): 
100      """Gradient of the function 
101   
102          f(x) = L1(mu, x), 
103   
104      where L1(mu, x) is the Nesterov smoothed L1-norm. 
105      """ 
106      alpha = (1.0 / mu) * beta 
107      asnorm = np.abs(alpha) 
108      i = asnorm > 1.0 
109      alpha[i] = np.divide(alpha[i], asnorm[i]) 
110   
111      return alpha 
112   
113   
114 -class L2(Function): 
115   
116 -    def __init__(self, l, rng=RandomUniform(0, 1)): 
117   
118          super(L2, self).__init__(l, rng=rng) 
119   
120 -    def grad(self, x): 
121          """Sub-gradient of the function 
122   
123              f(x) = |x|_2, 
124   
125          where |x|_2 is the L2-norm. 
126          """ 
127          norm_beta = norm2(x) 
128          if norm_beta > TOLERANCE: 
129              return x * (1.0 / norm_beta) 
130          else: 
131              D = x.shape[0] 
132              u = (self.rng(D, 1) * 2.0) - 1.0  # [-1, 1]^D 
133              norm_u = norm2(u) 
134              a = self.rng()  # [0, 1] 
135   
136              return (self.l * (a / norm_u)) * u 
137   
138   
139 -def grad_l2(beta, rng=RandomUniform(0, 1)): 
140      """Sub-gradient of the function 
141   
142          f(x) = |x|_2, 
143   
144      where |x|_2 is the L2-norm. 
145      """ 
146      norm_beta = norm2(beta) 
147      if norm_beta > TOLERANCE: 
148          return beta * (1.0 / norm_beta) 
149      else: 
150          D = beta.shape[0] 
151          u = (rng(D, 1) * 2.0) - 1.0  # [-1, 1]^D 
152          norm_u = norm2(u) 
153          a = rng()  # [0, 1] 
154   
155          return u * (a / norm_u) 
156   
157   
158 -class L2Squared(Function): 
159   
160 -    def __init__(self, l): 
161   
162          super(L2Squared, self).__init__(l) 
163   
164 -    def grad(self, x): 
165          """Gradient of the function 
166   
167              f(x) = (1 / 2) * |x|²_2, 
168   
169          where |x|²_2 is the squared L2-norm. 
170          """ 
171          return self.l * x 
172   
173   
174 -def grad_l2_squared(beta, rng=None): 
175      """Gradient of the function 
176   
177          f(x) = (1 / 2) * |x|²_2, 
178   
179      where |x|²_2 is the squared L2-norm. 
180      """ 
181      return beta 
182   
183   
184 -class NesterovFunction(Function): 
185   
186      __metaclass__ = abc.ABCMeta 
187   
188 -    def __init__(self, l, A, mu=TOLERANCE, rng=RandomUniform(-1, 1), 
189                   norm=L2.grad, **kwargs): 
190   
191          super(NesterovFunction, self).__init__(l, rng=rng, norm=norm, **kwargs) 
192   
193          self.A = A 
194          self.mu = mu 
195   
196 -    def grad(self, x): 
197   
198          grad_Ab = 0 
199          for i in xrange(len(self.A)): 
200              Ai = self.A[i] 
201              Ab = Ai.dot(x) 
202              grad_Ab += Ai.T.dot(self.norm(Ab, self.rng)) 
203   
204          return self.l * grad_Ab 
205   
206 -    def smoothed_grad(self, x): 
207   
208          alpha = self.alpha(x) 
209   
210          Aa = self.A[0].T.dot(alpha[0]) 
211          for i in xrange(1, len(self.A)): 
212              Aa += self.A[i].T.dot(alpha[i]) 
213   
214          return self.l * Aa 
215   
216 -    def alpha(self, x): 
217          """ Dual variable of the Nesterov function. 
218          """ 
219          alpha = [0] * len(self.A) 
220          for i in xrange(len(self.A)): 
221              alpha[i] = self.A[i].dot(x) * (1.0 / self.mu) 
222   
223          # Apply projection 
224          alpha = self.project(alpha) 
225   
226          return alpha 
227   
228 -    def project(self, alpha): 
229   
230          for i in xrange(len(alpha)): 
231              astar = alpha[i] 
232              normas = np.sqrt(np.sum(astar ** 2.0)) 
233              if normas > 1.0: 
234                  astar *= 1.0 / normas 
235              alpha[i] = astar 
236   
237          return alpha 
238   
239   
240 -class TotalVariation(Function): 
241   
242 -    def __init__(self, l, A, rng=RandomUniform(0, 1)): 
243   
244          super(TotalVariation, self).__init__(l, A=A, rng=rng) 
245   
246 -    def grad(self, x): 
247          """Gradient of the function 
248   
249              f(x) = TV(x), 
250   
251          where TV(x) is the total variation function. 
252          """ 
253          beta_flat = x.ravel() 
254          Ab = np.vstack([Ai.dot(beta_flat) for Ai in self.A]).T 
255          Ab_norm2 = np.sqrt(np.sum(Ab ** 2.0, axis=1)) 
256   
257          upper = Ab_norm2 > TOLERANCE 
258          grad_Ab_norm2 = Ab 
259          grad_Ab_norm2[upper] = (Ab[upper].T / Ab_norm2[upper]).T 
260   
261          lower = Ab_norm2 <= TOLERANCE 
262          n_lower = lower.sum() 
263   
264          if n_lower: 
265              D = len(self.A) 
266              vec_rnd = (self.rng(n_lower, D) * 2.0) - 1.0 
267              norm_vec = np.sqrt(np.sum(vec_rnd ** 2.0, axis=1)) 
268              a = self.rng(n_lower) 
269              grad_Ab_norm2[lower] = (vec_rnd.T * (a / norm_vec)).T 
270   
271          grad = np.vstack([self.A[i].T.dot(grad_Ab_norm2[:, i]) \ 
272                            for i in xrange(len(self.A))]) 
273          grad = grad.sum(axis=0) 
274   
275          return self.l * grad.reshape(x.shape) 
276   
277   
278 -def grad_tv(beta, A, rng=RandomUniform(0, 1)): 
279      beta_flat = beta.ravel() 
280      Ab = np.vstack([Ai.dot(beta_flat) for Ai in A]).T 
281      Ab_norm2 = np.sqrt(np.sum(Ab ** 2.0, axis=1)) 
282   
283      upper = Ab_norm2 > TOLERANCE 
284      grad_Ab_norm2 = Ab 
285      grad_Ab_norm2[upper] = (Ab[upper].T / Ab_norm2[upper]).T 
286   
287      lower = Ab_norm2 <= TOLERANCE 
288      n_lower = lower.sum() 
289   
290      if n_lower: 
291          D = len(A) 
292          vec_rnd = (rng(n_lower, D) * 2.0) - 1.0 
293          norm_vec = np.sqrt(np.sum(vec_rnd ** 2.0, axis=1)) 
294          a = rng(n_lower) 
295          grad_Ab_norm2[lower] = (vec_rnd.T * (a / norm_vec)).T 
296   
297      grad = np.vstack([A[i].T.dot(grad_Ab_norm2[:, i]) for i in xrange(len(A))]) 
298      grad = grad.sum(axis=0) 
299   
300      return grad.reshape(beta.shape) 
301   
302   
303 -class GroupLasso(Function): 
304   
305 -    def __init__(self, l, A, rng=RandomUniform(-1, 1)): 
306   
307          super(GroupLasso, self).__init__(l, A, rng=rng) 
308   
309   
310 -def grad_gl(beta, A, rng=RandomUniform(-1, 1)): 
311   
312      return _Nesterov_grad(beta, A, rng, grad_l2) 
313   
314   
315 -class SmoothedTotalVariation(NesterovFunction): 
316   
317 -    def __init__(self, l, A, mu=TOLERANCE): 
318   
319          super(SmoothedTotalVariation, self).__init__(l, A, mu=mu) 
320   
321 -    def grad(self, x): 
322          """Gradient of the function 
323   
324              f(x) = TV(mu, x), 
325   
326          where TV(mu, x) is the Nesterov smoothed total variation function. 
327          """ 
328          return self.smoothed_grad(x) 
329   
330 -    def project(self, alpha): 
331          """ Projection onto the compact space of the smoothed TV function. 
332          """ 
333          ax = alpha[0] 
334          ay = alpha[1] 
335          az = alpha[2] 
336          anorm = ax ** 2.0 + ay ** 2.0 + az ** 2.0 
337          i = anorm > 1.0 
338   
339          anorm_i = anorm[i] ** 0.5  # Square root is taken here. Faster. 
340          ax[i] = np.divide(ax[i], anorm_i) 
341          ay[i] = np.divide(ay[i], anorm_i) 
342          az[i] = np.divide(az[i], anorm_i) 
343   
344          return [ax, ay, az] 
345   
346   
347 -def grad_tvmu(beta, A, mu): 
348   
349      alpha = _Nestetov_alpha(beta, A, mu, _Nesterov_TV_project) 
350   
351      return _Nesterov_grad_smoothed(A, alpha) 
352   
353   
354 -class SmoothedGroupLasso(NesterovFunction): 
355   
356 -    def __init__(self, l, A, mu=TOLERANCE): 
357   
358          super(SmoothedGroupLasso, self).__init__(l, A, mu=mu) 
359   
360 -    def grad(self, x): 
361          """Gradient of the function 
362   
363              f(x) = GL(mu, x), 
364   
365          where GL(mu, x) is the Nesterov smoothed group lasso function. 
366          """ 
367          return self.smoothed_grad(x) 
368   
369   
370 -def grad_glmu(beta, A, mu): 
371   
372      alpha = _Nestetov_alpha(beta, A, mu, _Nesterov_project) 
373   
374      return _Nesterov_grad_smoothed(A, alpha) 
375   
376   
377 -class SmoothedGroupTotalVariation(NesterovFunction): 
378   
379 -    def __init__(self, l, A, mu=TOLERANCE): 
380   
381          super(SmoothedGroupTotalVariation, self).__init__(l, A, mu=mu) 
382   
383 -    def grad(self, x): 
384          """Gradient of the function 
385   
386              f(x) = GroupTV(mu, x), 
387   
388          where GroupTV(mu, x) is the Nesterov smoothed group total variation 
389          function. 
390          """ 
391          return self.smoothed_grad(x) 
392   
393 -    def project(self, a): 
394          """ Projection onto the compact space of the smoothed Group TV 
395          function. 
396          """ 
397          for g in xrange(0, len(a), 3): 
398   
399              ax = a[g + 0] 
400              ay = a[g + 1] 
401              az = a[g + 2] 
402              anorm = ax ** 2.0 + ay ** 2.0 + az ** 2.0 
403              i = anorm > 1.0 
404   
405              anorm_i = anorm[i] ** 0.5  # Square root is taken here. Faster. 
406              ax[i] = np.divide(ax[i], anorm_i) 
407              ay[i] = np.divide(ay[i], anorm_i) 
408              az[i] = np.divide(az[i], anorm_i) 
409   
410              a[g + 0] = ax 
411              a[g + 1] = ay 
412              a[g + 2] = az 
413   
414          return a 
415   
416   
417 -def grad_grouptvmu(beta, A, mu): 
418   
419      alpha = _Nestetov_alpha(beta, A, mu, _Nesterov_GroupTV_project) 
420   
421      return _Nesterov_grad_smoothed(A, alpha) 
422   
423   
424 -def _Nesterov_GroupTV_project(a): 
425      """ Projection onto the compact space of the smoothed Group TV function. 
426      """ 
427      for g in xrange(0, len(a), 3): 
428   
429          ax = a[g + 0] 
430          ay = a[g + 1] 
431          az = a[g + 2] 
432          anorm = ax ** 2.0 + ay ** 2.0 + az ** 2.0 
433          i = anorm > 1.0 
434   
435          anorm_i = anorm[i] ** 0.5  # Square root is taken here. Faster. 
436          ax[i] = np.divide(ax[i], anorm_i) 
437          ay[i] = np.divide(ay[i], anorm_i) 
438          az[i] = np.divide(az[i], anorm_i) 
439   
440          a[g + 0] = ax 
441          a[g + 1] = ay 
442          a[g + 2] = az 
443   
444      return a 
445   
446   
447 -def _Nesterov_grad(beta, A, rng=RandomUniform(-1, 1), grad_norm=grad_l2): 
448   
449      grad_Ab = 0 
450      for i in xrange(len(A)): 
451          Ai = A[i] 
452          Ab = Ai.dot(beta) 
453          grad_Ab += Ai.T.dot(grad_norm(Ab, rng)) 
454   
455      return grad_Ab 
456   
457   
458 -def _Nesterov_grad_smoothed(A, alpha): 
459   
460      Aa = A[0].T.dot(alpha[0]) 
461      for i in xrange(1, len(A)): 
462          Aa += A[i].T.dot(alpha[i]) 
463   
464      return Aa 
465   
466   
467 -def _Nestetov_alpha(beta, A, mu, proj): 
468      """ Dual variable of the Nesterov function. 
469      """ 
470      alpha = [0] * len(A) 
471      for i in xrange(len(A)): 
472          alpha[i] = A[i].dot(beta) * (1.0 / mu) 
473   
474      # Apply projection. 
475      alpha = proj(alpha) 
476   
477      return alpha 
478   
479   
480 -def _Nesterov_project(alpha): 
481   
482      for i in xrange(len(alpha)): 
483          astar = alpha[i] 
484          normas = np.sqrt(np.sum(astar ** 2.0)) 
485          if normas > 1.0: 
486              astar *= 1.0 / normas 
487          alpha[i] = astar 
488   
489      return alpha 
490   
491   
492 -def _Nesterov_TV_project(alpha): 
493      """ Projection onto the compact space of the smoothed TV function. 
494      """ 
495      ax = alpha[0] 
496      ay = alpha[1] 
497      az = alpha[2] 
498      anorm = ax ** 2.0 + ay ** 2.0 + az ** 2.0 
499      i = anorm > 1.0 
500   
501      anorm_i = anorm[i] ** 0.5  # Square root is taken here. Faster. 
502      ax[i] = np.divide(ax[i], anorm_i) 
503      ay[i] = np.divide(ay[i], anorm_i) 
504      az[i] = np.divide(az[i], anorm_i) 
505   
506      return [ax, ay, az] 
507   
508   
509  if __name__ == "__main__": 
510      import doctest 
511      doctest.testmod() 
512