# Backfitting algorithm

In statistics, the backfitting algorithm is a simple iterative procedure used to fit a generalized additive model. It was introduced in 1985 by Leo Breiman and Jerome Friedman along with generalized additive models. In most cases, the backfitting algorithm is equivalent to the Gauss–Seidel method algorithm for solving a certain linear system of equations

## 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:

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]

## Motivation

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.

## Explicit derivation for two dimensions

For the two dimensional case, we can formulate the backfitting algorithm explicitly. 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-1)}]}$

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 assume our constant ${\displaystyle \alpha }$ is zero and we set ${\displaystyle {\hat {f}}_{2}^{(0)}=0}$ then we get

${\displaystyle {\hat {f}}_{1}^{(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}$

This converges if ${\displaystyle \|S_{1}S_{2}\|<1}$.

## Issues

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.

## Modified algorithm

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 Si 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}}}}$


## References

1. Hastie, Trevor, Robert Tibshirani and Jerome Friedman (2001). The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer, ISBN 0-387-95284-5.
• {{#invoke:Citation/CS1|citation

|CitationClass=journal }}

• {{#invoke:Citation/CS1|citation

|CitationClass=journal }}