Trend - Polynomial Regression - Theory

Preface

Goal: Goes further from linear regression to polynomial regression. The theoretical foundations, explains why the math works.

The statistical properties for polynomial regression is different especially the utilization of gram matrix. From just basic math, we need to know linear algebra to do matrix operation.

Warning: This article contain basic matrix operations, such as transpose and inverse. Please be prepared.

Cheatsheet

I combined different things into a single coherent diagram.

• Regression model
• Error and residual analysis
• Hypothesis testing
• R² and R²(adjusted)

You can see how the polynomial regression differs from linear regressions.

You can obtain both the SVG source here:

• 📎 github.com/…/math/trend/equflow/equflow-04.svg

Note that I only cover samples, not the population. Population has a only slight difference. We will get into this later, still in this article.

---

1: Equation for Polynomial

Polynomial Regression (not the least square method)

Instead of focusing on Covariance, the equation is started with gram matrix.

Big Picture

Most of the equation in the left side part has been covered in polynomial calculation article.

• Trend - Polynomial Interpolation

I don’t think that I need to refresh the equation over and over again. But let’s dexribe the equation notation in other forms.

Observed Data

Let’s say we have some pairs of points from data observations.

$$\begin{align*} && {x_i}_\text{ (observed)} &= [\ldots] \quad , \\ && {y_i}_\text{ (observed)} &= [\ldots] \\ \end{align*}$$

We are going to find the curve fittings, and get the prediction results: $$\begin{align*} && {\hat{y}_i}_\text{ (prediction)} &= [\ldots] \\ \end{align*}$$

Let’s find the solution for this using polynomial coefficients.

Fit Model

Polynomial Prediction Model

We have the generic theoritical model with this explicit form:

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

Most of the time I only consider these three curve fitting model for predictions:

1. Linear Fit (first order)
2. Quadratic Fit (second order)
3. Cubic Fit (third order)

$$\begin{array}{|l|l|l|} \hline \scriptstyle \text{1st order} & \scriptstyle \text{linear} & \hat{y}_i = \beta_0 + \beta_1{x_i} \rule{0pt}{2ex} \\ \hline \scriptstyle \text{2nd order} & \scriptstyle \text{quadratic} & \hat{y}_i = \beta_0 + \beta_1{x_i} + \beta_2{x_i}^2 \rule{0pt}{2ex} \\ \hline \scriptstyle \text{3rd order} & \scriptstyle \text{cubic} & \hat{y}_i = \beta_0 + \beta_1{x_i} + \beta_2{x_i}^2 + \beta_3{x_i}^3 \rule{0pt}{2ex} \\ \hline \end{array}$$

Where order is the degree of the polynomial. And β (beta) denotes the coefficients for specific polynomial degrees.

Notes: higher-order terms help model curvature, but also risk overfitting. This can be problematic in real world situations.

Matrix Form

From here we have the famous vandermonde matrix, in general it is looks like this.

$$\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 & \ddots & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^k \\ \end{bmatrix} \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \vdots \\ \beta_k \\ \end{bmatrix} = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \\ \end{bmatrix} \\ $$

Where

• n: Number of observations (rows in X)

For the three curve fittings, the actual matrix can be described as following:

$$\begin{array}{|l|l|l|} \hline \scriptstyle \text{linear} & \mathbf{X} = \begin{bmatrix} 1 & x_1 \\ 1 & x_2 \\ \vdots & \vdots \\ 1 & x_n \\ \end{bmatrix} & \mathbf{y} = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \\ \end{bmatrix} \\ \hline \scriptstyle \text{quadratic} & \mathbf{X} = \begin{bmatrix} 1 & x_1 & x_1^2 \\ 1 & x_2 & x_2^2 \\ \vdots & \vdots & \vdots \\ 1 & x_n & x_n^2 \\ \end{bmatrix} & \mathbf{y} = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \\ \end{bmatrix} \\ \hline \scriptstyle \text{cubic} & \mathbf{X} = \begin{bmatrix} 1 & x_1 & x_1^2 & x_1^3 \\ 1 & x_2 & x_2^2 & x_2^3 \\ \vdots & \vdots & \vdots & \vdots \\ 1 & x_n & x_n^2 & x_n^3 \\ \end{bmatrix} & \mathbf{y} = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \\ \end{bmatrix} \\ \hline \end{array}$$

Solution

This is the solution, where X is the predictor variables, and y is dependent variables.

$$\begin{aligned} && X \mathbf{\beta} &= \mathbf{y} \\ \Rightarrow && X^\top X \mathbf{\beta} &= X^\top \mathbf{y} \\ \Rightarrow && \mathbf{\beta} &= (X^\top X)^{-1} X^\top \mathbf{y} \end{aligned}$$

For each curve fittings, the coefficients solutions can be described as following:

