How Neural Networks Learn

A model that got worse before it got better

Picture a fraud-detection team at a payments company operating at UPI scale, scoring tens of millions of transactions a day. Every quarter they retrain their neural scorer on fresh data, and every quarter they widen the hidden layers a little, because more transaction volume means more patterns to represent. One quarter, at a particular width, the held-out fraud-catch rate collapses: worse than the smaller model they shipped last quarter, worse than random guessing on some slices. The instinct is to assume a bug, a leaked feature, a broken pipeline. But suppose the pipeline is clean. Suppose the model at that exact width is behaving exactly as the mathematics says it should, and the next width up, wider still, fixes the collapse and beats every earlier model. That non-monotonic shape, worse-then-better as capacity grows, is not an anomaly to be engineered away. It has a name, it has a precise mechanism, and understanding it is the difference between a practitioner who panics at a bad retrain and a researcher who knows exactly where on the curve they are standing.

This chapter is not about the mechanics of how a network's weights get updated; the chain rule and the backward pass are covered in the companion chapter Backpropagation: How Neural Networks Learn. This chapter asks a different, research-literacy question: once you have a network large enough to represent almost anything, and an optimizer that can drive its training error to exactly zero, why does it generalize at all, and why does the relationship between "how big" and "how good" refuse to behave the way a first course in statistics says it must?

What "learning" means, formally

Strip away the neural-network framing for a moment. Learning, in the statistical sense, is choosing a function from a hypothesis class H that minimizes expected loss on unseen data, using only a finite training sample. You cannot minimize the true expected loss directly, because you do not have access to the true data distribution, only n samples from it. So you minimize a proxy, the empirical risk (average loss on the training set), and hope the result generalizes. Two separate questions hide inside that hope. First, is the function you want even representable inside H? That is a question about approximation, about the expressive power of your hypothesis class. Second, given that it is representable, will minimizing empirical risk actually find a function that performs well on new data, rather than one that has simply memorized the training sample? That is a question about generalization. Neural network theory keeps these two questions distinct because, as the next two sections show, a positive answer to the first tells you almost nothing about the second.

The Universal Approximation Theorem, and its silences

The oldest positive result about neural networks answers the approximation question. Cybenko (1989) proved that a single-hidden-layer network with a sigmoid activation can approximate any continuous function on a compact subset of R^n to arbitrary accuracy, given enough hidden units. Hornik, Stinchcombe, and White (1989) generalized this to a broad class of "squashing" activation functions, and Hornik (1991) extended it further, showing that it is the multilayer feedforward architecture itself, not any specific activation function, that supplies universality. Leshno et al. (1993) then pinned down the exact necessary and sufficient condition: a single-hidden-layer network is a universal approximator if and only if its activation function is not a polynomial. This is why ReLU, tanh, sigmoid, GELU, and dozens of other activations all work for building universal approximators, while a network built entirely from linear (or polynomial) activations never will, no matter how deep or wide, because a composition of linear maps is still just a linear map.

A separate strand asks whether depth can substitute for width. Lu et al. (NeurIPS 2017) showed that ReLU networks with width bounded by a small constant multiple of the input dimension are still universal approximators, provided you are allowed unbounded depth. Width and depth trade off against each other as resources for expressive power; neither one alone is the "real" source of universality.

Read carelessly, the Universal Approximation Theorem sounds like a complete license: if the function you want exists, a big enough network can represent it, so training should find it, and it should work. Every clause after "represent it" is where the theorem stays silent, and each silence is worth naming precisely, because this is the first misconception this topic invites.

• It is an existence proof, not a construction. The theorem guarantees weights exist that achieve the approximation; it gives you no algorithm for finding them and no guarantee that gradient descent, which only has access to local gradient information on a non-convex loss surface, will ever reach them.
• It does not bound the width you need. For a generic continuous function on a d-dimensional domain, the number of hidden units required to hit a fixed accuracy can grow exponentially in d, the same curse of dimensionality that makes naive nearest-neighbor methods intractable in high dimensions. "A network exists" and "a network of a size you can afford to train exists" are different claims.
• It says nothing about data outside the domain you approximated on. The theorem is about matching a target function on a compact set; it makes no statement about a finite training sample drawn from that set, nor about behavior on points the training sample never covered. Approximation theory and statistical generalization are formally unrelated results that happen to both live inside the phrase "neural networks can learn anything."

