parsimony.functions.nesterov.grouptv

1 # -*- coding: utf-8 -*- 2 """ 3 The :mod:`parsimony.functions.nesterov.grouptv` module contains the loss 4 function and helper functions for group Total variation, Group TV, smoothed 5 using Nesterov's smoothing technique. 6 7 Created on Mon May 5 11:46:45 2014 8 9 Copyright (c) 2013-2014, CEA/DSV/I2BM/Neurospin. All rights reserved. 10 11 @author: Tommy Löfstedt 12 @email: lofstedt.tommy@gmail.com 13 @license: BSD 3-clause. 14 """ 15 import scipy.sparse as sparse 16 import numpy as np 17 18 #from .properties import NesterovFunction 19 from .. import properties 20 import parsimony.utils.consts as consts 21 import parsimony.utils as utils 22 23 __all__ = ["GroupTotalVariation", 24 "linear_operator_from_masks", "A_from_masks", 25 "linear_operator_from_rects", "A_from_rects"]

26 27 28 -class GroupTotalVariation(properties.NesterovFunction, 29 properties.Penalty, 30 properties.Constraint):

31 """The smoothed Group total variation (Group TV) function 32 33 f(beta) = l * (GroupTV(beta) - c), 34 35 where GroupTV(beta) is the smoothed group total variation function. The 36 constrained version has the form 37 38 GroupTV(beta) <= c. 39 """

40 - def __init__(self, l, c=0.0, A=None, mu=0.0, penalty_start=0):

41 """ 42 Parameters 43 ---------- 44 l : Non-negative float. The Lagrange multiplier, or regularisation 45 constant, of the function. 46 47 c : Float. The limit of the constraint. The function is feasible if 48 f(beta) <= c. The default value is c=0, i.e. the default is a 49 regularised formulation. 50 51 A : Numpy array (usually sparse). The linear operator for the Nesterov 52 formulation. Will have length 3 * number of groups, and the 53 group A matrices are assumed to be next to eachother in the 54 list. A may not be None! 55 56 mu : Non-negative float. The regularisation constant for the smoothing. 57 58 penalty_start : Non-negative integer. The number of columns, variables 59 etc., to exempt from penalisation. Equivalently, the first 60 index to be penalised. Default is 0, all columns are included. 61 """ 62 super(GroupTotalVariation, self).__init__(l, A=A, mu=mu, 63 penalty_start=penalty_start) 64 self._p = A[0].shape[1] 65 66 self.c = float(c) 67 68 self.reset()

69

70 - def reset(self):

71 72 self._lambda_max = None

73

74 - def f(self, beta):

75 """ Function value. 76 """ 77 if self.l < consts.TOLERANCE: 78 return 0.0 79 80 if self.penalty_start > 0: 81 beta_ = beta[self.penalty_start:, :] 82 else: 83 beta_ = beta 84 85 A = self.A() 86 87 f = 0.0 88 for g in xrange(0, len(A), 3): 89 f += np.sum(np.sqrt(A[g + 0].dot(beta_) ** 2.0 + 90 A[g + 1].dot(beta_) ** 2.0 + 91 A[g + 2].dot(beta_) ** 2.0)) 92 93 return self.l * (f - self.c)

94

95 - def phi(self, alpha, beta):

96 """Function value with known alpha. 97 98 From the interface "NesterovFunction". 99 """ 100 if self.l < consts.TOLERANCE: 101 return 0.0 102 103 if self.penalty_start > 0: 104 beta_ = beta[self.penalty_start:, :] 105 else: 106 beta_ = beta 107 108 Aa = self.Aa(alpha) 109 110 alpha_sqsum = 0.0 111 for a in alpha: 112 alpha_sqsum += np.sum(a ** 2.0) 113 114 mu = self.get_mu() 115 116 return self.l * ((np.dot(beta_.T, Aa)[0, 0] 117 - (mu / 2.0) * alpha_sqsum) - self.c)

