Polynomial Regression With Matrices

In polynomial regression, we model the response $y$ as a polynomial of degree $k$ in the predictor $x{:}$

$$y = \beta_0 + \beta_1 x + \beta_2 x^2 + \cdots + \beta_k x^k.$$

Although this model is nonlinear in $x,$ it is linear in the coefficients $\beta_j.$ So we can write it in the same matrix form as ordinary linear regression, $\mathbf{y} = X\boldsymbol{\beta},$ by stacking the powers of each $x_i$ into a design matrix:

$$X = \begin{bmatrix} 1 & x_1 & x_1^2 & \cdots & x_1^k \\ 1 & x_2 & x_2^2 & \cdots & x_2^k \\ \vdots & \vdots & \vdots & & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^k \end{bmatrix}, \qquad \mathbf{y} = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}.$$

The least-squares coefficient vector is then given by the same normal equation:

$$\boldsymbol{\beta} = (X^T X)^{-1} X^T \mathbf{y}.$$

The only thing that changes when we move from linear to polynomial regression is the design matrix.

Quick illustration. To fit a quadratic to the data $(1, 4), (2, 7), (3, 5),$ the design matrix and response are

$$X = \begin{bmatrix} 1 & 1 & 1 \\ 1 & 2 & 4 \\ 1 & 3 & 9 \end{bmatrix}, \qquad \mathbf{y} = \begin{bmatrix} 4 \\ 7 \\ 5 \end{bmatrix}.$$