Source Code for Module parsimony.functions.nesterov.l1

  1  # -*- coding: utf-8 -*- 
  2  """ 
  3  The :mod:`parsimony.functions.nesterov.L1` module contains the loss function 
  4  for the L1 penalty, smoothed using Nesterov's technique. 
  5   
  6  Created on Mon Feb  3 17:00:56 2014 
  7   
  8  Copyright (c) 2013-2014, CEA/DSV/I2BM/Neurospin. All rights reserved. 
  9   
 10  @author:  Tommy Löfstedt, Vincent Guillemot, Edouard Duchesnay and 
 11            Fouad Hadj-Selem 
 12  @email:   lofstedt.tommy@gmail.com, edouard.duchesnay@cea.fr 
 13  @license: BSD 3-clause. 
 14  """ 
 15  import numpy as np 
 16  import scipy.sparse as sparse 
 17   
 18  #from .properties import NesterovFunction 
 19  from .. import properties 
 20  import parsimony.utils.consts as consts 
 21  import parsimony.utils.maths as maths 
 22  import parsimony.utils as utils 
 23   
 24  __all__ = ["L1", "linear_operator_from_variables", "A_from_variables"] 
 25   
 26   
 27 -class L1(properties.NesterovFunction, 
 28           properties.Penalty, 
 29           properties.Constraint): 
 30      """The proximal operator of the smoothed L1 function 
 31   
 32          f(beta) = l * (L1mu(beta) - c), 
 33   
 34      where L1mu(\beta) is the smoothed L1 function. The constrained version has 
 35      the form 
 36   
 37          L1mu(beta) <= c. 
 38      """ 
 39 -    def __init__(self, l, c=0.0, A=None, mu=0.0, penalty_start=0): 
 40          """ 
 41          Parameters 
 42          ---------- 
 43          l : Non-negative float. The Lagrange multiplier, or regularisation 
 44                  constant, of the function. 
 45   
 46          c : Float. The limit of the constraint. The function is feasible if 
 47                  ||beta||_1 <= c. The default value is c=0, i.e. the default is 
 48                  a regularisation formulation. 
 49   
 50          A : A (usually sparse) matrix. The linear operator for the Nesterov 
 51                  formulation. May not be None. 
 52   
 53          mu : Non-negative float. The regularisation constant for the smoothing. 
 54   
 55          penalty_start : Non-negative integer. The number of columns, variables 
 56                  etc., to exempt from penalisation. Equivalently, the first 
 57                  index to be penalised. Default is 0, all columns are included. 
 58          """ 
 59          super(L1, self).__init__(l, A=A, mu=mu, penalty_start=penalty_start) 
 60   
 61          self.c = float(c) 
 62   
 63 -    def f(self, beta): 
 64          """ Function value. 
 65          """ 
 66          if self.l < consts.TOLERANCE: 
 67              return 0.0 
 68   
 69          if self.penalty_start > 0: 
 70              beta_ = beta[self.penalty_start:, :] 
 71          else: 
 72              beta_ = beta 
 73   
 74          return self.l * (maths.norm1(beta_) - self.c) 
 75   
 76 -    def fmu(self, beta, mu=None): 
 77          """Returns the smoothed function value. 
 78   
 79          Parameters: 
 80          ---------- 
 81          beta : A weight vector. 
 82   
 83          mu : The regularisation constant for the smoothing. 
 84          """ 
 85          fmu = super(L1, self).fmu(beta, mu=mu) 
 86   
 87          return fmu - self.l * self.c 
 88   
 89  #        if mu is None: 
 90  #            mu = self.get_mu() 
 91  # 
 92  #        alpha = self.alpha(beta) 
 93  #        alpha_sqsum = 0.0 
 94  #        for a in alpha: 
 95  #            alpha_sqsum += np.sum(a ** 2.0) 
 96  # 
 97  #        Aa = self.Aa(alpha) 
 98  # 
 99  #        return self.l * ((np.dot(beta.T, Aa)[0, 0] 
