Classical statistics has a well-known rule of thumb: a model with more parameters than training examples will OVERFIT catastrophically, perfectly memorizing the training set (including its noise) while generalizing terribly to new data. This intuition is genuinely correct for classical models, and it predicts that test error should get WORSE and worse as a model's parameter count approaches, and then exceeds, the number of training examples. Modern deep learning routinely violates this rule in a specific, reproducible, and genuinely surprising way: massively overparameterized networks (far more parameters than training examples, sometimes by many orders of magnitude) often generalize EXTREMELY well, frequently better than a "properly-sized," classically-recommended model would.
The "double descent" phenomenon (Belkin et al., 2019, building on earlier statistical learning theory) reconciles these two facts by showing that test error, plotted against parameter count, doesn't behave the way classical theory predicts (monotonically increasing once you cross the point where parameters equal training samples). Instead, it follows a genuine DOUBLE-DESCENT curve: test error decreases as parameters increase (as classical theory expects), then spikes sharply right around the "interpolation threshold" (where the model has JUST enough capacity to fit the training data exactly, and is forced to do so in a fragile, noise-sensitive way), and then, surprisingly, DECREASES again as parameters continue to grow well past that threshold, because a sufficiently overparameterized model has enough freedom to find a SMOOTH interpolating solution, not just any interpolating solution.
Implement build_features(x, weight, bias) (the same sigmoid random-feature construction from [01-universal-approximation], but taking an already-drawn weight/bias so the identical features can be applied to both a training set and a test set), fit_min_norm(hidden, y) (solve for output weights via np.linalg.pinv, the MINIMUM-NORM solution, essential once num_features exceeds the sample count and the system becomes underdetermined), and train_and_test_mse(x_train, y_train, x_test, y_test, num_features, rng), which draws one shared random feature set, fits it on the training data, and reports both training and test MSE.
build_features computes sigmoid(np.outer(x, weight) + bias), taking weight/bias as arguments rather than drawing them itself.fit_min_norm must use np.linalg.pinv, not np.linalg.lstsq: the pseudoinverse specifically returns the MINIMUM-NORM solution in the overparameterized (underdetermined) case, which is the numerically well-behaved choice this whole demonstration depends on.train_and_test_mse must draw weight/bias ONCE and reuse the exact same values for both build_features(x_train, ...) and build_features(x_test, ...), so training and test features live in the same space.0 once num_features >= len(x_train) (the interpolation threshold has been crossed).build_features is a one-liner: sigmoid(np.outer(x, weight) + bias), identical to [01-universal-approximation]'s random_hidden_features, just with weight/bias passed in rather than generated internally.
fit_min_norm(hidden, y) = np.linalg.pinv(hidden) @ y. Unlike np.linalg.lstsq, which can behave ambiguously or error on some underdetermined systems, np.linalg.pinv always returns a well-defined result, and specifically the minimum-Euclidean-norm one when multiple exact solutions exist.
In train_and_test_mse: weight = rng.randn(num_features); bias = rng.randn(num_features), drawn ONCE. Build hidden_train and hidden_test using those SAME weight/bias values (not two separate draws). Fit output_weight via fit_min_norm(hidden_train, y_train) (training data only), then compute both MSEs using hidden_train @ output_weight and hidden_test @ output_weight respectively.
Click "Run Tests" to test your implementation