Numerical Conditioning in Local Polynomial Fits

Why Conditioning Matters in a Local Polynomial Fit

A local polynomial estimate is obtained from a weighted least-squares problem. Even when the intended polynomial model is appropriate, the numerical calculation can be unreliable if the weighted design matrix is nearly singular or has columns on very different scales. The resulting coefficients may change substantially under small rounding errors, and predictions can be inaccurate.

Two causes should be separated. Poor scaling can make a well-defined fit difficult to compute in the chosen coordinates. Insufficient local information can make the fit singular, or nearly so, regardless of the computational method. Centering and rescaling help with the first problem; they cannot create information that the data do not contain.

The Weighted Design Matrix and Its Condition

For a univariate local polynomial of degree $p$, evaluate the fit at a point $x$ using bandwidth $h>0$. Define the local coordinate $u_i=(X_i-x)/h$, the polynomial basis $v_p(u)=(1,u,\ldots,u^p)^\mathsf{T}$, and nonnegative weights $w_i=K(u_i)$. Observations with zero weight contribute nothing. The weighted least-squares criterion is

\sum_{i=1}^n w_i\{Y_i-v_p(u_i)^\mathsf{T}\beta\}^2, \qquad A=\sum_{i=1}^n w_i v_p(u_i)v_p(u_i)^\mathsf{T}, \qquad b=\sum_{i=1}^n w_i v_p(u_i)Y_i.

When $A$ is invertible, the coefficient vector is $\hat\beta=A^{-1}b$, and the fitted value at $x$ is its intercept, $\hat m(x)=e_1^\mathsf{T}\hat\beta$, where $e_1=(1,0,\ldots,0)^\mathsf{T}$. Equivalently, let $Z$ have rows $\sqrt{w_i}v_p(u_i)^\mathsf{T}$; then $A=Z^\mathsf{T}Z$.

Definition: If $A$ is positive definite, its spectral condition number is $\kappa_2(A)=\lambda_{\max}(A)/\lambda_{\min}(A)$, where $\lambda_{\max}$ and $\lambda_{\min}$ are its largest and smallest eigenvalues. A large condition number means that some directions in the coefficient vector are much less informed by the weighted data than others. If $A$ is singular, its condition number is infinite and the weighted least-squares coefficients are not uniquely determined.

The weighted design has full column rank exactly when $A$ is positive definite. In one dimension, this requires at least $p+1$ distinct predictor values with positive weights. Indeed, if there are $p+1$ such values, the corresponding rows contain a square Vandermonde matrix with nonzero determinant; if there are fewer, the matrix has fewer than $p+1$ distinct rows and cannot have full column rank. Having enough distinct points is a rank condition, not a guarantee of good conditioning: those points may be extremely close together.

Centering and Rescaling the Polynomial Basis

The basis should reflect the local fitting problem. Using raw powers of $X_i$, such as $1,X_i,X_i^2$, can be numerically awkward when $x$ is far from zero. Centering replaces these columns with powers of $X_i-x$, or more commonly with powers of $u_i=(X_i-x)/h$. Centering puts the evaluation point at coordinate zero; rescaling expresses distances in bandwidth units.

Changing coordinates does not change the polynomial being fitted if the transformation is invertible. To see this in general, suppose the original basis vector is $\phi(u)$, and a new basis is $\psi(u)=T\phi(u)$ for an invertible matrix $T$. If the new coefficient vector is $\gamma$, then $\psi(u)^\mathsf{T}\gamma=\phi(u)^\mathsf{T}T^\mathsf{T}\gamma$. Thus the same fitted polynomial has $\beta=T^\mathsf{T}\gamma$. The least-squares criterion, which depends on fitted values and not on the labels assigned to coefficients, is unchanged. The matrix and coefficients may be better scaled in the new basis, but a nearly unsupported direction in the data remains nearly unsupported.

Worked Example: Centering Predictors Near a Large Evaluation Point

Fit a local linear polynomial at $x=1001$ to the three observations $(X_i,Y_i)=(1000,1),(1001,2),(1002,4)$, using equal positive weights. Compare the raw basis $(1,X_i)$ with the centered basis $(1,X_i-1001)$.

