1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
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
30 from math import sqrt
31 from numpy import dot, transpose
32 from numpy.linalg import cholesky, eig, inv, LinAlgError, solve
33
34
35 -def gmw(dfk, d2fk, I, n, mach_acc, print_prefix, print_flag, return_matrix=0):
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
44
45
46
47
48
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
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
64 a = 1.0 * d2fk
65 r = 0.0 * d2fk
66 e = 0.0 * dfk
67 P = 1.0 * I
68
69
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
76 for j in range(n):
77
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
84 if q != j:
85
86 p = 1.0 * I
87
88
89 row_P = 1.0*P[:, q]
90 P[:, q] = P[:, j]
91 P[:, j] = row_P
92
93
94 row_p = 1.0*p[q]
95 p[q] = p[j]
96 p[j] = row_p
97
98
99 a = dot(p, dot(a, p))
100 r = dot(r, p)
101
102
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
110 if print_flag >= 3:
111 e[j] = dj - a[j, j]
112
113
114 r[j, j] = sqrt(dj)
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]
119
120
121 L = dot(P, transpose(r))
122
123
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
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