So the Universal Approximation Theorem tells you neural networks are expressive enough to be dangerous, not that any particular training run will end up somewhere useful. The rest of this chapter is about what actually governs that second question.

What classical theory predicts, and gets visibly wrong

Classical statistical learning theory, built on tools like Vapnik-Chervonenkis (VC) dimension, gives a clean prediction about the relationship between model complexity and test error. As you increase model capacity, training error falls (a more expressive hypothesis class can fit the training sample more closely), while test error traces a U-shape: it falls at first as the model escapes underfitting, reaches a minimum at some "just right" capacity, then rises again as the model starts fitting noise specific to the training sample rather than signal that generalizes, the classical bias-variance trade-off. Crucially, this theory predicts that once model capacity is large enough to drive training error to exactly zero (the model can perfectly interpolate every training point, including its noise), test error should be at its worst, and should stay bad or get worse as you add still more capacity beyond that point. A model with far more parameters than training examples, by this account, is a model that has simply memorized noise and has no business generalizing.

Modern deep networks routinely violate this prediction. Production language models and vision networks have parameter counts orders of magnitude larger than their training sets in the relevant statistical sense, achieve zero or near-zero training error, and still generalize well, often better than smaller models tuned to sit at the classical "sweet spot." Belkin, Hsu, Ma, and Mandal (PNAS, 2019) were the first to name and formalize this pattern across model families, from decision trees to random-feature models to small neural networks, calling it "double descent." Nakkiran et al. (ICLR, 2020) then showed the same shape holds for real, large deep networks (ResNets, CNNs, and transformer-style models), as a function of model size, of training time, and of dataset size, and introduced the language of an "interpolation threshold": the model complexity at which training error first reaches (near) zero.

The double descent curve

The shape itself is what makes the phenomenon worth a name rather than a footnote. Below the interpolation threshold, test risk behaves exactly as classical theory predicts, a U-shaped curve with a genuine "classical sweet spot." Right at the threshold, where the model has just enough capacity to fit the training data exactly, test risk spikes, often dramatically, because the interpolating solution at that exact capacity is poorly conditioned: there is exactly one way to fit the data exactly, and that one way is extremely sensitive to the specific noise in the training sample. Past the threshold, in the heavily overparameterized regime where many different sets of weights could all achieve zero training error, test risk falls again, a second descent, because the optimizer no longer has to pick the one fragile interpolating solution; it has room to pick among many, and (under conditions this chapter returns to) it tends to pick a well-behaved one.

A worked example: watching the spike happen

The cleanest place to see double descent is not a deep network at all, it is ordinary linear regression, because there the entire curve can be computed in closed form and reproduced exactly. Fix n = 20 training points. The true relationship uses only 5 informative features, with true coefficients beta_true = [3.0, -2.0, 1.5, 0.5, -1.0], plus Gaussian label noise with standard deviation sigma = 0.5. Now sweep the number of features the model is allowed to use, p, from 1 up to 3000, where features beyond the first 5 are pure noise, i.i.d. Gaussian, with no relationship to the label at all. For every value of p, fit ordinary least squares; when p > n, the system is underdetermined and has infinitely many exact solutions, so we take the minimum-norm one, which is exactly what numpy.linalg.lstsq returns via its SVD-based solver.

import numpy as np

np.random.seed(0)

n_train = 20
n_test = 200
p_max = 3000
k_true = 5
beta_true = np.array([3.0, -2.0, 1.5, 0.5, -1.0])
sigma = 0.5

X_train_full = np.random.randn(n_train, p_max)
X_test_full = np.random.randn(n_test, p_max)

