The Annals of Statistics

2001, Vol. 29, No. 5, 1189–1232

1999 REITZ LECTURE

GREEDY FUNCTION APPROXIMATION:
A GRADIENT BOOSTING MACHINE1

By Jerome H. Friedman
Stanford University
Function estimation/approximation is viewed from the perspective of
numerical optimization in function space, rather than parameter space. A
connection is made between stagewise additive expansions and steepest descent
minimization. A general gradient descent “boosting” paradigm is developed for
additive expansions based on any fitting criterion. Specific algorithms are
presented for least-squares, least absolute deviation, and Huber-M loss functions
for regression, and multiclass logistic likelihood for classification. Special
enhancements are derived for the particular case where the individual additive
components are regression trees, and tools for interpreting such “TreeBoost”
models are presented. Gradient boost ing ofregression trees produces
competitive, highly robust, interpretable procedures for both regression and
classification, especially appropriate for mining less than clean data. Connections
between this approach and the boosting methods ofFreund and Shapire and
Friedman, Hastie and Tib shirani are discussed.

1. Function estimation. In the function estimation or “predictive learn ing” problem,

one has a system consisting ofa random “output” or “response” variable y and a set
ofrandom “input” or “explanatory” variables x = x1 xn . Using a “training” sample yi
xi N1 ofknown y x -values, the goal is to obtain an estimate or approximation F x ,
ofthe function F∗ x mapping x to y, that minimizes the expected value ofsome
specified loss function L y F x over the joint distribution ofall y x -values,

F∗ = arg minFEy xL y F x = arg minFEx Ey (1) L y F x x

Frequently employed loss functions L y F include squared-error y − F 2

and absolute error y − F for y ∈ R 1 (regression) and negative binomial log
likelihood, log 1 + e−2yF , when y ∈ −1 1 (classification).
A common procedure is to restrict F x to be a member ofa parameterized class
offunctions F x P , where P = P 1 P2 is a finite set ofparameters whose joint values
identify individual class members. In this article we focus

Received May 1999; revised April 2001.

1
Supported in part by CSIRO Mathematical and Information Science, Australia; Department ofEnergy
Contract DE-AC03-76SF00515; and NSF Grant DMS-97-64431. AMS 2000 subject classifications. 62-02,
62-07, 62-08, 62G08, 62H30, 68T10. Key words and phrases. Function estimation, boosting, decision trees,
robust nonparametric regression.
1189
1190 J. H. FRIEDMAN

on “additive” expansions ofthe form

F x βm am M1 = M

(2) βmh x am
m=1
 The (generic) function h x a in (2) is usually a simple parameterized func tion
ofthe input variables x, characterized by parameters a = a1 a2 . The individual
terms differ in the joint values am chosen for these parameters. Such expansions (2)
are at the heart ofmany function approximation meth ods such as neural networks
[Rumelhart, Hinton, and Williams (1986)], radial basis functions [Powell (1987)],
MARS [Friedman (1991)], wavelets [Donoho (1993)] and support vector machines
[Vapnik (1995)]. Ofspecial interest here is the case where each ofthe functions h x a m
is a small regression tree, such as those produced by CART TM [Breiman, Friedman,
Olshen and Stone (1983)]. For a regression tree the parameters a m are the splitting
variables, split locations and the terminal node means ofthe individual trees.

1.1. Numerical optimization. In general, choosing a parameterized model F x P

changes the function optimization problem to one of parameter opti mization,
P∗ = arg minP
(3) P

where

P = Ey xL y F x P
and then
F∗ x = F x P ∗

For most F x P and L, numerical optimization methods must be applied to solve (3).
This often involves expressing the solution for the parameters in the form

P∗ = M

(4) pm
m=0

where p0 is an initial guess and p m M1are successive increments (“steps” or

“boosts”), each based on the sequence ofpreceding steps. The prescription for
computing each step pm is defined by the optimization method.

1.2. Steepest-descent. Steepest-descent is one ofthe simplest ofthe frequ ently

used numerical minimization methods. It defines the increments pm M1 (4) as follows.
First the current gradient gm is computed:

P=Pm−1
gm = gjm = ∂ P ∂Pj

GREEDY FUNCTION APPROXIMATION 1191

where

pi The step is taken to be Pm−1 =

m −1
i=0

pm =−ρmgm
where

ρm = arg minρ (5) Pm−1 − ρgm

The negative gradient −gm is said to define the “steepest-descent” direction and
(5) is called the “line search” along that direction.

2. Numerical optimization in function space. Here we take a “non parametric”

approach and apply numerical optimization in function space. That is, we consider F
x evaluated at each point x to be a “parameter” and seek to minimize

F = Ey xL y F x = Ex Ey L y F x x
or equivalently,

φ F x = Ey L y F x x
at each individual x, directly with respect to F x . In function space there are an
infinite number ofsuch parameters, but in data sets (discussed below) only a finite
number F xi N1 are involved. Following the numerical optimization paradigm we take
the solution to be

F∗ x = M
m=0 fm x

where f0 x is an initial guess, and fm x M1 are incremental functions (“steps” or

“boosts”) defined by the optimization method.
For steepest-descent,

(6) fm x =−ρmgm x
with
and ∂φ F x ∂F x F x =Fm−1 x = x ∂F x
gm x = ∂Ey L y F x
F x =Fm−1 x

m −1
Fm−1 x = i=0 fi x
1192 J. H. FRIEDMAN

Assuming sufficient regularity that one can interchange differentiation and

integration, this becomes

gm x = E y ∂L y F x

(7)
F x =Fm−1 x
∂F x x
The multiplier ρm in (6) is given by the line search
ρm = arg minρ Ey xL y Fm−1 (8) x − ρgm x
 3. Finitedata. This nonparametric approach breaks down when the joint
distribution of y x is estimated by a finite data sample y i xi N1 . In this case Ey · x
cannot be estimated accurately by its data value at each x i, and even ifit could, one
would like to estimate F∗ x at x values other than the training sample points.
Strength must be borrowed from nearby data points by imposing smoothness on the
solution. One way to do this is to assume a parameterized form such as (2) and do
parameter optimization as discussed in Section 1.1 to minimize the corresponding
data based estimate ofexpected loss,

m=1
βm am M1 =
M M
1 L yi
arg min β m a m N
β mh xi a m
i=1

In situations where this is infeasible one can try a “greedy stagewise” approach. For
m = 1 2 M,
βm am = arg min β a N
L yi Fm−1 xi + βh xi (9) a i=1

and then

Fm x = Fm−1 (10) x + βmh x am

Note that this stagewise strategy is different from stepwise approaches that readjust
previously entered terms when new ones are added. In signal processing this
stagewise strategy is called “matching pursuit” [Mallat and Zhang (1993)] where L y
F is squared-error loss and the h x am M1 are called basis functions, usually
taken from an overcomplete waveletlike dictionary. In machine learning,
(9), (10) is called “boosting” where y ∈ −1 1 and L y F is either an
exponential loss criterion e−yF [Freund and Schapire (1996), Schapire and Singer
(1998)] or negative binomial log likelihood [Friedman, Hastie and Tibshirani (2000)
(here after reffered to as FHT00)]. The function h x a is called a “weak learner” or
“base learner” and is usually a classification tree.
Suppose that for a particular loss L y F and/or base learner h x a the solution to
(9) is difficult to obtain. Given any approximator F m−1 x , the function βmh x am (9),
(10) can be viewed as the best greedy step toward the data-based estimate of F ∗ x
(1), under the constraint that the step “direc tion” h x a m be a member ofthe
parameterized class offunctions h x a . It can thus be regarded as a steepest descent
step (6) under that constraint. By
GREEDY FUNCTION APPROXIMATION 1193 construction, the data-

based analogue ofthe unconstrained negative gradi

ent (7), −gm xi =−

∂L yi F xi ∂F xi
F x =Fm−1 x
in the N
gives the best steepest-descent step direction −g m = −gm xi N1
dimensional data space at Fm−1 x . However, this gradient is defined only at the data
and cannot be generalized to other x-values. One possi
points xi N1 bility for
generalization is to choose that member of the parameterized class h x a m that
produces hm = h xi am N1 most parallel to −gm ∈ RN. This is the h x a most
highly correlated with −gm x over the data distribution. It can be obtained from
the solution
am = arg min a β N
−gm xi − βh xi a 2 (11) i=1

is used in place ofthe uncon

This constrained negative gradient h x am strained
one −gm x (7) in the steepest-descent strategy. Specifically, the line search (8) is
performed

ρm = arg minρ N
L yi Fm−1 xi + ρh xi (12) am i=1

and the approximation updated,

Fm x = Fm−1 x + ρmh x am
Basically, instead ofobtaining the solution under a smoothness constraint (9), the
constraint is applied to the unconstrained (rough) solution by fit ting h x a to the
“pseudoresponses” ˜yi =−gm xi Ni=1 (7). This permits the replacement of the difficult
function minimization problem (9) by least-squares function minimization (11),
followed by only a single parameter optimization based on the original criterion (12).
Thus, for any h x a for which a feasible least-squares algorithm exists for solving
(11), one can use this approach to minimize any differentiable loss L y F in
conjunction with forward stage wise additive modeling. This leads to the following
(generic) algorithm using steepest-descent.

Algorithm 1 (Gradient Boost).

1. F0 x = arg minρ
N
i=1 L yi ρ
2. For m = 1 to M do:
∂F xi
3. y˜i =− ∂L yi F xi
F x =Fm−1 x i=1N
4. am = arg mina β
N 2
i=1 ˜yi − βh xi a 5. ρm = arg minρ
N
i=1 L yi Fm−1 xi + ρh xi am 6. Fm x = Fm−1 x + ρmh x am
1194 J. H. FRIEDMAN

7. endFor
end Algorithm

Note that any fitting criterion that estimates conditional expectation (given x) could
in principle be used to estimate the (smoothed) negative gradient (7) at line 4
ofAlgorithm 1. Least-squares (11) is a natural choice owing to the superior
computational properties ofmany least-squares algorithms.
In the special case where y ∈ −1 1 and the loss function L y F
depends on y and F only through their product L y F = L yF , the analogy
ofboost ing (9), (10) to steepest-descent minimization has been noted in
the machine learning literature [Ratsch, Onoda and Muller (1998), Breiman
(1999)]. Duffy and Helmbold (1999) elegantly exploit this analogy to
motivate their GeoLev and GeoArc procedures. The quantity yF is called
the “margin” and the steepest-descent is performed in the space of
margin values, rather than the space offunction values F. The latter
approach permits application to more general loss functions where the
notion of margins is not apparent. Drucker (1997) employs a different
strategy of casting regression into the framework of classification in the
context ofthe AdaBoost algorithm [Freund and Schapire (1996)].

