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

Source Code for Module minfx.hessian_mods.se99

  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 Schnabel and Eskow (SE) revised modified Cholesky factorisation 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 numpy import array, dot, float64, identity, sort, sqrt, take, transpose, zeros 
 31  from numpy.linalg import cholesky, eig, LinAlgError, solve 
 32   
 33   
34 -def se99(dfk, d2fk, I, n, tau, tau_bar, mu, print_prefix, print_flag, return_matrix=0):
35 """A revised modified cholesky factorisation algorithm. 36 37 Schnabel, R. B. and Eskow, E. 1999, A revised modifed cholesky factorisation algorithm. SIAM J. Optim. 9, 1135-1148. 38 """ 39 40 # Test matrix. 41 d2fk = array([[4, 2, 1], [2, 6, 3], [1, 3, -0.004]], float64) 42 d2fk = array([[ 1890.3, -1750.6, -315.8, 3000.3], 43 [ -1705.6, 1538.3, 284.9, -2706.6], 44 [ -315.8, 284.9, 52.5, -501.2], 45 [ 3000.3, -2706.6, -501.2, 4760.8]], float64) 46 n = len(d2fk) 47 I = identity(n, float64) 48 49 # Create the matrices A and L. 50 A = 1.0 * d2fk 51 L = 0.0 * d2fk 52 53 # Permutation matrix. 54 P = 1.0 * I 55 56 # Calculate gamma. 57 gamma = abs(d2fk[0, 0]) 58 for i in range(n): 59 gamma = max(abs(d2fk[i, i]), gamma) 60 61 # Debugging. 62 E = 0.0 * d2fk 63 if print_flag >= 3: 64 print(print_prefix + "tau: " + repr(tau)) 65 print(print_prefix + "tau_bar: " + repr(tau_bar)) 66 print(print_prefix + "mu: " + repr(mu)) 67 print(print_prefix + "d2fk:\n" + repr(d2fk)) 68 print(print_prefix + "Orig A:\n" + repr(A)) 69 print(print_prefix + "Orig L:\n" + repr(L)) 70 print(print_prefix + "gamma: " + repr(gamma)) 71 72 73 # Phase one, A potentially positive definite. 74 ############################################# 75 76 for j in range(n): 77 if print_flag >= 3: 78 print("\n" + print_prefix + "Iteration j = " + repr(j)) 79 print(print_prefix + "Init A:\n" + repr(A)) 80 print(print_prefix + "Init L:\n" + repr(L)) 81 print(print_prefix + "Phase 1.") 82 83 # Calculate max_Aii and min_Aii 84 max_Aii = A[j, j] 85 min_Aii = A[j, j] 86 for i in range(j, n): 87 max_Aii = max(A[i, i], max_Aii) 88 min_Aii = min(A[i, i], min_Aii) 89 90 # Test for phase 2, A not positive definite. 91 if max_Aii < tau_bar * gamma or min_Aii < - mu * max_Aii: 92 exec_phasetwo(A, L, P, I, E, j, n, tau, tau_bar, gamma, print_prefix, print_flag) 93 break 94 95 else: 96 # Pivot on maximum diagonal of remaining submatrix. 97 i = j 98 for k in range(j, n): 99 if A[k, k] >= A[i, i]: 100 i = k 101 102 # Interchange row and column i and j. 103 p = 1.0 * I 104 if i != j: 105 if print_flag >= 3: 106 print(print_prefix + "Switching rows and columns of " + repr(i) + " and " + repr(j) + " of A and L.") 107 108 # Modify the permutation matrices by swaping rows. 109 row_p, row_P = 1.0*p[i], 1.0*P[i] 110 p[i], P[i] = p[j], P[j] 111 p[j], P[j] = row_p, row_P 112 113 # Permute A. 114 A = dot(p, dot(A, transpose(p))) 115 L = dot(p, dot(L, transpose(p))) 116 117 if print_flag >= 3: 118 print(print_prefix + "Permuted A:\n" + repr(A)) 119 print(print_prefix + "Permuted L:\n" + repr(L)) 120 121 # Test for phase 2 again. 122 min_num = 1e99 123 for i in range(j+1, n): 124 min_num = min(min_num, A[i, i] - A[i, j]**2 / A[j, j]) 125 if j+1 <= n and min_num < - mu * gamma: 126 exec_phasetwo(A, L, P, I, E, j, n, tau, tau_bar, gamma, print_prefix, print_flag) 127 break 128 129 # Perform jth iteration of factorisation. 130 else: 131 jiter_factor(A, L, j, n) 132 133 # Debugging. 134 if print_flag >= 3: 135 print(print_prefix + "Factorised A:\n" + repr(A)) 136 print(print_prefix + "Factorised L:\n" + repr(L)) 137 #import sys 138 #sys.exit() 139 140 # The Cholesky factor. 141 L = dot(P, L) 142 143 # Debugging. 144 if print_flag >= 3: 145 print("\n" + print_prefix + "Fin") 146 print(print_prefix + "d2fk:\n" + repr(d2fk)) 147 print(print_prefix + "A:\n" + repr(A)) 148 print(print_prefix + "L:\n" + repr(L)) 149 try: 150 print(print_prefix + "chol:\n" + repr(cholesky(d2fk))) 151 except LinAlgError: 152 print(print_prefix + "Matrix is not positive definite - Cholesky decomposition cannot be computed.") 153 print(print_prefix + "E:\n" + repr(E)) 154 print(print_prefix + "P:\n" + repr(P)) 155 print(print_prefix + "Reconstructed d2fk:\n" + repr(dot(L, transpose(L)))) 156 print(print_prefix + "d2fk + E:\n" + repr(d2fk + dot(P, dot(E, P)))) 157 158 import sys 159 sys.exit() 160 161 # Calculate the Newton direction. 162 y = solve(L, dfk) 163 if return_matrix: 164 return -solve(transpose(L), y), dot(L, transpose(L)) 165 else: 166 return -solve(transpose(L), y)
167 168
169 -def exec_phasetwo(A, L, P, I, E, j, n, tau, tau_bar, gamma, print_prefix, print_flag):
170 """Phase 2 of the algorithm, A not positive definite.""" 171 172 # Debugging. 173 if print_flag >= 3: 174 print(print_prefix + "Phase 2.") 175 176 if j == n: 177 # delta = Enn. 178 delta = -A[n-1, n-1] + max(-tau*A[n-1, n-1] / (1.0 - tau), tau_bar*gamma) 179 if print_flag >= 3: 180 print(print_prefix + "delta: " + repr(delta)) 181 E[n-1, n-1] = delta 182 A[n-1, n-1] = A[n-1, n-1] + delta 183 L[n-1, n-1] = sqrt(A[n-1, n-1]) 184 185 else: 186 # Number of iterations performed in phase one (less 1). 187 k = j - 1 188 189 # Calculate the lower Gerschgorin bounds of Ak+1. 190 g = zeros(n, float64) 191 for i in range(k+1, n): 192 sum_Aij = 0.0 193 for s in range(k+1, i-1): 194 sum_Aij = sum_Aij + abs(A[i, j]) 195 sum_Aji = 0.0 196 for s in range(i+1, n): 197 sum_Aji = sum_Aji + abs(A[j, i]) 198 g[i] = A[i, i] - sum_Aij - sum_Aji 199 if print_flag >= 3: 200 print(print_prefix + "g: " + repr(g)) 201 202 # Modified Cholesky decomposition. 203 delta_prev = 0.0 204 for j in range(k+1, n-2): 205 # Pivot on the maximum lower Gerschgorin bound estimate. 206 i = j 207 for k in range(j, n): 208 if g[k] >= g[i]: 209 i = k 210 211 # Interchange row and column i and j. 212 p = 1.0 * I 213 if i != j: 214 if print_flag >= 3: 215 print(print_prefix + "Switching rows and columns of " + repr(i) + " and " + repr(j) + " of A and L.") 216 217 # Modify the permutation matrices by swaping rows. 218 row_p, row_P = 1.0*p[i], 1.0*P[i] 219 p[i], P[i] = p[j], P[j] 220 p[j], P[j] = row_p, row_P 221 222 # Permute A. 223 A = dot(p, dot(A, transpose(p))) 224 L = dot(p, dot(L, transpose(p))) 225 226 if print_flag >= 3: 227 print(print_prefix + "Permuted A:\n" + repr(A)) 228 print(print_prefix + "Permuted L:\n" + repr(L)) 229 230 # Calculate Ejj and add to diagonal. 231 normj = 0 232 for i in range(j+1, n): 233 normj = normj + abs(A[i, j]) 234 delta = max(0.0, -A[j, j] + max(normj, tau_bar*gamma), delta_prev) # delta = Enn. 235 if delta > 0.0: 236 A[j, j] = A[j, j] + delta 237 delta_prev = delta 238 if print_flag >= 3: 239 print(print_prefix + "delta: " + repr(delta)) 240 E[j, j] = delta 241 242 # Update Gerschgorin bound estimates. 243 if A[j, j] != normj: 244 temp = 1.0 - normj / A[j, j] 245 for i in range(j+1, n): 246 g[i] = g[i] + abs(A[i, j]) * temp 247 248 # Perform jth iteration of factorisation. 249 jiter_factor(A, L, j, n) 250 251 # Final 2*2 submatrix. 252 mini = take(take(A, [n-2, n-1], axis=0), [n-2, n-1], axis=1) 253 mini[0, 1] = mini[1, 0] 254 eigenvals = sort(eig(mini)) 255 delta = max(0.0, -eigenvals[0] + max(tau*(eigenvals[1] - eigenvals[0])/(1.0 - tau), tau_bar*gamma), delta_prev) 256 if delta > 0.0: 257 A[n-2, n-2] = A[n-2, n-2] + delta 258 A[n-1, n-1] = A[n-1, n-1] + delta 259 delta_prev = delta 260 if print_flag >= 3: 261 print(print_prefix + "delta: " + repr(delta)) 262 E[n-2, n-2] = E[n-1, n-1] = delta 263 264 L[n-2, n-2] = sqrt(A[n-2, n-2]) 265 L[n-1, n-2] = A[n-1, n-2] / L[n-2, n-2] 266 L[n-1, n-1] = sqrt(A[n-1, n-1] - L[n-1, n-2]**2) 267 268 if print_flag >= 3: 269 print(print_prefix + "mini:\n" + repr(mini)) 270 print(print_prefix + "eigenvals: " + repr(eigenvals)) 271 print(print_prefix + "Factorised A:\n" + repr(A)) 272 print(print_prefix + "Factorised L:\n" + repr(L))
273 274
275 -def jiter_factor(A, L, j, n):
276 """Perform jth iteration of factorisation.""" 277 278 L[j, j] = sqrt(A[j, j]) 279 for i in range(j+1, n): 280 L[i, j] = A[i, j] / L[j, j] 281 for k in range(j+1, i+1): 282 A[i, k] = A[i, k] - L[i, j]*L[k, j]
283