Coverage for frank/minimizer.py: 74%

1# Frankenstein: 1D disc brightness profile reconstruction from Fourier data 

2# using non-parametric Gaussian Processes 

3# 

4# Copyright (C) 2019-2020 R. Booth, J. Jennings, M. Tazzari 

5# 

6# This program is free software: you can redistribute it and/or modify 

7# it under the terms of the GNU General Public License as published by 

8# the Free Software Foundation, either version 3 of the License, or 

9# (at your option) any later version. 

10# 

11# This program is distributed in the hope that it will be useful, 

12# but WITHOUT ANY WARRANTY; without even the implied warranty of 

13# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

14# GNU General Public License for more details. 

15# 

16# You should have received a copy of the GNU General Public License 

17# along with this program. If not, see <https://www.gnu.org/licenses/> 

18# 

19"""This module contains methods solving non-linear minimization problems. 

20""" 

21import numpy as np 

22import scipy.linalg 

23 

24class BaseLineSearch(object): 

25 """Base class for back-tracking line searches. 

26  

27 Parameters 

28 ---------- 

29 reduce_step: function (dx, x) --> dx, optional 

30 Can be provided to limit the maximum step allowed in newton iterations 

31 given the current position (x) and suggested step (dx). No limiting is 

32 done if this function is not provided.  

33 """ 

34 

35 def __init__(self, reduce_step=None): 

36 self.reduction = None 

37 

38 if reduce_step is None: 

39 self.reduce_step = lambda dx, _: dx 

40 else: 

41 self.reduce_step = reduce_step 

42 

43 

44class LineSearch(BaseLineSearch): 

45 """Back-tracking line search for Newton's method and gradient descent. 

46 

47 Iteratively changes the step size used in the update until an acceptable 

48 reduction in the cost function is found. 

49 

50 Parameters 

51 ---------- 

52 armijo : float, default = 1e-4 

53 Coefficient used when deciding whether to accept a new solution.  

54 Smaller is more tolerant. 

55 min_step_frac : foat, default = 0.1 

56 Minimum amount that the step length can be reduced by in a single 

57 iteration when trying to find an improved guess. 

58 reduce_step: function (dx, x) --> dx, optional 

59 Can be provided to limit the maximum step allowed in newton iterations 

60 given the current position (x) and suggested step (dx). No limiting is 

61 done if this function is not provided.  

62 

63 Notes 

64 ----- 

65 This implementation is baeed on numerical recipes. 

66 """ 

67 

68 def __init__(self, armijo=1e-4, min_step_frac=0.1, reduce_step=None): 

69 super(LineSearch, self).__init__(reduce_step) 

70 

71 self._armijo = armijo 

72 self._l_min = min_step_frac 

73 

74 def __call__(self, func, jac, x0, p, f0=None, root=True): 

75 """Find a good step using backtracking. 

76 

77 Parameters 

78 ---------- 

79 func : function, 

80 The function that we are trying to find the root of 

81 jac : Jacobian object, 

82 Must provide a "dot" method that returns the dot-product of the 

83 Jacobian with a vector. 

84 x0 : array 

85 Current best guess for the solution 

86 p : array 

87 Step direction. 

88 f0 : array, optional. 

89 Evaluation of func at x0, i.e. func(x0). If not provided then this 

90 will be evaluated internally. 

91 root :  

92 

93 Returns 

94 ------- 

95 x_new : array 

96 Recommended new point 

97 f_new : array 

98 func(x_new) 

99 nfev : int 

100 Number of function evaluations 

101 failed : bool 

102 Whether the line-search failed to produce a guess 

103 """ 

104 def eval(x): 

105 fvec = func(x) 

106 if root: 

107 f = 0.5*np.dot(fvec, fvec) 

108 return fvec, f 

109 else: 

110 return fvec, fvec 

111 

112 nfev = 0 

113 if f0 is None: 

114 f0, cost = eval(x0) 

115 nfev += 1 

116 else: 

117 if root: 

118 cost = 0.5*np.dot(f0, f0) 

119 else: 

120 cost = f0 

121 

122 # First make sure the step isn't trying to change x by too much 

123 p = self.reduce_step(p, x0) 

124 

125 # Compute the expected change in f due to the step p assuming f is 

126 # exactly described by its linear expansion about the root: 

127 if root: 

128 delta_f = np.dot(f0, np.dot(jac, p)) 

129 else: 

130 delta_f = np.dot(jac, p) 

131 

132 if delta_f > 0: 

133 raise ValueError("Round off in slope calculation") 

134 

135 # Start with the full newton step 

136 lam = 1.0 

137 cost_save = lam_save = None 

138 while True: 

139 

140 x_new = x0 + lam*p 

141 if np.all(x_new == x0): 

142 return x0, f0, nfev, True 

143 

144 f_new, cost_new = eval(x_new) 

145 nfev += 1 

146 

147 # Have we got an acceptable step? 

148 if cost_new <= (cost + self._armijo*lam*delta_f): 

