Universal approximation: why one wide hidden layer can fit any function, in principle

The problem, from first principles

A striking, genuinely surprising mathematical fact underlies everything else in this curriculum: a neural network with just ONE hidden layer, no depth at all, can approximate essentially any reasonable function to arbitrary precision, given enough hidden units. This is the Universal Approximation Theorem (Cybenko, 1989; Hornik, 1991), and it's worth sitting with why it's surprising: [02-layers/01-linear-forward]'s linear_forward on its own can only ever represent a STRAIGHT LINE (or a flat hyperplane, in higher dimensions), no matter how you set its weights. Stack a nonlinearity (like sigmoid, used here) between two linear layers, though, and each hidden unit becomes a smoothly-adjustable "bump" or "step" that can be placed, scaled, and combined with every other hidden unit's own bump, and a large enough COLLECTION of these adjustable bumps can be combined to approximate arbitrarily complex curves.

The theorem is an EXISTENCE proof, not a practical training recipe: it says a wide-enough network with the right weights exists, not that gradient descent will easily find those weights, or that "wide and shallow" is actually a good architecture in practice (it usually isn't; Representation learning: why depth learns hierarchical features, not one big lookup, immediately following this question, is exactly about why DEPTH, not just width, tends to work much better in practice). This question demonstrates the existence claim directly and numerically: build a wide hidden layer with RANDOM (untrained) weights, and show that simply choosing the right OUTPUT layer weights (a cheap, closed-form least-squares fit, no gradient descent needed) is enough to approximate a genuinely nonlinear target function, and that approximation gets dramatically better as the hidden layer gets WIDER.

From theory to code

Implement random_hidden_features(x, num_hidden, rng), fit_output_weights(hidden, y), and approximate_function(x, y, num_hidden, rng). random_hidden_features builds num_hidden sigmoid neurons, each with its own RANDOM weight and bias (already generated for you, weight = rng.randn(num_hidden) * 5.0 and similarly for bias), applied to every point in x, producing a (len(x), num_hidden) matrix. fit_output_weights solves, via np.linalg.lstsq, for the linear combination of those hidden-layer columns that best matches target y. approximate_function ties both together and reports the mean squared error of the fit.

Constraints

• random_hidden_features must return a matrix of shape (len(x), num_hidden): one column per hidden neuron, applying THAT neuron's own weight and bias to every point in x.
• fit_output_weights must use np.linalg.lstsq(hidden, y, rcond=None) (least squares), not gradient descent, this is a closed-form, exact best fit given the FIXED (untrained) hidden layer.
• approximate_function must reuse both functions above rather than reimplementing their logic.
• The returned mse should decrease substantially as num_hidden increases, for the same target function, demonstrating the "wider approximates better" intuition.

Hints