$$\begin{array}{|l|l|} \hline \scriptstyle \text{linear} & \left(\mathbf{X}^T \mathbf{X}\right) \begin{bmatrix} \beta_0 \\ \beta_1 \\ \end{bmatrix} = \mathbf{X}^T \mathbf{y} \\ \hline \scriptstyle \text{quadratic} & \left(\mathbf{X}^T \mathbf{X}\right) \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \end{bmatrix} = \mathbf{X}^T \mathbf{y} \\ \hline \scriptstyle \text{cubic} & \left(\mathbf{X}^T \mathbf{X}\right) \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \beta_3 \\ \end{bmatrix} = \mathbf{X}^T \mathbf{y} \\ \hline \end{array}$$

Explicit Gram Matrix

Now we also the generalized form the gram matrix (Xᵗ.X).

$$X^\top X = \begin{bmatrix} n & \sum x_m & \sum x_m^2 & \cdots & \sum x_m^k \\ \sum x_m & \sum x_m^2 & \sum x_m^3 & \cdots & \sum x_m^{k+1} \\ \sum x_m^2 & \sum x_m^3 & \sum x_m^4 & \cdots & \sum x_m^{k+2} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \sum x_m^k & \sum x_m^{k+1} & \sum x_m^{k+2} & \cdots & \sum x_m^{2k} \end{bmatrix}$$

Where

• m: Index used in sums for (Xᵗ.X) entries

The matrix is always symmetric.

We are going use computation. But for clarity, we can show the actual gram matrix for each curve fittings can be described as following:

$$\begin{array}{|l|l|} \hline \scriptstyle \text{linear} & \begin{bmatrix} n & \sum x_i \\ \sum x_i & \sum x_i^2 \\ \end{bmatrix} \begin{bmatrix} a \\ b \\ \end{bmatrix} = \begin{bmatrix} \sum y_i \\ \sum x_i y_i \\ \end{bmatrix} \\ \hline \scriptstyle \text{quadratic} & \begin{bmatrix} n & \sum x_i & \sum x_i^2 \\ \sum x_i & \sum x_i^2 & \sum x_i^3 \\ \sum x_i^2 & \sum x_i^3 & \sum x_i^4 \\ \end{bmatrix} \begin{bmatrix} a \\ b \\ c \\ \end{bmatrix} = \begin{bmatrix} \sum y_i \\ \sum x_i y_i \\ \sum x_i^2 y_i \\ \end{bmatrix} \\ \hline \scriptstyle \text{cubic} & \begin{bmatrix} n & \sum x_i & \sum x_i^2 & \sum x_i^3 \\ \sum x_i & \sum x_i^2 & \sum x_i^3 & \sum x_i^4 \\ \sum x_i^2 & \sum x_i^3 & \sum x_i^4 & \sum x_i^5 \\ \sum x_i^3 & \sum x_i^4 & \sum x_i^5 & \sum x_i^6 \\ \end{bmatrix} \begin{bmatrix} a \\ b \\ c \\ d \\ \end{bmatrix} = \begin{bmatrix} \sum y_i \\ \sum x_i y_i \\ \sum x_i^2 y_i \\ \sum x_i^3 y_i \\ \end{bmatrix} \\ \hline \end{array}$$

Inverse Gram Matrix

To calculate each Coefficients, we need the inverse Gram Matrix (Xᵗ.X)ˉ¹ from the x observed values. There is no simple closed-form expression for inverse gram matrix (Xᵗ.X)ˉ¹. We will use Excel computation instead.

Fit Prediction

The Estimated Coefficients (Calculated from Data), can be described as below:

$$\hat{\beta} = (X^\top X)^{-1} X^\top \mathbf{y}$$

Predicted Values

By Model Type

Using the estimated coefficients β, we substitute into the explicit polynomial forms. For each observed value xᵢ, we compute predictions yᵢ for each model:

$$\begin{array}{|l|l|} \hline \scriptstyle \textbf{Model Type} & \textbf{Prediction Equation} \\ \hline \scriptstyle \text{Observed predictors} & \mathbf{x} = [x_1, x_2, \ldots, x_n]^\top \\ \hline \scriptstyle \text{Linear fit} & \hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i \\ \hline \scriptstyle \text{Quadratic fit} & \hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i + \hat{\beta}_2 x_i^2 \\ \hline \scriptstyle \text{Cubic fit} & \hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i + \hat{\beta}_2 x_i^2 + \hat{\beta}_3 x_i^3 \\ \hline \end{array}$$

Direct Comparison

Predicted vs Observed Values

The actual calculations use the same observed xᵢvalues, enabling direct comparison with observed yᵢ:

