Linear Regression and Ridge Regression

Linear regression predicts a real-valued label as a weighted sum of the features:

ŷ = w₀ + w₁x₁ + w₂x₂ + … + w_dx_d

w₀ is the intercept (bias) and w₁ … w_d are the weights (slopes). Training chooses the weights that minimise the sum of squared errors (SSE), also called the residual sum of squares:

SSE = Σᵢ (yᵢ − ŷᵢ)²

This is the method of least squares.

Simple linear regression (one feature)

With one feature, ŷ = w₀ + w₁x. Setting the partial derivatives of SSE to zero gives:

• w₁ = S_xy / S_xx, where S_xy = Σ (xᵢ − x̄)(yᵢ − ȳ) and S_xx = Σ (xᵢ − x̄)².
• w₀ = ȳ − w₁ x̄, so the fitted line always passes through the point (x̄, ȳ).

Worked example. Fit a line to the points (1, 2), (2, 4), (3, 5), (4, 4), (5, 5).

x̄ = 3 and ȳ = 4, so w₁ = 6 / 10 = 0.6 and w₀ = 4 − 0.6 × 3 = 2.2. The fitted line is ŷ = 2.2 + 0.6x.

How good is the fit? R²

• SST = Σ (yᵢ − ȳ)² is the total variation of y around its mean.
• SSE = Σ (yᵢ − ŷᵢ)² is the variation the line does not explain.
• R² = 1 − SSE / SST is the fraction of the variation explained by the model. R² = 1 is a perfect fit; R² = 0 means the line is no better than predicting ȳ.

For the example, the predictions are 2.8, 3.4, 4.0, 4.6, 5.2, so SSE = 0.64 + 0.36 + 1 + 0.36 + 0.04 = 2.4, SST = 4 + 0 + 1 + 0 + 1 = 6 and R² = 1 − 2.4/6 = 0.6. For simple linear regression, R² equals the square of the correlation coefficient r between x and y.

xs = [1, 2, 3, 4, 5]
ys = [2, 4, 5, 4, 5]
n = len(xs)
mx, my = sum(xs) / n, sum(ys) / n
sxy = sum((x - mx) * (y - my) for x, y in zip(xs, ys))
sxx = sum((x - mx) ** 2 for x in xs)
syy = sum((y - my) ** 2 for y in ys)
w1 = sxy / sxx
w0 = my - w1 * mx
pred = [w0 + w1 * x for x in xs]
sse = sum((y - p) ** 2 for y, p in zip(ys, pred))
print(f"w1 = {w1:.2f}, w0 = {w0:.2f}")
print("predictions:", [round(p, 2) for p in pred])
print(f"SSE = {sse:.2f}, SST = {syy:.2f}, R^2 = {1 - sse / syy:.2f}")
print(f"r = {sxy / (sxx * syy) ** 0.5:.4f}, r^2 = {sxy ** 2 / (sxx * syy):.2f}")
print(f"prediction at x = 6: {w0 + w1 * 6:.2f}")

w1 = 0.60, w0 = 2.20
predictions: [2.8, 3.4, 4.0, 4.6, 5.2]
SSE = 2.40, SST = 6.00, R^2 = 0.60
r = 0.7746, r^2 = 0.60
prediction at x = 6: 5.80

Multiple linear regression and the normal equations

With d features, put each sample as a row of a matrix X with a leading column of 1s (for the intercept), and the labels in a vector y. The model is ŷ = Xw and

SSE(w) = ‖y − Xw‖².

Setting the gradient −2Xᵀ(y − Xw) to zero gives the normal equations:

XᵀX w = Xᵀy, so w = (XᵀX)⁻¹ Xᵀy when XᵀX is invertible.

XᵀX is invertible exactly when the columns of X are linearly independent. If one feature is a copy (or a linear combination) of others, called multicollinearity, XᵀX is singular and the least-squares weights are not unique.

Geometrically, Xw is the projection of y onto the column space of X, and the residual y − Xw is orthogonal to every column of X. That is why the residuals of a model with an intercept always sum to zero.

from fractions import Fraction as F

# Two features; the labels follow y = 1 + 2*x1 + 3*x2 exactly.
X = [[1, 0, 1], [1, 1, 0], [1, 1, 1], [1, 2, 1], [1, 0, 2]]
y = [1 + 2 * r[1] + 3 * r[2] for r in X]

def matmul(A, B):
    return [[sum(A[i][k] * B[k][j] for k in range(len(B))) for j in range(len(B[0]))] for i in range(len(A))]

