1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24 """Cholesky Hessian modification with added multiple of the identity.
25
26 This file is part of the U{minfx optimisation library<https://sourceforge.net/projects/minfx>}.
27 """
28
29
30 from numpy import dot, sqrt, trace, transpose
31 from numpy.linalg import LinAlgError, cholesky, solve
32
33
34 -def cholesky_mod(dfk, d2fk, I, n, print_prefix, print_flag, return_matrix=0):
35 """Cholesky with added multiple of the identity.
36
37 Algorithm 6.3 from page 145 of 'Numerical Optimization' by Jorge Nocedal and Stephen J. Wright, 1999, 2nd ed.
38
39 Returns the modified Newton step.
40 """
41
42
43 min_aii = 1e99
44 for i in range(n):
45 min_aii = min(d2fk[i, i], min_aii)
46
47
48 norm = sqrt(trace(dot(d2fk, d2fk)))
49 half_norm = norm / 2.0
50
51
52 if min_aii > 0.0:
53 tk = 0.0
54 else:
55 tk = half_norm
56
57
58 if print_flag >= 3:
59 print(print_prefix + "Frobenius norm: " + repr(norm))
60 print(print_prefix + "min aii: " + repr(min_aii))
61 print(print_prefix + "tk: " + repr(tk))
62
63
64 while True:
65 if print_flag >= 3:
66 print(print_prefix + "Iteration")
67
68
69 matrix = d2fk + tk * I
70
71
72 try:
73 L = cholesky(matrix)
74 if print_flag >= 3:
75 print(print_prefix + "\tCholesky matrix L:")
76 for i in range(n):
77 print(print_prefix + "\t\t" + repr(L[i]))
78 except LinAlgError:
79 if print_flag >= 3:
80 print(print_prefix + "\tLinearAlgebraError, matrix is not positive definite.")
81
82
83 tk = max(2.0*tk, half_norm)
84 else:
85
86 break
87
88
89 y = solve(L, dfk)
90 if return_matrix:
91 return -solve(transpose(L), y), matrix
92 else:
93 return -solve(transpose(L), y)
94