Coverage for gpaw/directmin/ls_etdm.py: 77%

190 statements  

« prev     ^ index     » next       coverage.py v7.7.1, created at 2025-07-08 00:17 +0000

1""" 

2Class for finding an optimal 

3step alpha length in the optimization problem: 

4x = x + alpha * p. 

5It is for lcao 

6""" 

7 

8import numpy as np 

9from gpaw.directmin.tools import minimum_cubic_interpol 

10 

11 

12def is_descent(phi_0, phi_j, eps=1.0e-2): 

13 return phi_j <= phi_0 + eps * abs(phi_0) 

14 

15 

16def is_descent_and_approximate_wolfe_conditions( 

17 der_phi_0, phi_0, der_phi_j, phi_j, eps=1.0e-6, delta=0.1, sigma=0.9): 

18 """ 

19 William W. Hager and Hongchao Zhang 

20 SIAM J. optim., 16(1), 170-192. 

21 """ 

22 

23 descent = is_descent(phi_0, phi_j, eps) 

24 awc = (2.0 * delta - 1.0) * der_phi_0 >= der_phi_j >= sigma * der_phi_0 

25 

26 return descent and awc 

27 

28 

29def get_slength(p, wfs, mode=None): 

30 if mode is None: 

31 mode = wfs.mode 

32 if mode == 'lcao': 

33 p_all_kpts = np.hstack([p[k] for k in p]) 

34 return np.linalg.norm(p_all_kpts) 

35 else: 

36 ret = 0.0 

37 for k in p: 

38 for val in p[k]: 

39 ret += np.real(wfs.integrate(val, val, global_integral=False)) 

40 ret = wfs.world.sum_scalar(ret) 

41 return np.sqrt(ret) 

42 

43 

44class MaxStep: 

45 

46 def __init__(self, evaluate_phi_and_der_phi, max_step=0.2): 

47 """ 

48 

49 :param evaluate_phi_and_der_phi: 

50 """ 

51 

52 self.evaluate_phi_and_der_phi = evaluate_phi_and_der_phi 

53 self.max_step = max_step 

54 self.name = 'max-step' 

55 

56 def todict(self): 

57 return {'name': self.name, 

58 'max_step': self.max_step} 

59 

60 def step_length_update(self, x, p, wfs, *args, mode=None, **kwargs): 

61 

62 kd = kwargs['kpdescr'] 

63 

64 slength = get_slength(p, wfs, mode) 

65 slength = kd.comm.max_scalar(slength) 

66 

67 a_star = self.max_step / slength if slength > self.max_step else 1.0 

68 

69 phi_star, der_phi_star, g_star = self.evaluate_phi_and_der_phi( 

70 x, p, a_star, wfs, *args) 

71 

72 return a_star, phi_star, der_phi_star, g_star 

73 

74 

75class StrongWolfeConditions(MaxStep): 

76 """ 

77 From a book of Jorge Nocedal and Stephen J. Wright 'Numerical 

78 Optimization' Second Edition, 2006 (p. 56) 

79 

80 This call should return a_star, phi_star, der_phi_star, g_star, 

81 where a_star is step length satisfied the strong Wolfe condts: 

82 

83 f(x_k + a_k p_k) <= f(x_k) + c_1 a_k grad f_k cdot p_k, 

84 

85 |grad f(x_k + a_k p_k) cdot p_k | <= c_2 |grad f_k cdot p_k|, 

86 

87 or descent conditions and approximate Wolfe conditions from 

88 

89 William W. Hager and Hongchao Zhang 

90 SIAM J. optim., 16(1), 170-192. 

91 

92 phi = f (x_k + a_k*p_k) 

93 der_phi = grad f(x_k + a_k p_k) cdot p_k 

94 g = grad f(x_k + a_k p_k) 

95 """ 

96 

97 def __init__(self, evaluate_phi_and_der_phi, 

98 c1=1.0e-4, c2=0.9, 

99 searchdirtype=None, max_iter=3, eps_dx=1.0e-10, 

100 eps_df=1.0e-10, use_descent_and_awc=True): 

101 """ 

102 :param evaluate_phi_and_der_phi: function which calculate 

103 phi, der_phi and g for given A_s, P_s, n_dim and alpha 

104 A_s[s] is skew-hermitian matrix, P_s[s] is matrix direction 

105 :param searchdirtype: used only in initial guess for alpha 

106 :param max_iter: maximum number of iterations 

107 :param eps_dx: length of minimal interval where alpha can 

108 be found 

109 :param eps_df: minimal change of function 

110 :param c1: see above 

111 :param c2: see above 

112 :param use_descent_and_awc: check descent and approximate Wolfe 

113 conditions 

114 

115 this class works fine, but these parameters eps_dx, eps_df 

116 might not be needed 

117 """ 

