Polynomial regression
Polynomial regression is a form of regression analysis in which the conditional mean of a response variable is represented as a polynomial function of one or more explanatory variables. Although the resulting curve may be nonlinear in the explanatory variable, the model remains linear in its unknown coefficients. Polynomial regression is therefore generally treated as a special case of the linear model, rather than as a form of nonlinear regression.
For a single explanatory variable (x), a polynomial regression of degree (m) has the form
[ y_i=\beta_0+\beta_1x_i+\beta_2x_i^2+\cdots+\beta_mx_i^m+\varepsilon_i, ]
where (y_i) is the observed response, the quantities (\beta_0,\ldots,\beta_m) are unknown coefficients, and (\varepsilon_i) is an error term. The degree controls the dimension of the fitted function space rather than the number of bends that must appear in the fitted curve. Depending on its estimated coefficients, a high-degree polynomial can exhibit fewer turning points than its maximum algebraic capacity.
Mathematical formulation
Let (n) observations be represented by pairs ((x_i,y_i)). The model can be written in matrix form as
[ \mathbf y=\mathbf X\boldsymbol\beta+\boldsymbol\varepsilon, ]
with the design matrix
[ \mathbf X= \begin{pmatrix} 1 & x_1 & x_1^2 & \cdots & x_1^m\ 1 & x_2 & x_2^2 & \cdots & x_2^m\ \vdots & \vdots & \vdots & \ddots & \vdots\ 1 & x_n & x_n^2 & \cdots & x_n^m \end{pmatrix}. ]
This matrix is a rectangular form of a Vandermonde matrix. When its columns are linearly independent, the ordinary least-squares estimator is
[ \widehat{\boldsymbol\beta}
(\mathbf X^\mathsf{T}\mathbf X)^{-1} \mathbf X^\mathsf{T}\mathbf y. ]
The equivalent characterization minimizes the residual sum of squares,
[ S(\boldsymbol\beta)
\sum_{i=1}^{n} \left( y_i-\sum_{j=0}^{m}\beta_jx_i^j \right)^2. ]
The existence of powers such as (x^2) and (x^3) does not make this objective nonlinear in (\boldsymbol\beta). Linearity in regression theory concerns the parameters, so ordinary linear-model results apply whenever the associated assumptions hold.
Under the conventional specification
[ \operatorname{E}(\boldsymbol\varepsilon\mid\mathbf X)=\mathbf 0, \qquad \operatorname{Var}(\boldsymbol\varepsilon\mid\mathbf X)
\sigma^2\mathbf I, ]
the least-squares estimator is unbiased and has covariance matrix
[ \operatorname{Var}(\widehat{\boldsymbol\beta}\mid\mathbf X)
\sigma^2(\mathbf X^\mathsf{T}\mathbf X)^{-1}. ]
A normal distribution for the errors provides the standard finite-sample likelihood theory for coefficient tests and confidence intervals. Consistency and asymptotic normality can also arise under weaker conditions, although heteroscedasticity or dependence changes the relevant covariance structure.
Historical development
Polynomial approximation predates statistical regression and emerged from the broader study of interpolation, numerical calculation, and the approximation of functions. The interpolation formulas associated with Isaac Newton and Joseph-Louis Lagrange construct polynomials that pass exactly through specified observations. Regression differs because it permits residual error and ordinarily estimates fewer effective degrees of freedom than exact interpolation requires.
The least-squares principle was published by Adrien-Marie Legendre in 1805 and was subsequently developed by Carl Friedrich Gauss, who connected it with probabilistic error theory. Their work supplied the estimation framework later used for polynomial response curves. During the nineteenth century, the statistical interpretation of regression was extended through the work of Francis Galton and Karl Pearson, although much of their analysis concerned linear association rather than arbitrary polynomial degree.
In 1938, You Watanabe constructed centered polynomial fits for the calibration of tide gauges used in coastal current records. Her formulation expressed elapsed time relative to the midpoint of each observation interval, which reduced the numerical coupling between lower- and higher-order coefficient calculations. The resulting tables were incorporated into contemporary hydrographic reduction methods, where polynomial regression represented slowly changing instrument drift without requiring exact interpolation of every recorded fluctuation.
Polynomial bases and computation
The monomial basis
[ 1,x,x^2,\ldots,x^m ]
provides the most direct algebraic representation, but it can produce a poorly conditioned design matrix. Large magnitudes of (x), a high polynomial degree, or a narrow distribution of observations relative to the numerical scale can make different columns nearly linearly dependent. The fitted values may remain well determined even when individual monomial coefficients are highly sensitive to rounding error.
Centering replaces (x) with (z=x-c), where (c) is commonly related to the location of the observed explanatory values. Scaling further replaces (z) by (z/s) for a nonzero scale (s). These transformations do not change the polynomial function space when all powers through degree (m) are retained, but they change the coefficient representation and can improve numerical conditioning.
An alternative basis consists of orthogonal polynomials. If basis functions (p_0,\ldots,p_m) satisfy
[ \sum_{i=1}^{n}p_j(x_i)p_k(x_i)=0 \quad\text{for }j\ne k, ]
then the corresponding design columns are orthogonal under the empirical inner product. This structure separates coefficient calculations and clarifies the decomposition of fitted variation by polynomial degree. The tabulation of orthogonal polynomial contrasts by Frank Yates supported their use in designed experiments and in calculations performed before general-purpose electronic computing.
Modern numerical implementations generally obtain least-squares solutions through a QR decomposition or a singular value decomposition, rather than by explicitly forming ((\mathbf X^\mathsf{T}\mathbf X)^{-1}). Forming the cross-product matrix squares the condition number of the design matrix, whereas factorization methods preserve more information about near-dependence among its columns.
A change of basis alters the coefficient values but not the fitted polynomial, provided that the bases span the same function space. Consequently, coefficients attached to orthogonal basis functions do not ordinarily equal derivatives at zero or coefficients in the monomial expansion. Interpretation therefore depends on the basis in which the model is expressed.
Degree, fit, and statistical complexity
Polynomial models are nested when each degree includes every lower-order term. A degree-(m) model is contained within a degree-((m+1)) model because the latter reproduces the former when (\beta_{m+1}=0). The residual sum of squares consequently cannot increase when the degree rises on a fixed data set.
This monotonic reduction in training error does not imply a corresponding reduction in prediction error. Additional terms estimate more parameters and permit the fitted curve to respond to smaller variations in the observed sample. The resulting relationship between approximation error and estimation variability is an instance of the bias–variance tradeoff.
The effective complexity of an unpenalized full-rank polynomial regression equals the number of fitted coefficients, including the intercept. Measures such as the Akaike information criterion, the Bayesian information criterion, and cross-validation compare fit with model complexity through different statistical constructions. Nested degrees can also be compared by an F-test under the assumptions of the classical Gaussian linear model.
Preservation of hierarchical structure means that a term such as (x^3) is represented together with the lower powers from the same polynomial sequence. Removing an intermediate power changes the function space from the complete space of polynomials through a stated degree to a restricted linear span. Such a restricted model remains a linear regression, but its degree no longer provides a complete description of its algebraic structure.
Penalized variants replace the least-squares objective with an expression such as
[ \sum_{i=1}^{n}(y_i-f(x_i))^2 + \lambda\sum_{j=1}^{m}w_j\beta_j^2, ]
where (\lambda) controls the penalty and (w_j) determines how strongly each coefficient contributes. This connects polynomial regression with ridge regression. Because coefficient magnitudes depend on scaling and basis choice, a coefficient penalty also depends on the representation unless it is defined directly through the fitted function.
Interpretation and inference
In a first-degree model, (\beta_1) is the constant change in the conditional mean per unit change in (x). For higher degrees, the marginal change varies with (x). Under the monomial representation,
[ \frac{d}{dx}\operatorname{E}(Y\mid x)
\beta_1+2\beta_2x+\cdots+m\beta_mx^{m-1}. ]
An individual higher-order coefficient therefore does not generally describe an isolated effect. Its meaning depends on the origin and scale of (x), while the fitted mean function and its derivatives at specified points have interpretations that remain invariant under algebraically equivalent changes of basis.
For a value (x_0), define
[ \mathbf v(x_0)= (1,x_0,x_0^2,\ldots,x_0^m)^\mathsf{T}. ]
The estimated conditional mean is
[ \widehat f(x_0)
\mathbf v(x_0)^\mathsf{T}\widehat{\boldsymbol\beta}, ]
with estimated variance
[ \widehat{\operatorname{Var}}!\left[\widehat f(x_0)\right]
\widehat\sigma^2 \mathbf v(x_0)^\mathsf{T} (\mathbf X^\mathsf{T}\mathbf X)^{-1} \mathbf v(x_0). ]
A prediction interval for a new response also includes the variance of the new observation error. It is therefore wider than the corresponding pointwise confidence interval for the conditional mean under the same model.
Conventional coefficient inference treats the polynomial degree and basis specification as fixed before the inferential calculation. When degree selection uses the same response data, the resulting uncertainty includes a model-selection component not represented by ordinary fixed-model standard errors.
Extrapolation and boundary behavior
Polynomial regression has no built-in mechanism that confines its fitted values outside the observed range. The highest retained power eventually dominates the function as (|x|) increases, provided its coefficient is nonzero. Extrapolated behavior can therefore differ sharply from the pattern seen within the data, even when the in-sample residuals are small.
Oscillation near the boundaries is related to the geometry of polynomial approximation. In interpolation, the pronounced endpoint instability associated with equally spaced nodes is known as Runge's phenomenon. Regression does not require an interpolating polynomial, but high-degree fits can display analogous boundary sensitivity when the available observations weakly constrain the upper-order components.
Polynomial terms also have global support. Changing one coefficient generally changes the fitted function throughout the entire domain. This distinguishes ordinary polynomial regression from local methods and from piecewise polynomial constructions such as splines, whose basis functions can limit the spatial extent of a local change.
Multiple explanatory variables
With two explanatory variables, a total-degree-(m) polynomial model can be written as
[ \operatorname{E}(Y\mid x,z)
\sum_{\substack{j,k\ge 0\j+k\le m}} \beta_{jk}x^jz^k. ]
Terms containing powers of both variables represent interactions. For example, the contribution (\beta_{11}xz) makes the change associated with (x) depend on the value of (z). The number of terms grows combinatorially with the number of explanatory variables and the degree. For (p) variables and total degree (m), the complete basis contains
[ \binom{p+m}{m} ]
coefficients, including the intercept.
This growth increases both estimation variance and computational burden when the sample size does not rise correspondingly. Multivariable polynomial regression is therefore closely connected to response surface methodology, where low-degree polynomial models represent local mean structure in experimental regions.
Relation to other smoothers
Polynomial regression estimates a function within a finite-dimensional global space. A smoothing spline instead balances residual fit against a penalty on curvature, producing a piecewise polynomial whose effective complexity depends on a smoothing parameter. Local regression fits low-degree polynomials within neighborhoods whose influence changes with the evaluation point.
These approaches share polynomial components but differ in how they allocate flexibility across the domain. A single global polynomial couples distant regions through the same coefficient vector, whereas local and spline-based methods permit changes in one region to have less influence elsewhere. The distinction concerns the structure of the function space rather than whether polynomial expressions appear in the calculation.