Package minfx :: Package hessian_mods :: Module gmw81
[hide private]
[frames] | no frames]

Source Code for Module minfx.hessian_mods.gmw81

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2003-2013 Edward d'Auvergne                                   # 
  4  #                                                                             # 
  5  # This file is part of the minfx optimisation library,                        # 
  6  # https://sourceforge.net/projects/minfx                                      # 
  7  #                                                                             # 
  8  # This program is free software: you can redistribute it and/or modify        # 
  9  # it under the terms of the GNU General Public License as published by        # 
 10  # the Free Software Foundation, either version 3 of the License, or           # 
 11  # (at your option) any later version.                                         # 
 12  #                                                                             # 
 13  # This program is distributed in the hope that it will be useful,             # 
 14  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 15  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 16  # GNU General Public License for more details.                                # 
 17  #                                                                             # 
 18  # You should have received a copy of the GNU General Public License           # 
 19  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 20  #                                                                             # 
 21  ############################################################################### 
 22   
 23  # Module docstring. 
 24  """The Gill, Murray, and Wright (GMW) modified Cholesky Hessian modification algorithm. 
 25   
 26  This file is part of the U{minfx optimisation library<https://sourceforge.net/projects/minfx>}. 
 27  """ 
 28   
 29  # Python module imports. 
 30  from math import sqrt 
 31  from numpy import dot, transpose 
 32  from numpy.linalg import cholesky, eig, inv, LinAlgError, solve 
 33   
 34   
35 -def gmw(dfk, d2fk, I, n, mach_acc, print_prefix, print_flag, return_matrix=0):
36 """The Gill, Murray, and Wright modified Cholesky algorithm. 37 38 Algorithm 6.5 from page 148 of 'Numerical Optimization' by Jorge Nocedal and Stephen J. Wright, 1999, 2nd ed. 39 40 Returns the modified Newton step. 41 """ 42 43 # Test matrix. 44 #d2fk = array([[4, 2, 1], [2, 6, 3], [1, 3, -0.004]], float64) 45 #n = len(d2fk) 46 #I = identity(n, float64) 47 48 # Calculate gamma(A) and xi(A). 49 gamma = 0.0 50 xi = 0.0 51 for i in range(n): 52 gamma = max(abs(d2fk[i, i]), gamma) 53 for j in range(i+1, n): 54 xi = max(abs(d2fk[i, j]), xi) 55 56 # Calculate delta and beta. 57 delta = mach_acc * max(gamma + xi, 1.0) 58 if n == 1: 59 beta = sqrt(max(gamma, mach_acc)) 60 else: 61 beta = sqrt(max(gamma, xi / sqrt(n**2 - 1.0), mach_acc)) 62 63 # Initialise data structures. 64 a = 1.0 * d2fk 65 r = 0.0 * d2fk 66 e = 0.0 * dfk 67 P = 1.0 * I 68 69 # Debugging. 70 if print_flag >= 3: 71 old_eigen = eig(d2fk) 72 print(print_prefix + "dfk: " + repr(dfk)) 73 print(print_prefix + "d2fk:\n" + repr(d2fk)) 74 75 # Main loop. 76 for j in range(n): 77 # Row and column swapping, find the index > j of the largest diagonal element. 78 q = j 79 for i in range(j+1, n): 80 if abs(a[i, i]) >= abs(a[q, q]): 81 q = i 82 83 # Interchange row and column j and q (if j != q). 84 if q != j: 85 # Temporary permutation matrix for swaping 2 rows or columns. 86 p = 1.0 * I 87 88 # Modify the permutation matrix P by swaping columns. 89 row_P = 1.0*P[:, q] 90 P[:, q] = P[:, j] 91 P[:, j] = row_P 92 93 # Modify the permutation matrix p by swaping rows (same as columns because p = pT). 94 row_p = 1.0*p[q] 95 p[q] = p[j] 96 p[j] = row_p 97 98 # Permute a and r (p = pT). 99 a = dot(p, dot(a, p)) 100 r = dot(r, p) 101 102 # Calculate dj. 103 theta_j = 0.0 104 if j < n-1: 105 for i in range(j+1, n): 106 theta_j = max(theta_j, abs(a[j, i])) 107 dj = max(abs(a[j, j]), (theta_j/beta)**2, delta) 108 109 # Calculate e (not really needed!). 110 if print_flag >= 3: 111 e[j] = dj - a[j, j] 112 113 # Calculate row j of r and update a. 114 r[j, j] = sqrt(dj) # Damned sqrt introduces roundoff error. 115 for i in range(j+1, n): 116 r[j, i] = a[j, i] / r[j, j] 117 for k in range(j+1, i+1): 118 a[i, k] = a[k, i] = a[k, i] - r[j, i] * r[j, k] # Keep matrix a symmetric. 119 120 # The Cholesky factor of d2fk. 121 L = dot(P, transpose(r)) 122 123 # Debugging. 124 if print_flag >= 3: 125 print(print_prefix + "r:\n" + repr(r)) 126 print(print_prefix + "e: " + repr(e)) 127 print(print_prefix + "e rot: " + repr(dot(P, e))) 128 print(print_prefix + "P:\n" + repr(P)) 129 print(print_prefix + "L:\n" + repr(L)) 130 print(print_prefix + "L.LT:\n" + repr(dot(L, transpose(L)))) 131 print(print_prefix + "L.LT - d2fk:\n" + repr(dot(L, transpose(L)) - d2fk)) 132 print(print_prefix + "dot(P, chol(L.LT)):\n" + repr(dot(P, cholesky(dot(L, transpose(L)))))) 133 print(print_prefix + "dot(P, chol(P.L.LT.PT)):\n" + repr(dot(P, cholesky(dot(P, dot(dot(L, transpose(L)), transpose(P))))))) 134 E = 0.0 * d2fk 135 for i in range(n): 136 E[i, i] = e[i] 137 E = dot(P, dot(E, transpose(P))) 138 print(print_prefix + "E:\n" + repr(E)) 139 print(print_prefix + "PT.d2fk+E.P:\n" + repr(dot(transpose(P), dot(d2fk+E, P)))) 140 try: 141 chol = cholesky(dot(transpose(P), dot(d2fk+E, P))) 142 print(print_prefix + "chol(PT.d2fk+E.P):\n" + repr(chol)) 143 print(print_prefix + "rT - chol(PT.d2fk+E.P):\n" + repr(transpose(r) - chol)) 144 chol = dot(P, chol) 145 print(print_prefix + "chol:\n" + repr(chol)) 146 print(print_prefix + "chol reconstructed:\n" + repr(dot(chol, transpose(chol)))) 147 except LinAlgError: 148 print(print_prefix + "Matrix is not positive definite - Cholesky decomposition cannot be computed.") 149 eigen = eig(dot(L, transpose(L))) 150 print(print_prefix + "Old eigenvalues: " + repr(old_eigen)) 151 print(print_prefix + "New eigenvalues: " + repr(eigen)) 152 print(print_prefix + "Newton dir: " + repr(-solve(transpose(L), solve(L, dfk)))) 153 print(print_prefix + "Newton dir using inverse: " + repr(-dot(inv(d2fk+E), dfk))) 154 155 # Calculate the Newton direction. 156 y = solve(L, dfk) 157 if return_matrix: 158 return -solve(transpose(L), y), dot(L, transpose(L)) 159 else: 160 return -solve(transpose(L), y)
161