118

119 - def feasible(self, beta):

120 """Feasibility of the constraint. 121 122 From the interface "Constraint". 123 """ 124 if self.penalty_start > 0: 125 beta_ = beta[self.penalty_start:, :] 126 else: 127 beta_ = beta 128 129 A = self.A() 130 131 f = 0.0 132 for g in xrange(0, len(A), 3): 133 f += np.sum(np.sqrt(A[g + 0].dot(beta_) ** 2.0 + 134 A[g + 1].dot(beta_) ** 2.0 + 135 A[g + 2].dot(beta_) ** 2.0)) 136 137 return f <= self.c

138

139 - def L(self):

140 """ Lipschitz constant of the gradient. 141 142 From the interface "LipschitzContinuousGradient". 143 """ 144 if self.l < consts.TOLERANCE: 145 return 0.0 146 147 lmaxA = self.lambda_max() 148 149 return self.l * lmaxA / self.mu

150

151 - def lambda_max(self):

152 """ Largest eigenvalue of the corresponding covariance matrix. 153 154 From the interface "Eigenvalues". 155 """ 156 # Note that we can save the state here since lmax(A) does not change. 157 # TODO: This only work if the elements of self._A are scipy.sparse. We 158 # should allow dense matrices as well. 159 if self._lambda_max is None: 160 161 from parsimony.algorithms.nipals import FastSparseSVD 162 163 A = sparse.vstack(self.A()) 164 # TODO: Add max_iter here! 165 v = FastSparseSVD().run(A) # , max_iter=max_iter) 166 us = A.dot(v) 167 self._lambda_max = np.sum(us ** 2.0) 168 169 return self._lambda_max

170

171 - def project(self, a):

172 """ Projection onto the compact space of the Nesterov function. 173 174 From the interface "NesterovFunction". 175 """ 176 for g in xrange(0, len(a), 3): 177 178 ax = a[g + 0] 179 ay = a[g + 1] 180 az = a[g + 2] 181 anorm = ax ** 2.0 + ay ** 2.0 + az ** 2.0 182 i = anorm > 1.0 183 184 anorm_i = anorm[i] ** 0.5 # Square root is taken here. Faster. 185 ax[i] = np.divide(ax[i], anorm_i) 186 ay[i] = np.divide(ay[i], anorm_i) 187 az[i] = np.divide(az[i], anorm_i) 188 189 a[g + 0] = ax 190 a[g + 1] = ay 191 a[g + 2] = az 192 193 return a

194

195 - def M(self):

196 """The maximum value of the regularisation of the dual variable. We 197 have 198 199 M = max_{alpha in K} 0.5*|alpha|²_2. 200 201 From the interface "NesterovFunction". 202 """ 203 A = self.A() 204 n = 0 205 for g in xrange(0, len(A), 3): 206 n += A[g].nnz # The number of non-zero rows of Ag, the gth group. 207 208 return float(n) / 2.0

209

210 - def estimate_mu(self, beta):

211 """Computes a "good" value of mu with respect to the given beta. 212 213 From the interface "NesterovFunction". 214 """ 215 if self.penalty_start > 0: 216 beta_ = beta[self.penalty_start:, :] 217 else: 218 beta_ = beta 219 220 A = self.A() 221 max_norm = 0.0 222 for g in xrange(0, len(A), 3): 223 224 ax = A[g + 0].dot(beta_) 225 ay = A[g + 1].dot(beta_) 226 az = A[g + 2].dot(beta_) 227 228 anorm = (ax ** 2.0 + ay ** 2.0 + az ** 2.0) ** 0.5 # Compute norm. 229 230 max_norm = max(max_norm, np.max(anorm)) # The overall maximum. 231 232 return max_norm

233

234 235 @utils.deprecated("linear_operator_from_masks") 236 -def A_from_masks(*args, **kwargs):