In the raw basis, the design cross-product is

A_{\mathrm{raw}}= \begin{pmatrix} 3&3003\\ 3003&3006005 \end{pmatrix}, \qquad \det(A_{\mathrm{raw}})=3(3006005)-3003^2=6.

The trace is $3006008$. Since the determinant is the product of the two positive eigenvalues and the larger eigenvalue is close to the trace, the smaller eigenvalue is approximately $6/3006008\approx 1.996\times 10^{-6}$. The larger eigenvalue is approximately $3006008$, giving $\kappa_2(A_{\mathrm{raw}})\approx 1.51\times 10^{12}$.

Now let $t_i=X_i-1001$, so the predictor values are $-1,0,1$. In the centered basis, $A_{\mathrm{cent}}=\begin{pmatrix}3&0\\0&2\end{pmatrix}$, with condition number $3/2$. The response vector in this basis is $b_{\mathrm{cent}}=(7,3)^\mathsf{T}$, so $\hat\beta=(7/3,3/2)^\mathsf{T}$. The fitted value at the evaluation point, where $t=0$, is $7/3$. Centering changes neither the local line nor its fitted value; it changes the numerical scale of the calculation substantially.

When Rescaling Helps, and When It Does Not

Even bandwidth-standardized coordinates can be tightly clustered if the positive-weight observations happen to lie very close to the evaluation point. In a polynomial basis, higher powers of small coordinates are especially small. This makes columns nearly redundant in the weighted design. Dividing the local coordinates by a representative local spread can improve their scale, provided the change is applied consistently to the basis and fitted value.

Worked Example: Rescaling a Tightly Clustered Quadratic Fit

Consider degree two, equal weights, and local coordinates $u=(-0.02,0,0.02)$, with responses $Y=(1,2,5)$. In the basis $(1,u,u^2)$, the cross-product matrix is

A_u= \begin{pmatrix} 3&0&0.0008\\ 0&0.0008&0\\ 0.0008&0&0.00000032 \end{pmatrix}.

The determinant of the block for the constant and quadratic terms is $3(0.00000032)-(0.0008)^2=0.00000032$. Its large eigenvalue is about $3$, and its small eigenvalue is about $0.00000032/3=1.067\times10^{-7}$. The linear-term eigenvalue is $0.0008$, so the full matrix has condition number about $2.81\times10^7$.

Set $t=u/0.02$, giving $t=(-1,0,1)$. In the basis $(1,t,t^2)$, the cross-product matrix is

A_t= \begin{pmatrix} 3&0&2\\ 0&2&0\\ 2&0&2 \end{pmatrix}.

Its eigenvalues are $2$ and $(5\pm\sqrt{17})/2$, so its condition number is $\{(5+\sqrt{17})/(5-\sqrt{17})\}\approx10.4$. The rescaling improves the conditioning while retaining the same quadratic polynomial: $2+2t+t^2=2+100u+2500u^2$. It does not change the fitted value at $u=0$, which is $2$. The three distinct locations make the fit identifiable, but the original coordinates made it unnecessarily difficult to compute.

The next distinction is essential: a singular design is not a scaling problem. A polynomial basis transformation that is invertible preserves the rank of the design. If the data do not identify all $p+1$ coefficients, rescaling cannot supply the missing direction.

Worked Example: Detecting Exact Rank Failure

Suppose a degree-two local fit has positive weights only at $u=-1$ and $u=1$, with responses $Y=1$ and $Y=3$. Its design rows are $(1,-1,1)$ and $(1,1,1)$. The coefficient vector $c=(-1,0,1)^\mathsf{T}$ is nonzero, but both rows multiplied by $c$ equal zero:

(1,-1,1)c=0,\qquad (1,1,1)c=0.

Consequently, $Zc=0$ and $Ac=Z^\mathsf{T}Zc=0$; $A$ is singular. In fact, the polynomial $u^2-1$ vanishes at both observed locations, so adding any multiple of it to a fitted polynomial leaves both fitted values unchanged. The data identify no unique quadratic coefficient vector. A broader neighborhood might include additional distinct locations, or one might fit a lower-degree polynomial if appropriate. Merely changing units or centering cannot resolve this exact rank failure.