$$\begin{array}{|l|l|} \hline \scriptstyle \text{samples} & {x_i}_\text{ (observed)} = [\ldots] \\ \hline \scriptstyle \text{samples} & {y_i}_\text{ (observed)} = [\ldots] \\ \hline \scriptstyle \text{linear} & {\hat{y}_i}_\text{ (prediction)} = [\ldots] \\ \hline \scriptstyle \text{quadratic} & {\hat{y}_i}_\text{ (prediction)} = [\ldots] \\ \hline \scriptstyle \text{cubic} & {\hat{y}_i}_\text{ (prediction)} = [\ldots] \\ \hline \end{array}$$

The model comparison reveals how polynomial complexity affects fit quality.

$$\begin{array}{|l|l|l|} \hline \text{Sample} & x_i \text{ (observed)} & y_i \text{ (observed)} \\ \hline 1 & x_1 & y_1 \\ 2 & x_2 & y_2 \\ \vdots & \vdots & \vdots \\ n & x_n & y_n \\ \hline \end{array}$$ $$\begin{array}{|l|l|l|l|} \hline \text{Prediction} & \text{Linear} & \text{Quadratic} & \text{Cubic} \\ \hline \hat{y}_1 & \hat{\beta}_0 + \hat{\beta}_1 x_1 & \hat{\beta}_0 + \hat{\beta}_1 x_1 + \hat{\beta}_2 x_1^2 & \hat{\beta}_0 + \cdots + \hat{\beta}_3 x_1^3 \\ \hat{y}_2 & \hat{\beta}_0 + \hat{\beta}_1 x_2 & \hat{\beta}_0 + \hat{\beta}_1 x_2 + \hat{\beta}_2 x_2^2 & \hat{\beta}_0 + \cdots + \hat{\beta}_3 x_2^3 \\ \vdots & \vdots & \vdots & \vdots \\ \hat{y}_n & \hat{\beta}_0 + \hat{\beta}_1 x_n & \hat{\beta}_0 + \hat{\beta}_1 x_n + \hat{\beta}_2 x_n^2 & \hat{\beta}_0 + \cdots + \hat{\beta}_3 x_n^3 \\ \hline \end{array}$$

The interesting parts comes in the next section below.

---

2: Equation for Regression

Most of the equation in the top side part alos has been covered in polynomial calculation article.

• Trend - Properties - Cheatsheet

No Variance?

Where’s Waldo?

So where is the famous Standard Deviation? Pearson, and covariance? Actually these properties are only relevant for linear regression. For polynomial regression we will utilize other properties. No slope, no intercept. Dang! Unlearn process happened.

We are still using SST (total variance) calculate SSR, with mean inside the equation.

$$\sum (y_i - \bar{y})^2$$

We still have covariance, but not in the simple form as used in previous linear regression.

Degrees of Freedom

In polynomial regression, the degrees of freedom (df) for the error term (residuals) is calculated as:

$$\begin{align*} & df = n - p - 1 \\ & \scriptstyle \text{(degree of freedom)} \end{align*}$$

Where the equation itself can be interpreted as:

• n = Total number of observations (data points).
• p = Number of predictors (independent variables) in the model.
• -1 = Accounts for the intercept term (β0​).

Now we can rewrite for each model as:

$$\begin{array}{|l|l|l|} \hline \scriptstyle \text{linear} & p = 1: \text{only } x & df = n - 2 \\ \hline \scriptstyle \text{quadratic} & p = 2: x \text{ and } x^2 & df = n - 3 \\ \hline \scriptstyle \text{cubic} & p = 3: x, x^2, x^3 & df = n - 4 \\ \hline \end{array}$$

The use of the predictors is pretty clear.

Residual

The sum of squared errors (SSE) is calculated the same way for all models, but the residuals differ based on the model’s complexity:

$$\begin{align*} SSE &= \sum\limits_{i=1}^{n} (\epsilon_i)^2 \\ &= \sum (y_i - \hat{y}_i)^2 \\ \end{align*}$$

Where the predicted values of y vary by model:

$$\begin{array}{|l|l|} \hline \scriptstyle \text{linear} & \sum (y_i - \left(\beta_0 + \beta_1{x_i})\right)^2 \\ \hline \scriptstyle \text{quadratic} & \sum (y_i - \left(\beta_0 + \beta_1{x_i} + \beta_2{x_i}^2)\right)^2 \\ \hline \scriptstyle \text{cubic} & \sum (y_i - \left(\beta_0 + \beta_1{x_i} + \beta_2{x_i}^2 + \beta_3{x_i}^3)\right)^2 \\ \hline \end{array}$$

Don’t worry, we are not going to derive these equation. We’ll compute these directly in Excel using tabular calculations.

Coefficient of Determination

The R Square (R²)

In general, the R Square (R²) is still defined as follows.

$$R^2 = 1 - \frac{\sum (y_i - \hat{y}_i)^2}{\sum (y_i - \bar{y})^2}$$

Using tabular in Excel (or python), we can just put the residual calculation above into above equation.

