diff --git a/article/sim_numerics.pdf b/article/sim_numerics.pdf index 8fe3ae2..616b12a 100644 Binary files a/article/sim_numerics.pdf and b/article/sim_numerics.pdf differ diff --git a/article/sim_numerics.tex b/article/sim_numerics.tex index c2d5ff0..ebddead 100644 --- a/article/sim_numerics.tex +++ b/article/sim_numerics.tex @@ -628,16 +628,35 @@ It is not the guarantee that caps us. With the safeguard removed entirely the same problems cost $46$, $57$, $61$ and $68$ against bisection's $57$, $77$, $45$ and $45$ --- better on the wide brackets, \emph{worse than bisection} on \texttt{fun3.1} and \texttt{fun5.1} --- and still far from Chandrupatla's $24$ and -$44$. The cap is the step rule itself: $|f|/|f'|$ is an absolute length, and an -absolute length cannot march across decades. That is the same wall the fixed -fraction (\texttt{cbrt}) and the logarithmic change of variable each hit from a -different side, both buying the wide brackets at the price of the canonical set. +$44$. -So the limit is stated exactly: the placement of the similarity reading improves -the method up to the point where the \emph{form} of the derived step becomes the -binding constraint, and on brackets spanning many decades it does. Removing that -constraint means replacing the step the contraction law derives, which is a -different method rather than a better reading of this one. +Nor is it that the step fails to be relative to the bracket. With the ratio +similarity $1-S=\min(1,|f|/(|f'|D))$, so the step is exactly +$\min(D,|f|/|f'|)$: a \emph{fraction} of the distance to the far end, and the +similarity is that fraction. Nor does the iterate creep along one side without +moving the proven bracket: measured probe by probe, $91.5\%$ of our evaluations +on the $2\times10^4$-wide bracket do shrink the bracket, slightly \emph{more} +often than Chandrupatla's $87.5\%$. + +What differs is how much each probe removes. Ours shrinks the bracket by +$0.93$ halvings per evaluation there, Chandrupatla's by $2.26$ --- a factor of +about two against a factor of about five. On a wide bracket of a convex function +the Newton step systematically undershoots, because a linear model of the value +axis cannot see the curvature that separates the iterate from the root, so our +probe lands well short and the new bracket is about half the old one; a model +that carries the curvature --- inverse quadratic interpolation --- lands far +closer and removes about four fifths. + +So the limit is stated exactly, and it is the \emph{order of the model the +similarity is read against}, not the placement of the reading and not the +guarantee. Three parameter-free ways round it have now been measured and all +three fail in the same way, buying the wide brackets at the price of the +canonical set: a fixed fraction of the bracket (\texttt{cbrt}), a logarithmic +change of variable, and a curvature-corrected scale (\texttt{derived3\_m}, +whose derivation holds but whose estimator of the local order fires on $27.9\%$ +of steps where the truth is a simple root). The fourth way is to model the +curvature directly, which is inverse quadratic interpolation --- a different +method rather than another reading of this one. On that set our shape reading ($2159$ re-scored) is \emph{worse} than bisection ($2136$); only the dissimilarity reading ($1740$) is clearly ahead, which is the diff --git a/diag_sign_flips.py b/diag_sign_flips.py new file mode 100644 index 0000000..fd33727 --- /dev/null +++ b/diag_sign_flips.py @@ -0,0 +1,104 @@ +"""Correcting the diagnosis: the step IS a fraction of the bracket, so what stalls? + +With S = ratio, 1 - S = min(1, |f|/(|slope| D)), hence step = min(D, |f|/|slope|): the step is +a FRACTION of the distance to the far end and the fraction is the similarity. So "an absolute +length cannot march" was wrong. + +The candidate mechanism instead: when the Newton length is much shorter than the bracket the +iterate creeps from ONE SIDE, the sign never changes, and the PROVEN bracket does not shrink +at all — only the guard's bisections shrink it. Chandrupatla, bisecting on 72% of steps and +placing its accepted probes by inverse quadratic interpolation, may instead CROSS the root +and collapse the bracket in one step. + +Measured here, per problem: the share of evaluations that flip the sign of f (i.e. actually +shrink the proven bracket), and the geometric bracket reduction per evaluation. + + python -u diag_sign_flips.py +""" + +from __future__ import annotations + +import math + +from scipy.optimize._tstutils import _CHANDRUPATLA_TESTS_DICTS + +from simroot.baselines import chandrupatla, chandrupatla_scipy +from simroot.pure import pure_root +from sweep_full import ACC, MAXFEV, XTOL + + +def traced_bracket(f, a, b): + """Replay a run's probes and count how many flipped the sign, plus the bracket history.""" + xs = [] + + def g(x): + v = f(x) + xs.append((x, v)) + return v + return g, xs + + +def replay(xs, a, b, fa, fb): + """Walk the probes in order, maintaining the proven bracket; count sign flips.""" + A, FA, B, FB = a, fa, b, fb + flips = 0 + widths = [abs(B - A)] + for x, v in xs: + if not (min(A, B) < x < max(A, B)): + continue + if FA * v < 0.0: + B, FB = x, v + flips += 1 if (FB * FA < 0.0) else 0 + else: + A, FA = x, v + widths.append(abs(B - A)) + return flips, widths + + +def run(method, f, a, b, root): + g, xs = traced_bracket(f, a, b) + try: + x, n, conv = method(g, a, b, XTOL, MAXFEV) + except Exception: + return None + x = float(x) + if not (conv and math.isfinite(x) and (abs(x - root) <= ACC or f(x) == 0.0)): + return None + fa, fb = f(a), f(b) + _flips, widths = replay(xs[2:], a, b, fa, fb) + # how often did a probe actually reduce the proven bracket, and by how much on average + reduced = sum(1 for i in range(1, len(widths)) if widths[i] < widths[i - 1] * 0.999) + halvings = (math.log2(widths[0] / max(widths[-1], 1e-300)) / max(len(xs), 1)) + return len(xs), reduced, 100.0 * reduced / max(len(xs), 1), halvings + + +def main(): + print(f"{'problem':<9} {'method':<14} {'evals':>6} {'probes that shrank':>19} " + f"{'share':>7} {'halvings/eval':>14}") + for d in _CHANDRUPATLA_TESTS_DICTS: + pid = str(d["ID"]) + if pid not in ("fun1.3", "fun1.4", "fun1.5", "fun3.1", "fun5.1"): + continue + f = d["f"] + a, b = (float(v) for v in d["bracket"]) + root = float(d["root"]) + ours = run(lambda ff, aa, bb, xt, mf: pure_root(ff, aa, bb, xt, mf, similarity="ratio", + scale="derived3", guard="safe", + prestep="slope", + direction="dissimilar"), + f, a, b, root) + best = None + for impl in (chandrupatla, chandrupatla_scipy): + r = run(impl, f, a, b, root) + if r and (best is None or r[0] < best[0]): + best = r + for nm, r in (("ours", ours), ("chandrupatla", best)): + if r is None: + print(f"{pid:<9} {nm:<14} (not scorable)") + continue + n, red, share, halv = r + print(f"{pid:<9} {nm:<14} {n:6d} {red:19d} {share:6.1f}% {halv:14.3f}") + + +if __name__ == "__main__": + main() diff --git a/results_sign_flips.txt b/results_sign_flips.txt new file mode 100644 index 0000000..ecd2636 --- /dev/null +++ b/results_sign_flips.txt @@ -0,0 +1,11 @@ +problem method evals probes that shrank share halvings/eval +fun1.3 ours 19 15 78.9% 2.498 +fun1.3 chandrupatla 15 13 86.7% 3.099 +fun1.4 ours 59 54 91.5% 0.930 +fun1.4 chandrupatla 24 21 87.5% 2.256 +fun1.5 ours 63 60 95.2% 1.179 +fun1.5 chandrupatla 44 41 93.2% 1.684 +fun3.1 ours 64 61 95.3% 0.672 +fun3.1 chandrupatla 44 42 95.5% 0.955 +fun5.1 ours 57 47 82.5% 0.741 +fun5.1 chandrupatla 44 42 95.5% 0.955