import json, math, random import numpy as np from scipy.optimize import linprog import torch import torch.nn as nn SEED = 2056 np.random.seed(SEED); random.seed(SEED); torch.manual_seed(SEED) torch.set_num_threads(4) def fresh_seed(s=SEED): np.random.seed(s); random.seed(s); torch.manual_seed(s) class ICNN(nn.Module): def __init__(self, d=2, h=24, L=2): super().__init__(); self.L=L self.A=nn.ParameterList([nn.Parameter(torch.randn(h, d if k==0 else h)*.35) for k in range(L)]) self.U=nn.ParameterList([nn.Parameter(torch.randn(h,d)*.35) for _ in range(L)]) self.b=nn.ParameterList([nn.Parameter(torch.zeros(h)) for _ in range(L)]) self.aw=nn.Parameter(torch.zeros(h)); self.u=nn.Parameter(torch.zeros(d)); self.c=nn.Parameter(torch.zeros(1)) def forward(self,x): z=x for k in range(self.L): W=torch.nn.functional.softplus(self.A[k])+1e-6 z=torch.relu(z@W.T + x@self.U[k].T + self.b[k]) return z@torch.nn.functional.softplus(self.aw)+x@self.u+self.c class MLP(nn.Module): def __init__(self,d=2,h=24,L=2): super().__init__(); layers=[]; q=d for _ in range(L): layers += [nn.Linear(q,h),nn.ReLU()]; q=h layers += [nn.Linear(q,1)] self.net=nn.Sequential(*layers) def forward(self,x): return self.net(x).squeeze(-1) def convexity_sweep(): # Prediction 1: for nonnegative W and output w, midpoint Jensen violations are zero. # Prediction 2: introducing negative hidden-to-hidden entries causes violations at a # positive rate; this is the exact architectural boundary (sign >= 0). rng=np.random.default_rng(SEED); n=30000; d=2; h=4 x=rng.uniform(-2,2,(n,d)); y=rng.uniform(-2,2,(n,d)); t=rng.uniform(0,1,n)[:,None] xm=t*x+(1-t)*y rows=[] for negfrac in [0.0, .05, .10, .25, .50, 1.0]: W1=np.abs(rng.normal(size=(h,d))); W2=np.abs(rng.normal(size=(h,h))) if negfrac: mask=rng.random(W2.shape)1e-8)), 'mean_positive_violation':float(np.maximum(gap,0).mean())}) return rows def controlled_sign_sweep(): # Controlled prediction: f_gamma(x)=ReLU(1+gamma*ReLU(x)) is convex for # gamma>=0; for gamma<0 it has a downward slope after x=0 and must violate # Jensen convexity. The pair (-1,1) at t=1/2 gives violation # f(0)-(f(-1)+f(1))/2 = -gamma/2 for -1<=gamma<0. xs=np.linspace(-1,1,401) rows=[] for gamma in [-2.0,-1.0,-.5,-.1,0.0,.1,.5,1.0,2.0]: def f(x): return np.maximum(0,1+gamma*np.maximum(0,x)) maxv=0.0 for a in xs: for b in xs[::8]: for t in (.25,.5,.75): gap=float(f(t*a+(1-t)*b)-t*f(a)-(1-t)*f(b)) maxv=max(maxv,gap) predicted=max(0.0,-gamma/2.0) if -1<=gamma<0 else (0.5 if gamma < -1 else 0.0) rows.append({'gamma':gamma,'predicted_max_violation':predicted,'observed_max_violation':maxv}) return rows def epigraph_check(): # Prediction 3: LP epigraph minimum equals forward ReLU output when all downstream # coefficients are positive; gap should be numerical tolerance. rng=np.random.default_rng(SEED+1); d=2; h=3; N=80 W=np.abs(rng.normal(size=(h,d)))*.7; U=rng.normal(size=(h,d)); b=rng.normal(size=h); w=np.abs(rng.normal(size=h))+.2 gaps=[] for _ in range(N): x=rng.uniform(-2,2,d); pre=W@x+U@x+b; z=np.maximum(0,pre); forward=w@z # min w'z subject z >= pre and z >=0 res=linprog(w,A_ub=-np.eye(h),b_ub=-np.maximum(pre,0),bounds=[(0,None)]*h,method='highs') gaps.append(abs(res.fun-forward)) # A two-layer version, with fixed x and nonnegative W2, is also an LP. W1=np.abs(rng.normal(size=(h,d))); U1=rng.normal(size=(h,d)); b1=rng.normal(size=h) W2=np.abs(rng.normal(size=(h,h))); U2=rng.normal(size=(h,d)); b2=rng.normal(size=h); w2=np.abs(rng.normal(size=h))+.2 gaps2=[] for _ in range(N): x=rng.uniform(-2,2,d); p1=W1@x+U1@x+b1; z1=np.maximum(0,p1); p2=W2@z1+U2@x+b2; z2=np.maximum(0,p2); forward=w2@z2 # Variables z1,z2. Constraints -z1<=-p1; -z2+W2z1<=-(U2x+b2), plus nonnegative bounds. Aub=np.zeros((2*h,2*h)); bub=np.zeros(2*h) Aub[:h,:h]=-np.eye(h); bub[:h]=-p1 Aub[h:,:h]=W2; Aub[h:,h:]=-np.eye(h); bub[h:]=-(U2@x+b2) res=linprog(np.r_[np.zeros(h),w2],A_ub=Aub,b_ub=bub,bounds=[(0,None)]*(2*h),method='highs') gaps2.append(abs(res.fun-forward)) return {'one_layer_max_abs_gap':float(max(gaps)), 'two_layer_max_abs_gap':float(max(gaps2)), 'predicted_gap':0.0} def fit_compare(): # Convex target: quadratic. Compare ICNN with a width-matched ordinary ReLU MLP. fresh_seed(); rng=np.random.default_rng(SEED) X=rng.uniform(-2,2,(600,2)).astype('float32'); y=(.8*X[:,0]**2+1.4*X[:,1]**2+.35*X[:,0]*X[:,1]+.2*X[:,0]).astype('float32') tr=slice(0,450); va=slice(450,None); xt=torch.tensor(X); yt=torch.tensor(y) out={} for name, model in [('icnn',ICNN()),('relu_mlp',MLP())]: fresh_seed(SEED+ (0 if name=='icnn' else 1)); model=model.to('cpu'); opt=torch.optim.AdamW(model.parameters(),lr=.01,weight_decay=1e-5) for step in range(1200): opt.zero_grad(); pred=model(xt[tr]); loss=((pred-yt[tr])**2).mean(); loss.backward(); opt.step() with torch.no_grad(): train=float(((model(xt[tr])-yt[tr])**2).mean()); val=float(((model(xt[va])-yt[va])**2).mean()) out[name]={'train_mse':train,'validation_mse':val,'steps':1200} out['metric']='validation MSE (lower is better)'; return out def main(): result={'device':'cpu','predictions':{ 'P1':'nonnegative W,w => Jensen violation <= numerical tolerance (predicted 0)', 'P2':'negative W entries create positive Jensen violations; violation rate rises with negative fraction', 'P3':'positive output coefficients make ReLU epigraph LP exact (predicted objective gap 0)'}, 'convexity_sweep':convexity_sweep(),'controlled_sign_sweep':controlled_sign_sweep(),'epigraph_check':epigraph_check(),'fit_comparison':fit_compare()} with open('results.json','w') as f: json.dump(result,f,indent=2) print(json.dumps(result,indent=2)) main()