Source Code for Module linear_operators.iterative.optimization

 1  """ 
 2  Use the optimization scikits to define a set of optimizer and models 
 3  using LinearOperators. 
 4  """ 
 5  import numpy as np 
 6  from scikits.optimization import step, line_search, criterion, optimizer 
 7  from scikits.optimization.optimizer import StandardOptimizer 
 8  import criterions 
 9   
10  # defaults 
11  default_step = step.PRPConjugateGradientStep() 
12  default_line_search = line_search.QuadraticInterpolationSearch(0) 
13  default_stop_criterion = criterion.criterion() 
14   
15  # a class wrapping StandardOptimizer 
16 -class Optimizer(object): 
17 -    def __init__(self, criterion, 
18                   x0=None, 
19                   step=default_step, 
20                   line_search=default_line_search, 
21                   stop_criterion=default_stop_criterion): 
22          self.criterion = criterion 
23          self.n_variables = criterion.n_variables 
24          self.step = step 
25          self.line_search = line_search 
26          self.stop_criterion = stop_criterion 
27          self.first_guess() 
28          self.optimizer = StandardOptimizer(function=self.criterion, 
29                                             step=self.step, 
30                                             line_search=self.line_search, 
31                                             criterion=self.stop_criterion, 
32                                             x0=self.current_solution) 
33          self.optimizer.recordHistory = self.callback 
34 -    def callback(self, **kwargs): 
35          state = self.optimizer.state 
36          iter_ = state['iteration'] 
37          crit = state['new_value'] 
38          step = state['alpha_step'] 
39          if state['iteration'] == 0: 
40              print("Iteration \t Step \t \Criterion") 
41          print("\t%i \t %e \t %e" % (iter_, step, crit)) 
42 -    def first_guess(self, x0=None): 
43          """ 
44          Sets current_solution attribute to initial value. 
45          """ 
46          if x0 is None: 
47              self.current_solution = np.zeros(self.n_variables) 
48          else: 
49              self.current_solution = copy(x0) 
50 -    def iterate(self): 
51          self.optimizer.iterate() 
52 -    def __call__(self): 
53          self.optimizer.optimize() 
54   
55 -class QuadraticOptimizer(Optimizer): 
56 -    def __init__(self, model, data, priors, hypers, **kwargs): 
57          criterion = criterions.QuadraticCriterion(model, data, hypers=hypers, 
58                                                    priors=priors) 
59          kwargs['line_search'] = self.optimal_step 
60          Optimizer.__init__(self, criterion, **kwargs) 
61   
62 -    def optimal_step(self, origin, state): 
63          d = state['direction'] 
64          g = state['gradient'] 
65          H = state['function'].model 
66          Ds = state['function'].priors 
67          hypers = state['function'].hypers 
68          algo_norms = state['function'].norms 
69          # replace norms by Norm2 if not a Norm2 instance 
70          # to handle properly Norm2 with C covariance matrices ... 
71          algo_norms = [n if isinstance(n, norms.Norm2) else norms.Norm2() for n in algo_norms] 
72          # compute quadratic optimal step 
73          a = -.5 * np.dot(d.T, g) 
74          a /= algo_norms[0](H * d) + np.sum([h * n(D * d) for h, D, n in zip(hypers, Ds, algo_norms)]) 
75          return a 
76