1import json
  2import math
  3import numpy as np
  4
  5
  6def metric(x, R=1.0, delta=1e-10):
  7    x = np.asarray(x, dtype=float)
  8    q = 1.0 - float(x @ x) / (R * R)
  9    if q <= 0:
 10        raise ValueError("point is outside the open ball")
 11    # Hessian[-log(1-|x|^2/R^2)] = a I + b xx^T
 12    a = 2.0 / (R * R * q)
 13    b = 4.0 / (R**4 * q*q)
 14    return a * np.eye(len(x)) + b * np.outer(x, x) + delta * np.eye(len(x))
 15
 16
 17def barrier_value(x, R=1.0):
 18    q = 1.0 - float(x @ x) / (R * R)
 19    return -math.log(q) if q > 0 else float("inf")
 20
 21
 22def intrinsic_grad_sq(x, grad, R=1.0):
 23    G = metric(x, R)
 24    return float(grad @ np.linalg.solve(G, grad))
 25
 26
 27def project_ball(x, R=1.0, tol=1e-8):
 28    n = np.linalg.norm(x)
 29    lim = R * (1.0 - tol)
 30    return x if n < lim else x * (lim / max(n, 1e-30))
 31
 32
 33def verification():
 34    rng = np.random.default_rng(123)
 35    max_ratio = 0.0
 36    min_eig = float("inf")
 37    ratios = []
 38    # The theoretical bound is |grad g|_G^2 <= 1.
 39    for _ in range(5000):
 40        direction = rng.normal(size=2)
 41        direction /= np.linalg.norm(direction)
 42        # deliberately include points extremely near the boundary
 43        r = 10 ** rng.uniform(-4, -0.00005)
 44        x = r * direction
 45        q = 1 - x @ x
 46        grad_g = 2 * x / q
 47        val = intrinsic_grad_sq(x, grad_g)
 48        ratios.append(val)
 49        max_ratio = max(max_ratio, val)
 50        min_eig = min(min_eig, np.linalg.eigvalsh(metric(x)).min())
 51    # Check exact radial formula at several radii: 2u/(1+u).
 52    radial_errors = []
 53    for r in np.linspace(0.0, 0.999999, 100):
 54        x = np.array([r, 0.0])
 55        u = r*r
 56        observed = intrinsic_grad_sq(x, 2*x/(1-u))
 57        radial_errors.append(abs(observed - 2*u/(1+u)))
 58    return {
 59        "max_intrinsic_grad_g_squared": max_ratio,
 60        "max_bound_violation": max(0.0, max_ratio - 1.0),
 61        "minimum_metric_eigenvalue": min_eig,
 62        "max_radial_formula_error": max(radial_errors),
 63    }
 64
 65
 66def boundary_scaling_check():
 67    rows = []
 68    grad = np.array([1.0, 0.0])
 69    for r in [0.0, 0.5, 0.9, 0.99, 0.9999, 0.999999]:
 70        x = np.array([r, 0.0])
 71        euclidean_step = np.linalg.norm(grad)
 72        natural_step = np.linalg.norm(np.linalg.solve(metric(x), grad))
 73        rows.append({"radius": r, "euclidean_step": euclidean_step,
 74                     "barrier_step": natural_step,
 75                     "ratio": natural_step / euclidean_step})
 76    return rows
 77
 78
 79def run_optimizer(kind, target, steps=300, eta=0.18, R=1.0, seed=7):
 80    rng = np.random.default_rng(seed)
 81    x = np.array([0.0, 0.0])
 82    losses, radii, qvals, grad_norms = [], [], [], []
 83    # Smooth quadratic, plus small fixed noise to expose boundary overshoot behavior.
 84    for t in range(steps):
 85        grad = x - target
 86        noisy_grad = grad + 0.015 * rng.normal(size=2)
 87        if kind == "barrier-natural":
 88            G = metric(x, R)
 89            step = np.linalg.solve(G, noisy_grad)
 90            trial = x - eta * step
 91            # Backtracking is the specified projection safeguard.
 92            local_eta = eta
 93            while np.linalg.norm(trial) >= R * (1 - 1e-9):
 94                local_eta *= 0.5
 95                trial = x - local_eta * step
 96                if local_eta < 1e-12:
 97                    break
 98            x = trial
 99        elif kind == "projected-euclidean":
100            x = project_ball(x - eta * noisy_grad, R)
101        else:
102            raise ValueError(kind)
103        loss = 0.5 * float(np.sum((x - target) ** 2))
104        losses.append(loss)
105        radii.append(float(np.linalg.norm(x)))
106        qvals.append(float(1 - x @ x))
107        grad_norms.append(float(np.linalg.norm(noisy_grad)))
108    return {
109        "final_loss": losses[-1],
110        "best_loss": min(losses),
111        "final_radius": radii[-1],
112        "minimum_boundary_margin_q": min(qvals),
113        "boundary_hits_or_backtracks": int(sum(q <= 1e-8 for q in qvals)),
114        "max_gradient_norm": max(grad_norms),
115        "loss_first_20": float(np.mean(losses[:20])),
116        "loss_last_20": float(np.mean(losses[-20:])),
117    }
118
119
120def main():
121    check = verification()
122    # Target outside the ball: the constrained infimum is at the boundary.
123    target = np.array([1.20, 0.15])
124    results = {
125        "verification": check,
126        "boundary_step_scaling": boundary_scaling_check(),
127        "setup": {"domain": "||x|| < 1", "target": target.tolist(), "steps": 300, "eta": 0.18},
128        "projected_euclidean": run_optimizer("projected-euclidean", target),
129        "barrier_natural": run_optimizer("barrier-natural", target),
130        "repeat_summary": {
131            kind: {
132                "mean_final_loss": float(np.mean([run_optimizer(kind, target, seed=k)["final_loss"] for k in range(5)])),
133                "mean_final_radius": float(np.mean([run_optimizer(kind, target, seed=k)["final_radius"] for k in range(5)])),
134                "mean_min_q": float(np.mean([run_optimizer(kind, target, seed=k)["minimum_boundary_margin_q"] for k in range(5)])),
135            } for kind in ["projected-euclidean", "barrier-natural"]
136        },
137    }
138    print(json.dumps(results, indent=2))
139
140
141if __name__ == "__main__":
142    main()