1"""Toy verification of position-only active-noise cancellation.
  2
  3The quadratic is f(theta)=k theta^2/2.  Observed minibatch gradient is
  4 g=k*theta+a+epsilon, where a is an OU process.  The optimizer estimates a
  5from the residual g-k*theta (position and observed gradient only), then uses
  6 theta <- theta-eta*(g-c*ahat).
  7
  8This is the minimal implementable part of the proposed controller.  The
  9script also computes the exact deterministic augmented transition matrix
 10for [theta,a,ahat] and checks rho<1 against simulated stability.
 11"""
 12import json
 13import numpy as np
 14
 15
 16def kf_run(seed, tau, eta, c, n=12000, k=1.0, sigma_a=1.0,
 17           sigma_g=0.20, burn=2000):
 18    rng = np.random.default_rng(seed)
 19    q = np.exp(-1.0 / tau)
 20    Q = sigma_a**2 * (1.0-q*q)
 21    # Stationary scalar Kalman covariance for measurement a + epsilon.
 22    P = sigma_a**2
 23    h = 0.0
 24    theta = 0.0
 25    xs, hs, aa = [], [], []
 26    a = rng.normal(0.0, sigma_a)
 27    for t in range(n):
 28        a = q*a + rng.normal(0.0, np.sqrt(Q))
 29        eps = rng.normal(0.0, sigma_g)
 30        g = k*theta + a + eps
 31        # Prediction and measurement update using residual g - k*theta.
 32        hp = q*h
 33        Pp = q*q*P + Q
 34        R = sigma_g*sigma_g
 35        gain = Pp/(Pp+R)
 36        h = hp + gain*((g-k*theta)-hp)
 37        P = (1.0-gain)*Pp
 38        theta = theta - eta*(g-c*h)
 39        if t >= burn:
 40            xs.append(theta); hs.append(h); aa.append(a)
 41    return float(np.mean(np.square(xs))), float(np.mean(np.square(np.array(aa)-np.array(hs)))), float(np.mean(np.array(xs)))
 42
 43
 44def transition(eta, k, tau, c):
 45    """Noise-free expected transition for [theta,a,ahat]."""
 46    q = np.exp(-1.0/tau)
 47    # In the zero-observation-noise filter, the update gain is 1.
 48    # For the actual noisy Kalman filter use its stationary gain below.
 49    R = .20**2
 50    Q = 1-q*q
 51    P = 1.0
 52    # converge scalar Riccati to get fixed gain
 53    for _ in range(10000):
 54        Pp=q*q*P+Q
 55        L=Pp/(Pp+R)
 56        P=(1-L)*Pp
 57    # theta' = (1-eta*k)theta -eta*a + eta*c*h
 58    # h' = L*a + q*(1-L)*h (residual is a in expectation)
 59    return np.array([[1-eta*k, -eta + eta*c*L, eta*c*q*(1-L)],
 60                     [0, q, 0],
 61                     [0, L, q*(1-L)]], dtype=float)
 62
 63
 64def feedback_transition(dt, k, G, Ga, tau, L):
 65    """A_cl for e=theta-lambda, a, ahat under the stated lambda feedback."""
 66    q=np.exp(-dt/tau)
 67    return np.array([[1+dt*(G-k), dt, dt*Ga],
 68                     [0,q,0], [0,L, q*(1-L)]], float)
 69
 70
 71def empirical_boundary():
 72    # Deterministic trajectories establish the observed boundary, while the
 73    # matrix supplies the claimed prediction.  Noise is deliberately absent.
 74    tau, k, c = 5.0, 1.0, 0.7
 75    last_stable = None; first_unstable = None
 76    for eta in np.linspace(.01, 2.5, 500):
 77        A=transition(eta,k,tau,c); rho=max(abs(np.linalg.eigvals(A)))
 78        x=np.array([1., 1., 0.])
 79        for _ in range(300): x=A@x
 80        stable=(rho < 1.0)
 81        if stable: last_stable=(eta,rho)
 82        elif first_unstable is None: first_unstable=(eta,rho)
 83    return last_stable, first_unstable
 84
 85
 86def main():
 87    taus=[0.2,1.0,5.0,20.0]
 88    etas=[0.15,0.30,0.50]
 89    rows=[]
 90    for tau in taus:
 91        for eta in etas:
 92            base=kf_run(10,tau,eta,0.0)
 93            idea=kf_run(10,tau,eta,0.85)
 94            rows.append({'tau':tau,'eta':eta,'sgd_mse':base[0],
 95                         'controller_mse':idea[0],
 96                         'ratio':idea[0]/base[0],
 97                         'residual_mse':idea[1]})
 98    lb, fu=empirical_boundary()
 99    # A representative feedback-loop stability check, matching the formula.
100    fb=[]
101    for G in [0.1,0.5,1.0,2.0]:
102        A=feedback_transition(0.1,1.0,G,0.5,5.0,0.7)
103        fb.append({'G':G,'rho':float(max(abs(np.linalg.eigvals(A))))})
104    out={'rows':rows,'predicted_vs_empirical_boundary':{
105        'last_stable_eta_rho':lb,'first_unstable_eta_rho':fu},
106        'feedback_rho_scan':fb}
107    print(json.dumps(out,indent=2))
108
109if __name__=='__main__': main()