Linear Regression Explained: Normal Equation to Gradient Descent

Here is a strange fact. You can find the best fitting line through a cloud of points in a single step, with one formula. Or you can start with a terrible line and improve it slowly, taking thousands of tiny steps. Both land on the same answer. Today we find out why, and when each one wins.

• Normal equation: one shot
• Gradient descent: step by step

Last episode we saw that a model must generalize, and that too much flexibility leads to overfitting. Now we meet the simplest model that actually learns from data: linear regression. It is simple enough to solve exactly, which makes it the perfect place to understand every moving part.

• Last time: overfitting and bias-variance
• Today: our first real model

We will define the model, measure its errors with squares you can actually see, and compare three loss measures. Then we derive the normal equation using matrix calculus, check it with three points by hand, derive the gradient, write NumPy code, and finish with scaling, polynomials, and the assumptions that can break everything.

• Hypothesis and residuals
• Normal equation derivation
• Gradient descent and NumPy
• Scaling, polynomials, assumptions

Suppose x is the size of a house and y is its price. Our hypothesis says the predicted price, written y hat, equals w times x plus b. The weight w is the slope, how much the price rises per extra unit of size. The bias b is the intercept, the prediction when x is zero. Learning means choosing w and b.

• y-hat = w x + b
• w: slope, b: intercept

For each data point, the residual is the actual value minus the predicted value. Graphically it is the vertical gap between the point and the line. Points above the line give positive residuals, points below give negative ones. A good line makes these gaps small overall, but we need one number to say how small.

• Residual = actual minus predicted
• Vertical gaps, not perpendicular

Adding raw residuals fails, because positives and negatives cancel. A flat line through the middle of a rising cloud can have a residual sum of exactly zero while being a terrible fit. We need every error to count as positive. We could take absolute values, or we could square them.

• Positives cancel negatives
• Bad lines can score zero

Here is the picture to remember. Draw a square on each residual, with the residual as its side. The area of each square is the squared error. Linear regression finds the line that makes the total area of all those squares as small as possible. Big misses make huge squares, so they dominate the total.

• Area of square = error squared
• Minimize total area

Mean squared error is one over n, times the sum from i equals one to n of y sub i minus y hat sub i, squared. Here n is the number of examples, y sub i is the true value, and y hat sub i is the prediction. We divide by n so the score does not grow just because we have more data.

• MSE = (1/n) sum (y_i - y-hat_i)^2
• n = number of points

Try three errors: one, two, and three. Squared they become one, four, and nine. The sum is fourteen. Divide by three and the mean squared error is about four point six seven. Notice the error of three contributes nine, more than the other two combined. Squaring punishes big mistakes much harder.

• Errors: 1, 2, 3
• Squares: 1, 4, 9
• MSE = 14/3 = 4.67

Mean absolute error averages the absolute residuals. For our errors that is two. Root mean squared error takes the square root of the MSE, giving about two point one six. RMSE brings the number back into the original units, dollars instead of dollars squared, which makes it easier to interpret. It is always at least as large as MAE.

• MAE: average absolute error
• RMSE: square root of MSE
• RMSE has original units

The choice matters. Squared loss is smooth, so calculus works nicely, but a single outlier can drag the line toward it. Absolute loss is more robust, yet it has a sharp corner at zero, which makes it harder to solve in closed form. Squared loss also corresponds to assuming Gaussian noise, which we saw in the maximum likelihood episode.

• MSE: smooth, outlier-sensitive
• MAE: robust, kink at zero
• Squared loss matches Gaussian noise

Real problems have many features: size, bedrooms, age of the house. Multiple linear regression predicts with b plus w one x one plus w two x two, and so on. With two features the model is a flat plane through three dimensional space. With a hundred features it is a flat hyperplane we cannot draw, but the math is identical.

• y-hat = b + w1 x1 + w2 x2 + ...
• Fit a plane, not a line

