Backfitting algorithm

Additive models are a class of non-parametric regression models of the form:

${\displaystyle Y_{i}=\alpha +\sum _{j=1}^{p}f_{j}(X_{ij})+\epsilon _{i}}$

where each ${\displaystyle X_{1},X_{2},\ldots ,X_{p}}$ is a variable in our ${\displaystyle p}$-dimensional predictor ${\displaystyle X}$, and ${\displaystyle Y}$ is our outcome variable. ${\displaystyle \epsilon }$ represents our inherent error, which is assumed to have mean zero. The ${\displaystyle f_{j}}$ represent unspecified smooth functions of a single ${\displaystyle X_{j}}$. Given the flexibility in the ${\displaystyle f_{j}}$, we typically do not have a unique solution: ${\displaystyle \alpha }$ is left unidentifiable as one can add any constants to any of the ${\displaystyle f_{j}}$ and subtract this value from ${\displaystyle \alpha }$. It is common to rectify this by constraining

${\displaystyle \sum _{i=1}^{N}f_{j}(X_{ij})=0}$ for all ${\displaystyle j}$

leaving

${\displaystyle \alpha =1/N\sum _{i=1}^{N}y_{i}}$

necessarily.

The backfitting algorithm is then:

   Initialize ${\displaystyle {\hat {\alpha }}=1/N\sum _{i=1}^{N}y_{i},{\hat {f_{j}}}\equiv 0}$,${\displaystyle \forall j}$
   Do until ${\displaystyle {\hat {f_{j}}}}$ converge:
       For each predictor j:
           (a) ${\displaystyle {\hat {f_{j}}}\leftarrow {\text{Smooth}}[\lbrace y_{i}-{\hat {\alpha }}-\sum _{k\neq j}{\hat {f_{k}}}(x_{ik})\rbrace _{1}^{N}]}$ (backfitting step)
           (b) ${\displaystyle {\hat {f_{j}}}\leftarrow {\hat {f_{j}}}-1/N\sum _{i=1}^{N}{\hat {f_{j}}}(x_{ij})}$ (mean centering of estimated function)

where ${\displaystyle {\text{Smooth}}}$ is our smoothing operator. This is typically chosen to be a cubic spline smoother but can be any other appropriate fitting operation, such as:

• local polynomial regression
• kernel smoothing methods
• more complex operators, such as surface smoothers for second and higher-order interactions

In theory, step (b) in the algorithm is not needed as the function estimates are constrained to sum to zero. However, due to numerical issues this might become a problem in practice.[1]

If we consider the problem of minimizing the expected squared error:

${\displaystyle \min E[Y-(\alpha +\sum _{j=1}^{p}f_{j}(X_{j}))]^{2}}$

There exists a unique solution by the theory of projections given by:

${\displaystyle f_{i}(X_{i})=E[Y-(\alpha +\sum _{j\neq i}^{p}f_{j}(X_{j}))|X_{i}]}$

for i = 1, 2, ..., p.

This gives the matrix interpretation:

${\displaystyle {\begin{pmatrix}I&P_{1}&\cdots &P_{1}\\P_{2}&I&\cdots &P_{2}\\\vdots &&\ddots &\vdots \\P_{p}&\cdots &P_{p}&I\end{pmatrix}}{\begin{pmatrix}f_{1}(X_{1})\\f_{2}(X_{2})\\\vdots \\f_{p}(X_{p})\end{pmatrix}}={\begin{pmatrix}P_{1}Y\\P_{2}Y\\\vdots \\P_{p}Y\end{pmatrix}}}$

where ${\displaystyle P_{i}(\cdot )=E(\cdot |X_{i})}$. In this context we can imagine a smoother matrix, ${\displaystyle S_{i}}$, which approximates our ${\displaystyle P_{i}}$ and gives an estimate, ${\displaystyle S_{i}Y}$, of ${\displaystyle E(Y|X)}$

${\displaystyle {\begin{pmatrix}I&S_{1}&\cdots &S_{1}\\S_{2}&I&\cdots &S_{2}\\\vdots &&\ddots &\vdots \\S_{p}&\cdots &S_{p}&I\end{pmatrix}}{\begin{pmatrix}f_{1}\\f_{2}\\\vdots \\f_{p}\end{pmatrix}}={\begin{pmatrix}S_{1}Y\\S_{2}Y\\\vdots \\S_{p}Y\end{pmatrix}}}$

or in abbreviated form

${\displaystyle {\hat {S}}f=QY\,}$

An exact solution of this is infeasible to calculate for large np, so the iterative technique of backfitting is used. We take initial guesses ${\displaystyle f_{i}^{(0)}}$ and update each ${\displaystyle f_{i}^{(j)}}$ in turn to be the smoothed fit for the residuals of all the others:

${\displaystyle {\hat {f_{i}}}^{(j)}\leftarrow {\text{Smooth}}[\lbrace y_{i}-{\hat {\alpha }}-\sum _{k\neq j}{\hat {f_{k}}}(x_{ik})\rbrace _{1}^{N}]}$

Looking at the abbreviated form it is easy to see the backfitting algorithm as equivalent to the Gauss–Seidel method for linear smoothing operators S.

Following,[2] we can formulate the backfitting algorithm explicitly for the two dimensional case. We have:

${\displaystyle f_{1}=S_{1}(Y-f_{2}),f_{2}=S_{2}(Y-f_{1})}$

If we denote ${\displaystyle {\hat {f}}_{1}^{(i)}}$ as the estimate of ${\displaystyle f_{1}}$ in the ith updating step, the backfitting steps are