100  #                          - (mu / 2.0) * alpha_sqsum) - self.c) 
101   
102 -    def phi(self, alpha, beta): 
103          """ Function value with known alpha. 
104   
105          From the interface "NesterovFunction". 
106          """ 
107          if self.l < consts.TOLERANCE: 
108              return 0.0 
109   
110          if self.penalty_start > 0: 
111              beta_ = beta[self.penalty_start:, :] 
112          else: 
113              beta_ = beta 
114   
115          return self.l * ((np.dot(alpha[0].T, beta_)[0, 0] 
116                           - (self.mu / 2.0) * np.sum(alpha[0] ** 2.0)) - self.c) 
117   
118 -    def grad(self, beta): 
119          """ Gradient of the function at beta. 
120   
121          From the interface "Gradient". Overloaded since we can do it faster 
122          than the default. 
123          """ 
124          alpha = self.alpha(beta) 
125   
126          return self.l * alpha[0] 
127   
128 -    def L(self): 
129          """ Lipschitz constant of the gradient. 
130   
131          From the interface "LipschitzContinuousGradient". 
132          """ 
133          return self.l / self.mu 
134   
135 -    def alpha(self, beta): 
136          """ Dual variable of the Nesterov function. 
137   
138          From the interface "NesterovFunction". Overloaded since we can do it 
139          faster than the default. 
140          """ 
141          if self.penalty_start > 0: 
142              beta_ = beta[self.penalty_start:, :] 
143          else: 
144              beta_ = beta 
145   
146          alpha = self.project([beta_ * (1.0 / self.mu)]) 
147   
148          return alpha 
149   
150      @staticmethod 
151 -    def project(a): 
152          """ Projection onto the compact space of the Nesterov function. 
153   
154          From the interface "NesterovFunction". 
155          """ 
156          a = a[0] 
157          anorm = np.abs(a) 
158          i = anorm > 1.0 
159          anorm_i = anorm[i] 
160          a[i] = np.divide(a[i], anorm_i) 
161   
162          return [a] 
163   
164 -    def M(self): 
165          """ The maximum value of the regularisation of the dual variable. We 
166          have 
167   
168              M = max_{alpha in K} 0.5*|alpha|²_2. 
169   
170          From the interface "NesterovFunction". 
171          """ 
172          A = self.A() 
173          return A[0].shape[0] / 2.0 
174   
175 -    def estimate_mu(self, beta): 
176          """ Computes a "good" value of mu with respect to the given \beta. 
177   
178          From the interface "NesterovFunction". 
179          """ 
180          if self.penalty_start > 0: 
181              beta_ = beta[self.penalty_start:, :] 
182          else: 
183              beta_ = beta 
184   
185          return np.max(np.absolute(beta_)) 
186   
187 -    def feasible(self, beta): 
188          """Feasibility of the constraint. 
189   
190          From the interface "Constraint". 
191          """ 
192          if self.penalty_start > 0: 
193              beta_ = beta[self.penalty_start:, :] 
194          else: 
195              beta_ = beta 
196   
197          return maths.norm1(beta_) <= self.c 
198   
199   
200  @utils.deprecated("linear_operator_from_variables") 
201 -def A_from_variables(*args, **kwargs): 
202   
203      return linear_operator_from_variables(*args, **kwargs) 
204   
205   
206 -def linear_operator_from_variables(num_variables, penalty_start=0): 
207      """Generates the linear operator for the L1 Nesterov function from number 
208      of variables. 
209   
210      Parameters: 
211      ---------- 
212      num_variables : Positive integer. The total number of variables, including 
213              the intercept variable(s). 
214   
215      penalty_start : Non-negative integer. The number of variables to exempt 
216              from penalisation. Equivalently, the first index to be penalised. 
217              Default is 0, all variables are included. 
218      """ 
219      A = sparse.eye(num_variables - penalty_start, 
220                     num_variables - penalty_start, 
221                     format="csr") 
222   
223      return [A] 
224