1import json
  2import numpy as np
  3
  4
  5def cubic_tensor(c, g, H, C):
  6    """Return A such that A[(1,s)]^3 = c+g.s+.5*s'Hs+C[s,s,s]/6."""
  7    g, H, C = np.asarray(g), np.asarray(H), np.asarray(C)
  8    n = len(g)
  9    A = np.zeros((n + 1,) * 3, dtype=float)
 10    A[0, 0, 0] = c
 11    for i in range(n):
 12        for p in ((0, 0, i + 1), (0, i + 1, 0), (i + 1, 0, 0)):
 13            A[p] = g[i] / 3.0
 14    for i in range(n):
 15        for j in range(n):
 16            # Three placements of the zero index sum to (1/2) H_ij s_i s_j.
 17            for p in ((0, i + 1, j + 1), (i + 1, 0, j + 1),
 18                      (i + 1, j + 1, 0)):
 19                A[p] = (0.5 * H[i, j]) / 3.0
 20    A[1:, 1:, 1:] = C / 6.0
 21    return A
 22
 23
 24def model(c, g, H, C, s):
 25    return c + g @ s + .5 * s @ H @ s + np.einsum('ijk,i,j,k', C, s, s, s) / 6.
 26
 27
 28def contract(A, u, v):
 29    return np.einsum('ijk,j,k->i', A, u, v)
 30
 31
 32def pam(A, sweeps=30, beta=0.8, seed=0):
 33    rng = np.random.default_rng(seed)
 34    d = A.shape[0]
 35    U = [rng.normal(size=d) for _ in range(3)]
 36    U = [u / np.linalg.norm(u) for u in U]
 37    history = []
 38    def F():
 39        return np.einsum('ijk,i,j,k', A, *U)
 40    for _ in range(sweeps):
 41        for block in range(3):
 42            old = U[block].copy()
 43            before = F()  # proximal term is zero at the old block
 44            others = [U[i] for i in range(3)]
 45            a = contract(A, others[(block + 1) % 3], others[(block + 2) % 3])
 46            z = beta * old - a
 47            if np.linalg.norm(z) < 1e-14:
 48                continue
 49            U[block] = z / np.linalg.norm(z)
 50            after = F() + .5 * beta * np.sum((U[block] - old) ** 2)
 51            history.append((after - before, before, after))
 52    return U, np.asarray(history)
 53
 54
 55def run():
 56    rng = np.random.default_rng(12)
 57    n = 4
 58    c = .37
 59    g = rng.normal(size=n)
 60    H = rng.normal(size=(n, n)); H = (H + H.T) / 2
 61    C0 = rng.normal(size=(n, n, n))
 62    C = sum(C0.transpose(p) for p in ((0,1,2),(1,0,2),(2,1,0),(0,2,1),(1,2,0),(2,0,1))) / 6
 63    A = cubic_tensor(c, g, H, C)
 64
 65    # Prediction 1: exact homogeneous representation.
 66    rep_errors = []
 67    for _ in range(1000):
 68        s = rng.normal(size=n)
 69        rep_errors.append(abs(np.einsum('ijk,i,j,k', A, np.r_[1., s], np.r_[1., s], np.r_[1., s]) - model(c, g, H, C, s)))
 70    representation_max_error = float(max(rep_errors))
 71
 72    # Prediction 2 sweep: exact block solves never increase their proximal objective,
 73    # for every beta; larger beta is the predicted stronger proximal stabilization.
 74    beta_sweep = {}
 75    for beta in (0.0, 0.1, 0.8, 5.0):
 76        _, ch = pam(A, sweeps=40, beta=beta, seed=4)
 77        beta_sweep[str(beta)] = {'max_increase': float(np.max(ch[:, 0])),
 78                                 'final_objective': float(ch[-1, 2])}
 79    max_block_increase = beta_sweep['0.8']['max_increase']
 80    final_F = beta_sweep['0.8']['final_objective']
 81
 82    # Prediction 3 sweep: for a quartic mismatch q*r^4/4, rho decreases with
 83    # radius and its acceptance boundary is the positive root of rho=eta.
 84    sdir = -g / np.linalg.norm(g)
 85    eta = .25
 86    radii = np.geomspace(.03, 2.0, 400)
 87    aa = -(g @ sdir); bb = -.5 * (sdir @ H @ sdir)
 88    dd = -np.einsum('ijk,i,j,k', C, sdir, sdir, sdir) / 6
 89    q_sweep = {}
 90    for q in (.2, .8, 1.7, 4.0):
 91        observed = []
 92        for r in radii:
 93            dec = aa*r + bb*r*r + dd*r**3
 94            rho = (dec - q*r**4/4) / dec
 95            observed.append(rho >= eta)
 96        roots = np.roots([q/4, -(1-eta)*dd, -(1-eta)*bb, -(1-eta)*aa])
 97        pos = sorted(float(z.real) for z in roots if abs(z.imag)<1e-8 and z.real>1e-12)
 98        rstar = pos[0] if pos else float('nan')
 99        accepted = np.where(np.asarray(observed))[0]
