Locally Weighted Regression

a straight line is too rigid for a wiggly world, and a global degree-9 polynomial is a hostage negotiation. 𐃏 locally weighted regression (LWR — and its robust cousin LOWESS) takes the diplomatic route: fit the simplest possible model, but fit it freshly at every query point, paying attention only to the training points nearby.

motivation#

• linear regression commits to one \(\theta\) for the whole input space. if the true \(f\) bends, the residue of that commitment is bias everywhere.
• the fix need not be a fancier global family. any smooth function is locally linear — taylor says so — so a linear fit weighted toward a neighbourhood of \(x_0\) can track an arbitrary smooth \(f\).
• the price: there is no longer a “trained model”. LWR is memory-based and non-parametric — like k-nearest neighbours, it keeps the entire training set and does all of its work at prediction time. training is \(O(1)\); every query costs a fresh weighted least-squares solve. 𐃏

the estimator#

weighted least squares at a query point#

fix a query \(x_0\). assign each training point a weight \(w_i(x_0) \ge 0\) that decays with distance from \(x_0\), then solve the weighted least-squares problem

\begin{equation} \hat\theta(x_0) = \arg\min_{\theta} \sum_{i=1}^{n} w_i(x_0)\,\big(y_i - \theta^\top x_i\big)^2 . \end{equation}

stack the weights into \(W = \operatorname{diag}(w_1(x_0), \dots, w_n(x_0))\). setting the gradient \(-2X^\top W(y - X\theta)\) to zero gives the closed form

\begin{equation} \boxed{\;\hat\theta(x_0) = (X^\top W X)^{-1} X^\top W y\;} \end{equation}

and the prediction is \(\hat f(x_0) = x_0^\top \hat\theta(x_0)\). ordinary least squares is the special case \(W = I\) — LWR is just OLS wearing distance-tinted glasses.

• in practice one augments \(x\) with a constant (and, for local polynomial regression, powers of \(x - x_0\)): degree 0 is a weighted average, degree 1 a weighted line, degree 2 a weighted parabola.
• note \(\hat\theta\) depends on \(x_0\) through \(W\): a different regression for every query. the final curve \(x_0 \mapsto \hat f(x_0)\) is smooth even though no single global model produced it.
• each prediction is linear in \(y\): \(\hat f(x_0) = \ell(x_0)^\top y\) for a weight vector \(\ell(x_0)\) that depends only on the inputs. LWR is a linear smoother, so an effective degrees of freedom exists: \(\mathrm{df} = \operatorname{tr}(L)\), where row \(i\) of \(L\) is \(\ell(x_i)^\top\) (Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome, 2009).

kernels#

the weight function is a kernel \(w_i(x_0) = K_\lambda(x_0, x_i)\). two standard choices:

• gaussian: \(K(x_0, x_i) = \exp\!\big(-\lVert x_i - x_0\rVert^2 / 2\tau^2\big)\). every point gets some weight; \(\tau\) is the bandwidth.
• tricube (the LOWESS default): with \(u = \lvert x_i - x_0\rvert / \lambda\),

\begin{equation} K(u) = \begin{cases} (1 - u^3)^3 & u < 1 \\ 0 & u \ge 1, \end{cases} \end{equation}

compactly supported (points beyond \(\lambda\) contribute exactly nothing — cheaper, and the fit at \(x_0\) is immune to far-away outliers) and smooth at the boundary of its support (Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome, 2009).

• LOWESS makes the bandwidth adaptive: \(\lambda(x_0)\) is set to the distance to the \(k\)-th nearest neighbour, so a fixed fraction (the “span”) of the data is in play everywhere — wide windows where data are sparse, narrow where dense.¹

bandwidth is the bias-variance dial#

everything you know from the bias-variance decomposition applies with \(\tau\) as the knob:

the sweet spot minimises estimated test error, and since each fit is cheap and the model never “trains”, leave-one-out cross validation is the natural selector — for linear smoothers it even has a closed form, no refitting required.

connection to kernel smoothing#

LWR is the polished member of a family (Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome, 2009):

