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
« 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"""
8import numpy as np
9from gpaw.directmin.tools import minimum_cubic_interpol
12def is_descent(phi_0, phi_j, eps=1.0e-2):
13 return phi_j <= phi_0 + eps * abs(phi_0)
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 """
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
26 return descent and awc
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)
44class MaxStep:
46 def __init__(self, evaluate_phi_and_der_phi, max_step=0.2):
47 """
49 :param evaluate_phi_and_der_phi:
50 """
52 self.evaluate_phi_and_der_phi = evaluate_phi_and_der_phi
53 self.max_step = max_step
54 self.name = 'max-step'
56 def todict(self):
57 return {'name': self.name,
58 'max_step': self.max_step}
60 def step_length_update(self, x, p, wfs, *args, mode=None, **kwargs):
62 kd = kwargs['kpdescr']
64 slength = get_slength(p, wfs, mode)
65 slength = kd.comm.max_scalar(slength)
67 a_star = self.max_step / slength if slength > self.max_step else 1.0
69 phi_star, der_phi_star, g_star = self.evaluate_phi_and_der_phi(
70 x, p, a_star, wfs, *args)
72 return a_star, phi_star, der_phi_star, g_star
75class StrongWolfeConditions(MaxStep):
76 """
77 From a book of Jorge Nocedal and Stephen J. Wright 'Numerical
78 Optimization' Second Edition, 2006 (p. 56)
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:
83 f(x_k + a_k p_k) <= f(x_k) + c_1 a_k grad f_k cdot p_k,
85 |grad f(x_k + a_k p_k) cdot p_k | <= c_2 |grad f_k cdot p_k|,
87 or descent conditions and approximate Wolfe conditions from
89 William W. Hager and Hongchao Zhang
90 SIAM J. optim., 16(1), 170-192.
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 """
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
115 this class works fine, but these parameters eps_dx, eps_df
116 might not be needed
117 """
119 super().__init__(
120 evaluate_phi_and_der_phi)
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'
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}
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
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)
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
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)
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
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
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
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
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
216 if np.abs(alpha_max - alpha[i]) < 1.0e-6:
217 alpha_max = 1.5 * alpha[i]
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)
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
232 a_new = minimum_cubic_interpol(
233 alpha[i], alpha_max, phi_i, phi_max,
234 der_phi_i, der_phi_max)
236 alpha.append(a_new)
237 phi_i_1 = phi_i
238 der_phi_i_1 = der_phi_i
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
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
252 i += 1
254 return a_star, phi_star, der_phi_star, g_star
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):
259 max_iter = self.max_iter
260 i = 0
262 while True:
263 a_j = minimum_cubic_interpol(
264 a_lo, a_hi, f_lo, f_hi, df_lo, df_hi)
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
272 phi_j, der_phi_j, g_j = \
273 self.evaluate_phi_and_der_phi(x, p, a_j, wfs, *args)
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
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
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
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
302 a_lo = a_j
303 f_lo = phi_j
304 df_lo = der_phi_j
306 i += 1
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
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
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
332 return a_star, phi_star, der_phi_star, g_star
334 def init_guess(self, phi_0, der_phi_0, phi_old, der_phi_old,
335 alpha_old=1.0):
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
354 return alpha_1