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