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

Source Code for Module minimise.hessian_mods.gmw81

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2003, 2004 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   
 24  from LinearAlgebra import LinAlgError, cholesky_decomposition, eigenvalues, inverse, solve_linear_equations 
 25  from math import sqrt 
 26  from Numeric import Float64, array, dot, identity, matrixmultiply, transpose 
 27   
 28   
29 -def gmw(dfk, d2fk, I, n, mach_acc, print_prefix, print_flag, return_matrix=0):
30 """The Gill, Murray, and Wright modified Cholesky algorithm. 31 32 Algorithm 6.5 from page 148 of 'Numerical Optimization' by Jorge Nocedal and Stephen J. 33 Wright, 1999, 2nd ed. 34 35 Returns the modified Newton step. 36 """ 37 38 # Test matrix. 39 #d2fk = array([[4, 2, 1], [2, 6, 3], [1, 3, -0.004]], Float64) 40 #n = len(d2fk) 41 #I = identity(n, Float64) 42 43 # Calculate gamma(A) and xi(A). 44 gamma = 0.0 45 xi = 0.0 46 for i in xrange(n): 47 gamma = max(abs(d2fk[i, i]), gamma) 48 for j in xrange(i+1, n): 49 xi = max(abs(d2fk[i, j]), xi) 50 51 # Calculate delta and beta. 52 delta = mach_acc * max(gamma + xi, 1.0) 53 if n == 1: 54 beta = sqrt(max(gamma, mach_acc)) 55 else: 56 beta = sqrt(max(gamma, xi / sqrt(n**2 - 1.0), mach_acc)) 57 58 # Initialise data structures. 59 a = 1.0 * d2fk 60 r = 0.0 * d2fk 61 e = 0.0 * dfk 62 P = 1.0 * I 63 64 # Debugging. 65 if print_flag >= 3: 66 old_eigen = eigenvalues(d2fk) 67 print print_prefix + "dfk: " + `dfk` 68 print print_prefix + "d2fk:\n" + `d2fk` 69 70 # Main loop. 71 for j in xrange(n): 72 # Row and column swapping, find the index > j of the largest diagonal element. 73 q = j 74 for i in xrange(j+1, n): 75 if abs(a[i, i]) >= abs(a[q, q]): 76 q = i 77 78 # Interchange row and column j and q (if j != q). 79 if q != j: 80 # Temporary permutation matrix for swaping 2 rows or columns. 81 p = 1.0 * I 82 83 # Modify the permutation matrix P by swaping columns. 84 row_P = 1.0*P[:, q] 85 P[:, q] = P[:, j] 86 P[:, j] = row_P 87 88 # Modify the permutation matrix p by swaping rows (same as columns because p = pT). 89 row_p = 1.0*p[q] 90 p[q] = p[j] 91 p[j] = row_p 92 93 # Permute a and r (p = pT). 94 a = dot(p, dot(a, p)) 95 r = dot(r, p) 96 97 # Calculate dj. 98 theta_j = 0.0 99 if j < n-1: 100 for i in xrange(j+1, n): 101 theta_j = max(theta_j, abs(a[j, i])) 102 dj = max(abs(a[j, j]), (theta_j/beta)**2, delta) 103 104 # Calculate e (not really needed!). 105 if print_flag >= 3: 106 e[j] = dj - a[j, j] 107 108 # Calculate row j of r and update a. 109 r[j, j] = sqrt(dj) # Damned sqrt introduces roundoff error. 110 for i in xrange(j+1, n): 111 r[j, i] = a[j, i] / r[j, j] 112 for k in xrange(j+1, i+1): 113 a[i, k] = a[k, i] = a[k, i] - r[j, i] * r[j, k] # Keep matrix a symmetric. 114 115 # The Cholesky factor of d2fk. 116 L = dot(P, transpose(r)) 117 118 # Debugging. 119 if print_flag >= 3: 120 print print_prefix + "r:\n" + `r` 121 print print_prefix + "e: " + `e` 122 print print_prefix + "e rot: " + `dot(P, e)` 123 print print_prefix + "P:\n" + `P` 124 print print_prefix + "L:\n" + `L` 125 print print_prefix + "L.LT:\n" + `dot(L, transpose(L))` 126 print print_prefix + "L.LT - d2fk:\n" + `dot(L, transpose(L)) - d2fk` 127 print print_prefix + "dot(P, chol(L.LT)):\n" + `dot(P, cholesky_decomposition(dot(L, transpose(L))))` 128 print print_prefix + "dot(P, chol(P.L.LT.PT)):\n" + `dot(P, cholesky_decomposition(dot(P, dot(dot(L, transpose(L)), transpose(P)))))` 129 E = 0.0 * d2fk 130 for i in xrange(n): 131 E[i, i] = e[i] 132 E = dot(P, dot(E, transpose(P))) 133 print print_prefix + "E:\n" + `E` 134 print print_prefix + "PT.d2fk+E.P:\n" + `dot(transpose(P), dot(d2fk+E, P))` 135 try: 136 chol = cholesky_decomposition(dot(transpose(P), dot(d2fk+E, P))) 137 print print_prefix + "chol(PT.d2fk+E.P):\n" + `chol` 138 print print_prefix + "rT - chol(PT.d2fk+E.P):\n" + `transpose(r) - chol` 139 chol = dot(P, chol) 140 print print_prefix + "chol:\n" + `chol` 141 print print_prefix + "chol reconstructed:\n" + `dot(chol, transpose(chol))` 142 except LinAlgError: 143 print print_prefix + "Matrix is not positive definite - Cholesky decomposition cannot be computed." 144 eigen = eigenvalues(dot(L, transpose(L))) 145 print print_prefix + "Old eigenvalues: " + `old_eigen` 146 print print_prefix + "New eigenvalues: " + `eigen` 147 print print_prefix + "Newton dir: " + `-solve_linear_equations(transpose(L), solve_linear_equations(L, dfk))` 148 print print_prefix + "Newton dir using inverse: " + `-matrixmultiply(inverse(d2fk+E), dfk)` 149 150 # Calculate the Newton direction. 151 y = solve_linear_equations(L, dfk) 152 if return_matrix: 153 return -solve_linear_equations(transpose(L), y), dot(L, transpose(L)) 154 else: 155 return -solve_linear_equations(transpose(L), y)
156