From here we can calculate the adjusted R Square with following equation:

$$R^2 _\text{(adjusted)} = 1 - \frac{{(1 - R^2)(n - 1)}}{{n - p - 1}}$$

MSE

MSE stand for Mean Standard Errors. The mean word refers to the average of squared errors, adjusted for degrees of freedom.

$$MSE = \frac{SSE}{df} = \frac{\sum (y_i - \hat{y}_i)^2 }{n - p - 1}$$

And for further calculation we can have this helper.

$$RMSE = \sqrt{MSE}$$

Gram Matrix

To calculate Standard Error of each Coefficients, we need Gram Matrix (Xᵗ.X) from the x observed values.

From the matrix form:

$$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 & \ddots & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^k \\ \end{bmatrix}$$

Covariance Matrix of β

β = coefficient estimates

The theoretical covariance matrix is:

$$\text{Cov}(\hat{\mathbf{\beta}}) = \sigma^2 (X^\top X)^{-1}$$

where:

• σ² is the true error variance (unknown in practice)
• (Xᵗ.X)ˉ¹ is the inverse Gram matrix

In practice, we estimate σ² using the Mean Squared Error (MSE):

$$\text{Cov}(\hat{\mathbf{\beta}}) \approx \text{MSE} \cdot (X^\top X)^{-1}$$

To find inverse gram matrix (Xᵗ.X)ˉ¹, we have to rely on computation instead of finding generalization equation.

Variance of β

β = coefficient estimates

The true variance of the coefficient estimates is

$$\text{Var}(\hat{\beta}_j) = \sigma^2 \cdot \left(X^\top X\right)^{-1}_{jj}$$

Where

• σ² is the unknown error variance

In practice, we estimate σ² using the Mean Squared Error (MSE):

$$\text{Var}(\hat{\beta}_j) \approx \text{MSE} \cdot \left(X^\top X\right)^{-1}_{jj}$$

Diagonal Matrix

The diagonal of the inverse Gram matrix, the diags|(Xᵗ.X)ˉ¹|, represents the scaled variance of β (variance per unit of σ²):

$$\left(X^\top X\right)^{-1}_{jj} \approx \text{Var}_{scaled}(\hat{\beta}_j)\cdot$$

Or in other notation can be written as

$$\text{diag}|(X^\top X)^{-1}|_j$$

The same applied for diags|(Xᵗ.X)ˉ¹|ⱼ. There is no generalization either, so we have to rely on computation instead.

Standard Error of Coefficient

Standard deviation is the square root of the variance right? Then the standard error (SE) is indeed the square root of the variance of the coefficient estimates.

For any coefficient βⱼ​ in polynomial regression:

$$\begin{align*} SE(\hat{\beta}_j) &= \sqrt{\text{Var}(\hat{\beta}_j) } \\ &= \sqrt{MSE \cdot (X^\top X)^{-1}_{jj}} \end{align*}$$

We can break down for each case. The standard Errors by Polynomial Degree are:

$$\begin{array}{|l|l|} \hline \scriptstyle \text{linear} & SE(\hat{\beta}_0) , SE(\hat{\beta}_1) \\ \hline \scriptstyle \text{quadratic} & SE(\hat{\beta}_0) , SE(\hat{\beta}_1) , SE(\hat{\beta}_2) \\ \hline \scriptstyle \text{cubic} & SE(\hat{\beta}_0) , SE(\hat{\beta}_1) , SE(\hat{\beta}_2) , SE(\hat{\beta}_3) \\ \hline \end{array}$$

t-value

In Polynomial Regression

The t-value for each coefficient β​ tests the null hypothesis H₀: βj=0. It is calculated as:

$$t_{\text{value}, j} = \frac{\hat{\beta}_j}{SE(\hat{\beta}_j)}$$

Where:

• Numerator β​ = Estimated coefficient
• Denumerator SE(β) = Standard error of β​ (from previous section)

H₀: βj=0 asserts that the predictor xⱼ (e.g., x, x², etc.) has no effect on the response variable y. If H₀​ is true, the term xⱼ should be dropped from the model.

p-value

No need any gamma function complexity for practical calculation. For beginner like us, it is better to use excel built-in fomula instead.

Confidence Intervals

We can continue to other statistical properties beyond just coefficients, but I gues I need to stop before this article become to complex. Let’s keep this article to be statistics for beginner.

Now it is a good time to write down these equations to tabular excel.

---

What Lies Ahead 🤔?

Math is fun, right?

From theoretical foundations we are going shift to practical implementation in daily basis, using spreadsheet formula, python tools, and visualizations.

While the theory explains why the math works, the practice shows how to execute it with real tools. With the flow, from theory to tabular spreadsheets, then from we are going to build visualization using python.

Consider diving into the next step by exploring [ Trend - Polynomial Regression - Formula ].