Linear Regression and the Least Squares Objective

A linear model predicts a scalar target $y \in \mathbb{R}$ from a feature vector $x \in \mathbb{R}^d$ as

where $w \in \mathbb{R}^d$ is the weight vector and $b \in \mathbb{R}$ is the bias (or intercept). It is common to absorb the bias into $w$ by appending a constant feature $x_0 = 1$ to every input, giving the compact form $\hat{y} = w^\top x$ with $w, x \in \mathbb{R}^{d+1}$. The model is called linear because $\hat{y}$ is a linear (affine) function of the parameters $w$ — not necessarily of the raw inputs, since features can include polynomials, logs, or interactions of the original variables. Given a training set $\{(x^{(i)}, y^{(i)})\}_{i=1}^n$, the natural way to judge a candidate $w$ is to look at the residuals $r_i = y^{(i)} - \hat{y}^{(i)} = y^{(i)} - w^\top x^{(i)}$. The least-squares objective (also called the sum of squared errors, SSE, or residual sum of squares, RSS) is

Dividing by $n$ gives the mean squared error (MSE), which is scale-independent across dataset sizes and is the version most often reported and optimized in practice:

Fitting the model means solving

Why squares, rather than absolute residuals or some other penalty? Two independent reasons converge on the same answer. Geometrically, minimizing the sum of squared residuals is equivalent to finding the orthogonal projection of the target vector $y \in \mathbb{R}^n$ onto the subspace spanned by the feature columns — the unique point in that subspace closest to $y$ in Euclidean distance, which makes squared error the natural companion to linear algebra. Statistically, if we assume $y^{(i)} = w^\top x^{(i)} + \varepsilon_i$ with independent noise $\varepsilon_i \sim \mathcal{N}(0, \sigma^2)$, the likelihood of the data is

Taking the log and dropping constants that don't depend on $w$, maximizing the log-likelihood is exactly the same as minimizing $\sum_i (y^{(i)} - w^\top x^{(i)})^2$. So least squares is the maximum-likelihood estimator under Gaussian noise — not an ad hoc choice but the loss implied by the simplest and most common noise assumption.

Fitting a linear regression model has three conceptual steps: represent the data as matrices, write down the loss as a function of $w$, and minimize it. The minimization itself can be done in closed form or iteratively; this article treats the setup and the shape of the objective, while the actual solvers are developed in their own right in The Normal Equations and Closed-Form Solution and Gradient Descent for Parameter Estimation. Stack the $n$ feature vectors (each with a leading 1 for the bias) into a design matrix $X \in \mathbb{R}^{n \times (d+1)}$, and stack the targets into $y \in \mathbb{R}^n$. The predictions for the whole dataset are $\hat{y} = Xw$, and the least-squares objective becomes

This is a quadratic function of $w$. Expanding it,

which is a sum of a constant, a linear term in $w$, and a quadratic term with matrix $X^\top X$. Because $X^\top X$ is positive semi-definite (it equals $A^\top A$ for some $A$, and $v^\top X^\top X v = \|Xv\|^2 \ge 0$ for all $v$), $J$ is a convex function of $w$: its graph is a bowl, possibly with a flat direction but never a saddle or a local bump. This convexity is the single most important structural fact about linear regression — it is what guarantees that "find the minimum" is a well-posed, solvable problem rather than a search over an unpredictable landscape. To find the minimizer, take the gradient of $J$ with respect to $w$ and set it to zero:

