Gradient Descent & Regularization

The closed-form solution of ordinary least squares (OLS) is elegant, but it does not scale to huge feature counts, does not exist for most other models, and says nothing about controlling overfitting. This lesson adds the two tools that do: gradient descent — the optimization engine behind nearly all of modern ML — and regularization — the standard brake on model complexity.

Gradient descent

The idea is the one that closed the previous lesson. To find the minimum of a loss \(J(w)\), instead of solving \(\nabla J = 0\) on paper, walk downhill: measure the slope where you stand, take a step against it, repeat.

\(\eta\) is the learning rate, the size of the step. The gradient \(\nabla_w J\) points in the direction where the loss grows fastest, so it is the minus sign that makes the method descend rather than climb.

Only one thing is missing before this becomes code: knowing what \(\nabla_w J\) actually is. It is a list of partial derivatives, one per parameter, and it is worth assembling that list by hand before accepting the matrix form ready-made.

Where the gradient comes from

Start with the smallest case, the two-parameter line. Write the loss as a mean, so the scale does not depend on \(n\):

Each term is the square of a residual, \(e_i^2\). By the chain rule, the derivative of \(e_i^2\) with respect to any parameter is \(2e_i\) times the derivative of \(e_i\) with respect to that parameter. So everything reduces to knowing how the residual reacts to each weight — and that is easy, because \(e_i\) depends on the weights linearly:

The \(-1\) appears because raising the intercept by one unit lowers every residual by one unit. The \(-x_i\) appears because raising the slope by one unit lowers the residual by \(x_i\) — the further a point sits from the origin, the more it feels a change in slope.

Substituting:

And so on: parameter \(j\)

With \(d\) features the model is \(\hat{y}_i = w_0 + w_1x_{i1} + \dots + w_dx_{id}\). Adopt the convention \(x_{i0} = 1\) — the column of ones already sitting in the design matrix — and the intercept stops being a special case: it becomes just the weight on a constant column. The derivative of the residual then always has the same shape,

and a single formula covers every parameter:

Notice the first row is not an exception to the last: it is the last, with \(x_{i0}=1\). It is the same column-of-ones trick that turned \(\hat{y} = w_0 + w_1x\) into \(\hat{y} = Xw\).

Stacking them: the matrix form

Now collect those \(d+1\) derivatives into a vector. Entry \(j\) is \(\sum_i e_i x_{ij}\), which is the \(j\)-th column of \(X\) dotted with the residual vector. Doing that for every column at once is multiplying by \(X^\top\):

There is no magic in this formula: it is the table above, written compactly. That is why one line of numpy computes the gradient of a model with a million parameters — the structure is exactly the one for \(w_0\) and \(w_1\), repeated.

Two readings come out of this, and both travel far.

First: set the gradient to zero and you recover the normal equations, \(X^\top(y - Xw) = 0\) — which is the closed-form solution from last lesson. The closed form and the gradient are the same statement by two routes: one solves \(\nabla J = 0\) in a single move, the other walks there.

Second: look at what the gradient is made of — a residual-weighted sum. A point the model already predicts well enters with \(e_i \approx 0\) and barely votes; a point it badly misses dominates the step. So the rule \(w \leftarrow w - \eta\,\nabla J\) reads, in plain words: move each weight in the direction the errors point, by an amount proportional to how large those errors are. That sentence stays true, word for word, when the model is a 100-layer network and \(\nabla J\) comes from backpropagation.

One step, by hand

Three points — \((1,2), (2,4), (3,7)\) — starting from \(w_0 = w_1 = 0\), with \(\eta = 0.1\).

With both weights at zero the line predicts 0 everywhere, so the residual is \(y\) itself:

Drop these into the two formulas we just derived:

Both came out negative, which makes sense: both weights are too low, and raising either one reduces the loss. The step therefore raises both:

In that one step the loss falls from \(J = 23\) to \(J = 0.625\). Repeating the procedure:

Read the table slowly, because it shows three things at once.

The first step is enormous: the loss drops from 23 to 0.63, nearly all the work. Then progress crawls — 0.344, then 0.327. That shape is the rule, not the exception, and it is why watching \(J\) fall fast early on says almost nothing about having converged.

The gradients change sign between iterations 1 and 2. That is the step having overshot the minimum and being walked back — the same bouncing the learning-rate simulator shows below, here in miniature.

And the two parameters do not move at the same pace: at iteration 0 the derivative in \(w_1\) is 2.4 times the one in \(w_0\), because each residual enters the \(w_1\) sum multiplied by \(x_i\), which reaches 3. Hold on to that asymmetry — it is the subject of the landscape just below.

