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