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 Warning: This code is not currently functional.
27
28 This file is part of the U{minfx optimisation library<https://gna.org/projects/minfx>}.
29 """
30
31
32 from numpy import dot, sort, sqrt, transpose
33 from numpy.linalg import eig, inv
34
35
36 -def gmw_old(dfk, d2fk, I, n, mach_acc, print_prefix, print_flag, return_matrix=0):
37 """Modified Cholesky Hessian modification.
38
39 Algorithm 6.5 from page 148.
40
41 Somehow the data structures l, d, and c can be stored in d2fk!
42 """
43
44
45
46
47 if print_flag >= 3:
48 print("d2fk: " + repr(d2fk))
49 eigen = eig(d2fk)
50 eigenvals = sort(eigen[0])
51 print("Eigenvalues: " + repr(eigenvals))
52
53
54 gamma = 0.0
55 xi = 0.0
56 for i in range(n):
57 gamma = max(abs(d2fk[i, i]), gamma)
58 for j in range(i+1, n):
59 xi = max(abs(d2fk[i, j]), xi)
60
61
62 delta = mach_acc * max(gamma + xi, 1)
63 beta = sqrt(max(gamma, xi / sqrt(n**2 - 1.0), mach_acc))
64
65
66 d2fk_orig = 1.0 * d2fk
67 c = 0.0 * d2fk
68 d = 0.0 * dfk
69 l = 1.0 * I
70 e = 0.0 * dfk
71 P = 1.0 * I
72
73
74 for k in range(n):
75 c[k, k] = d2fk[k, k]
76
77
78 for j in range(n):
79 if print_flag >= 3:
80 print("\n<j: " + repr(j) + ">")
81
82
83 p = 1.0 * I
84 q = j
85 for i in range(j, n):
86 if abs(c[q, q]) <= abs(c[i, i]):
87 q = i
88 if print_flag >= 3:
89 print("Row and column swapping.")
90 print("i range: " + repr(range(j, n)))
91 print("q: " + repr(q))
92 print("a: " + repr(d2fk))
93 print("c: " + repr(c))
94 print("d: " + repr(d))
95 print("l: " + repr(l))
96 if q != j:
97
98 temp_p, temp_P = 1.0*p[:, q], 1.0*P[:, q]
99 p[:, q], P[:, q] = p[:, j], P[:, j]
100 p[:, j], P[:, j] = temp_p, temp_P
101
102
103 c = dot(p, dot(c, p))
104 d2fk = dot(p, dot(d2fk, p))
105
106 if print_flag >= 3:
107 print("p: " + repr(p))
108 print("a(mod): " + repr(d2fk))
109 print("c(mod): " + repr(c))
110 print("d(mod): " + repr(d))
111 print("l(mod): " + repr(l))
112
113
114 if print_flag >= 3:
115 print("\nCalculate the elements of l.")
116 print("s range: " + repr(range(j)))
117 for s in range(j):
118 if print_flag >= 3:
119 print("s: " + repr(s))
120 l[j, s] = c[j, s] / d[s]
121 if print_flag >= 3:
122 print("l: " + repr(l))
123
124
125 if print_flag >= 3:
126 print("\nCalculate c[i, j].")
127 print("i range: " + repr(range(j+1, n)))
128 for i in range(j+1, n):
129 sum = 0.0
130 for s in range(j):
131 if print_flag >= 3:
132 print("s range: " + repr(range(j)))
133 sum = sum + l[j, s] * c[i, s]
134 c[i, j] = d2fk[i, j] - sum
135 if print_flag >= 3:
136 print("c: " + repr(c))
137
138
139 if print_flag >= 3:
140 print("\nCalculate d.")
141 theta_j = 0.0
142 if j < n-1:
143 if print_flag >= 3:
144 print("j < n, " + repr(j) + " < " + repr(n-1))
145 print("i range: " + repr(range(j+1, n)))
146 for i in range(j+1, n):
147 theta_j = max(theta_j, abs(c[i, j]))
148 else:
149 if print_flag >= 3:
150 print("j >= n, " + repr(j) + " >= " + repr(n-1))
151 d[j] = max(abs(c[j, j]), (theta_j/beta)**2, delta)
152 if print_flag >= 3:
153 print("theta_j: " + repr(theta_j))
154 print("d: " + repr(d))
155
156
157 if print_flag >= 3:
158 print("\nCalculate c[i, i].")
159 if j < n-1:
160 if print_flag >= 3:
161 print("j < n, " + repr(j) + " < " + repr(n-1))
162 print("i range: " + repr(range(j+1, n)))
163 for i in range(j+1, n):
164 c[i, i] = c[i, i] - c[i, j]**2 / d[j]
165 else:
166 if print_flag >= 3:
167 print("j >= n, " + repr(j) + " >= " + repr(n-1))
168 if print_flag >= 3:
169 print("c: " + repr(c))
170
171
172 e[j] = d[j] - c[j, j]
173
174 for i in range(n):
175 d2fk[i, i] = d2fk[i, i] + e[i]
176 d2fk = dot(P, dot(d2fk, transpose(P)))
177
178
179 if print_flag >= 3:
180 print("\nFin:")
181 print("gamma(A): " + repr(gamma))
182 print("xi(A): " + repr(xi))
183 print("delta: " + repr(delta))
184 print("beta: " + repr(beta))
185 print("c: " + repr(c))
186 print("d: " + repr(d))
187 print("l: " + repr(l))
188 temp = 0.0 * I
189 for i in range(len(d)):
190 temp[i, i] = d[i]
191 print("m: " + repr(dot(l, sqrt(temp))))
192 print("mmT: " + repr(dot(dot(l, sqrt(temp)), transpose(dot(l, sqrt(temp))))))
193 print("P: " + repr(P))
194 print("e: " + repr(e))
195 print("d2fk: " + repr(d2fk))
196 eigen = eig(d2fk)
197 eigenvals = sort(eigen[0])
198 print("Eigenvalues: " + repr(eigenvals))
199 import sys
200 sys.exit()
201
202
203 if return_matrix:
204 return -dot(inv(d2fk), dfk), d2fk
205 else:
206 return -dot(inv(d2fk), dfk)
207