4. Applications: additivemodeling. In this section the gradient boost ing strategy is

applied to several popular loss criteria: least-squares (LS), least absolute deviation
(LAD), Huber (M), and logistic binomial log-likelihood (L). The first serves as a “reality
check”, whereas the others lead to new boosting algorithms.

4.1. Least-squares regression. Here L y F = y − F 2/2. The pseudore sponse in

line 3 ofAlgorithm 1 is y˜i = yi − Fm−1 xi . Thus, line 4 simply fits the current residuals
and the line search (line 5) produces the result ρm = βm, where βm is the minimizing β
ofline 4. Therefore, gradient boosting on squared-error loss produces the usual
stagewise approach ofiteratively fitting the current residuals.

Algorithm 2 (LS Boost).

F0 x = ¯y
For m = 1 to M do:
y˜i = yi − Fm−1 xi i = 1 N
ρm am = arg mina ρ
N 2
i=1 ˜yi − ρh xi a
Fm x = Fm−1 x + ρmh x am
endFor
end Algorithm

4.2. Least absolute deviation (LAD) regression. For the loss function
L y F = y − F , one has
y˜i =− ∂L yi F xi
= sign yi − Fm−1 xi (13)
F x =Fm−1 x
∂F xi
GREEDY FUNCTION APPROXIMATION 1195

This implies that h x a is fit (by least-squares) to the sign ofthe current residuals in
line 4 ofAlgorithm 1. The line search (line 5) becomes

ρm = arg minρ N i=1 yi − Fm−1 xi − ρh xi am

(14) = arg minρ N i=1 yi − Fm−1 xi

h xi a m · h x i am − ρ

= medianW yi − Fm−1 xi h xi am N1 wi = h x i a m

Here medianW · is the weighted median with weights wi. Inserting these results [(13),
(14)] into Algorithm 1 yields an algorithm for least absolute deviation boosting, using
any base learner h x a .
 4.3. Regression trees. Here we consider the special case where each base
learner is an J-terminal node regression tree [Breiman, Friedman, Olshen and Stone
(1983)]. Each regression tree model itselfhas the additive form

h x bj Rj J1 = J

bj1 x ∈ Rj (15)
j=1

Here Rj J1 are disjoint regions that collectively cover the space ofall joint values ofthe
predictor variables x. These regions are represented by the ter minal nodes ofthe
corresponding tree. The indicator function 1 · has the value 1 ifits argument is true,
and zero otherwise. The “parameters” ofthis base learner (15) are the coefficients b j
J J
1 , and the quantities that define the boundaries ofthe regions R j 1 . These are the
splitting variables and the values ofthose variables that represent the
splits at the nonterminal nodes of the tree. Because the regions are
disjoint, (15) is equivalent to the prediction rule: if x ∈ R j then h x = bj.
For a regression tree, the update at line 6 ofAlgorithm 1 becomes Fm x =

Fm−1 x + ρm J

(16) bjm1 x ∈ Rjm j=1

Here Rjm J1 are the regions defined by the terminal nodes ofthe tree at the mth
iteration. They are constructed to predict the pseudoresponses ˜y i N1 (line 3) by
least-squares (line 4). The bjm are the corresponding least-squares coefficients,

bjm = avexi∈Rjm y˜i

The scaling factor ρm is the solution to the “line search” at line 5.

1196 J. H. FRIEDMAN

The update (16) can be alternatively expressed as

J
Fm x = Fm−1 x +
(17) γjm1 x ∈ Rjm
j=1

with γjm = ρmbjm. One can view (17) as adding J separate basis functions at
each step 1 x ∈ Rjm J1 , instead ofa single additive one as in (16). Thus, in this
case one can further improve the quality of the fit by using the opti mal coefficients
for each of these separate basis functions (17). These optimal coefficients are the
solution to

γjm J1 = arg N
i=1 yi Fm−1 xi + J γj1 x ∈ Rjm
min γj J1 L j=1

Owing to the disjoint nature ofthe regions produced by regression trees, this reduces
to

γjm = arg min

L yi Fm−1 xi (18) + γ
γ xi∈Rjm

This is just the optimal constant update in each terminal node region, based on the
loss function L, given the current approximation Fm−1 x . For the case ofLAD
regression (18) becomes

γjm = medianxi∈Rjm yi − Fm−1 xi

which is simply the median ofthe current residuals in the jth terminal node at the mth
iteration. At each iteration a regression tree is built to best predict the sign ofthe
current residuals yi −Fm−1 xi , based on a least-squares criterion. Then the
approximation is updated by adding the median ofthe residuals in each ofthe derived
terminal nodes.

Algorithm 3 (LAD TreeBoost).

F0 x = median yi N1
For m = 1 to M do:
y˜i = sign yi − Fm−1 xi i = 1 N
Rjm J1 = J-terminal node tree ˜yi xi N1
γjm = medianxi∈Rjm yi − Fm−1 xi j = 1 J
Fm x = Fm−1 x +
J
j=1 γjm1 x ∈ Rjm
endFor
end Algorithm

This algorithm is highly robust. The trees use only order information on the
individual input variables xj, and the pseudoresponses y˜i (13) have only two values,
y˜i ∈ −1 1 . The terminal node updates are based on medians.
GREEDY FUNCTION APPROXIMATION 1197

An alternative approach would be to build a tree to directly minimize the loss

criterion,
i=1
treem x = arg
N yi − Fm−1 xi − tree xi
min J-node tree

Fm x = Fm−1 xi + treem x
However, Algorithm 3 is much faster since it uses least-squares to induce the trees.
Squared-error loss is much more rapidly updated than mean absolute deviation when
searching for splits during the tree building process.

4.4. M-Regression. M-regression techniques attempt resistance to long tailed

error distributions and outliers while maintaining high efficiency for normally
distributed errors. We consider the Huber loss function [Huber (1964)]

Ly F= 1 2
2 y−F y−F≤δ
δ y − F − δ/2 y − F >δ. (19) Here the pseudoresponse is

y˜i =− ∂L yi F xi ∂F xi F x =Fm−1 x

yi − Fm−1 xi yi − Fm−1 xi ≤ δ, =
δ · sign yi − Fm−1 xi yi − Fm−1 xi >δ,
and the line search becomes
ρm = arg minρ N
L yi Fm−1 xi + ρh xi (20) am i=1

with L given by (19). The solution to (19), (20) can be obtained by standard iterative
methods [see Huber (1964)].
The value ofthe transition point δ defines those residual values that are
considered to be “outliers,” subject to absolute rather than squared-error
loss. An optimal value will depend on the distribution of y − F ∗ x , where F∗
is the true target function (1). A common practice is to choose the value of
δ to be the α-quantile ofthe distribution of y − F ∗ x , where 1 − α
controls the breakdown point ofthe procedure. The “breakdown point” is
the fraction of observations that can be arbitrarily modified without
seriously degrading the quality ofthe result. Since F ∗ x is unknown one uses
the current estimate Fm−1 x as an approximation at the mth iteration. The
distribution of y − Fm−1 x is estimated by the current residuals, leading to
N
δm = quantileα yi − Fm−1 xi 1

With regression trees as base learners we use the strategy ofSection 4.3, that is,
a separate update (18) in each terminal node Rjm. For the Huber loss
1198 J. H. FRIEDMAN

(19) the solution to (18) can be approximated by a single step ofthe standard iterative
procedure [Huber (1964)] starting at the median

r˜jm = medianxi∈Rjm rm−1 xi

N
where rm−1 xi 1 are the current residuals
rm−1 xi = yi − Fm−1 xi
The approximation is

γjm =˜rjm +1 Njm sign rm−1 xi − ˜rjm ·min δm

xi∈Rjm
abs rm−1 xi − ˜rjm

where Njm is the number ofobservations in the jth terminal node. This gives the
following algorithm for boosting regression trees based on Huber loss (19).

Algorithm 4 M TreeBoost.
F0 x = median yi N1
For m = 1 to M do:
rm−1 xi = yi − Fm−1 xi i = 1 N
 N
δm =quantileα rm−1 xi 1

y˜i = rm−1 xi >δm i = 1 N

rm−1 xi rm−1 xi ≤ δm δm · sign rm−1 xi
Rjm J1 = J-terminal node tree ˜yi xi N1 r˜jm =
medianxi∈Rjm rm−1 xi , j = 1 J
γjm =˜rjm + 1
Njm
j=1J

xi∈Rjm signrm−1 xi − ˜rjm · min δm abs rm−1

xi − ˜rjm ,
 Fm x = Fm−1 x +
J
j=1 γjm1 x ∈ Rjm
endFor
end Algorithm

According to the motivations underlying robust regression, this algorithm should

have properties similar to that ofleast-squares boosting (Algorithm 2) for normally
distributed errors, and similar to that of least absolute deviation regression (Algorithm
3) with very long-tailed distributions. For error distri butions with only moderately long
tails it can have performance superior to both (see Section 6.2).

4.5. Two-class logistic regression and classification. Here the loss function is
negative binomial log-likelihood (FHT00)
L y F = log 1 + exp −2yF y ∈ −1 1
where 1
F x = 2 log Pr y = 1 x

(21) Pr y =−1 x
GREEDY FUNCTION APPROXIMATION 1199 The pseudoresponse is

∂L yi F xi
= 2yi/ 1 + exp 2yiFm−1 xi (22)
F x =Fm−1 x
y˜i =− ∂F xi
i=1
The line search becomes ρm = arg minρ N
log 1 + exp −2yi Fm−1 xi + ρh xi am

With regression trees as base learners we again use the strategy (Section 4.3)
ofseparate updates in each terminal node Rjm:

log 1 + exp −2yi Fm−1 xi (23) + γ γjm = arg min

γ xi∈Rjm

There is no closed-form solution to (23). Following FHT00, we approximate it by a

single Newton–Raphson step. This turns out to be
xi∈Rjm
γjm = xi∈Rjm y˜i ˜yi 2 − ˜yi

with y˜i given by (22). This gives the following algorithm for likelihood gradient
boosting with regression trees.

Algorithm 5 (LK TreeBoost).

F0 x = 12 log 1+¯y
1−¯y
For m = 1 to M do:
y˜i = 2yi/ 1 + exp 2yiFm−1 xi , i = 1 N
 Rjm J1 = J-terminal node tree ˜yi xi N1
γjm =

xi∈Rjm y˜i/

xi∈Rjm ˜yi 2 − ˜yi , j = 1 J

Fm x = Fm−1 x +
J
j=1 γjm1 x ∈ Rjm
endFor
end Algorithm

The final approximation FM x is related to log-odds through (21). This can be

inverted to yield probability estimates

p+ x = Pr y = 1 x = 1/ 1 + e−2FM x

