1import json
  2import numpy as np
  3from scipy.optimize import minimize_scalar
  4
  5
  6def bdd_q(M, dims):
  7    starts = np.cumsum([0] + list(dims))
  8    qs = []
  9    for i in range(len(dims)):
 10        a, b = starts[i], starts[i + 1]
 11        D = M[a:b, a:b]
 12        off = np.concatenate([M[a:b, starts[j]:starts[j+1]]
 13                              for j in range(len(dims)) if j != i], axis=1)
 14        try:
 15            z = np.linalg.solve(D, off)
 16            qs.append(float(np.max(np.sum(np.abs(z), axis=1))))
 17        except np.linalg.LinAlgError:
 18            qs.append(float('inf'))
 19    return max(qs), qs
 20
 21
 22def fixed_point(r, b, max_steps=2000, tol=1e-8):
 23    x = np.zeros(2)
 24    for k in range(1, max_steps + 1):
 25        xn = b + r * np.array([x[1], x[0]])
 26        if np.max(np.abs(xn - x)) < tol:
 27            return True, k, xn
 28        x = xn
 29    return False, max_steps, x
 30
 31
 32def toy_sweep():
 33    rows = []
 34    b = np.array([1.0, -0.4])
 35    for r in np.linspace(0, 1.2, 49):
 36        M = np.array([[1., -r], [-r, 1.]])
 37        q, block_q = bdd_q(M, [1, 1])
 38        sv = np.linalg.svd(M, compute_uv=False)
 39        cond = float(sv[0] / sv[-1]) if sv[-1] > 1e-14 else float('inf')
 40        ok, steps, _ = fixed_point(r, b)
 41        rows.append(dict(r=float(r), q=float(q), block_q=block_q,
 42                         determinant=float(np.linalg.det(M)), cond2=cond,
 43                         converged=ok, iterations=steps))
 44    below = [z['r'] for z in rows if z['converged']]
 45    above = [z['r'] for z in rows if not z['converged']]
 46    fit = [(z['q'], z['cond2']) for z in rows if z['q'] <= .9]
 47    scaled = [c * (1-q) for q, c in fit if np.isfinite(c)]
 48    return {
 49        'rows': rows,
 50        'predictions': {
 51            'q_equals_coupling_max_abs_error': max(abs(z['q'] - z['r']) for z in rows),
 52            'predicted_singularity_q': 1.0,
 53            'observed_zero_determinant_r': next(z['r'] for z in rows if abs(z['determinant']) < 1e-12),
 54            'observed_fixed_point_boundary_midpoint': (max(below) + min(above)) / 2,
 55            'predicted_fixed_point_boundary_q': 1.0,
 56            'condition_scaled_c_times_1_minus_q_median_q_le_0.9': float(np.median(scaled)),
 57            'condition_scaled_range_q_le_0.9': [float(min(scaled)), float(max(scaled))]
 58        }
 59    }
 60
 61
 62def constrained_fit(seed=7, lam=0.0, delta=.2):
 63    rng = np.random.default_rng(seed)
 64    r_true = .95
 65    b = np.array([1.0, -.4])
 66    y = np.linalg.solve(np.array([[1., -r_true], [-r_true, 1.]]), b)
 67    y = y + .001 * rng.normal(size=2)
 68
 69    def parts(r):
 70        M = np.array([[1., -r], [-r, 1.]])
 71        pred = np.linalg.solve(M, b)
 72        task = .5 * np.sum((pred-y)**2)
 73        q = abs(r)
 74        penalty = lam * max(0., q-(1-delta))**2
 75        return task + penalty, task, q
 76
 77    opt = minimize_scalar(lambda r: parts(r)[0], bounds=(-.99, .99), method='bounded',
 78                          options={'xatol': 1e-12, 'maxiter': 1000})
 79    total, task, q = parts(opt.x)
 80    return {'r_learned': float(opt.x), 'q': float(q), 'task_loss': float(task),
 81            'total_loss': float(total), 'target_q': r_true,
 82            'certificate_limit': 1-delta, 'optimizer_success': bool(opt.success)}
 83
 84
 85def random_certificate_check(seed=19, n=1000):
 86    rng = np.random.default_rng(seed)
 87    worst_q = 0.0
 88    violations = 0
 89    for _ in range(n):
 90        # Two 2x2 diagonal blocks, with off-diagonal blocks scaled to q<0.9.
 91        D1 = np.eye(2) + 0.2 * rng.normal(size=(2, 2))
 92        D2 = np.eye(2) + 0.2 * rng.normal(size=(2, 2))
 93        C12 = rng.normal(size=(2, 2))
 94        C21 = rng.normal(size=(2, 2))
 95        # Scale both cross blocks until the exact BDD q is below 0.9.
 96        M = np.block([[D1, C12], [C21, D2]])
 97        q0, _ = bdd_q(M, [2, 2])
 98        scale = 0.85 / max(q0, 1e-12)
 99        M = np.block([[D1, scale*C12], [scale*C21, D2]])
100        q, _ = bdd_q(M, [2, 2])
101        worst_q = max(worst_q, q)
102        if abs(np.linalg.det(M)) < 1e-8:
103            violations += 1
104    return {'samples': n, 'worst_q': float(worst_q),
105            'near_singular_count': violations}
106
107
108def main():
109    out = {'toy': toy_sweep(),
110           'random_block_certificate': random_certificate_check(),
111           'fit_unconstrained': constrained_fit(lam=0.0),
112           'fit_bdd_lam_10': constrained_fit(lam=10.0),
113           'fit_bdd_lam_1000': constrained_fit(lam=1000.0),
114           'fit_bdd_lam_100000': constrained_fit(lam=100000.0)}
115    with open('results.json', 'w') as f:
116        json.dump(out, f, indent=2)
117    p = out['toy']['predictions']
118    print('q/coupling max error:', p['q_equals_coupling_max_abs_error'])
119    print('observed singularity r:', p['observed_zero_determinant_r'],
120          'predicted:', p['predicted_singularity_q'])
121    print('observed fixed-point boundary:', p['observed_fixed_point_boundary_midpoint'],
122          'predicted:', p['predicted_fixed_point_boundary_q'])
123    print('median cond*(1-q):', p['condition_scaled_c_times_1_minus_q_median_q_le_0.9'])
124    for k in ('fit_unconstrained', 'fit_bdd_lam_10', 'fit_bdd_lam_1000'):
125        print(k + ':', out[k])
126
127
128if __name__ == '__main__':
129    main()
130