149 self.reduction = lam 

150 return x_new, f_new, nfev, False 

151 

152 # Try back tracking: 

153 if lam == 1.0: 

154 # First attempt, make a second order model of the cost 

155 # against lam. 

156 lam_new = - 0.5*delta_f / (cost_new - cost - delta_f) 

157 else: 

158 # Use a third order model 

159 r1 = (cost_new - cost - lam*delta_f)/(lam*lam) 

160 r2 = (cost_save - cost - lam_save*delta_f)/(lam_save*lam_save) 

161 

162 a = (r1 - r2) / (lam - lam_save) 

163 b = (lam*r2-lam_save*r1) / (lam - lam_save) 

164 

165 if a == 0: 

166 lam_new = - 0.5*delta_f / b 

167 else: 

168 d = b*b - 3*a*delta_f 

169 if d < 0: 

170 lam_new = 0.5*lam 

171 elif b <= 0: 

172 lam_new = (-b+np.sqrt(d))/(3*a) 

173 else: 

174 lam_new = -1 * delta_f / (b + np.sqrt(d)) 

175 

176 lam_new = min(0.5*lam, lam_new) 

177 

178 if np.isnan(lam_new): 

179 #raise ValueError("Nan in line search") 

180 lam_new = self._l_min*lam 

181 

182 lam_save = lam 

183 cost_save = cost_new 

184 lam = max(lam_new, self._l_min*lam) 

185 

186 

187def MinimizeNewton(fun, jac, hess, guess, line_search, 

188 max_step=10**5, max_hev=1000, tol=1e-5): 

189 """Minimize a function using Newton's method with back-tracking 

190 

191 Note, if Newton's method fails to find improvement a gradient descent step 

192 with back-tracking will be tried instead. 

193 

194 Convergence is assumed when the r.m.s. jacobian is below the requested 

195 theshold. MinimizeNewton will guess if convergence fails due to too many  

196 iterations, hessian calculations, or failure to improve the solution. 

197 

198 Parameters 

199 ---------- 

200 fun : function, of N variables returning a scalar 

201 Function to minimize 

202 jac : function, of N variables returning an array of N variables 

203 Gradient (jacobian) of the function 

204 hess : function, of N variables returning an array of NxN variables 

205 Hession of the function 

206 guess : array, lent N 

207 Initial guess for the solution 

208 line_search : LineSearch object 

209 Back-tracking line_search object 

210 max_step : int 

211 Maximum number of steps allowed 

212 max_hev : int 

213 Maximum number of Hessian evaluations allowed 

214 tol : float, 

215 Tolerance parameter for convergence. Convergence occurs when: 

216 np.sqrt(np.mean(jac(x)**2)) < tol 

217  

218 Returns 

219 ------- 

220 x : array, 

221 Best guess of the solution 

222 status : int, 

223 0 : Success 

224 1 : Iteration failed to improve estimate 

225 2 : Too many iterations 

226 3 : Too many hessian evaluations 

227 """ 

228 need_hess = True 

229 nfev = 1 

230 nhess = 0 

231 x = guess 

232 fx = fun(x) 

233 for nstep in range(max_step): 

234 if need_hess: 

235 if nhess == max_hev: 

236 return x, (3, nstep, nfev, nhess) 

237 

238 j_sol = (scipy.linalg.lu_factor(hess(x)), scipy.linalg.lu_solve) 

239 nhess += 1 

240 

241 jx = jac(x) 

242 dx = j_sol[1](j_sol[0], -jx) 

243 

244 if np.dot(jx, dx) < 0: 

245 x, fx, fev, failed = line_search(fun, jx, x, dx, fx, False) 

246 nfev += fev 

247 else: 

248 failed = True 

249 

250 # Use gradient descent when Newton's method fails 

251 if failed: 

252 x, fx, fev, failed_descent = line_search(fun, jx, x, -jx, fx, False) 

253 nfev += fev 

254 

255 

256 # Try again once more with a different back tracking algorithm  

257 if failed_descent: 

258 dx = line_search.reduce_step(-jx, x) 

259 alpha = 2**-4 

260 for _ in range(10): 

261 xn = x + dx 

262 fn = fun(xn) 

263 nfev += 1 

264 if fn < fx: 

265 break 

266 else: 

267 dx *= alpha 

268 else: 

269 # Neither gradient descent nor Newton's method can improve 

270 # the solution 

271 return x, (1, nstep, nfev, nhess) 

272 

273 fx = fn 

274 x = xn 

275 

276 need_hess = failed or (line_search.reduction != 1.0) 

277 

278 # Small enough gradient 

279 #if np.sqrt(np.mean(jac(x)**2)) < tol: 

280 # return x, 0 

281 if (np.abs(jac(x))*np.abs(x)).max() < tol * max(np.abs(fx), 1): 

282 return x, (0, nstep, nfev, nhess) 

283 

284 return x, (2, nstep, nfev, nhess)