p− x = Pr y =−1 x = 1/ 1 + e2FM x
These in turn can be used for classification,

yˆ x = 2 · 1 c −1 1 p+ x >c 1 −1 p− x − 1
where c ˆy y is the cost associated with predicting yˆ when the truth is y.
1200 J. H. FRIEDMAN

4.5.1. Influence trimming. The empirical loss function for the two-class logistic
regression problem at the mth iteration is

φm ρ a = N
log 1 + exp −2yiFm−1 xi · exp −2yiρh xi (24) a i=1

If yiFm−1 xi is very large, then (24) has almost no dependence on ρh x i a for small to
moderate values near zero. This implies that the ith observation y i xi has almost no
influence on the loss function, and therefore on its solution

ρm am = arg min ρ a φm ρ a

This suggests that all observations y i xi for which yiFm−1 xi is relatively very large
can be deleted from all computations of the mth iteration without having a substantial
effect on the result. Thus,

wi = exp −2yiFm−1 xi (25)

can be viewed as a measure ofthe “influence” or weight ofthe ith observation on the
estimate ρmh x am .
More generally, from the nonparametric function space perspective of Section 2,
the parameters are the observation function values F x i N1 . The influence on an
estimate to changes in a “parameter” value F x i (holding all the other parameters
fixed) can be gauged by the second derivative ofthe loss function with respect to that
parameter. Here this second derivative at the mth iteration is ˜y i 2 − ˜yi with y˜i
given by (22). Thus, another measure of the influence or “weight” ofthe ith
observation on the estimate ρmh x am at the mth iteration is

wi = ˜yi 2 − ˜yi (26)

Influence trimming deletes all observations with wi-values less than wl α , where l α
is the solution to

l α
N
wi = α
wi (27)
i=1 i=1

Here w i N1 are the weights wi N1 arranged in ascending order. Typical values

are α ∈ 0 05 0 2 . Note that influence trimming based on (25), (27) is
identical to the “weight trimming” strategy employed with Real AdaBoost,
whereas (26), (27) is equivalent to that used with LogitBoost, in FHT00.
There it was seen that 90% to 95% ofthe observations were often deleted
without sacrificing accuracy ofthe estimates, using either influence
measure. This results in a corresponding reduction in computation by
factors of 10 to 20.
GREEDY FUNCTION APPROXIMATION 1201

4.6. Multiclass logistic regression and classification. Here we develop a gradient-

descent boosting algorithm for the K-class problem. The loss func tion is
K K
L yk F k x 1 =−
yk log pk (28) x
k=1

where yk = 1(class = k) ∈ 0 1 , and pk x = Pr yk = 1 x . Following FHT00, we

use the symmetric multiple logistic transform
1
Fk x = log pk x − K K
log pl (29) x
l=1

or equivalently

K
pk x = exp Fk x
exp Fl (30) x l=1

Substituting (30) into (28) and taking first derivatives one has
∂L yil Fl xi Kl=1

= yik − pk m−1 xi (31)

K
Fl x =Fl x
y˜ik =− ∂Fk xi m−1 1

where pk m−1 x is derived from Fk m−1 x through (30). Thus, K-trees are induced at
each iteration m to predict the corresponding current residuals for each class on the
probability scale. Each ofthese trees has J-terminal nodes, with corresponding
regions Rjkm Jj=1. The model updates γjkm corresponding to these regions are the
solution to
N K φ
γjkm = arg min γjk i=1 k=1
 J
yik Fk m−1 xi + j=1 γjk1 xi ∈ Rjm

where φ yk Fk =−yk log pk from (28), with Fk related to pk through (30). This has no
closed form solution. Moreover, the regions corresponding to the different class trees
overlap, so that the solution does not reduce to a separate calculation within each
region ofeach tree in analogy with (18). Following FHT00, we approximate the
solution with a single Newton–Raphson step, using a diagonal approximation to the
Hessian. This decomposes the problem into a separate calculation for each terminal
node of each tree. The result is

xi∈Rjkm y˜ik
γjkm = K − 1
xi∈Rjkm ˜yik 1 − ˜yik (32)
K

This leads to the following algorithm for K-class logistic gradient boosting.

Algorithm 6 (LK TreeBoost).

Fk0 x = 0, k = 1 K
For m = 1 to M do:
pk x = exp Fk x /
K
l=1 exp Fl x , k = 1 K
1202 J. H. FRIEDMAN

For k = 1 to K do:
y˜ik = yik − pk xi , i = 1 N
Rjkm Jj=1 = J-terminal node tree ˜yik xi N1
γjkm = K−1 K
x ∈R ik
i jkm y˜

x ∈R
i jkm ˜yik 1− ˜yik ,j=1J
 Fkm x = Fk m−1 x+
J
j=1 γjkm1 x ∈ Rjkm
endFor
endFor
end Algorithm

can be used to obtain corresponding prob

The final estimates FkM x K1 ability
estimates pkM x K1 through (30). These in turn can be used for classification
ˆ K c k k pk M x
k x = arg min 1≤k≤K k =1

where c k k is the cost associated with predicting the kth class when the truth is k .
Note that for K = 2, Algorithm 6 is equivalent to Algorithm 5. Algorithm 6 bears a
close similarity to the K-class LogitBoost procedure ofFHT00, which is based on
Newton–Raphson rather than gradient descent in function space. In that algorithm K
trees were induced, each using corre sponding pseudoresponses

y˜ik = K − 1 yik − pk xi

pk xi 1 − pk xi (33)
K
and a weight

wk xi = pk xi 1 − pk xi (34) applied to each observation ˜yik xi . The

terminal node updates were
xi∈Rjkmwk xi ˜yik

γjkm =

xi∈Rjkmwk xi
which is equivalent to (32). The difference between the two algorithms is the splitting
criterion used to induce the trees and thereby the terminal regions Rjkm J1 .
The least-squares improvement criterion used to evaluate potential splits ofa
currently terminal region R into two subregions Rl Rr is i2 Rl Rr = wlwr
wl + wr ¯yl −¯yr 2 (35)
where y¯l, y¯r are the left and right daughter response means respectively, and w l, wr
are the corresponding sums ofthe weights. For a given split, using (31) with unit
weights, or (33) with weights (34), give the same val ues for y¯ l, y¯r. However, the
weight sums wl, wr are different. Unit weights (L K TeeBoost) favor splits that are
symmetric in the number of observations in
GREEDY FUNCTION APPROXIMATION 1203

each daughter node, whereas (34) (LogitBoost) favors splits for which the sums ofthe
currently estimated response variances var yik = pk xi 1 − pk xi are more equal.
LK TreeBoost has an implementation advantage in numerical stability. LogitBoost
becomes numerically unstable whenever the value of(34) is close to zero for any
, which happens quite frequently. This is a con
observation xi sequence ofthe difficulty
that Newton–Raphson has with vanishing second derivatives. Its performance is
strongly affected by the way this problem is handled (see FHT00, page 352). L K
TreeBoost has such difficulties only when (34) is close to zero for all observations in
a terminal node. This happens much less frequently and is easier to deal with when it
does happen.
Influence trimming for the multiclass procedure is implemented in the same way
as that for the two-class case outlined in Section 4.5.1. Associated with each
“observation” yik xi is an influence wik = ˜yik 1 − ˜yik which is used for deleting
observations (27) when inducing the kth tree at the current iteration m.

5. Regularization. In prediction problems, fitting the training data too closely can
be counterproductive. Reducing the expected loss on the training data beyond some
point causes the population-expected loss to stop decreas ing and often to start
increasing. Regularization methods attempt to prevent such “overfitting” by
constraining the fitting procedure. For additive expan sions (2) a natural
regularization parameter is the number ofcomponents M. This is analogous to
are consid
stepwise regression where the h x a m M1 ered explanatory variables that
are sequentially entered. Controlling the value of M regulates the degree to which
expected loss on the training data can be minimized. The best value for M can be
estimated by some model selection method, such as using an independent “test” set,
or cross-validation.
Regularizing by controlling the number ofterms in the expansion places an implicit
prior beliefthat “sparse” approximations involving fewer terms are likely to provide
better prediction. However, it has often been found that reg ularization through
shrinkage provides superior results to that obtained by restricting the number
ofcomponents [Copas (1983)]. In the context ofaddi tive models (2) constructed in a
forward stagewise manner (9), (10), a simple shrinkage strategy is to replace line 6
ofthe generic algorithm (Algorithm 1) with

Fm x = Fm−1 (36) x + ν · ρmh x am 0 <ν ≤ 1

and making the corresponding equivalent changes in all ofthe specific algo rithms
(Algorithms 2–6). Each update is simply scaled by the value ofthe “learning rate”
parameter ν.
 Introducing shrinkage into gradient boosting (36) in this manner provides two
regularization parameters, the learning rate ν and the number ofcom ponents M.
Each one can control the degree of fit and thus affect the best value for the other
one. Decreasing the value of ν increases the best value for M. Ideally one should
estimate optimal values for both by minimizing a
1204 J. H. FRIEDMAN

model selection criterion jointly with respect to the values ofthe two param eters.
There are also computational considerations; increasing the size of M produces a
proportionate increase in computation.
We illustrate this ν–M trade-off through a simulation study. The training sample
consists of5000 observations yi xi with

yi = F∗ xi + εi

The target function F∗ x , x ∈ R 10, is randomly generated as described in Section

6.1. The noise ε was generated from a normal distribution with zero mean, and
variance adjusted so that

E ε = 12Ex F∗ x − medianxF∗ x
giving a signal-to-noise ratio of2/1. For this illustration the base learner h x a is taken
to be an 11-terminal node regression tree induced in a best first manner (FHT00). A
general discussion oftree size choice appears in Section 7.
Figure 1 shows the lack offit (LOF) ofLS TreeBoost, LAD TreeBoost, and L2
TreeBoost as a function ofnumber ofterms (iterations) M, for several
values ofthe shrinkage parameter ν ∈ 1 0 0 25 0 125 0 06 . For the first
two methods, LOF is measured by the average absolute error ofthe
estimate F M x relative to that ofthe optimal constant solution

