Source Code for Module parsimony.functions.nesterov.l1tv

  1  # -*- coding: utf-8 -*- 
  2  """ 
  3  The :mod:`parsimony.functions.nesterov.L1TV` module contains the loss function 
  4  for the L1 + TV penalty, smoothed together using Nesterov's technique. 
  5   
  6  Created on Mon Feb  3 17:04:14 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 math 
 16   
 17  import scipy.sparse as sparse 
 18  import numpy as np 
 19   
 20  #from .properties import NesterovFunction 
 21  from .. import properties 
 22  import parsimony.utils.consts as consts 
 23  import parsimony.utils.maths as maths 
 24  import parsimony.utils as utils 
 25  import tv 
 26  import l1 
 27   
 28  __all__ = ["L1TV", 
 29             "linear_operator_from_mask", "A_from_mask", 
 30             "linear_operator_from_shape", "A_from_shape"] 
 31   
 32   
 33 -class L1TV(properties.NesterovFunction, 
 34             properties.Penalty, 
 35             properties.Eigenvalues): 
 36      """The proximal operator of the smoothed sum of the TV and L1 functions 
 37   
 38          f(beta) = (l1 * L1(beta) + tv * TV(beta))_mu, 
 39   
 40      where (...)_mu means that what's within parentheses is smoothed. 
 41      """ 
 42 -    def __init__(self, l1, tv, A=None, mu=0.0, penalty_start=0): 
 43          """ 
 44          Parameters 
 45          ---------- 
 46          l1 : Non-negative float. The Lagrange multiplier, or regularisation 
 47                  constant, of the smoothed L1 part of the function. 
 48   
 49          tv : Non-negative float. The Lagrange multiplier, or regularisation 
 50                  constant, of the smoothed total variation part of the function. 
 51   
 52          A : A list or tuple with 4 elements of (usually sparse) arrays. The 
 53                  linear operator for the smoothed L1+TV. The first element must 
 54                  be the linear operator for L1 and the following three for TV. 
 55                  May not be None. 
 56   
 57          mu : Non-negative float. The regularisation constant for the smoothing. 
 58   
 59          penalty_start : Non-negative integer. The number of columns, variables 
 60                  etc., to exempt from penalisation. Equivalently, the first 
 61                  index to be penalised. Default is 0, all columns are included. 
 62          """ 
 63          self.g = float(tv) 
 64   
 65          # WARNING: Number of non-zero rows may differ from p. 
 66          self._p = A[0].shape[1] 
 67          A = [l1 * A[0], 
 68               tv * A[1], 
 69               tv * A[2], 
 70               tv * A[3]] 
 71   
 72          super(L1TV, self).__init__(l1, A=A, mu=mu, penalty_start=penalty_start) 
 73   
 74          self.reset() 
 75   
 76 -    def reset(self): 
 77   
 78          self._lambda_max = None 
 79   
 80 -    def f(self, beta): 
 81          """ Function value. 
 82          """ 
 83          if self.l < consts.TOLERANCE and self.g < consts.TOLERANCE: 
 84              return 0.0 
 85   
 86          if self.penalty_start > 0: 
 87              beta_ = beta[self.penalty_start:, :] 
 88          else: 
 89              beta_ = beta 
 90   
 91          # lambda and gamma are in A. 
 92          A = self.A() 
 93          return maths.norm1(A[0].dot(beta_)) + \ 
 94                 np.sum(np.sqrt(A[1].dot(beta_) ** 2.0 + 
 95                                A[2].dot(beta_) ** 2.0 + 
 96                                A[3].dot(beta_) ** 2.0)) 
 97   
 98 -    def fmu(self, beta, mu=None): 
 99          """Returns the smoothed function value. 
