Source Code for Module minimise.trust_region

 1  from Numeric import dot, sqrt 
 2   
 3   
 4 -def trust_region(delta_max=1e5, delta0=1.0, eta=0.2): 
 5          """An algorithm for trust region radius selection. 
 6   
 7          Page 68 from 'Numerical Optimization' by Jorge Nocedal and Stephen J. Wright, 1999 
 8          """ 
 9   
10          while 1: 
11                  # Optain pk. 
12                  pk = get_pk() 
13   
14                  # Evaluate rho. 
15                  rho = calc_rho(func, xk, pk, dfk, Bk) 
16   
17                  # Calculate the Euclidean norm of pk. 
18                  norm_pk = sqrt(dot(pk, pk)) 
19   
20                  # Choose the trust region radius for the next iteration. 
21                  # Rho is close to zero or negative, therefore the trust region is shrunk. 
22                  if rho < 0.25: 
23                          delta_new = 0.25 * norm_pk 
24                  # Rho is close to one and pk has reached the boundary of the trust region, therefore the trust region is expanded. 
25                  elif rho > 0.75 and norm_pk == delta: 
26                          delta_new = min(2.0*delta, delta_max) 
27                  # Rho is positive but not close to one, therefore the trust region is unaltered. 
28                  else: 
29                          delta_new = delta 
30   
31                  # Choose the position for the next iteration. 
32                  if rho > eta: 
33                          xk_new = xk + pk 
34                  else: 
35                          xk_new = xk 
36   
37                  # Update. 
38                  delta = delta_new 
39                  xk = xk_new 
40   
41   
42 -def calc_rho(func, xk, pk, dfk, Bk): 
43          """Function to calculate the ratio rho used to choose the trust region radius. 
44   
45          The ratio is defined as: 
46   
47                          f(xk) - f(xk + pk) 
48                  rho  =  ------------------ 
49                            mk(0) - mk(pk) 
50   
51          Where the numerator is called the actual reduction and the denominator is the predicted reduction. 
52          """ 
53   
54          # Actual reduction. 
55          fk = apply(func, (xk,)+args) 
56          f_pk = apply(func, (xk + pk,)+args) 
57          act_red = f - f_p 
58   
59          # Predicted reduction. 
60          mk_pk = fk + dot(dfk, pk) + 0.5 * dot(pk, dot(Bk, pk)) 
61          pred_red = f - mk_pk 
62   
63          # Rho. 
64          rho = act_red / pred_red 
65          return rho 
66