∗
A F M x = Ex F x − F M x
Ex F∗ x − medianxF∗ x (37)
For logistic regression the y-values were obtained by thresholding at the median of
F∗ x over the distribution of x-values; F∗ xi values greater than the median were
assigned yi = 1; those below the median were assigned y i =−1. The Bayes error
rate is thus zero, but the decision boundary is fairly complicated. There
are two LOF measures for L2 TreeBoost; minus twice log-likelihood (“deviance”)
and the misclassification error rate Ex 1 y = sign F M x . The values ofall LOF
measures were computed by using an independent validation data set of10,000
observations.
As seen in Figure 1, smaller values ofthe shrinkage parameter ν (more shrinkage)
are seen to result in better performance, although there is a dimin ishing return for
the smallest values. For the larger values, behavior charac teristic ofoverfitting is
observed; performance reaches an optimum at some value of M and thereafter
diminishes as M increases beyond that point. This effect is much less pronounced
with LAD TreeBoost, and with the error rate criterion ofL 2 TreeBoost. For smaller
values of ν there is less overfitting, as would be expected.
Although difficult to see except for ν = 1, the misclassification error rate (lower
right panel) continues to decrease well after the logistic likelihood has reached its
optimum (lower left panel). Thus, degrading the likelihood by overfitting actually
improves misclassification error rate. Although perhaps
GREEDY FUNCTION APPROXIMATION 1205
Fig. 1. Performance of three gradient boosting algorithms as a function of number of iterations M. The four
curves correspond to shrinkage parameter values of ν ∈ 1 0 0 25 0 125 0 06 and are in that order
(top to bottom) at the extreme right of each plot.

counterintuitive, this is not a contradiction; likelihood and error rate measure different

aspects of fit quality. Error rate depends only on the sign of F M x whereas likelihood
is affected by both its sign and magnitude. Apparently, overfitting degrades the
quality of the magnitude estimate without affecting (and sometimes improving) the
sign. Thus, misclassification error is much less sensitive to overfitting.
Table 1 summarizes the simulation results for several values of ν including those
shown in Figure 1. Shown for each ν-value (row) are the iteration number at which
the minimum LOF was achieved and the corresponding minimizing value (pairs
ofcolumns).
1206 J. H. FRIEDMAN

Table 1
Iteration number giving the best fit and the best fit value for several shrinkage parameter ν-values, with
three boosting methods
 LS A FM x LAD A FM x L2 −2 log like L2 error rate
1.0 15 0.48 19 0.57 20 0.60 436 0.111 0.5 43 0.40 19 0.44 80 0.50 371 0.106 0.25 77 0.34 84 0.38
310 0.46 967 0.099 0.125 146 0.32 307 0.35 570 0.45 580 0.098 0.06 326 0.32 509 0.35 1000 0.44
994 0.094 0.03 855 0.32 937 0.35 1000 0.45 979 0.097

The ν–M trade-off is clearly evident; smaller values of ν give rise to larger optimal
M-values. They also provide higher accuracy, with a diminishing return for ν< 0 125.
The misclassification error rate is very flat for M 200, so that optimal M-values for it
are unstable.
Although illustrated here for just one target function and base learner (11- terminal
node tree), the qualitative nature ofthese results is fairly universal. Other target
functions and tree sizes (not shown) give rise to the same behav ior. This suggests
that the best value for ν depends on the number ofiterations M. The latter should be
made as large as is computationally convenient or feasible. The value of ν should
then be adjusted so that LOF achieves its min imum close to the value chosen for M.
IfLOF is still decreasing at the last iteration, the value of ν or the number ofiterations
M should be increased, preferably the latter. Given the sequential nature of the
algorithm, it can eas ily be restarted where it finished previously, so that no
computation need be repeated. LOF as a function of iteration number is most
conveniently esti mated using a left-out test sample.
As illustrated here, decreasing the learning rate clearly improves perfor mance,
usually dramatically. The reason for this is less clear. Shrinking the model update
(36) at each iteration produces a more complex effect than direct proportional
shrinkage ofthe entire model

F ν x = ¯y + ν · F (38) M x − ¯y

where F M x is the model induced without shrinkage. The update ρ mh x am at each

iteration depends on the specific sequence ofupdates at the previous iterations.
Incremental shrinkage (36) produces very different models than global shrinkage
(38). Empirical evidence (not shown) indicates that global shrinkage (38) provides at
best marginal improvement over no shrinkage, far from the dramatic effect of
incremental shrinkage. The mystery underlying the success ofincremental shrinkage
is currently under investigation.

6. Simulation studies. The performance of any function estimation method

depends on the particular problem to which it is applied. Important characteristics of
problems that affect performance include training sample size N, true underlying
“target” function F∗ x (1), and the distribution of
GREEDY FUNCTION APPROXIMATION 1207

the departures, ε, of y x from F∗ x . For any given problem, N is always known and
sometimes the distribution of ε is also known, for example when y is binary
(Bernoulli). When y is a general real-valued variable the distribution of ε is seldom
known. In nearly all cases, the nature of F∗ x is unknown.
In order to gauge the value ofany estimation method it is necessary to accu rately
evaluate its performance over many different situations. This is most conveniently
accomplished through Monte Carlo simulation where data can be generated
according to a wide variety ofprescriptions and resulting perfor mance accurately
calculated. In this section several such studies are presented in an attempt to
understand the properties ofthe various Gradient TreeBoost procedures developed in
the previous sections. Although such a study is far more thorough than evaluating
the methods on just a few selected examples, real or simulated, the results ofeven a
large study can only be regarded as suggestive.

6.1. Randomfunction generator. One ofthe most important character istics of any
problem affecting performance is the true underlying target function F ∗ x (1). Every
method has particular targets for which it is most appropriate and others for which it
is not. Since the nature of the target func tion can vary greatly over different
problems, and is seldom known, we com pare the merits ofregression tree gradient
boosting algorithms on a variety of different randomly generated targets. Each one
takes the form

F∗ x = 20

algl zl (39)
l=1

The coefficients al 201 are randomly generated from a uniform distribution a l U −1

1 . Each gl zl is a function ofa randomly selected subset, ofsize n l, ofthe n-input
variables x. Specifically,
n
zl = xPl j l

j=1

where each Pl is a separate random permutation ofthe integers 1 2 n . The size

ofeach subset nl is itselftaken to be random, n l = 1 5 + r , with r being drawn from an
exponential distribution with mean λ = 2. Thus, the expected number ofinput
variables for each gl zl is between three and four. However, most often there will be
fewer than that, and somewhat less often, more. This reflects a bias against strong
very high-order interaction effects. However, for any realized F ∗ x there is a good
chance that at least a few ofthe 20 functions g l zl will involve higher-order
interactions. In any case, F∗ x will be a function ofall, or nearly all, ofthe input
variables.
Each gl zl is an nl-dimensional Gaussian function
gl zl = exp −12 zl − µl TVl zl − µl (40)
where each ofthe mean vectors µl 201 is randomly generated from the same
distribution as that ofthe input variables x. The nl × nl covariance matrix Vl
1208 J. H. FRIEDMAN

is also randomly generated. Specifically,

Vl = UlDlUTl

where Ul is a random orthonormal matrix (uniform on Haar measure) and Dl = diag

d1l ···dnll . The square roots ofthe eigenvalues are randomly gen

erated from a uniform distribution djl U a b , where the limits a b depend

on the distribution ofthe input variables x.
For all ofthe studies presented here, the number ofinput variables was taken to be
n = 10, and their joint distribution was taken to be standard nor mal x N 0 I . The
eigenvalue limits were a = 0 1 and b = 2 0. Although the tails ofthe normal distribution
are often shorter than that ofdata encountered in practice, they are still more realistic
than uniformly distributed inputs often used in simulation studies. Also, regression
trees are immune to the effects of long-tailed input variable distributions, so shorter
tails gives a relative advan tage to competitors in the comparisons.
In the simulation studies below, 100 target functions F ∗ x were randomly
generated according to the above prescription (39), (40). Performance is evalu ated
in terms ofthe distribution ofapproximation inaccuracy [relative approx imation error
(37) or misclassification risk] over these different targets. This approach allows a
wide variety of quite different target functions to be gen erated in terms ofthe shapes
oftheir contours in the ten-dimensional input space. Although lower order interactions
are favored, these functions are not especially well suited to additive regression
trees. Decision trees produce ten sor product basis functions, and the components g l
zl ofthe targets F∗ x are not tensor product functions. Using the techniques
described in Section 8, visualizations ofthe dependencies ofthe first randomly
generated function on some ofits more important arguments are shown in Section
8.3.
Although there are only ten input variables, each target is a function of all ofthem.
In many data mining applications there are many more than ten inputs. However, the
relevant dimensionalities are the intrinsic dimensional ity ofthe input space, and the
number ofinputs that actually influence the output response variable y. In problems
with many input variables there are usually high degrees ofcollinearity among many
ofthem, and the number of roughly independent variables (approximate intrinsic
dimensionality) is much smaller. Also, target functions often strongly depend only on
a small subset of all ofthe inputs.

6.2. Error distribution. In this section, LS TreeBoost, LAD TreeBoost, and M

TreeBoost are compared in terms oftheir performance over the 100 target functions
for two different error distributions. Best-first regression trees with 11 terminal nodes
were used with all algorithms. The breakdown param eter for the M TreeBoost was
set to its default value α = 0 9. The learning rate parameter (36) was set to ν = 0 1 for
all TreeBoost procedures in all of the simulation studies.
GREEDY FUNCTION APPROXIMATION 1209

One hundred data sets yi xi N1 were generated according to yi = F∗

xi + εi

where F∗ x represents each ofthe 100 target functions randomly generated as

described in Section 6.1. For the first study, the errors ε i were generated from a
normal distribution with zero mean, and variance adjusted so that

E ε = Ex F∗ x − medianxF∗ (41) x
giving a 1/1 signal-to-noise ratio. For the second study the errors were gen erated
from a “slash” distribution, εi = s · u/v , where u N 0 1 and v U 0 1 . The scale
factor s is adjusted to give a 1/1 signal-to-noise ratio (41). The slash distribution has
very thick tails and is often used as an extreme to test robustness. The training
sample size was taken to be N = 7500, with 5000 used for training, and 2500 left out
as a test sample to estimate the optimal number ofcomponents M. For each ofthe
100 trials an additional validation sample of5000 observations was generated
(without error) to evaluate the approximation inaccuracy (37) for that trial.
The left panels ofFigure 2 show boxplots ofthe distribution ofapproxima tion
inaccuracy (37) over the 100 targets for the two error distributions for each ofthe
three methods. The shaded area ofeach boxplot shows the interquar tile range ofthe
distribution with the enclosed white bar being the median.
Fig. 2. Distribution of absolute approximation error (left panels) and error relative to the best (right panels)
for LS TreeBoost, LAD TreeBoost and M TreeBoost for normal and slash error dis tributions. LS TreeBoost,
performs best with the normal error distribution. LAD TreeBoost and M TreeBoost both perform well with
slash errors. M TreeBoost is very close to the best for both error distributions. Note the use of logarithmic
scale in the lower right panel.
1210 J. H. FRIEDMAN

