1import json, time
 2import numpy as np
 3from active_set_router import active_set_simplex, projected_gradient, objective
 4
 5np.set_printoptions(precision=4, suppress=True)
 6rng = np.random.default_rng(1263)
 7E = 16
 8# Orthogonal feature construction makes the exact Gram spectrum known while
 9# retaining matrix-free Z/Z.T products in the solver.
10Q, _ = np.linalg.qr(rng.normal(size=(E,E)))
11U, _ = np.linalg.qr(rng.normal(size=(E,E)))
12xtrue = np.zeros(E); xtrue[[0,3,7,12]] = [0.05, 0.15, 0.30, 0.50]
13rows=[]
14for lam in [1e-4, 1e-3, 1e-2, 1e-1, 1.0]:
15    # eigenvalues of Z'Z range from 1 to 1e4; ridge controls kappa
16    eig = np.geomspace(1.0, 1e4, E)
17    Z = Q @ np.diag(np.sqrt(eig)) @ U.T
18    A = Z.T@Z + lam*np.eye(E)
19    b = A@xtrue
20    cond = np.linalg.eigvalsh(A)[-1]/np.linalg.eigvalsh(A)[0]
21    t0=time.perf_counter(); xa, ma = active_set_simplex(Z,b,lam=lam,cg_tol=1e-11,max_pivots=300); ta=time.perf_counter()-t0
22    t0=time.perf_counter(); xp, pg = projected_gradient(Z,b,lam=lam,steps=100000,tol=1e-10); tp=time.perf_counter()-t0
23    # KKT residual after choosing equality multiplier from free coordinates
24    g0 = A@xa-b
25    free = xa > 1e-8
26    nu = -np.mean(g0[free]) if np.any(free) else 0.0
27    station = np.max(np.abs((g0+nu)[free])) if np.any(free) else np.nan
28    dual = np.min((g0+nu)[~free]) if np.any(~free) else 0.0
29    rows.append(dict(lambda_=lam, condition=cond, sqrt_condition=np.sqrt(cond),
30                     active_cg_matvecs=ma['cg_matvecs'], pivots=ma['pivots'],
31                     active_free=ma['free'], pg_steps=pg,
32                     active_objective=objective(Z,b,xa,lam), pg_objective=objective(Z,b,xp,lam),
33                     simplex_error=float(abs(xa.sum()-1)), min_x=float(xa.min()),
34                     kkt_free=float(station), kkt_bound_min=float(dual),
35                     active_time=ta, pg_time=tp))
36# Separate spectrum sweep: predicted CG dependence on sqrt(kappa), using the
37# same simplex problem but varying the raw Gram dynamic range.
38scaling=[]
39lam=1e-3
40for ratio in [1, 10, 100, 1000, 10000]:
41    eig=np.geomspace(1., ratio, E)
42    Z=Q@np.diag(np.sqrt(eig))@U.T
43    A=Z.T@Z+lam*np.eye(E); b=A@xtrue
44    cond=np.linalg.eigvalsh(A)[-1]/np.linalg.eigvalsh(A)[0]
45    xa,ma=active_set_simplex(Z,b,lam=lam,cg_tol=1e-11,max_pivots=300)
46    xp,pg=projected_gradient(Z,b,lam=lam,steps=100000,tol=1e-10)
47    scaling.append(dict(raw_ratio=ratio, condition=cond, sqrt_condition=np.sqrt(cond),
48                        cg_matvecs=ma['cg_matvecs'], pg_steps=pg,
49                        objective_gap=objective(Z,b,xp,lam)-objective(Z,b,xa,lam)))
50result={'ridge_sweep':rows,'spectrum_sweep':scaling,
51        'predictions':[
52          'Increasing ridge decreases kappa(A) monotonically and should reduce CG work.',
53          'Across spectra, CG work should track sqrt(kappa) more closely than projected-gradient steps track kappa.',
54          'The active-set output should have exact simplex feasibility and nonnegative coefficients, with sparse support.'
55        ]}
56with open('results.json','w') as f: json.dump(result,f,indent=2)
57print(json.dumps(result,indent=2))