${\displaystyle {\hat {f}}_{1}^{(i)}=S_{1}[Y-{\hat {f}}_{2}^{(i-1)}],{\hat {f}}_{2}^{(i)}=S_{2}[Y-{\hat {f}}_{1}^{(i)}]}$

By induction we get

${\displaystyle {\hat {f}}_{1}^{(i)}=Y-\sum _{\alpha =0}^{i-1}(S_{1}S_{2})^{\alpha }(I-S_{1})Y-(S_{1}S_{2})^{i-1}S_{1}{\hat {f}}_{2}^{(0)}}$

and

${\displaystyle {\hat {f}}_{2}^{(i)}=S_{2}\sum _{\alpha =0}^{i-1}(S_{1}S_{2})^{\alpha }(I-S_{1})Y+S_{2}(S_{1}S_{2})^{i-1}S_{1}{\hat {f}}_{2}^{(0)}}$

If we set ${\displaystyle {\hat {f}}_{2}^{(0)}=0}$ then we get

${\displaystyle {\hat {f}}_{1}^{(i)}=Y-S_{2}^{-1}{\hat {f}}_{2}^{(i)}=[I-\sum _{\alpha =0}^{i-1}(S_{1}S_{2})^{\alpha }(I-S_{1})]Y}$

${\displaystyle {\hat {f}}_{2}^{(i)}=[S_{2}\sum _{\alpha =0}^{i-1}(S_{1}S_{2})^{\alpha }(I-S_{1})]Y}$

Where we have solved for ${\displaystyle {\hat {f}}_{1}^{(i)}}$ by directly plugging out from ${\displaystyle f_{2}=S_{2}(Y-f_{1})}$.

We have convergence if ${\displaystyle \|S_{1}S_{2}\|<1}$. In this case, letting ${\displaystyle {\hat {f}}_{1}^{(i)},{\hat {f}}_{2}^{(i)}{\xrightarrow {}}{\hat {f}}_{1}^{(\infty )},{\hat {f}}_{2}^{(\infty )}}$:

${\displaystyle {\hat {f}}_{1}^{(\infty )}=Y-S_{2}^{-1}{\hat {f}}_{2}^{(\infty )}=Y-(I-S_{1}S_{2})^{-1}(I-S_{1})Y}$

${\displaystyle {\hat {f}}_{2}^{(\infty )}=S_{2}(I-S_{1}S_{2})^{-1}(I-S_{1})Y}$

We can check this is a solution to the problem, i.e. that ${\displaystyle {\hat {f}}_{1}^{(i)}}$ and ${\displaystyle {\hat {f}}_{2}^{(i)}}$ converge to ${\displaystyle f_{1}}$ and ${\displaystyle f_{2}}$ correspondingly, by plugging these expressions into the original equations.

The choice of when to stop the algorithm is arbitrary and it is hard to know a priori how long reaching a specific convergence threshold will take. Also, the final model depends on the order in which the predictor variables ${\displaystyle X_{i}}$ are fit.

As well, the solution found by the backfitting procedure is non-unique. If ${\displaystyle b}$ is a vector such that ${\displaystyle {\hat {S}}b=0}$ from above, then if ${\displaystyle {\hat {f}}}$ is a solution then so is ${\displaystyle {\hat {f}}+\alpha b}$ is also a solution for any ${\displaystyle \alpha \in \mathbb {R} }$. A modification of the backfitting algorithm involving projections onto the eigenspace of S can remedy this problem.

We can modify the backfitting algorithm to make it easier to provide a unique solution. Let ${\displaystyle {\mathcal {V}}_{1}(S_{i})}$ be the space spanned by all the eigenvectors of S_i that correspond to eigenvalue 1. Then any b satisfying ${\displaystyle {\hat {S}}b=0}$ has ${\displaystyle b_{i}\in {\mathcal {V}}_{1}(S_{i})\forall i=1,\dots ,p}$ and ${\displaystyle \sum _{i=1}^{p}b_{i}=0.}$ Now if we take ${\displaystyle A}$ to be a matrix that projects orthogonally onto ${\displaystyle {\mathcal {V}}_{1}(S_{1})+\dots +{\mathcal {V}}_{1}(S_{p})}$, we get the following modified backfitting algorithm:

   Initialize ${\displaystyle {\hat {\alpha }}=1/N\sum _{1}^{N}y_{i},{\hat {f_{j}}}\equiv 0}$,${\displaystyle \forall i,j}$, ${\displaystyle {\hat {f_{+}}}=\alpha +{\hat {f_{1}}}+\dots +{\hat {f_{p}}}}$
   Do until ${\displaystyle {\hat {f_{j}}}}$ converge:
       Regress ${\displaystyle y-{\hat {f_{+}}}}$ onto the space ${\displaystyle {\mathcal {V}}_{1}(S_{i})+\dots +{\mathcal {V}}_{1}(S_{p})}$, setting ${\displaystyle a=A(Y-{\hat {f_{+}}})}$
       For each predictor j:
           Apply backfitting update to ${\displaystyle (Y-a)}$ using the smoothing operator ${\displaystyle (I-A_{i})S_{i}}$, yielding new estimates for ${\displaystyle {\hat {f_{j}}}}$

• Breiman, L. & Friedman, J. H. (1985). "Estimating optimal transformations for multiple regression and correlations (with discussion)". Journal of the American Statistical Association. 80 (391): 580–619. doi:10.2307/2288473. JSTOR 2288473.
• Hastie, T. J. & Tibshirani, R. J. (1990). "Generalized Additive Models". Monographs on Statistics and Applied Probability. 43.
• Härdle, Wolfgang; et al. (June 9, 2004). "Backfitting". Archived from the original on 2015-05-10. Retrieved 2015-08-19.

• R Package for GAM backfitting
• R Package for BRUTO backfitting