y_train = X_train_full[:, :k_true] @ beta_true + sigma * np.random.randn(n_train)
y_test = X_test_full[:, :k_true] @ beta_true + sigma * np.random.randn(n_test)

checkpoints = [1,5,10,15,18,19,20,21,22,25,30,50,100,300,800,1800,3000]

for p in checkpoints:
    Xtr = X_train_full[:, :p]
    Xte = X_test_full[:, :p]
    beta_hat, *_ = np.linalg.lstsq(Xtr, y_train, rcond=None)
    train_pred = Xtr @ beta_hat
    test_pred = Xte @ beta_hat
    train_mse = np.mean((train_pred - y_train) ** 2)
    test_mse = np.mean((test_pred - y_test) ** 2)
    print(f"p={p:4d}  train_MSE={train_mse:8.4f}  test_MSE={test_mse:9.4f}  ||beta_hat||={np.linalg.norm(beta_hat):8.3f}")

Running this produces the following (actual output, this seed, unaltered):

p=   1  train_MSE=  6.1133  test_MSE=   9.2554  ||beta_hat||=   2.424
p=   5  train_MSE=  0.1396  test_MSE=   0.4566  ||beta_hat||=   3.814
p=  10  train_MSE=  0.1016  test_MSE=   0.4997  ||beta_hat||=   3.811
p=  15  train_MSE=  0.0600  test_MSE=   0.6405  ||beta_hat||=   3.767
p=  18  train_MSE=  0.0347  test_MSE=   1.2028  ||beta_hat||=   3.784
p=  19  train_MSE=  0.0281  test_MSE=   5.6334  ||beta_hat||=   5.187
p=  20  train_MSE=  0.0000  test_MSE= 974.3024  ||beta_hat||=  27.675
p=  21  train_MSE=  0.0000  test_MSE=  11.7183  ||beta_hat||=   5.138
p=  22  train_MSE=  0.0000  test_MSE=   7.4218  ||beta_hat||=   3.315
p=  25  train_MSE=  0.0000  test_MSE=   6.6445  ||beta_hat||=   3.210
p=  30  train_MSE=  0.0000  test_MSE=   7.7229  ||beta_hat||=   3.048
p=  50  train_MSE=  0.0000  test_MSE=  13.1752  ||beta_hat||=   2.512
p= 100  train_MSE=  0.0000  test_MSE=  18.2114  ||beta_hat||=   1.662
p= 300  train_MSE=  0.0000  test_MSE=  18.1980  ||beta_hat||=   0.929
p= 800  train_MSE=  0.0000  test_MSE=  18.6500  ||beta_hat||=   0.587
p=1800  train_MSE=  0.0000  test_MSE=  18.7313  ||beta_hat||=   0.403
p=3000  train_MSE=  0.0000  test_MSE=  18.9618  ||beta_hat||=   0.312

Trace what happens step by step. From p = 5 to p = 18, adding features barely helps and slowly hurts, test MSE drifts from 0.46 up to 1.20, a mild version of the classical overfitting slope, while train MSE keeps shrinking. At p = 20 = n, the design matrix is exactly square: there is only one solution that fits all 20 points exactly, and that solution is catastrophic, test MSE explodes to 974, roughly 2,100 times the error at the best classical model (p = 5). Immediately past the threshold, at p = 21, the blow-up nearly vanishes, down to 11.7, and by p = 100 onward test MSE has settled into a plateau around 18-19 that barely moves as p climbs all the way to 3000. That plateau, spike, and recovery, all in a model with a closed-form solution, is the same qualitative shape production teams see when they scale a neural network's width past the point where it starts to memorize its training set.

The spike has an exact explanation, not just a description. At p = n = 20, with X square, the OLS solution is beta_hat = X^{-1} y = beta_true + X^{-1} eps, so the estimation error is entirely X^{-1} eps. Averaging the squared test error over a fresh isotropic test point gives exactly E[||beta_hat - beta_true||^2] = sigma^2 · trace((X^T X)^{-1}) = sigma^2 · sum(1 / s_i^2), where s_i are the singular values of the training design matrix. The smaller the smallest singular value, the more that one term dominates the sum and the more the estimator's variance explodes.

