Source Code for Module mvpa.clfs.model_selector

  1  # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- 
  2  # vi: set ft=python sts=4 ts=4 sw=4 et: 
  3  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 
  4  # 
  5  #   Copyright (c) 2008 Emanuele Olivetti <emanuele@relativita.com> 
  6  #   See COPYING file distributed along with the PyMVPA package for the 
  7  #   copyright and license terms. 
  8  # 
  9  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 
 10  """Model selction.""" 
 11   
 12  __docformat__ = 'restructuredtext' 
 13   
 14   
 15  import numpy as N 
 16  from mvpa.base import externals 
 17  from mvpa.misc.exceptions import InvalidHyperparameterError 
 18   
 19  if externals.exists("scipy", raiseException=True): 
 20      import scipy.linalg as SL 
 21   
 22  # no sense to import this module if openopt is not available 
 23  if externals.exists("openopt", raiseException=True): 
 24      try: 
 25          from openopt import NLP 
 26      except ImportError: 
 27          from scikits.openopt import NLP 
 28   
 29  if __debug__: 
 30      from mvpa.base import debug 
 31   
 32 -def _openopt_debug(): 
 33      # shut up or make verbose OpenOpt 
 34      # (-1 = no logs, 0 = small log, 1 = verbose) 
 35      if __debug__: 
 36          da = debug.active 
 37          if 'OPENOPT' in da: 
 38              return 1 
 39          elif 'MOD_SEL' in da: 
 40              return 0 
 41      return -1 
 42   
 43   
 44 -class ModelSelector(object): 
 45      """Model selection facility. 
 46   
 47      Select a model among multiple models (i.e., a parametric model, 
 48      parametrized by a set of hyperparamenters). 
 49      """ 
 50   
 51 -    def __init__(self, parametric_model, dataset): 
 52          self.parametric_model = parametric_model 
 53          self.dataset = dataset 
 54          self.hyperparameters_best = None 
 55          self.log_marginal_likelihood_best = None 
 56          self.problem = None 
 57          pass 
 58   
 59 -    def max_log_marginal_likelihood(self, hyp_initial_guess, maxiter=1, 
 60              optimization_algorithm="scipy_cg", ftol=1.0e-3, fixedHypers=None, 
 61              use_gradient=False, logscale=False): 
 62          """ 
 63          Set up the optimization problem in order to maximize 
 64          the log_marginal_likelihood. 
 65   
 66          :Parameters: 
 67   
 68            parametric_model : Classifier 
 69              the actual parameteric model to be optimized. 
 70   
 71            hyp_initial_guess : numpy.ndarray 
 72              set of hyperparameters' initial values where to start 
 73              optimization. 
 74   
 75            optimization_algorithm : string 
 76              actual name of the optimization algorithm. See 
 77              http://scipy.org/scipy/scikits/wiki/NLP 
 78              for a comprehensive/updated list of available NLP solvers. 
 79              (Defaults to 'ralg') 
 80   
 81            ftol : float 
 82              threshold for the stopping criterion of the solver, 
 83              which is mapped in OpenOpt NLP.ftol 
 84              (Defaults to 1.0e-3) 
 85   
 86            fixedHypers : numpy.ndarray (boolean array) 
 87              boolean vector of the same size of hyp_initial_guess; 
 88              'False' means that the corresponding hyperparameter must 
 89              be kept fixed (so not optimized). 
 90              (Defaults to None, which during means all True) 
 91   
 92          NOTE: the maximization of log_marginal_likelihood is a non-linear 
 93          optimization problem (NLP). This fact is confirmed by Dmitrey, 
 94          author of OpenOpt. 
 95          """ 
 96          self.problem = None 
 97          self.use_gradient = use_gradient 
 98          self.logscale = logscale # use log-scale on hyperparameters to enhance numerical stability 
 99          self.optimization_algorithm = optimization_algorithm 