100        observed_boundary = float(radii[accepted[-1]]) if len(accepted) else 0.0
101        q_sweep[str(q)] = {'predicted_radius': rstar,
102                           'observed_grid_boundary': observed_boundary,
103                           'relative_grid_error': abs(observed_boundary-rstar)/rstar,
104                           'pattern_matches': bool(np.all(np.asarray(observed)==(radii<=rstar)))}
105    rstar = q_sweep['1.7']['predicted_radius']
106    first_rejected = float(radii[np.where(~np.asarray([((aa*r+bb*r*r+dd*r**3)-1.7*r**4/4)/(aa*r+bb*r*r+dd*r**3)>=eta for r in radii]))[0][0]])
107    transition_matches = q_sweep['1.7']['pattern_matches']
108
109    # Secondary same-objective comparison: fixed-step gradient descent vs guarded cubic-direction steps.
110    def f(x): return model(c, g, H, C, x) + q * np.sum(x**2)**2 / 4
111    def gradient(x): return g + H @ x + .5 * np.einsum('ijk,j,k->i', C, x, x) + q * np.sum(x*x) * x
112    def safeguarded(x, radius):
113        gg = gradient(x); ss = radius * (-gg / (np.linalg.norm(gg) + 1e-12))
114        if model(c, g, H, C, x + ss) <= model(c, g, H, C, x) and f(x + ss) < f(x): return x + ss, True
115        return x, False
116    xg = np.zeros(n); xs = np.zeros(n); lr = .08; rejects = 0
117    for _ in range(80):
118        xg -= lr * gradient(xg)
119        xs, ok = safeguarded(xs, .18); rejects += int(not ok)
120    return {'prediction_1_representation_max_error': representation_max_error,
121            'prediction_2_max_proximal_block_increase': max_block_increase,
122            'prediction_2_final_recorded_block_objective': final_F,
123            'prediction_2_beta_sweep': beta_sweep,
124            'prediction_3_eta': eta, 'prediction_3_quartic_q': 1.7,
125            'prediction_3_predicted_transition_radius': rstar,
126            'prediction_3_first_rejected_grid_radius': first_rejected,
127            'prediction_3_quartic_strength_sweep': q_sweep,
128            'prediction_3_acceptance_patterns_match': transition_matches,
129            'mini_final_objective_gradient': float(f(xg)),
130            'mini_final_objective_safeguarded': float(f(xs)),
131            'mini_safeguarded_rejections': rejects, 'mini_steps': 80}
132
133if __name__ == '__main__':
134    result = run()
135    with open('results.json', 'w') as fp: json.dump(result, fp, indent=2)
136    print(json.dumps(result, indent=2))