import numpy as np
np.random.seed(0)
n_train = 20; p_max = 3000
X_train_full = np.random.randn(n_train, p_max)
X20 = X_train_full[:, :20]
s = np.linalg.svd(X20, compute_uv=False)
print(np.round(s, 3))
print('smallest sv:', s.min(), ' condition number:', s.max() / s.min())
print('sum 1/s^2:', np.sum(1 / s**2))

Output:

[8.417 7.755 6.931 6.541 5.865 5.649 4.77  4.403 4.291 3.978 3.702 2.916
 2.727 2.131 1.869 1.716 1.441 0.709 0.288 0.024]
smallest sv: 0.023647828425396644  condition number: 355.9417269078388
sum 1/s^2: 1804.2818316473918

The smallest singular value of this particular 20-by-20 matrix is 0.0236, giving a condition number of 356 and sum(1/s_i^2) = 1804.28. Plugging in, the expected excess test MSE at this threshold is 0.25 × 1804.28 ≈ 451, plus the irreducible noise floor of 0.25, about 451.3 in total. The observed value, 974.3, is roughly 2.2 times that expectation, which is exactly what a theory built on the smallest singular value predicts should happen: that quantity is heavy-tailed and dominated by a single near-zero singular value, so any one random seed will land far from the mean. This is also why published double-descent figures are almost always averaged over many random seeds near the threshold, a single run there is close to meaningless in magnitude, only in shape.

The plateau has its own explanation, and it is the more instructive one. As p grows far past n, the minimum-norm interpolating solution has to spread its "explaining power" across an ever-larger number of directions, most of which carry no real signal in this experiment (features 6 through 3000 are pure noise by construction). The theoretical floor for a model that gives up entirely on identifying the true 5-dimensional signal and instead predicts the unconditional mean is ||beta_true||^2 + sigma^2 = 16.5 + 0.25 = 16.75. The observed plateau, 18.2 to 19.0, sits just above that floor, the gap explained by the fact that n = 20 is a small sample, so even the asymptotic minimum-norm estimator retains some residual estimation noise. The interpolator has, in effect, given up on recovering the true coefficients precisely and settled for a safe, low-variance, near-featureless prediction, escaping the catastrophic variance of the threshold but not fully recovering the accuracy available at the classical sweet spot (0.46). That distinction, escaping catastrophe versus fully recovering accuracy, is exactly the hinge the next section turns on.

Why overparameterization doesn't have to fail: four mechanisms

The toy experiment above used features with a flat, isotropic spectrum, meaning every direction in feature space is equally likely to look informative to the estimator, and no direction is genuinely favored. That turns out to be close to the worst possible setting for a full recovery in the second descent. Bartlett, Long, Lugosi, and Tsigler (PNAS, 2020) worked out precisely when a minimum-norm interpolator can fit training noise exactly and still generalize almost as well as the best possible model, a phenomenon they named "benign overfitting." Their condition is about the spectrum of the feature covariance: benign overfitting needs a small number of directions that carry most of the signal's energy, plus a very large number of directions, individually almost negligible, that can each absorb a tiny sliver of the noise essentially for free. An isotropic spectrum, every direction equally weighted, has no such structure to exploit, which is precisely why the toy model's second descent stalled well above the classical sweet spot instead of beating it. Real data, and real trained representations, are rarely isotropic: image and language features concentrate most of their energy in a comparatively low-dimensional subspace, which is exactly the kind of decaying spectrum benign overfitting requires.

