Local regression

The subsets of data used for each weighted least squares fit in LOESS are determined by a nearest neighbors algorithm. A user-specified input to the procedure called the "bandwidth" or "smoothing parameter" determines how much of the data is used to fit each local polynomial. The smoothing parameter, ${\displaystyle \alpha }$, is the fraction of the total number n of data points that are used in each local fit. The subset of data used in each weighted least squares fit thus comprises the ${\displaystyle n\alpha }$ points (rounded to the next largest integer) whose explanatory variables' values are closest to the point at which the response is being estimated.[7]

Since a polynomial of degree k requires at least k + 1 points for a fit, the smoothing parameter ${\displaystyle \alpha }$ must be between ${\displaystyle \left(\lambda +1\right)/n}$ and 1, with ${\displaystyle \lambda }$ denoting the degree of the local polynomial.

${\displaystyle \alpha }$ is called the smoothing parameter because it controls the flexibility of the LOESS regression function. Large values of ${\displaystyle \alpha }$ produce the smoothest functions that wiggle the least in response to fluctuations in the data. The smaller ${\displaystyle \alpha }$ is, the closer the regression function will conform to the data. Using too small a value of the smoothing parameter is not desirable, however, since the regression function will eventually start to capture the random error in the data.

The local polynomials fit to each subset of the data are almost always of first or second degree; that is, either locally linear (in the straight line sense) or locally quadratic. Using a zero degree polynomial turns LOESS into a weighted moving average. Higher-degree polynomials would work in theory, but yield models that are not really in the spirit of LOESS. LOESS is based on the ideas that any function can be well approximated in a small neighborhood by a low-order polynomial and that simple models can be fit to data easily. High-degree polynomials would tend to overfit the data in each subset and are numerically unstable, making accurate computations difficult.

As mentioned above, the weight function gives the most weight to the data points nearest the point of estimation and the least weight to the data points that are furthest away. The use of the weights is based on the idea that points near each other in the explanatory variable space are more likely to be related to each other in a simple way than points that are further apart. Following this logic, points that are likely to follow the local model best influence the local model parameter estimates the most. Points that are less likely to actually conform to the local model have less influence on the local model parameter estimates.

The traditional weight function used for LOESS is the tri-cube weight function,

${\displaystyle w(d)=(1-|d|^{3})^{3}}$

where d is the distance of a given data point from the point on the curve being fitted, scaled to lie in the range from 0 to 1.[7]

However, any other weight function that satisfies the properties listed in Cleveland (1979) could also be used. The weight for a specific point in any localized subset of data is obtained by evaluating the weight function at the distance between that point and the point of estimation, after scaling the distance so that the maximum absolute distance over all of the points in the subset of data is exactly one.

Consider the following generalisation of the linear regression model with a metric ${\displaystyle w(x,z)}$ on the target space ${\displaystyle \mathbb {R} ^{m}}$ that depends on two parameters, ${\displaystyle x,z\in \mathbb {R} ^{p}}$. Assume that the linear hypothesis is based on ${\displaystyle p}$ input parameters and that, as customary in these cases, we embed the input space ${\displaystyle \mathbb {R} ^{p}}$ into ${\displaystyle \mathbb {R} ^{p+1}}$ as ${\displaystyle x\mapsto {\hat {x}}:=(1,x)}$, and consider the following loss function

${\displaystyle \operatorname {RSS} _{x}(A)=\sum _{i=1}^{N}(y_{i}-A{\hat {x}}_{i})^{T}w_{i}(x)(y_{i}-A{\hat {x}}_{i}).}$

Here, ${\displaystyle A}$ is an ${\displaystyle m\times (p+1)}$ real matrix of coefficients, ${\displaystyle w_{i}(x):=w(x_{i},x)}$ and the subscript i enumerates input and output vectors from a training set. Since ${\displaystyle w}$ is a metric, it is a symmetric, positive-definite matrix and, as such, there is another symmetric matrix ${\displaystyle h}$ such that ${\displaystyle w=h^{2}}$. The above loss function can be rearranged into a trace by observing that ${\displaystyle y^{T}wy=(hy)^{T}(hy)=\operatorname {Tr} (hyy^{T}h)=\operatorname {Tr} (wyy^{T})}$. By arranging the vectors ${\displaystyle y_{i}}$ and ${\displaystyle {\hat {x}}_{i}}$ into the columns of a ${\displaystyle m\times N}$ matrix ${\displaystyle Y}$ and an ${\displaystyle (p+1)\times N}$ matrix ${\displaystyle {\hat {X}}}$ respectively, the above loss function can then be written as

${\displaystyle \operatorname {Tr} (W(x)(Y-A{\hat {X}})^{T}(Y-A{\hat {X}}))}$

where ${\displaystyle W}$ is the square diagonal ${\displaystyle N\times N}$ matrix whose entries are the ${\displaystyle w_{i}(x)}$s. Differentiating with respect to ${\displaystyle A}$ and setting the result equal to 0 one finds the extremal matrix equation

${\displaystyle A{\hat {X}}W(x){\hat {X}}^{T}=YW(x){\hat {X}}^{T}.}$

Assuming further that the square matrix ${\displaystyle {\hat {X}}W(x){\hat {X}}^{T}}$ is non-singular, the loss function ${\displaystyle \operatorname {RSS} _{x}(A)}$ attains its minimum at

${\displaystyle A(x)=YW(x){\hat {X}}^{T}({\hat {X}}W(x){\hat {X}}^{T})^{-1}.}$

A typical choice for ${\displaystyle w(x,z)}$ is the Gaussian weight

${\displaystyle w(x,z)=\exp \left(-{\frac {\|x-z\|^{2}}{2\alpha ^{2}}}\right)}$

• Fortran implementation
• C implementation (from the R project)
• Lowess implementation in Cython by Carl Vogel
• Python implementation (in Statsmodels)
• LOESS Smoothing in Excel
• LOESS implementation in pure Julia
• JavaScript implementation
• Java implementation

 This article incorporates public domain material from the National Institute of Standards and Technology.