def solve(A, b):
    # Gauss-Jordan elimination with exact fractions.
    n = len(A)
    M = [[F(v) for v in row] + [F(bv)] for row, bv in zip(A, b)]
    for c in range(n):
        p = next(r for r in range(c, n) if M[r][c] != 0)
        M[c], M[p] = M[p], M[c]
        M[c] = [v / M[c][c] for v in M[c]]
        for r in range(n):
            if r != c:
                M[r] = [a - M[r][c] * b for a, b in zip(M[r], M[c])]
    return [row[-1] for row in M]

Xt = [list(col) for col in zip(*X)]
XtX = matmul(Xt, X)
Xty = [sum(Xt[i][k] * y[k] for k in range(len(y))) for i in range(3)]
print("X^T X =", XtX)
print("X^T y =", Xty)
print("w =", [str(v) for v in solve(XtX, Xty)])

X^T X = [[5, 4, 5], [4, 6, 3], [5, 3, 7]]
X^T y = [28, 25, 32]
w = ['1', '2', '3']

The normal equations recover the true weights (1, 2, 3) exactly because the labels have no noise.

Gradient descent

For large data sets, inverting XᵀX (cost about d³) is replaced by gradient descent: start with any w and repeat

w ← w − η ∇SSE(w), where ∇ of the mean squared error is −(2/n) Xᵀ(y − Xw).

η is the learning rate. Too small and convergence is slow; too large and the updates overshoot and diverge. The squared-error loss of linear regression is convex, so gradient descent with a suitable η reaches the global minimum.

xs = [1, 2, 3, 4, 5]
ys = [2, 4, 5, 4, 5]
n = len(xs)
w0, w1, eta = 0.0, 0.0, 0.05
for step in range(1, 2001):
    g0 = -2 / n * sum(y - (w0 + w1 * x) for x, y in zip(xs, ys))
    g1 = -2 / n * sum((y - (w0 + w1 * x)) * x for x, y in zip(xs, ys))
    w0, w1 = w0 - eta * g0, w1 - eta * g1
    if step in (1, 10, 100, 2000):
        print(f"step {step:4}: w0 = {w0:.4f}, w1 = {w1:.4f}")

step    1: w0 = 0.4000, w1 = 1.3200
step   10: w0 = 0.6072, w1 = 1.0412
step  100: w0 = 1.8566, w1 = 0.6951
step 2000: w0 = 2.2000, w1 = 0.6000

Ridge regression (L2 regularization)

When there are many features, or features are nearly collinear, least squares can give huge weights that fit noise. Ridge regression adds a penalty on the size of the weights:

minimise Σ (yᵢ − ŷᵢ)² + λ Σⱼ wⱼ² (the intercept w₀ is usually not penalised)

• λ ≥ 0 is the regularization strength, a hyperparameter chosen by validation.
• λ = 0 gives ordinary least squares; as λ → ∞ the weights shrink towards 0 (but never become exactly 0).
• The solution is w = (XᵀX + λI)⁻¹ Xᵀy. Adding λI makes the matrix invertible for any λ > 0, so ridge regression always has a unique solution, even with collinear features.
• Larger λ means higher bias, lower variance.

For one feature with centred data, the ridge slope is w₁ = S_xy / (S_xx + λ). In the worked example, λ = 2 gives w₁ = 6 / (10 + 2) = 0.5 instead of 0.6: the slope shrinks towards zero.

(Lasso uses the penalty λ Σ |wⱼ| instead, which can set weights exactly to zero. It is not in the GATE DA syllabus by name, but the contrast with ridge is a common question.)

xs = [1, 2, 3, 4, 5]
ys = [2, 4, 5, 4, 5]
mx, my = sum(xs) / 5, sum(ys) / 5
sxy = sum((x - mx) * (y - my) for x, y in zip(xs, ys))
sxx = sum((x - mx) ** 2 for x in xs)
for lam in [0, 2, 10, 100]:
    w1 = sxy / (sxx + lam)
    print(f"lambda = {lam:3}: w1 = {w1:.4f}, w0 = {my - w1 * mx:.4f}")

lambda =   0: w1 = 0.6000, w0 = 2.2000
lambda =   2: w1 = 0.5000, w0 = 2.5000
lambda =  10: w1 = 0.3000, w0 = 3.1000
lambda = 100: w1 = 0.0545, w0 = 3.8364

Practice questions