Multivariate Local Polynomial Regression

Why Multivariate Local Polynomial Regression?

When the predictor $X$ has $d$ coordinates, a local fit must use nearby observations in a $d$-dimensional neighborhood. The central construction is to scale each nearby displacement by a bandwidth $h$, then fit a polynomial in those scaled coordinates by kernel-weighted least squares. This keeps the local design on a stable scale as $h$ shrinks and makes the roles of the polynomial basis and kernel moments explicit.

The discussion below concerns an interior point $x$, where the predictor density is positive and the kernel neighborhood does not meet the boundary of the predictor domain. A single scalar bandwidth $h$ is used for all coordinates; unequal coordinate bandwidths are a separate issue.

Scaled Coordinates, Basis, and Moment Matrix

Let $X_i\in\mathbb R^d$ and $Y_i\in\mathbb R$, for $i=1,\ldots,n$, be independent observations from the same distribution. Write $m(z)=E[Y\mid X=z]$, and estimate $m(x)$ at a fixed point $x$. For an observation near $x$, define its scaled displacement $U_i=(X_i-x)/h$.

A multi-index is a vector $\alpha=(\alpha_1,\ldots,\alpha_d)$ of nonnegative integers. Write $|\alpha|=\sum_{j=1}^d\alpha_j$, $\alpha!=\prod_{j=1}^d\alpha_j!$, and $u^\alpha=\prod_{j=1}^d u_j^{\alpha_j}$. The degree-$p$ polynomial basis $r_p(u)$ contains one coordinate $u^\alpha$ for every $|\alpha|\le p$, with the constant coordinate first. Its dimension is $q_p=\binom{d+p}{p}$. Let $e_0$ denote the vector selecting that first, constant coordinate.

Definition: The degree-$p$ multivariate local polynomial estimator minimizes $$ \sum_{i=1}^n K(U_i)\{Y_i-r_p(U_i)^\mathsf{T}b\}^2 $$ over $b\in\mathbb R^{q_p}$, and sets $\hat m(x)=e_0^\mathsf{T}\hat b$. Here $K:\mathbb R^d\to[0,\infty)$ is a kernel and $h>0$ is the bandwidth. The factor $h^{-d}$, sometimes included in the kernel weights, does not change the minimizer.

The coordinates of $b$ approximate the Taylor coefficients in scaled coordinates: the coefficient for $u^\alpha$ is approximately $h^{|\alpha|}D^\alpha m(x)/\alpha!$. The matrix that controls the local fit is the kernel moment matrix

S=\int_{\mathbb R^d}r_p(u)r_p(u)^\mathsf{T}K(u)\,du.

For a stable population fit, assume $S$ is finite and positive definite. This means no nonzero polynomial in the chosen basis has zero squared integral under the kernel weights. The sample weighted least-squares matrix must also be invertible for the estimator to be defined.

Worked Example: A Bivariate Local Linear Basis

Take $d=2$ and $p=1$, with the product uniform kernel $K(u_1,u_2)=\tfrac14$ on $[-1,1]^2$ and zero elsewhere. The basis is $r_1(u)=(1,u_1,u_2)^\mathsf{T}$. Since the kernel is symmetric, the first moments and the mixed moment vanish. Also, $$ \int K(u)\,du=1,\qquad \int u_1^2K(u)\,du=\int u_2^2K(u)\,du=\frac13. $$ Thus $S=\operatorname{diag}(1,1/3,1/3)$, which is positive definite. The fit estimates a local intercept and two local slopes. There are three coefficients because the bivariate linear basis has one constant and one term for each coordinate.

Interior Bias: A Moment-Matrix Derivation

Let $f$ be the density of $X$, with $f(x)>0$, and define the population weighted least-squares coefficient $\beta_h$ by minimizing $$ \int K(u)\{m(x+hu)-r_p(u)^\mathsf{T}b\}^2f(x+hu)\,du. $$ This population target separates the smoothing bias from the sampling variation. Suppose $m$ has continuous derivatives through order $p+1$ near $x$, $f$ is continuous and positive at $x$, and $K$ is bounded with compact support. For each multi-index $\alpha$, put $$ R_\alpha=\int r_p(u)u^\alpha K(u)\,du. $$ The generic leading bias expansion is as follows.

Theorem: Under the conditions just stated, if $S$ is positive definite, then $$ e_0^\mathsf{T}\beta_h-m(x) =h^{p+1}\sum_{|\alpha|=p+1} \frac{D^\alpha m(x)}{\alpha!}\, e_0^\mathsf{T}S^{-1}R_\alpha +o(h^{p+1}). $$ The coefficient can be zero for a particular kernel or symmetry; in that case the first nonzero bias term is of smaller order than the generic $h^{p+1}$ term.