The outer hinges represent the points closest to (plus/minus) 1.5 interquar tile range
units from the (upper/lower) quartiles. The isolated bars represent individual points
outside this range (outliers).
These plots allow the comparison ofthe overall distributions, but give no
information concerning relative performance for individual target functions. The right
two panels ofFigure 2 attempt to provide such a summary. They show distributions
oferror ratios, rather than the errors themselves. For each target function and
method, the error for the method on that target is divided by the smallest error
obtained on that target, over all ofthe methods (here three) being compared. Thus,
for each of the 100 trials, the best method receives a value of1.0 and the others
receive a larger value. Ifa particu lar method was best (smallest error) for all 100
target functions, its resulting distribution (boxplot) would be a point mass at the value
1.0. Note that the logarithm ofthis ratio is plotted in the lower right panel.
From the left panels of Figure 2 one sees that the 100 targets represent a fairly
wide spectrum of difficulty for all three methods; approximation errors vary by over a
factor of two. For normally distributed errors LS TreeBoost is the superior performer,
as might be expected. It had the smallest error in 73 ofthe trials, with M TreeBoost
best the other 27 times. On average LS TreeBoost was 0.2% worse than the best, M
TreeBoost 0.9% worse, and LAD TreeBoost was 7.4% worse than the best.
With slash-distributed errors, things are reversed. On average the approxi mation
error for LS TreeBoost was 0.95, thereby explaining only 5% target variation. On
individual trials however, it could be much better or much worse. The performance of
both LAD TreeBoost and M TreeBoost was much better and comparable to each
other. LAD TreeBoost was best 32 times and M TreeBoost 68 times. On average
LAD TreeBoost was 4.1% worse than the best, M TreeBoost 1.0% worse, and LS
TreeBoost was 364.6% worse that the best, over the 100 targets.
The results suggest that ofthese three, M TreeBoost is the method of choice. In
both the extreme cases ofvery well-behaved (normal) and very badly behaved (slash)
errors, its performance was very close to that of the best. By comparison, LAD
TreeBoost suffered somewhat with normal errors, and LS TreeBoost was disastrous
with slash errors.

6.3. LS TreeBoost versus MARS. All Gradient TreeBoost algorithms pro duce
piecewise constant approximations. Although the number ofsuch pieces is generally
much larger than that produced by a single tree, this aspect ofthe approximating

function F M x might be expected to represent a disadvantage with respect to

methods that provide continuous approximations, especially when the true underlying
target F∗ x (1) is continuous and fairly smooth. All ofthe randomly generated target
functions (39), (40) are continuous and very smooth. In this section we investigate the
extent ofthe piecewise constant disadvantage by comparing the accuracy ofGradient
TreeBoost with that of MARS [Friedman (1991)] over these 100 targets. Like
TreeBoost, MARS pro duces a tensor product based approximation. However, it uses
continuous func-
GREEDY FUNCTION APPROXIMATION 1211

tions as the product factors, thereby producing a continuous approximation. It also

uses a more involved (stepwise) strategy to induce the tensor products. Since MARS
is based on least-squares fitting, we compare it to LS Tree Boost using normally
distributed errors, again with a 1/1 signal-to-noise ratio (41). The experimental setup
is the same as that in Section 6.2. It is interesting to note that here the performance
of MARS was considerably enhanced by using the 2500 observation test set for
model selection, rather than its default
generalized cross-validation (GCV) criterion [Friedman (1991)]. The top left panel
ofFigure 3 compares the distribution ofMARS average absolute approximation errors,
over the 100 randomly generated target func tions (39), (40), to that ofLS TreeBoost
from Figure 2. The MARS distribution is seen to be much broader, varying by almost
a factor of three. There were many targets for which MARS did considerably better
than LS TreeBoost, and many for which it was substantially worse. This further
illustrates the fact that the nature of the target function strongly influences the relative
performance of different methods. The top right panel of Figure 3 shows the
distribution oferrors, relative to the best for each target. The two methods exhibit
similar performance based on average absolute error. There were a number
oftargets where each one substantially outperformed the other. The bottom two
panels ofFigure 3 show corresponding plots based on root mean squared error. This
gives proportionally more weight to larger errors in assessing lack ofperformance.
For LS TreeBoost the two error measures have close to the same values for all of the
100 targets. However with MARS, root mean squared error is typically 30% higher
than average absolute error. This indicates that MARS predictions tend to be either
very close to, or far from, the target. The errors from LS TreeBoost are more evenly
distributed. It tends to have fewer very large errors or very small errors. The latter
may be a consequence ofthe piecewise constant nature ofthe approximation which
makes it difficult to get arbitrarily close to very smoothly varying targets with
approximations offinite size. As Figure 3 illustrates, relative performance can be quite
sensitive to the criterion used to measure it.
These results indicate that the piecewise constant aspect ofTreeBoost
approximations is not a serious disadvantage. In the rather pristine environ ment
ofnormal errors and normal input variable distributions, it is competi tive with MARS.
The advantage ofthe piecewise constant approach is robust ness; specifically, it
provides immunity to the adverse effects of wide tails and outliers in the distribution
ofthe input variables x. Methods that produce continuous approximations, such as
MARS, can be extremely sensitive to such problems. Also, as shown in Section 6.2,
M TreeBoost (Algorithm 4) is nearly as accurate as LS TreeBoost for normal errors
while, in addition, being highly resistant to output y-outliers. Therefore in data mining
applications where the cleanliness ofthe data is not assured and x- and/or y-outliers
may be present, the relatively high accuracy, consistent performance and robustness
of M TreeBoost may represent a substantial advantage.
1212 J. H. FRIEDMAN

Fig. 3. Distribution of approximation error (left panels) and error relative to the best (right panels) for LS
TreeBoost and MARS. The top panels are based on average absolute error, whereas the bottom ones use
root mean squared error. For absolute error the MARS distribution is wider, indicating more frequent better
and worse performance than LS TreeBoost. MARS performance as measured by root mean squared error
is much worse, indicating that it tends to more frequently make both larger and smaller errors than LS
TreeBoost.

6.4. LK TreeBoost versus K-class LogitBoost and [Link]. In this section the
TreeBoost is compared to that of K-class Log
performance of LK itBoost (FHT00) and
[Link] [Schapire and Singer (1998)] over the 100 randomly generated targets
(Section 6.1). Here K = 5 classes are generated by thresholding each target at its
0.2, 0.4, 0.6 and 0.8 quantiles over the dis tribution ofinput x-values. There are N =
7500 training observations for each trial (1500 per class) divided into 5000 for
training and 2500 for model selec tion (number ofiterations, M). An independently
generated validation sample of5000 observations was used to estimate the error rate
for each target. The
GREEDY FUNCTION APPROXIMATION 1213

Bayes error rate is zero for all targets, but the induced decision boundaries can
become quite complicated, depending on the nature ofeach individual target function
F∗ x . Regression trees with 11 terminal nodes were used for each method.
Figure 4 shows the distribution oferror rate (left panel), and its ratio to the smallest
(right panel), over the 100 target functions, for each of the three methods. The error
rate ofall three methods is seen to vary substantially over these targets. L K TreeBoost
is seen to be the generally superior performer. It had the smallest error for 78 of the
trials and on average its error rate was 0.6% higher than the best for each trial.
LogitBoost was best on 21 of the targets and there was one tie. Its error rate was
3.5% higher than the best on average. [Link] was never the best performer,
and on average it was 15% worse than the best.
Figure 5 shows a corresponding comparison, with the LogitBoost and
[Link] procedures modified to incorporate incremental shrinkage (36), with
the shrinkage parameter set to the same (default) value ν = 0 1 used with L K
TreeBoost. Here one sees a somewhat different picture. Both LogitBoost and
[Link] benefit substantially from shrinkage. The performance of all three
procedures is now nearly the same, with LogitBoost perhaps hav ing a slight
advantage. On average its error rate was 0.5% worse that the best; the
corresponding values for LK TreeBoost and [Link] were 2.3% and 3.9%,
respectively. These results suggest that the relative performance of these methods is
more dependent on their aggressiveness, as parameterized by learning rate, than on
their structural differences. LogitBoost has an addi

Fig. 4. Distribution of error rate on a five-class problem (left panel) and error rate relative to the best (right
panel) for LK TreeBoost, LogitBoost, and [Link]. LK TreeBoost exhibits superior performance.
1214 J. H. FRIEDMAN
Fig. 5. Distribution of error rate on a five-class problem (left panel), and error rate relative to the best (right
panel), for LK TreeBoost, and with proportional shrinkage applied to LogitBoost and RealAdaBoost. Here the
performance of all three methods is similar.

tional internal shrinkage associated with stabilizing its pseudoresponse (33) when the
denominator is close to zero (FHT00, page 352). This may account for its slight
superiority in this comparison. In fact, when increased shrink age is applied to L K
TreeBoost (ν = 0 05) its performance improves, becoming identical to that
ofLogitBoost shown in Figure 5. It is likely that when the shrinkage parameter is
carefully tuned for each of the three methods, there would be little performance
differential between them.

7. Tree boosting. The GradientBoost procedure (Algorithm 1) has two primary

metaparameters, the number ofiterations M and the learning rate parameter ν (36).
These are discussed in Section 5. In addition to these, there are the metaparameters
associated with the procedure used to estimate the base learner h x a . The primary
focus of this paper has been on the use of best-first induced regression trees with a
fixed number ofterminal nodes, J. Thus, J is the primary metaparameter ofthis base
learner. The best choice for its value depends most strongly on the nature ofthe
target function, namely the highest order ofthe dominant interactions among the
variables.
Consider an ANOVA expansion ofa function

Fx= fj xj + fjk xj xk +
fjkl xj xk xl (42) +···
j j k j k l

The first sum is called the “main effects” component of F x . It consists ofa sum
offunctions that each depend on only one input variable. The particular functions f j xj
N
1 are those that provide the closest approximation to F x
GREEDY FUNCTION APPROXIMATION 1215

under this additive constraint. This is sometimes referred to as an “additive” model

