"""The pure similarity root finder: one similarity function, one step rule. direction : the side of the bracket that contains the sign change length : (1 - S(z)) * (distance to the far end), z = f(x) / scale Everything about the method is in S and in where `scale` comes from. scale CHOSEN (a fixed reference value) -> linear convergence at a rate that is a property of the problem, not of the method, and therefore fails on some problems (stage 0b, 5-8 of 8) scale DERIVED (the local slope of f, read off the two most recent evaluations) -> the step becomes |f|/|f'| and the convergence is superlinear (stage 0d/0e) The derivation, from the contraction law measured in stage 0c: e_{n+1} = e_n * |1 - c*|f'(r)|*D/scale| , c = slope of (1 - S) at 0 rate = 0 <=> c = scale / (|f'(r)|*D) so the one quantity the pure rule needs is the LOCAL SCALE OF THE VALUE AXIS. Choosing it is what makes the method brittle; reading it off the data is what makes it parameter free - the same move the SIM control paper makes when it normalises the input space from the data instead of choosing membership widths by hand. Safeguards ---------- An interpolation-only rule has no worst-case bound. Two parameter-free safeguards are provided and are measured against each other: "bracket" bisect whenever a similarity step failed to halve the BRACKET. Sound, but expensive: a secant step routinely collapses the error while barely moving the proven bracket, so the guard fires on steps that were in fact excellent (measured: 150 vs 86 evaluations over the stage-0 problem set). "step" bisect whenever the proposed step is not at most half the step taken two evaluations ago. Step lengths are then forced to decrease geometrically, and a converging secant sequence - whose steps shrink far faster than by half - never triggers it. This is the criterion Brent's method uses internally, stated without his tuning constants. It does NOT on its own bound the worst case, and an earlier version of this docstring claimed that it did. A sequence that creeps towards a non-root from one side with steps that halve every time has small steps by construction, so the guard stays silent, the sign never changes, and the proven bracket never moves. Cheap, but not sound alone. "both" the step guard, with the bracket guard underneath it as the backstop that actually carries the worst-case bound. The step guard does the work on ordinary problems, so the bracket guard almost never fires and almost never costs anything. "project" do not REPLACE the similarity step at all - PROJECT it onto the largest interval still consistent with the bisection bound: allowed_j = 2^(n_max - j) * (xtol/2) x = clamp(x_sim, hi - allowed_j, lo + allowed_j) with n_max = 2*ceil(log2(w0 / xtol)) fixed by the initial bracket. The next bracket is at most max(x - lo, hi - x), so the clamp gives w_{j+1} <= allowed_j directly, and induction bounds the run by n_max steps in exact arithmetic - one more when the bracket lands exactly on the boundary in floating point. The guarantee is therefore 2*ceil(log2(w0/xtol)) + 3 evaluations: twice bisection plus one step. Zero slack (n_max = ceil(log2(w0/xtol)), a worst case of exactly bisection) is degenerate and was measured as such: it demands a halving at every single step, leaving the projection no freedom, and the guard behaves as bisection. The three guards above discard a computed step; this one keeps as much of it as the bound allows, which is why it does not flatten the order. The construction is ITP's projection stage. ITP's three constants (kappa1, kappa2, n0) all live in its TRUNCATION stage, which exists to repair regula falsi; the similarity step needs no repair, so the truncation - and with it every constant - is dropped. "safe" the step guard for cheapness, with the projection carrying the bound. This is the configuration the paper proposes: a proven worst case AND the full order, which the three earlier guards could not deliver together. The only constant in any of them is 1/2, which is the definition of bisection rather than a tuned value. """ from __future__ import annotations import math from .baselines import Counter, _Budget # -------------------------------------------------------------------------- # [0,1] similarity functions of a normalised residual z: S(0) = 1 (identical # to the expectation), S -> 0 far from it. # -------------------------------------------------------------------------- PURE_SIMS = { # the ratio form S = 1 - |f|/scale. Note that the symmetric ratio # min(|x|,|y|)/max(|x|,|y|) is identically 0 when one argument is the # target, so this is the only usable reading of "ratio similarity to 0". "ratio": lambda z: max(0.0, 1.0 - abs(z)), "sqrt": lambda z: max(0.0, 1.0 - math.sqrt(abs(z))), "cbrt": lambda z: max(0.0, 1.0 - abs(z) ** (1.0 / 3.0)), "rational": lambda z: 1.0 / (1.0 + abs(z)), "ratio_soft": lambda z: 1.0 / (1.0 + abs(z)) ** 2, "exp_abs": lambda z: math.exp(-abs(z)), # similarities that approach 1 quadratically: c = 0, predicted rate 1, # kept because the law has to be falsifiable in both directions "rat_sq": lambda z: 1.0 / (1.0 + z * z), "exp_sq": lambda z: math.exp(-z * z), "cosine": lambda z: 1.0 / math.sqrt(1.0 + z * z), } GUARDS = ("none", "bracket", "step", "both", "project", "safe", "credit") class _Safeguard: """Parameter-free worst-case protection. See the module docstring.""" __slots__ = ("mode", "guard_width", "pending", "step_2", "step_1", "n_half", "eps", "j", "allow", "credit") def __init__(self, mode, width, xtol=None, budget_factor=2.0, budget_extra=None): if mode not in GUARDS: raise ValueError(f"unknown guard {mode!r}, expected one of {GUARDS}") self.mode = mode self.guard_width = width self.pending = False self.step_2 = width # step two evaluations ago self.step_1 = width # step one evaluation ago # projection state: the bisection budget, fixed by the initial bracket self.eps = 0.5 * xtol if xtol else 0.0 self.j = 0 # set for every mode: the step loop decrements it unconditionally, and # leaving it to the projection branch made every other guard raise self.allow = 0.0 self.n_half = 0 self.credit = 1.0 # one bit of endowment; see project_credit if mode in ("project", "safe"): if not xtol or width <= 0.0: raise ValueError("the projection guard needs xtol and a " "positive initial width") half = math.ceil(math.log2(width / xtol)) if width > xtol else 0 # The budget is TWICE the bisection budget. Zero slack (budget = # half) is degenerate: a worst case of exactly bisection demands a # halving at every single step, so the projection leaves no # freedom at all - measured, r_0 was 0.05 of a unit bracket and # the guard behaved as bisection. Doubling is not a tuned # constant: it is the same factor 2 already claimed as the bound, # and it buys total freedom until the budget is nearly spent. # Two ways to give the projection room, and the difference # matters for whether this method contains a constant at all. # # budget_extra = k -> n_max = half + k. The requirement is # only r_0 >= 0, i.e. n_max >= half; at n_max = half there # is no room and the guard degenerates to bisection. The # MINIMAL non-degenerate choice is therefore k = 1 - not a # tuned value but the smallest integer for which the # construction works at all. # budget_factor = c -> n_max = ceil(c * half). This is what # the paper used, at c = 2, and the justification given for # it ("the factor already in the bound") was circular: the # bound is 2x because 2 was chosen. Kept as an ablation # knob so the circularity can be measured rather than # argued about. if budget_extra is not None: self.n_half = half + int(budget_extra) else: self.n_half = max(half, int(math.ceil(budget_factor * half))) # The allowance 2^(n_max - j).eps is halved once per step, so it # is carried forward rather than recomputed: an exponentiation in # the inner loop was costing more than the projection itself. # Halving is exact in binary floating point, so this changes no # iterate and no evaluation count - only the arithmetic. try: self.allow = 2.0 ** self.n_half * self.eps except OverflowError: self.allow = math.inf def project_credit(self, x_sim, lo, hi): """The projection with the budget EARNED rather than scheduled. The scheduled version fixes the allowance in advance from the initial width, so the multiplier in n_max = c*half is a free constant, and the argument that c = 2 is "the factor already in the bound" is circular: the bound is 2x because 2 was chosen. Measured, c = 1.25 is 8% better, so the constant is not neutral either. Here the method pays its own way. Let s_j = log2(w_j / w_{j+1}) be the shrink a step delivers in bits; bisection delivers exactly 1. Keep C <- C + s_j - 1, C >= 0, the surplus banked so far, and allow a step to under-deliver only by what has been banked: the width after the step may be at most 2^(C-1) * w. Since C never goes negative, the shrinks over j steps sum to at least j - C_0, so the run costs at most half + C_0 steps. C_0 = 1 - one bit of endowment - is the smallest integer for which the rule is not degenerate: at C_0 = 0 the allowance is exactly w/2 on the first step, which forces the midpoint and makes the method bisection. So the guarantee is bisection plus one step, and the constant is the promise, not a tuning knob. """ w = hi - lo # The allowance is w * 2^(C-1), so it exceeds w exactly when C >= 1, # and then it restricts nothing. Testing that FIRST removes the # exponentiation from the common path and, more importantly, removes # a crash: C accumulates the log of every shrink and is unbounded, so # 2**(C-1) overflows a double at C ~ 1025. That is reachable - sixty # nine-decade steps reach C = 1012, and a bracket of width 1e300 with # a tolerance of 1e-8 has about 1023 bits of shrink available - and # the OverflowError would escape, since only the budget exception is # caught here. if self.credit >= 1.0: return x_sim allowed = w * (2.0 ** (self.credit - 1.0)) left, right = hi - allowed, lo + allowed if left > right: return 0.5 * (lo + hi) return min(max(x_sim, left), right) def note_shrink(self, w_before, w_after): """Bank what the step delivered beyond a halving, or draw on the bank for what it fell short.""" if w_before <= 0.0 or w_after <= 0.0: return s = math.log2(w_before / w_after) self.credit = max(0.0, self.credit + s - 1.0) def project(self, x_sim, lo, hi): """Pull the similarity step back to the edge of the interval that the bisection budget still allows. Returns the admissible iterate. The budget says the next bracket may be at most `allowed = 2^(n_max - j) * eps`. Whichever side the sign change turns out to be on, that bracket is at most max(x - lo, hi - x), so the admissible set is exactly [hi - allowed, lo + allowed]. Writing it this way rather than as a radius r = allowed - (hi-lo)/2 around the midpoint matters numerically: r is the difference of two nearly equal numbers, and the cancellation made the computed radius slightly too large. Measured, that admitted points which broke the invariant by a relative 8.9e-5 and cost exactly one extra evaluation on eleven Chandrupatla problems. The clamped form below performs no cancelling subtraction and cannot overshoot.""" allowed = self.allow if not math.isfinite(allowed): return x_sim left, right = hi - allowed, lo + allowed if left > right: # the budget no longer permits anything but a halving return 0.5 * (lo + hi) return min(max(x_sim, left), right) def bisect_before_step(self, width): """Called before the similarity step is even computed. Only the bracket guard can decide at this point.""" if self.mode not in ("bracket", "both"): return False if width <= 0.5 * self.guard_width: self.guard_width = width # The similarity step that is about to be taken is itself the # pending one. Clearing the flag here instead granted TWO # consecutive similarity steps after every halving, which makes # the guarantee one bisection in three evaluations (3x bisection) # rather than one in two (2x) - measured at 139 evaluations # against a 2x bound of 114 on Chandrupatla fun4.3. self.pending = True return False if self.pending: self.pending = False return True self.pending = True return False def reject_step(self, step): """Called with the proposed displacement. Only the step guard can decide at this point.""" if self.mode not in ("step", "both", "safe"): return False return abs(step) > 0.5 * self.step_2 def note_step(self, step): self.step_2 = self.step_1 self.step_1 = abs(step) def _slope3(x2, f2, x1, f1, x0, f0): """Derivative at x2 of the parabola through the three most recent points, as Newton divided differences: f'(x2) ~ f[x1,x2] + f[x0,x1,x2] * (x2 - x1) Falls back to the two-point slope when a point is repeated or missing, so the first iteration of a run behaves exactly like scale="derived".""" if x2 == x1: return 0.0 d12 = (f2 - f1) / (x2 - x1) # An overflowing difference must report "no usable slope" (0.0), which # sends pure_root to its bisection fallback. Returning the infinity # instead would make the scale infinite, hence z = 0, hence "perfectly # similar to the target", hence a zero step - a stall dressed up as # success. if not math.isfinite(d12): return 0.0 if x0 is None or x0 == x1 or x0 == x2: return d12 d01 = (f1 - f0) / (x1 - x0) d012 = (d12 - d01) / (x2 - x0) s = d12 + d012 * (x2 - x1) return s if math.isfinite(s) else d12 def _local_order(x2, f2, x1, f1, x0, f0): """The local order m of f, from three evaluations and nothing else. Redo the contraction law WITHOUT assuming |f| ~ |f'(r)|.|x-r|. Take instead |f| ~ K.e^m near the current point, m = 1 being a simple root and m > 1 a locally flat approach. Then for a power similarity S = 1 - z^p, z = |f| / (|f'|.D) = e / (m.D) (K cancels) step = z^p . D = D^(1-p) . (e/m)^p and the step lands on the root when (e/D)^(p-1) = m^p, i.e. p = L / (L - ln m), L = ln(e/D) < 0. m = 1 gives p = 1, the earlier law. m > 1 drives p below 1 on its own - at e/D = 0.1, m = 10 it gives p = 0.50, and at e/D = 0.5 it gives 0.23, which is the neighbourhood of the cube-root similarity that stage 7 found winning on the wide-bracket families. Substituting that p back collapses the whole expression to step = e = m.|f|/|f'|, so the exponent is a reparameterisation and the substance is the multiplier m. m is measurable without knowing the root: with theta = f.f''/(f')^2, f ~ K e^m, f' ~ K m e^(m-1), f'' ~ K m(m-1) e^(m-2) => theta = m(m-1)/m^2 = 1 - 1/m => m = 1/(1 - theta). Both derivatives are divided differences of the three most recent points, so nothing is chosen or fitted. Returns 1.0 whenever the power-law model does not apply - theta outside (0,1) means the data does not describe a flat approach, and the honest response is to decline the correction rather than to clamp it into range. """ if x0 is None or len({x0, x1, x2}) < 3: return 1.0 d01 = (f1 - f0) / (x1 - x0) d12 = (f2 - f1) / (x2 - x1) d012 = (d12 - d01) / (x2 - x0) fp = d12 + d012 * (x2 - x1) # f'(x2) on the parabola fpp = 2.0 * d012 # f''(x2) on the parabola if fp == 0.0 or not math.isfinite(fp) or not math.isfinite(fpp): return 1.0 theta = f2 * fpp / (fp * fp) if not math.isfinite(theta) or not (0.0 < theta < 1.0): return 1.0 m = 1.0 / (1.0 - theta) return m if math.isfinite(m) and m >= 1.0 else 1.0 def pure_root(f, a, b, xtol, maxfev, similarity="ratio", scale="derived", guard="step", stats=None, budget_factor=2.0, budget_extra=None, endowment=1.0): """The pure similarity root finder. similarity : key into PURE_SIMS scale : "derived" - local slope of f from the two most recent evaluations (parameter free, order ~1.618) "derived3_m" - REFUTED as an improvement. The three-point scale divided by the local order m of f, which is what the contraction law gives once the linearity assumption is dropped. The derivation is sound and the order is preserved; the estimator of m is not sharp enough to pay for itself. See the note at the use site. "derived_rho", "derived3_rho" - REFUTED. The same, with the scale multiplied by the fraction of its own promise the previous step kept. Retained for reproducibility only; see the note at the use site. Do not adopt. "derived3" - local slope from the three most recent evaluations, i.e. the derivative at x of the parabola through them. Still parameter free - a divided difference is a computation, not a fitted coefficient - but the slope error drops from O(h) to O(h^2), which lifts the order from 1.618 towards 1.84. This is the direct consequence of the stage-4 finding that the whole deficit against Brent is the endgame order, and nothing else. "start" - |f| at the starting point (a chosen reference) "bracket" - largest |f| at the current bracket ends guard : one of GUARDS endowment : ABLATION KNOB for guard="credit", not a setting. The method is defined at 1 bit, and 1 is not a chosen value: the guarantee is bisection + ceil(C0) steps, so every C0 in (0,1] promises the same thing - one step - and among those 1 is the most permissive and therefore the cheapest. Measured over the 199 problems: 0.1 -> 8395, 0.5 -> 7657, 1 -> 6588, all with the same one-step guarantee; 2 -> 6404 but two steps, 5 -> 5621 but five. Nothing is free to choose once the promise is fixed, and the promise is the specification. stats : optional dict; on return it carries 'forced' = how many steps the safeguard turned into bisections. Comparable with the same key from baselines.brent, so the two methods' fallback rates can be put side by side. Returns (x, nfev, converged). """ if similarity not in PURE_SIMS: raise ValueError(f"unknown similarity {similarity!r}") if scale not in ("derived", "derived3", "derived_rho", "derived3_rho", "derived3_m", "start", "bracket"): raise ValueError(f"unknown scale {scale!r}") S = PURE_SIMS[similarity] g = Counter(f, maxfev) x = 0.5 * (a + b) try: fa, fb = g(a), g(b) if fa == 0.0: return a, g.n, True if fb == 0.0: return b, g.n, True if fa * fb > 0.0: raise ValueError("f(a) and f(b) must straddle zero") A, FA, B, FB = a, fa, b, fb # the proven bracket if abs(fa) <= abs(fb): x, fx, x_prev, f_prev = a, fa, b, fb else: x, fx, x_prev, f_prev = b, fb, a, fa # the third point does not exist yet; _slope3 falls back to two points x_prev2, f_prev2 = None, None start_scale = abs(fx) sg = _Safeguard(guard, abs(B - A), xtol, budget_factor, budget_extra) sg.credit = endowment forced = 0 if stats is not None: stats['forced'] = 0 while abs(B - A) > xtol: if sg.bisect_before_step(abs(B - A)): x_new = 0.5 * (A + B) forced += 1 # keep the step history honest under guard="both": a bisection # forced by the bracket backstop is still a step that was taken sg.note_step(x_new - x) else: far = B if abs(B - x) >= abs(A - x) else A if scale in ("derived", "derived_rho"): slope = ((fx - f_prev) / (x - x_prev)) if x != x_prev else 0.0 sc = abs(slope) * abs(far - x) elif scale in ("derived3", "derived3_rho", "derived3_m"): slope = _slope3(x, fx, x_prev, f_prev, x_prev2, f_prev2) sc = abs(slope) * abs(far - x) if scale == "derived3_m": # REFUTED as an improvement - kept for reproducibility. # The DERIVATION holds: dropping the linearity # assumption from the contraction law does produce the # march exponent, and it collapses to the multiplier m. # The ESTIMATOR does not: theta from three points # cannot tell a genuinely flat approach from curvature # noise. Measured, it reports m > 1 on 27.9% of steps # across Alefeld-Potra-Shi, where almost every root is # simple and the truth is m = 1, with excursions to # m = 34.8; on the flat fun3/4/5 families it fires on # 35.8%. A separation of 27.9 against 35.8 is no # separation, and the false positives cost 273 # evaluations: 3079 against 2806 on APS, 6304 against # 6026 over both sets. # # It does pass the two tests rho failed - the order # stays at 2.001 and the eight standard problems # improve to 82 - which is exactly why it must be # labelled: a small hand-picked set makes it look like # a win. # # step = m.|f|/|f'| is step = e for a locally flat # approach; dividing the scale by m is the same thing # written in the similarity's own terms. m >= 1, so # like rho this can only lengthen a step - but unlike # rho it is read from the CURVATURE, which stays # meaningful as |f| goes to zero, and it returns to 1 # at a simple root, leaving the endgame untouched. sc /= _local_order(x, fx, x_prev, f_prev, x_prev2, f_prev2) elif scale == "start": sc = start_scale else: sc = max(abs(FA), abs(FB)) if scale.endswith("_rho"): # REFUTED - kept only so the experiment is reproducible. # Measured (stage 8): the fun3/4/5 families it was built # for do not improve (1364 against 1367, where cbrt gets # 677); the eight standard problems get twice as # expensive (84 -> 161); and the observed order collapses # from 2.000 to 0.904. # # Why the derivation was wrong: rho equates "how much # value the step removed" with "how trustworthy the local # model is". Near the root |f_n| and |f_{n-1}| are both # negligible and their ratio is noise, so rho thrashes in # exactly the place where the correction has to be inert. # The alternative credibility - the ratio of the assumed # slope to the realised one - is closed too: the realised # slope is what the derived scale already computes, so # correcting by it would double-count. # # The phase detector. A step of |f|/|f'| is a promise to # land on the root; rho is the fraction of that promise # the previous step actually kept: # # rho = 1 - |f_n| / |f_{n-1}| # # Near the root the local model is right and rho -> 1, so # the correction vanishes and the order is untouched. Far # from it - the approach phase on a very wide bracket - # the model under-delivers, rho is small, and the step # must be 1/rho longer to remove what is left. Dividing # the scale by rho does exactly that. # # rho <= 1 always, so the correction can only LENGTHEN the # step, never shorten it below the Newton step. rho <= 0 # means the value did not improve at all: there is no # credibility to read, so the correction is dropped rather # than guessed, and the bracket projection keeps the step # safe either way. if f_prev != 0.0: rho = 1.0 - abs(fx) / abs(f_prev) if 0.0 < rho <= 1.0: sc *= rho if sc == 0.0 or not math.isfinite(sc): x_new = 0.5 * (A + B) else: step = (1.0 - S(abs(fx) / sc)) * (far - x) if sg.reject_step(step): x_new = 0.5 * (A + B) step = x_new - x forced += 1 else: x_new = x + step if guard in ("project", "safe", "credit"): projected = (sg.project_credit if guard == "credit" else sg.project)(x_new, min(A, B), max(A, B)) if projected != x_new: forced += 1 x_new = projected step = x_new - x sg.note_step(step) if not math.isfinite(x_new) or not (min(A, B) < x_new < max(A, B)): x_new = 0.5 * (A + B) if x_new == x or x_new == x_prev: x_new = 0.5 * (A + B) if x_new == x or x_new == x_prev: # the bracket can no longer be split in floating point if stats is not None: stats['forced'] = forced return x_new, g.n, abs(B - A) <= xtol w_before = abs(B - A) fnew = g(x_new) sg.j += 1 # one step of the bisection budget is spent sg.allow *= 0.5 x_prev2, f_prev2 = x_prev, f_prev x_prev, f_prev = x, fx x, fx = x_new, fnew if fnew == 0.0: if stats is not None: stats['forced'] = forced return x_new, g.n, True if FA * fnew < 0.0: B, FB = x_new, fnew else: A, FA = x_new, fnew if guard == "credit": sg.note_shrink(w_before, abs(B - A)) if stats is not None: stats['forced'] = forced return 0.5 * (A + B), g.n, True except _Budget: if stats is not None: stats['forced'] = forced return x, g.n, False