To write this compactly, stack the data into a matrix X with one row per example. Add a column of ones on the left so the bias becomes just another weight. Collect all parameters into one vector called theta. Now every prediction at once is simply X times theta. One short expression replaces a loop over every example.

• Add a column of ones
• X is n by d plus one
• theta stacks bias and weights

The residual vector is r equals y minus X theta. The sum of squared errors is r transpose times r, because multiplying a row of residuals by a column of residuals gives each residual times itself, added up. So we minimize J of theta equals y minus X theta, transposed, times y minus X theta.

• Residual vector r = y - X theta
• SSE = r transpose r

The loss, as a function of theta, is a smooth bowl. At the bottom of a bowl the slope is flat in every direction, meaning the gradient equals the zero vector. So our plan is to compute the gradient of J with respect to theta, set it to zero, and solve for theta. We just need matrix calculus.

• Find theta where gradient is zero
• Bowl-shaped loss has one minimum

First expand the product using the rule that the transpose of a product reverses order. We get y transpose y, minus theta transpose X transpose y, minus y transpose X theta, plus theta transpose X transpose X theta. The two middle terms are single numbers and are transposes of each other, so they are equal. That gives minus two theta transpose X transpose y.

• (y - X theta)^T (y - X theta)
• = y^T y - 2 theta^T X^T y + theta^T X^T X theta

Now differentiate term by term. The first term does not involve theta, so its gradient is zero. For the second, the gradient of theta transpose a is just a, giving minus two X transpose y. For the third, the gradient of theta transpose A theta is two A theta when A is symmetric, and X transpose X is symmetric. That gives two X transpose X theta.

• d/dtheta of y^T y = 0
• d/dtheta of theta^T a = a
• d/dtheta of theta^T A theta = 2 A theta

Set the gradient to zero: two X transpose X theta minus two X transpose y equals zero. Move the second term across and cancel the twos to get X transpose X theta equals X transpose y. These are the normal equations. Multiply both sides by the inverse of X transpose X, and we get theta equals the inverse of X transpose X, times X transpose y.

• Gradient = 2 X^T X theta - 2 X^T y = 0
• X^T X theta = X^T y
• theta = (X^T X)^-1 X^T y

There is a beautiful picture behind this. The columns of X span a flat subspace. The vector y usually sits outside it. The best prediction is the projection of y onto that subspace, its shadow, and the residual points straight out at a right angle. The normal equations simply say X transpose times the residual equals zero.

• Prediction = projection of y
• Residual is orthogonal to columns of X

How do we know this point is a minimum and not a maximum? The second derivative, the Hessian, is two X transpose X. This matrix is positive semi-definite, as we saw in the eigenvalue episode, so the loss curves upward everywhere. The surface is convex, with no local minima to get stuck in.

• Hessian = 2 X^T X
• Positive semi-definite
• Convex: no local traps

Let us run the formula by hand on three points: one comma one, two comma two, and three comma two. The design matrix X has rows one one, one two, and one three, with the first column of ones for the bias. The target vector y is one, two, two.

• (1, 1), (2, 2), (3, 2)
• Build X and y

X transpose X has entries: three, which is the count of points; six, the sum of the x values; and fourteen, the sum of x squared, which is one plus four plus nine. X transpose y has two entries. The first is one plus two plus two, which is five. The second is one times one, plus two times two, plus three times two, which is eleven.

• X^T X = [[3, 6], [6, 14]]
• X^T y = [5, 11]

The determinant of the two by two matrix is three times fourteen minus six times six, which is six. The inverse swaps the diagonal, negates the off diagonal, and divides by six. Multiply by X transpose y: the first entry is seventy minus sixty six, over six, giving two thirds. The second is minus thirty plus thirty three, over six, giving one half.

• det = 3x14 - 6x6 = 6
• Inverse = (1/6)[[14, -6], [-6, 3]]
• theta = [2/3, 1/2]

