Polynomial & Basis-Function Regression
Why fitting curves is still "linear" regression once you expand the basis.
On this page
Beginner: "linear regression" sounds like it should only ever produce a straight line, so a curvy fitted line showing up in a "linear model" can feel like a contradiction. It isn't one. Linear regression doesn't mean the fitted curve has to be straight — it means the model is a straight-line (linear) function of its parameters. If you take your raw input x and, before fitting anything, build a whole set of new features out of it — x, x², x³, ... for polynomials, or sin(x), cos(x) for oscillating signals, or a handful of "bump" shapes centered at different points — and then run ordinary linear regression on those transformed features, the fitted curve can bend and wiggle as much as the transformations allow, while the fitting procedure itself is completely unchanged plain linear regression underneath.
Intermediate: this trick has a name — basis expansion. Pick a fixed set of functions φ₁, φ₂, ..., φ_K (called basis functions), each one a specific, computable transformation of x, and model the target as a weighted sum of them: ŷ = θ₀ + θ₁φ₁(x) + θ₂φ₂(x) + ⋯ + θ_Kφ_K(x). Every basis function is evaluated the same way for every data point and never changes during fitting — only the weights θ_k are learned. That is exactly the structure ordinary least squares was built for: a linear combination of columns of a design matrix. Swap the raw feature columns for basis-function columns and every tool from ordinary linear regression (2.3.1) — the normal equations, the Gauss-Markov guarantee, and the ridge/lasso penalties from earlier in this module — applies to the result completely unchanged, because none of that machinery ever looks at what a feature column means, only that it's a fixed number attached to each row.
Advanced: the general recipe is "choose the features, then run linear regression," and different choices of basis specialize to different familiar techniques. Polynomial powers φ_k(x) = x^k give ordinary polynomial regression. Gaussian "bumps" centered at fixed points, φ_k(x) = exp(-(x-c_k)²/2s²), give radial-basis-function (RBF) regression — the diagram below builds and animates exactly this case, deliberately choosing a different basis than the polynomial one Module 1's overfitting lesson (2.1.4) already used, to make the point that basis expansion is a general idea, not a synonym for "add polynomial terms." Sine/cosine pairs at a few frequencies give Fourier-feature regression, well suited to periodic signals. Push the idea to its limit — an infinite-dimensional basis, computed only implicitly through a kernel function rather than ever building the columns explicitly — and you arrive at kernel methods, previewed briefly in the real-world examples below and developed properly in a later Kernel Methods module.
Worked intuition — why this doesn't add a new fitting algorithm. Suppose you have ten (x, y) points that clearly follow a gentle curve, not a line. You could reach for some exotic nonlinear solver — or you could build a new table with columns 1, x, x², x³ from those same ten x values, hand that table and y to the exact same normal-equations solver used for plain OLS, and get back four numbers θ₀, θ₁, θ₂, θ₃. Plot θ₀+θ₁x+θ₂x²+θ₃x³ and it curves through your ten points. Nothing about the solver changed; only the table handed to it did. That is the entire idea this lesson formalizes — feature engineering, not a new model.
A general basis expansion: a fixed set of basis functions φ₀, ..., φ_K (φ₀ ≡ 1 supplies the intercept), combined with learned weights θ_k. The model is linear in θ for any choice of the φ's, however nonlinear each φ is in x. Named special cases used throughout this lesson and elsewhere:
Polynomial regression fixes c_k at nothing (powers of x instead of bumps); RBF regression fixes a set of centers c_k and a bandwidth s; Fourier regression fixes a fundamental frequency ω and takes harmonics of it. In every case, once the φ's are decided, fitting θ is the same least-squares problem.
Step 1 — the polynomial case, made concrete. Suppose you want to fit ŷ = θ₀ + θ₁x + θ₂x² + ⋯ + θ_d x^d to n data points (x_i, y_i). Build a matrix Φ — a Vandermonde matrix — whose i-th row is every power of that row's x_i, from 0 up to d:
The predictions for every row, all at once, are exactly the matrix-vector product Φθ — row i of that product is θ₀ + θ₁x_i + θ₂x_i² + ⋯ + θ_d x_i^d, precisely the polynomial evaluated at x_i. Fitting by least squares means minimizing the sum of squared errors between predictions and targets:
Compare this to plain OLS's objective from 2.3.1, ‖y − Xθ‖², where X is the ordinary design matrix of raw features. The two expressions are identical in form — Φ is simply standing in the exact slot X occupied. Differentiating and setting the gradient to zero, exactly as in the plain-OLS derivation:
— the exact normal equations from 2.3.1, with Φ written where X used to be. No new algebra, no new solver, no new assumption was needed to get here: the Gauss-Markov theorem's guarantee of the best linear unbiased estimator, and every ridge/lasso penalty developed elsewhere in this module, transfer over unchanged simply by reading "design matrix" as "matrix of basis-function values" instead of "matrix of raw features."
Step 2 — nothing above used powers of x. Re-read the argument and check what property of Φ it actually leaned on. It only used that Φ is some fixed n×(K+1) matrix of numbers, that θ enters the prediction as Φθ (a linear combination of Φ's columns), and that the loss is ‖y − Φθ‖². Nothing in the gradient computation or the resulting normal equations referenced "power of x" anywhere. So replace the Vandermonde columns with any fixed, computable transformations φ₁(x), φ₂(x), …, φ_K(x) of the input — Gaussian bumps, sine/cosine pairs, indicator functions for which region x falls into, anything at all — and build:
— the identical closed form solves it. The only place the choice of basis ever enters the computation is in how the numbers inside Φ get computed in the first place; once that one matrix is built, every downstream step (normal equations, ridge's added λI, lasso's coordinate descent) proceeds exactly as it does for plain OLS on raw features. This is the precise sense in which polynomial, RBF, and Fourier regression are not three different algorithms — they are one algorithm, OLS, run on three different feature tables.
Where this is used: this is exactly why radial-basis-function networks, Fourier-feature regression, and (as the real-world examples below expand on) even certain kernel methods can all be fit with the same closed-form linear algebra as plain OLS — nothing about the derivation above cared what the columns of Φ represent, so any method that reduces to "build a fixed feature table, then fit weights linearly" inherits the entire OLS toolkit for free: closed-form solutions, the ridge/lasso penalties from this module, and the statistical guarantees this chapter has built up around linear estimators.
Fixed amber-ticked Gaussian bumps are centered at evenly spaced points; the intro reveals the fit as it's actually computed — one bump's own scaled contribution appearing at a time, summing into the solid purple running total. Afterward, drag the slider from 2 to 20 bumps and watch the same underfit-to-overfit pattern Module 1 showed with polynomial degree, driven this time by how many basis functions the model is allowed.
Toggle between three basis families from the Formula section above and watch how differently their φ_k(x) shapes behave away from the center of the domain: polynomial powers are global and grow fastest near the edges (the root cause of Runge's phenomenon in the Pitfall section below); Gaussian bumps stay local, decaying to zero away from their own center; Fourier terms oscillate everywhere but stay bounded.
All three tabs do the same two-step recipe the derivation describes: build a design matrix whose columns are basis-function values (Gaussian bumps here, not polynomial powers), then solve for θ exactly as OLS does. The library tab shows the same pattern also holds for the polynomial special case — PolynomialFeatures followed by LinearRegression is precisely "build Φ, then run OLS," expressed as an off-the-shelf pipeline.
#include <bits/stdc++.h>
using namespace std;
double trueFn(double x) { return 1.3 * sin(1.15 * x) + 0.22 * x - 1.1; }
struct Rng {
uint64_t state;
explicit Rng(uint64_t seed) : state(seed) {}
double next() {
state = state * 6364136223846793005ULL + 1442695040888963407ULL;
uint32_t xorshifted = (uint32_t)(((state >> 18u) ^ state) >> 27u);
uint32_t rot = (uint32_t)(state >> 59u);
uint32_t out = (xorshifted >> rot) | (xorshifted << ((32 - rot) & 31));
return (double)out / 4294967295.0;
}
double uniform(double lo, double hi) { return lo + next() * (hi - lo); }
double gaussian() {
double u1 = max(next(), 1e-12), u2 = next();
return sqrt(-2.0 * log(u1)) * cos(2.0 * M_PI * u2);
}
};
// Column 0 is the intercept; columns 1..K are Gaussian bumps -- the same Phi
// matrix from the derivation, just with radial-basis columns instead of powers of x.
vector<vector<double>> rbfDesignMatrix(const vector<double>& xs,
const vector<double>& centers, double width) {
vector<vector<double>> Phi(xs.size(), vector<double>(centers.size() + 1));
for (size_t i = 0; i < xs.size(); i++) {
Phi[i][0] = 1.0;
for (size_t k = 0; k < centers.size(); k++) {
double d = xs[i] - centers[k];
Phi[i][k + 1] = exp(-(d * d) / (2 * width * width));
}
}
return Phi;
}
vector<double> solveLinear(vector<vector<double>> A, vector<double> b) {
int n = (int)b.size();
for (int col = 0; col < n; col++) {
int pivot = col;
for (int r = col + 1; r < n; r++)
if (fabs(A[r][col]) > fabs(A[pivot][col])) pivot = r;
swap(A[col], A[pivot]);
swap(b[col], b[pivot]);
double d = fabs(A[col][col]) < 1e-12 ? 1e-12 : A[col][col];
for (int r = col + 1; r < n; r++) {
double f = A[r][col] / d;
for (int c = col; c < n; c++) A[r][c] -= f * A[col][c];
b[r] -= f * b[col];
}
}
vector<double> x(n);
for (int r = n - 1; r >= 0; r--) {
double sum = b[r];
for (int c = r + 1; c < n; c++) sum -= A[r][c] * x[c];
x[r] = sum / A[r][r];
}
return x;
}
int main() {
Rng rng(2024);
int n = 40;
vector<double> xs(n), ys(n);
for (int i = 0; i < n; i++) {
xs[i] = rng.uniform(0.4, 9.6);
ys[i] = trueFn(xs[i]) + 0.28 * rng.gaussian();
}
sort(xs.begin(), xs.end());
int K = 8;
vector<double> centers(K);
for (int k = 0; k < K; k++) centers[k] = 0.4 + k * (9.6 - 0.4) / (K - 1);
double width = max(0.35, (9.6 - 0.4) / (K - 1) * 0.65);
auto Phi = rbfDesignMatrix(xs, centers, width);
int p = K + 1;
vector<vector<double>> PhiTPhi(p, vector<double>(p, 0.0));
vector<double> PhiTy(p, 0.0);
for (int i = 0; i < p; i++) {
for (int j = 0; j < p; j++)
for (int r = 0; r < n; r++) PhiTPhi[i][j] += Phi[r][i] * Phi[r][j];
PhiTPhi[i][i] += 1e-6;
for (int r = 0; r < n; r++) PhiTy[i] += Phi[r][i] * ys[r];
}
vector<double> theta = solveLinear(PhiTPhi, PhiTy);
double sse = 0.0;
for (int i = 0; i < n; i++) {
double pred = theta[0];
for (int k = 0; k < K; k++) {
double d = xs[i] - centers[k];
pred += theta[k + 1] * exp(-(d * d) / (2 * width * width));
}
double err = pred - ys[i];
sse += err * err;
}
cout << "training MSE: " << sse / n << "\n";
cout << "Same normal-equations solver as plain OLS -- only Phi's columns changed.\n";
return 0;
}- Polynomial trend fitting in physics and engineering measurements — fitting a low-degree polynomial (quadratic, cubic) to calibration curves, thermal expansion data, or projectile-motion measurements is a direct, everyday use of the exact Vandermonde-matrix derivation above: the physical relationship is often genuinely well-approximated by a low-degree polynomial, and the fit is produced by the same OLS solver used for any other linear regression.
- Radial basis function networks in geospatial interpolation — estimating a continuous surface (temperature, elevation, pollutant concentration) from scattered sensor or survey readings often places RBF centers at or near the sensor locations and fits weights exactly as the diagram above does, just in two spatial dimensions instead of one. This basis-function view is closely related to (though not identical in its statistical assumptions to) geostatistical kriging.
- Fourier-feature regression for periodic or seasonal signals — modeling electricity demand, retail sales, or vibration/audio measurements that repeat on a known cycle (daily, weekly, yearly) works well with sine/cosine features at the known frequency and its harmonics; the seasonal component becomes literally a handful of columns fit by the same linear regression as everything else in this module, often combined with a polynomial trend term in the same design matrix.
- The kernel trick, previewed briefly. Some very high- or even infinite-dimensional basis expansions are computationally infeasible to build explicitly column by column — but for certain algorithms, the fitting math only ever needs dot products between basis vectors, never the basis vectors themselves, and a kernel function can compute that dot product directly without ever forming
Φ. This is the core idea behind kernel ridge regression and support vector machines, introduced properly in the Linear Algebra chapter's kernel-methods lesson and developed in full in a later Kernel Methods module — worth flagging here only because it is, at its core, the exact same basis-expansion idea pushed to an extreme where the basis is never written down. - RBF networks as an early alternative to standard neural networks for function approximation and control — before deep learning's rise, radial-basis-function networks (a hidden layer of Gaussian-bump units, exactly as diagrammed above, feeding a linear output layer) were a popular universal-function-approximation architecture, precisely because the output layer's weights could be fit by ordinary linear regression once the bump centers and widths were chosen.
- Time-series trend-plus-seasonality decomposition — many classical forecasting pipelines fit a single linear regression whose design matrix concatenates a low-degree polynomial trend basis with a Fourier seasonal basis in one design matrix, estimating both components' weights in one normal-equations solve rather than as two separate models.
- Extrapolating a high-degree polynomial outside the training range is a much worse failure mode than the interpolation-region wiggling Module 1's overfitting lesson (2.1.4) covers, and it deserves its own name: Runge's phenomenon for the oscillation that appears near the edges of the fitted range, and simple algebra for what happens beyond it — since a degree-
dpolynomial's leading termθ_d x^deventually dominates every other term as|x|grows, any high-degree polynomial fit races off toward±∞just past the edge of the data it was trained on, often within a small step outside the training range. RBF and Fourier bases don't remove this risk entirely, but they fail more gently near the edges of the data — a Gaussian bump decays to (near) zero far from its center rather than exploding, which is one reason basis choice matters beyond just "fits the training data well." - Choosing a basis with no domain reasoning about the true relationship's shape. Fitting a polynomial basis to an obviously periodic signal, or fitting a handful of widely-spaced RBF bumps to a signal with sharp discontinuities, wastes basis functions representing shapes the truth doesn't have and starves the model of the shapes it does need. The choice of basis is itself a modeling decision that deserves the same domain thinking as choosing which raw features to collect in the first place — it is not a purely mechanical step to skip past.
- Forgetting that more basis functions is still more capacity. Adding more polynomial terms, more RBF bumps, or more Fourier harmonics increases the model's flexibility exactly the way increasing polynomial degree did in Module 1's overfitting lesson (2.1.4) — because, per the derivation above, it is the identical mechanism: more columns in
Φ, more free parameters inθ, the same bias-variance trade-off. Basis expansion does not let you escape that trade-off; it only relocates the decision from "how curvy should the fitted line be" to "how many, and which, fixed feature functions do I supply" — and the same held-out validation discipline from 2.1.4 is exactly as necessary here as it is for raw polynomial degree.
Going deeper
Basis-function regression's biggest structural weakness is that every polynomial or Fourier basis function is global — it's nonzero across the entire domain, so a coefficient chosen to fit the data well near one point can distort the fitted curve everywhere else too, which is exactly the mechanism behind Runge's phenomenon flagged above. Splines, the very next lesson (2.3.7), are best understood as basis-function regression with one specific, carefully engineered choice of basis: piecewise polynomials defined only on small local intervals between fixed "knot" points, glued together with continuity constraints (matching value, and usually matching first and second derivatives) exactly at each knot. Each piece only has influence near its own interval, so a stray wiggle forced by one region's data can no longer propagate across the whole curve the way a single misbehaving global polynomial term can. Nothing about the fitting machinery changes — it is still, underneath, the same θ̂ = (Φᵗ Φ)⁻¹Φᵗy this lesson derived — only the basis functions filling Φ's columns get smarter about being local instead of global. That's the whole motivation for the next lesson existing at all: a better basis, not a different algorithm.
A colleague fits y = θ0 + θ1·sin(x) + θ2·cos(x) + θ3·sin(2x) + θ4·cos(2x) to some data using ordinary least squares, and the resulting curve clearly oscillates. They say: 'this can't be linear regression, look how curvy it is.' What's wrong with that reasoning, and what would you actually call this model?
The reasoning conflates two different meanings of 'linear.' The FITTED CURVE is certainly not a straight line in x -- but 'linear regression' refers to the model being linear in its PARAMETERS theta, not in x. Here, once sin(x), cos(x), sin(2x), and cos(2x) are treated as fixed, precomputed feature columns (exactly like the Phi matrix in the derivation above), the prediction is a plain linear combination of those columns with weights theta0..theta4 -- so this is fit by, and is, ordinary least squares, with the closed-form normal equations applying completely unchanged. The right name for it is Fourier-feature (basis-expansion) regression: a specific, well-motivated choice of basis for a signal expected to have periodic structure, not a departure from linear regression at all.
"Linear" in linear regression describes the parameters, not the shape of the fitted curve. Basis expansion exploits that gap: replace raw x with any fixed set of computable features φ_k(x) — powers, Gaussian bumps, sinusoids, or anything else a domain calls for — and every tool built for plain OLS (closed-form solving, Gauss-Markov, ridge, lasso) carries over unchanged, because the underlying optimization problem, ‖y − Φθ‖², is unchanged. What does change is capacity, and basis expansion doesn't sidestep the bias-variance trade-off from Module 1 — it only moves where you make the capacity decision, from "how high a polynomial degree" to "how many, and which, basis functions." The next lesson, splines, is exactly this idea again with one deliberately smarter choice of local basis.