100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380 | def make_step_size(
rule: str | float | Callable = "constant",
*,
f: Callable | None = None,
gradient: Callable | None = None,
**options: Any,
) -> Callable:
"""Build ``step(n, x=None, d=None, *, gradient_value=None, f_value=None)``.
Schedules accept a positive number, a callable of ``n``, or the names
``constant``, ``not_summable``, ``square_summable_not_summable``,
``geometric_series``, and ``user_defined``. Their formulas are respectively
``a``, ``a/sqrt(n)``, ``a/(b+n)``, and ``a*r**n``; ``n >= 1``.
Line searches are ``armijo``, ``wolfe``, ``strong_wolfe``, and ``exact``
(case-insensitive). They require ``f``; Armijo/Wolfe additionally require
a descent direction. A supplied
``gradient`` callback always takes precedence. Without one, ordinary
centered differences are used, or specular differences for names prefixed
with ``specular_``. Cached values must match that raw gradient and objective.
Common line-search options are ``t_0=1``, ``max_alpha=1e8``,
``max_iter=100``, and ``h=None``. Armijo/Wolfe use ``c_1=1e-4`` and
``rho=0.5``; Wolfe also uses ``c_2=0.9`` (``c_3`` is a legacy alias for
strong Wolfe). ``exact`` uses golden-section minimization on
``[0, max_alpha]`` with ``tol=1e-8`` and at most ``max_iter`` objective
evaluations, including the initial value when uncached. It is derivative
free; ``specular_exact`` is an identical alias. The name describes
numerical bounded minimization, not a global or symbolic guarantee.
Failed searches raise :class:`LineSearchError`; no unchecked fallback step
is returned. Specular slopes do not imply classical smooth Wolfe guarantees.
"""
name = rule.strip().lower() if isinstance(rule, str) else None
schedules = {
"constant", "not_summable", "square_summable_not_summable",
"geometric_series", "user_defined",
}
if name is None or name in schedules:
if callable(rule):
_unknown(options, set())
schedule = rule
elif name is None:
_unknown(options, set())
value = _positive(rule, "step size")
schedule = lambda n: value
elif name == "user_defined":
_unknown(options, {"user_defined_rule"})
schedule = options.get("user_defined_rule")
if not callable(schedule):
raise TypeError("user_defined_rule must be callable")
else:
allowed = {"a"}
if name == "square_summable_not_summable":
allowed.add("b")
if name == "geometric_series":
allowed.add("r")
_unknown(options, allowed)
a = _positive(options.get("a", 1.0), "a")
if name == "constant":
schedule = lambda n: a
elif name == "not_summable":
schedule = lambda n: a / math.sqrt(n)
elif name == "square_summable_not_summable":
b = _positive(options.get("b", 1.0), "b", zero=True)
schedule = lambda n: a / (b + n)
else:
r = _positive(options.get("r", 0.5), "r")
if r >= 1:
raise ValueError("r must be strictly between zero and one")
schedule = lambda n: a * r**n
def scheduled(n, x=None, d=None, *, gradient_value=None, f_value=None):
return _positive(schedule(_integer(n, "n")), "scheduled step")
return scheduled
specular = name.startswith("specular_")
base = name.removeprefix("specular_")
if base not in {"armijo", "wolfe", "strong_wolfe", "exact"}:
raise ValueError(f"Unknown step-size rule: {rule!r}")
if not callable(f):
raise TypeError("Line searches require a callable f")
if gradient is not None and not callable(gradient):
raise TypeError("gradient must be callable")
allowed = {"t_0", "max_alpha", "max_iter", "h"}
if base == "exact":
allowed.add("tol")
else:
allowed.update({"c_1", "rho"})
if base != "armijo":
allowed.add("c_2")
if base == "strong_wolfe":
allowed.add("c_3")
_unknown(options, allowed)
initial = _positive(options.get("t_0", 1.0), "t_0")
bound = _positive(options.get("max_alpha", 1e8), "max_alpha")
budget = _integer(options.get("max_iter", 100), "max_iter")
h = options.get("h")
if h is not None:
h = _positive(h, "h")
c1 = _positive(options.get("c_1", 1e-4), "c_1")
rho = _positive(options.get("rho", 0.5), "rho")
if c1 >= 1 or rho >= 1:
raise ValueError("c_1 and rho must be strictly between zero and one")
if "c_2" in options and "c_3" in options:
raise TypeError("Use only one of c_2 and its legacy alias c_3")
c2 = _positive(options.get("c_2", options.get("c_3", 0.9)), "c_2")
if base in {"wolfe", "strong_wolfe"} and not c1 < c2 < 1:
raise ValueError("Wolfe parameters must satisfy 0 < c_1 < c_2 < 1")
tolerance = _positive(options.get("tol", 1e-8), "tol")
if gradient is None:
if specular:
def raw_gradient(point):
calc = _specular_derivative if np.asarray(point).ndim == 0 else _specular_gradient
return calc(f, point, h=h)
else:
def raw_gradient(point):
return _centered_gradient(f, point, h)
else:
raw_gradient = gradient
def search(n, x=None, d=None, *, gradient_value=None, f_value=None):
_integer(n, "n")
point, direction = _point(x, "x"), _point(d, "d")
if point.shape != direction.shape:
raise ValueError("x and d must have the same shape")
if not np.any(direction):
raise LineSearchError("Line search requires a nonzero descent direction")
def slope(at, cached=None):
try:
value = raw_gradient(_argument(at)) if cached is None else cached
except ArithmeticError:
return math.nan
array = np.asarray(value)
if array.shape != point.shape or array.dtype.kind not in "iuf":
raise ValueError("Raw gradient must be real and have the same shape as x")
if not np.all(np.isfinite(array)):
return math.nan
with np.errstate(over="ignore", invalid="ignore"):
return float(np.sum(array.astype(np.longdouble) * direction.astype(np.longdouble)))
try:
phi0 = _scalar(f(_argument(point)) if f_value is None else f_value, "objective value")
except ArithmeticError as exc:
raise LineSearchError("Initial objective value must be finite") from exc
if not math.isfinite(phi0):
raise LineSearchError("Initial objective value must be finite")
if base != "exact":
slope0 = slope(point, gradient_value)
if not math.isfinite(slope0):
raise LineSearchError("Initial slope must be finite")
if slope0 >= 0:
raise LineSearchError("Line search requires a descent direction (negative raw slope)")
def trial(alpha):
try:
at = np.array([
math.fma(alpha, float(di), float(xi))
for xi, di in zip(point.reshape(-1), direction.reshape(-1))
]).reshape(point.shape)
except OverflowError:
return np.full_like(point, math.inf), math.inf
if not np.all(np.isfinite(at)) or np.array_equal(at, point):
return at, math.inf
try:
value = _scalar(f(_argument(at)), "objective value")
except ArithmeticError:
return at, math.inf
return at, value if math.isfinite(value) else math.inf
def armijo(alpha, phi):
threshold = np.longdouble(phi0) + np.longdouble(c1) * alpha * slope0
return math.isfinite(phi) and phi <= threshold
def curvature(value):
if not math.isfinite(value):
return False
if base == "strong_wolfe":
return abs(value) <= -c2 * slope0
return value >= c2 * slope0
if base == "exact":
evaluations = 0 if f_value is not None else 1
best_alpha, best_value = 0.0, phi0
def evaluate(alpha):
nonlocal evaluations, best_alpha, best_value
if evaluations >= budget:
raise LineSearchError("Bounded minimization evaluation budget exhausted")
evaluations += 1
_, value = trial(alpha)
if math.isfinite(value) and (value < best_value or (value == best_value and best_alpha == 0)):
best_alpha, best_value = alpha, value
return value
left, right = 0.0, bound
evaluate(right)
if initial < bound:
evaluate(initial)
ratio = (math.sqrt(5.0) - 1.0) / 2.0
c, e = right - ratio * (right - left), left + ratio * (right - left)
fc, fe = evaluate(c), evaluate(e)
while right - left > tolerance * max(1.0, abs(left + (right - left) / 2.0)):
previous = (left, right)
if fc <= fe:
right, e, fe = e, c, fc
c = right - ratio * (right - left)
fc = evaluate(c)
else:
left, c, fc = c, e, fe
e = left + ratio * (right - left)
fe = evaluate(e)
if (left, right) == previous:
raise LineSearchError("Bounded minimization stagnated at floating-point precision")
if best_alpha <= 0:
raise LineSearchError("Bounded minimization found no positive improving step")
return best_alpha
alpha = min(initial, bound)
if base == "armijo":
for _ in range(budget):
at, phi = trial(alpha)
if np.array_equal(at, point):
raise LineSearchError("Armijo step stagnated at floating-point precision")
if armijo(alpha, phi):
return alpha
next_alpha = alpha * rho
if next_alpha <= 0 or next_alpha == alpha:
break
alpha = next_alpha
raise LineSearchError("Armijo search could not accept a step")
previous_alpha, previous_phi = 0.0, phi0
low = high = low_phi = None
for iteration in range(budget):
if low is not None:
alpha = low + (high - low) / 2.0
if alpha == low or alpha == high or alpha <= 0:
raise LineSearchError("Wolfe bracket stagnated at floating-point precision")
at, phi = trial(alpha)
if np.array_equal(at, point):
next_alpha = min(alpha / rho, bound)
if low is not None or next_alpha <= alpha:
raise LineSearchError("Wolfe step stagnated at floating-point precision")
alpha = next_alpha
continue
if not armijo(alpha, phi) or (low is not None and phi >= low_phi):
if low is None:
low, low_phi = previous_alpha, previous_phi
high = alpha
continue
if low is None and iteration > 0 and phi >= previous_phi:
low, low_phi, high = previous_alpha, previous_phi, alpha
continue
value = slope(at)
if curvature(value):
return alpha
if not math.isfinite(value):
if low is None:
low, low_phi = previous_alpha, previous_phi
high = alpha
continue
if low is not None:
if (value >= 0) == (high >= low):
high = low
low, low_phi = alpha, phi
elif value >= 0:
low, low_phi, high = alpha, phi, previous_alpha
else:
previous_alpha, previous_phi = alpha, phi
next_alpha = min(alpha / rho, bound)
if next_alpha <= alpha:
raise LineSearchError("Wolfe search reached max_alpha without acceptable curvature")
alpha = next_alpha
raise LineSearchError("Wolfe search exhausted its iteration budget")
return search
|