100          self.hyp_initial_guess = N.array(hyp_initial_guess) 
101          self.hyp_initial_guess_log = N.log(self.hyp_initial_guess) 
102          if fixedHypers is None: 
103              fixedHypers = N.zeros(self.hyp_initial_guess.shape[0],dtype=bool) 
104              pass 
105          self.freeHypers = -fixedHypers 
106          if self.logscale: 
107              self.hyp_running_guess = self.hyp_initial_guess_log.copy() 
108          else: 
109              self.hyp_running_guess = self.hyp_initial_guess.copy() 
110              pass 
111          self.f_last_x = None 
112   
113          def f(x): 
114              """ 
115              Wrapper to the log_marginal_likelihood to be 
116              maximized. 
117              """ 
118              # XXX EO: since some OpenOpt NLP solvers does not 
119              # implement lower bounds the hyperparameters bounds are 
120              # implemented inside PyMVPA: (see dmitrey's post on 
121              # [SciPy-user] 20080628). 
122              # 
123              # XXX EO: OpenOpt does not implement optimization of a 
124              # subset of the hyperparameters so it is implemented here. 
125              # 
126              # XXX EO: OpenOpt does not implement logrithmic scale of 
127              # the hyperparameters (to enhance numerical stability), so 
128              # it is implemented here. 
129              self.f_last_x = x.copy() 
130              self.hyp_running_guess[self.freeHypers] = x 
131              # REMOVE print "guess:",self.hyp_running_guess,x 
132              try: 
133                  if self.logscale: 
134                      self.parametric_model.set_hyperparameters(N.exp(self.hyp_running_guess)) 
135                  else: 
136                      self.parametric_model.set_hyperparameters(self.hyp_running_guess) 
137                      pass 
138              except InvalidHyperparameterError: 
139                  if __debug__: debug("MOD_SEL","WARNING: invalid hyperparameters!") 
140                  return -N.inf 
141              try: 
142                  self.parametric_model.train(self.dataset) 
143              except (N.linalg.linalg.LinAlgError, SL.basic.LinAlgError, ValueError): 
144                  # Note that ValueError could be raised when Cholesky gets Inf or Nan. 
145                  if __debug__: debug("MOD_SEL", "WARNING: Cholesky failed! Invalid hyperparameters!") 
146                  return -N.inf 
147              log_marginal_likelihood = self.parametric_model.compute_log_marginal_likelihood() 
148              # REMOVE print log_marginal_likelihood 
149              return log_marginal_likelihood 
150   
151          def df(x): 
152              """ 
153              Proxy to the log_marginal_likelihood first 
154              derivative. Necessary for OpenOpt when using derivatives. 
155              """ 
156              self.hyp_running_guess[self.freeHypers] = x 
157              # REMOVE print "df guess:",self.hyp_running_guess,x 
158              # XXX EO: Most of the following lines can be skipped if 
159              # df() is computed just after f() with the same 
160              # hyperparameters. The partial results obtained during f() 
161              # are what is needed for df(). For now, in order to avoid 
162              # bugs difficult to trace, we keep this redunundancy. A 
163              # deep check with how OpenOpt works or using memoization 
164              # should solve this issue. 
165              try: 
166                  if self.logscale: 
167                      self.parametric_model.set_hyperparameters(N.exp(self.hyp_running_guess)) 
168                  else: 
169                      self.parametric_model.set_hyperparameters(self.hyp_running_guess) 
170                      pass 
171              except InvalidHyperparameterError: 
172                  if __debug__: debug("MOD_SEL", "WARNING: invalid hyperparameters!") 
173                  return -N.inf 
174              # Check if it is possible to avoid useless computations 
175              # already done in f(). According to tests and information 
176              # collected from OpenOpt people, it is sufficiently 
177              # unexpected that the following test succeed: 
178              if N.any(x!=self.f_last_x): 
179                  if __debug__: debug("MOD_SEL","UNEXPECTED: recomputing train+log_marginal_likelihood.") 
180                  try: 
181                      self.parametric_model.train(self.dataset) 
182                  except (N.linalg.linalg.LinAlgError, SL.basic.LinAlgError, ValueError): 
183                      if __debug__: debug("MOD_SEL", "WARNING: Cholesky failed! Invalid hyperparameters!") 
184                      # XXX EO: which value for the gradient to return to 
185                      # OpenOpt when hyperparameters are wrong? 
186                      return N.zeros(x.size) 
187                  log_marginal_likelihood = self.parametric_model.compute_log_marginal_likelihood() # recompute what's needed (to be safe) REMOVE IN FUTURE! 
188                  pass 
189              if self.logscale: 
190                  gradient_log_marginal_likelihood = self.parametric_model.compute_gradient_log_marginal_likelihood_logscale() 
191              else: 
192                  gradient_log_marginal_likelihood = self.parametric_model.compute_gradient_log_marginal_likelihood() 
193                  pass 
194              # REMOVE print "grad:",gradient_log_marginal_likelihood 
195              return gradient_log_marginal_likelihood[self.freeHypers] 
196   
197   
198          if self.logscale: 
199              # vector of hyperparameters' values where to start the search 
200              x0 = self.hyp_initial_guess_log[self.freeHypers] 
201          else: 
202              x0 = self.hyp_initial_guess[self.freeHypers] 
203              pass 
204          self.contol = 1.0e-20 # Constraint tolerance level 
205          # XXX EO: is it necessary to use contol when self.logscale is 
206          # True and there is no lb? Ask dmitrey. 
207          if self.use_gradient: 
208              # actual instance of the OpenOpt non-linear problem 
209              self.problem = NLP(f, x0, df=df, contol=self.contol, goal='maximum') 
210          else: 
211              self.problem = NLP(f, x0, contol=self.contol, goal='maximum') 
212              pass 
213          self.problem.name = "Max LogMargLikelihood" 
214          if not self.logscale: 
215               # set lower bound for hyperparameters: avoid negative 
216               # hyperparameters. Note: problem.n is the size of 
217               # hyperparameters' vector 
218              self.problem.lb = N.zeros(self.problem.n)+self.contol 
219              pass 
220          # max number of iterations for the optimizer. 
221          self.problem.maxiter = maxiter 
222          # check whether the derivative of log_marginal_likelihood converged to 
223          # zero before ending optimization 
224          self.problem.checkdf = True 
225           # set increment of log_marginal_likelihood under which the optimizer stops 
226          self.problem.ftol = ftol 
227          self.problem.iprint = _openopt_debug() 
228          return self.problem 
229   
230   
231 -    def solve(self, problem=None): 
232          """Solve the maximization problem, check outcome and collect results. 
233          """ 
234          # XXX: this method can be made more abstract in future in the 
235          # sense that it could work not only for 
236          # log_marginal_likelihood but other measures as well 
237          # (e.g. cross-valideted error). 
238   
239          if N.all(self.freeHypers==False): # no optimization needed 
240              self.hyperparameters_best = self.hyp_initial_guess.copy() 
241              try: 
242                  self.parametric_model.set_hyperparameters(self.hyperparameters_best) 
243              except InvalidHyperparameterError: 
244                  if __debug__: debug("MOD_SEL", "WARNING: invalid hyperparameters!") 
245                  self.log_marginal_likelihood_best = -N.inf 
246                  return self.log_marginal_likelihood_best 
247              self.parametric_model.train(self.dataset) 
248              self.log_marginal_likelihood_best = self.parametric_model.compute_log_marginal_likelihood() 
249              return self.log_marginal_likelihood_best 
250   
251          result = self.problem.solve(self.optimization_algorithm) # perform optimization! 
252          if result.stopcase == -1: 
253              # XXX: should we use debug() for the following messages? 
254              # If so, how can we track the missing convergence to a 
255              # solution? 
256              print "Unable to find a maximum to log_marginal_likelihood" 
257          elif result.stopcase == 0: 
258              print "Limits exceeded" 
259          elif result.stopcase == 1: 
260              self.hyperparameters_best = self.hyp_initial_guess.copy() 
261              if self.logscale: 
262                  self.hyperparameters_best[self.freeHypers] = N.exp(result.xf) # best hyperparameters found # NOTE is it better to return a copy? 
263              else: 
264                  self.hyperparameters_best[self.freeHypers] = result.xf 
265                  pass 
266              self.log_marginal_likelihood_best = result.ff # actual best vuale of log_marginal_likelihood 
267              pass 
268          self.stopcase = result.stopcase 
269          return self.log_marginal_likelihood_best 
270   
271      pass 
272