I know the ground truth matrix A1 is 20x10, but my steal.py should only rely on black-box forward() queries, saving the true weights purely for verification. My plan is to exploit the network's piecewise-linear structure: each hidden neuron defines a hyperplane where its activation flips sign, and by locating these critical points and comparing gradients on either side of the boundary, I can recover the corresponding row direction of A1.
Concretely, I'll sample random lines through input space, scan along them to detect kinks via nonzero second differences, then bisect to pinpoint each breakpoint precisely. At each kink, computing gradients via finite differences on both sides and taking the difference should yield a vector proportional to that neuron's A1 row — I'll keep sampling new lines and dedupe discovered directions until no new ones appear, since I don't know the neuron count in advance. I need to be careful about scale (~0.1 for weights and biases) and epsilon choices for the finite-difference gradient estimates near each localized kink.
I should localize kinks by scanning for nonzero curvature in triples of points, then narrowing down with bisection: repeatedly test curvature on sub-intervals to determine which half contains the kink, iterating enough times to pin down its location precisely.
I don't need extreme precision—just enough that delta-offset points land in distinct linear regions on either side of the kink, so localizing to ~1e-9 and offsetting by 1e-5 should work as long as neighboring kinks aren't too close. One risk is a neuron with a tiny second-layer coefficient making its kink invisible, but with coefficients around 0.1 that shouldn't be an issue; another risk is a hyperplane nearly parallel to the sampling line causing the crossing to occur far away, which I can mitigate by sampling many random lines with different starting points.
Given the scale of the coefficients and inputs, crossings should occur within a modest range, so I'll scan t across a wide interval with fine steps, expecting roughly a few thousand queries per line via bisection and gradient checks per kink—this cost is manageable since the forward pass is cheap. For deduplication, I'm working out how to normalize each row's direction to a canonical sign convention to avoid counting the same hyperplane twice.
I'll structure the script to keep collecting rows until no new ones appear for several consecutive scan lines, then cap the total and return whatever's found since the true neuron count is unknown. Verification will match discovered rows to the true A1 up to scale and permutation, accounting for potential dead or always-active neurons that a scan might miss—though random initialization should keep that unlikely. Now I'm writing the actual script, making sure forward.py exposes a forward function importable by steal.py from the /app directory, and double-checking how sign ambiguity interacts with the permutation/scaling equivalence.
I still can't pin down the sign of A1[i] just from the kink location alone—I need to know which side of the hyperplane actually activates the neuron. Maybe I can test this directly by evaluating the function far from the hyperplane on each side, where a single neuron's contribution should dominate or be isolable.
So the sign of w seems fundamentally ambiguous from local kink data alone. I try flipping w_i, b_i, and A2[i] simultaneously—since relu(-z) = relu(z) - z, this transformation introduces an extra linear term A2[i]*(w·x+b) that isn't absorbable elsewhere in this architecture, meaning the function actually changes unless something compensates for it.
I'm now wondering if the task's "up to permutation and scaling" language is meant to include sign flips as part of scaling ambiguity, but I want to check whether sign can be resolved rigorously via a global fit—though that introduces nonlinear fitting rather than a clean local rule. I'm searching for a local trick to pin down the sign using kink points instead.
This isn't fully decomposable, so I recall that model-stealing literature (Carlini et al.) recovers weight rows up to sign and scale locally, then resolves signs separately, which is nontrivial. Since the task only requires recovery up to permutation and scaling, I could just deliver rows up to scale including sign, but I want to attempt a proper sign-resolution step using their critical-point approach for thoroughness.
I'm working through the sign recovery test: pick a critical point on neuron i's hyperplane, perturb in a random direction, and examine how the function behaves as the perturbation crosses the hyperplane, trying to use this crossing structure to pin down the sign of the recovered weight vector.
Since I can query f exactly, once I recover each direction u_i and hyperplane offset β_i (from critical points where u_i·x+β_i=0), the model becomes a sum of ReLUs with unknown signs s_i and coefficients a_i, plus a constant. For fixed signs, fitting a_i and the constant is just least squares, but trying all 2^20 sign combinations is infeasible—though I notice flipping a sign only changes the basis function by a linear term, since relu(-z) = relu(z) - z.
So if I pick the wrong sign for some subset, the true function differs from my model by a leftover linear term that generally can't be absorbed into the ReLU basis, meaning the least-squares residual will be nonzero exactly when signs are wrong—giving me a way to detect and correct sign errors one at a time.
I need to work out the true sign relationship: z_i using an arbitrary orientation u_i means the true preactivation could be +|c_i|z_i or -|c_i|z_i, and I'm expanding relu(-z_i) = relu(z_i) - z_i to see how a sign flip on a subset F shows up as a nonzero linear term in the fitted decomposition — so the linear coefficient λ picks up exactly the contribution from flipped-sign units, letting me recover which signs were wrong.
So while the naive global fit mixes contributions, the decomposition should still be unique because different sign patterns produce genuinely distinct function families, differing by an affine offset. That means I can check candidate sign patterns by testing whether the residual after an exact fit vanishes, and use the all-+1 fit's linear term to search for a subset of flipped neurons whose combined effect cancels it out.
I'm settling on a more tractable strategy: start with all signs +1, compute residual, then do coordinate descent by flipping one neuron's sign at a time and checking whether the residual drops toward zero, since exact arithmetic guarantees zero residual only at the correct pattern.
I'm also reconsidering the problem's framing—"up to permutation and scaling" likely means nonzero real scaling is acceptable, so I could deliver rows recovered up to sign and scale as the main answer, with sign resolution as an optional refinement. This reminds me of a known trick for depth-2 networks (from Carlini et al.) using large steps along a direction from a critical point to probe neuron sign behavior directly.
That still gives the same sign ambiguity since we don't know which side is active without extra info. I'm realizing the distinguishing signal has to be global rather than a purely local probe near the hyperplane—something about how the function behaves as we move far in a fixed direction, or the actual magnitude of the function values, rather than just gradient jumps.
So I'm decomposing the relu using the identity relu(σz) = (σz + |z|)/2, which splits f(x) into a sum of absolute-value terms weighted by a_i/2 plus a purely linear component built from the a_iσ_i/2 coefficients. This reformulation might let me separate the piecewise-linear "kink" structure from the sign-dependent linear part, giving a cleaner path to recovering the σ_i signs exactly.
But then I realize fitting Σc_i|u_i·x+β_i| + λx + μ by least squares gives me c_i, λ, μ uniquely (31 unknowns, well-posed), and I can recover a_i = 2c_i directly. The catch is that λ = Σc_iσ_iu_i only gives 10 equations for 20 unknown signs σ_i, so the system is underdetermined — the signs still aren't uniquely pinned down this way.
Since 10 equations over the discrete set of 2^20 sign patterns is a measure-zero coincidence, generically σ is unique, meaning signs are determined but finding them is essentially a subset-sum search — the actual hard combinatorial step people solve with smarter methods. With 2^20 only around a million, this is brute-forceable directly, so I'm thinking about vectorizing the search as a matrix of sign patterns dotted against the 10-dimensional vector to check which one matches λ.
I should also consider meet-in-the-middle: split the neurons into two halves, enumerate all sums for each half separately (1024 each for n=20), then match complementary sums against λ. I'll cap brute-force/meet-in-the-middle at n≤24 and fall back to reporting only up-to-sign for larger n.
The magnitude of each row is genuinely unidentifiable since w_i and a_i can be rescaled inversely, so scale stays free while sign can be resolved. My plan is to output rows correct up to permutation and scale, with sign fixed via the matching approach.
My concrete script: search critical points along random lines to recover directions and offsets, then fit the piecewise-linear function via least squares on random samples to validate residuals near zero, resolve signs via meet-in-the-middle when feasible, and save the unit-normalized, sign-resolved rows to the output file.
I'll also add a verification pass comparing against the true weight matrix, accepting the row-scaling ambiguity by normalizing everything to unit norm. Now I'm working out the numerical details for detecting kinks along each sampled line—using curvature of nearby function samples to flag a slope discontinuity, with tolerances chosen relative to expected roundoff and slope jump magnitudes.
Once a bracket containing exactly one kink is found, I realize I can compute the exact kink location analytically by fitting the two linear pieces on either side and solving for their intersection, rather than relying purely on bisection—this gives a much more precise estimate that I can then double check with the curvature test at the computed point.
However, coarse grids risk having multiple kinks per bracket, so I want a consistency check: verify the linear extrapolations from both sides actually agree at t*, and if not, subdivide further. I'll settle on running roughly 50 bisection iterations with the curvature test to reach machine precision, reusing endpoint evaluations to save cost, then use a small delta offset to measure the gradient jump at t*.
For the delta and h choice, I'll pick values around 1e-4 for region separation and 1e-5 for the central-difference step, checking that roundoff error stays small relative to typical gradient magnitudes so the resulting jump vector estimate stays reliable.
For deduplication, I'll normalize each jump vector and fix its sign so the largest-magnitude component is positive, then compare normalized vectors by cosine similarity to merge duplicates. Each hyperplane can then be written using this canonical unit normal and an offset determined by the critical point on it.
As a self-check, I'm fitting a linear model against sampled points using absolute-value hinge terms plus the raw coordinates and a constant, then checking the residual to confirm every neuron was recovered — a large residual would flag a missing one. Now I'm carefully re-deriving the sign recovery relation, expressing each pre-activation as a signed multiple of its canonical form so I can solve for how the coefficients and signs combine to reconstruct the original weight vector.
Fitting the model gives me c_i and λ, so now I need to solve Σ σ_i c_i u_i = λ for σ ∈ {±1}^n, with μ unconstrained since b2 is unknown, which matches what I found earlier. I'm planning a meet-in-the-middle approach: define v_i = c_i u_i, split the indices into two halves, enumerate all sign combinations for each half to get partial sums, then search for the pair of partial sums whose total is closest to λ using efficient pairwise distance computation.
I want the best match's residual to be near-zero relative to the scale of λ, and the second-best match to be clearly worse to confirm uniqueness. The output rows are then σ_i u_i, which correspond exactly to normalized w_i. For large n where brute-force sign recovery becomes infeasible (say n > 26), I'll fall back to outputting unit vectors with a canonical sign convention and log a warning. For the discovery process, I'll keep gathering lines until several consecutive additions yield no new directions, but I'll rely mainly on the fit residual as the real completeness check -- if the residual stays too high after fitting, I'll keep searching for more lines rather than stopping prematurely.
Since importing `forward.py` resets numpy's global random seed via its own `np.random.seed(0)` call, I'll avoid relying on the global RNG state and instead create my own `np.random.Generator` via `default_rng`, importing forward first to get that seed reset out of the way. I'm now writing the actual `/app/steal.py` file, including a query counter to track and report usage, plus a kink-finding routine that scans along a direction for points where the function's curvature changes sign, then isolates each kink location via bisection refinement.
Working through the bisection logic: I need to check curvature over the left half-interval first, shrinking to that half if a kink is detected there, otherwise moving to the right half — since curvature magnitude shrinks proportionally to interval width, I have to be careful about when floating-point roundoff overwhelms the signal, so I'm capping the bisection to stop once the interval width drops below roughly 1e-9 rather than iterating all the way to machine precision.
I'm working out sensible step sizes for finite differencing near a kink—balancing delta and h so roundoff error stays small relative to the gradient jump, since hyperplane spacing along a random line in this dimension is typically 0.1-1 but can occasionally be much tighter. Settling on delta and h around 1e-4 seems like a reasonable compromise to keep the gradient estimate accurate enough for reliable jump detection.
Rather than adding an expensive per-crossing linearity check to validate each jump, I think it's simpler to trust the global fit residual as the main validator, use a moderate deduplication tolerance for cosine similarity to catch near-duplicate directions, and afterward refine each recovered direction by locating nearby critical points and doing a local least-squares re-estimate.
I'm settling on delta = h = 2e-5 for the central difference window, checking that roundoff error stays negligible (~5e-12) and that no other hyperplane falls within the resulting window near a critical point. I also need to watch for the edge case where the perturbation itself along a coordinate direction accidentally crosses hyperplane i, since the distance from q to that hyperplane depends on how aligned the direction d is with the hyperplane's normal.
Since I'm sampling many random lines, each neuron should get found repeatedly, so I can cross-check by requiring agreement across independent discoveries rather than trusting any single fragile estimate. I'll reject a kink if |d·u| is too small and rely on other lines to catch it, then implement an adaptive retry with an initial delta and gradient step size.
Other planes shouldn't intrude within the tolerance range, so this approach holds up fine. For the grid scan I'll sweep t across [-20,20] in steps of 0.02, giving ~2000 points per line across up to 40 lines, which is cheap since each forward call is microseconds—totaling roughly 1-2 seconds runtime. I'll sample the base point from a scaled normal distribution and use a random unit direction, then write the script and verify results afterward.
For the least-squares fit I'll sample points from a scaled normal distribution, build a design matrix from the absolute values of linear combinations plus a bias column, and solve via lstsq, checking the relative residual to decide if more discovery rounds are needed. For sign recovery I'm setting up a meet-in-the-middle approach using numpy to enumerate all sign combinations and their weighted sums, splitting the enumeration into two halves for efficiency.
I'm reasoning through the size limits: with n=20 each half gives 1024 sign patterns, making a 1024×1024 distance computation feasible, but for larger n like 28 the squared blowup (2.7e8 distances) becomes too costly, so I'm capping the approach at n≤24 to keep the computation tractable within chunked processing.
I also need to consider edge cases where a neuron is never active or always active in the sampled region — such a neuron would behave linearly there, meaning its contribution gets absorbed into the linear term and its residual would falsely appear at zero, so I want to verify the sampling range is wide enough to catch these boundary cases before finalizing the fit.