Second, the optimizer itself is not neutral among the many solutions that fit the training data exactly. Soudry et al. (2018) showed that gradient descent on separable classification data converges, in direction, to the maximum-margin solution, not an arbitrary one, even though many other solutions fit the training labels equally well. Gunasekar et al. (2018) generalized this into a broader picture of "implicit regularization": the geometry of the optimization algorithm itself, not any explicit penalty term, selects a specific, usually simple, member of the (huge) set of zero-training-error solutions. Overparameterization is not a liability here, it is what creates the large set of equally-fitting solutions for the optimizer's implicit bias to choose among; a model with barely enough capacity to interpolate, by contrast, has only one such solution available, fragile or not.

Third, the Neural Tangent Kernel framework (Jacot, Gabriel, and Hongler, NeurIPS 2018) shows that in the limit of infinite width, a neural network's training dynamics become linear in its parameters and equivalent to kernel regression against a fixed, architecture-determined kernel. This reframes "why do enormous networks generalize" as a question already well understood in classical kernel methods, generalization is governed by the smoothness properties of that kernel and by an implicit norm-minimization very similar to the minimum-norm mechanism in the worked example above, not by counting parameters at all.

Fourth, Frankle and Carbin (ICLR 2019) proposed the lottery ticket hypothesis: a randomly initialized dense network contains a much smaller subnetwork that, trained in isolation from the same initialization, matches the full network's test accuracy. On this view, part of what extra width buys you is not "more capacity used," it is more chances for a lucky, well-initialized sparse subnetwork to exist inside the larger one, which the optimizer can then find and effectively train while leaving the rest near-inert. Consistent with this, Dziugaite and Roy (UAI, 2017) produced the first non-vacuous PAC-Bayes generalization bounds for real trained deep networks by exploiting exactly this idea in reverse, that trained solutions tend to sit in wide, flat regions of the loss landscape, which a PAC-Bayes analysis can certify generalize well, independent of the raw parameter count.

The common thread across all four mechanisms is the same: parameter count alone was never the right measure of "complexity" for a network trained by gradient-based optimization. The relevant complexity is some norm, margin, or spectral property of the specific solution the optimizer actually finds, and overparameterization is beneficial precisely because it gives the optimizer more low-complexity solutions to land on, not because it gives the model more raw capacity to memorize with.

Common misconception

The misconception this chapter set out to correct: "if a model has more parameters than training examples, it must have memorized noise and will generalize poorly." Classical statistical learning theory does predict exactly this, and for models fit without any implicit or explicit bias toward simple solutions, it can be true, the catastrophic spike at p = n = 20 in the worked example is a real instance of it. But it is not a fact about parameter counts, it is a fact about what happens at the specific, fragile capacity where exactly one solution exists that fits the training data. Push capacity further past that point, so that many solutions fit the data exactly, and whether the model generalizes well becomes a question about which of those many solutions the optimizer is biased toward finding, not a question about how many parameters it has. "Overparameterized" is not synonymous with "overfit"; it is, more often in modern practice, the precondition for the optimizer's implicit regularization to have room to operate at all.

Active recall

Attempt each question before reading its answer.

1. State what the Universal Approximation Theorem guarantees, and name one thing it explicitly does not guarantee.
2. In the worked linear-regression experiment, why did test MSE spike specifically at p = 20 rather than gradually worsening as p approached 20 from below?
3. Why did train MSE hit exactly 0.0000 for every p ≥ 20, while test MSE kept changing substantially across that same range?
4. State, qualitatively, the covariance-spectrum condition Bartlett et al. (2020) identify for "benign overfitting," and explain why the worked example's noise features violate it.
5. What does the infinite-width limit in the Neural Tangent Kernel framework simplify about training dynamics, and why does that simplification help explain generalization?
6. Ripple-effect: suppose the label noise in the worked example is increased from sigma = 0.5 to sigma = 2.0, with the same design matrices X_train_full and X_test_full (unchanged singular values). Trace the effect on (a) the expected height of the spike at p = 20, (b) the plateau value for very large p, and (c) the location of the classical sweet spot.

Answers.