because the contributions ofeach xj, fj xj , add to the contributions of the others. This
is a different and more restrictive definition of “additive” than (2). The second sum
consists offunctions ofpairs ofinput variables. They are called the two-variable
“interaction effects.” They are chosen so that along with the main effects they provide
the closest approximation to F x under the limitation ofno more than two-variable
interactions. The third sum represents three-variable interaction effects, and so on.
The highest interaction order possible is limited by the number ofinput variables n.
However, especially for large n, many target functions F ∗ x encountered in practice
can be closely approximated by ANOVA decomposi tions ofmuch lower order. Only
the first few terms in (42) are required to capture the dominant variation in F ∗ x . In
fact, considerable success is often achieved with the additive component alone
[Hastie and Tibshirani (1990)]. Purely additive approximations are also produced by
the “naive” -Bayes method [Warner, Toronto, Veasey and Stephenson (1961)], which
is often highly successful in classification. These considerations motivated the bias
toward lower-order interactions in the randomly generated target functions (Section
6.1) used for the simulation studies.
The goal offunction estimation is to produce an approximation F x that closely
matches the target F∗ x . This usually requires that the dominant interaction order of
∗
F x be similar to that of F x . In boosting regression trees, the interaction
order can be controlled by limiting the size ofthe indi vidual trees induced
at each iteration. A tree with J terminal nodes produces a function with
interaction order at most min J − 1 n . The boosting pro cess is additive, so
the interaction order ofthe entire approximation can be no larger than the
largest among its individual components. Therefore, with any ofthe
TreeBoost procedures, the best tree size J is governed by the effec tive
interaction order ofthe target F ∗ x . This is usually unknown so that J becomes a
metaparameter ofthe procedure to be estimated using a model selection criterion
such as cross-validation or on a left-out subsample not used in training. However, as
discussed above, it is unlikely that large trees would ever be necessary or desirable.
Figure 6 illustrates the effect of tree size on approximation accuracy for
the 100 randomly generated functions (Section 6.1) used in the simulation
studies. The experimental set-up is the same as that used in Section 6.2.
Shown is the distribution ofabsolute errors (37) (left panel), and errors
relative to the lowest for each target (right panel), for J ∈ 2 3 6 11 21 .
The first value J = 2 produces additive main effects components only; J = 3
produces additive and two-variable interaction terms, and so on. A J
terminal node tree can produce interaction levels up to a maximum ofmin
J−1 n , with typical values being less than that, especially when J − 1 n.
As seen in Figure 6 the smallest trees J ∈ 2 3 produce lower accuracy
on average, but their distributions are considerably wider than the others.
This means that they produce more very accurate, and even more very
inaccurate,
1216 J. H. FRIEDMAN
Fig. 6. Distribution of absolute approximation error (left panel) and error relative to the best (right panel) for
LS TreeBoost with different sized trees, as measured by number of terminal nodes J. The distribution using
the smallest trees J ∈ 2 3 is wider, indicating more frequent better and worse performance than with the
larger trees, all of which have similar performance.

approximations. The smaller trees, being restricted to low-order

interactions, are better able to take advantage oftargets that happen to be
oflow interaction level. However, they do quite badly when trying to
approximate the high order interaction targets. The larger trees J ∈ 6 11
21 are more consistent. They sacrifice some accuracy on low-order
interaction targets, but do much better on the higher-order functions.
There is little performance difference among the larger trees, with perhaps
some slight deterioration for J = 21. The J = 2 trees produced the most
accurate approximation eight times; the corresponding numbers for J ∈ 3
6 11 21 were 2, 30, 31, 29, respectively. On average the J = 2 trees had
errors 23.2% larger than the lowest for each target, while the others had
corresponding values of16.4%, 2.4%, 2.2% and 3.7%, respectively. Higher
accuracy should be obtained when the best tree size J is individually
estimated for each target. In practice this can be accomplished by
evaluating the use of different tree sizes with an independent test data
set, as illustrated in Section 9.

8. Interpretation. In many applications it is useful to be able to inter pret the

derived approximation F x . This involves gaining an understanding ofthose

particular input variables that are most influential in contributing to its variation, and

the nature ofthe dependence of F x on those influential inputs. To the extent that F
x at least qualitatively reflects the nature ofthe target function F ∗ x (1), such tools
can provide information concerning the underlying relationship between the inputs x
and the output variable y. In this section, several tools are presented for interpreting
TreeBoost approxima-
GREEDY FUNCTION APPROXIMATION 1217

tions. Although they can be used for interpreting single decision trees, they tend to
be more effective in the context of boosting (especially small) trees. These
interpretative tools are illustrated on real data examples in Section 9.

8.1. Relative importance of input variables. Among the most useful des criptions
ofan approximation F x are the relative influences I j, of the individual inputs xj, on
the variation of F x over the joint input variable distribution. One such measure is

Ij =  2
1/2
Ex ∂F x
 
(43)
∂xj · varx xj

For piecewise constant approximations produced by decision trees, (43) does not
strictly exist and it must be approximated by a surrogate measure that reflects its
properties. Breiman, Friedman, Olshen and Stone (1983) proposed

I
2
j T =J −1
2
ıˆ (44) t 1 vt = j
t=1

where the summation is over the nonterminal nodes t ofthe J-terminal node tree T, v t
is the cor
is the splitting variable associated with node t, and ıˆ 2t responding empirical
improvement in squared error (35) as a result ofthe split. The right-hand side of(44) is
associated with squared influence so that its units correspond to those of(43).
Breiman, Friedman, Olshen and Stone (1983) used (44) directly as a measure
ofinfluence, rather than squared influ ence. For a collection ofdecision trees T m M1 ,
obtained through boosting, (44) can be generalized by its average over all ofthe
trees,

2 1 M
I j= M
2
I j (45) Tm
m=1

in the sequence.
The motivation for (44), (45) is based purely on heuristic arguments. As a partial
justification we show that it produces expected results when applied in the simplest
context. Consider a linear target function

F∗ x = a 0 + n

ajxj (46)
j=1

where the covariance matrix ofthe inputs is a multiple ofthe identity E x x −¯x
x −¯x T= cIn
In this case the influence measure (43) produces

Ij = aj (47)
1218 J. H. FRIEDMAN
 Table 2 shows the results ofa small simulation study similar to those in Section 6,
but with F∗ x taken to be linear (46) with coefficients

aj = −1 j (48) j

and a signal-to-noise ratio of1/1 (41). Shown are the mean and standard deviation
ofthe values of(44), (45) over ten random samples, all with F ∗ x given by (46), (48).
The influence ofthe estimated most influential variable x j∗ is arbitrarily assigned the
value Ij∗ = 100, and the estimated values of the others scaled accordingly. The
estimated importance ranking ofthe input variables was correct on every one ofthe
ten trials. As can be seen in Table 2, the estimated relative influence values are
consistent with those given by (47) and (48).
In Breiman, Friedman, Olshen and Stone 1983, the influence measure (44) is
augmented by a strategy involving surrogate splits intended to uncover the masking
ofinfluential variables by others highly associated with them. This strategy is most
helpful with single decision trees where the opportunity for variables to participate in
splitting is limited by the size J ofthe tree in (44). In the context ofboosting, however,
the number ofsplitting opportunities is vastly increased (45), and surrogate
unmasking is correspondingly less essential.
In K-class logistic regression and classification (Section 4.6) there are K (logistic)
regression functions FkM x Kk=1, each described by a sequence of M trees. In this
case (45) generalizes to

2 1 M
I jk = M
2
I j (49) Tkm
m=1

where Tkm is the tree induced for the kth class at iteration m. The quantity I jk can be
interpreted as the relevance ofpredictor variable xj in separating class k from the
other classes. The overall relevance of xj can be obtained by

Table 2
Estimated mean and standard deviation of input variable
relative influence for a linear target function

VariableMean Standard

10 100.0 0.0
9 90.3 4.3
8 80.0 4.1
7 69.8 3.9
6 62.1 2.3
5 51.7 2.0
4 40.3 4.2
3 31.3 2.9
2 22.2 2.8
1 13.0 3.2
GREEDY FUNCTION APPROXIMATION 1219 averaging over all

classes

1 K
I j = K k=1 I jk
However, the individual I jk themselves can be quite useful. It is often the case that
different subsets of variables are highly relevant to different subsets of classes. This
more detailed knowledge can lead to insights not obtainable by examining only
overall relevance.

8.2. Partial dependence plots. Visualization is one ofthe most powerful

interpretational tools. Graphical renderings ofthe value of F x as a func tion ofits
arguments provides a comprehensive summary ofits dependence on the joint values
ofthe input variables. Unfortunately, such visualization is limited to low-dimensional
arguments. Functions ofa single real-valued vari able x F x can be plotted as a
graph ofthe values of F x against each corresponding value of x. Functions ofa
single categorical variable can be represented by a bar plot, each bar representing
one ofits values, and the bar height the value ofthe function. Functions oftwo real-
valued variables can be pictured using contour or perspective mesh plots. Functions
ofa categorical variable and another variable (real or categorical) are best
summarized by a sequence of(“trellis”) plots, each one showing the dependence of F
x on the second variable, conditioned on the respective values ofthe first variable
[Becker and Cleveland (1996)].
Viewing functions of higher-dimensional arguments is more difficult. It is therefore
useful to be able to view the partial dependence of the approximation F x on
selected small subsets ofthe input variables. Although a collection of such plots can
seldom provide a comprehensive depiction ofthe approximation, it can often produce

helpful clues, especially when F x is dominated by low order interactions (Section

7).
Let zl be a chosen “target” subset, ofsize l, ofthe input variables x, zl = z1
zl ⊂ x1 xn

and z\l be the complement subset

z\l ∪ zl = x

The approximation F x in principle depends on variables in both subsets F x = F

zl z\l

Ifone conditions on specific values for the variables in z\l, then F x can be
considered as a function only of the variables in the chosen subset zl,

F z\l zl = F zl z\l (50)

1220 J. H. FRIEDMAN

In general, the functional form of F z\l zl will depend on the particular values chosen
for z\l. If, however, this dependence is not too strong then the average function

F l zl = Ez\l F x = F zl z\l p\l z\l dz\l (51)

can represent a useful summary of the partial dependence of F x on the chosen

variable subset zl. Here p\l z\l is the marginal probability density
of z\l,
p x dzl (52)
p\l z\l =
where p x is the joint density ofall ofthe inputs x. This complement marginal density
(52) can be estimated from the training data, so that (51) becomes
1 N
F l zl = N
F zl zi \l (53)
i=1

In the special cases where the dependence of F x on zl is additive, F x =

F l zl + F \l z\l (54)
or multiplicative,

F x = F l zl · F \l z\l (55)
the form of F z\l zl (50) does not depend on the joint values ofthe complement
variables z\l. Then F l zl (51) provides a complete description ofthe nature ofthe

variation of F x on the chosen input variable subset zl.

An alternative way ofsummarizing the dependence of F x on a subset zl is to
directly model F x as a function of zl on the training data F l zl = Ex F x zl = F x p
z\l zl dz\l (56)

However, averaging over the conditional density in (56), rather than the marginal
density in (51), causes F l zl to reflect not only the dependence of F x on the
, but in addition, apparent depen
selected variable subset zl dencies induced solely
by the associations between them and the complement variables z \l. For example,
(56)
ifthe contribution of zl happens to be additive (54) or multiplicative (55), F l zl
would not evaluate to the correspond
ing term or factor F l zl , unless the joint density
p x happened to be the product

p x = pl zl · p\l z\l (57)

