Linear regression as a least-squares linear algebra problem (bridge to stats course) · Linear Algebra

§20.5 built least squares as a linear-algebra tool; the statistics course's regression chapter derived the same formulas from calculus and probability instead. This lesson closes the loop explicitly: linear regression is §20.5's least squares, with A built from the predictor data — no new mathematics, only a new name for the columns.

The design matrix

For n data points (x_i,y_i), simple linear regression fits \hat y=\beta_0+\beta_1x. Writing this for all n points as one matrix equation:

$\begin{pmatrix}y_1\\y_2\\\vdots\\y_n\end{pmatrix}\approx\begin{pmatrix}1&x_1\\1&x_2\\\vdots&\vdots\\1&x_n\end{pmatrix}\begin{pmatrix}\beta_0\\\beta_1\end{pmatrix}$

— exactly \vec y\approx A\vec\beta, with A (the design matrix) having a column of 1s (for the intercept) and a column of x-values. This system is essentially always inconsistent (n equations, 2 unknowns, real data never fits a line exactly) — precisely the setup §20.5 built least squares to handle.

Normal equations, translated

§20.5's A^TA\hat{\vec x}=A^T\vec b becomes, expanding A^TA directly:

$A^TA=\begin{pmatrix}n&\sum x_i\\\sum x_i&\sum x_i^2\end{pmatrix},\qquad A^T\vec y=\begin{pmatrix}\sum y_i\\\sum x_iy_i\end{pmatrix}$

