1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24 from LinearAlgebra import LinAlgError, cholesky_decomposition, eigenvalues, inverse, solve_linear_equations
25 from math import sqrt
26 from Numeric import Float64, array, dot, identity, matrixmultiply, transpose
27
28
29 -def gmw(dfk, d2fk, I, n, mach_acc, print_prefix, print_flag, return_matrix=0):
30 """The Gill, Murray, and Wright modified Cholesky algorithm.
31
32 Algorithm 6.5 from page 148 of 'Numerical Optimization' by Jorge Nocedal and Stephen J.
33 Wright, 1999, 2nd ed.
34
35 Returns the modified Newton step.
36 """
37
38
39
40
41
42
43
44 gamma = 0.0
45 xi = 0.0
46 for i in xrange(n):
47 gamma = max(abs(d2fk[i, i]), gamma)
48 for j in xrange(i+1, n):
49 xi = max(abs(d2fk[i, j]), xi)
50
51
52 delta = mach_acc * max(gamma + xi, 1.0)
53 if n == 1:
54 beta = sqrt(max(gamma, mach_acc))
55 else:
56 beta = sqrt(max(gamma, xi / sqrt(n**2 - 1.0), mach_acc))
57
58
59 a = 1.0 * d2fk
60 r = 0.0 * d2fk
61 e = 0.0 * dfk
62 P = 1.0 * I
63
64
65 if print_flag >= 3:
66 old_eigen = eigenvalues(d2fk)
67 print print_prefix + "dfk: " + `dfk`
68 print print_prefix + "d2fk:\n" + `d2fk`
69
70
71 for j in xrange(n):
72
73 q = j
74 for i in xrange(j+1, n):
75 if abs(a[i, i]) >= abs(a[q, q]):
76 q = i
77
78
79 if q != j:
80
81 p = 1.0 * I
82
83
84 row_P = 1.0*P[:, q]
85 P[:, q] = P[:, j]
86 P[:, j] = row_P
87
88
89 row_p = 1.0*p[q]
90 p[q] = p[j]
91 p[j] = row_p
92
93
94 a = dot(p, dot(a, p))
95 r = dot(r, p)
96
97
98 theta_j = 0.0
99 if j < n-1:
100 for i in xrange(j+1, n):
101 theta_j = max(theta_j, abs(a[j, i]))
102 dj = max(abs(a[j, j]), (theta_j/beta)**2, delta)
103
104
105 if print_flag >= 3:
106 e[j] = dj - a[j, j]
107
108
109 r[j, j] = sqrt(dj)
110 for i in xrange(j+1, n):
111 r[j, i] = a[j, i] / r[j, j]
112 for k in xrange(j+1, i+1):
113 a[i, k] = a[k, i] = a[k, i] - r[j, i] * r[j, k]
114
115
116 L = dot(P, transpose(r))
117
118
119 if print_flag >= 3:
120 print print_prefix + "r:\n" + `r`
121 print print_prefix + "e: " + `e`
122 print print_prefix + "e rot: " + `dot(P, e)`
123 print print_prefix + "P:\n" + `P`
124 print print_prefix + "L:\n" + `L`
125 print print_prefix + "L.LT:\n" + `dot(L, transpose(L))`
126 print print_prefix + "L.LT - d2fk:\n" + `dot(L, transpose(L)) - d2fk`
127 print print_prefix + "dot(P, chol(L.LT)):\n" + `dot(P, cholesky_decomposition(dot(L, transpose(L))))`
128 print print_prefix + "dot(P, chol(P.L.LT.PT)):\n" + `dot(P, cholesky_decomposition(dot(P, dot(dot(L, transpose(L)), transpose(P)))))`
129 E = 0.0 * d2fk
130 for i in xrange(n):
131 E[i, i] = e[i]
132 E = dot(P, dot(E, transpose(P)))
133 print print_prefix + "E:\n" + `E`
134 print print_prefix + "PT.d2fk+E.P:\n" + `dot(transpose(P), dot(d2fk+E, P))`
135 try:
136 chol = cholesky_decomposition(dot(transpose(P), dot(d2fk+E, P)))
137 print print_prefix + "chol(PT.d2fk+E.P):\n" + `chol`
138 print print_prefix + "rT - chol(PT.d2fk+E.P):\n" + `transpose(r) - chol`
139 chol = dot(P, chol)
140 print print_prefix + "chol:\n" + `chol`
141 print print_prefix + "chol reconstructed:\n" + `dot(chol, transpose(chol))`
142 except LinAlgError:
143 print print_prefix + "Matrix is not positive definite - Cholesky decomposition cannot be computed."
144 eigen = eigenvalues(dot(L, transpose(L)))
145 print print_prefix + "Old eigenvalues: " + `old_eigen`
146 print print_prefix + "New eigenvalues: " + `eigen`
147 print print_prefix + "Newton dir: " + `-solve_linear_equations(transpose(L), solve_linear_equations(L, dfk))`
148 print print_prefix + "Newton dir using inverse: " + `-matrixmultiply(inverse(d2fk+E), dfk)`
149
150
151 y = solve_linear_equations(L, dfk)
152 if return_matrix:
153 return -solve_linear_equations(transpose(L), y), dot(L, transpose(L))
154 else:
155 return -solve_linear_equations(transpose(L), y)
156