1import numpy as np
 2from scipy.linalg import expm
 3
 4
 5def discretize(A, B, dt):
 6    n, k = A.shape[0], B.shape[1]
 7    aug = np.zeros((n + k, n + k)); aug[:n, :n] = A; aug[:n, n:] = B
 8    E = expm(aug * dt)
 9    return E[:n, :n], E[:n, n:]
10
11
12def simulate(m=2, p=2, N=5000, dt=.03, noise=0., seed=0, warmup=1000):
13    rng = np.random.default_rng(seed); n = 2
14    A = np.array([[-.35, .8], [-.6, -.55]])
15    Bbase = np.array([[1.0, .25], [.15, 1.0]])
16    B = Bbase[:, :m]
17    Cbase = np.array([[1., .25], [-.35, 1.], [.4, -.7]])
18    C = Cbase[:p]
19    Af = np.array([[0., 1.], [-2.25, -1.5]]); bf = np.array([[0.], [1.]])
20    Ad, Bd = discretize(A, B, dt); Df, Gf = discretize(Af, bf, dt)
21    x = np.zeros(n); zu = np.zeros((m, n)); zy = np.zeros((p, n)); rows=[]; ys=[]
22    for t in range(N + warmup):
23        u = .8*rng.normal(size=m) + .35*np.sin(.011*t + np.arange(m))
24        y = C @ x + noise*rng.normal(size=p)
25        zu = (Df @ zu.T).T + u[:, None]*Gf[:, 0][None, :]
26        zy = (Df @ zy.T).T + y[:, None]*Gf[:, 0][None, :]
27        x = Ad @ x + Bd @ u
28        if t >= warmup: rows.append(np.r_[zu.ravel(), zy.ravel()]); ys.append(y)
29    return np.asarray(rows), np.asarray(ys)
30
31
32def svd_info(Z, tau=1e-3):
33    s = np.linalg.svd(Z-Z.mean(0), compute_uv=False)
34    return s, int(np.sum(s/s[0] > tau))
35
36
37def features(H):
38    return np.column_stack([H[:,i]*H[:,j] for i in range(H.shape[1]) for j in range(i,H.shape[1])])
39
40
41def condition(X):
42    s=np.linalg.svd(X,compute_uv=False)
43    return float(s[0]/max(s[-1],1e-30))
44
45
46def main():
47    np.set_printoptions(precision=4,suppress=True)
48    print('PREDICTION 1: SIMO 2n=4; MIMO n(m+1)=6')
49    for m,p,label,expected in [(1,2,'SIMO',4),(2,2,'MIMO',6)]:
50        Z,_=simulate(m=m,p=p,seed=3); s,r=svd_info(Z)
51        print(f'{label}: ambient={Z.shape[1]}, predicted={expected}, observed={r}, ratios={s/s[0]}')
52
53    print('\nPREDICTION 2: filter transient decays and rank approaches theorem')
54    for warmup in [0,10,100,300,1000,3000]:
55        Z,_=simulate(m=2,p=2,seed=4,warmup=warmup); s,r=svd_info(Z)
56        print(f'warmup={warmup:4d}, rank(tau=1e-3)={r}, s7/s1={s[6]/s[0]:.3e}, s8/s1={s[7]/s[0]:.3e}')
57
58    print('\nPREDICTION 2: observation-noise tail scales linearly with noise')
59    for noise in [0., 1e-5, 1e-4, 1e-3, 1e-2]:
60        Z,_=simulate(m=2,p=2,seed=5,noise=noise); s,_=svd_info(Z)
61        tail=float(np.mean(s[6:]/s[0]))
62        print(f'noise={noise:.0e}, tail_mean/s1={tail:.3e}, tail/noise={tail/max(noise,1e-30):.3e}')
63
64    print('\nPREDICTION 3: cutoff transition when tau crosses noisy tail')
65    Z,_=simulate(m=2,p=2,seed=5,noise=1e-3); s,_=svd_info(Z)
66    print('singular_ratios=',s/s[0])
67    for tau in [1e-1,1e-2,1e-3,1e-4,1e-5,1e-6,1e-8]:
68        _,r=svd_info(Z,tau); print(f'tau={tau:.0e}, rank={r}')
69
70    print('\nCRITIC REGRESSION: quadratic value target, full vs reduced vs under-compressed')
71    Z,Y=simulate(m=2,p=2,seed=7,noise=1e-3); split=len(Z)//2; Zc=Z-Z[:split].mean(0)
72    U=np.linalg.svd(Zc[:split],full_matrices=False)[2].T
73    target=np.sum(Y**2,axis=1)
74    for r in [3,6,8]:
75        H=Zc@U[:,:r]; X=features(H); tr,te=X[:split],X[split:]
76        beta=np.linalg.solve(tr.T@tr+1e-5*np.eye(tr.shape[1]),tr.T@target[:split])
77        mse=np.mean((te@beta-target[split:])**2)
78        print(f'r={r:2d}, features={X.shape[1]:3d}, test_mse={mse:.5e}, feature_cond={condition(tr):.3e}')
79
80if __name__=='__main__': main()