1import json
  2import math
  3import numpy as np
  4
  5
  6def decimated(A, Q, d):
  7    n = A.shape[0]
  8    Ad = np.eye(n)
  9    Qd = np.zeros_like(Q)
 10    # A^j Q (A^j)^T, j=0,...,d-1
 11    Aj = np.eye(n)
 12    for j in range(d):
 13        Qd += Aj @ Q @ Aj.T
 14        Aj = A @ Aj
 15    return Aj, (Qd + Qd.T) * 0.5
 16
 17
 18def riccati(A, H, Q, R, d, max_iter=10000, tol=1e-11):
 19    Ad, Qd = decimated(A, Q, d)
 20    P = Qd.copy()
 21    for _ in range(max_iter):
 22        S = H @ P @ H.T + R
 23        # Joseph-equivalent form of the stated prior Riccati update
 24        K = np.linalg.solve(S, H @ P @ Ad.T).T
 25        Pn = Ad @ P @ Ad.T - K @ S @ K.T + Qd
 26        Pn = (Pn + Pn.T) * 0.5
 27        if np.max(np.abs(Pn - P)) < tol:
 28            P = Pn
 29            break
 30        P = Pn
 31    return P, Ad, Qd
 32
 33
 34def choose_stride(A, H, Q, R, tau, dmax):
 35    vals = []
 36    for d in range(1, dmax + 1):
 37        P, _, _ = riccati(A, H, Q, R, d)
 38        vals.append(float(np.linalg.eigvalsh(P).max()))
 39    feasible = [d + 1 for d, v in enumerate(vals) if v <= tau]
 40    return (max(feasible) if feasible else 1), vals
 41
 42
 43def math_check(seed=7):
 44    rng = np.random.default_rng(seed)
 45    # Stable, observable 2-D system; empirical covariance is measured over many
 46    # independent d-step prior rollouts initialized from the filtered covariance.
 47    A = np.array([[0.88, 0.12], [-0.04, 0.82]])
 48    H = np.eye(2)
 49    Q = np.diag([0.015, 0.008])
 50    R = np.diag([0.025, 0.025])
 51    P1, _, _ = riccati(A, H, Q, R, 1)
 52    rows = []
 53    for d in (1, 2, 4, 8):
 54        P, Ad, Qd = riccati(A, H, Q, R, d)
 55        # P is the steady-state prior at an observation interval of d.
 56        # The d-step prior starts immediately after an observation, hence
 57        # sample from the corresponding posterior covariance P_plus.
 58        S = H @ P @ H.T + R
 59        K = P @ H.T @ np.linalg.inv(S)
 60        P_plus = (P - K @ S @ K.T + (P - K @ S @ K.T).T) * 0.5
 61        x = rng.multivariate_normal(np.zeros(2), P_plus, size=60000)
 62        noise = rng.multivariate_normal(np.zeros(2), Qd, size=60000)
 63        e = (x @ Ad.T) + noise
 64        empirical = np.cov(e, rowvar=False, bias=True)
 65        pred = float(np.linalg.eigvalsh(P).max())
 66        emp = float(np.linalg.eigvalsh(empirical).max())
 67        rows.append({"d": d, "predicted_lambda_max": pred,
 68                     "empirical_lambda_max": emp,
 69                     "relative_error": abs(pred - emp) / max(pred, 1e-12)})
 70    # monotonic signal expected for this stable system
 71    monotone = all(rows[i]["predicted_lambda_max"] <= rows[i+1]["predicted_lambda_max"] + 1e-10
 72                   for i in range(len(rows)-1))
 73    tau = 0.13
 74    selected = choose_stride(A, H, Q, R, tau, 8)[0]
 75    return {"rows": rows, "predicted_monotone": monotone,
 76            "tau": tau, "selected_stride": selected,
 77            "max_relative_error": max(r["relative_error"] for r in rows)}
 78
 79
 80class ToyWorld:
 81    # Mild nonlinear latent dynamics; observation is the latent vector plus noise.
 82    def __init__(self, A, q=0.006, r=0.04, seed=0):
 83        self.A = A
 84        self.q = q
 85        self.r = r
 86        self.rng = np.random.default_rng(seed)
 87
 88    def step(self, z):
 89        return self.A @ z + np.array([0.035 * np.tanh(z[1]),
 90                                      -0.025 * np.tanh(z[0])]) + self.rng.normal(0, math.sqrt(self.q), 2)
 91
 92    def observe(self, z):
 93        return z + self.rng.normal(0, math.sqrt(self.r), 2)
 94
 95
 96def rollout(strategy, A, H, Q, R, tau, seed, T=100, d_fixed=4):
 97    # A scheduler decides when to read an expensive observation. Between reads,
 98    # the model rolls forward using its local linear dynamics (here the known toy model).
 99    world = ToyWorld(A, q=float(Q[0, 0]), r=float(R[0, 0]), seed=seed)
