5.1 Linear Regression, Design Matrices, and Least Squares: Linearity Refers to Parameters, Not Input Being Limited to Straight Lines

A fitted line is a starting point for predicting task duration from resource records. Adding categories, interactions, and numeric fields requires a design matrix that defines each column and reveals redundant representations. Otherwise, the model can diligently calculate several interchangeable sets of coefficients.

Linear models are well-suited as the first class of supervised learning because they're simple; more precisely, because they clearly expose the relationships between feature representation, loss function, optimization, and interpretability.

Learning Objectives

• Use a design matrix to formulate a multivariate linear model;
• Distinguish between "linearity in parameters" and "linearity in original inputs";
• Understand the least squares projection perspective, rank, and numerical solution methods;
• Separate prediction, association interpretation, and causal inference.

1. From a Single Line to a Design Matrix

For $n$ samples and $p$ features:

$$y=X\beta +\epsilon .$$

If the intercept is included in the all-ones column of $X$, then $X\in {\mathbb{R}}^{n\times (p+1)}$ and $\beta$ represents the coefficients.

A single row can contain:

1                         intercept
team_size                 original numeric value
log(resources)            nonlinear transformation
team_size²                polynomial term
team_size × mission_type interaction term
region_north              one-hot indicator

The model remains linear in $\beta$, making it a linear model; however, it can produce curved relationships and interactions with the original inputs. The shape of the decision boundary or response is determined by the feature mapping.

2. The Goal of Ordinary Least Squares

OLS minimizes the squared residual along the response axis:

$$\hat{\beta }=\arg \underset{\beta }{min}\Vert y-X\beta {\Vert }_2^2.$$

This is not the perpendicular distance from data points to a geometric line. If X also has measurement error, consider an errors-in-variables model. Total least squares is one option under specific error and scaling assumptions, rather than an automatic fix.

The squared loss amplifies large residuals, making the method sensitive to outliers. Choosing it may stem from a Gaussian conditional-noise model, computational convenience, or an approximation of business loss, never because "linear regression is default" should we overlook the cost of error.

3. The Normal Equation Does Not Equal Explicit Inversion

The first-order condition is:

$$X^{T}X\hat{\beta }=X^{T}y.$$

If $X$ has full column rank, the solution can be written as $(X^{T}X{)}^{-1}{X}^{T}y$; however, numerical computation should generally avoid explicitly constructing the inverse matrix. Methods such as QR decomposition, SVD, or specialized least-squares solvers are more numerically stable and can handle rank-deficient cases.

Complexity depends on $n$, $p$, sparsity, and the solver. For dense matrices with $n\gg p$, a common QR solve costs about $O(n{p}^2)$; forming and factoring normal equations costs about $O(n{p}^2+p^3)$ and squares the condition number. Sparse high-dimensional problems require a cost model based on sparse operations and iteration counts, rather than a blanket $O(n^3)$ claim.

4. Viewing Through the Lens of Projection

The fitted values $\hat{y}=X\hat{\beta }$ represent the projection of $y$ onto the column space of $X$. The residuals $r=y-\hat{y}$ are orthogonal to the column space:

$$X^{T}r=0.$$

With an intercept, residuals sum to zero at the exact unweighted OLS solution. This runnable Python 3 + NumPy example checks orthogonality; save it as ols_demo.py and run python3 ols_demo.py.

import numpy as np

x = np.array([0.0, 1.0, 2.0])
y = np.array([1.0, 3.0, 4.0])
X = np.column_stack([np.ones_like(x), x])
beta, _, rank, _ = np.linalg.lstsq(X, y, rcond=None)
residual = y - X @ beta
np.testing.assert_allclose(X.T @ residual, 0.0, atol=1e-12)
print(np.round(beta, 6), rank)

Expect [1.166667 1.5 ] 2: intercept, slope, and matrix rank. Spacing can vary slightly by NumPy version. The assertion checks residual orthogonality, without establishing accuracy on new data. See NumPy's lstsq reference for rank and minimum-norm behavior.

