Tangential Bellman Tie Resolver / tie_resolver_experiment.py
Failed on benchmark
1import json
2import numpy as np
3from pathlib import Path
4
5SEED = 2137
6rng = np.random.default_rng(SEED)
7
8
9def H(z, delta, P, beta):
10 # z, delta: [B,S]; P: [B,S,S]
11 continuation = np.min(z + delta, axis=0)
12 return delta + beta * np.einsum("bxy,y->bx", P, continuation)
13
14
15def solve(delta, P, beta, tol=1e-13, max_iter=10000):
16 z = np.zeros_like(delta)
17 for k in range(max_iter):
18 zn = H(z, delta, P, beta)
19 if np.max(np.abs(zn-z)) < tol:
20 return zn, k + 1
21 z = zn
22 raise RuntimeError("fixed point did not converge")
23
24
25def quantile_width(xs, ps, lo=.1, hi=.9):
26 # Width between interpolated crossings; NaN if sweep is too narrow.
27 def cross(q):
28 ind = np.where((ps[:-1]-q) * (ps[1:]-q) <= 0)[0]
29 if len(ind) == 0:
30 return np.nan
31 i = ind[0]
32 if ps[i+1] == ps[i]: return xs[i]
33 return xs[i] + (q-ps[i])*(xs[i+1]-xs[i])/(ps[i+1]-ps[i])
34 return abs(cross(hi)-cross(lo))
35
36
37def contraction_check():
38 B, S, beta = 4, 7, .73
39 raw = rng.random((B,S,S)); P = raw/raw.sum(axis=2, keepdims=True)
40 delta = rng.random((B,S))
41 z = rng.normal(size=(B,S)); w = rng.normal(size=(B,S))
42 ratio = np.max(np.abs(H(z,delta,P,beta)-H(w,delta,P,beta))) / np.max(np.abs(z-w))
43 # Use a generic instance to avoid accidental finite-state exactness.
44 zstar, _ = solve(delta,P,beta)
45 zk = np.zeros_like(delta)
46 errors = []
47 bounds = []
48 initial = np.max(np.abs(zk-zstar))
49 for k in range(9):
50 errors.append(float(np.max(np.abs(zk-zstar))))
51 bounds.append(float(beta**k * initial))
52 zk = H(zk,delta,P,beta)
53 positive = [errors[i+1]/errors[i] for i in range(8) if errors[i] > 1e-14]
54 beta_sweep=[]
55 for bb in [.2, .5, .8, .95]:
56 zz,_=solve(delta,P,bb); cur=np.zeros_like(delta); es=[]
57 for _ in range(12):
58 es.append(np.max(np.abs(cur-zz))); cur=H(cur,delta,P,bb)
59 rr=es[-1]/es[-2] if es[-2] > 1e-14 else 0.0
60 beta_sweep.append({"beta":bb,"late_error_ratio":float(rr),"predicted_upper":bb})
61 return {"beta": beta, "operator_ratio": float(ratio),
62 "max_error_over_bound": float(max(e/b for e,b in zip(errors,bounds))),
63 "iteration_ratios": positive, "predicted_ratio_upper": beta,
64 "beta_sweep":beta_sweep}
65
66
67def tied_mdp():
68 # State 0 is the decision state. Each branch reaches a different
69 # continuation state; continuation states self-loop.
70 B, S, beta = 3, 4, .8
71 P = np.zeros((B,S,S))
72 destinations = [1,2,3]
73 for b in range(B):
74 P[b,0,destinations[b]] = 1.
75 for x in range(1,S): P[b,x,x] = 1.
76 delta = np.zeros((B,S))
77 future_cost = np.array([0., .20, .35, .27])
78 delta[:,1:] = future_cost[1:]
79 # Set root deficits so the three fixed-point scores tie, then add a
80 # positive common offset so signed perturbations remain valid deficits.
81 z,_ = solve(delta,P,beta)
82 base_scores = z + delta
83 target = np.max(base_scores[:,0])
84 delta[:,0] += (target-base_scores[:,0])/2 + 1.0
85 return P, delta, beta
86
87
88def scaling_check():
89 P, delta0, beta = tied_mdp()
90 B = 3
91 # Verify exact tie before perturbation.
92 z0,_ = solve(delta0,P,beta); s0=z0+delta0
93 tie_spread=float(np.ptp(s0[:,0]))
94 Ns=[20,50,100,200]
95 widths=[]; curves={}
96 # kappa=d=1: predicted critical temperature tau ~ N^-1.
97 for N in Ns:
98 tau=.35/N
99 ts=np.linspace(-8*tau,8*tau,501)
100 ps=[]
101 for t in ts:
102 d=delta0.copy(); d[0,0]+=t
103 z,_=solve(d,P,beta)
104 scores=z[:,0]+d[:,0]
105 a=np.exp(-(scores-scores.min())/tau); p=a/a.sum()
106 ps.append(p[0])
107 ps=np.asarray(ps)
108 widths.append(float(quantile_width(ts,ps)))
109 curves[N]=(ts/tau,ps)
110 # Collapse error against N=100 on common normalized coordinates.
111 grid=np.linspace(-8,8,1001)
112 ref=np.interp(grid,*curves[100])
113 collapse=max(float(np.max(np.abs(np.interp(grid,*curves[N])-ref))) for N in Ns)
114 scaled_widths=[w*N for w,N in zip(widths,Ns)]
115
116 # Ordinary noisy argmin: fixed critic noise does not inherit N^-1 scaling.
117 noise_rng=np.random.default_rng(SEED+1)
118 sigma=.08; widths_arg=[]
119 ts=np.linspace(-.8,.8,301)
120 noises=noise_rng.normal(0,sigma,size=(30000,3))
121 for N in Ns:
122 # N is only a label for the sampled-pool regime here; fixed noise
123 # demonstrates the baseline's lack of critical N scaling.
124 probs=[]
125 for t in ts:
126 scores=np.zeros((len(noises),3)); scores[:,0]=2*t
127 probs.append(np.mean(np.argmin(scores+noises,axis=1)==0))
128 widths_arg.append(float(quantile_width(ts,np.asarray(probs))))
129 # Under critic noise, compare soft resolver probabilities and hard argmin
130 # against the known lowest continuation-cost branch (branch 0 at root).
131 rr=np.random.default_rng(SEED+9)
132 M=50000; noise=rr.normal(0,.06,size=(M,3)); perturb=.03
133 noisy=np.zeros((M,3))+noise; noisy[:,0]+=perturb
134 hard=np.mean(np.argmin(noisy,axis=1)==0)
135 tau=.03
136 logits=-noisy/tau; logits-=logits.max(axis=1,keepdims=True)
137 soft=np.exp(logits); soft/=soft.sum(axis=1,keepdims=True)
138 resolver_prob=float(np.mean(soft[:,0]))
139 return {"tie_spread":tie_spread, "Ns":Ns, "tau_rule":"0.35/N (kappa=d=1)",
140 "resolver_widths":widths, "resolver_width_times_N":scaled_widths,
141 "predicted_width_times_N_constant":True,
142 "normalized_curve_max_deviation":collapse,
143 "baseline_noisy_argmin_widths":widths_arg,
144 "baseline_width_times_N":[w*N for w,N in zip(widths_arg,Ns)],
145 "predicted_baseline_width_N_scaling":"constant width, hence width*N grows linearly",
146 "direct_noisy_selection": {"perturbation":perturb,"noise_sigma":.06,
147 "baseline_hard_best_branch_frequency":float(hard),
148 "resolver_soft_probability_best_branch":resolver_prob}}
149
150
151def main():
152 out={"seed":SEED, "contraction":contraction_check(), "scaling":scaling_check()}
153 Path("results.json").write_text(json.dumps(out,indent=2))
154 print(json.dumps(out,indent=2))
155
156if __name__ == "__main__": main()