Solving this 2\times2 system (by §16.3's explicit inverse formula) reproduces the statistics course's slope and intercept formulas exactly — \hat\beta_1=S_{xy}/S_{xx} and \hat\beta_0=\bar y-\hat\beta_1\bar x are just the closed-form solution of these particular normal equations, worked out by hand once and memorized, rather than solved generically each time.

Multiple regression is the same equation, wider A

Adding more predictors costs nothing structurally: for k predictors, A gains more columns (n\times(k+1), still one column of 1s plus one per predictor), and \hat{\vec\beta}=(A^TA)^{-1}A^T\vec y is identical code to the simple case — this is precisely why §20.5 built the theory in general matrix form instead of specializing to two variables from the start. Everything else transfers too: R^2=1-\|\vec r\|^2/\|\vec y-\bar y\mathbb1\|^2 (§20.5's residual, normalized), and \vec r\perp\operatorname{Col}(A) (§20.5's orthogonality condition) is exactly the statistics-course fact that "residuals are orthogonal to every predictor," now recognized as the defining property of least squares itself, not a separate coincidence.

Doing it in Python

import numpy as np

x = np.array([1., 2., 3., 4., 5.])
y = np.array([2.1, 3.9, 6.2, 7.8, 10.1])

A = np.column_stack([np.ones_like(x), x])
beta = np.linalg.solve(A.T @ A, A.T @ y)
print(f"intercept={beta[0]:.4f}, slope={beta[1]:.4f}")

# Multiple regression: add x^2 as a second predictor
A2 = np.column_stack([np.ones_like(x), x, x**2])
beta2 = np.linalg.solve(A2.T @ A2, A2.T @ y)
print(f"quadratic fit: {beta2[0]:.4f} + {beta2[1]:.4f}x + {beta2[2]:.4f}x^2")

intercept=0.0500, slope=1.9900
quadratic fit: 0.2000 + 1.8614x + 0.0214x^2

Confirming R^2 and residual orthogonality — the exact same computation as §20.5, just with regression-specific names:

import numpy as np

x = np.array([1., 2., 3., 4., 5.])
y = np.array([2.1, 3.9, 6.2, 7.8, 10.1])
A = np.column_stack([np.ones_like(x), x])

beta = np.linalg.solve(A.T @ A, A.T @ y)
fitted = A @ beta
residual = y - fitted

ss_res = (residual**2).sum()
ss_tot = ((y - y.mean())**2).sum()
r_squared = 1 - ss_res / ss_tot

print("R^2:", round(float(r_squared), 4))
print("residual . column1 (intercept):", round(float(residual @ A[:, 0]), 10) or 0.0)
print("residual . column2 (x):", round(float(residual @ A[:, 1]), 10) or 0.0)

R^2: 0.9973
residual . column1 (intercept): 0.0
residual . column2 (x): 0.0

Worked example

Fit \hat y=\beta_0+\beta_1x to (1,3),(2,5),(3,8) via the normal equations, matching this to the statistics course's formulas.

A=\begin{pmatrix}1&1\\1&2\\1&3\end{pmatrix}, \vec y=(3,5,8).

$A^TA=\begin{pmatrix}3&6\\6&14\end{pmatrix},\qquad A^T\vec y=\begin{pmatrix}16\\37\end{pmatrix}$

Solve \begin{pmatrix}3&6\\6&14\end{pmatrix}\hat\beta=\begin{pmatrix}16\\37\end{pmatrix}: R_2\leftarrow R_2-2R_1: \begin{pmatrix}3&6\\0&2\end{pmatrix}\hat\beta=\begin{pmatrix}16\\5\end{pmatrix}\Rightarrow\hat\beta_1=2.5; 3\hat\beta_0+6(2.5)=16\Rightarrow\hat\beta_0=1/3.

$\boxed{\hat y=\tfrac13+2.5x}$

Sanity check against \hat\beta_1=S_{xy}/S_{xx}: \bar x=2,\bar y=16/3. S_{xx}=(1-2)^2+0+(3-2)^2=2. S_{xy}=(-1)(3-16/3)+0+(1)(8-16/3)=(-1)(-7/3)+(1)(8/3)=7/3+8/3=5. \hat\beta_1=5/2=2.5 ✓ — matches exactly, confirming the general normal-equations route and the statistics course's specialized formula agree, as they must.

Your turn

1. Why is the system A\vec\beta=\vec y almost always inconsistent for real regression data?

2. What does "residuals orthogonal to \operatorname{Col}(A)" mean concretely for a regression with an intercept column?

3. True or false: adding more predictor columns to A can only increase (never decrease) R^2 on the training data.

Solutions

1. A has n rows (one per data point) but only k+1 columns (predictors plus intercept), with n\gg k+1 typically — a genuinely tall matrix (§17.5's terminology). §17.5 already established that a tall matrix's column space is a proper subspace of \mathbb{R}^n, so \vec y essentially never happens to land exactly inside it: real, noisy y-values almost surely aren't an exact linear combination of the predictor columns.

2. Since one column of A is \vec1 (the intercept column), \vec r\cdot\vec1=0 means the residuals sum to exactly zero — exactly the statistics course's stated fact "residuals always sum to zero, by construction," now identified as a direct consequence of §20.5's general orthogonality condition rather than a separate algebraic coincidence.

3. True. Adding a column to A enlarges \operatorname{Col}(A) (§17.4: more spanning vectors can only increase, never shrink, a span), so the closest point in the bigger column space to \vec y can only be at least as close as before (§20.5's projection is always the closest point available) — \|\vec r\| can only shrink or stay the same, so R^2 can only rise or stay flat. (This is exactly why statistics courses need a separate penalized R^2 or held-out validation to detect overfitting — R^2 on the training data alone can't ever flag it.)

Check yourself in code

Fit \hat y=\beta_0+\beta_1x to (0,1),(1,4),(2,4),(3,7) via the normal equations and report \beta_0,\beta_1 and R^2.

Print exactly this:

beta0=1.3, beta1=1.8
R^2: 0.9

import numpy as np

x = np.array([0., 1., 2., 3.])
y = np.array([1., 4., 4., 7.])
A = np.column_stack([np.ones_like(x), x])

beta = np.linalg.solve(A.T @ A, A.T @ y)
print(f"beta0={round(beta[0], 4)}, beta1={round(beta[1], 4)}")
# compute R^2 and print "R^2: <value rounded to 4dp>"

import numpy as np

x = np.array([0., 1., 2., 3.])
y = np.array([1., 4., 4., 7.])
A = np.column_stack([np.ones_like(x), x])

beta = np.linalg.solve(A.T @ A, A.T @ y)
print(f"beta0={round(beta[0], 4)}, beta1={round(beta[1], 4)}")

fitted = A @ beta
ss_res = ((y - fitted)**2).sum()
ss_tot = ((y - y.mean())**2).sum()
print("R^2:", round(float(1 - ss_res / ss_tot), 4))

Linear regression is §20.5's least-squares problem with A built from predictor columns (plus a 1s column for the intercept) — the normal equations, the closed-form slope/intercept formulas, "residuals sum to zero," and R^2 are all direct instances of general facts already proved: orthogonal projection onto \operatorname{Col}(A), residual orthogonality, and monotonic growth of the column space under more predictors.

Next, closing this course: solving systems of differential equations — Calculus IV's phase-plane systems (§13.8), revisited now that eigenvalues and eigenvectors are fully developed rather than introduced as a self-contained prerequisite note.

50. Linear regression as a least-squares linear algebra problem (bridge to stats course)