import json from pathlib import Path import numpy as np SEED = 2844 rng = np.random.default_rng(SEED) OUT = Path('results.json') def relu(x): return np.maximum(x, 0.0) def rollout(A, B, b, x0, U): xs = [x0.copy()] zs = [] x = x0.copy() for u in U: z = A @ x + B @ u + b zs.append(z.copy()) x = relu(z) xs.append(x.copy()) return np.asarray(xs), np.asarray(zs) def observability(A, C, T): blocks = [] P = np.eye(A.shape[0]) for _ in range(T + 1): blocks.append(C @ P) P = P @ A return np.vstack(blocks) def rank_tol(M): s = np.linalg.svd(M, compute_uv=False) return int(np.sum(s > max(M.shape) * np.finfo(float).eps * s[0])) if s[0] else 0 def prediction_positive_regime(): # Positive A, input, bias, and initial state imply every preactivation is positive. n, p, steps = 8, 3, 30 A = 0.82 * np.eye(n) + 0.06 * np.ones((n, n)) / n B = 0.12 * np.ones((n, p)) b = np.full(n, 0.03) x0 = np.full(n, 0.4) U = rng.uniform(0.0, 1.0, size=(steps, p)) xr, z = rollout(A, B, b, x0, U) # Explicit affine recurrence, which must equal ReLU recurrence. xa = [x0.copy()] x = x0.copy() for u in U: x = A @ x + B @ u + b xa.append(x.copy()) xa = np.asarray(xa) return { 'predicted_max_difference': 0.0, 'observed_max_difference': float(np.max(np.abs(xr - xa))), 'negative_preactivation_fraction': float(np.mean(z < 0)), 'pass': bool(np.max(np.abs(xr - xa)) < 1e-12 and np.all(z >= 0)), } def prediction_rank_transition(): # A diagonal system with distinct eigenvalues and one coordinate observed. # For C=[1,0,...], T+1 observations span k coordinates only when T+1 >= k. n = 6 eig = np.array([0.92, 0.78, 0.64, 0.50, 0.36, 0.22]) A = np.diag(eig) C = np.ones((1, n)) # all modes visible; Vandermonde rows become full rank at T=n-1 rows = [] for T in range(0, n + 3): Q = observability(A, C, T) s = np.linalg.svd(Q, compute_uv=False) rows.append({'T': T, 'rank': rank_tol(Q), 'sigma_min': float(s[-1])}) observed_transition = next(r['T'] for r in rows if r['rank'] == n) predicted_transition = n - 1 return { 'predicted_full_rank_first_T': predicted_transition, 'observed_full_rank_first_T': observed_transition, 'sweep': rows, 'pass': observed_transition == predicted_transition, } def prediction_noise_scaling(): # In the linear positive regime y=Qx0+noise. Least-squares error is governed by 1/sigma_min(Q). n = 5 eig = np.array([0.95, 0.80, 0.65, 0.50, 0.35]) A = np.diag(eig) C = np.ones((1, n)) noise = 1e-3 x_trials = rng.normal(size=(1500, n)) result = [] for T in range(n - 1, 13): Q = observability(A, C, T) s = np.linalg.svd(Q, compute_uv=False) # Keep initial conditions positive while preserving a linear identification problem. x = 0.5 + 0.15 * x_trials y = x @ Q.T + rng.normal(scale=noise, size=(len(x), Q.shape[0])) est = np.linalg.lstsq(Q, y.T, rcond=None)[0].T rmse = float(np.sqrt(np.mean((est - x) ** 2))) result.append({'T': T, 'sigma_min': float(s[-1]), 'rmse': rmse, 'rmse_times_sigma': rmse * float(s[-1])}) # Compare the monotonic qualitative prediction and a fitted log slope. sig = np.array([r['sigma_min'] for r in result]) err = np.array([r['rmse'] for r in result]) slope = float(np.polyfit(np.log(sig), np.log(err), 1)[0]) monotonic = bool(np.all(np.diff(sig) > 0) and np.all(np.diff(err) < 0)) return { 'predicted': 'sigma_min increases with T and reconstruction RMSE decreases approximately as 1/sigma_min', 'observed_log_error_vs_log_sigma_slope': slope, 'sweep': result, 'pass': monotonic and slope < -0.5, } def positivity_penalty_sweep(): # A fourth practical check: shifting bias upward removes masks and penalty. n, p, steps = 8, 2, 20 A = 0.65 * np.eye(n) + 0.03 * np.ones((n, n)) B = 0.08 * np.ones((n, p)) x0 = np.full(n, 0.2) U = rng.uniform(-1, 1, size=(steps, p)) out = [] for bias in [-0.20, -0.10, 0.0, 0.10, 0.20]: xs, z = rollout(A, B, np.full(n, bias), x0, U) neg = np.maximum(-z, 0) out.append({'bias': bias, 'negative_fraction': float(np.mean(z < 0)), 'positivity_penalty': float(np.sum(neg * neg))}) return out def baseline_vs_positive(): # Same observable system dimensions and inputs; compare signed ReLU dynamics # against a nonnegative system. The latter is the proposed regime. n, m, steps = 8, 4, 12 C = np.eye(n)[:m] x0 = np.full(n, 0.3) U = rng.uniform(-0.05, 0.05, size=(steps, 2)) cases = { 'baseline_signed_relu': (0.55 * rng.normal(size=(n, n)) / np.sqrt(n), 0.05 * rng.normal(size=(n, 2)), np.zeros(n)), 'positive_regime_relu': (0.55 * (0.5 + rng.random((n, n))) / np.sqrt(n), 0.05 * (0.5 + rng.random((n, 2))), np.full(n, 0.12)), } out = {} for name, (A, B, b) in cases.items(): xs, zs = rollout(A, B, b, x0, U) # Finite differences test the linear perturbation claim along the realized path. d = rng.normal(size=n); d *= 1e-6 xp, _ = rollout(A, B, b, x0 + d, U) actual = (xp[1:] - xs[1:]) predicted = np.asarray([np.linalg.matrix_power(A, t + 1) @ d for t in range(steps)]) rel_err = np.linalg.norm(actual - predicted) / max(np.linalg.norm(actual), 1e-15) Q = observability(A, C, steps - 1) out[name] = { 'negative_preactivation_fraction': float(np.mean(zs < 0)), 'linearized_perturbation_relative_error': float(rel_err), 'sigma_min_Q': float(np.linalg.svd(Q, compute_uv=False)[-1]), 'observability_rank': rank_tol(Q), } return out def main(): results = { 'seed': SEED, 'positive_regime_exactness': prediction_positive_regime(), 'observability_rank_transition': prediction_rank_transition(), 'noise_reconstruction_scaling': prediction_noise_scaling(), 'positivity_bias_sweep': positivity_penalty_sweep(), 'baseline_vs_positive': baseline_vs_positive(), } results['mechanism_pass'] = all([ results['positive_regime_exactness']['pass'], results['observability_rank_transition']['pass'], results['noise_reconstruction_scaling']['pass'], ]) OUT.write_text(json.dumps(results, indent=2)) print(json.dumps(results, indent=2)) if __name__ == '__main__': main()