These are the normal equations; when $X^\top X$ is invertible, $w^\star = (X^\top X)^{-1} X^\top y$ solves the problem exactly in one step. When $X^\top X$ is singular or the data is very large, iterative methods like gradient descent are preferred instead. Both routes are explored in depth in their dedicated articles — the point here is that the objective itself, being a convex quadratic, is what makes both routes work. ``` ALGORITHM FitLinearRegressionBySSE(X, y) Input: X, an n x d matrix of feature rows (bias column optionally prepended) y, a length-n vector of targets Output: w, the least-squares weight vector 1. Augment X with a column of ones -> X' (n x (d+1)), unless already done 2. Define J(w) = (y - X' w)^T (y - X' w) # the objective to minimize 3. Solve for the stationary point of J: gradient(w) = -2 X'^T y + 2 X'^T X' w set gradient(w) = 0 ==> X'^T X' w = X'^T y # normal equations 4. If X'^T X' is invertible: w <- (X'^T X')^{-1} X'^T y # closed form, see Normal Equations article else: w <- solve via gradient descent or pseudo-inverse # see Gradient Descent article 5. return w ``` ``` ALGORITHM PredictAndEvaluate(w, X_new, y_new) Input: fitted weights w, new feature matrix X_new, true targets y_new (for evaluation) Output: predictions, and MSE 1. y_hat <- X_new w 2. residuals <- y_new - y_hat 3. MSE <- (1/n) sum(residuals[i]^2 for i in 1..n) 4. return y_hat, MSE ```

Objective construction. Building the design matrix $X$ costs $O(nd)$ to lay out; computing $X^\top X$ costs $O(nd^2)$ multiply-adds (each of the $d \times d$ entries is an inner product over $n$ rows), and $X^\top y$ costs $O(nd)$. Evaluating the loss $J(w)$ for a given $w$ costs $O(nd)$ — one pass over the data to compute residuals and square them.

Solving for $w$ depends on the chosen method:

• Closed form via the normal equations requires inverting (or factorizing) the $(d+1) \times (d+1)$ matrix $X^\top X$, which costs $O(d^3)$ using standard Gaussian elimination or Cholesky decomposition, on top of the $O(nd^2)$ to form $X^\top X$. Total: $O(nd^2 + d^3)$. This is efficient when $d$ is small to moderate (up to a few thousand) even for very large $n$, since the cost grows only linearly in $n$.
• Iterative solvers (gradient descent) cost $O(nd)$ per iteration and typically run for $T$ iterations, giving $O(Tnd)$ total — attractive when $d$ is huge or $X^\top X$ is ill-conditioned or too large to form explicitly.

Space. Storing $X$ costs $O(nd)$; storing $X^\top X$ costs $O(d^2)$, independent of $n$, which is why the normal-equation route scales gracefully to huge $n$ as long as $d$ stays modest.

Correctness — why minimizing $J$ is the right thing to do. The argument has two independent legs. Geometric leg. The set of vectors reachable as $Xw$ for some $w$ is exactly the column space of $X$, a subspace $\mathcal{C}(X) \subseteq \mathbb{R}^n$. Minimizing $\|y - Xw\|^2$ over all $w$ is minimizing distance from $y$ to a point in $\mathcal{C}(X)$. Basic Euclidean geometry says the closest point in a subspace to an external point is its orthogonal projection, and the residual vector $y - Xw^\star$ at that closest point must be perpendicular to every vector in the subspace — in particular perpendicular to every column of $X$:

which is precisely the normal equations derived from calculus above. The two derivations (orthogonality and vanishing gradient) agreeing is not a coincidence — the same minimum has both a geometric and a calculus certificate. Convexity leg. $J(w) = w^\top(X^\top X)w - 2 w^\top X^\top y + y^\top y$ has Hessian $\nabla^2_w J = 2 X^\top X$, which is positive semi-definite for any $X$ (shown above). A twice-differentiable function with PSD Hessian everywhere is convex, and for a convex function any stationary point (where the gradient vanishes) is a global minimum, not merely a local one. If $X^\top X$ is strictly positive definite (full column rank, i.e., $d+1 \le n$ and no exact linear dependencies among features), the minimum is additionally unique. If $X^\top X$ is singular (more features than independent data directions, or collinear features), infinitely many $w$ achieve the same minimal loss, and regularization — see Ridge Regression and L2 Regularization — is used to pick out one well-behaved solution.