Partial dependence functions (51) can be used to help interpret models pro duced
by any “black box” prediction method, such as neural networks, support vector
machines, nearest neighbors, radial basis functions, etc. When there are a large
number ofpredictor variables, it is very useful to have a measure of
GREEDY FUNCTION APPROXIMATION 1221

relevance (Section 8.1) to reduce the potentially large number variables and variable
combinations to be considered. Also, a pass over the data (53) is required to
evaluate each F l zl for each set of joint values zl ofits argument. This can be time-
consuming for large data sets, although subsampling could help somewhat.
 For regression trees based on single-variable splits, however, the partial
(51) is straight
dependence of F x on a specified target variable subset zl forward to
evaluate given only the tree, without reference to the data itself (53). For a specific
set ofvalues for the variables zl, a weighted traversal of the tree is performed. At the
root ofthe tree, a weight value of1 is assigned. For each nonterminal node visited, ifits
split variable is in the target subset z l, the appropriate left or right daughter node is
visited and the weight is not modified. Ifthe node’s split variable is a member ofthe
complement subset z\l, then both daughters are visited and the current weight is
multiplied by the fraction of training observations that went left or right, respectively,
at that node.
Each terminal node visited during the traversal is assigned the current value ofthe

weight. When the tree traversal is complete, the value of F l zl is the corresponding
weighted average ofthe F x values over those termi nal nodes visited during the
tree traversal. For a collection of M regression trees, obtained through boosting, the
results for the individual trees are simply averaged.
For purposes ofinterpretation through graphical displays, input variable
subsets oflow cardinality (l ≤ 2) are most useful. The most informative of
such subsets would likely be comprised ofthe input variables deemed to be
among the most influential (44), (45) in contributing to the variation of F x .
Illustrations are provided in Sections 8.3 and 9.
The closer the dependence of F x on the subset zl is to being additive (54) or
multiplicative (55), the more completely the partial dependence function F l zl (51)
captures the nature ofthe influence ofthe variables in z l on the derived approximation
F x . Therefore, subsets z l that group together those influential inputs that have
complex [nonfactorable (55)] interactions between them will provide the most
revealing partial dependence plots. As a diagnostic, both F l zl and F l z\l can be
separately computed for candidate subsets. The value ofthe multiple correlation over
the training data between F x and F l zl F \l z\l and/or F l zl · F \l z\l can be used to
gauge the degree of additivity and/or factorability of F x with respect to a chosen
subset zl. As an additional diagnostic, F z\l zl (50) can be computed for a small
number of z\l-values randomly selected from the training data. The resulting functions
of zl can be compared to F l zl to judge the variability ofthe partial dependence of F x
on zl, with respect to changing values of z\l.
In K-class logistic regression and classification (Section 4.6) there are K (logistic)
regression functions Fk x Kk=1. Each is logarithmically related to p k x = Pr y = k x
through (29). Larger values of Fk x imply higher
1222 J. H. FRIEDMAN

probability ofobserving class k at x. Partial dependence plots ofeach F k x on variable

subsets zl most relevant to that class (49) provide information on how the input
variables influence the respective individual class probabilities.

8.3. Randomly generated function. In this section the interpretational tools

described in the preceding two sections are applied to the first (ofthe 100) randomly
generated functions (Section 6.1) used for the Monte Carlo studies ofSection 6.
 Figure 7 shows the estimated relative importance (44), (45) ofthe 10 input
predictor variables. Some are seen to be more influential than others, but no small
subset appears to dominate. This is consistent with the mechanism used to generate
these functions.
Figure 8 displays single variable (l = 1) partial dependence plots (53) on the six
most influential variables. The hash marks at the base ofeach plot represent the
deciles ofthe corresponding predictor variable distribution. The piecewise constant
nature ofthe approximation is evident. Unlike most approximation methods, there is
no explicit smoothness constraint imposed upon TreeBoost models. Arbitrarily sharp
discontinuities can be accommo dated. The generally smooth trends exhibited in
these plots suggest that a smooth approximation best describes this target. This is
again consistent with the way these functions were generated.

Fig. 7. Relative importance of the input predictor variables for the first randomly generated function used in
the Monte Carlo studies.
GREEDY FUNCTION APPROXIMATION 1223
Fig. 8. Single-variable partial dependence plots for the six most influential predictor variables for the first
randomly generated function used in the simulation studies.

Figure 9 displays two-variable (l = 2) partial dependence plots on some of the

more influential variables. Interaction effects of varying degrees are indicated among
these variable pairs. This is in accordance with the way in which these target
functions were actually generated (39), (40).
Given the general complexity ofthese generated targets as a function of their
arguments, it is unlikely that one would ever be able to uncover their complete
detailed functional form through a series of such partial dependence plots. The goal
is to obtain an understandable description ofsome ofthe impor tant aspects of the
functional relationship. In this example the target function was generated from a
known prescription, so that at least qualitatively we can verify that this is the case
here.

9. Real data. In this section the TreeBoost regression algorithms are illustrated on
two moderate-sized data sets. The results in Section 6.4 suggest that the properties
TreeBoost are very sim
ofthe classification algorithm LK ilar to those ofLogitBoost,
which was extensively applied to data in FHT00. The first (scientific) data set
consists ofchemical concentration measurements on rock samples, and the second
(demographic) is sample survey questionnaire data. Both data sets were partitioned
into a learning sample consisting oftwo thirds ofthe data, with the remaining data
being used as a test sample for choosing the model size (number ofiterations M).
The shrinkage parameter (36) was set to ν = 0 1.
1224 J. H. FRIEDMAN
Fig. 9. Two-variable partial dependence plots on a few of the important predictor variables for the first
randomly generated function used in the simulation studies.

9.1. Garnet data. This data set consists ofa sample of N = 13317 garnets collected
from around the world [Griffin, Fisher, Friedman, Ryan and O’ Reilly (1997)]. A garnet
is a complex Ca–Mg–Fe–Cr silicate that commonly occurs as a minor phase in rocks
making up the earth’s mantle. The variables associated with each garnet are the
concentrations ofvarious chemicals and the tectonic plate setting where the rock was
collected:

TiO2 Cr2O3 FeO MnO MgO CaO Zn Ga Sr Y Zr tec

The first eleven variables representing concentrations are real-valued. The last
variable (tec) takes on three categorical values: “ancient stable shields,” “Proterozoic
shield areas,” and “young orogenic belts.” There are no missing values in these data,
but the distribution ofmany ofthe variables tend to be highly skewed toward larger
values, with many outliers.
The purpose ofthis exercise is to estimate the concentration oftitanium (TiO 2) as a
function ofthe joint concentrations ofthe other chemicals and the tectonic plate index.
GREEDY FUNCTION APPROXIMATION 1225

Table 3
Average absolute error of LS TreeBoost, LAD TreeBoost, and M TreeBoost on the garnet
data for varying numbers of terminal nodes in the individual trees

Terminal nodes LS LAD M

2 0.58 0.57 0.57

 3 0.48 0.47 0.46
4 0.49 0.45 0.45
6 0.48 0.44 0.43
11 0.47 0.44 0.43 21 0.46 0.43 0.43

Table 3 shows the average absolute error in predicting the output y-variable,
relative to the optimal constant prediction,

A y F x =Ey x y−F x
Ey y − median y
(58)

based on the test sample, for LS TreeBoost, LAD TreeBoost, and M TreeBoost for
several values ofthe size (number ofterminal nodes) J ofthe constituent trees. Note
that this prediction error measure (58) includes the additive irre ducible error
associated with the (unknown) underlying target function F ∗ x (1). This irreducible
error adds same amount to all entries in Table 3. Thus, differences in those entries
reflect a proportionally greater improvement in approximation error (37) on the target
function itself.
For all three methods the additive (J = 2) approximation is distinctly infe rior to that
using larger trees, indicating the presence of interaction effects (Section 7) among
the input variables. Six terminal node trees are seen to be adequate and using only
three terminal node trees is seen to provide accuracy within 10% ofthe best. The
errors ofLAD TreeBoost and M TreeBoost are smaller than those ofLS TreeBoost
and similar to each other, with perhaps M TreeBoost having a slight edge. These
results are consistent with those obtained in the simulation studies as shown in
Figures 2 and 6.
Figure 10 shows the relative importance (44), (45) ofthe 11 input variables in
predicting TiO2 concentration based on the M TreeBoost approximation using six
terminal node trees. Results are very similar for the other models in Table 3 with
similar errors. Ga and Zr are seen to be the most influential with MnO being
somewhat less important. The top three panels ofFigure 11 show the partial

dependence (51) ofthe approximation F x on these three most influential variables.

The bottom three panels show the partial dependence of F x on the three pairings
of these variables. A strong interaction effect between Ga and Zr is clearly evident. F
x has very little dependence on either variable when the other takes on its smallest
values. As the value ofone ofthem is increased, the dependence of F x on the other
is correspondingly amplified. A somewhat smaller interaction effect is seen between
MnO and Zr.
1226 J. H. FRIEDMAN
Fig. 10. Relative influence of the eleven input variables on the target variation for the garnet data. Ga and
Zr are much more influential that the others.

Fig. 11. Partial dependence plots for the three most influential input variables in the garnet data. Note the
different vertical scales for each plot. There is a strong interaction effect between Zr and Ga, and a
somewhat weaker one between Zr and MnO.
GREEDY FUNCTION APPROXIMATION 1227

Table 4
Variables for the demographic data

Variable Demographic Number values Type

 1 sex 2 cat
2 martial status 5 cat
3 age 7 real
4 education 6 real
5 occupation 9 cat
6 income 9 real
7 years in Bay Area 5 real
8 dual incomes 2 cat
9 number in household 9 real
10 number in household<18 9 real
11 householder status 3 cat
12 type ofhome 5 cat
13 ethnic classification 8 cat
14 language in home 3 cat

9.2. Demographic data. This data set consists of N = 9409 questionnaires filled
out by shopping mall customers in the San Francisco Bay Area [Impact Resources,
Inc, Columbus, Ohio (1987)]. Here we use answers to the first 14 questions, relating
to demographics, for illustration. These questions are listed in Table 4. The data are
seen to consist ofa mixture ofreal and categorical variables, each with a small
numbers ofdistinct values. There are many miss ing values.
We illustrate TreeBoost on these data by modeling income as a function of the
other 13 variables. Table 5 shows the average absolute error in predicting income,
relative to the best constant predictor (58), for the three regression TreeBoost
algorithms.
There is little difference in performance among the three methods. Owing to the
highly discrete nature ofthese data, there are no outliers or long-tailed distributions
among the real-valued inputs or the output y. There is also very little reduction in
error as the constituent tree size J is increased, indicating

