How to reconcile these two matrix equations for obtaining the coefficients for a linear least squares fit?

Question

In ordinary least squares linear regression, given a set of data points $(x_1,y_1),(x_2,y_2),...(x_N,y_N)$, that we want to fit to the function $y=\beta_0 + \beta_1 x$, we would usually write the linear model as the matrix equation $$ \begin{align} \mathbf{Y} &= \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \\ \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{pmatrix} &= \begin{pmatrix} 1 & x_1 \\ 1 & x_2 \\ \vdots & \vdots \\ 1 & x_N \end{pmatrix} \cdot \begin{pmatrix} \beta_0 \\ \beta_1\end{pmatrix} + \begin{pmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_N \end{pmatrix} \end{align} \tag{1} $$ where $\beta_0$ and $\beta_1$ are the intercept and slope of the model, respectively, and $\varepsilon_k$ are the residuals between the model and each data point.

Minimizing the sum of the squares of these residuals (by differentiating with respect to $\boldsymbol{\beta}$ and setting the result equal to zero) then gives the equation $$ \mathbf{X}^\textrm{T} \mathbf{Y} = (\mathbf{X}^\textrm{T}\mathbf{X})\boldsymbol{\beta} \tag{2} $$ and so the typical equation to solve to find the best fit parameters is $$ \boldsymbol{\beta} = (\mathbf{X}^\textrm{T}\mathbf{X})^{-1}\mathbf{X}^\textrm{T}\mathbf{Y} \tag{3} $$ (for example, see this link).

This is fine for me so far.

However, I am confused that it seems like the exact same result is obtained by solving instead $$ \mathbf{Y} = \mathbf{X}\boldsymbol{\beta} \tag{4} $$ The following Matlab script illustrates this, where the outputs for beta1 and beta2 are identical:

x_data = [0.1 ; 3.6 ; 5.2 ; 8.1];
y_data = [0.15 ; 3.5  ; 5.4  ; 7.6];

X = [ones(length(x_data),1) x_data];
Y = y_data;

% First method
beta1 = (X.'*X) \ (X.'*Y);        % instead of backslash, can also do beta = inv(X.'*X)*(X.')*Y;

% Second method
beta2 = X \ Y;

Equation (2) arises from minimizing the residuals, whereas Eq. (4) seems to have just neglected the residuals altogether?

Is it as simple as the fact that the $\mathbf{X}^\textrm{T}$ can simply cancel on both sides of Eq. (2)? In which case, why do we see the form given in Eq. (2) used so much?

(It is obviously much more computationally expensive to do method 1, compared with method 2.)

What am I missing here?

Thanks

Answer 1

Under the condition that $X \in \mathbb{R}^{n \times p}$ is of full column rank and $Y = X\beta$ has a solution (which is implied by the consistency condition $\operatorname{rank}(\begin{bmatrix}X & Y \end{bmatrix}) = \operatorname{rank}(X)$), you do make a good point: the solution $\hat{\beta}$ to the linear system $Y = X\beta$ indeed coincides with the closed-form least-squares estimate expression $\hat{\beta} = (X^TX)^{-1}X^TY$, which is readily from Eq. (2). However, the reason is slightly more involved than what you proposed "Is it as simple as the fact that the $X^T$ can simply cancel on both sides of Eq. (2)?" Simply cancelling $X^T$ from both sides of Eq. (2) as treating it as a non-zero real number is not rigorous enough. You will need some linear algebra weapon here.

One argument goes as follows: another way of viewing the normal equation $X^TX\beta = X^TY$ is \begin{align} X^T(Y - X\beta) = 0, \tag{$*$} \end{align} hence if $\beta$ solves Eq. (2), then $Y - X\beta$ must fall in the null space of $X^T$, whose dimensionality is $p - \operatorname{rank}(X^T) = p - p = 0$ (this is where you need to use linear algebra theory). This implies that if $Y - X\beta$ solves $(*)$, then $Y - X\beta = 0$, i.e., $\beta$ must be a solution to $Y = X\beta$. The converse direction is trivial.

An alternative argument is by using QR decomposition of $X$ -- this is the default numerical recipe for R to get least-squares estimate. Suppose the QR decomposition of $X$ is $X = QR$, where $Q \in \mathbb{R}^{n \times p}$ is a column-orthogonal matrix (i.e., $Q^TQ = I_{(p)}$), and $R$ is an invertible (as $\operatorname{rank}(X) = p$) upper-triangular matrix. Eq. (2) can then be rewritten as $R^TQ^TQR\beta = R^TQ^TY$, which in turn by the orthogonality of $Q$ and invertibility of $R$ reduces to $R\beta = Q^TY$. On the other hand, in terms of $Q$ and $R$, Eq. (4) is $Y = QR\beta$, which becomes $Q^TY = R\beta$ after multiplying $Q^T$ on both sides. Therefore both Eq. (2) and Eq. (4) are essentially the equation $R\beta = Q^TY$, hence must yield the same numerical solution.

Finally, Eq. (3) is used so much because it gives a explicit expression of the least-squares estimate of $\beta$ in terms of the raw inputs/observations $X$ and $Y$, which is valuable in theory and perhaps in teaching/communication (of course, as you mentioned, in implementation, no sensible software would really try to invert the matrix $X^TX$). By comparison, it is less convenient to call the LSE as, say, "a solution to the linear system $Y = X\beta$". In addition, Eq. (2) helps to reveal (or is a natural consequence of) the nature of LSE: the LSE $\hat{\beta}$ aims to minimize the Euclidean distance $\|Y - X\beta\|^2 = (Y - X\beta)^T(Y - X\beta)$, and you can see that the term $X^TX$ naturally emerges from this objective function.