1. The theorem (Cybenko, 1989; Hornik et al., 1989; Hornik, 1991) guarantees that a sufficiently wide (or, per Lu et al., 2017, sufficiently deep) network with a non-polynomial activation can approximate any continuous function on a compact domain to arbitrary accuracy. It does not guarantee that gradient descent will find the weights that achieve that approximation (it is an existence proof, not a constructive algorithm), nor does it bound how large the network needs to be to hit a given accuracy for a generic function (that bound can be exponential in the input dimension), nor does it say anything about generalization to data outside the approximated domain.

2. At exactly p = n = 20, the design matrix is square, so there is exactly one exact solution to the 20 training equations, and that unique solution equals beta_true + X^{-1}eps: whatever the specific noise realization was, it gets fully absorbed into the coefficient estimates, amplified by X^{-1}, whose size is controlled by the smallest singular value of X (0.0236 here, giving a condition number of 356). Below p = 20, the system is overdetermined and least squares averages over the noise rather than absorbing it point-for-point, so there is no equivalent amplification; the degradation is gradual, matching the mild rise seen from p = 5 to p = 19 (0.46 to 5.63).

3. Once p ≥ n, the system is underdetermined (or exactly determined at p = 20): there exists at least one exact solution to the 20 training equations, and numpy.linalg.lstsq returns one such solution (the minimum-norm one), so training residuals are exactly zero by construction, regardless of how badly that solution generalizes. Test error, evaluated on 200 fresh points the fit never saw, has no such guarantee: it depends on how well the specific minimum-norm solution recovers the true 5-dimensional signal amid an ever-larger number of pure-noise directions, which is precisely what varies from 974 down to a plateau near 18-19 as p grows.

4. Benign overfitting requires the feature covariance's spectrum to have a small number of large eigenvalues, capturing most of the signal's energy cheaply, and a very large number of small eigenvalues, each able to absorb an individually negligible amount of noise. The worked example uses i.i.d. Gaussian features for every one of the up to 3000 columns, an isotropic (flat) spectrum with no such large/small split, so there is no cheap place for the interpolator to hide the noise without also diluting the true signal, which is exactly why the second descent plateaus around 18-19 rather than dropping back toward the 0.46 achieved by the classical sweet spot.

5. In the infinite-width limit, the network's output becomes a linear function of its parameters throughout training (the tangent kernel itself stops changing), which turns the generally non-convex neural network training problem into an equivalent convex kernel regression problem. Because kernel regression's generalization behavior is already well understood analytically, this equivalence gives a rigorous account of when and why very wide networks generalize, in terms of the smoothness of the induced kernel, without needing to solve the original non-convex optimization problem directly.

6. The design matrices and hence all 20 singular values are unchanged, only sigma changes, from 0.5 to 2.0, so sigma^2 scales by (2.0/0.5)^2 = 16. (a) The expected excess risk at the threshold is sigma^2 · sum(1/s_i^2), so it scales by the same factor 16: from an expected ~451 to an expected ~7,217; any single realized run, given the roughly 2x variance seen at the original noise level, could plausibly land well into the thousands. (b) The plateau value is ||beta_true||^2 + sigma^2 = 16.5 + sigma^2: it rises only modestly, from 16.75 to 20.5, because the fixed signal term (16.5) dominates the noise term at both noise levels. (c) The classical sweet spot shifts toward a smaller p: with noisier labels, each additional feature buys less signal recovery for the variance it costs, so the bias-variance-optimal capacity is lower than before, and its achievable risk floor rises, since there is now more noise for even the best-tuned classical model to average away with the same 20 training points.

Think About It

Think about this: How would you explain how neural networks learn to a friend who has never seen a computer? What real-world analogy would you use? Imagine you had to build a system using these concepts — what would be your first step? Try this: before moving on, write down three things you learned and one question you still have.

Key Takeaways — Summary and Recap

Let us recap what we covered: the core ideas behind how neural networks learn, how they connect to real-world applications, and why they matter for your journey in computer science. Remember these key points as you move forward. For competitive exam preparation (CBSE, JEE, BITSAT), focus on understanding the WHY behind each concept, not just the WHAT.