118 

119 super().__init__( 

120 evaluate_phi_and_der_phi) 

121 

122 self.max_iter = max_iter 

123 self.searchdirtype = searchdirtype 

124 self.eps_dx = eps_dx 

125 self.eps_df = eps_df 

126 self.c1 = c1 

127 self.c2 = c2 

128 self.use_descent_and_awc = use_descent_and_awc 

129 self.name = 'swc-awc' 

130 

131 def todict(self): 

132 return {'name': self.name, 

133 'max_iter': self.max_iter, 

134 'searchdirtype': self.searchdirtype, 

135 'eps_dx': self.eps_dx, 

136 'eps_df': self.eps_df, 

137 'c1': self.c1, 

138 'c2': self.c2} 

139 

140 def step_length_update(self, x, p, wfs, *args, **kwargs): 

141 c1 = self.c1 

142 c2 = self.c2 

143 phi_0 = kwargs['phi_0'] 

144 der_phi_0 = kwargs['der_phi_0'] 

145 phi_old = kwargs['phi_old'] 

146 der_phi_old = kwargs['der_phi_old'] 

147 alpha_old = kwargs['alpha_old'] 

148 alpha_max = kwargs['alpha_max'] 

149 alpha_1 = self.init_guess(phi_0=phi_0, der_phi_0=der_phi_0, 

150 phi_old=phi_old, 

151 der_phi_old=der_phi_old, 

152 alpha_old=alpha_old) 

153 i = 1 

154 

155 if phi_0 is None or der_phi_0 is None: 

156 phi_0, der_phi_0, g_0 = \ 

157 self.evaluate_phi_and_der_phi(x, p, 0.0, wfs, *args) 

158 

159 alpha = [0.0, alpha_1] 

160 phi_i_1 = phi_0 

161 der_phi_i_1 = der_phi_0 

162 max_iter = self.max_iter 

163 phi_max = None 

164 der_phi_max = None 

165 g_max = None 

166 

167 while True: 

168 if np.abs(alpha[-1] - alpha_max) < 1.0e-6 and phi_max is not None: 

169 phi_i, der_phi_i, g_i = phi_max, der_phi_max, g_max 

170 phi_max = None 

171 der_phi_max = None 

172 else: 

173 phi_i, der_phi_i, g_i = \ 

174 self.evaluate_phi_and_der_phi(x, p, alpha[i], wfs, *args) 

175 

176 if self.use_descent_and_awc: 

177 if is_descent_and_approximate_wolfe_conditions( 

178 der_phi_0, phi_0, der_phi_i, phi_i): 

179 a_star = alpha[i] 

180 phi_star = phi_i 

181 der_phi_star = der_phi_i 

182 g_star = g_i 

183 break 

184 

185 if phi_i > phi_0 + c1 * alpha[i] * der_phi_0 or \ 

186 (phi_i >= phi_i_1 and i > 1): 

187 a_star, phi_star, der_phi_star, g_star = \ 

188 self.zoom(alpha[i - 1], alpha[i], 

189 phi_i_1, der_phi_i_1, 

190 phi_i, der_phi_i, x, p, 

191 phi_0, der_phi_0, c1, c2, wfs, *args) 

192 break 

193 

194 if np.fabs(der_phi_i) <= -c2 * der_phi_0: 

195 a_star = alpha[i] 

196 phi_star = phi_i 

197 der_phi_star = der_phi_i 

198 g_star = g_i 

199 break 

200 

201 if der_phi_i >= 0.0: 

202 a_star, phi_star, der_phi_star, g_star = \ 

203 self.zoom(alpha[i], alpha[i - 1], 

204 phi_i, der_phi_i, 

205 phi_i_1, der_phi_i_1, 

206 x, p, phi_0, der_phi_0, c1, c2, wfs, *args) 

207 break 

208 

209 if i == max_iter: 

210 a_star = alpha[i] 

211 phi_star = phi_i 

212 der_phi_star = der_phi_i 

213 g_star = g_i 

214 break 

215 

216 if np.abs(alpha_max - alpha[i]) < 1.0e-6: 

217 alpha_max = 1.5 * alpha[i] 

218 

219 if phi_max is None or der_phi_max is None: 

220 phi_max, der_phi_max, g_max = \ 

221 self.evaluate_phi_and_der_phi(x, p, alpha_max, wfs, *args) 

222 

223 if self.use_descent_and_awc: 

224 if is_descent_and_approximate_wolfe_conditions( 

225 der_phi_0, phi_0, der_phi_max, phi_max): 

226 a_star = alpha_max 

227 phi_star = phi_max 

228 der_phi_star = der_phi_max 

229 g_star = g_max 

