1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
|
"""Deterministic numerical checks for the claims in THEORY.md."""
import itertools
import math
import numpy as np
def simultaneous_mse(rng, depth, width, directions, trials=12000):
gradients = rng.normal(size=(depth, width))
gradients /= np.linalg.norm(gradients, axis=1, keepdims=True)
estimate = np.zeros((trials, width), dtype=np.float64)
for _ in range(directions):
xi = rng.integers(0, 2, size=(trials, depth, width), dtype=np.int8)
xi = 2.0 * xi - 1.0
directional = np.einsum("tdw,dw->t", xi, gradients)
estimate += xi[:, 0, :] * directional[:, None] / directions
empirical = np.square(estimate - gradients[0]).sum(axis=1).mean()
total_energy = np.square(gradients).sum()
theoretical = (width * total_energy - np.square(gradients[0]).sum()) / directions
return empirical, theoretical
def check_simultaneous_variance():
print("SIMULTANEOUS RADEMACHER MSE")
print("depth width K empirical theory ratio")
cases = ((2, 8, 1), (4, 8, 1), (8, 8, 1),
(4, 32, 1), (4, 32, 4), (4, 32, 16))
for index, (depth, width, directions) in enumerate(cases):
empirical, theoretical = simultaneous_mse(
np.random.default_rng(100 + index), depth, width, directions)
ratio = empirical / theoretical
print(f"{depth:5d} {width:5d} {directions:2d} "
f"{empirical:9.4f} {theoretical:9.4f} {ratio:6.3f}")
assert abs(ratio - 1.0) < 0.06
def cubic_loss(z, linear, cubic):
return linear @ z + cubic * np.power(z, 3).sum() / 6.0
def check_sigma_bias():
# Exhaustive directions remove Monte Carlo error, exposing only the
# centered-difference bias. For this separable cubic it is exactly O(sigma^2).
dim = 8
z = np.linspace(-0.4, 0.5, dim)
linear = np.linspace(0.2, 0.9, dim)
cubic = 1.7
true_gradient = linear + 0.5 * cubic * np.square(z)
directions = np.asarray(list(itertools.product((-1.0, 1.0), repeat=dim)))
sigmas = np.asarray((0.005, 0.01, 0.02, 0.04))
biases = []
print("\nCENTERED FINITE-DIFFERENCE BIAS")
print("sigma bias_norm")
for sigma in sigmas:
estimates = []
for xi in directions:
coefficient = (cubic_loss(z + sigma * xi, linear, cubic)
- cubic_loss(z - sigma * xi, linear, cubic)) / (2.0 * sigma)
estimates.append(xi * coefficient)
bias = np.linalg.norm(np.mean(estimates, axis=0) - true_gradient)
biases.append(bias)
print(f"{sigma:5.3f} {bias:.9e}")
slope = np.polyfit(np.log(sigmas), np.log(biases), 1)[0]
print(f"log-log slope={slope:.6f}")
assert abs(slope - 2.0) < 1e-6
def check_descent_threshold():
rng = np.random.default_rng(301)
dim = 20
gradient = rng.normal(size=dim)
orthogonal = rng.normal(size=dim)
orthogonal -= orthogonal.dot(gradient) * gradient / gradient.dot(gradient)
orthogonal /= np.linalg.norm(orthogonal)
neg_gradient = -gradient / np.linalg.norm(gradient)
cosine = 0.35
update = cosine * neg_gradient + math.sqrt(1.0 - cosine ** 2) * orthogonal
beta = 2.3
threshold = 2.0 * cosine * np.linalg.norm(gradient) / (beta * np.linalg.norm(update))
def quadratic_change(step):
return step * gradient.dot(update) + 0.5 * beta * step ** 2 * update.dot(update)
below = quadratic_change(0.9 * threshold)
above = quadratic_change(1.1 * threshold)
print("\nSMOOTH DESCENT THRESHOLD")
print(f"cos={cosine:.3f} eta_max={threshold:.6f} "
f"delta_below={below:+.6f} delta_above={above:+.6f}")
assert below < 0 < above
def check_predictor_timescale():
rng = np.random.default_rng(401)
examples = 50000
cells = 12
h = rng.normal(size=(examples, cells))
slopes = np.exp(rng.normal(scale=0.25, size=cells))
traffic = h * slopes
print("\nNEUTRAL PREDICTOR TIMESCALE")
print("eta steps residual_power empirical theory")
for eta in (0.002, 0.01, 0.05):
for steps in (0, 20, 100, 500):
predictor = slopes * (1.0 - (1.0 - eta) ** steps)
residual = traffic - h * predictor
empirical = np.square(residual).mean() / np.square(traffic).mean()
theoretical = (1.0 - eta) ** (2 * steps)
print(f"{eta:5.3f} {steps:5d} {empirical:14.8f} {theoretical:14.8f}")
assert abs(empirical - theoretical) < 1e-12
def check_conditional_projection():
"""Finite-sample L2 projection realizes the innovation identities exactly."""
rng = np.random.default_rng(351)
examples = 4096
h = rng.uniform(-1.5, 1.5, size=(examples, 2))
basis = np.column_stack((
np.ones(examples), h[:, 0], h[:, 1], np.square(h[:, 0]),
np.square(h[:, 1]), h[:, 0] * h[:, 1], np.sin(h[:, 0]),
np.cos(h[:, 1])))
coefficients = rng.normal(size=(basis.shape[1], 3))
conditional_mean = basis @ coefficients
innovation = rng.normal(scale=0.4, size=conditional_mean.shape)
# Project the finite-sample noise off every soma-measurable basis vector.
innovation -= basis @ np.linalg.lstsq(basis, innovation, rcond=None)[0]
nuisance = conditional_mean + innovation
linear_basis = basis[:, :3]
restricted = linear_basis @ np.linalg.lstsq(
linear_basis, nuisance, rcond=None)[0]
lhs = np.square(nuisance - restricted).mean()
irreducible = np.square(innovation).mean()
approximation = np.square(conditional_mean - restricted).mean()
orthogonality = np.abs(basis.T @ innovation / examples).max()
teaching = rng.normal(size=conditional_mean.shape)
raw = teaching + nuisance
residual = teaching + innovation
alpha = (np.linalg.norm(residual, axis=1)
/ np.linalg.norm(raw, axis=1).clip(min=1e-12))
matched = alpha[:, None] * raw
def row_cosine(left, right):
numerator = np.sum(left * right, axis=1)
denominator = (np.linalg.norm(left, axis=1)
* np.linalg.norm(right, axis=1)).clip(min=1e-12)
return numerator / denominator
direction_difference = np.abs(
row_cosine(raw, teaching) - row_cosine(matched, teaching)).max()
norm_difference = np.abs(
np.linalg.norm(matched, axis=1) - np.linalg.norm(residual, axis=1)).max()
print("\nCONDITIONAL INNOVATION PROJECTION")
print(f"orthogonality={orthogonality:.3e} pythagorean_error="
f"{abs(lhs - irreducible - approximation):.3e}")
print(f"norm_match_error={norm_difference:.3e} "
f"direction_change={direction_difference:.3e}")
assert orthogonality < 2e-14
assert abs(lhs - irreducible - approximation) < 2e-14
assert norm_difference < 2e-14
assert direction_difference < 2e-14
def squared_cosine(x, y):
return float((x.ravel() @ y.ravel()) ** 2
/ ((x.ravel() @ x.ravel()) * (y.ravel() @ y.ravel())))
def check_innovation_identification():
rng = np.random.default_rng(501)
examples = 100000
cells = 8
h = rng.normal(size=(examples, cells))
traffic_slopes = np.linspace(0.8, 1.5, cells)
teaching_slopes = np.linspace(0.3, 0.7, cells)
innovation = rng.normal(scale=0.6, size=(examples, cells))
noise = rng.normal(scale=0.15, size=(examples, cells))
teaching = h * teaching_slopes + innovation
traffic = h * traffic_slopes
apical = teaching + traffic + noise
neutral_residual = apical - traffic
task_coeff = (apical * h).mean(axis=0) / np.square(h).mean(axis=0)
task_residual = apical - h * task_coeff
raw_alignment = squared_cosine(apical, teaching)
neutral_alignment = squared_cosine(neutral_residual, teaching)
task_alignment = squared_cosine(task_residual, teaching)
predictable_fraction = np.square(h * teaching_slopes).mean() / np.square(teaching).mean()
retained_fraction = np.square(teaching - h * teaching_slopes).mean() / np.square(teaching).mean()
print("\nINNOVATION IDENTIFICATION")
print(f"squared cosine raw={raw_alignment:.4f} neutral={neutral_alignment:.4f} "
f"task_fit={task_alignment:.4f}")
print(f"teaching predictable={predictable_fraction:.4f} "
f"retained_after_task_fit={retained_fraction:.4f}")
assert neutral_alignment > raw_alignment + 0.15
assert task_alignment < neutral_alignment - 0.15
assert abs(retained_fraction - (1.0 - predictable_fraction)) < 0.01
def main():
check_simultaneous_variance()
check_sigma_bias()
check_descent_threshold()
check_conditional_projection()
check_predictor_timescale()
check_innovation_identification()
print("\nALL THEORY CHECKS PASSED")
if __name__ == "__main__":
main()
|