5. Rank Deficiency and Collinearity

If columns are linearly dependent, $X^{T}X$ is singular and thus not invertible, leading to non-unique coefficient estimates. For example, including an intercept along with all one-hot columns for a categorical variable creates a dummy-variable trap.

SVD or pseudoinverse methods select a minimum-norm solution. Standard OLS fitted values $X\hat{\beta }$ on the training matrix are a unique projection even when coefficients are nonunique. On new inputs that break the training column relationships, different minimizing coefficient vectors can predict different values. The approach to resolving this depends on the goal:

• Remove redundant encoding;
• Merge or reparameterize the variables;
• Apply ridge regularization;
• Collect more data that spans the correlated directions;
• Interpret only identifiable linear combinations.

High correlation does not imply that any column is useless, nor does it establish a causal direction.

6. Weighted least squares

If the sample has known differing precision or sampling weights, minimize:

$$\sum _i{w}_i(y_i-x_i^T\beta {)}^2.$$

Frequency weights, inverse variance weights, and sampling weights have distinct meanings and should not be interchanged. Weights also alter the target population; training, metrics, and standard errors must all be computed and interpreted according to their design.

7. Ridge is a Stabilizer for Pathological Problems

$${\hat{\beta }}_{ridge}=\arg \underset{\beta }{min}\Vert y-X\beta {\Vert }_2^2+\lambda \sum _{j=1}^p{\beta }_j^2.$$

Ridge adds constraints to correlated or weakly identifiable directions, reducing variance at the cost of increased bias. It typically does not penalize the intercept, and the relative impact of features under the same $\lambda$ is determined by their scale.

Regularization alters the estimand; ridge coefficients cannot be directly used to compute ordinary OLS p-values. Predictive performance and inferential validity must be treated separately.

8. Prediction Boundaries

The model extrapolates using its feature map: raw numeric features give linear extrapolation, while polynomial features can grow much faster. Training resource values of 10–100 do not stop the model returning a value at 10,000, even when that value violates physical constraints.

Check the following:

• Training and deployment ranges;
• Impossible negative predictions or proportions greater than 1;
• Interaction terms and units of measurement;
• Behavior of feature transformations at extreme values;
• Prediction intervals and grouped prediction errors.

Transforming the target variable, using generalized linear models (GLMs), or applying constraints may better align with the underlying data generation process.

9. Coefficients Are Not Causal Effects

Given other columns in the design matrix, ${\beta }_j$ describes a conditional association or predictive parameter. Omitting confounders, selection bias, measurement error, or post-treatment controls all undermine causal interpretation.

Volumes 12, Chapters 6 and 8 discuss regression diagnostics and causal identification. High out-of-sample $R^2$ in machine learning does not guarantee intervention effectiveness.

Common Misconceptions

• Linear models can only fit straight lines: Feature maps can include polynomials and interactions.
• OLS minimizes the perpendicular distance from points to the line: It actually minimizes the residual in the y-direction.
• Implementing the normal equation means you must compute the inverse: QR or SVD methods are more numerically stable.
• A significant or large coefficient means the feature is important: Scale, correlation, and causal structure all influence interpretability.

Exercise

1. Construct the design matrix using x, x², category, x×category.
2. Verify that for OLS with an intercept, $X^{T}r\approx 0$.
3. Generate a perfectly collinear column and compare QR/SVD, pseudoinverse, and ridge regression.
4. Evaluate model extrapolation beyond the training range and check against business domain boundaries.

Summary

Linear regression connects feature representation with coefficient estimation. Least squares is a projection problem under a well-defined loss function; its stability and reliability depend on numerical rank, scale, and interpretability of the target variables.

The next lesson replaces the response with a Bernoulli distribution: logistic regression does not directly model 0/1 outcomes, but instead models the log-odds of the event occurring.