1import json
  2import math
  3import numpy as np
  4
  5
  6def euclidean_fix(s, y):
  7    c = float(s @ y)
  8    d = abs(c) - c
  9    return y + d * s / float(s @ s)
 10
 11
 12def metric_fix(s, y, B):
 13    c = float(s @ y)
 14    d = abs(c) - c
 15    return y + d * (B @ s) / float(s @ B @ s)
 16
 17
 18def bfgs_update(B, s, y):
 19    Bs = B @ s
 20    return B - np.outer(Bs, Bs) / (s @ Bs) + np.outer(y, y) / (s @ y)
 21
 22
 23# A two-loop L-BFGS inverse-Hessian application. Pairs are (s,y), with y.T s > 0.
 24def two_loop(g, pairs, H0_scale=1.0):
 25    q = g.copy()
 26    alphas = []
 27    for s, y in reversed(pairs):
 28        rho = 1.0 / (s @ y)
 29        a = rho * (s @ q)
 30        alphas.append(a)
 31        q = q - a * y
 32    if pairs:
 33        s, y = pairs[-1]
 34        gamma = (s @ y) / max(y @ y, 1e-30)
 35    else:
 36        gamma = H0_scale
 37    r = gamma * q
 38    for (s, y), a in zip(pairs, reversed(alphas)):
 39        rho = 1.0 / (s @ y)
 40        b = rho * (y @ r)
 41        r = r + s * (a - b)
 42    return r
 43
 44
 45def objective(x):
 46    # Smooth bounded nonconvex quartic in coordinate zero plus a stiff quadratic block.
 47    return -0.5*x[0]**2 + 0.10*x[0]**4 + 0.05*np.sum(x[1:]**2)
 48
 49
 50def grad(x):
 51    g = 0.1*x.copy()
 52    g[0] = -x[0] + 0.4*x[0]**3
 53    return g
 54
 55
 56def run_lbfgs(correct=False, seed=7, n=12, steps=180, memory=7):
 57    rng = np.random.default_rng(seed)
 58    x = np.r_[0.35, 0.1*rng.normal(0, 1.0, n-1)]
 59    pairs = []
 60    f0 = objective(x)
 61    corrected = 0
 62    negative = 0
 63    backtracks = 0
 64    accepted = []
 65    for _ in range(steps):
 66        g = grad(x)
 67        p = -two_loop(g, pairs, H0_scale=1.0)
 68        if not np.isfinite(p).all() or g @ p >= 0:
 69            p = -g
 70        f = objective(x)
 71        slope = g @ p
 72        t = 1.0
 73        for _bt in range(35):
 74            xn = x + t*p
 75            if objective(xn) <= f + 1e-4*t*slope:
 76                break
 77            t *= 0.5
 78        backtracks += int(round(-math.log2(t)))
 79        if t < 2**-34:
 80            break
 81        gn = grad(xn)
 82        s = xn - x
 83        y = gn - g
 84        c = float(s @ y)
 85        scale = np.linalg.norm(s)*np.linalg.norm(y)
 86        if c < -1e-4*scale:
 87            negative += 1
 88            if correct:
 89                y = euclidean_fix(s, y)
 90                corrected += 1
 91        if s @ y > 1e-12*max(np.linalg.norm(s)*np.linalg.norm(y), 1.0):
 92            pairs.append((s, y))
 93            pairs = pairs[-memory:]
 94        x = xn
 95        accepted.append(objective(x))
 96    return {
 97        'initial_loss': f0, 'final_loss': float(objective(x)),
 98        'best_loss': float(np.min(accepted)) if accepted else f0,
 99        'negative_pairs': negative, 'corrected_pairs': corrected,
100        'backtracking_halvings': backtracks, 'steps': len(accepted),
101        'finite': bool(np.isfinite(x).all() and np.isfinite(objective(x)))
102    }
103
104
105def main():
106    rng = np.random.default_rng(123)
107    # Prediction 1/2: for c=-rho||s||||y||, Euclidean correction has
108    # s.T ytilde=|c| and ||delta||/||y||=2rho exactly.
109    sweep = []
110    for rho in [0.001, 0.01, 0.05, 0.1, 0.25, 0.5, 0.9]:
111        s = rng.normal(size=9)
112        y0 = rng.normal(size=9)
113        # Construct y with prescribed normalized negative curvature.
114        sh = s / np.linalg.norm(s)
115        perp = y0 - sh*(sh @ y0)
116        perp /= np.linalg.norm(perp)
117        y = -rho*np.linalg.norm(y0)*sh + math.sqrt(max(1-rho*rho, 0))*np.linalg.norm(y0)*perp
118        c = s @ y
119        yt = euclidean_fix(s, y)
120        delta_ratio = np.linalg.norm(yt-y)/np.linalg.norm(y)
121        sweep.append({'rho': rho, 'observed_secant_over_abs_c': float((s@yt)/abs(c)),
122                      'predicted_secant_over_abs_c': 1.0,
123                      'observed_delta_ratio': float(delta_ratio),
124                      'predicted_delta_ratio': 2*rho})
125
126    # Prediction 3: the BFGS update is SPD for every corrected pair.
127    pd_sweep = []
128    for dim in [3, 8, 20]:
129        mins = []
130        for trial in range(40):
131            A = rng.normal(size=(dim, dim))
132            B = A.T@A + 0.2*np.eye(dim)
133            s = rng.normal(size=dim)
134            y = -abs(rng.normal())*s + rng.normal(size=dim)*0.2
135            if s@y >= 0: y = y - (abs(s@y)+0.3)*s/(s@s)
136            yt = euclidean_fix(s, y)
137            Bn = bfgs_update(B, s, yt)
138            mins.append(float(np.linalg.eigvalsh((Bn+Bn.T)/2).min()))
139        pd_sweep.append({'dimension': dim, 'minimum_observed_eigenvalue': min(mins),
140                         'predicted_lower_bound': 0.0})
141
142    baseline = run_lbfgs(correct=False)
143    idea = run_lbfgs(correct=True)
144    result = {
145        'math_verification': {
146            'prediction_1_secant_exactness': 's.T ytilde / abs(s.T y) = 1 for every rho',
147            'prediction_2_correction_scaling': '||ytilde-y||/||y|| = 2 rho for prescribed c=-rho||s||||y||',
148            'prediction_3_positive_definiteness': 'BFGS corrected update has strictly positive eigenvalues when B is SPD',
149            'sweep': sweep, 'spd_sweep': pd_sweep
150        },
151        'mini_experiment': {'baseline_ordinary_reject': baseline, 'idea_euclidean_correction': idea}
152    }
153    with open('results.json', 'w') as f:
154        json.dump(result, f, indent=2)
155    print(json.dumps(result, indent=2))
156
157
158if __name__ == '__main__':
159    main()