Multiple Linear Regression and Least Squares

Definition Multiple linear regression predicts a numerical target from $p$ features $x_1, \dots, x_p$ with

For $n$ training cases, collect the parameters in $\boldsymbol{\theta} = (b, w_1, \dots, w_p)$ and build the design matrix $X$ from a leading column of ones, which carries the intercept, followed by one column per feature. Then $X$ is $n \times (p + 1)$, $\boldsymbol{\theta}$ has $p + 1$ entries, and the vector of predictions $\hat{\mathbf{y}} = X\boldsymbol{\theta}$ has $n$ entries. Because the ones column already supplies the intercept, no separate intercept is added. Least squares chooses $\boldsymbol{\theta}$ to minimize $\lVert X\boldsymbol{\theta} - \mathbf{y} \rVert_2^2$, the sum of the squared residuals $r_i = \hat{y}_i - y_i$, which on this page are prediction minus observed as in Simple Linear Regression.

Definition For a real $n \times k$ matrix $X$ and a target vector $\mathbf{y}$ in $\mathbb{R}^n$, the normal equations are

They come from the objective $J(\boldsymbol{\theta}) = \tfrac12 \lVert X\boldsymbol{\theta} - \mathbf{y} \rVert_2^2$, whose gradient is

The function $J$ is convex, because its Hessian $X^{\top}X$ is positive semidefinite, so $\boldsymbol{\theta}$ minimizes $J$ exactly when the gradient is zero. The least-squares solutions are therefore exactly the solutions of the normal equations, and at least one always exists. Entry $j$ of $X^{\top}(X\boldsymbol{\theta} - \mathbf{y})$ is the dot product of column $j$ with the residual vector, so the equations say that the residual is orthogonal to every column of $X$. With an intercept column, the first equation says that the residuals add to zero.

Definition The column space $\operatorname{Col} X$ is the set of all vectors $X\boldsymbol{\theta}$, the prediction vectors the model can produce for the $n$ training cases. The least-squares fitted vector $\hat{\mathbf{y}} = X\hat{\boldsymbol{\theta}}$ is the orthogonal projection $\mathbf{p}$ of $\mathbf{y}$ onto $\operatorname{Col} X$. For every $\boldsymbol{\theta}$, the vector $X\boldsymbol{\theta} - \mathbf{p}$ lies in $\operatorname{Col} X$ and $\mathbf{p} - \mathbf{y}$ is orthogonal to it, so

Hence $\boldsymbol{\theta}$ is a least-squares solution exactly when $X\boldsymbol{\theta} = \mathbf{p}$, and the smallest sum of squared residuals is $\lVert \mathbf{p} - \mathbf{y} \rVert_2^2$. The fitted vector is unique even when several coefficient vectors produce it. This picture lives in observation space $\mathbb{R}^n$, with one axis per training case, not in feature space, where the fitted model is a plane over the features.

Definition An $n \times k$ matrix $X$ has full column rank when $\operatorname{rank} X = k$, that is, when its columns are linearly independent, so that $X\mathbf{v} = \mathbf{0}$ only for $\mathbf{v} = \mathbf{0}$. This requires $n \ge k$, so a design with an intercept and $p$ features needs at least $p + 1$ cases. The matrix $X^{\top}X$ is invertible if and only if $X$ has full column rank, because $X^{\top}X\mathbf{v} = \mathbf{0}$ gives $\lVert X\mathbf{v} \rVert_2^2 = \mathbf{v}^{\top}X^{\top}X\mathbf{v} = 0$ and hence $X\mathbf{v} = \mathbf{0}$. With full column rank the least-squares coefficients are unique,

Without it, if $\hat{\boldsymbol{\theta}}$ is one least-squares solution, the solutions are exactly the vectors $\hat{\boldsymbol{\theta}} + \mathbf{h}$ with $X\mathbf{h} = \mathbf{0}$. They all give the same fitted vector, so the predictions for the training cases stay unique while the coefficients do not. A feature that is constant across the cases, for example, repeats the ones column up to a factor and breaks full column rank.

Definition In a fitted multiple regression, the coefficient $w_j$ is the partial effect of $x_j$ in the model, the change in the predicted value when $x_j$ increases by one unit while every other feature in the model stays fixed. It is measured in target units per unit of $x_j$. It describes the fitted model, and it reflects a comparison the data support only when they contain cases in which $x_j$ varies while the other features stay roughly the same. It is not a causal effect unless the study design or a separate causal argument justifies that reading. If another feature is computed from $x_j$, such as $x_j^2$, the other features cannot stay fixed and $w_j$ alone has no such reading.

Definition Let $X$ be a real $n \times k$ matrix with full column rank. Its reduced QR decomposition is

where $Q$ is $n \times k$ with orthonormal columns, so that $Q^{\top}Q = I_k$, and $R$ is $k \times k$, upper triangular and invertible. Requiring positive diagonal entries in $R$ makes the decomposition unique. Substituting it into the normal equations gives $R^{\top}R\boldsymbol{\theta} = R^{\top}Q^{\top}\mathbf{y}$, and since $R^{\top}$ is invertible the least-squares coefficients solve

a triangular system solved by back substitution, while the fitted vector is $QQ^{\top}\mathbf{y}$. This route never forms $X^{\top}X$. For full column rank the condition number $\kappa_2(X)$, the ratio of the largest to the smallest singular value, satisfies $\kappa_2(X^{\top}X) = \kappa_2(X)^2$, so working through $X^{\top}X$ can lose up to about twice as many correct digits to rounding as working with $X$ directly. The singular value decomposition is another stable route. It also handles rank deficiency by treating singular values below a stated tolerance as zero, and it is the method behind NumPy’s lstsq. Neither method can create information that the columns do not contain.