100   
101          Parameters: 
102          ---------- 
103          beta : Numpy array. The weight vector. 
104   
105          mu : Non-negative float. The regularisation constant for the smoothing. 
106          """ 
107          if mu is None: 
108              mu = self.get_mu() 
109   
110          alpha = self.alpha(beta) 
111          alpha_sqsum = 0.0 
112          for a in alpha: 
113              alpha_sqsum += np.sum(a ** 2.0) 
114   
115          Aa = self.Aa(alpha) 
116   
117          return np.dot(beta.T, Aa)[0, 0] - (mu / 2.0) * alpha_sqsum 
118   
119 -    def phi(self, alpha, beta): 
120          """ Function value with known alpha. 
121   
122          From the interface "NesterovFunction". 
123          """ 
124          if self.l < consts.TOLERANCE and self.g < consts.TOLERANCE: 
125              return 0.0 
126   
127          Aa = self.Aa(alpha) 
128   
129          alpha_sqsum = 0.0 
130          for a in alpha: 
131              alpha_sqsum += np.sum(a ** 2.0) 
132   
133          if self.penalty_start > 0: 
134              beta_ = beta[self.penalty_start:, :] 
135          else: 
136              beta_ = beta 
137   
138          return np.dot(beta_.T, Aa)[0, 0] - (self.mu / 2.0) * alpha_sqsum 
139   
140 -    def lambda_max(self): 
141          """ Largest eigenvalue of the corresponding covariance matrix. 
142   
143          From the interface "Eigenvalues". 
144          """ 
145          # Note that we can save the state here since lmax(A) does not change. 
146          if len(self._A) == 4 and self._A[2].nnz == 0 and self._A[3].nnz == 0: 
147  #        if len(self._shape) == 3 \ 
148  #            and self._shape[0] == 1 and self._shape[1] == 1: 
149              # TODO: Instead of p, this should really be the number of non-zero 
150              # rows of A. 
151              p = self._A[1].shape[0] 
152              lmaxTV = 2.0 * (1.0 - math.cos(float(p - 1) * math.pi 
153                                             / float(p))) 
154              self._lambda_max = lmaxTV * self.g ** 2.0 + self.l ** 2.0 
155   
156          elif self._lambda_max is None: 
157   
158              from parsimony.algorithms.nipals import FastSparseSVD 
159   
160              A = sparse.vstack(self.A()[1:]) 
161              # TODO: Add max_iter here!! 
162              v = FastSparseSVD().run(A)  # , max_iter=max_iter) 
163              us = A.dot(v) 
164              self._lambda_max = np.sum(us ** 2.0) + self.l ** 2.0 
165   
166          return self._lambda_max 
167   
168  #    """ Linear operator of the Nesterov function. 
169  # 
170  #    From the interface "NesterovFunction". 
171  #    """ 
172  #    def A(self): 
173  # 
174  #        return self._A 
175   
176  #    """ Computes A^\T\alpha. 
177  # 
178  #    From the interface "NesterovFunction". 
179  #    """ 
180  #    def Aa(self, alpha): 
181  # 
182  #        A = self.A() 
183  #        Aa = A[0].T.dot(alpha[0]) 
184  #        for i in xrange(1, len(A)): 
185  #            Aa += A[i].T.dot(alpha[i]) 
186  # 
187  #        return Aa 
188   
189 -    def lA(self): 
190          """ Linear operator of the Nesterov function multiplied by the 
191          corresponding Lagrange multipliers. 
192   
193          Note that in this case, the A matrices are already multiplied by the 
194          Lagrange multipliers. 
195          """ 
196          A = self.A() 
197   
198          return A 
199   
200 -    def alpha(self, beta): 
201          """ Dual variable of the Nesterov function. 
202   
203          From the interface "NesterovFunction". Overloaded since we need to do 
204          more than the default. 
205          """ 
206          A = self.A() 
207   
208          if self.penalty_start > 0: 
209              beta_ = beta[self.penalty_start:, :] 
210          else: 
211              beta_ = beta 
212   
213          a = [0] * len(A) 
214          a[0] = (1.0 / self.mu) * A[0].dot(beta_) 
215          a[1] = (1.0 / self.mu) * A[1].dot(beta_) 
216          a[2] = (1.0 / self.mu) * A[2].dot(beta_) 
217          a[3] = (1.0 / self.mu) * A[3].dot(beta_) 
218          # Remember: lambda and gamma are already in the A matrices. 
219   
220          return self.project(a) 
221   
222 -    def project(self, a): 
223          """ Projection onto the compact space of the Nesterov function. 
224   
225          From the interface "NesterovFunction". 
226          """ 
227          # L1 
228          al1 = a[0] 
229          anorm_l1 = np.abs(al1) 
230          i_l1 = anorm_l1 > 1.0 
231          anorm_l1_i = anorm_l1[i_l1] 
232          al1[i_l1] = np.divide(al1[i_l1], anorm_l1_i) 
233   
234          # TV 
235          ax = a[1] 
236          ay = a[2] 
237          az = a[3] 
238          anorm_tv = ax ** 2.0 + ay ** 2.0 + az ** 2.0 
239          i_tv = anorm_tv > 1.0 
240   
241          anorm_tv_i = anorm_tv[i_tv] ** 0.5  # Square root taken here. Faster. 
242          ax[i_tv] = np.divide(ax[i_tv], anorm_tv_i) 
243          ay[i_tv] = np.divide(ay[i_tv], anorm_tv_i) 
244          az[i_tv] = np.divide(az[i_tv], anorm_tv_i) 
245   
246          return [al1, ax, ay, az] 
247   
248 -    def estimate_mu(self, beta): 
249          """Computes a "good" value of mu with respect to the given beta. 
250   
251          From the interface "NesterovFunction". 
252          """ 
253          raise NotImplementedError("We do not use this here!") 
254   
255 -    def M(self): 
256          """ The maximum value of the regularisation of the dual variable. We 
257          have 
258   
259              M = max_{alpha in K} 0.5*|alpha|²_2. 
260   
261          From the interface "NesterovFunction". 
262          """ 
263          A = self.A() 
264   
265          # A[0] is L1, A[1-3] is TV. 
266          return (A[0].shape[0] / 2.0) \ 
267               + (A[1].shape[0] / 2.0) 
268   
269   
270  @utils.deprecated("linear_operator_from_mask") 
271 -def A_from_mask(*args, **kwargs): 
272   
273      return linear_operator_from_mask(*args, **kwargs) 
274   
275   
276  # TODO: Do we need to take the number of variables here? 
277  # Why not use np.prod(shape) + penalty_start instead and save a parameter? 
278 -def linear_operator_from_mask(mask, num_variables, penalty_start=0): 
279      """Generates the linear operator for the total variation Nesterov function 
280      from a mask for a 3D image. 
281   
282      Parameters 
283      ---------- 
284      mask : Numpy array. The mask. The mask does not involve any intercept 
285              variables. 
286   
287      num_variables : Positive integer. The total number of variables, including 
288              the intercept variable(s). 
289   
290      penalty_start : Non-negative integer. The number of variables to exempt 
291              from penalisation. Equivalently, the first index to be penalised. 
292              Default is 0, all variables are included. 
293      """ 
294      Atv, _ = tv.A_from_mask(mask) 
295      Al1 = l1.A_from_variables(num_variables, penalty_start=penalty_start) 
296   
297      return Al1[0], Atv[0], Atv[1], Atv[2] 
298   
299   
300  @utils.deprecated("linear_operator_from_shape") 
301 -def A_from_shape(*args, **kwargs): 
302   
303      return linear_operator_from_shape(*args, **kwargs) 
304   
305   
306  # TODO: Do we need to take the number of variables here? 
307  # Why not use np.prod(shape) + penalty_start instead and save a parameter? 
308 -def linear_operator_from_shape(shape, num_variables, penalty_start=0): 
309      """Generates the linear operator for the total variation Nesterov function 
310      from the shape of a 3D image. 
311   
312      Parameters 
313      ---------- 
314      shape : List or tuple with 1, 2 or 3 elements. The shape of the 1D, 2D or 
315              3D image. shape has the form (Z, Y, X), where Z is the number of 
316              "layers", Y is the number of rows and X is the number of columns. 
317              The shape does not involve any intercept variables. 
318   
319      num_variables : Positive integer. The total number of variables, including 
320              the intercept variable(s). 
321   
322      penalty_start : Non-negative integer. The number of variables to exempt 
323              from penalisation. Equivalently, the first index to be penalised. 
324              Default is 0, all variables are included. 
325      """ 
326      Atv, _ = tv.A_from_shape(shape) 
327      Al1 = l1.A_from_variables(num_variables, penalty_start=penalty_start) 
328   
329      return Al1[0], Atv[0], Atv[1], Atv[2] 
330