minimise.hessian

1 ############################################################################### 2 # # 3 # Copyright (C) 2003 Edward d'Auvergne # 4 # # 5 # This file is part of the program relax. # 6 # # 7 # relax is free software; you can redistribute it and/or modify # 8 # it under the terms of the GNU General Public License as published by # 9 # the Free Software Foundation; either version 2 of the License, or # 10 # (at your option) any later version. # 11 # # 12 # relax is distributed in the hope that it will be useful, # 13 # but WITHOUT ANY WARRANTY; without even the implied warranty of # 14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 15 # GNU General Public License for more details. # 16 # # 17 # You should have received a copy of the GNU General Public License # 18 # along with relax; if not, write to the Free Software # 19 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA # 20 # # 21 ############################################################################### 22 23 from LinearAlgebra import LinAlgError, cholesky_decomposition, eigenvalues, solve_linear_equations 24 from MLab import tril 25 from Numeric import Float64, array, dot, identity, sort, sqrt, take, transpose, zeros 26 27

28 -def se99(dfk, d2fk, I, n, tau, tau_bar, mu, print_prefix, print_flag, return_matrix=0):

29 """A revised modified cholesky factorisation algorithm. 30 31 Schnabel, R. B. and Eskow, E. 1999, A revised modifed cholesky factorisation algorithm. 32 SIAM J. Optim. 9, 1135-1148. 33 """ 34 35 # Test matrix. 36 d2fk = array([[4, 2, 1], [2, 6, 3], [1, 3, -0.004]], Float64) 37 d2fk = array([[ 1890.3, -1750.6, -315.8, 3000.3], 38 [ -1705.6, 1538.3, 284.9, -2706.6], 39 [ -315.8, 284.9, 52.5, -501.2], 40 [ 3000.3, -2706.6, -501.2, 4760.8]], Float64) 41 n = len(d2fk) 42 I = identity(n, Float64) 43 44 # Create the matrices A and L. 45 A = 1.0 * d2fk 46 L = 0.0 * d2fk 47 48 # Permutation matrix. 49 P = 1.0 * I 50 51 # Calculate gamma. 52 gamma = abs(d2fk[0, 0]) 53 for i in xrange(n): 54 gamma = max(abs(d2fk[i, i]), gamma) 55 56 # Debugging. 57 E = 0.0 * d2fk 58 if print_flag >= 3: 59 print print_prefix + "tau: " + `tau` 60 print print_prefix + "tau_bar: " + `tau_bar` 61 print print_prefix + "mu: " + `mu` 62 print print_prefix + "d2fk:\n" + `d2fk` 63 print print_prefix + "Orig A:\n" + `A` 64 print print_prefix + "Orig L:\n" + `L` 65 print print_prefix + "gamma: " + `gamma` 66 67 68 # Phase one, A potentially positive definite. 69 ############################################# 70 71 for j in xrange(n): 72 if print_flag >= 3: 73 print "\n" + print_prefix + "Iteration j = " + `j` 74 print print_prefix + "Init A:\n" + `A` 75 print print_prefix + "Init L:\n" + `L` 76 print print_prefix + "Phase 1." 77 78 # Calculate max_Aii and min_Aii 79 max_Aii = A[j, j] 80 min_Aii = A[j, j] 81 for i in xrange(j, n): 82 max_Aii = max(A[i, i], max_Aii) 83 min_Aii = min(A[i, i], min_Aii) 84 85 # Test for phase 2, A not positive definite. 86 if max_Aii < tau_bar * gamma or min_Aii < - mu * max_Aii: 87 exec_phasetwo(A, L, P, I, E, j, n, tau, tau_bar, gamma, print_prefix, print_flag) 88 break 89 90 else: 91 # Pivot on maximum diagonal of remaining submatrix. 92 i = j 93 for k in xrange(j, n): 94 if A[k, k] >= A[i, i]: 95 i = k 96 97 # Interchange row and column i and j. 98 p = 1.0 * I 99 if i != j: 100 if print_flag >= 3: 101 print print_prefix + "Switching rows and columns of " + `i` + " and " + `j` + " of A and L." 102 103 # Modify the permutation matrices by swaping rows. 104 row_p, row_P = 1.0*p[i], 1.0*P[i] 105 p[i], P[i] = p[j], P[j] 106 p[j], P[j] = row_p, row_P 107 108 # Permute A. 109 A = dot(p, dot(A, transpose(p))) 110 L = dot(p, dot(L, transpose(p))) 111 112 if print_flag >= 3: 113 print print_prefix + "Permuted A:\n" + `A` 114 print print_prefix + "Permuted L:\n" + `L` 115 116 # Test for phase 2 again. 117 min_num = 1e99 118 for i in xrange(j+1, n): 119 min_num = min(min_num, A[i, i] - A[i, j]**2 / A[j, j]) 120 if j+1 <= n and min_num < - mu * gamma: 121 exec_phasetwo(A, L, P, I, E, j, n, tau, tau_bar, gamma, print_prefix, print_flag) 122 break 123 124 # Perform jth iteration of factorisation. 125 else: 126 jiter_factor(A, L, j, n) 127 128 # Debugging. 129 if print_flag >= 3: 130 print print_prefix + "Factorised A:\n" + `A` 131 print print_prefix + "Factorised L:\n" + `L` 132 #import sys 133 #sys.exit() 134 135 # The Cholesky factor. 136 L = dot(P, L) 137 138 # Debugging. 139 if print_flag >= 3: 140 print "\n" + print_prefix + "Fin" 141 print print_prefix + "d2fk:\n" + `d2fk` 142 print print_prefix + "A:\n" + `A` 143 print print_prefix + "L:\n" + `L` 144 try: 145 print print_prefix + "chol:\n" + `cholesky_decomposition(d2fk)` 146 except LinAlgError: 147 print print_prefix + "Matrix is not positive definite - Cholesky decomposition cannot be computed." 148 print print_prefix + "E:\n" + `E` 149 print print_prefix + "P:\n" + `P` 150 print print_prefix + "Reconstructed d2fk:\n" + `dot(L, transpose(L))` 151 print print_prefix + "d2fk + E:\n" + `d2fk + dot(P, dot(E, P))` 152 153 import sys 154 sys.exit() 155 156 # Calculate the Newton direction. 157 y = solve_linear_equations(L, dfk) 158 if return_matrix: 159 return -solve_linear_equations(transpose(L), y), dot(L, transpose(L)) 160 else: 161 return -solve_linear_equations(transpose(L), y)

