Linear Regression via the Normal Equations

Fit a line to data by solving a small linear system, and understand why the QR route is numerically safer.

The least-squares problem

Given data pairs, fit y = m x + c (or a general linear model) by minimizing the sum of squared residuals ||X beta - y||^2. Setting the gradient to zero gives the normal equations X^T X beta = X^T y, a small square system to solve for the coefficients beta.

import numpy as np
x=np.array([0,1,2,3,4.]); y=np.array([1,3,7,8,12.])
X=np.vstack([x,np.ones_like(x)]).T   # columns: slope, intercept
beta=np.linalg.solve(X.T@X, X.T@y)
print('slope,intercept:',np.round(beta,3))
print('via lstsq:',np.round(np.linalg.lstsq(X,y,rcond=None)[0],3))

What it computes

The solution beta = (X^T X)^-1 X^T y projects y onto the column space of X - the fitted values are the closest point in that space to the data, and the residual is orthogonal to it. This geometric picture (orthogonal projection) is the heart of least squares.

The conditioning trap

Forming X^T X squares the condition number of X, so if the columns are nearly collinear the normal equations lose accuracy fast. The numerically sound route is QR: factor X = Q R and solve R beta = Q^T y by back-substitution, which never forms X^T X. For rank-deficient or very ill-conditioned problems, the SVD-based pseudoinverse is safest.

Regularization

Adding a penalty gives ridge regression: solve (X^T X + lambda I) beta = X^T y. This stabilizes near-collinear problems, shrinks coefficients, and reduces overfitting - a small change with large practical benefit whenever features are correlated or data are scarce.