Honest tradeoffs. Squared error weights large residuals quadratically, so a single wild outlier can dominate the fit far more than it should — the least-squares line can be dragged noticeably toward one bad point. Squared error also implicitly assumes a symmetric, Gaussian-shaped noise process; if the true noise is heavy-tailed or skewed, the maximum-likelihood justification for least squares no longer holds, and alternatives like least-absolute-deviations (L1 loss) or robust regression are more appropriate. Linear regression is also, by construction, blind to nonlinear relationships in the raw features — it fits a linear function of whatever features you give it, so nonlinearity in the true relationship must be captured by feature engineering (polynomial terms, interactions) rather than by the model itself.

Consider predicting exam score $y$ from hours studied $x$, using four students, and fit a model $\hat{y} = w_1 x + w_0$ (one feature plus bias).

Build the design matrix (bias column first) and target vector:

Step 1 — form $X^\top X$ and $X^\top y$.

(Check: the (1,1) entry of $X^\top X$ is $n = 4$; the (1,2) entry is $\sum x_i = 1+2+3+5=11$; the (2,2) entry is $\sum x_i^2 = 1+4+9+25=39$. And $\sum y_i = 52+63+75+90=280$; $\sum x_i y_i = 52+126+225+450 = 853$ — let's recompute carefully: $1\cdot52=52$, $2\cdot63=126$, $3\cdot75=225$, $5\cdot90=450$; sum $=52+126+225+450=853$. So $X^\top y = \begin{bmatrix}280\\853\end{bmatrix}$.)

Step 2 — solve the normal equations $X^\top X\, w = X^\top y$:

Determinant of $X^\top X$: $4 \cdot 39 - 11 \cdot 11 = 156 - 121 = 35$. Inverse: $\dfrac{1}{35}\begin{bmatrix} 39 & -11 \\ -11 & 4 \end{bmatrix}$.

Compute each entry: $39 \times 280 = 10920$; $11 \times 853 = 9383$; so $w_0$ numerator $= 10920 - 9383 = 1537$, giving $w_0 = 1537/35 \approx 43.914$. $-11 \times 280 = -3080$; $4 \times 853 = 3412$; so $w_1$ numerator $= -3080 + 3412 = 332$, giving $w_1 = 332/35 \approx 9.486$. So the fitted model is approximately

Step 3 — check residuals and the loss.

Sum of squared residuals $J(w^\star) \approx 1.96 + 0.01 + 6.92 + 1.80 = 10.69$. Mean squared error $\approx 10.69/4 \approx 2.67$.

Sanity check via orthogonality. The correctness argument above claims $X^\top r = 0$ at the optimum, i.e., residuals sum to zero and are uncorrelated with $x$. Sum of residuals: $-1.40+0.11+2.63-1.34 = 0.00$ (up to rounding) — confirmed. This is not a coincidence; it is the geometric fact that the residual vector is orthogonal to the column space of $X$, and the all-ones column forces $\sum r_i = 0$ exactly at the true optimum.

The least-squares method predates modern machine learning by two centuries and was developed to solve an astronomical problem: reconciling noisy, overdetermined observations of planetary and cometary orbits.

Linear regression fit by least squares is the workhorse of quantitative fields far beyond machine learning. In economics and finance, it underlies hedonic pricing models, CAPM asset-pricing regressions, and demand forecasting. In the sciences, it is the default tool for quantifying how one measured quantity depends on another, and its assumptions (linearity, homoscedastic Gaussian noise) are taught as the baseline model against which more complex ones are judged. In engineering, it is used for calibration curves, sensor fusion, and system identification, where fast, interpretable, closed-form fits are valuable in real time. Inside machine learning systems, least-squares linear regression appears as a baseline model against which every fancier method is benchmarked, as the final layer of many regression networks (a linear readout on top of learned features), and as the core building block generalized by Ridge Regression and L2 Regularization (adding a penalty for numerical stability and to control overfitting) and by Logistic Regression and the Sigmoid (replacing the identity output with a sigmoid for classification). The two solution methods it sets up — direct linear-algebra solves via the Normal Equations and iterative first-order optimization via Gradient Descent — recur throughout the rest of supervised learning, from support vector machines to the training of deep neural networks, making this simple model an essential template rather than a historical curiosity.

Ask the tutor: If two features in my dataset are almost perfectly correlated, what happens to the least-squares weights, and how does ridge regression fix it?