162 163

164 -def exec_phasetwo(A, L, P, I, E, j, n, tau, tau_bar, gamma, print_prefix, print_flag):

165 """Phase 2 of the algorithm, A not positive definite.""" 166 167 # Debugging. 168 if print_flag >= 3: 169 print print_prefix + "Phase 2." 170 171 if j == n: 172 # delta = Enn. 173 delta = -A[n-1, n-1] + max(-tau*A[n-1, n-1] / (1.0 - tau), tau_bar*gamma) 174 if print_flag >= 3: 175 print print_prefix + "delta: " + `delta` 176 E[n-1, n-1] = delta 177 A[n-1, n-1] = A[n-1, n-1] + delta 178 L[n-1, n-1] = sqrt(A[n-1, n-1]) 179 180 else: 181 # Number of iterations performed in phase one (less 1). 182 k = j - 1 183 184 # Calculate the lower Gerschgorin bounds of Ak+1. 185 g = zeros(n, Float64) 186 for i in xrange(k+1, n): 187 sum_Aij = 0.0 188 for s in xrange(k+1, i-1): 189 sum_Aij = sum_Aij + abs(A[i, j]) 190 sum_Aji = 0.0 191 for s in xrange(i+1, n): 192 sum_Aji = sum_Aji + abs(A[j, i]) 193 g[i] = A[i, i] - sum_Aij - sum_Aji 194 if print_flag >= 3: 195 print print_prefix + "g: " + `g` 196 197 # Modified Cholesky decomposition. 198 delta_prev = 0.0 199 for j in xrange(k+1, n-2): 200 # Pivot on the maximum lower Gerschgorin bound estimate. 201 i = j 202 for k in xrange(j, n): 203 if g[k] >= g[i]: 204 i = k 205 206 # Interchange row and column i and j. 207 p = 1.0 * I 208 if i != j: 209 if print_flag >= 3: 210 print print_prefix + "Switching rows and columns of " + `i` + " and " + `j` + " of A and L." 211 212 # Modify the permutation matrices by swaping rows. 213 row_p, row_P = 1.0*p[i], 1.0*P[i] 214 p[i], P[i] = p[j], P[j] 215 p[j], P[j] = row_p, row_P 216 217 # Permute A. 218 A = dot(p, dot(A, transpose(p))) 219 L = dot(p, dot(L, transpose(p))) 220 221 if print_flag >= 3: 222 print print_prefix + "Permuted A:\n" + `A` 223 print print_prefix + "Permuted L:\n" + `L` 224 225 # Calculate Ejj and add to diagonal. 226 normj = 0 227 for i in xrange(j+1, n): 228 normj = normj + abs(A[i, j]) 229 delta = max(0.0, -A[j, j] + max(normj, tau_bar*gamma), delta_prev) # delta = Enn. 230 if delta > 0.0: 231 A[j, j] = A[j, j] + delta 232 delta_prev = delta 233 if print_flag >= 3: 234 print print_prefix + "delta: " + `delta` 235 E[j, j] = delta 236 237 # Update Gerschgorin bound estimates. 238 if A[j, j] != normj: 239 temp = 1.0 - normj / A[j, j] 240 for i in xrange(j+1, n): 241 g[i] = g[i] + abs(A[i, j]) * temp 242 243 # Perform jth iteration of factorisation. 244 jiter_factor(A, L, j, n) 245 246 # Final 2*2 submatrix. 247 mini = take(take(A, [n-2, n-1]), [n-2, n-1], 1) 248 mini[0, 1] = mini[1, 0] 249 eigenvals = sort(eigenvalues(mini)) 250 delta = max(0.0, -eigenvals[0] + max(tau*(eigenvals[1] - eigenvals[0])/(1.0 - tau), tau_bar*gamma), delta_prev) 251 if delta > 0.0: 252 A[n-2, n-2] = A[n-2, n-2] + delta 253 A[n-1, n-1] = A[n-1, n-1] + delta 254 delta_prev = delta 255 if print_flag >= 3: 256 print print_prefix + "delta: " + `delta` 257 E[n-2, n-2] = E[n-1, n-1] = delta 258 259 L[n-2, n-2] = sqrt(A[n-2, n-2]) 260 L[n-1, n-2] = A[n-1, n-2] / L[n-2, n-2] 261 L[n-1, n-1] = sqrt(A[n-1, n-1] - L[n-1, n-2]**2) 262 263 if print_flag >= 3: 264 print print_prefix + "mini:\n" + `mini` 265 print print_prefix + "eigenvals: " + `eigenvals` 266 print print_prefix + "Factorised A:\n" + `A` 267 print print_prefix + "Factorised L:\n" + `L`

268 269

270 -def jiter_factor(A, L, j, n):

271 """Perform jth iteration of factorisation.""" 272 273 L[j, j] = sqrt(A[j, j]) 274 for i in xrange(j+1, n): 275 L[i, j] = A[i, j] / L[j, j] 276 for k in xrange(j+1, i+1): 277 A[i, k] = A[i, k] - L[i, j]*L[k, j]

278

Source Code for Module minimise.hessian_mods.se99