The best line is y hat equals zero point six six seven plus point five x. Predictions are about one point one six seven, one point six six seven, and two point one six seven. The residuals are minus point one six seven, plus point three three three, and minus point one six seven. They sum to zero, as the normal equations promise, because the bias column is orthogonal to the residual.

• Line: y-hat = 0.667 + 0.5 x
• Residuals: -0.167, 0.333, -0.167
• Residuals sum to zero

Now score it. The squared errors are point zero two eight, point one one one, and point zero two eight, totaling about point one six seven. Divide by three for an MSE of point zero five five six. The square root gives an RMSE of about point two three six. The MAE is point two two two. No other line can beat that total area of squares.

• SSE = 0.1667
• MSE = 0.0556
• RMSE = 0.236, MAE = 0.222

The normal equation is exact, but it has costs. Inverting a matrix with d features takes roughly d cubed operations. If features are duplicated or perfectly correlated, X transpose X is singular and has no inverse. And with huge datasets, forming that matrix may not even be practical. This is where gradient descent comes in.

• Inverting is O(d^3)
• X^T X may be singular
• Data may not fit in memory

The derivation is the one we just did, with a one over n in front. The gradient of the mean squared error is two over n, times X transpose, times X theta minus y. Read it simply: take the prediction errors, weight them by each feature, and average. The gradient points uphill, so we step in the opposite direction.

• J = (1/n)(X theta - y)^T (X theta - y)
• Gradient = (2/n) X^T (X theta - y)

The update rule is theta becomes theta minus alpha times the gradient. Alpha is the learning rate, the size of each step. Repeat many times and theta slides down the bowl to the same minimum the normal equation finds directly. Because the loss is convex, there is only one bottom to reach.

• theta = theta - alpha x gradient
• alpha: learning rate

Use our three points. Start at theta equals zero zero, so every prediction is zero and the errors, prediction minus truth, are minus one, minus two, minus two. The MSE is three. The gradient is two thirds times X transpose times those errors: minus three point three three and minus seven point three three. With alpha point one, theta moves to point three three and point seven three.

• Start theta = [0, 0]
• Gradient = [-3.33, -7.33]
• alpha = 0.1: theta = [0.33, 0.73]

Choosing alpha is the classic tuning problem. Too small, and training crawls. Just right, and the loss falls smoothly. Too large, and each step overshoots the valley, bouncing higher and higher until the loss explodes. A good habit is to plot loss against iteration. It should fall steadily. If it rises, shrink the learning rate.

• Too small: crawls
• Just right: smooth
• Too large: diverges

The gradient we wrote uses every example, called batch gradient descent. Stochastic gradient descent uses one random example per step, which is noisy but cheap. Mini batch uses a small group, such as thirty two or sixty four examples, balancing speed and stability. This is what scales linear models to enormous datasets, and what neural networks use later.

• Batch: all data per step
• Stochastic: one example
• Mini-batch: a small group

In NumPy, first build X by stacking a column of ones next to the feature. Then solve the normal equations with np dot linalg dot solve, passing X transpose X and X transpose y. Prefer solve over explicitly inverting, because it is faster and numerically more stable. On our three points it prints point six six seven and point five.

• Add ones column
• Use solve, not inv

Gradient descent is just as short. Start theta at zeros. In a loop, compute the gradient as two over n times X transpose times X theta minus y, then subtract alpha times it. There is no inner loop over examples, because matrix multiplication handles them all at once. A few thousand iterations later, theta matches the closed form.

• Vectorized: no loops over rows
• Loop only over iterations

Now compare with scikit-learn. LinearRegression fits the model and exposes intercept underscore and coef underscore. It does not need a ones column, since it handles the intercept itself. Run it on our data and you get point six six seven and point five, matching both of our from scratch versions. Understanding the math means the library is no longer a black box.

• LinearRegression: closed-form
• SGDRegressor: gradient-based
• Same coefficients

