SabaPivot's picture
download
raw
19.1 kB
"""Claim 3 -- Theorem 3.5 / Propositions D.5-D.6 of arXiv:2601.20180.
"Theorem 3.5 shows a phase transition to tractability: when rho <= 1 + O_eps(eps^4),
an eps-performatively stable point can be computed in poly(d, log(1/eps)) time via
an ellipsoid-method algorithm exploiting hypomonotonicity."
Executable programme:
A. the exact algebraic identity (11) behind Theorem D.1, and the exact
hypomonotonicity modulus of F = id - T for a (1+sigma)-expansive T.
B. symbolic check that sigma <= eps in Prop D.6 yields eps' = Theta(eps^{1/4}),
and, inverted, that a target accuracy delta needs rho <= 1 + Theta(delta^4).
C. a numerically stable (square-root form) ellipsoid method for the VI, run on
monotone instances; iteration counts fitted against C * d^a * log(1/eps)^b.
D. the phase transition: certificate scaling in sigma, the mechanism by which
hypomonotonicity breaks the ellipsoid cut, a nonlinear expansive family, and
a head-to-head against repeated risk minimisation at rho >= 1.
"""
import numpy as np
import sympy as sp
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from common import FIGS, dump
SEED = 20260725
rng = np.random.default_rng(SEED)
res = {"seed": SEED, "claim": "Theorem 3.5 phase transition"}
# =====================================================================
# A. hypomonotonicity identity and its exact modulus
# =====================================================================
ux, uy, wx, wy = sp.symbols("ux uy wx wy", real=True)
U, W = sp.Matrix([ux, uy]), sp.Matrix([wx, wy])
lhs = (U - W).dot(U) # <F(x)-F(x'), x-x'> with F = id - T, w = T(x)-T(x')
rhs = sp.Rational(1, 2) * (U.dot(U) - W.dot(W) + (U - W).dot(U - W))
res["identity_11"] = {
"statement": "<F(x)-F(x'),x-x'> = 1/2(||u||^2 - ||T(x)-T(x')||^2 + ||u-(T(x)-T(x'))||^2)",
"symbolically_exact": bool(sp.simplify(lhs - rhs) == 0),
"simplifies_to": "<u, u - w>",
}
print(
"A. identity (11) exact:", res["identity_11"]["symbolically_exact"], "-> <u, u-w>"
)
worst = []
for sigma in [0.0, 1e-4, 1e-2, 0.1, 0.5, 1.0]:
best = np.inf
for _ in range(100000):
uu = rng.normal(size=3)
uu /= np.linalg.norm(uu)
ww = rng.normal(size=3)
ww = ww / np.linalg.norm(ww) * (1 + sigma) * rng.uniform(0, 1)
best = min(best, float(np.dot(uu, uu - ww)))
# the analytic worst case w = (1+sigma)u, which the random search only
# approaches; including it shows the minimum is attained and equals -sigma
# (fixed vector, so the shared rng stream is untouched)
uu = np.array([1.0, 0.0, 0.0])
best = min(best, float(np.dot(uu, uu - (1 + sigma) * uu)))
worst.append(
{
"sigma": sigma,
"empirical_worst_modulus": -best,
"exact_worst_modulus": sigma,
"paper_bound_sigma_plus_sigma_sq_over_2": sigma + sigma**2 / 2,
"paper_bound_valid": bool(-best <= sigma + sigma**2 / 2 + 1e-9),
"empirical_reaches_exact_within": abs(-best - sigma),
}
)
res["hypomonotonicity_modulus"] = {
"note": "1/2(||u||^2-||w||^2+||u-w||^2) = <u,u-w>, whose minimum over "
"||w|| <= (1+sigma)||u|| is exactly -sigma||u||^2, so F = id - T is "
"exactly sigma-hypomonotone. The paper's sigma+sigma^2/2 is valid but "
"loose by sigma^2/2.",
"sweep": worst,
}
print(
"A. paper's bound sigma+sigma^2/2 valid in %d/%d; exact modulus is sigma "
"(max empirical deviation %.3e)"
% (
sum(w["paper_bound_valid"] for w in worst),
len(worst),
max(w["empirical_reaches_exact_within"] for w in worst),
)
)
# =====================================================================
# B. the eps^{1/4} <-> eps^4 exponent, symbolically
# =====================================================================
e, D, s, delta = sp.symbols("epsilon D sigma delta", positive=True)
eps_prime_expr = sp.sqrt(2 * D * sp.sqrt((2 + s) * (e + (s + s**2 / 2) * D**2)))
lim = sp.simplify(sp.limit(eps_prime_expr.subs(s, e) / e ** sp.Rational(1, 4), e, 0))
sol = sp.solve(sp.Eq(sp.sqrt(2 * D * sp.sqrt(2 * (e + e * D**2))), delta), e)
res["exponent"] = {
"prop_D6_formula": str(eps_prime_expr),
"limit_eps_prime_over_eps_to_the_quarter": str(lim),
"exponent_is_one_quarter": bool(lim.is_finite and lim != 0),
"inverted_eps_as_function_of_target_delta": [str(x) for x in sol],
"inverted_exponent_is_4": bool(
any(sp.degree(sp.simplify(sp.expand(x)), delta) == 4 for x in sol)
),
}
print(
"B. eps'/eps^{1/4} ->",
str(lim),
"| inversion gives eps ~ delta^4:",
res["exponent"]["inverted_exponent_is_4"],
)
# =====================================================================
# C. a numerically stable ellipsoid method for the VI
# =====================================================================
def ball_gap(x, Fx, R):
"""max_{||y||<=R} <Fx, x-y> = <Fx,x> + R||Fx|| -- the VI gap on the ball."""
return float(np.dot(Fx, x) + R * np.linalg.norm(Fx))
def ellipsoid_vi(F, d, R=1.0, iters=20000, target=None, x_true=None):
"""Square-root-form ellipsoid method on X = {||x||_2 <= R}.
E_k = {c + B u : ||u|| <= 1}. Cut: <F(c), x - c> <= 0, valid for every VI
solution when F is monotone (Theorem D.1 / Prop D.5). Also records how often
a known solution x_true violates the cut -- the mechanism hypomonotonicity
attacks.
"""
c = np.zeros(d)
B = R * np.sqrt(d + 1.0) * np.eye(d)
tau = 1.0 - np.sqrt((d - 1.0) / (d + 1.0))
scale = d / np.sqrt(d * d - 1.0)
best_gap, best_x = np.inf, c.copy()
cuts, bad_cuts = 0, 0
for k in range(iters):
nc = np.linalg.norm(c)
if nc > R: # feasibility cut for the ball
a = c / nc
else:
Fx = F(c)
gap = ball_gap(c, Fx, R)
if gap < best_gap:
best_gap, best_x = gap, c.copy()
if target is not None and best_gap <= target:
return best_x, best_gap, k + 1, cuts, bad_cuts
na = np.linalg.norm(Fx)
if na < 1e-300:
return best_x, best_gap, k + 1, cuts, bad_cuts
a = Fx / na
cuts += 1
if x_true is not None and float(np.dot(a, x_true - c)) > 1e-15:
bad_cuts += 1 # the true solution is cut away -> invalid cut
Bta = B.T @ a
n = np.linalg.norm(Bta)
if not np.isfinite(n) or n < 1e-300:
break
bt = Bta / n
c = c - (B @ bt) / (d + 1.0)
B = scale * (B - tau * np.outer(B @ bt, bt))
if not np.isfinite(B).all():
break
return best_x, best_gap, iters, cuts, bad_cuts
def monotone_instance(d, rng, mono_shift=0.0):
"""F(x) = A(x - x0) with A skew(+mono_shift*I): monotone, unique solution x0."""
Bm = rng.normal(size=(d, d))
A = (Bm - Bm.T) / np.sqrt(2.0 * d) + mono_shift * np.eye(d)
x0 = rng.uniform(-0.3, 0.3, size=d)
return (lambda x, A=A, x0=x0: A @ (x - x0)), A, x0
scaling = []
for d in [2, 4, 8, 16, 32]:
for eps in [1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8]:
for trial in range(3):
F, A, x0 = monotone_instance(d, rng)
cap = int(8 * d * (d + 1) * np.log(1 / eps)) + 2000
_, gap, it, cuts, bad = ellipsoid_vi(
F, d, R=1.0, iters=cap, target=eps, x_true=x0
)
scaling.append(
{
"d": d,
"eps": eps,
"trial": trial,
"iterations": it,
"final_gap": gap,
"reached": bool(gap <= eps),
"invalid_cuts": bad,
}
)
n_reached = sum(r["reached"] for r in scaling)
X = np.array(
[
[1.0, np.log(r["d"]), np.log(np.log(1.0 / r["eps"]))]
for r in scaling
if r["reached"]
]
)
y = np.array([np.log(r["iterations"]) for r in scaling if r["reached"]])
coef, *_ = np.linalg.lstsq(X, y, rcond=None)
r2 = 1 - ((y - X @ coef) ** 2).sum() / ((y - y.mean()) ** 2).sum()
res["ellipsoid_scaling"] = {
"model": "iterations = C * d^a * (log(1/eps))^b",
"C": float(np.exp(coef[0])),
"exponent_a_in_d": float(coef[1]),
"exponent_b_in_log_1_over_eps": float(coef[2]),
"R2": float(r2),
"predicted_a": 2.0,
"predicted_b": 1.0,
"n_runs": len(scaling),
"n_reached_target": int(n_reached),
"total_invalid_cuts_on_monotone_instances": int(
sum(r["invalid_cuts"] for r in scaling)
),
"measurements": scaling,
}
print(
"C. %d/%d runs hit target; fit iterations ~ %.3f * d^%.3f * (log 1/eps)^%.3f "
"(R^2 = %.4f; theory a=2, b=1); invalid cuts on monotone instances: %d"
% (
n_reached,
len(scaling),
np.exp(coef[0]),
coef[1],
coef[2],
r2,
sum(r["invalid_cuts"] for r in scaling),
)
)
# =====================================================================
# D. the phase transition
# =====================================================================
def rotation_expansive(d, sigma, theta, rng):
"""g(x) = M(x-x0)+x0 with M = (1+sigma)*blockdiag(R(theta)).
Lipschitz constant of g is exactly 1+sigma, so rho = L*beta/alpha = 1+sigma;
the hypomonotonicity modulus of F = id - g is max(0, (1+sigma)cos(theta)-1).
"""
Q = np.zeros((d, d))
for i in range(d // 2):
ct, st = np.cos(theta), np.sin(theta)
Q[2 * i : 2 * i + 2, 2 * i : 2 * i + 2] = np.array([[ct, -st], [st, ct]])
M = (1 + sigma) * Q
x0 = rng.uniform(-0.3, 0.3, size=d)
return (lambda x, M=M, x0=x0: M @ (x - x0) + x0), M, x0
def rrm_best_gap(g, d, R=1.0, iters=5000, rng=None):
x = rng.uniform(-0.3, 0.3, size=d)
best = np.inf
for _ in range(iters):
best = min(best, float(np.linalg.norm(x - g(x))))
x = g(x)
n = np.linalg.norm(x)
if n > R:
x = x / n * R
return best
d, Ddiam, eps_evi = 4, 2.0, 1e-12
phase = []
for theta_deg in [0.0, 30.0, 60.0, 90.0]:
for sigma in [0.0, 1e-6, 1e-4, 1e-3, 1e-2, 0.05, 0.1, 0.3, 0.5, 1.0]:
g, M, x0 = rotation_expansive(d, sigma, np.deg2rad(theta_deg), rng)
F = lambda x, g=g: x - g(x)
sym = (np.eye(d) - M + (np.eye(d) - M).T) / 2
sigma_hypo = float(max(0.0, -np.linalg.eigvalsh(sym).min()))
bx, gap, it, cuts, bad = ellipsoid_vi(F, d, R=1.0, iters=6000, x_true=x0)
fp_gap = float(np.linalg.norm(bx - g(bx)))
# Prop D.5 / D.6 certificates, using the generic Lipschitz bound L <= 2+sigma
vi_cert = 2 * Ddiam * np.sqrt((2 + sigma) * (eps_evi + sigma_hypo * Ddiam**2))
fp_cert = np.sqrt(
2
* Ddiam
* np.sqrt((2 + sigma) * (eps_evi + (sigma + sigma**2 / 2) * Ddiam**2))
)
phase.append(
{
"theta_deg": theta_deg,
"sigma": sigma,
"rho": 1 + sigma,
"sigma_hypomonotone": sigma_hypo,
"ellipsoid_vi_gap": gap,
"ellipsoid_fixed_point_gap": fp_gap,
"dist_to_true_fixed_point": float(np.linalg.norm(bx - x0)),
"rrm_best_fixed_point_gap": rrm_best_gap(g, d, rng=rng),
"prop_D5_vi_certificate": float(vi_cert),
"prop_D6_fixed_point_certificate": float(fp_cert),
"achieved_within_vi_certificate": bool(gap <= vi_cert + 1e-12),
"achieved_within_fp_certificate": bool(fp_gap <= fp_cert + 1e-12),
"cuts_applied": cuts,
"invalid_cuts": bad,
"invalid_cut_fraction": float(bad / max(cuts, 1)),
"iterations": it,
}
)
ok_vi = sum(p["achieved_within_vi_certificate"] for p in phase)
ok_fp = sum(p["achieved_within_fp_certificate"] for p in phase)
res["phase_transition"] = phase
print(
"D. Prop D.5 VI certificate holds in %d/%d; Prop D.6 fixed-point certificate "
"holds in %d/%d" % (ok_vi, len(phase), ok_fp, len(phase))
)
sel = [p for p in phase if p["theta_deg"] == 0.0 and p["sigma"] > 0]
sig = np.array([p["sigma"] for p in sel])
slope_vi = float(
np.polyfit(np.log(sig), np.log([p["prop_D5_vi_certificate"] for p in sel]), 1)[0]
)
slope_fp = float(
np.polyfit(
np.log(sig), np.log([p["prop_D6_fixed_point_certificate"] for p in sel]), 1
)[0]
)
res["certificate_slopes"] = {
"vi_certificate_slope_vs_sigma": slope_vi,
"predicted_vi_slope": 0.5,
"fixed_point_certificate_slope_vs_sigma": slope_fp,
"predicted_fixed_point_slope": 0.25,
"interpretation": "a target fixed-point accuracy delta therefore needs "
"sigma = rho - 1 ~ delta^4, which is Theorem 3.5.",
}
print(
"D. certificate slopes vs sigma: VI %.4f (theory 0.5), fixed point %.4f "
"(theory 0.25)" % (slope_vi, slope_fp)
)
# --- mechanism: how often does hypomonotonicity invalidate the ellipsoid cut?
mech = {}
for theta_deg in [0.0, 30.0, 60.0, 90.0]:
rows = [p for p in phase if p["theta_deg"] == theta_deg]
mech["theta_%d" % int(theta_deg)] = [
{
"sigma": p["sigma"],
"sigma_hypo": p["sigma_hypomonotone"],
"invalid_cut_fraction": p["invalid_cut_fraction"],
}
for p in rows
]
res["cut_validity_mechanism"] = mech
print(
"D. invalid-cut fraction at theta=0: "
+ ", ".join(
"sigma=%g:%.2f" % (r["sigma"], r["invalid_cut_fraction"])
for r in mech["theta_0"]
)
)
# --- a nonlinear expansive family, where the ellipsoid can actually degrade
def nonlinear_expansive(d, sigma, amp, rng, theta=np.pi / 2):
Q = np.zeros((d, d))
for i in range(d // 2):
ct, st = np.cos(theta), np.sin(theta)
Q[2 * i : 2 * i + 2, 2 * i : 2 * i + 2] = np.array([[ct, -st], [st, ct]])
M = (1 + sigma) * Q
om = 9.0
x0 = rng.uniform(-0.2, 0.2, size=d)
def g(x):
return M @ (x - x0) + x0 + amp * np.sin(om * x)
return g, M, x0
nl = []
for sigma in [0.0, 1e-3, 1e-2, 0.05, 0.1, 0.3, 0.5]:
for amp in [0.0, 0.02, 0.05]:
g, M, x0 = nonlinear_expansive(4, sigma, amp, rng)
F = lambda x, g=g: x - g(x)
# empirical Lipschitz constant of g (i.e. rho) over random secants
Lemp = 0.0
for _ in range(4000):
a1 = rng.uniform(-1, 1, size=4)
a2 = a1 + rng.normal(size=4) * 0.05
Lemp = max(
Lemp, float(np.linalg.norm(g(a1) - g(a2)) / np.linalg.norm(a1 - a2))
)
bx, gap, it, cuts, bad = ellipsoid_vi(F, 4, R=1.0, iters=8000)
nl.append(
{
"sigma": sigma,
"sine_amplitude": amp,
"empirical_rho": Lemp,
"ellipsoid_vi_gap": gap,
"ellipsoid_fixed_point_gap": float(np.linalg.norm(bx - g(bx))),
"rrm_best_fixed_point_gap": rrm_best_gap(g, 4, rng=rng),
"invalid_cut_fraction": float(bad / max(cuts, 1)),
"iterations": it,
}
)
res["nonlinear_family"] = nl
print("D. nonlinear expansive family (rho measured empirically):")
for r in nl:
print(
" sigma=%-5g amp=%-4g rho_emp=%.3f | ellipsoid fp gap %.3e | "
"RRM fp gap %.3e"
% (
r["sigma"],
r["sine_amplitude"],
r["empirical_rho"],
r["ellipsoid_fixed_point_gap"],
r["rrm_best_fixed_point_gap"],
)
)
# --- RRM vs ellipsoid head to head
comparison = []
for theta_deg, sigma in [
(90.0, 0.0),
(90.0, 0.01),
(90.0, 0.1),
(60.0, 0.0),
(30.0, 0.01),
(0.0, 0.01),
]:
g, M, x0 = rotation_expansive(4, sigma, np.deg2rad(theta_deg), rng)
F = lambda x, g=g: x - g(x)
bx, gap, it, cuts, bad = ellipsoid_vi(F, 4, R=1.0, iters=6000)
comparison.append(
{
"theta_deg": theta_deg,
"rho": 1 + sigma,
"ellipsoid_vi_gap": gap,
"ellipsoid_iterations": it,
"ellipsoid_fixed_point_gap": float(np.linalg.norm(bx - g(bx))),
"rrm_best_fixed_point_gap_5000_iters": rrm_best_gap(g, 4, rng=rng),
}
)
res["rrm_vs_ellipsoid"] = comparison
for c_ in comparison:
print(
"D. theta=%3.0f deg rho=%.2f : ellipsoid fp gap %.2e (%d iters) | RRM fp gap %.2e"
% (
c_["theta_deg"],
c_["rho"],
c_["ellipsoid_fixed_point_gap"],
c_["ellipsoid_iterations"],
c_["rrm_best_fixed_point_gap_5000_iters"],
)
)
# ------------------------------------------------------------------ figures
fig, axes = plt.subplots(1, 3, figsize=(15, 4.2))
ax = axes[0]
for dd in [2, 4, 8, 16, 32]:
pts = sorted(
[(r["eps"], r["iterations"]) for r in scaling if r["d"] == dd and r["reached"]]
)
if pts:
xs = sorted(set(p[0] for p in pts))
ys = [np.median([p[1] for p in pts if p[0] == x]) for x in xs]
ax.plot([np.log(1 / x) for x in xs], ys, "o-", label="d=%d" % dd)
ax.set_xlabel("log(1/eps)")
ax.set_ylabel("ellipsoid iterations")
ax.set_title("iterations grow linearly in log(1/eps)")
ax.legend(fontsize=8)
ax.grid(alpha=0.3)
ax = axes[1]
ds = [2, 4, 8, 16, 32]
for eps in [1e-3, 1e-5, 1e-8]:
ys = [
np.median(
[
r["iterations"]
for r in scaling
if r["d"] == dd and r["eps"] == eps and r["reached"]
]
)
for dd in ds
]
ax.loglog(ds, ys, "o-", label="eps=%g" % eps)
ax.loglog(
ds, [ys[0] * (np.array(ds) / ds[0]) ** 2][0], "k--", alpha=0.5, label="slope 2"
)
ax.set_xlabel("dimension d")
ax.set_ylabel("ellipsoid iterations")
ax.set_title("iterations ~ d^%.2f (theory 2)" % coef[1])
ax.legend(fontsize=8)
ax.grid(alpha=0.3, which="both")
ax = axes[2]
ax.loglog(
sig,
[p["prop_D6_fixed_point_certificate"] for p in sel],
"o-",
label="Prop D.6 certificate (slope %.2f)" % slope_fp,
)
ax.loglog(
sig,
[p["prop_D5_vi_certificate"] for p in sel],
"s-",
label="Prop D.5 VI certificate (slope %.2f)" % slope_vi,
)
ax.loglog(
sig,
[max(p["ellipsoid_fixed_point_gap"], 1e-17) for p in sel],
"^--",
label="achieved fixed-point gap",
)
ax.set_xlabel("sigma (rho = 1 + sigma)")
ax.set_ylabel("accuracy")
ax.set_title("guaranteed accuracy vs expansiveness")
ax.legend(fontsize=8)
ax.grid(alpha=0.3, which="both")
plt.tight_layout()
p = FIGS + "/claim3_ellipsoid_phase.png"
plt.savefig(p, dpi=130)
print("wrote", p)
res["verdict"] = {
"identity_exact": res["identity_11"]["symbolically_exact"],
"eps_quarter_exponent_confirmed": res["exponent"]["exponent_is_one_quarter"],
"poly_d_log_eps_confirmed": bool(coef[1] < 3.0 and coef[2] < 2.0 and r2 > 0.9),
"certificates_hold": bool(ok_vi == len(phase) and ok_fp == len(phase)),
"paper_constant_loose_by": "sigma^2/2 (the exact modulus is sigma)",
}
dump("claim3_ellipsoid_phase.json", res)

Xet Storage Details

Size:
19.1 kB
·
Xet hash:
9117971f2796df05506967c983016c9cd479e2606ba84ac966210bf33b1ac01c

Xet efficiently stores files, intelligently splitting them into unique chunks and accelerating uploads and downloads. More info.