Proof. Taylor expansion on the compact kernel support gives $$ m(x+hu)=\sum_{|\alpha|\le p}\frac{h^{|\alpha|}D^\alpha m(x)}{\alpha!}u^\alpha +h^{p+1}\sum_{|\alpha|=p+1}\frac{D^\alpha m(x)}{\alpha!}u^\alpha +o(h^{p+1}), $$ uniformly over the support. The first sum is exactly $r_p(u)^\mathsf{T}\beta^{(p)}$ for a coefficient vector $\beta^{(p)}$, whose constant coordinate is $m(x)$. The population normal equations have matrix $$ M_h=\int r_p(u)r_p(u)^\mathsf{T}K(u)f(x+hu)\,du $$ and right-hand side $\int r_p(u)K(u)m(x+hu)f(x+hu)\,du$. Subtracting $M_h\beta^{(p)}$ from that right-hand side leaves $$ h^{p+1}f(x)\sum_{|\alpha|=p+1} \frac{D^\alpha m(x)}{\alpha!}R_\alpha+o(h^{p+1}). $$ Indeed, continuity of $f$ gives $f(x+hu)=f(x)+o(1)$ uniformly on the support, and the Taylor remainder is uniform there. Also $M_h\to f(x)S$, which is invertible because $f(x)>0$ and $S$ is positive definite. Multiplying by $M_h^{-1}$, the factor $f(x)$ cancels in the leading term. Finally, $e_0^\mathsf{T}\beta^{(p)}=m(x)$. This proves the expansion. $\square$

Worked Example: Bivariate Local Linear Bias

Use the kernel and basis from the first example, and suppose the design density is constant near $x=(0,0)$. Let $$ m(z_1,z_2)=2+3z_1-2z_2+4z_1^2+6z_1z_2+2z_2^2. $$ For $p=1$, the leading omitted terms have degree two. The required second derivatives are $D^{(2,0)}m=8$, $D^{(1,1)}m=6$, and $D^{(0,2)}m=4$. Since $S^{-1}=\operatorname{diag}(1,3,3)$, the intercept coordinate of $S^{-1}R_{(2,0)}$ is $\int u_1^2K(u)\,du=1/3$; that of $S^{-1}R_{(0,2)}$ is also $1/3$. For $R_{(1,1)}$, it is $\int u_1u_2K(u)\,du=0$. Therefore $$ e_0^\mathsf{T}\beta_h-m(0) =h^2\left\{\frac{8}{2}\frac13+6(0)+\frac{4}{2}\frac13\right\} +o(h^2) =2h^2+o(h^2). $$ The mixed-curvature term contributes nothing here because the product kernel is symmetric in each coordinate.

Variance and Dimension

The variance depends on squared kernel weights, so it involves a second moment matrix distinct from $S$: $$ S_2=\int r_p(u)r_p(u)^\mathsf{T}K(u)^2\,du. $$ Assume additionally that $Y=m(X)+\varepsilon$, with $E[\varepsilon\mid X]=0$, conditional variance $\operatorname{Var}(\varepsilon\mid X=z)=\sigma^2(z)$ continuous at $x$, and $\sigma^2(x)>0$. With the kernel and moment conditions above, and $h\to0$, $nh^d\to\infty$, the leading variance is

Theorem: At an interior point, $$ \operatorname{Var}\{\hat m(x)\}\ \sim\ \frac{\sigma^2(x)}{nh^d f(x)} e_0^\mathsf{T}S^{-1}S_2S^{-1}e_0. $$ In particular, the variance has order $(nh^d)^{-1}$ when the displayed constant is positive.

Derivation. The sample normal equations use the matrix $\sum_i K(U_i)r_p(U_i)r_p(U_i)^\mathsf{T}$. Dividing it by $nh^d$, its expectation tends to $f(x)S$: this follows by the change of variables $X=x+hu$, which contributes the volume factor $h^d$, and continuity of $f$. Thus the inverse sample matrix is approximately $(nh^df(x))^{-1}S^{-1}$. The leading random part of the normal-equation right-hand side is $\sum_i K(U_i)r_p(U_i)\varepsilon_i$. Its covariance is asymptotic to $nh^d f(x)\sigma^2(x)S_2$, by the same change of variables and continuity of the conditional variance. Combining the inverse matrix on both sides of this covariance gives $$ \frac{\sigma^2(x)}{nh^df(x)}S^{-1}S_2S^{-1}. $$ Taking the intercept coordinate yields the stated expression. The Taylor residual in the conditional mean is of order $h^{p+1}$ locally and contributes a smaller variance term than the noise term when $\sigma^2(x)>0$. The approximation uses a growing effective local sample size $nh^d$; requiring suitable finite moments, for example a conditional $2+\delta$ moment for some $\delta>0$, also supports the usual asymptotic normal approximation. $\square$