100    true = np.zeros(2)
101    estimate = np.zeros(2)
102    P = np.eye(2) * 0.05
103    calls = 0
104    sqerr = []
105    t = 0
106    strides = []
107    while t < T:
108        if strategy == "adaptive":
109            d, _ = choose_stride(A, H, Q, R, tau, 8)
110        else:
111            d = d_fixed
112        d = min(d, T - t)
113        # Observation update at current time (the expensive encoder call).
114        y = world.observe(true)
115        calls += 1
116        S = H @ P @ H.T + R
117        K = P @ H.T @ np.linalg.inv(S)
118        estimate = estimate + K @ (y - H @ estimate)
119        P = (np.eye(2) - K @ H) @ P
120        # Predict d transitions, recording prediction error at every frame.
121        for _ in range(d):
122            true = world.step(true)
123            estimate = A @ estimate + np.array([0.035 * np.tanh(estimate[1]),
124                                                -0.025 * np.tanh(estimate[0])])
125            P = A @ P @ A.T + Q
126            sqerr.append(float(np.sum((estimate - true) ** 2)))
127            t += 1
128        strides.append(d)
129    return {"rmse": float(np.sqrt(np.mean(sqerr))), "encoder_calls": calls,
130            "mean_stride": float(np.mean(strides)), "strides": strides}
131
132
133def experiment():
134    A = np.array([[0.88, 0.12], [-0.04, 0.82]])
135    H = np.eye(2)
136    Q = np.diag([0.006, 0.006])
137    R = np.diag([0.04, 0.04])
138    # Chosen from the analytical covariance curve: d=4 is feasible, d=5 is not.
139    tau = 0.105
140    check = math_check()
141    out = {"math_check": check, "runs": {}}
142    for strategy, d in (("baseline_stride_1", 1), ("fixed_stride_4", 4), ("fixed_stride_8", 8), ("adaptive", 0)):
143        vals = [rollout(strategy, A, H, Q, R, tau, seed=100+i, T=240, d_fixed=d)
144                for i in range(12)]
145        out["runs"][strategy] = {
146            "rmse_mean": float(np.mean([v["rmse"] for v in vals])),
147            "rmse_std": float(np.std([v["rmse"] for v in vals])),
148            "calls_mean": float(np.mean([v["encoder_calls"] for v in vals])),
149            "mean_stride": float(np.mean([v["mean_stride"] for v in vals])),
150            "example_strides": vals[0]["strides"][:12]
151        }
152    # Claimed stability signal: increasing the dominant eigenvalue increases
153    # predicted uncertainty and can force a shorter admissible interval.
154    probe = []
155    for radius in (0.90, 0.95, 0.98, 0.99, 0.995):
156        Ap = np.diag([radius, 0.75])
157        selected, curve = choose_stride(Ap, H, Q, R, tau, 16)
158        probe.append({"dominant_radius": radius, "selected_stride": selected,
159                      "lambda_max_by_d": curve})
160    out["near_unit_radius_probe"] = probe
161    return out
162
163
164if __name__ == "__main__":
165    print(json.dumps(experiment(), indent=2))