Suppose one feature is house size in the thousands and another is bedrooms from one to five. The loss valley becomes a long narrow canyon, and gradient descent zig zags across it painfully slowly. Standardizing each feature, subtracting its mean and dividing by its standard deviation, makes the contours nearly circular and the path almost straight.

• Features on different scales
• Stretched, narrow valley
• Zig-zag descent

One common mistake: computing the mean and standard deviation on the whole dataset before splitting. That leaks information from the test set. Fit the scaler on training data only, then apply those same numbers to validation and test data. Note the normal equation itself does not care about scale, but gradient descent very much does.

• Fit scaler on training data only
• Apply same transform to test
• Normal equation does not need it

What if the relationship is curved? Add columns for x squared, x cubed, and so on. The model draws a curve, yet it is still linear regression, because it is linear in the parameters. We simply changed the design matrix, and the same normal equation or gradient descent works unchanged. Scaling now matters even more, since x cubed is huge.

• Add x squared, x cubed as features
• Still linear in the parameters

Remember the bias variance trade off. A degree one line underfits. A cubic captures the shape. A degree fifteen polynomial threads through every training point and swings wildly in between, memorizing noise. Training error keeps falling with degree, but test error forms a U. Choose complexity using validation data, never training error.

• Degree 1: underfit
• Degree 3: good
• Degree 15: overfit

Linear regression leans on assumptions. The relationship should be linear in the features. Errors should be independent of each other. The spread of errors should be constant, called homoscedasticity. And for confidence intervals and tests, errors should be roughly normal. Predictions can still work when these bend, but inference becomes unreliable.

• Linearity
• Independent errors
• Constant variance
• Roughly normal errors

The best diagnostic is the residual plot, residuals on the vertical axis against predictions on the horizontal. A healthy plot is a shapeless cloud around zero. A curve means you missed nonlinearity, so add polynomial terms. A funnel means variance changes with the prediction, so consider transforming the target, for example with a logarithm.

• Residuals vs predictions
• Random cloud = good
• Curve or funnel = trouble

To check normality, plot a histogram of the residuals, or better, a Q Q plot. It compares residual quantiles to those of a perfect normal distribution. Points hugging the diagonal line mean a good match. Heavy tails curving away at both ends warn of outliers that could be distorting your fit.

• Histogram of residuals
• Q-Q plot near diagonal line

When two features carry almost the same information, say height in centimeters and height in inches, X transpose X becomes nearly singular. The valley turns into a flat trough, and the weights can swing to huge opposite values with tiny changes in data. Predictions may remain fine, but individual coefficients become meaningless. Regularization, coming in episode seventeen, cures this.

• Correlated features
• X^T X nearly singular
• Unstable, huge weights

Two more limits. Because errors are squared, a single extreme point can pull the whole line toward it, so inspect your data. And a line keeps going forever. Predicting far outside the range you trained on, extrapolating, assumes the trend continues, which is often false. Also remember that a strong fit shows association, not causation.

• Squared loss amplifies outliers
• Do not predict far outside data range

So which should you use? With a modest number of features, the normal equation is exact and needs no tuning. With very many features or enormous data, gradient descent and its stochastic variants scale far better, at the price of choosing a learning rate and scaling features. Both minimize the same convex bowl.

• Normal: exact, cubic cost
• Gradient descent: scalable, needs alpha

Let us recap. Linear regression minimizes the total area of squared residuals. Setting the gradient to zero gives the normal equation, theta equals X transpose X inverse times X transpose y. The gradient of MSE, two over n times X transpose times errors, powers gradient descent. Scale your features, add polynomial terms carefully, and always inspect residuals.

• Residuals, squares, MSE, MAE, RMSE
• theta = (X^T X)^-1 X^T y
• Gradient: (2/n) X^T (X theta - y)
• Scale, check residuals

But what if we want to predict a category, like spam or not spam, instead of a number? A straight line fails badly there. Next episode we bend the line with the sigmoid function to get logistic regression, and meet softmax for multiple classes.

• Predict categories, not numbers
• Sigmoid squeezes output to 0-1