Left running to 2000 iterations it reaches \(w_0 = -0.667\) and \(w_1 = 2.5\), exactly what the normal equations give — the same answer, arrived at by a slower road.

The landscape you are walking down

So far these have been numbers in a table. It is worth seeing what they are geometrically, because the picture explains both the slowness and its cure at once.

Think of \(J(w_0, w_1)\) as a landscape: every point of the plane is a pair of weights, that is, a candidate line, and the height there is the error that line makes. Minimising means hunting for the bottom of the valley. The gradient is the direction of steepest ascent at that point, so \(-\nabla J\) is the direction of steepest descent.

The left panel shows the contours of that landscape, with the closed-form optimum marked. On the right, the same weights drawn as a line through the data. Click on the left to drop a starting guess, then keep pressing Step:

Now switch on centre x and start again from a similar point. Same data, same algorithm, completely different behaviour.

Uncentred, the contours form a long, narrow valley. The reason is the same asymmetry you saw in the table: \(w_0\) and \(w_1\) become coupled, because tilting a line whose \(x\) values all sit far from zero also raises or lowers its height. The gradient points across the valley instead of along it, so the path zigzags and drags. Centred, the valley becomes a round bowl and the same steps go nearly straight to the answer.

The learning rate

• \(\eta\) too small → tiny steps, painfully slow convergence;
• \(\eta\) too large → steps overshoot the minimum; the loss oscillates or diverges;
• practical recipe: try \(\eta \in \{10^{-3}, 10^{-2}, 10^{-1}\}\), monitor the training-loss curve — it should fall smoothly.

Feel it yourself. On the loss below the threshold sits exactly at \(\eta = 1\), which makes the whole story visible in one place: at \(\eta = 0.5\) the step lands on the minimum in one move; below that it creeps in; between 0.5 and 1 it overshoots but still closes in, bouncing from side to side; at exactly \(\eta = 1\) it bounces between the same two points for ever, never getting closer; past 1 it explodes. Set \(\eta = 1.05\) and step through it:

Batch, stochastic, and mini-batch

from sklearn.linear_model import SGDRegressor
model = SGDRegressor(loss='squared_error', penalty='l2', alpha=1e-4,
                     learning_rate='invscaling', max_iter=1000)

The same loop — with different losses — trains logistic regression, SVMs, and neural networks. Learn it once, reuse it everywhere.

From lines to curves: polynomial features

Linear regression is linear in the parameters, not necessarily in the inputs. Expanding features to powers, \(x \mapsto (x, x^2, \dots, x^p)\), fits polynomials with the same OLS machinery:

from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures

model = make_pipeline(PolynomialFeatures(degree=3), LinearRegression())

But flexibility cuts both ways:

Degree 1 underfits — too rigid to follow the sine. Degree 15 overfits — 16 parameters chase 25 noisy points, producing wild oscillations. The right panel keeps all 16 parameters but adds a Ridge penalty: the curve relaxes back to the signal. That is regularization at work.

Underfitting and overfitting, measured

The figure convinces the eye. It is worth measuring, because the numbers show something the picture hides.

The data are 100 points of \(y = 0.5x^2 + x + 2\) with noise of standard deviation 3, for \(x\) between \(-3\) and \(3\). Three models: degree 1, degree 2 — the right degree — and degree 30. For each, the error on its own training set beside the error under five-fold cross-validation:

Degree 30 has the lowest training error of the three. Anyone picking a model by that column would pick the worst one.

The most instructive column, though, is the last. Look at the five folds of degree 30 one by one:

On three of the five it does as well as degree 2. On one, it explodes. Overfitting is not being wrong all the time — it is being wrong unpredictably, which is why the spread across folds gives the problem away before the mean does.

That yields a pocket diagnostic:

There is also a floor. The noise has standard deviation 3, and that part of \(y\) does not depend on \(x\), so no model can remove it. Degree 2 lands near that floor, which is all anyone can ask of a model.

"Too complex" is relative to how much data you have

Complexity is not a property of the model alone: it is the relationship between the number of parameters and the number of points.

A degree-10 polynomial has 11 free parameters — ten powers plus the intercept. As long as the training points do not outnumber those, the curve can pass exactly through every one of them:

Zero training error is not a sign of success: it is a sign that nothing was left over. With ten points and eleven parameters the model never had to learn anything about the shape of the curve — memorizing was enough. From the eleventh point on it has to start choosing, and the training error lifts off the floor.