• degree 0 is the nadaraya–watson kernel smoother: minimising \(\sum_i w_i (y_i - \theta_0)^2\) over a constant \(\theta_0\) gives the kernel-weighted average

\begin{equation} \hat f(x_0) = \frac{\sum_i K_\lambda(x_0, x_i)\, y_i}{\sum_i K_\lambda(x_0, x_i)}. \end{equation}

• but a locally-constant fit has a defect at the edges: near a boundary the window is one-sided, so if \(f\) is sloping, the average of the neighbours sits systematically above or below the truth. this is boundary bias, and it is first-order — it does not vanish quickly with more data.
• local linear regression fixes it automatically: fitting a slope lets the estimator extrapolate the trend across the asymmetric window. local linear fits are exact for linear \(f\) everywhere, boundaries included, killing the first-order term in the bias expansion (Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome, 2009). the code below makes this concrete — watch the two estimators disagree at \(x_0 = 0\) and \(x_0 = 1\).
• local quadratic fits go one better in curved interior regions (locally-linear fits trim the hills and fill the valleys — curvature bias) at the price of higher variance. degree 1 is the everyday compromise.
• the general kernel view of learning — feature spaces, gram matrices, representer theorems — lives on the kernel methods page; here the kernel is doing a humbler job, defining “nearby”.

one warning from the curse of dimensionality: local methods lean entirely on neighbourhoods being small and populated. in high dimensions neighbourhoods that contain enough points stop being local, boundary points become the majority, and LWR quietly degrades. it is a low-dimensional specialist.

running means, loess and local likelihood#

the estimator above is the polished form. the classical route to it runs in steps — fixed-count neighbourhoods first, then cleveland’s loess, then the likelihood reading that carries local fitting beyond gaussian responses — and each step makes one of the claims above precise. 𐃏

running means and running lines#

before kernels, the crudest local methods use a fixed count of neighbours. with the data sorted, the symmetric \(k\)-nearest neighbourhood of \(x_i\) is the index set \(N_k(x_i) = \{\max(1, i-k), \dots, \min(n, i+k)\}\), of size \(n_{ik}\): \(2k + 1\) in the interior, truncated at the ends (for \(i = 2, k = 2\) it is \(\{1,2,3,4\}\)).

the running mean is a box kernel: it jumps whenever a point enters the window, is never differentiable, and suffers boundary bias because its window is one-sided at the ends. the running line removes the first-order boundary bias but is still only piecewise linear.

loess in detail#

fixed bandwidth versus fixed neighbour count#

a kernel smoother fixes the window width \(h\); loess fixes the number of points \(k\). so loess has roughly constant variance across \(x\) but becomes biased where data are sparse (far points are pulled in to make up the count), while a fixed-\(h\) kernel smoother keeps its bias local but its variance grows where few points fall in the window. on Wage, ksmooth(age, wage, kernel = "normal") with \(h = 5\) (gaussian sd about 1.9 years in ksmooth’s scaling) and \(h = 2\) (sd about 0.7, nearly the mean at a single integer age) give 99.1 and 93.9 at age 70, the second visibly jagged. a model global in some variables and local in another, \(y = \beta_0(z) + \beta_1(z)x + \varepsilon\), is a varying coefficient model.

boundary bias, proved#

local likelihood#

with gaussian errors,

\begin{equation} -2\,l(y_i \mid \boldsymbol\beta) = \sigma^{-2}(y_i - \mathbf x_i^\top\boldsymbol\beta)^2 + \log(2\pi\sigma^2)\,, \end{equation}

so the local weighted least-squares criterion is, up to constants, \(-2\sum_i w_i(x_0)\,l(y_i \mid \boldsymbol\beta)\). written this way the gaussian log-likelihood can be swapped for a binomial or poisson one with a link: the same kernel weights give local logistic or local poisson regression, the route from smoothing to non-parametric GLMs and additive models.

proofs#

in code#

LWR from scratch on a noisy sine — gaussian kernel, local polynomials of degree 0 (nadaraya–watson) and degree 1, and a bandwidth sweep:

import numpy as np

rng = np.random.default_rng(7)
f = lambda x: np.sin(2 * np.pi * x)
n, sigma = 60, 0.25
x = np.sort(rng.uniform(0, 1, n))
y = f(x) + rng.normal(0, sigma, n)

def lwr(x0, x, y, tau, degree=1):
    """fit a weighted polynomial at query x0, return prediction there."""
    w = np.exp(-((x - x0) ** 2) / (2 * tau ** 2))     # gaussian kernel
    A = np.vander(x - x0, degree + 1, increasing=True) # [1, (x-x0), ...]
    W = np.diag(w)
    theta = np.linalg.solve(A.T @ W @ A, A.T @ W @ y)  # (A'WA)^{-1} A'Wy
    return theta[0]                                    # value at x0

xg = np.linspace(0, 1, 101)
print(f"{'tau':>6} {'rmse (deg 0)':>13} {'rmse (deg 1)':>13}")
for tau in [0.02, 0.05, 0.1, 0.2, 0.5, 2.0]:
    for deg in (0, 1):
        yhat = np.array([lwr(x0, x, y, tau, deg) for x0 in xg])
        rmse = np.sqrt(((yhat - f(xg)) ** 2).mean())
        if deg == 0:
            row = f"{tau:>6.2f} {rmse:>13.4f}"
        else:
            row += f" {rmse:>13.4f}"
    print(row)

# a few point predictions at the best bandwidth, incl. the boundary
tau = 0.1
print("\n x0    truth   deg-0   deg-1")
for x0 in [0.00, 0.25, 0.50, 0.75, 1.00]:
    p0, p1 = lwr(x0, x, y, tau, 0), lwr(x0, x, y, tau, 1)
    print(f"{x0:.2f} {f(x0):>8.3f} {p0:>7.3f} {p1:>7.3f}")

   tau  rmse (deg 0)  rmse (deg 1)
  0.02        0.1327        0.1268
  0.05        0.0773        0.0674
  0.10        0.1620        0.1231
  0.20        0.3455        0.2812
  0.50        0.5892        0.4374
  2.00        0.6956        0.4523

 x0    truth   deg-0   deg-1
0.00    0.000   0.253  -0.052
0.25    1.000   0.843   0.843
0.50    0.000  -0.052  -0.022
0.75   -1.000  -0.808  -0.806
1.00   -0.000  -0.459  -0.043

three lessons in one table:

• bandwidth u-shape. rmse bottoms out near \(\tau = 0.05\) and climbs in both directions: at \(\tau = 0.02\) the curve chases individual noisy points (variance), by \(\tau = 0.5\) the “local” window spans half a period of the sine (bias), and \(\tau = 2\) is effectively a global fit to a sine wave.
• boundary bias, live. at the interior points the two degrees agree closely. at \(x_0 = 0\) and \(x_0 = 1\) the degree-0 smoother misses by \(0.25\) and \(0.46\) — its one-sided window averages points that all sit above (respectively below) the truth. degree 1 extrapolates the local trend and lands within \(0.06\). exactly the theory’s prediction.
• interior peaks are trimmed. both estimators read \(0.84\) at the crest \(x_0 = 0.25\) where the truth is \(1\) — smoothing flattens local maxima (curvature bias); a local quadratic would recover most of it.

complexity check: every query re-solves a small \((d{+}1)\times(d{+}1)\) system built from all \(n\) points — \(O(n)\) work per prediction, versus \(O(1)\) for a fitted parametric model. with a compact kernel and a k-d tree the constant improves, but the memory-based character remains.

see also#

• ordinary least squares — the \(W = I\) special case, and the engine inside every local fit
• k-nearest neighbours — the same memory-based bet with a cruder (rectangular) kernel
• bias-variance decomposition — what the bandwidth is really adjusting
• cross validation — choosing \(\tau\), with a free closed form for linear smoothers
• curse of dimensionality — why this page lives happily in 1–3 dimensions and dies in 30

References#

Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome (2009). The Elements of Statistical Learning, Springer.