Table 5
Average absolute error of LS TreeBoost, LAD TreeBoost, and M TreeBoost on the
demographic data for varying numbers of terminal nodes in the individual trees

Terminal nodes LS LAD M

2 0.60 0.63 0.61

3 0.60 0.62 0.59
4 0.59 0.59 0.59
6 0.59 0.58 0.59
11 0.59 0.57 0.58 21 0.59 0.58 0.58
1228 J. H. FRIEDMAN
Fig. 12. Relative influence of the 13 input variables on the target variation for the demographic data. No
small group of variables dominate.

lack ofinteractions among the input variables; an approximation additive in the

individual input variables (J = 2) seems to be adequate. Figure 12 shows the relative
importance ofthe input variables in predicting income, based on the (J = 2) LS
TreeBoost approximation. There is no small subset ofthem that dominates. Figure 13
shows partial dependence plots on the six most influential variables. Those for the
categorical variables are rep resented as bar plots, and all plots are centered to have
zero mean over the data. Since the approximation consists of main effects only [first
sum in (42)], these plots completely describe the corresponding contributions fj xj
ofeach ofthese inputs.
There do not appear to be any surprising results in Figure 13. The depen dencies
for the most part confirm prior suspicions and suggest that the approx imation is
intuitively reasonable.

10. Data mining. As “off the shelf ” tools for predictive data mining, the TreeBoost
procedures have some attractive properties. They inherit the favor able
characteristics oftrees while mitigating many ofthe unfavorable ones. Among the
most favorable is robustness. All TreeBoost procedures are invari ant under all
(strictly) monotone transformations of the individual input vari ables. For example,
using xj log xj exj , or xaj as the jth input variable yields the same result. Thus, the need
for considering input variable transformations is eliminated. As a consequence ofthis
invariance, sensitivity to long-tailed
GREEDY FUNCTION APPROXIMATION 1229
Fig. 13. Partial dependence plots for the six most influential input variables in the demographic data. Note
the different vertical scales for each plot. The abscissa values for age and education are codes
representing consecutive equal intervals. The dependence of income on age is nonmonotonic reaching a
maximum at the value 5, representing the interval 45–54 years old.

distributions and outliers is also eliminated. In addition, LAD TreeBoost is completely

robust against outliers in the output variable y as well. M Tree Boost also enjoys a
fair measure of robustness against output outliers.
Another advantage ofdecision tree induction is internal feature selection. Trees
tend to be quite robust against the addition ofirrelevant input variables. Also, tree-
based models handle missing values in a unified and elegant manner [Breiman,
Friedman, Olshen and Stone (1983)]. There is no need to consider external
imputation schemes. TreeBoost clearly inherits these properties as well.
The principal disadvantage ofsingle tree models is inaccuracy. This is a
consequence ofthe coarse nature oftheir piecewise constant approximations,
especially for smaller trees, and instability, especially for larger trees, and the fact
that they involve predominately high-order interactions. All of these are mitigated by
boosting. TreeBoost procedures produce piecewise constant approximations, but the
granularity is much finer. TreeBoost enhances sta bility by using small trees and by
the effect of averaging over many of them. The interaction level of TreeBoost
approximations is effectively controlled by limiting the size ofthe individual constituent
trees.
Among the purported biggest advantages ofsingle tree models is inter pretability,
whereas boosted trees are thought to lack this feature. Small trees
1230 J. H. FRIEDMAN

can be easily interpreted, but due to instability such interpretations should be treated
with caution. The interpretability oflarger trees is questionable [Ripley (1996)].
TreeBoost approximations can be interpreted using partial dependence plots in
conjunction with the input variable relative importance measure, as illustrated in
Sections 8.3 and 9. While not providing a complete description, they at least offer
some insight into the nature of the input– output relationship. Although these tools
can be used with any approxima tion method, the special characteristics oftree-based
models allow their rapid calculation. Partial dependence plots can also be used with
single regression trees, but as noted above, more caution is required owing to
greater instability.
After sorting ofthe input variables, the computation ofthe regression Tree Boost
procedures (LS , LAD , and M TreeBoost) scales linearly with the num ber
ofobservations N, the number ofinput variables n and the number of iterations M. It
scales roughly as the logarithm ofthe size ofthe constituent trees J. In addition, the
classification algorithm LK TreeBoost scales linearly with the number ofclasses K; but
it scales highly sublinearly with the num ber ofiterations M, ifinfluence trimming
(Section 4.5.1) is employed. As a point ofreference, applying M TreeBoost to the
garnet data ofSection 9.1 (N = 13317 n = 11 J = 6 M = 500) required 20 seconds on a
933Mh Pen tium III computer.
As seen in Section 5, many boosting iterations (M 500) can be required to obtain
optimal TreeBoost approximations, based on small values ofthe shrinkage parameter
ν (36). This is somewhat mitigated by the very small size ofthe trees induced at each
iteration. However, as illustrated in Figure 1, improvement tends to be very rapid
initially and then levels off to slower increments. Thus, nearly optimal approximations
can be achieved quite early (M 100) with correspondingly much less computation.
These near-optimal approximations can be used for initial exploration and to provide
an indication ofwhether the final approximation will be ofsufficient accuracy to warrant
continuation. Iflack offit improves very little in the first few iterations (say 100), it is
unlikely that there will be dramatic improvement later on. Ifcontin uation is judged to
be warranted, the procedure can be restarted where it left off previously, so that no
computational investment is lost. Also, one can use larger values ofthe shrinkage
parameter to speed initial improvement for this purpose. As seen in Figure 1, using ν
= 0 25 provided accuracy within 10% of the optimal (ν = 0 1) solution after only 20
iterations. In this case however, boosting would have to be restarted from the
beginning if a smaller shrinkage parameter value were to be subsequently employed.
The ability ofTreeBoost procedures to give a quick indication ofpotential
predictability, coupled with their extreme robustness, makes them a useful
preprocessing tool that can be applied to imperfect data. If sufficient pre dictability is
indicated, further data cleaning can be invested to render it suitable for more
sophisticated, less robust, modeling procedures.
Ifmore data become available after modeling is complete, boosting can be
continued on the new data starting from the previous solution. This usually improves
accuracy provided an independent test data set is used to monitor
GREEDY FUNCTION APPROXIMATION 1231

improvement to prevent overfitting on the new data. Although the accuracy increase
is generally less than would be obtained by redoing the entire analysis on the
combined data, considerable computation is saved.
Boosting on successive subsets ofdata can also be used when there is insuf
ficient random access main memory to store the entire data set. Boosting can be
applied to “arcbites” ofdata [Breiman (1997)] sequentially read into main memory,
each time starting at the current solution, recycling over previous subsets as time
permits. Again, it is crucial to use an independent test set to stop training on each
individual subset at that point where the estimated accuracy ofthe combined
approximation starts to diminish.

Acknowledgments. Helpful discussions with Trevor Hastie, Bogdan Popescu and

Robert Tibshirani are gratefully acknowledged.

REFERENCES
Becker, R. A. and Cleveland, W. S (1996). The design and control ofTrellis display. J. Comput. Statist.
Graphics 5 123–155.
Breiman, L. (1997). Pasting bites together for prediction in large data sets and on-line. Technical report,
Dept. Statistics, Univ. California, Berkeley.
Breiman, L. (1999). Prediction games and arcing algorithms. Neural Comp. 11 1493–1517. Breiman, L.,
Friedman, J. H., Olshen, R. and Stone, C. (1983). Classification and Regression Trees. Wadsworth,
Belmont, CA.
Copas, J. B. (1983). Regression, prediction, and shrinkage (with discussion). J. Roy. Statist. Soc. Ser. B 45
311–354.
Donoho, D. L. (1993). Nonlinear wavelete methods for recovery of signals, densities, and spec tra from
indirect and noisy data. In Different Perspectives on Wavelets. Proceedings of Symposium in
Applied Mathematics (I. Daubechies, ed.) 47 173–205. Amer. Math. Soc., Providence RI.
Drucker, H. (1997). Improving regressors using boosting techniques. Proceedings of Fourteenth
International Conference on Machine Learning (D. Fisher, Jr., ed.) 107–115. Morgan Kaufmann,
San Francisco.
Duffy, N. and Helmbold, D. (1999). A geometric approach to leveraging weak learners. In Computational
Learning Theory. Proceedings of 4th European Conference EuroCOLT99 (P. Fischer and H. U.
Simon, eds.) 18–33. Springer, New York.
Freund, Y. and Schapire, R. (1996). Experiments with a new boosting algorithm. In Machine Learning:
Proceedings of the Thirteenth International Conference 148–156. Morgan Kaufman, San
Francisco.
Friedman, J. H. (1991). Multivariate adaptive regression splines (with discussion). Ann. Statist. 19 1–141.
Friedman J. H., Hastie, T. and Tibshirani, R. (2000). Additive logistic regression: a statistical view ofboosting
(with discussion). Ann. Statist. 28 337–407.
Griffin, W. L., Fisher, N. I., Friedman J. H., Ryan, C. G. and O’Reilly, S. (1999). Cr-Pyrope garnets in
lithospheric mantle. J. Petrology. 40 679–704.
Hastie, T. and Tibshirani, R. (1990). Generalized Additive Models. Chapman and Hall, London. Huber, P.
(1964). Robust estimation ofa location parameter. Ann. Math. Statist. 35 73–101. Mallat, S. and Zhang, Z.
(1993). Matching pursuits with time frequency dictionaries. IEEE Trans. Signal Processing 41 3397–3415.
Powell, M. J. D. (1987). Radial basis functions for multivariate interpolation: a review. In Algo rithms for
Approximation (J. C. Mason and M. G. Cox, eds.) 143–167. Clarendon Press, Oxford.
Ratsch, G., Onoda, T. and Muller, K. R. (1998). Soft margins for AdaBoost. NeuroCOLT Tech nical Report
NC-TR-98-021.
1232 J. H. FRIEDMAN

Ripley, B. D. (1996). Pattern Recognition and Neural Networks. Cambridge Univ. Press. Rumelhart, D. E.,
Hinton, G. E. and Williams, R. J. (1986). Learning representations by back propagating errors. Nature 323
533–536.
Schapire, R. and Singer, Y. (1998). Improved boosting algorithms using confidence-rated predic tions. In
Proceedings of the Eleventh Annual Conference on Computational Learning Theory. ACM, New
York.
Vapnik, V. N. (1995). The Nature of Statistical Learning Theory. Springer, New York. Warner, J. R., Toronto,
A. E., Veasey, L. R. and Stephenson, R. (1961). A mathematical model for medical diagnosis—application
to congenital heart disease. J. Amer. Med. Assoc. 177 177–184.

Department of Statistics
Sequoia Hall
Stanford University
Stanford, California 94305
E-mail: jhf@[Link]