Notice what this implies: the same degree-10 model is too complex for 10 points and perfectly reasonable for 100. There is no "right" degree in the abstract, only a right degree for a quantity of data. The tool that makes this dependence visible is the learning curve — training and validation error against training-set size — covered in model selection.

Regularization

Instead of restricting the number of parameters, penalize their magnitude — add a complexity term to the loss:

Ridge (L2) — Tikhonov, 1943; Hoerl & Kennard, 1970

• shrinks all coefficients smoothly toward zero (never exactly zero);
• spreads weight across correlated features — the standard cure for multicollinearity;
• closed form still exists: \(\hat{w} = (X^\top X + \alpha I)^{-1} X^\top y\) — the \(\alpha I\) makes the matrix invertible even with collinear features.

Lasso (L1) — Tibshirani, 1996

• the absolute-value penalty has corners at zero: solutions land exactly at zero for weak features;
• performs automatic feature selection — the surviving nonzero coefficients name the features that matter;
• among a group of highly correlated features, it tends to keep one arbitrarily and zero the rest.

Elastic Net blends both penalties (l1_ratio) — a robust default when features are many and correlated.

from sklearn.linear_model import Ridge, Lasso, ElasticNet

Ridge(alpha=1.0)
Lasso(alpha=0.1)
ElasticNet(alpha=0.1, l1_ratio=0.5)

Why L1 zeroes and L2 does not

Both bullet lists above assert the difference. The reason is geometric, and it is worth seeing rather than believing.

Write the penalty as a budget instead: minimise the squared error subject to \(\sum |w_j| \le t\) (lasso) or \(\sum w_j^2 \le t^2\) (ridge). The two formulations are equivalent — every \(\alpha\) has a matching \(t\). Now the picture is two objects: the elliptical contours of the error, centred on the OLS solution, and the feasible region.

The constrained solution is where the ellipse, growing outward from the OLS optimum, first touches the region:

The circle has no corners, so the touch point can be anywhere on it — and lands on an axis only by coincidence. The diamond has four corners, and they sit exactly on the axes. A smooth curve expanding toward a pointed region tends to meet a point first, and meeting a corner is a coefficient being exactly zero.

Drag the OLS optimum around and watch how often the diamond is touched at a corner while the circle almost never is. Lasso does not select features by thresholding — it selects them because of the shape of its constraint.

The one-dimensional case, in closed form

With a single standardised feature the two penalties have explicit solutions, and they say the same thing in algebra:

With \(\alpha = 1\):

Ridge divides; division never reaches zero. Lasso subtracts a fixed amount and clips at zero — so every coefficient below the threshold \(\alpha/2\) does not shrink, it disappears. That subtraction is the corner of the diamond, written as a formula.

Watching the path

Run the penalty across a range of \(\alpha\) and plot every coefficient. Eight features here: three carry real signal, three are pure noise, and two are near-copies of each other.

Lasso removes the noise features first — exactly the selection you hoped for — and the count keeps dropping. Ridge never drops below eight at any \(\alpha\): the noise coefficients get small and stay.

Then switch to the correlated pair \(c_1, c_2\), whose true coefficients are both 2. At \(\alpha = 2\) lasso reports (1.89, 0.98) — it is starting to pick one and discard the other, and which one is close to arbitrary. Ridge reports (1.02, 0.97): it splits the weight evenly. That is the whole reason ridge is the standard answer to multicollinearity and lasso is not.

The knob α

\(\alpha\) trades data fidelity against coefficient size:

• \(\alpha \to 0\): plain OLS (no brake);
• \(\alpha \to \infty\): all coefficients crushed to ~0, model predicts the mean (full brake);
• the right \(\alpha\) is not known in advance — it is chosen by cross-validation (RidgeCV, LassoCV, or a grid search).

Class materials

References

• Tikhonov, A. N. "On the stability of inverse problems." C. R. (Doklady) Acad. Sci. URSS, n. Ser. 39 (1943), 176–179. zbMATH record
• Robbins, H.; Monro, S. "A Stochastic Approximation Method." Annals of Math. Statistics 22 (1951). DOI
• Hoerl, A. E.; Kennard, R. W. "Ridge Regression: Biased Estimation for Nonorthogonal Problems." Technometrics 12 (1970). DOI
• Tibshirani, R. "Regression Shrinkage and Selection via the Lasso." JRSS B 58 (1996). DOI

The full course bibliography is on the references page.

Quiz