1import json
  2import numpy as np
  3
  4SEED = 2759
  5rng = np.random.default_rng(SEED)
  6
  7
  8def make_problem(gamma):
  9    # H = I + gamma*(ones-I): SPD for gamma < 1, while Jacobi rho=2*gamma.
 10    E = gamma * (np.ones((3, 3)) - np.eye(3))
 11    H = np.eye(3) + E
 12    r = np.array([1.0, -0.7, 0.45])
 13    return H, E, r
 14
 15
 16def spectral_radius(A):
 17    return float(np.max(np.abs(np.linalg.eigvals(A))))
 18
 19
 20def run_jacobi(H, E, r, steps, alpha=1.0):
 21    z = np.zeros(3)
 22    residuals = []
 23    Dinv = np.diag(1.0 / np.diag(H))
 24    M = (1-alpha)*np.eye(3) - alpha*Dinv@E
 25    for _ in range(steps):
 26        z = (1-alpha)*z + alpha*Dinv@(r-E@z)
 27        residuals.append(float(np.linalg.norm(H@z-r)))
 28    return z, np.array(residuals), spectral_radius(M)
 29
 30
 31def run_gs(H, r, steps):
 32    z = np.zeros(3)
 33    residuals = []
 34    for _ in range(steps):
 35        for i in range(3):
 36            z[i] = (r[i] - (H[i] @ z - H[i, i]*z[i])) / H[i, i]
 37        residuals.append(float(np.linalg.norm(H@z-r)))
 38    return z, np.array(residuals)
 39
 40
 41def run_gd(H, r, steps, lr=0.25):
 42    z = np.zeros(3)
 43    residuals = []
 44    for _ in range(steps):
 45        z -= lr*(H@z-r)
 46        residuals.append(float(np.linalg.norm(H@z-r)))
 47    return z, np.array(residuals)
 48
 49
 50def run_adam(H, r, steps, lr=0.15):
 51    z = np.zeros(3); m = np.zeros(3); v = np.zeros(3)
 52    residuals = []
 53    for t in range(1, steps+1):
 54        g = H@z-r
 55        m = .9*m + .1*g; v = .999*v + .001*g*g
 56        z -= lr*(m/(1-.9**t))/(np.sqrt(v/(1-.999**t))+1e-8)
 57        residuals.append(float(np.linalg.norm(H@z-r)))
 58    return z, np.array(residuals)
 59
 60
 61def main():
 62    gammas = [0.10, 0.30, 0.49, 0.51, 0.70, 0.90]
 63    rows = []
 64    for g in gammas:
 65        H,E,r = make_problem(g)
 66        rho = spectral_radius(E)
 67        _, res, rho_measured = run_jacobi(H,E,r,40)
 68        # empirical asymptotic ratio, avoiding initial transient
 69        ratio = float(np.median(res[-8:-1]/res[-9:-2]))
 70        _, damp_res, damp_rho = run_jacobi(H,E,r,200,alpha=.5)
 71        _, gs_res = run_gs(H,r,40)
 72        _, gd_res = run_gd(H,r,40,lr=.25)
 73        _, adam_res = run_adam(H,r,40)
 74        if rho < 0.8:
 75            policy, _, policy_res, policy_rho = 'parallel', *run_jacobi(H,E,r,40,alpha=1.0)
 76        elif rho < 1.0:
 77            policy, _, policy_res, policy_rho = 'damped', *run_jacobi(H,E,r,40,alpha=.5)
 78        else:
 79            policy, _, policy_res = 'sequential', *run_gs(H,r,40)
 80            policy_rho = None
 81        rows.append({
 82            'gamma': g, 'predicted_rho_2gamma': 2*g, 'policy': policy, 'policy_residual_40': float(policy_res[-1]), 'measured_rho': rho,
 83            'jacobi_matrix_rho': rho_measured, 'jacobi_ratio_last': ratio,
 84            'jacobi_residual_40': float(res[-1]), 'damped_alpha_.5_rho': damp_rho,
 85            'damped_residual_200': float(damp_res[-1]), 'gs_residual_40': float(gs_res[-1]),
 86            'gd_residual_40': float(gd_res[-1]), 'adam_residual_40': float(adam_res[-1]),
 87            'spd_min_eigenvalue': float(np.min(np.linalg.eigvalsh(H)))
 88        })
 89
 90    boundary = []
 91    for g in np.arange(0.45, 0.551, 0.01):
 92        H,E,r = make_problem(float(g))
 93        _, rr, mr = run_jacobi(H,E,r,100)
 94        boundary.append({'gamma': float(g), 'predicted_rho': float(2*g),
 95                         'matrix_rho': mr, 'residual_growth_100': float(rr[-1]/max(rr[0],1e-300))})
 96
 97    # Direct numerical checks of the three claims.
 98    stable_below = [x for x in rows if x['gamma'] < .5 and x['jacobi_residual_40'] < 1e-6]
 99    unstable_above = [x for x in rows if x['gamma'] > .5 and x['jacobi_residual_40'] > 1e2]
100    ratio_check = [x for x in rows if .2 <= x['gamma'] <= .49 and abs(x['jacobi_ratio_last']-x['predicted_rho_2gamma']) < .03]
101    damp_check = [x for x in rows if x['gamma'] == .9 and x['damped_alpha_.5_rho'] < 1 and x['damped_residual_200'] < 0.2]
102    report = {
103        'seed': SEED,
104        'predictions': {
105            'boundary': 'rho=2*gamma, transition at gamma=0.5',
106            'contraction': 'stable Jacobi residual ratio approaches rho',
107            'damping': 'alpha=0.5 stabilizes gamma=0.9 despite undamped rho=1.8'
108        },
109        'rows': rows, 'boundary_sweep': boundary,
110        'checks': {
111            'stable_below_count': len(stable_below), 'unstable_above_count': len(unstable_above),
112            'ratio_match_count': len(ratio_check), 'damping_check': bool(damp_check)
113        }
114    }
115    with open('results.json','w') as f: json.dump(report,f,indent=2)
116    print(json.dumps(report, indent=2))
117
118if __name__ == '__main__':
119    main()