vvsandClaude Opus 5
68df1d25db
Impact, every sweep regenerated: the paper's key configurations are unchanged to the evaluation (proposed 2806/3220, credit 3568/3020, Chandrupatla, TOMS 748, Brent, bisection); 46 stated numbers move (ITP CHA 2857 -> 2896, ITP oracle 2611 -> 2647, record CHA figures by 14-18, the grid means by ~1%) and one finding weakens: the measured-order exponent of section 7.1 was cheaper in 29 of 60 far cells and 7 of 7 safeguards; with the fix 12 of 61 and 6 of 7, and on the wide-bracket set it now loses in every safeguard and drops 5 problems. PREREG G1-G3 outcome on the fixed code: G2 (repeating step pattern, period from the data) PASSES Y1-Y4 and is ACCEPTED by the stop rule - proposed rule APS 2806 -> 2782 (below Chandrupatla's 2787), CHA 3220 -> 3169, both 6026 -> 5951, per problem vs Chandrupatla on APS 95/20/38 p = 8e-7. G1 fails Y1/Y2 (its Y3 miss was the sign bug), G3 fails Y2, A fails Y2. Regression 887 pass; the two red tests are the paper-number check (awaiting the paper update) and test_blend (fixed since, copy in stage16). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
391 lines
16 KiB
Python
391 lines
16 KiB
Python
"""Stage 4 - WHERE the similarity method wins, and WHY.
|
|
|
|
The totals say the parameter-free similarity rule loses to Brent and TOMS748
|
|
on Alefeld-Potra-Shi and beats them on Chandrupatla. A total does not explain
|
|
anything, so this file takes the runs apart.
|
|
|
|
Every root search has two phases with completely different economics:
|
|
|
|
APPROACH getting from the initial bracket into the neighbourhood of the
|
|
root. Cost here is governed by how boldly a method is willing
|
|
to jump, and by how often its safeguard drags it back to the
|
|
midpoint. Bracket width drives this phase.
|
|
|
|
ENDGAME going from the neighbourhood down to the tolerance. Cost here
|
|
is governed by the ORDER of convergence and nothing else:
|
|
secant 1.618, Brent's inverse quadratic 1.839, TOMS748's cubic
|
|
up to 2.7.
|
|
|
|
Prediction, if the method is what the derivation says it is:
|
|
- its step is |f|/|f'| with a secant slope, so its order is ~1.618 and it
|
|
must LOSE the endgame to both Brent and TOMS748, by a roughly constant
|
|
number of evaluations per problem;
|
|
- its approach carries no cautious truncation, so it must WIN the
|
|
approach, by an amount that GROWS with the width of the bracket.
|
|
- therefore: narrow brackets -> the endgame dominates -> we lose;
|
|
wide brackets -> the approach dominates -> we win.
|
|
|
|
All three are measured below. If the width correlation is absent, the
|
|
explanation is wrong and must not go into the paper.
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
from simroot.baselines import _opposite, _same_sign # 2026-09-21 fix
|
|
|
|
import math
|
|
|
|
from scipy.optimize._tstutils import (_APS_TESTS_DICTS,
|
|
_CHANDRUPATLA_TESTS_DICTS)
|
|
|
|
from simroot.baselines import bisection, brent, toms748
|
|
from simroot.pure import pure_root
|
|
|
|
XTOL = 1e-12
|
|
MAXFEV = 200
|
|
ACC = 1e-6
|
|
|
|
EPS = 2.220446049250313e-16
|
|
|
|
|
|
def phase_boundary(width, scale=0.0, xtol=XTOL):
|
|
"""Where the approach ends and the endgame begins.
|
|
|
|
A fixed relative distance does not work: at 1e-4 of a bracket of width
|
|
1e12 the boundary sits at 1e8, which the very first evaluation clears, so
|
|
every wide problem would report its entire cost as endgame.
|
|
|
|
The boundary is the geometric mean of the initial width and the tolerance
|
|
that is actually ATTAINABLE. The attainable one is not xtol: near
|
|
abscissae of magnitude 1e12 the spacing of floating point is about 2e-4,
|
|
so no method can drive the bracket to 1e-12 there and the journey is
|
|
shorter than xtol suggests. Using xtol regardless put the boundary too
|
|
low and charged the whole tail to the approach - measured on a width of
|
|
1e12, bisection split 40 / 15 instead of evenly.
|
|
|
|
`scale` is the largest |abscissa| in the bracket; pass it whenever the
|
|
bracket is far from the origin."""
|
|
attainable = max(xtol, EPS * abs(scale))
|
|
return math.sqrt(width * attainable)
|
|
|
|
|
|
def ours(f, a, b, xtol, maxfev):
|
|
return pure_root(f, a, b, xtol, maxfev, similarity="ratio",
|
|
scale="derived", guard="step")
|
|
|
|
|
|
METHODS = [("OURS (similarity)", ours), ("brent", brent),
|
|
("toms748", toms748), ("bisection", bisection)]
|
|
|
|
|
|
def traced(f):
|
|
"""Wrap f so the sequence of abscissae it is asked about is recorded."""
|
|
xs = []
|
|
|
|
def g(x):
|
|
xs.append(x)
|
|
return f(x)
|
|
|
|
return g, xs
|
|
|
|
|
|
def phases(method, f, a, b, root):
|
|
"""Split a run into approach and endgame evaluations."""
|
|
g, xs = traced(f)
|
|
try:
|
|
x, n, conv = method(g, a, b, XTOL, MAXFEV)
|
|
except Exception:
|
|
return None
|
|
if not conv or abs(x - root) > ACC:
|
|
return None
|
|
target = phase_boundary(abs(b - a), max(abs(a), abs(b)))
|
|
approach = len(xs)
|
|
for i, xi in enumerate(xs):
|
|
if abs(xi - root) <= target:
|
|
approach = i + 1
|
|
break
|
|
return approach, len(xs) - approach, len(xs)
|
|
|
|
|
|
def observed_order(method, f, a, b, root):
|
|
"""Order p from the tail of the running-best error sequence:
|
|
e_{n+1} ~ C e_n^p => p = log(e_{n+1}/e_n) / log(e_n/e_{n-1})."""
|
|
g, xs = traced(f)
|
|
try:
|
|
method(g, a, b, XTOL, MAXFEV)
|
|
except Exception:
|
|
return float("nan")
|
|
best, errs = float("inf"), []
|
|
for xi in xs:
|
|
e = abs(xi - root)
|
|
if e < best:
|
|
best = e
|
|
errs.append(e)
|
|
errs = [e for e in errs if e > 1e-14]
|
|
if len(errs) < 4:
|
|
return float("nan")
|
|
ps = []
|
|
for i in range(1, len(errs) - 1):
|
|
d1, d2 = errs[i] / errs[i - 1], errs[i + 1] / errs[i]
|
|
if 0 < d1 < 1 and 0 < d2 < 1:
|
|
ps.append(math.log(d2) / math.log(d1))
|
|
if not ps:
|
|
return float("nan")
|
|
return sorted(ps)[len(ps) // 2]
|
|
|
|
|
|
def replay_bracket(f, a, b, xs):
|
|
"""Reconstruct the bracket a bracketing method maintained, from the
|
|
sequence of abscissae it probed, and score every probe against it.
|
|
|
|
Returns (positions, shrinks, skipped) where
|
|
positions[i] where probe i fell inside the bracket, 0..1
|
|
(0.5 is the midpoint, i.e. a bisection step)
|
|
shrinks[i] the fraction of the bracket that probe i removed
|
|
skipped probes that fell outside the proven bracket and so
|
|
cannot be scored against it. Counting them matters:
|
|
dropping them silently would flatter any method that
|
|
probes outside, and scipy's toms748 does.
|
|
|
|
The reconstruction is exact for a bracketing method because the bracket
|
|
is determined by the signs of f, which are recomputed here.
|
|
"""
|
|
lo, hi = (a, b) if a < b else (b, a)
|
|
flo = f(lo)
|
|
positions, shrinks, skipped = [], [], 0
|
|
for xi in xs:
|
|
if not (lo < xi < hi):
|
|
skipped += 1
|
|
continue
|
|
w0 = hi - lo
|
|
positions.append((xi - lo) / w0)
|
|
fi = f(xi)
|
|
if _opposite(flo, fi):
|
|
hi = xi
|
|
else:
|
|
lo, flo = xi, fi
|
|
shrinks.append(1.0 - (hi - lo) / w0)
|
|
return positions, shrinks, skipped
|
|
|
|
|
|
def load():
|
|
out = []
|
|
for d in _APS_TESTS_DICTS:
|
|
a, b = d["bracket"]
|
|
args = d.get("args", ())
|
|
out.append(("APS", d["ID"], (lambda ff, aa: lambda x: ff(x, *aa))(d["f"], args),
|
|
float(a), float(b), float(d["root"])))
|
|
for d in _CHANDRUPATLA_TESTS_DICTS:
|
|
a, b = d["bracket"]
|
|
out.append(("CHA", d["ID"], (lambda ff: lambda x: ff(x))(d["f"]),
|
|
float(a), float(b), float(d["root"])))
|
|
return out
|
|
|
|
|
|
def main():
|
|
probs = load()
|
|
|
|
# ------------------------------------------------------------------
|
|
# 1. approach vs endgame, summed over the problems every method solves
|
|
# ------------------------------------------------------------------
|
|
print("=" * 78)
|
|
print("1. WHERE THE EVALUATIONS GO (only problems all four methods solve)")
|
|
print("=" * 78)
|
|
for setname in ("APS", "CHA"):
|
|
rows = {}
|
|
common = []
|
|
for s, pid, f, a, b, root in probs:
|
|
if s != setname:
|
|
continue
|
|
r = {nm: phases(m, f, a, b, root) for nm, m in METHODS}
|
|
if all(v is not None for v in r.values()):
|
|
common.append((pid, abs(b - a), r))
|
|
for nm, _m in METHODS:
|
|
app = sum(c[2][nm][0] for c in common)
|
|
end = sum(c[2][nm][1] for c in common)
|
|
rows[nm] = (app, end, app + end)
|
|
print(f"\n{setname}: {len(common)} problems solved by all four")
|
|
print(f"{'method':22s}{'approach':>11s}{'endgame':>10s}{'total':>9s}"
|
|
f"{'approach %':>12s}")
|
|
print("-" * 64)
|
|
for nm, _m in METHODS:
|
|
app, end, tot = rows[nm]
|
|
print(f"{nm:22s}{app:>11d}{end:>10d}{tot:>9d}"
|
|
f"{100.0 * app / tot:>11.1f}%")
|
|
|
|
# ------------------------------------------------------------------
|
|
# 2. observed order of convergence
|
|
# ------------------------------------------------------------------
|
|
print("\n" + "=" * 78)
|
|
print("2. OBSERVED ORDER OF CONVERGENCE (median over solved problems)")
|
|
print("=" * 78)
|
|
print(f"{'method':22s}{'median order':>14s}{'problems':>11s} theory")
|
|
print("-" * 66)
|
|
theory = {"OURS (similarity)": "1.618 (secant step |f|/|f'|)",
|
|
"brent": "1.839 (inverse quadratic)",
|
|
"toms748": "up to 2.7 (cubic)",
|
|
"bisection": "1.0 (linear, rate 1/2)"}
|
|
for nm, m in METHODS:
|
|
ps = [observed_order(m, f, a, b, root) for _s, _pid, f, a, b, root in probs]
|
|
ps = [p for p in ps if math.isfinite(p) and 0.5 < p < 5.0]
|
|
med = sorted(ps)[len(ps) // 2] if ps else float("nan")
|
|
print(f"{nm:22s}{med:>14.3f}{len(ps):>11d} {theory[nm]}")
|
|
|
|
# ------------------------------------------------------------------
|
|
# 3. the prediction: our margin over brent grows with bracket width
|
|
# ------------------------------------------------------------------
|
|
print("\n" + "=" * 78)
|
|
print("3. MARGIN OVER BRENT BY BRACKET WIDTH (negative = we spend fewer)")
|
|
print("=" * 78)
|
|
buckets = [(0, 1e1), (1e1, 1e3), (1e3, 1e6), (1e6, 1e15)]
|
|
labels = ["width < 10", "10 - 1e3", "1e3 - 1e6", "> 1e6"]
|
|
print(f"{'bracket width':16s}{'problems':>10s}{'OURS':>9s}{'brent':>8s}"
|
|
f"{'diff':>8s}{'diff/problem':>14s}")
|
|
print("-" * 66)
|
|
for (lo, hi), lab in zip(buckets, labels):
|
|
n_ours = n_brent = cnt = 0
|
|
for _s, _pid, f, a, b, root in probs:
|
|
w = abs(b - a)
|
|
if not (lo <= w < hi):
|
|
continue
|
|
ro = phases(ours, f, a, b, root)
|
|
rb = phases(brent, f, a, b, root)
|
|
if ro is None or rb is None:
|
|
continue
|
|
cnt += 1
|
|
n_ours += ro[2]
|
|
n_brent += rb[2]
|
|
if cnt:
|
|
d = n_ours - n_brent
|
|
print(f"{lab:16s}{cnt:>10d}{n_ours:>9d}{n_brent:>8d}{d:>8d}"
|
|
f"{d / cnt:>14.2f}")
|
|
print("\nThe trend is not monotone, so bracket width is a correlate and "
|
|
"not the cause.\nSection 4 tests the mechanism directly.")
|
|
|
|
# ------------------------------------------------------------------
|
|
# 4. the mechanism: how often does each method's safeguard fire?
|
|
# An interpolation step that gets thrown away and replaced by a
|
|
# midpoint is a wasted evaluation. If our advantage is real, our
|
|
# safeguard must fire LESS often than Brent's exactly on the set
|
|
# where we win, and not much less where we lose.
|
|
# ------------------------------------------------------------------
|
|
print("\n" + "=" * 78)
|
|
print("4. HOW OFTEN THE SAFEGUARD THROWS AN INTERPOLATION AWAY")
|
|
print("=" * 78)
|
|
print(f"{'set':6s}{'method':22s}{'steps':>9s}{'forced to midpoint':>21s}"
|
|
f"{'share':>9s}")
|
|
print("-" * 68)
|
|
for setname in ("APS", "CHA"):
|
|
for nm in ("OURS (similarity)", "brent"):
|
|
steps = forced = 0
|
|
for s, _pid, f, a, b, root in probs:
|
|
if s != setname:
|
|
continue
|
|
st = {}
|
|
if nm == "brent":
|
|
x, n, conv = brent(f, a, b, XTOL, MAXFEV, stats=st)
|
|
fired = st.get("bisect", 0)
|
|
else:
|
|
x, n, conv = pure_root(f, a, b, XTOL, MAXFEV,
|
|
similarity="ratio", scale="derived",
|
|
guard="step", stats=st)
|
|
fired = st.get("forced", 0)
|
|
if not conv or abs(x - root) > ACC:
|
|
continue
|
|
steps += n - 2 # the two endpoint evaluations are free
|
|
forced += fired
|
|
share = 100.0 * forced / steps if steps else float("nan")
|
|
print(f"{setname:6s}{nm:22s}{steps:>9d}{forced:>21d}"
|
|
f"{share:>8.1f}%")
|
|
print("\nREFUTED: Brent throws away 53% of its interpolations on APS, "
|
|
"where it WINS,\nand only 34% on Chandrupatla, where it LOSES. "
|
|
"The safeguard rate is not the\nmechanism. Section 5 tests the "
|
|
"next candidate.")
|
|
|
|
# ------------------------------------------------------------------
|
|
# 5. how BOLD are the steps?
|
|
#
|
|
# Brent accepts an interpolation only if 2|p| < 3*xm*q, i.e. only if
|
|
# it lands within three quarters of the way to the far end; otherwise
|
|
# the midpoint is used. The similarity rule has no such cap - its step
|
|
# is |f|/|f'| and may cross almost the whole bracket in one go. On a
|
|
# wide bracket whose root sits far from the middle, that cap is the
|
|
# difference between arriving in one step and bisecting towards it.
|
|
#
|
|
# Measured as the position t of each probe inside the bracket it was
|
|
# proposed in, t = (x - lo)/(hi - lo). t near 0 or 1 is a bold step;
|
|
# t = 0.5 is bisection. The bracket is reconstructed from the sign of
|
|
# f at each probe, which is exact for a bracketing method.
|
|
# ------------------------------------------------------------------
|
|
print("\n" + "=" * 78)
|
|
print("5. HOW BOLD IS EACH STEP (position inside the bracket it was "
|
|
"proposed in)")
|
|
print("=" * 78)
|
|
print(f"{'set':6s}{'method':22s}{'probes':>9s}{'bold (t<.1 or t>.9)':>22s}"
|
|
f"{'share':>9s}")
|
|
print("-" * 70)
|
|
for setname in ("APS", "CHA"):
|
|
for nm, m in (("OURS (similarity)", ours), ("brent", brent)):
|
|
probes = bold = 0
|
|
for s, _pid, f, a, b, root in probs:
|
|
if s != setname:
|
|
continue
|
|
g, xs = traced(f)
|
|
try:
|
|
x, _n, conv = m(g, a, b, XTOL, MAXFEV)
|
|
except Exception:
|
|
continue
|
|
if not conv or abs(x - root) > ACC:
|
|
continue
|
|
positions, _shrinks, _skipped = replay_bracket(f, a, b, xs[2:])
|
|
probes += len(positions)
|
|
bold += sum(1 for t in positions if t < 0.1 or t > 0.9)
|
|
share = 100.0 * bold / probes if probes else float("nan")
|
|
print(f"{setname:6s}{nm:22s}{probes:>9d}{bold:>22d}{share:>8.1f}%")
|
|
print("\nREFUTED: on CHA Brent is BOLDER than we are (56.3% vs 41.8%) and "
|
|
"still loses.\nBoldness is not the mechanism either. Section 6 "
|
|
"measures step QUALITY rather\nthan step size.")
|
|
|
|
# ------------------------------------------------------------------
|
|
# 6. step quality: how much of the bracket does one probe remove?
|
|
#
|
|
# Boldness is how far a probe goes; quality is whether it went the
|
|
# right way. A probe at t = 0.95 that lands on the wrong side of the
|
|
# root removes 5% of the bracket; the same probe on the right side
|
|
# removes 95%. This is the one quantity that must differ if a method
|
|
# with a LOWER order of convergence nevertheless needs fewer probes.
|
|
# ------------------------------------------------------------------
|
|
print("\n" + "=" * 78)
|
|
print("6. STEP QUALITY (fraction of the bracket removed per probe)")
|
|
print("=" * 78)
|
|
print(f"{'set':6s}{'method':22s}{'probes':>9s}{'median removed':>17s}"
|
|
f"{'mean removed':>15s}{'probes skipped':>16s}")
|
|
print("-" * 86)
|
|
for setname in ("APS", "CHA"):
|
|
for nm, m in (("OURS (similarity)", ours), ("brent", brent),
|
|
("toms748", toms748), ("bisection", bisection)):
|
|
shrinks, skipped = [], 0
|
|
for s, _pid, f, a, b, root in probs:
|
|
if s != setname:
|
|
continue
|
|
g, xs = traced(f)
|
|
try:
|
|
x, _n, conv = m(g, a, b, XTOL, MAXFEV)
|
|
except Exception:
|
|
continue
|
|
if not conv or abs(x - root) > ACC:
|
|
continue
|
|
_pos, sh, sk = replay_bracket(f, a, b, xs[2:])
|
|
shrinks.extend(sh)
|
|
skipped += sk
|
|
if shrinks:
|
|
srt = sorted(shrinks)
|
|
med = srt[len(srt) // 2]
|
|
print(f"{setname:6s}{nm:22s}{len(shrinks):>9d}{med:>17.3f}"
|
|
f"{sum(shrinks) / len(shrinks):>15.3f}{skipped:>16d}")
|
|
print("\nBisection removes exactly 0.500 by construction and is the "
|
|
"reference row.")
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|