230 break 

231 

232 a_new = minimum_cubic_interpol( 

233 alpha[i], alpha_max, phi_i, phi_max, 

234 der_phi_i, der_phi_max) 

235 

236 alpha.append(a_new) 

237 phi_i_1 = phi_i 

238 der_phi_i_1 = der_phi_i 

239 

240 if np.abs(a_new - alpha_max) < 1.0e-6: 

241 phi_i = phi_max 

242 der_phi_i = der_phi_max 

243 g_i = g_max 

244 

245 if abs(alpha[-1] - alpha[-2]) < 1.0e-5: 

246 a_star = alpha[i] 

247 phi_star = phi_i 

248 der_phi_star = der_phi_i 

249 g_star = g_i 

250 break 

251 

252 i += 1 

253 

254 return a_star, phi_star, der_phi_star, g_star 

255 

256 def zoom(self, a_lo, a_hi, f_lo, df_lo, f_hi, df_hi, x, p, phi_0, 

257 der_phi_0, c1, c2, wfs, *args): 

258 

259 max_iter = self.max_iter 

260 i = 0 

261 

262 while True: 

263 a_j = minimum_cubic_interpol( 

264 a_lo, a_hi, f_lo, f_hi, df_lo, df_hi) 

265 

266 if a_j < 0.01: 

267 a_j = 0.1 

268 phi_j, der_phi_j, g_j = \ 

269 self.evaluate_phi_and_der_phi(x, p, a_j, wfs, *args) 

270 return a_j, phi_j, der_phi_j, g_j 

271 

272 phi_j, der_phi_j, g_j = \ 

273 self.evaluate_phi_and_der_phi(x, p, a_j, wfs, *args) 

274 

275 if self.use_descent_and_awc: 

276 if is_descent_and_approximate_wolfe_conditions( 

277 der_phi_0, phi_0, der_phi_j, phi_j): 

278 a_star = a_j 

279 phi_star = phi_j 

280 der_phi_star = der_phi_j 

281 g_star = g_j 

282 break 

283 

284 if phi_j > phi_0 + c1 * a_j * der_phi_0 or phi_j >= f_lo: 

285 a_hi = a_j 

286 f_hi = phi_j 

287 df_hi = der_phi_j 

288 

289 else: 

290 if abs(der_phi_j) <= -c2 * der_phi_0: 

291 a_star = a_j 

292 phi_star = phi_j 

293 der_phi_star = der_phi_j 

294 g_star = g_j 

295 break 

296 

297 if der_phi_j * (a_hi - a_lo) >= 0.0: 

298 a_hi = a_lo 

299 f_hi = f_lo 

300 df_hi = df_lo 

301 

302 a_lo = a_j 

303 f_lo = phi_j 

304 df_lo = der_phi_j 

305 

306 i += 1 

307 

308 if np.fabs(a_lo - a_hi) < self.eps_dx and a_lo < \ 

309 self.eps_dx: 

310 a_star = a_lo 

311 phi_star, der_phi_star, g_star = \ 

312 self.evaluate_phi_and_der_phi(x, p, a_star, wfs, *args) 

313 break 

314 

315 elif np.fabs(a_lo - a_hi) < self.eps_dx: 

316 a_star = a_lo 

317 phi_star, der_phi_star, g_star = \ 

318 self.evaluate_phi_and_der_phi(x, p, a_star, wfs, *args) 

319 break 

320 

321 if i == max_iter: 

322 if a_lo > self.eps_dx: 

323 a_star = a_lo 

324 phi_star, der_phi_star, g_star = \ 

325 self.evaluate_phi_and_der_phi(x, p, a_star, wfs, *args) 

326 else: 

327 a_star = a_hi 

328 phi_star, der_phi_star, g_star = \ 

329 self.evaluate_phi_and_der_phi(x, p, a_star, wfs, *args) 

330 break 

331 

332 return a_star, phi_star, der_phi_star, g_star 

333 

334 def init_guess(self, phi_0, der_phi_0, phi_old, der_phi_old, 

335 alpha_old=1.0): 

336 

337 # chose initial alpha 

338 if self.searchdirtype in ['quasi-newton', 'newton']: 

339 alpha_1 = 1.0 

340 else: 

341 if phi_old is not None and der_phi_old is not None: 

342 try: 

343 alpha_1 = 2.0 * (phi_0 - phi_old) / der_phi_old 

344 if alpha_1 < 0.1: 

345 if alpha_old < 0.1: 

346 alpha_1 = 10.0 

347 else: 

348 alpha_1 = alpha_old 

349 except ZeroDivisionError: 

350 alpha_1 = 1.0 

351 else: 

352 alpha_1 = 1.0 

353 

354 return alpha_1