237 238 return linear_operator_from_masks(*args, **kwargs)

239

240 241 -def linear_operator_from_masks(masks, weights=None):

242 """Generates the linear operator for the group total variation Nesterov 243 function from a mask for a 3D image. 244 245 Parameters 246 ---------- 247 masks : List of numpy arrays. The mask for each group. Each mask is an 248 integer (0 or 1) or boolean numpy array or the same shape as the 249 actual data. The mask does not involve any intercept variables. 250 251 weights : List of floats. The weights account for different group sizes, 252 or incorporates some prior knowledge about the importance of the 253 groups. Default value is the square roots of the group sizes. 254 """ 255 import parsimony.functions.nesterov.tv as tv 256 257 if isinstance(masks, tuple): 258 masks = list(masks) 259 260 A = [] 261 262 G = len(masks) 263 for g in xrange(G): 264 mask = masks[g] 265 266 if weights is None: 267 weight = np.sqrt(np.sum(mask)) 268 else: 269 weight = weights[g] 270 271 # Compute group A matrix 272 Ag, _ = tv.A_from_subset_mask(mask) 273 274 # Include the weights 275 if weight != 1.0 and weight != 1: 276 for A_ in Ag: 277 A_ *= weight 278 279 A += Ag 280 281 return A

282

283 284 @utils.deprecated("linear_operator_from_rects") 285 -def A_from_rects(*args, **kwargs):

286 287 return linear_operator_from_rects(*args, **kwargs)

288

289 290 -def linear_operator_from_rects(rects, shape, weights=None):

291 """Generates the linear operator for the group total variation Nesterov 292 function from the rectange of a 3D image. 293 294 Parameters 295 ---------- 296 rects : List of lists or tuples with 2-, 4- or 6-tuple elements. The shape 297 of the patch of the 1D, 2D or 3D image to smooth. The elements of 298 rects has the form ((x1, x2),), ((y1, y2), (x1, x2)) or ((z1, z2), 299 (y1, y2), (x1, x2)), where z is the "layers", y rows and x is the 300 columns and x1 means the first column to include, x2 is one beyond 301 the last column to include, and similarly for y and z. The rect 302 does not involve any intercept variables. 303 304 shape : List or tuple with 1, 2 or 3 integers. The shape of the 1D, 2D or 305 3D image. shape has the form (X,), (Y, X) or (Z, Y, X), where Z is 306 the number of "layers", Y is the number of rows and X is the number 307 of columns. The shape does not involve any intercept variables. 308 309 weights : List of floats. The weights account for different group sizes, 310 or incorporates some prior knowledge about the importance of the 311 groups. Default value is the square roots of the group sizes. 312 """ 313 import parsimony.functions.nesterov.tv as tv 314 315 A = [] 316 G = len(rects) 317 for g in xrange(G): 318 rect = rects[g] 319 if len(rect) == 1: 320 rect = [(0, 1), (0, 1), rect[0]] 321 elif len(rect) == 2: 322 rect = [(0, 1), rect[0], rect[1]] 323 324 while len(shape) < 3: 325 shape = tuple([1] + list(shape)) 326 327 mask = np.zeros(shape, dtype=bool) 328 z1 = rect[0][0] 329 z2 = rect[0][1] 330 y1 = rect[1][0] 331 y2 = rect[1][1] 332 x1 = rect[2][0] 333 x2 = rect[2][1] 334 mask[z1:z2, y1:y2, x1:x2] = True 335 336 if weights is None: 337 weight = np.sqrt(np.sum(mask)) 338 else: 339 weight = weights[g] 340 341 # Compute group A matrix 342 Ag, _ = tv.A_from_subset_mask(mask) 343 344 # Include the weights 345 if weight != 1.0 and weight != 1: 346 for A_ in Ag: 347 A_ *= weight 348 349 A += Ag 350 351 return A

352

Source Code for Module parsimony.functions.nesterov.grouptv