import json import numpy as np from scipy.linalg import expm, eigvalsh from scipy.optimize import linear_sum_assignment SEED = 1502 rng = np.random.default_rng(SEED) def reachable(MA, MB): n = MA.shape[0] seen = np.any(MB != 0, axis=1) changed = True while changed: changed = False newly = (MA.astype(bool) @ seen.astype(int)) > 0 add = newly & ~seen if np.any(add): seen |= add; changed = True return seen def boolean_core(MA, MB): """Boolean sparsity of [B, AB, ..., A^(n-1)B], retaining each block.""" n, m = MB.shape P = MB.astype(bool) blocks = [P.copy()] for _ in range(1, n): P = (MA.astype(int) @ P.astype(int)) > 0 blocks.append(P.copy()) return np.concatenate(blocks, axis=1) def max_matching_rows(pattern): r, c = pattern.shape if not np.any(pattern): return 0, [] # Add a large penalty for forbidden entries; assignment then maximizes allowed pairs. cost = np.where(pattern, 0.0, 1e6) ri, ci = linear_sum_assignment(cost) pairs = [(int(i), int(j)) for i, j in zip(ri, ci) if pattern[i, j]] return len(pairs), pairs def ctrb(A, B): n = A.shape[0]; P = B.copy(); blocks = [] for _ in range(n): blocks.append(P); P = A @ P return np.concatenate(blocks, axis=1) def gramian(A, B, T=8., steps=240): dt = T / steps; W = np.zeros((A.shape[0], A.shape[0])) for q in range(steps): E = expm(A * ((q + .5) * dt)); X = E @ B W += dt * (X @ X.T) return (W + W.T) / 2 def normalized_design(kind, n=10, m=2): MA = np.zeros((n, n), int); MB = np.zeros((n, m), int) if kind == 'matched': MB[0, 0] = 1; MB[5, 1] = 1 for i in range(1, 5): MA[i, i-1] = 1 for i in range(6, 10): MA[i, i-1] = 1 elif kind == 'unmatched': # Every state is reachable, but all descendants have the same predecessor/input history. MB[0, 0] = 1 for i in range(1, n): MA[i, 0] = 1 elif kind == 'inaccessible': MB[0, 0] = 1 for i in range(1, n-1): MA[i, i-1] = 1 MA[n-1, n-1] = 1 elif kind == 'random': MA = (rng.random((n, n)) < .22).astype(int); np.fill_diagonal(MA, 0) MB = (rng.random((n, m)) < .5).astype(int) if not np.any(MB): MB[0, 0] = 1 return MA, MB def weighted(MA, MB, seed, stable=True): r = np.random.default_rng(seed) A = MA * r.normal(0, .35, MA.shape); B = MB * r.normal(0, 1, MB.shape) rad = max(abs(np.linalg.eigvals(A))) if np.any(A) else 0 if stable and rad > 0: A *= .65 / rad return A, B def rank_tol(M): s = np.linalg.svd(M, compute_uv=False) return int(np.sum(s > 1e-9 * max(s[0], 1e-30))) def structural_report(): out = [] for kind in ['matched', 'unmatched', 'inaccessible', 'random']: MA, MB = normalized_design(kind); core = boolean_core(MA, MB); mm, _ = max_matching_rows(core) out.append({'design': kind, 'reachable': int(reachable(MA, MB).sum()), 'n': len(MA), 'core_nonzero': int(core.sum()), 'matching': mm, 'row_saturating': mm == len(MA)}) return out def gramian_sweep(): rows = [] for kind in ['matched', 'unmatched', 'inaccessible', 'random']: MA, MB = normalized_design(kind); ranks = []; mins = []; ctranks = [] for seed in range(8): A, B = weighted(MA, MB, seed); ctranks.append(rank_tol(ctrb(A, B))) W = gramian(A, B, T=8., steps=180); ev = eigvalsh(W) ranks.append(int(np.sum(ev > 1e-8 * max(ev[-1], 1e-30)))); mins.append(float(ev[0])) rows.append({'design': kind, 'mean_ctrb_rank': float(np.mean(ctranks)), 'mean_gramian_rank': float(np.mean(ranks)), 'min_gramian_rank': min(ranks), 'mean_lambda_min': float(np.mean(mins)), 'median_lambda_min': float(np.median(mins))}) return rows def horizon_sweep(): MA, MB = normalized_design('matched'); A, B = weighted(MA, MB, 7); result = [] for T in [.25, .5, 1, 2, 4, 8]: W = gramian(A, B, T=T, steps=180); ev = eigvalsh(W) result.append({'T': T, 'lambda_min': float(ev[0]), 'lambda_max': float(ev[-1]), 'normalized_lambda_min': float(ev[0] / ev[-1]), 'rank': int(np.sum(ev > 1e-8 * ev[-1]))}) return result def horizon_controllability_sweep(): # Discrete-time analogue: rank grows when additional input history reaches # successive chain states; the unmatched star saturates at its input rank. rows=[] for kind in ["matched", "unmatched", "inaccessible"]: MA,MB=normalized_design(kind); A,B=weighted(MA,MB,7) for H in [1,2,3,4,5,6,8,10]: P=B.copy(); blocks=[] for _ in range(H): blocks.append(P); P=A@P C=np.concatenate(blocks,axis=1) sv=np.linalg.svd(C,compute_uv=False) rows.append({"design":kind,"horizon":H,"rank":rank_tol(C), "min_nonzero_singular":float(sv[rank_tol(C)-1]) if rank_tol(C)>0 else 0.0}) return rows def baseline_comparison(): # Standard sparse random initialization versus matching mask, same n,m and # one fixed linear recurrent setup; report the controllability metric. MA,MB=normalized_design("matched"); RA,RB=normalized_design("random") vals=[] for seed in range(8): for name,ma,mb in [("matching",MA,MB),("random_sparse",RA,RB)]: A,B=weighted(ma,mb,seed); W=gramian(A,B,T=8.,steps=180); ev=eigvalsh(W) vals.append({"design":name,"rank":int(np.sum(ev>1e-8*max(ev[-1],1e-30))), "lambda_min":float(max(ev[0],0.0))}) return vals def main(): result = {'seed': SEED, 'structural': structural_report(), 'gramian_sweep': gramian_sweep(), 'matched_horizon_sweep': horizon_sweep(), 'horizon_controllability_sweep': horizon_controllability_sweep(), 'baseline_comparison': baseline_comparison(), 'predictions': [ 'accessibility plus row-saturating matching predicts generic full controllability rank', 'the matched finite-horizon Gramian has positive minimum eigenvalue for generic weights', 'increasing horizon increases the minimum Gramian eigenvalue for the matched chain, while inaccessible remains singular']} with open('results.json', 'w') as f: json.dump(result, f, indent=2) print(json.dumps(result, indent=2)) if __name__ == '__main__': main()