Two Results for Diagnosing Numerical Sensitivity

Result: If $A$ is positive definite and a perturbed calculation uses $A+\Delta$ and $b+\delta b$, suppose $\eta=\|A^{-1}\Delta\|_2<1$. For $\hat\beta=A^{-1}b$ and $\tilde\beta=(A+\Delta)^{-1}(b+\delta b)$, $$ \|\tilde\beta-\hat\beta\|_2 \leq \frac{\|A^{-1}\|_2}{1-\eta} \left(\|\delta b\|_2+\|\Delta\|_2\|\hat\beta\|_2\right). $$

Proof. Since $\eta<1$, the matrix $I+A^{-1}\Delta$ is invertible. Its inverse has norm at most $1/(1-\eta)$: the Neumann series $\sum_{k=0}^{\infty}(-A^{-1}\Delta)^k$ converges in operator norm, and its norm is bounded by the geometric sum $\sum_{k=0}^{\infty}\eta^k=1/(1-\eta)$. Subtracting the equations for the two coefficient vectors gives

(I+A^{-1}\Delta)(\tilde\beta-\hat\beta) =A^{-1}(\delta b-\Delta\hat\beta).

Taking norms, applying the inverse bound and the triangle inequality yields the stated inequality. It shows why small perturbations can matter more when the matrix is poorly conditioned: $\|A^{-1}\|_2=1/\lambda_{\min}(A)$ is large when the least-informed direction has a small eigenvalue. $\square$

A related computational fact is that the condition number of the weighted design $Z$ is the square root of that of $A=Z^\mathsf{T}Z$, when $Z$ has full column rank. Indeed, the eigenvalues of $Z^\mathsf{T}Z$ are the squared singular values of $Z$, so $\kappa_2(A)=\kappa_2(Z)^2$. This is one reason to solve weighted least squares with a numerically stable QR or singular-value decomposition method rather than explicitly forming $A^{-1}$. Forming normal equations squares the condition number of the design.

A Practical Diagnostic Sequence

Check the positive-weight locations.
Count distinct predictor values with $w_i>0$. Fewer than $p+1$ means exact rank failure for a univariate degree-$p$ fit.

Inspect the coordinate scale.
Center at the evaluation point and express distances in local units. For multivariate fits, use the coordinate-specific bandwidth scaling established for the fit.

Assess rank and conditioning.
Use a QR or singular-value decomposition of the weighted design to inspect its numerical rank and smallest singular value. A large condition number signals sensitivity, not automatically a mathematical failure.

Choose a remedy consistent with the problem.
Rescale a poorly scaled basis. If observations do not support the requested degree, consider a justified change to the local data or polynomial degree rather than treating a coordinate transformation as a fix.

A common mistake is to see a large raw-basis condition number, rescale the columns, and conclude that the local fit is now well supported. Rescaling can reduce artificial numerical imbalance, but it cannot increase the information in the neighborhood. Conversely, treating every large condition number as proof of non-identifiability is also incorrect: the first worked example has a unique fit in both bases, although one representation is far less stable numerically. Diagnose the rank and the scale separately.

Key takeaway: Centering and rescaling are changes of coordinates that can stabilize a local polynomial calculation while preserving its fitted function. Rank checks answer a different question: whether the positively weighted observations identify the requested polynomial at all.

Check Your Understanding

Use the distinction between numerical scaling and local information to answer these questions.

For a univariate degree-three local polynomial, what condition on the positive-weight predictor values is necessary for full column rank?
Why can centering at the evaluation point improve conditioning when predictors are near a large constant?
In the clustered quadratic example, why did rescaling improve conditioning without changing the fitted value at the evaluation point?
Can an invertible change of polynomial basis repair an exactly singular weighted design? Explain.
Why can solving weighted least squares through the normal equations be less numerically stable than using QR or singular-value decomposition?