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 cholesky_decomposition, solve_linear_equations
25 from Numeric import dot, sqrt, trace, transpose
26
27
28 -def cholesky(dfk, d2fk, I, n, print_prefix, print_flag, return_matrix=0):
29 """Cholesky with added multiple of the identity.
30
31 Algorithm 6.3 from page 145 of 'Numerical Optimization' by Jorge Nocedal and Stephen J.
32 Wright, 1999, 2nd ed.
33
34 Returns the modified Newton step.
35 """
36
37
38 min_aii = 1e99
39 for i in xrange(n):
40 min_aii = min(d2fk[i, i], min_aii)
41
42
43 norm = sqrt(trace(dot(d2fk, d2fk)))
44 half_norm = norm / 2.0
45
46
47 if min_aii > 0.0:
48 tk = 0.0
49 else:
50 tk = half_norm
51
52
53 if print_flag >= 3:
54 print print_prefix + "Frobenius norm: " + `norm`
55 print print_prefix + "min aii: " + `min_aii`
56 print print_prefix + "tk: " + `tk`
57
58
59 while 1:
60 if print_flag >= 3:
61 print print_prefix + "Iteration"
62
63
64 matrix = d2fk + tk * I
65
66
67 try:
68 L = cholesky_decomposition(matrix)
69 if print_flag >= 3:
70 print print_prefix + "\tCholesky matrix L:"
71 for i in xrange(n):
72 print print_prefix + "\t\t" + `L[i]`
73 except:
74 if print_flag >= 3:
75 print print_prefix + "\tLinearAlgebraError, matrix is not positive definite."
76
77
78 tk = max(2.0*tk, half_norm)
79 else:
80
81 break
82
83
84 y = solve_linear_equations(L, dfk)
85 if return_matrix:
86 return -solve_linear_equations(transpose(L), y), matrix
87 else:
88 return -solve_linear_equations(transpose(L), y)
89