The factor $h^d$ is important: a neighborhood with radius proportional to $h$ has volume proportional to $h^d$. Thus, at the same sample size and bandwidth, increasing the predictor dimension reduces the number of observations contributing appreciable weight. This geometric effect is one reason multivariate nonparametric estimation can be statistically demanding.

Worked Example: Variance for a Bivariate Local Linear Fit

Continue with the product uniform kernel, now under a uniform design on a square, so $f(x)=1/4$ at an interior point. Let the conditional noise variance be $\sigma^2(x)=4$, and take $p=1$. Since $S^{-1}e_0=e_0$, the variance constant uses only the first entry of $S_2$: $$ e_0^\mathsf{T}S^{-1}S_2S^{-1}e_0 =e_0^\mathsf{T}S_2e_0 =\int K(u)^2\,du =\frac{1}{16}\times 4 =\frac14. $$ There are four unit-area pieces in the support calculation, so the integral is $4/16$. In dimension two, $h^d=h^2$; substituting the values gives $$ \operatorname{Var}\{\hat m(x)\}\sim \frac{4}{nh^2(1/4)}\frac14 =\frac{4}{nh^2}. $$ The local effective sample size is proportional to $nh^2/4$, consistent with the inverse-variance order $1/(nh^2)$.

Exact Polynomial Reproduction and a Common Pitfall

A useful check on the basis is exact reproduction. If the true regression function is a polynomial of degree at most $p$, it can be represented exactly as $r_p(u)^\mathsf{T}b$ after substituting $z=x+hu$. In a noiseless sample with an invertible weighted design matrix, weighted least squares then recovers that polynomial exactly. This is why the local fit’s coordinates must be scaled consistently with the bandwidth.

Worked Example: A Local Quadratic Reproduces a Quadratic Surface

At $x=(0,0)$, let $m(z_1,z_2)=1+2z_1-z_2+3z_1^2+4z_1z_2+5z_2^2$, and use the degree-two basis $(1,u_1,u_2,u_1^2,u_1u_2,u_2^2)^\mathsf{T}$. With $z=hu$, the function becomes $$ m(hu_1,hu_2)=1+(2h)u_1-hu_2+(3h^2)u_1^2+(4h^2)u_1u_2+(5h^2)u_2^2. $$ Thus its coefficient vector in the scaled basis is $(1,2h,-h,3h^2,4h^2,5h^2)^\mathsf{T}$. If the observations are noiseless and the positively weighted design points make the six-column design matrix full rank, this vector gives zero residual and is the unique weighted least-squares solution. In particular, the fitted intercept is exactly $m(0)=1$, for any such bandwidth.

Pitfall: treating a multivariate bandwidth as if it created a one-dimensional neighborhood. The local sample size is proportional to $nh^d$, not $nh$. A variance calculation that omits the dimension-dependent volume factor gives the wrong rate.
Correction: Count the $d$-dimensional volume through the change of variables $X=x+hu$; its Jacobian is $h^d$.
Pitfall: reading the coefficient of $u^\alpha$ as an unscaled derivative. The coefficient approximates $h^{|\alpha|}D^\alpha m(x)/\alpha!$, not the derivative alone.
Correction: Keep the scaled basis $r_p((X_i-x)/h)$ and interpret its coefficients with the corresponding powers of $h$.

Key takeaway: The multivariate fit is governed by the polynomial moment matrix $S$ and the squared-kernel moment matrix $S_2$. At an interior point, the generic bias is of order $h^{p+1}$, while the variance is of order $1/(nh^d)$; symmetry or other moment cancellations can make the leading bias smaller.

Check Your Understanding

Try these questions using the definitions and expansions above.

For $d=3$ and $p=1$, how many coordinates are in $r_p(u)$, and what are they?
Why does a positive-definite moment matrix $S$ matter for defining the population local fit?
For a local linear fit in two dimensions, which feature of a product symmetric kernel made the mixed-curvature term vanish in the worked bias calculation?
Explain why the leading variance involves $K(u)^2$, whereas the moment matrix $S$ involves $K(u)$.
If the bandwidth is halved while $n$ and $d=2$ stay fixed, what happens to the leading variance order, ignoring changes in its constant?