1"""Forward-intersection spectral latent dynamics: numerical verification.
  2
  3Toy: G=span(e0,...,e3) in R^5.  e0/e1/e2 are supported stable modes;
  4e3 is an unstable nuisance mode whose image leaks into e4.  We compare
  5exact nullspace intersections with a finite-tolerance principal-angle filter.
  6"""
  7import json
  8import numpy as np
  9
 10
 11def orth(A, tol=1e-10):
 12    if A.size == 0:
 13        return np.zeros((A.shape[0], 0))
 14    u, s, _ = np.linalg.svd(A, full_matrices=False)
 15    scale = s[0] if len(s) else 1.0
 16    return u[:, :int(np.sum(s > tol * scale))]
 17
 18
 19def intersection_basis(Q, R, tol=1e-8):
 20    """Basis for range(Q) intersect range(R) from ker([Q,-R])."""
 21    p, q = Q.shape[1], R.shape[1]
 22    _, s, vh = np.linalg.svd(np.concatenate([Q, -R], axis=1), full_matrices=True)
 23    scale = s[0] if len(s) else 1.0
 24    rank = int(np.sum(s > tol * scale))
 25    return orth(Q @ vh[rank:, :].T[:p, :], tol)
 26
 27
 28def principal_filter(Q, R, tau):
 29    """Approximate intersection: retain Q directions with sin(angle)<=tau."""
 30    # Principal angles require orthonormal bases on both sides.
 31    Rorth = orth(R)
 32    U, c, _ = np.linalg.svd(Q.T @ Rorth, full_matrices=False)
 33    keep = c >= np.sqrt(max(0., 1. - tau * tau))
 34    return Q @ U[:, keep], c
 35
 36
 37def make_K(alpha):
 38    K = np.zeros((5, 5))
 39    K[0, 0], K[1, 1], K[2, 2] = .90, .70, .80
 40    K[3, 3] = 1.30                 # unsupported unstable mode
 41    K[4, 3] = alpha                # forward incompatibility/leakage
 42    K[4, 4] = 1.0
 43    return K
 44
 45
 46def run(seed=7):
 47    np.random.default_rng(seed)       # fixed-seed API; this toy is deterministic
 48    Q0 = np.eye(5)[:, :4]
 49    tau = .20
 50    predicted_boundary = 1.30 * tau / np.sqrt(1 - tau * tau)
 51    alphas = np.array([0., .10, .20, .25, .263, .27, .35, .60, 1.0])
 52    rows = []
 53    for alpha in alphas:
 54        K = make_K(alpha)
 55        Q1, cosines = principal_filter(Q0, K @ Q0, tau)
 56        Af, A0 = Q1.T @ K @ Q1, Q0.T @ K @ Q0
 57        ef, e0 = np.linalg.eigvals(Af), np.linalg.eigvals(A0)
 58        supported = [.90, .70, .80]
 59        support_err = max(min(abs(z - x) for z in ef) for x in supported)
 60        sin_angle = alpha / np.sqrt(1.30**2 + alpha**2)
 61        rows.append({
 62            "alpha": float(alpha),
 63            "predicted_nuisance_retained": bool(alpha <= predicted_boundary),
 64            "observed_nuisance_retained": bool(np.any(abs(ef - 1.30) < 1e-6)),
 65            "predicted_dim": 4 if alpha <= predicted_boundary else 3,
 66            "observed_dim": int(Q1.shape[1]),
 67            "sin_principal_angle": float(sin_angle),
 68            "spectral_radius_baseline": float(max(abs(e0))),
 69            "spectral_radius_filtered": float(max(abs(ef)) if len(ef) else 0.),
 70            "supported_eigenvalue_error": float(support_err),
 71            "cosines": [float(x) for x in cosines],
 72        })
 73
 74    exact_dims = {str(a): int(intersection_basis(Q0, make_K(a) @ Q0).shape[1])
 75                  for a in [0., .1, .5]}
 76
 77    # Equal-step rollout from a nuisance-contaminated latent state.
 78    alpha = .60
 79    K = make_K(alpha)
 80    Q1, _ = principal_filter(Q0, K @ Q0, tau)
 81    A0, Af = Q0.T @ K @ Q0, Q1.T @ K @ Q1
 82    z0 = np.array([1., -1., .5, 1.])
 83    base, filt = z0.copy(), Q1.T @ (Q0 @ z0)
 84    for _ in range(30):
 85        base, filt = A0 @ base, Af @ filt
 86    rollout = {"steps": 30, "baseline_norm": float(np.linalg.norm(base)),
 87               "filtered_norm": float(np.linalg.norm(filt)),
 88               "norm_ratio_filtered_over_baseline": float(np.linalg.norm(filt) / np.linalg.norm(base))}
 89
 90    observed = {
 91        "P1_boundary_observed_transition": "see rows; dimension/nuisance switch",
 92        "P2_max_supported_error": max(r["supported_eigenvalue_error"] for r in rows),
 93        "P3_exact_dimensions": exact_dims, "rollout": rollout}
 94    passed = (observed["P2_max_supported_error"] < 1e-10
 95              and exact_dims == {"0.0": 4, "0.1": 3, "0.5": 3}
 96              and all(r["observed_dim"] == r["predicted_dim"] for r in rows))
 97    return {"tau": tau,
 98            "predicted": {"P1_boundary_alpha": predicted_boundary,
 99                           "P2_supported_eigenpairs_preserved_error": "0 (floating point)",
100                           "P3_exact_intersection_dimension": "4 at alpha=0; 3 for any nonzero alpha"},
101            "observed": observed, "rows": rows, "pass": bool(passed)}
102
103
104if __name__ == "__main__":
105    print(json.dumps(run(), indent=2))