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