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 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
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
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
50 A = 1.0 * d2fk
51 L = 0.0 * d2fk
52
53
54 P = 1.0 * I
55
56
57 gamma = abs(d2fk[0, 0])
58 for i in range(n):
59 gamma = max(abs(d2fk[i, i]), gamma)
60
61
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
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
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
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
97 i = j
98 for k in range(j, n):
99 if A[k, k] >= A[i, i]:
100 i = k
101
102
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
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
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
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
130 else:
131 jiter_factor(A, L, j, n)
132
133
134 if print_flag >= 3:
135 print(print_prefix + "Factorised A:\n" + repr(A))
136 print(print_prefix + "Factorised L:\n" + repr(L))
137
138
139
140
141 L = dot(P, L)
142
143
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
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
173 if print_flag >= 3:
174 print(print_prefix + "Phase 2.")
175
176 if j == n:
177
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
187 k = j - 1
188
189
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
203 delta_prev = 0.0
204 for j in range(k+1, n-2):
205
206 i = j
207 for k in range(j, n):
208 if g[k] >= g[i]:
209 i = k
210
211
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
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
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
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)
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
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
249 jiter_factor(A, L, j, n)
250
251
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
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