Unsupervised Tree Boosting for Learning Probability Distributions
Naoki Awaya∗1 and Li Ma †2
1

arXiv:2101.11083v7 [stat.ME] 7 Jul 2023

2

Department of Statistics, Stanford University
Department of Statistical Science, Duke University
July 11, 2023

Abstract
We propose an unsupervised tree boosting algorithm for inferring the underlying sampling
distribution of an i.i.d. sample based on fitting additive tree ensembles in a fashion analogous tosupervised tree boosting. Integral to the algorithm is a new notion of “addition” on probability
distributions that leads to a coherent notion of “residualization”, i.e., subtracting a probability
distribution from an observation to remove the distributional structure from the sampling distribution of the latter. We show that these notions arise naturally for univariate distributions throughcumulative distribution function (CDF) transforms and compositions due to several “group-like”
properties of univariate CDFs. While the traditional multivariate CDF does not preserve theseproperties, a new definition of multivariate CDF can restore these properties, thereby allowing thenotions of “addition” and “residualization” to be formulated for multivariate settings as well. This
then gives rise to the unsupervised boosting algorithm based on forward-stagewise fitting of anadditive tree ensemble, which sequentially reduces the Kullback-Leibler divergence from the truth.
The algorithm allows analytic evaluation of the fitted density and outputs a generative model thatcan be readily sampled from. We enhance the algorithm with scale-dependent shrinkage and a
two-stage strategy that separately fits the marginals and the copula. The algorithm then performs
competitively to state-of-the-art deep-learning approaches in multivariate density estimation on
multiple benchmark data sets.

1

Introduction

In supervised learning such as classification and regression, boosting is acknowledged as one of themost powerful algorithms. It is acclaimed for the ability to overcome the curse of dimensionality
and achieve a desirable balance in bias-variance trade-off. The most popular boosting algorithms
can be thought of as sequentially fitting an additive ensemble of weak learners, often in the form∗†

nawaya@stanford.edu
li.ma@duke.edu

1

of regression or classification trees (e.g., Friedman, 2001; Hastie et al., 2009). The success of tree
boosting in supervised problems suggests that a similar strategy might also prevail in unsupervised
problems, where the ultimate objective involves learning the structures of some unknown probability
distribution based on a collection of training data from that distribution.
Our aim in this paper is to formulate a new additive tree model framework for probability distributions along with an unsupervised boosting algorithm that inherits the strength of the supervised
boosting. Our approach is motivated by the observation that some highly effective supervised tree
boosting algorithms are fit in each iteration based on a set of residuals rather than the original observations, thereby substantially simplifying the optimization task in each iteration. To realize thisstrategy in the unsupervised context, we introduce a notion of addition specialized for probability
measures which leads to a natural concept of the residual of an observation after “subtracting” a probability measure from it. It should be clarified that in the unsupervised setting, the natural operation of
addition—such as through taking weighted averages—does not work, as it is not straightforward to findsuch an embedding that renders a conceptually and computationally simple notion of “residualization”
of an observation, which removes a fitted measure from the underlying sampling distribution.
The notions of “addition” and “residualization” for probability distributions are formulated in termsof cumulative distribution function (CDF) transforms and compositions. We start from the case of
univariate distributions, for which the addition of two measures can be defined simply in terms ofa composition of their CDFs whereas the residual of an observation from subtracting a measure issimply the application of the corresponding CDF transform to that observation. In generalizing thisnotion of addition to multivariate distributions on Rd with d > 1, however, the classical notion of themultivariate CDF, which maps Rd to the interval (0, 1] is unsatisfactory, as easily seen, for example,
by the fact that one can neither define a composition of two such CDFs nor define residuals that stilllie in Rd . More fundamentally, multivariate CDFs do not preserve a set of “group-like” propertiesof 1D CDFs that underlie the notion of addition and residualization. Interestingly, it can be shownthat a proper notion of the CDF for multivariate distributions, which maps Rd to Rd , does exist for
probability measures defined on tree-partition structures (tree measures) and it naturally generalizes
the notions of addition and residuals to multivariate settings.
Based on these notions, we introduce an unsupervised tree boosting algorithm for learning probability measures based on forward-stagewise (FS) fitting of an additive tree ensemble. Our algorithm
in each iteration completes two operations that resemble those in supervised boosting: (i) computing
the current residuals by subtracting the fitted measure at the current iteration from the observationsand (ii) fitting a tree-based weak learner on the residuals and adding the estimated distribution to
the current fit. The algorithm enables straightforward analytical evaluation of the probability densityof the fitted distribution and produces a generative model for the fitted measure that can be directly
sampled from.
Because the notion of addition on probability measures in our boosting framework takes the formof function compositions, one can also view our approach as a normalizing flow (NF) (see i.e., Pa-2pamakarios, 2019; Papamakarios et al., 2021; Kobyzev et al., 2020). In NF approaches, one seeks tofind a sequence of transformations that moves the observed distribution into a baseline distribution
such as the uniform and the Gaussian. Most notably, Inouye and Ravikumar (2018) introduced atree-based NF method through composing a class of transforms that are equivalent to the generalized
CDF transforms we introduce. As such, their NF algorithm—called “tree density destructors” is inessence the same as our boosting algorithm aside from the difference in the choice of the base learner,
the strategies in regularization, and the other boosting-inspired specification strategies. These different practical choices we make—which are largely motivated from the boosting perspective—do lead to
substantial differences in empirical performance and computational efficiency, as we will demonstrate
in our numerical experiments. Besides, formulating the algorithm from the boosting perspective also
allows us to provide a more rigorous theoretical grounding for the algorithm. Therefore, our main
contribution is not in the novelty of the algorithm itself, but in connecting boosting and NF in boththeory and practice. In summary, our contributions are:
1. Reformulation of tree-based NFs as tree boosting. We introduce a formal notion of addition onprobability measures that leads to be group structure on (generalized) CDF transforms, based on
which we show that tree-based NFs can be understood as iterative fitting of an additive ensembleof probability measures in ways analogous to iterative fitting of weak learners to residuals underclassical tree boosting for supervised problems.
2. Theoretical justification. The boosting formulation allows us to justify the iterative algorithm formally from a decision-theoretic perspective in terms of sequential minimization of the KullbackLeibler divergence between the model and the unknown true measure. In addition, analogous tosupervised tree boosting (Breiman, 2004), we show that the unsupervised boosting is “highly expressive” in the sense that a wide class of distributions can be represented or well approximatedby a finite combination of highly constrained (or “weak”) tree-based density models.
3. Methodological improvement. The boosting formulation allows us to incorporate methodological
techniques—originally developed for supervised problems. These include guidelines for choosing the number of trees, setting the appropriate level of shrinkage/regularization, and choosing/specifying the base learner. We provide a comprehensive empirical evaluation of our proposed boosting by comparing it with the density destructor (Inouye and Ravikumar, 2018) andthe state-of-the-art NF algorithms such as MAF (Papamakarios et al., 2017) using simulation
examples and benchmark data sets. The results suggest that our new algorithm substantially improves performance over tree density destructors and is competitive to the state-of-the-art deeplearning based NF algorithms at a substantially small computation cost. The decision-theoretic
formulation also leads to a natural measure of variable importance based on the respective contribution of each dimension in reducing the overall KL divergence, which provides additional
insights on the relevance of each dimension in characterizing the underlying distribution andallows effective variable screening.
3We note that boosting for unsupervised learning has been considered by Ridgeway (2002); Rosset
and Segal (2002); Cui et al. (2021). These previous attempts however aim at constructing an ensemble
in the form of a weighted average of probability measures and fit such ensembles through gradientboosting (Mason et al., 1999) under various loss functions. These unsupervised gradient boostingalgorithms have yet to be demonstrated to be computationally efficient or perform well in highdimensional continuous sample spaces.
All proofs are given in Appendix A and B.

2

MethodIn this section, we start by defining notions of addition and residuals for one-dimensional settings andnext generalize them for multivariate distributions to introduce the new tree boosting algorithm.

2.1

CDF-based Addition and Residualization for Univariate Distributions

Without loss of generality, let (0, 1] represent the one-dimensional sample space. For ease of exploration, we shall assume that the distributions are absolutely continuous on the sample space withrespect to the Lebesgue measure and have full support.
We first make an observation that if a random variable X ∼ G, its sampling distribution, then
G(X) ∼ Unif(0, 1], where G denotes the CDF of G. (Throughout we will use bold font letters toindicate the CDFs of the corresponding distributions.) As such, the CDF transform “removes” thedistributional structure of G from the sampling distribution of X. Thus one can think of r = G(X)
as a “residual”. Moreover, Unif(0,1] serves as the notion of “zero”, whose CDF is the identity map, inthe space of probability distributions as it is the remaining distribution after “subtracting” the true
sampling distribution G from X.
Next we define a notion of “adding” two distributions G and H that is consistent with the above
notion of “subtraction” or “residualization”. Specifically, the addition of G1 and G2 , denoted as “G1 ⊕G2 ”, should satisfy the property that if X ∼ G1 ⊕ G2 , then r(1) = G1 (X) ∼ G2 . In other words, if thesampling distribution of X is the “sum” of G1 and G2 , then taking the “residual” of X with respect
to G1 should result in a random variable distributed as G2 . Such a notion of addition indeed exists:
G1 ⊕ G2 is the distribution whose CDF is G2 ◦ G1
where “◦” denotes function composition. Note that this notion of addition is not commutative. That
is, G1 ⊕ G2 ̸= G2 ⊕ G1 . Fortunately, as we will see, the operation of fitting an additive ensemble
as in supervised boosting requires only a non-abelian group-structure, which does not require thecommutativeness of the underlying addition. As such, the loss of commutativeness will pose nodifficulty in our construction of additive tree models and later an unsupervised boosting algorithm
based on the new notions of addition and residuals.

4Addition
Residual
Zero

Boosting for regression

Boosting for probability measures

h1 + · · · + hk
Py − k−1l=1 hl (x)0G 1 ⊕ · · · ⊕ Gk
Gk−1 ◦ · · · ◦ G1 (x)Uniform distribution on (0, 1]dTable 1: Key concepts in boosting
By iteratively applying such an addition, one can define the “sum” of k(≥ 1) probability measures
G1 , . . . , Gk . Specifically, the sum of G1 , . . . , Gk ,G1 ⊕ · · · ⊕ Gk is the distribution whose CDF is Gk ◦ · · · ◦ G1 .
The following property provides the basis for sequential addition and residualization, analogous tothose in supervised boosting.
Proposition 1. If G1 , . . . , Gk have full support on (0, 1], (i.e., when the CDFs G1 , . . . , Gk are strictlyincreasing), then for any i = 1, 2, . . . , k − 1
X ∼ G1 ⊕ · · · ⊕ Gk

if and only ifr(i) = Gi ◦ · · · ◦ G1 (X) ∼ Gi+1 ⊕ · · · ⊕ Gk .

Proposition 1 then implies that such residualization can be applied sequentially. That is, if r(k) is
the residual of x after subtracting G1 ⊕ · · · ⊕ Gk , then r(0) = x, and for k ≥ 1r(k) = Gk (r(k−1) ) = Gk ◦ Gk−1 (r(k−2) ) · · · = Gk ◦ Gk−1 ◦ · · · ◦ G1 (r(0) ).

(1)

The “additivity” induced by the composition of CDFs also induces an additivity on the corresponding
log-likelihood. Specifically, suppose gi = dGi /dµ is the probability density function (pdf) of Gi
for i = 1, 2, . . . , k with respect to Lebesgue measure µ. Then the density of the ensemble measure
Fk := G1 ⊕ · · · ⊕ Gk , fk = dFk /dµ, satisfies
fk (x) =

k
Y

gi (r

(i−1)

) or

log fk (x) =

i=1

k
X

log gi (r(i−1) ).

i=1

Table 1 summarizes the corresponding notions of addition, residuals, and zero in supervised boosting (in particular regression) and those in our unsupervised formulation. With these new notions, weare ready to introduce a boosting algorithm for learning one-dimensional distributions. However, themore interesting application involves multivarate (in fact high-dimensional) distributions. As such,
we first generalize these notions to multivariate cases, and then introduce a multivariate version of
our boosting algorithm that contains the (less interesting) univariate scenario as a special case.

5

2.2

Generalization to Multivariate Distributions

The above notions of addition and residuals do not find direct counterparts for multivariate measures
if one uses the traditional definition of CDFs for multivariate distributions. In particular, because
the traditional CDF is a mapping from (0, 1]d to (0, 1] instead of (0, 1]d , we cannot even take thecomposition of the CDFs or compute the residuals, which should remain in the same space as theoriginal observations. Beyond the minimal requirement that the appropriate notion of “CDF” should
map from (0, 1]d to (0, 1]d , it must also enjoy several group-like properties of univariate CDF’s.
We summarize four such properties that the “CDF” must satisfy to allow the definition of addition
and residualization to carry over into the multivariate setting:(C1) G is a mapping from (0, 1]d to (0, 1]d .
(C2) G is uniquely determined by G.
(C3) If X ∼ G, then G(X) ∼ Unif((0, 1]d ), the “zero”.
(C4) If X ∼ G1 ⊕G2 , the distribution is uniquely determined by its “CDF” G2 ◦G1 , then G1 (X) ∼ G2 .
Remark: (C1) and (C2) are needed for defining addition in terms of compositions. (C3) and (C4) areneeded for the proper notion of residuals.
For the purpose of constructing a tree additive ensemble model and a boosting algorithm, one
type of “CDFs” that satisfy these conditions are particularly useful as they are very easy to computefor probability distributions with piecewise constant densities defined on leafs of a recursive dyadic
partitioning of the sample space. As one can imagine, efficient computation of the “CDFs” for treebased models is critical as they will be computed many times during the fit to an additive ensemble.
2.2.1

Characterizing probability measures on a recursive dyadic partition tree

Next we describe the construction of this generalized notion of multivariate CDFs, which we call the“tree-CDF”, due to its connection to recursive bifurcating partition trees. We start by introducing
some additional notation related to recursive dyadic partitions.
A recursive dyadic partition of depth R is a sequence of nested dyadic partitions
1
A , A2 , . . . , AR on the sample space Ω. The first partition A1 only includes Ω, and for k = 2, . . . , R,
the partition Ak consists of all the sets generated by dividing each A ∈ Ak−1 into two children Al
and Ar , where Al ∪ Ar = A and Al ∩ Ar = ∅. (Throughout, we use subscripts l and r to indicatekleft and right children respectively.) We can denote the recursive partition using a tree T = ∪R
k=1 A .
As such, we refer to the sets in the partitions as “nodes”. We call the collection of nodes in AR the“terminal” nodes or “leafs” of T and denote it by L(T ); the nodes in other levels are the “non-leaf”
nodes or “interior” nodes, which we denote by N (T ) = T \L(T ).
We consider partition trees with axis-aligned partition lines. In this case, a node A ∈ T is of thefollowing rectangular form
A = (a1 , b1 ] × · · · × (ad , bd ].
6(2)

Figure 1: Visualization of the tree-based decomposition of a univariate CDF into three local moves.
The dotted lines indicate G(Al | A)/µ(Al | A) or G(Ar | A)/µ(Ar | A) on each node A.
For a non-leaf node A ∈ N (T ), the children Al and Ar are generated by dividing A in one of the ddimensions, say j ∗ ,Al = (a1 , b1 ] × · · · (aj ∗ , cj ∗ ] × · · · (ad , bd ] and Ar = (a1 , b1 ] × · · · (cj ∗ , bj ∗ ] × · · · (ad , bd ].

(3)

In the following, for each partition tree T , we let PT be the class of probability measures that areconditionally uniform on the leafs of T and have full support on Ω. That is,no
PT = G : G has full support on (0, 1]d and G(· | A) = µ(· | A) for every A ∈ L(T ) ,where µ is the uniform distribution, and G(·|A) and µ(·|A) are the corresponding conditional distributions on A.
To generalize the CDF transform from univariate to multivariate cases, first we note an interesting
multi-scale decomposition of the univariate CDF—a univariate CDF transform G for any distribution
G ∈ PT on an observation x can actually be computed sequentially in a fine-to-coarse fashion along
the branch in the partition tree T in which x falls. Specifically, suppose T has depth R, then forkx ∈ (0, 1], let {Ak }R
k=1 be a sequence of nodes in T such that Ak ∈ A andx ∈ AR ⊂ AR−1 ⊂ · · · ⊂ A1 = (0, 1].
The CDF transform G(x) can be decomposed into the composition of a sequence of “local movefunctions”
G(x) = GA1 ◦ · · · ◦ GAR−1 (x),(4)

where for any A = (a, b] ∈ N (T ) with two children Al = (a, c] and Ar = (c, b], the mapping GA :A → A is (up to a normalizing constant µ(A)) the CDF of a dyadic piecewise constant density equalto G(Al |A)/µ(Al ) on Al and G(Ar |A)/µ(Ar ) on Ar . More precisely, GA : A → A is given by
G(Al |A)GA (x) − a=x−aµ(Al |A)for x ∈ Aland7b − GA (x)
G(Ar |A)=b−xµ(Ar |A)for x ∈ Ar .

Note that the conditional measures and the input and output of GA have the following relationship
G(Al |A) > µ(Al |A)⇔

G(Ar |A) < µ(Ar |A)⇔

GA (x) > x,

G(Al |A) < µ(Al |A)⇔

G(Ar |A) > µ(Ar |A)⇔

GA (x) < x.

We call GA a “local move” function because it moves a point in A in the direction of the child node withless (conditional) probability mass than the (conditional) uniform measure as illustrated in Figure 1.
The amount of movement on A is proportional to the probability mass differential between the twochildren of A in G relative to µ.
If we think of applying the univariate CDF transform as “subtracting” the information containedin a probability measure from an observation, the decomposition in Equation 4 indicates that such
subtraction can be done sequentially through the local moves, each subtracting a piece of information
regarding the measure from the observation. This perspective leads to a generalization of the CDF
transform for the multivariate case as we describe below.
For a point x = (x1 , . . . , xd ) ∈ Ω = [0, 1)d , again let T be a recursive dyadic partition tree of depth
R on the sample space Ω, and {Ak }R
k=1 the sequence of nodes in T that contains x as before. Then we
define a mapping G : (0, 1]d → (0, 1]d in terms of a sequence of fine-to-coarse local moves along thatbranch in T . Specifically, for a node A ∈ N (T ) as in Equation 2 with children Al and Ar attained
from dividing A in the j ∗ th dimension as described in Equation 3, we define a local move mapping
GA : A → A such that for any x ∈ A, GA (x) = (GA,1 (x), . . . , GA,d (x)) where GA,j (x) = x for allj ̸= j ∗ , andGA,j ∗ (x) − aj ∗G(Al |A)=xj ∗ − aj ∗µ(Al |A)for x ∈ Alandbj ∗ − GA,j ∗ (x)G(Ar |A)=bj ∗ − xj ∗µ(Ar |A)for x ∈ Ar .

As illustrated in Figure 2, similar to the univariate case, the local move mapping is nothing but (upto a normalizing constant µ(A)) the CDF of a dyadic piecewise constant density on A, except that
now in the multivariate setting there are a total of d directions in which such a dyadic split can takeplace. As a transform, it moves x in the direction of the child node with less probability mass relative
to the uniform measure.
As before, we now define a mapping G : (0, 1]d → (0, 1]d , called a “tree-CDF”, as the compositionof these local move functions. That is,G(x) = GA1 ◦ · · · ◦ GAR−1 (x).
G is injective from (0, 1]d to (0, 1]d for any G with full support on (0, 1]d . Additionally, because GA
is surjective for every A, G is also surjective. One can also show that G is measurable. Hence wehave the following proposition that establishes Conditions (C1) for tree-CDFs, which is essential todefining the addition of multivariate distributions in terms of tree-CDF compositions.
Proposition 2. The tree-CDF mapping G : (0, 1]d 7→ (0, 1]d is bijective and measurable for any
G ∈ PT .
8

Figure 2: Top: Visualization of the local move functions in each level under R = 3. The nodes withthe darker color have higher conditional probabilities relative to the uniform measure. Bottom: An
example of GA with j ∗ = 2. The input and the output of GA are indicated by × and ◦, respectively.
The next two theorems show that our construction of G satisfies Conditions (C2) and (C3) as
well. That is, G uniquely determines G and applying the G mapping to an observation effectively“subtracts” the distributional structure in G from the sampling distribution of that observation.
Theorem 1. A measure G ∈ PT for some partition tree T can be determined by the tree-CDF mapping
G as follows
G(B) = µ({G(x) : x ∈ B})

for all B ∈ B(Ω).

Remark: Theorem 1 establishes (C2) and implies that G uniquely determines G, regardless of thetree from which G is defined. However, G is tree-specific, that is, to uniquely determine G we needa pair of the measure G and the finite tree T .
Theorem 2. If X ∼ G ∈ PT , then G(X) ∼ Unif((0, 1]d ). Conversely, if U ∼ Unif((0, 1]d ), thenG−1 (U ) ∼ G.
Remark: Theorem 2 establishes (C3) and shows that if one can compute the inverse map G−1 then one
essentially has a generative model, which allows generating samples from G based on “inverse-CDF”
sampling. More details on this will be given in Section 2.3.
2.2.2

Addition and residualization for multivariate settings

Let G1 , . . . , Gk be a collection of probability measures such that Gl ∈ PTl for l = 1, 2, . . . , k, and
let G1 , . . . , Gk be the corresponding tree-CDFs. As a generalization to the univariate case, next we
9define addition of distributions by composing their tree-CDFs. We first show that such a composition
indeed pins down a unique probability measure.
Lemma 1. For Gl ∈ PTl (l = 1, 2, . . . , k), the mapping Fk : B(Ω) 7→ (0, 1] defined as
Fk (B) = µ({Gk ◦ · · · ◦ G1 (x) : x ∈ B})

for B ∈ B(Ω).

(5)

is a probability measure.
Now we can define the sum of k distributions, G1 ⊕· · ·⊕Gk , as the measure Fk given in Equation 5.
This definition of addition contains the univariate case presented earlier as a special case. We notethat the addition implicitly involves the tree structures T1 , . . . , Tk . This dependency on the trees,
however, is suppressed in the “⊕” notation for simplicity without causing confusion.
Next we turn to the notion of residuals and generalize Proposition 1 to multivariate distributions,which establishes Condition (C4) for tree-CDFs.
Proposition 3. Let G1 , . . . , Gk be a collection of probability measures such that Gl ∈ PTl for l =1, 2, . . . , k. ThenX ∼ G1 ⊕ · · · ⊕ Gk

if and only ifr(i) = Gi ◦ · · · ◦ G1 (X) ∼ Gi+1 ⊕ · · · ⊕ Gk

for any i = 1, 2, . . . , k − 1.
Remark: This proposition implies Condition (C4) by setting k = 2 and i = 1.
Moreover, the sequential update of the residuals given in Equation 1 remains valid. The only
difference is that now the residualization in each step depends on an implicit partition tree structure,encapsulated in the corresponding tree-CDF.

2.3

An Unsupervised Boosting Algorithm based on Forward-stagewise (FS) Fitting

Equipped with the new notions of addition and residuals, we are ready to generalize our unsupervised
boosting algorithm to the multivariate setting based on forward-stagewise (FS) fitting. Suppose wehave an i.i.d. sample x1 , . . . , xn from an unknown distribution F , which we model as an additiveensemble of K probability measures
F = G1 ⊕ · · · ⊕ G K(6)

where each Gk is modeled as a member in PTk for some (unknown) Tk . We introduce an FS algorithm
in which we compute the residuals step-by-step and at the kth step, fit Gk to the current residuals.
The fit at the kth step produces an estimate for Gk along with a partition tree Tk , which is used todefine the tree-CDF in the next step for computing the new residuals.

10Initialization
Set r (0) = (x1 , . . . , xn ).
Forward-stagewise fitting
Repeat the following steps for k = 1, . . . , K:
1. Fit a weak learner that produces a pair of outputs (Gk , Tk ) to the residualized observationsr (k−1) , where Tk is an inferred partition tree and Gk is the tree-CDF for a measure Gk ∈
PTk .
(k)

(k)

(k)

2. Update the residuals r (k) = (r1 , . . . , rn ), where ri(k−1)= Gk (ri).

The output of the boosting algorithm in terms of the collection of pairs (Gk , Tk ) for k = 1, 2, . . . , K
contains all of the information from the data regarding the underlying distribution. (In fact the Gk ’s
alone contain all the relevant information, but the Tk ’s are indispensable for effectively representing
and storing the Gk ’s.)
Next we demonstrate two ways to extract such information. In particular, we show (i) how tocompute the density function of the fitted measure F at any point in the sample space analytically,and (ii) how to use the resulting generative model to draw Monte Carlo samples from the fitted
measure F based on “inverse-CDF” sampling.
2.3.1

Evaluating the Density Function of F .

Density estimation is a common objective in learning multivariate distributions. The next proposition
generalizes the additive decomposition of the log-likelihood for the univariate case and provides arecipe for evaluating the density for the fitted measure F analytically based on the output of the FS
algorithm.
Proposition 4. For any x ∈ (0, 1]d , the density f = dF/dµ for F of the additive form in Equation 6
is given as follows
f (x) =KY

gk (r(k−1) )or

log f (x) =

k=1

K
X

log gk (r(k−1) ),k=1

where gk = dGk /dµ is the density of Gk , r(k−1) = Gk−1 ◦ · · · ◦ G1 (x), is the residual for x after
subtracting G1 ⊕ · · · ⊕ Gk−1 , and in particular r(0) = x. In other words, the density f (x) is exactly the
product of the fitted density of each weak learner evaluated at the corresponding sequence of residuals.
2.3.2

A generative model for F .

It turns out that one can use the classical idea of inverse-CDF sampling to construct a generativemodel for F as a result of Theorem 2. Specifically, we can generate samples from F by first generating11

U ∼ Unif((0, 1]d ) and then compute the following transform−1
F−1 (U ) := G−1
1 ◦ · · · ◦ GK (U ),(7)where G−1
k is the corresponding inverse for the tree-CDF Gk for k = 1, 2, . . . , K. To implement the
sampler, we next obtain the analytic form of the inverse of a tree-CDF.
Recall that in Section 2.2.1 we showed that a tree-CDF G for a measure G ∈ PT can be expressedas the composition of a sequence of local move mappings GA : A → A along each subbranch of T .
−1−1
The inverse of the local move mapping GA can be expressed as G−1
A (y) = (GA,1 (y), . . . , GA,d (y)) for
any y = (y1 , . . . , yd ) ∈ A, whereyj(j ̸= j ∗ ),


G−1
(y
)
=
j
A,j
G′ −1 yj −aj(j = j ∗ ),bj −ajA,jand
a + cj −aj z′ −1jG(Al |A) j
GA,j (zj ) =c + bj −cj {z − G(A | A)}
j

G(Ar |A)jl

if yj ≤ aj + G(Al | A)(bj − aj ),if yj > aj + G(Al | A)(bj − aj ).

With the inverse local move function G−1
A available for all A ∈ N (T ), we can obtain the explicit form
for the inverse tree-CDF G−1 as
G−1 (y) = G(R−1) ◦ · · · ◦ G(1) (y),where for k = 1, . . . , R − 1,
G(k) (y) =X

G−1
A (y)1A (y).

A∈Ak

2.4

Decision-theoretic Considerations

In this subsection we show that our boosting algorithm can be interpreted as fitting the additive modelin Equation 6 by sequentially reducing the Kullback-Leibler divergence.
Let F ∗ be the true sampling distribution for the observations and f ∗ = dF ∗ /dµ its density function.
Again, let F = G1 ⊕ · · · ⊕ GK be the additive model for the distribution and f = dF/dµ its density.
We consider the entropy loss, i.e., the Kullback-Leibler (KL) divergence between F ∗ and F defined as
Z
f∗∗
KL(F ||F ) = log dF ∗ .
(8)f
The next lemma states that the entropy loss in Equation 8 can be decomposed into K components(ignoring a constant) each of which only depends on Gk .

12

Lemma 2. The Kullback-Leibler divergence can be written asZ∗

KL(F ||F ) =log f ∗ dF ∗ −KX{KL(F̃k ||µ) − KL(F̃k ||Gk )}.

(9)

k=1

where F̃k is the true distribution of the residualized observation after subtracting G1 , G2 , . . . , Gk−1 .
That is, F̃k is the true distribution of r(k−1) = Gk−1 ◦ · · · ◦ G1 (X), where X ∼ F ∗ .
Remark: Note that because X ∼ F ∗ , we have F̃1 = F ∗ , and for k = 2, . . . , K, F̃k is in the following
form
−1
d
F̃k (B) = F ∗ (G−1
1 ◦ · · · ◦ Gk−1 (B)) for all B ∈ B((0, 1] ).
The first term on the right-hand side of Equation (9) is a constant. The summand in the second
term is positive as long as the measure Gk is closer to F̃k than the uniform measure µ in terms of
KL divergence. Hence, unless F̃k = µ, the entropy loss could be reduced by adding an additional
measure Gk that is closer to F̃k than µ. In this way, fitting a measure Gk to the residuals r (k) in thekth step of our boosting algorithm can be understood as an operation to sequentially reduce the KL
divergence. Next we turn from the above insight at the population level to the practical strategy at
the finite-sample level for fitting F based on n i.i.d. observations {xi }ni=1 from F ∗ . First note thatR
minimizing the divergence KL(F ∗ ||F ) is equivalent to maximizing the average log-density log f dF ∗ .
Thus with a finite sample, we aim to maximize the sample (average) log-density of the training data,that is,n
1X
log f (xi ).
ni=1

It follows from Proposition 2 that the log-density can also be decomposed into the sum of K components, which we call “improvements”.
Lemma 3. The sample average log-density can be written as
n

K

i=1

k=1

X (n)
1X
log f (xi ) =Dk (Gk ),n(n)where for k = 1, 2, . . . , K, the improvement Dk (Gk ) is1(n)
Dk (Gk ) =
n

n
X

(k−1)

log gk (ri

)

with

(k−1)

ri

= Gk−1 ◦ · · · ◦ G1 (xi ).

i=1
(n)

Accordingly, the next proposition characterizes the “optimal” pair (Gk , Tk ) that maximizes Dk (Gk ).
(n)

Proposition 5. A pair of (Gk , Tk ) maximizes Dk (Gk ) if and only if(n)

Tk ∈ arg max
T ∈T

(n)

X

F̃k (A) log

A∈L(T )

13F̃k (A)
,
µ(A)

(10)

and
(n)

for all A ∈ L(Tk ),Gk (A) = F̃k (A)(n)

where F̃k(11)(k−1) n}i=1 . That is,is the empirical measure of the residuals r (k−1) = {rin(n)F̃k (B) =1X(k−1)
δB (ri)nfor B ∈ B((0, 1]d ).

i=1Remark 1: The summation in Equation 10 is the KL divergence between two discrete probability(n)measures with masses given by {F̃k (A)}A∈L(T ) and {µ(A)}A∈L(T ) respectively. Equation 10 implies
that the “optimal” tree Tk should allow maximal differentiation in KL divergence between the induced
discretizations of F̃ (n) and µ on its leafs. This proposition offers practical guidance on how to choose
a good weak learner, which will be detailed in Section 2.7.1.
Remark 2: As suggested in Lemma 3, the loss is reduced in each step as long as the improvement(n)
Dk (Gk ) is positive, and the improvement is maximized by the measure described in Proposition 5.
However, adopting the “optimal” base learner in fitting Gk as prescribed in Proposition 5 will generally lead to over-fitting. As in supervised boosting (Hastie et al., 2009), additional regularization isnecessary to reduce the variance of the weak learner, and this can be achieved in analogy to supervised
boosting through shrinkage toward the “zero”, here the uniform distribution. As will be detailed in
Section 2.7.2, one can still ensure the improvement to be positive when shrinkage is incorporated in
an appropriate way.

2.5

Group structure of tree-CDFs

In the previous sections we have defined the new operation to add two distributions, and discussed
the “group-like” structure on probability measures it induces. Here we make this notion more formal
by showing that the collection of tree-CDFs indeed form a group.
Proposition 6. Let G be a set of tree-CDFs defined as follows
G = {G | G is a tree-CDF of G ∈ PT for a finite tree T } .
Then, G generates a group under the composition ◦. Specifically, the identity map—which corresponds
to the uniform distribution—is the identity in the group. G is closed under ◦ and each element has aninverse in G.
This group structure is an example of the group-theoretic structure Inouye and Ravikumar (2018)
introduced to the family of density destructor transformations, though they did not show in theparticular case of tree density destructors the group structure exists. Also note that this group is notabelian, because ◦ is not commutative. This is distinct from the usual group structure defined on theclass of tree regressions in supervised settings, as the usual notion of addition is commutative. The
forward-stepwise fitting of an additive ensemble of elements in G does not require the group to beabelian.
14

2.6

Connection to Gradient Boosting

Many of the existing boosting algorithms can be regarded as iterative optimization of loss functions
with gradient descent (Mason et al., 1999; Friedman, 2001). In this subsection we discuss our new
boosting from this perspective and clarify the difference from existing gradient boosting methods.
First we note that as in gradient boosting, our new algorithm can be seen as fitting a linear additive
expansion. As shown in Proposition 4, for the ensemble measure F = G1 ⊕ · · · ⊕ GK , the log-densityof the ensemble measure f = dF/dµ evaluated at the observation xi (i = 1, . . . , n) can be decomposed
as follows:
K
X
(k−1)
log f (xi ) =log gk (ri
),
k=1
(k−1)
where gk = dgk /dµ is the density of Gk and ri= Gk−1 · · · ◦ G1 (xi ) is the residual.

From thisexpression we can see that log f bears resemblance to a “linear additive model” with which we estimate
the log-density function log f ∗ of the unknown measure F ∗ . The fitting of this estimate to the data isevaluated with the sum of the log densities
L(log f ) =nX− log f (xi ).

i=1

Note that Rosset and Segal (2002) also proposes an additive model for density estimation with thesame objective function, but their model is a weighted sum of density functions so is different from(k−1)
our model. We also note that the input of log gk is the residual ri
instead of xi itself due to thedefinition of our new addition rule that involves transformation with the tree-CDFs.
Suppose that we want to update log f , a current estimate of log f ∗ , by adding a new density
function log g, where g is the density function of the new tree measure G, to improve the fitting.
The standard approach in the gradient boosting is evaluating a gradient at the current estimate andapproximating its negative with the new function (Mason et al., 1999; Friedman, 2001). In our case,the gradient is constant: for all i,


∂L
= −1.
∂ log f (x) log f (x)=log f (xi )Approximating its negative, 1, with the new function log g is not reasonable. However, this result
implies that we can maximize the improvement in the loss function by maximizing the sum of the(K−1) nlog-densities evaluated at the current residuals {ri}i=1 ,nX(K−1)

log g(ri),

i=1or equivalently their average, which is exactly what our proposed boosting algorithm iteratively does
to fit the ensemble measure (see Section 2.4). Therefore, while our new boosting method constructs
the ensemble measure in a different way from standard gradient boosting, one can still justify ouralgorithm as a sequential optimization of the average log densities.
15

2.7

Practical Considerations

In this subsection we describe several practical considerations in implementing and applying the boosting algorithm. While they might first appear as technical details, we have found that they are criticalin achieving competitive performance and thus worth elaborating on. Several of these considerations
are drawn from similar considerations in supervised boosting.
2.7.1

Choice of a Weak Learner

Searching over all possible trees to solve Equation 10 in each step of the FS algorithm is computationally prohibitive. Nevertheless, Proposition 5 provides hints on how to choose good weak learners
that improve the KL divergence efficiently over the iterations. The simplest possible choice of a weaklearner, as is often implemented in supervised boosting is to implement a top-down greedy tree learning algorithm that maximizes Equation 10 one split at a time, as is done in fitting classification andregression trees (CART) (Hastie et al., 2009).
In our numerical examples and software, we adopt a weak learner based on a simplified version of an
unsupervised (Bayesian) CART model for probability distributions proposed in Awaya and Ma (2022).
Fitting this weak learner uses a stochastic one-step look-ahead strategy to choose splitting decisions
on each tree node, which generally produces closer approximation to the “optimal” tree splits thangreedy tree algorithms. See Theorem 4.1 in Awaya and Ma (2022) for an asymptotic justification—as
the sample size grows, it produces trees that satisfy Proposition 5 with probability increasing to 1.
Additional details about the weak learner can be found in Appendix C.
It is worth emphasizing that because we are only building “weak” learners that extract a small
fraction of the distributional structure in each iteration, one does not need to be precisely “optimal”
in each iteration. More importantly than being “optimal”, the weak learner should facilitate theappropriate shrinkage to avoid overfitting, which we elaborate in the next subsection.
2.7.2

Regularization through Scale-specific Shrinkage

Just as in supervised boosting, simply adopting the solution for Equation 10 (either exact or approximate) as the fit for Gk in each iteration will typically lead to overfitting even when the complexityof the tree Tk is restricted to be small. In particular, the fitted density will tend to have spikes ator near the training points. To avoid such overfitting, it is necessary to regularize or penalize thenon-smoothness in the fit for each Gk . This can be achieved through shrinkage toward “zero”, or the(n)uniform measure µ, thereby discounting the influence of the residuals (or its empirical measure F̃k )on fitting Gk . In supervised boosting it is typical to introduce a learning rate c0 ∈ (0, 1] that controls
how much shrinkage toward zero is applied in each iteration. In the current context, this traditional
strategy would correspond to setting(n)Gk = (1 − c0 )µ + c0 F̃k .

16

We found that in practice one can further improve upon this shrinkage strategy by allowing different
levels of shrinkage at different scales. The intuition is that depending on the smoothness of theunderlying function, overfitting can be more (or less) likely to happen in learning local details of the
distribution and thus one may benefit from enforcing a level of shrinkage that increases (or not) withthe depth in the tree Tk . Following this intuition, we specify a scale-dependent learning rate as follows
c(A) = c0 · (1 − log2 vol(A))−γ ,where A is a node in Tk , vol(A) is a volume of A. Then the shrinkage toward the uniform can be
specified on each node A ∈ N (Tk ) in terms of the conditional probability on the children of A(n)

Gk (Al | A) = (1 − c(A))µ(Al | A) + c(A)F̃k (Al | A) for A ∈ N (Tk ),(12)

Gk (· | A) = µ(· | A) for A ∈ L(Tk ),
(n)

(n)

(n)

(n)

where Al and Ar are the children nodes of A in Tk , F̃k (Al | A) = F̃k (Al )/F̃k (A) if F̃k (A) > 0(n)and F̃k (Al | A) = µ(Al | A) otherwise.
The node-specific learning rate c(A) controls how strongly one “pulls” the empirical measure F̃ (n)toward the uniform measure µ at the corresponding scale of A. It is specified with two tuning
parameters c0 ∈ (0, 1] and γ ≥ 0. The parameter c0 controls the global level of shrinkage, andwhen γ > 0 we introduce stronger shrinkage for small nodes, imposing stronger penalty on local
spikes. When γ = 0, this shrinkage reduces to the standard single learning rate specification describedabove. In practice, we recommend setting these tuning parameters by cross-validation.
Our next proposition shows that with shrinkage, the sample average log-density is steadily improved in each step of the FS algorithm until the residual distribution becomes the uniform measure.
Proposition 7. For any finite tree Tk , under the definition of Gk given in Equation 12, the improve(n)(n)
ment satisfies Dk (Gk ) ≥ 0 if c(A) ∈ (0, 1] for all A ∈ L(Tk ) unless F̃k (A) is indistinguishable from(n)
the uniform distribution on the tree Tk , that is, F̃k (A) = µ(A) for all A ∈ L(Tk ).
2.7.3

Evaluating Variable Importance

As in supervised learning, it is often desirable to evaluate the contribution of each dimension to theapproximation of the unknown measure F ∗ . Thus we provide a way to quantify variable importancein a conceptually similar manner to what is often used in supervised boosting (see Hastie et al.,2009). We note that Ram and Gray (2011) also introduced a notion of the variable importance indensity trees. While their definition is based on improvement in the L2 loss, ours is based on the KL
divergence, which is consistent with our earlier decision-theoretic discussion.
Specifically, because our boosting algorithm reduces the KL divergence from the unknown measure∗
F , a natural way of quantifying the importance of a variable is adding up the decrease in the KL
divergence due to splitting a tree node in the corresponding dimension. Lemma 3 shows that this

17

(n)

quantity can be expressed as the sum of the improvements Dk (Gk ). In particular, the improvement(n)
Dk (Gk ) can be further decomposed over the splits of the tree Tk as followsX
Gk (Al | A)
Gk (Ar | A)
(n)
(n)
(n)
(n)
F̃k (A) F̃k (Al | A) log
Dk (Gk ) =,+ F̃k (Ar | A) logµ(Al | A)µ(Ar | A)
A∈N (T )(n)where the empirical measure F̃kthe brackets can be written as

is as defined in Proposition 5. Note that the summation inside of

(n)

(n)

KL(F̃k (Al | A)||µ(Al | A)) − KL(F̃k (Al | A)||Gk (Al | A)),
where KL(p||q) = p log(p/q)+(1−p) log[(1−p)/(1−q)], and it quantifies the extent to which splitting A
makes Gk closer to the distribution of the residuals. Based on the decomposition, a natural definitionof the total contribution of dividing in the jth dimension isX
Gk (Al | A)
Gk (Ar | A)
(n)
(n)
(n)
IGk ,j =F̃k (A) F̃k (Al | A) log+ F̃k (Ar | A) log,µ(Al | A)µ(Ar | A)
A∈Nj (T )where Nj (T ) represents the collection of all nodes in T that are split in the jth dimension. Finally we
can define the importance of the jth variable in the additive measure F = G1 ⊕ · · · ⊕ GK by summingover the variable importance across the Gk ’s:
Ij =

K
X

IGk ,j .

k=1

2.7.4

Fitting the Margins and the Copula Separately and Addressing Technical Ties

In the density estimation literature, Lu et al. (2013) suggested a two-stage strategy for estimating
multivariate densities using tree-based models, which separately fits the marginal distributions andthen the dependence (or copula). From our experience, this strategy can often substantially improve
the fit of our unsupervised boosting algorithm.
This two-stage strategy is easy to realize in our algorithm. In the first stage, for each of the dimensions, one can adopt weak learners that are constrained to involving tree-CDFs based on partitions
along that single dimension. Computing the residuals with tree-CDFs defined on such a tree only
removes the marginal distributions from the observations. “Subtracting” all of the marginal distributions from the original observations results in a sample of residuals the remaining distribution withuniform marginals (i.e., the corresponding copula). Then in the second stage, the single-dimension
constraint on the partition trees is removed, and tree-CDFs are then fitted to the copula. The final
fit is simply the sum, in terms of tree-CDF compositions, of all of the marginals and the copula.
A related practical consideration regards tied values in the training data. (Ties in the margins
occur much more frequently than ties that occur simultaneously in all margins and thus the issue isparticularly relevant during the fitting of the marginal distributions in this two-stage strategy.) When
18

tied values occur, either from the actual data generative mechanism or due to technical reasons suchas rounding, the additive tree model itself will only assume that it is due to the actual data generative
mechanism and therefore there must be positive probability mass at those tied values, leading to spikesof estimated densities at those values. In practice, if the data generative mechanism is assumed tobe continuous and the ties are due to technical reasons such as rounding, one can avoid this issue
by a simple preprocessing step for the training data that “smooths out” those spikes by adding small
perturbations before fitting the model. We have found a simple strategy to be effective—when there
are ties in the training data at the same value x and the adjacent values are x− and x+ (x− < x < x+ ),we add uniform perturbation to the training data at x on the support (−(x − x− )/2, (x+ − x)/2).
2.7.5

Choosing the number of trees

For the learning rate and the number of trees, we use the standard strategy for supervised problems(Friedman, 2001) by setting the learning rate (c0 ) small—e.g., 0.1 or 0.01—along with a large numberof trees (e.g., hundreds to thousands).
While our unsupervised boosting algorithm is generally robust with respect to overspecificationin the number of trees, it is nevertheless beneficial to adopt an adaptive stopping strategy, which
terminates the boosting algorithm when substantitive improvement in the fitting is no longer expected,
to avoid excessive computation. We adopt a simple strategy to achieve this. In each iteration of thealgorithm, one can use a portion of the data, for example, 90% of the residuals, to fit the next
weak learner, and use the rest of the data to evaluate this tree measure by computing the average log
densities. We use this quantity to measure the improvement in the fitting. If the average improvement
given by, for example, the most recent 50 trees is non-positive, we terminate the algorithm. We havefound this adaptive stopping strategy highly effective in all of our numerical experiments—it did notaffect predictive performance noticeably while substantially reducing computing time.

2.8

Expressive Power of Additive Tree Ensembles

One interesting question is what kind of probability measures can be well approximated by the ensemble when relatively simple (e.g., shallow) tree-based weak learners are combined. This is often
referred to as the “expressive power” of the model in the machine learning literature, or the “support”
of the model in the statistical literature. The expressive power for several normalizing flows have been
analyzed (Huang et al., 2018; Jaini et al., 2019; Kong and Chaudhuri, 2020), and the large support
property of linear combinations of classification trees has been established in Breiman (2004). We
show next that a similar property holds for our unsupervised tree ensemble under the following set of
conditions.
Assumption 1. For k = 1, 2, ..., K, the pair of the tree Tk and the measure Gk that forms a componentof the tree ensemble satisfies the following conditions:19

1. The tree Tk can be any finite tree formed by the dyadic splitting rule described in Equation 2 andEquation 3. That is, it incorporates an axis-aligned splitting rule with flexible split points.
2. Each Tk has at least d + 1 leaf nodes, where d is the dimension of the sample space.
3. The measure Gk can be any conditionally uniform measure on Tk , namely, for every non-terminal
node A ∈ N (Tk ), the conditional probability Gk (Al | A) can be any value in (0, 1).
Most notable in the assumption is the second condition, which is in sharp contrast to theories on
single tree-based density models in the statistical literature. For instance, a popular tree-based density
model, the Pólya tree (PT) model is shown to have the large support under the assumption that asingle tree has infinite depth (Lavine, 1992; Ghosal and Van der Vaart, 2017). This is not surprising in
the context of additive trees, however, since Proposition 2 of Breiman (2004) for supervised boosting
essentially requires the same condition. With additive ensembles, we can combine small trees toexpress general continuous distributions, as formally stated in the next theorem.
Theorem 3. Let F ∗ be a probability measure that has a bounded density function. Then, under
Assumption 1, for any ϵ > 0, there exists a tree ensemble with a finite number of tree measures,
G1 ⊕ · · · ⊕ GK , that approximates F ∗ in terms of the KL divergence with this precision, i.e.,
KL(F ∗ ||G1 ⊕ · · · ⊕ GK ) < ϵ.

2.9Connection to Normalizing Flows

Under our definition of additive tree ensembles, the unknown distribution of the observation is modeled
as the transformation of the uniform distribution that takes the form
T1 ◦ · · · ◦ TK (U ), U ∼ Unif((0, 1]d ),where Tk = G−1
k . Given this expression, we can find a connection between our new boosting and thenormalizing flow (NF) methods, a class of machine learning algorithms for density estimation. For
comprehensive reviews of the NF, see Papamakarios (2019), Papamakarios et al. (2021), and Kobyzev
et al. (2020). In NF methods, one approximates the observation’s distribution with the transformationof known distributions such as the uniform and the Gaussian, and the transformation is represented
as a composition of multiple functions. From this viewpoint, our boosting method can be consideredan NF method in which we use the inverse tree-CDFs G−1
k as base transformations.
In NF methods, the estimated log-density is written in the form of the log-determinant of a flowtransformation (Papamakarios et al., 2021). We can show that in our case, the log-determinant is
identical to the sum of the log-densities that appears in our ensemble formulation. To see this, let
FK be a composition of tree CDFs Fk = GK ◦ · · · ◦ G1 , which uniquely corresponds to the ensemble
measure written as F = G1 ⊕ · · · ⊕ GK by Theorem 1. We can rewrite the log-determinant using the20

density function of the tree measures gk = dGk /dµ, as follows
log detF (x) =KXlog detGk (r(k−1) ) =

k=1

K
X

log gk (r(k−1) ),k=1

where r(k−1) = Gk−1 ◦ · · · ◦ G1 (x). The first equation follows the chain rule, and the second equation
is obtained by the fact that |detGk (x)| = gk (x) (almost everywhere). The last expression is the sumof the log-densities in our boosting algorithm (see Proposition 4).
From an algorithmic perspective, some NF methods are similar to our boosting algorithm in thatthey adopt an iterative fitting approach. That is, they sequentially transform the observations tomake the distribution close to the known distributions, as we sequentially residualize the observations
with Gk ’s. In particular, Inouye and Ravikumar (2018) introduces an NF method called “density
destructors”, which “substracts” information from i.i.d. samples until they are uniform samples. Notably, one particular type of density destructors—the “tree density destructors”—is defined based on
subtracting a tree-based transforms, which is equivalent to our notion of tree-CDF transform. The
sequential training on the residuals is feasible due to the underlying group structure over the treeCDF transforms. It substantially reduces the computational cost compared to methods that requires
stochastic gradient descent to train.
For most NF models, there tends to be a trade-off between the ease in evaluating the fitted
density or generating samples from the fitted distribution and that in achieving large expressive power(Papamakarios et al., 2021). It is worth noting that density evaluation and simulation given the fitted
additive model are both straightforward to implement under the boosting algorithm because thesetasks only require transforming inputs with the tree-CDF and its inverse function respectively, bothof which are available in closed forms. The computational cost of these tasks is O(RK) and in practice,
the cost is much smaller because the node splitting is often terminated in shallow levels on much of
the sample space. At the same time, the proposed additive tree ensemble is capable of expressing or
approximating general continuous distributions as shown in Section 2.8.

3

Numerical Experiments

In this section, we carry out numerical experiments to evaluate and demonstrate the behavior of our
method. We start with a set of simulations under several representative forms of density functions in
a 48-dimensional sample space, through which we examine the impact of the tuning parameters—thelearning rate, shrinkage parameters, and the number of trees—as well as the two-stage strategy on
the performance of the algorithm. In addition, we compare the predictive performance of our methodwith a state-of-the-art single tree learner as well as the closely related tree density destructor (Inouye
and Ravikumar, 2018). We then provide a comparative study based on several popular benchmark
data sets that pitches our method against several state-of-the-art NF methods including both deeplearning based NFs and the tree density destructor. Finally we demonstrate the computation of21

variable importance in both artificial datasets and the well-known MNIST handwritten digits data
(LeCun et al., 1998).
Unless otherwise noted, we adopt the two-stage strategy and the adaptive stopping discussed in
Section 2.7.4 and in Section 2.7.5, respectively, and set the number of measures for the first stage(estimation of the marginal distributions) to 100 per dimension and for the second stage (estimationof the dependence structures) to 5,000 unless otherwise statied. The evaluation of the computational
cost is done in a single-core AMD EPYC 7002 (2.50GHz) environment.

3.1

Simulation Study in 48D Sample Spaces

We consider three 48-dimensional scenarios whose true densities are illustrated in Figure 3, and set
the sample size n to 10,000. (The detailed information on the data generating process is provided inAppendix D.)
We evaluate the performance under two types of experimental settings. In the first experiment, we
set the value of c0 , the global shrinkage parameter, to 0.01, 0.1 and 0.99 while fixing γ, the shrinkage
parameter that penalizes small nodes, to 0. In this experiment we do not adopt the two-stage strategy
and nor the adaptive stopping in order to examine the effect of changing c0 and the number of trees.
The predictive performance given by the considered methods is visualized in the first row of Figure 4.
For completeness, we also report the predictive score for a state-of-the-art single tree learner (whichis essentially an unregularized version of the base learner in our boosting algorithm) to illustrate thedramatic improvement in the performance through fitting an addive ensemble. We can see that the
boosting substantially outperforms the single-tree method, and that the lower learning rate (c0 ) tends
to result in better predictive scores. Also, when c0 is small (0.01), the number of trees to needed tomaximize performance can be as large as 4000 to 5000.
In the second experiment, we compare the performance under different γ (the shrinkage parameterfor the penalty on small nodes) and evaluate the effect of introducing the two-stage strategy. In thisexperiment γ is set to 0.0, 0.5, . . . , 2.0, and c0 is fixed to 0.1. The results are reported in the second rowof Figure 4. Overall one can see that (i) a penalty for small nodes via a positive γ leads to improved
fit as long as γ is not too large (i.e., ≥ 1); and (ii) the two-stage strategy can further improve the fit.
To demonstrate the practical improvement through boosting inspired specification of our algorithm,
we also provide a comparison with the deep density destructor (DDD) (Inouye and Ravikumar, 2018)
in terms of predictive scores and computation cost. For the DDD algorithm, we use the random
tree destructor to fit tree measures to the residuals, which is a favorable choice in the numericalexperiment provided in Inouye and Ravikumar (2018), and we set the shrinkage parameter (α) to
0.2%, 2%, and 20% of the sample size (20, 200, 2000) following Inouye and Ravikumar (2018). The
results are provided in Figure 5. Our algorithm substantially outperforms the DDD and the amount ofcomputational time needed to achieve desirable performance is substantially less under our method.

22

Figure 3: The true marginal densities with respect to X1 and X2 .

3.2

Performance Comparison with State-of-the-art Density Estimators

We evaluate the performance of our boosting algorithm using seven popular benchmark data sets
recorded in the University of California, Irvine (UCI) machine learning repository (Dua and Graff,2017). We preprocessed the four data sets (“POWER”, “GAS”, “HEPMASS”, and “MINIBOONE”)
following Papamakarios et al. (2017) with the code provided at https://github.com/gpapamak/maf
and the three data sets (“AReM”, “CASP”, and “BANK”) with the code provided at https://zenodo.
org/record/4560982#.Yh4k_OiZOCo. The density functions are estimated based on the training data
sets, and the performance is measured by the predictive score, i.e., the average log-density evaluatedat the held-out testing sets, one for each data set.
We compare our approach with three normalizing flow (NF) methods using deep neural networks toconstruct transforms, which represent the state-of-the-art for density estimation in machine learning:
MADE (Germain et al., 2015), Real NVP (Dinh et al., 2017), and MAF (Papamakarios et al., 2017),
and with the DDD algorithm. Their predictive scores are taken from Papamakarios et al. (2017) and(Liu et al., 2021), where the detailed settings of the NF models are provided. For the DDD algorithm,
we evaluate the performance under three possible values of the shrinkage parameter (α), 0.2%, 2%,
and 20% of the sample size, and show the best predictive scores. For our boosting algorithm, thevalues of c0 , γ—the tuning parameters for regularization—are set to (0.1, 0.0) and (0.1, 0.5).
A comparison of the predictive scores given by the boosting and the other NF methods is provided in Table 2, where our method is labeled as “boostPM” (which stands for “boosting probability
measures”). Our unsupervised tree boosting is overall competitive with the NF methods and evenshows the best predictive performance for two data sets (“POWER” and “AReM”). Appendix E also
provides a visual comparison of the training data sets and replicated data sets simulated from thefitted generative model, and it confirms that the distributional structures are successfully captured.
It should be noted that among the considered methods, ours is the only one that is not based onneural networks but is a combination of the tree-based learners, and therefore requires only a tiny
fraction of the computational cost to train.
Table 3 presents the computation time on the fourlarge data sets measured in the same single-core environment and shows our method is substantiallyfaster to train in large n settings. Additional tables are provided in Appendix E which show that23

Figure 4: The comparison between the single tree method and the boosting method under different
tuning parameter settings. The first row provides the predictive performance of the single tree method
and the boosting under different c0 (global shrinkage parameter) and fixed γ = 0.0 (penalty for small
nodes). The second row compares the predictive scores under different γ values and fixed c0 = 0.1.
The shown predictive scores are averages given under 30 data sets generated with different random
seeds.
the predictive performance of our boosting scarcely changes under different random seeds though thetree-constructing algorithm involves random splitting, and that the computation cost for simulation
is much smaller compared to the cost for constructing the ensemble.

3.3

Evaluating Variable ImportanceIn this section we demonstrate the use of the variable importance measure defined in Section 2.7.3.
We first carry out an experiment using artificial 10-dimensional distributions of (X1 , . . . , X10 ) withn = 10, 000 using the following scenarios:
Scenario(1): Xj ∼ Beta(21−j , 21−j ) for j = 1, . . . , 10.
Scenario(2): Five pairs of random variables (Ym,1 , Ym,2 ) (m = 1, . . . , 5), each of which follows a
2-dimensional Gaussian distribution
" # "
#!
0
1 ρ
N
,
,
0
ρ 1
with ρ = 0.1, 0.3, 0.5, 0.7 and 0.9 respectively, and(X2(m−1)+1 , X2m ) = (Φ(Ym,1 ), Φ(Ym,2 )),24

Figure 5:
The comparison between the boosting algorithm and the DDD (First row: predictivescores, second row: computation time) based on different 30 data sets. The points correspond to themeans, and the intervals are made by adding/subtracting the standard deviations.
where Φ(·) is a CDF of the standard Gaussian distribution.
Scenario(3): Each of four groups of variables (Y1 ), (Y2 , Y3 ), (Y4 , Y5 , Y6 ), and (Y7 , Y8 , Y9 , Y10 ) follows
the uni-variate/multi-variate Gaussian distribution in which the marginal distribution is thestandard Gaussian and the correlation is 0.9, and Xj = Φ(Yj ) for j = 1, . . . , 10.
The computed variable importance is displayed in Figure 6. These results show the tendency of Xj withlarger j obtaining larger importance, and it implies that a variable acquires higher importance whenthe marginal distribution is substantially different from the uniform and/or it is strongly dependent
on the other variables.
We next compute the variable importance on the MNIST handwritten digits data (LeCun et al.,1998). The gray scale images contained in this data consist of 28 × 28 = 784 pixels, and each pixel
takes integer values ranging from 0 (black) to 255 (white). We obtain the data with the read_mnist
function in the R package dslabs (Irizarry and Gill, 2021) and as in Papamakarios et al. (2017), scale
them into [0, 1]. The tuning parameters c0 and γ are both set to 0.1.
Recall that our notion of variable importance characterizes how each variable contributes to thedeviation from the uniform measure in the underlying sampling distribution. For the particular
25

POWERGASHEPMASSMINIBOONE

MADE
RealNVP
MAF
DDD
BoostPM (0.0)
BoostPM (0.5)

0.40 (0.01)
0.17 (0.01)
0.30 (0.01)-1.91 (0.03)
1.20 (0.01)
1.12 (0.01)
AReM

8.47 (0.02)
8.33 (0.14)
10.08 (0.02)-4.40 (0.31)
9.56 (0.02)
9.42 (0.02)
CASP−15.15 (0.02)-18.71 (0.02)−17.39 (0.02)-28.47 (0.03)-20.16 (0.02)-19.43 (0.02)
BANK−12.27 (0.47)-13.55 (0.49)−11.68 (0.44)-54.78 (1.03)-20.38 (0.46)-17.28 (0.47)

MADE
RealNVP
MAF
DDD
BoostPM(0.0)
BoostPM(0.5)

6.00 (0.11)
9.52 (0.18)
9.49 (0.17)
5.73 (0.16)
12.28 (0.16)
12.10 (0.15)

21.82(0.23)
26.81 (0.15)
27.61 (0.13)
7.77 (0.19)
21.84 (0.15)
22.23 (0.15)

14.97 (0.53)
26.33 (0.22)
20.09 (0.20)
9.33 (0.19)
36.47 (0.17)
36.34 (0.17)

Table 2: The comparison of the predictive scores (MADE, RealNVP, MAF, DDD, and the proposed
boosting under γ = 0.0, 0.5). The means of the estimated log-densities and the two standard deviations
are provided. The two top-performing methods for each dataset is indicated in bold font.
application of zip-code digit recognization, a pixel is more informative about the underlying digit if
it has large variation over the range of intensities. As such, the practical meaning of “importance”
in this particular application is the opposite to the statistical importance—it is exactly those pixels
with intensities spread out over large ranges (and thus more uniform) that are informative about theunderlying digit. As such, we want to emphasize the difference between the “practical importance”
and that of the “distributional importance” in terms of KL as we defined before.
The computed “distributional importance” obtained for the ten different digits is visualized on
the left of Figure 7, and a sample of handwriting of 0 is provided on the right. We can see that thetb

Figure 6:

The computed variable importance for the three simulation examples.

26

POWERGASHEPMASSMINIBOONEnd1,659,917
6

852,1748

315,123
21

29,55643

MADE
RealNVP
MAF
BoostPM(0.0)
BoostPM(0.5)

3.43
43.72
15.37
0.61
2.25

19.33
24.13
77.53
0.75
2.06

8.88
60.80
34.03
0.92
2.31

0.90
1.42
0.80
0.40
1.01

Table 3:
The training time (hours) for the benchmark data sets, measured in minutes (standard
error in parentheses) in a single-core AMD EPYC 7002 (2.50GHz) environment. We show the time
given under the optimal settings for MADE and RealNVP, and the time under 5 autoregressive layers
for MAF. Also provided are the sample size and dimensionality of the data sets.
pixels with relatively low “distributional importance” and hence high practical importance lie along theoutlines of the digits. Hence in this case the “distributional importance” on the left side characterizes“the average shapes” of the handwritten numbers.

Figure 7: Left: The importance of pixels (variables) computed for the 10 different digits. Yellow/purple colors correspond to low/high importance, which indicates large/small difference in the
handwriting styles found at the pixels. Right: A sample of handwriting of 0.

4

Concluding Remarks

We have proposed an unsupervised boosting method for learning multivariate probability measures byintroducing new notions of addition and residuals based on tree-CDF transforms, and demonstrated
how one can carry out density estimation and simulate from the fitted measure based on the output
27

of the algorithm. Given its similarity to classical boosting for regression and classification, we expect
other techniques for the boosting in such contexts, for example subsampling (Friedman, 2002), could
further improve the performance of our boosting method. Due to the limited space, we could not
exploit all possible techniques in supervised boosting for improving the performance, but we expect
many of them may be effective.

Acknowledgments
LM’s research is partly supported by National Science Foundation grants DMS-1749789 and DMS2013930. NA is partly supported by a fellowship from the Nakajima Foundation.

References
Naoki Awaya and Li Ma. Hidden markov Pólya trees for high-dimensional distributions. Journal of
the American Statistical Association, pages 1–13, 2022.
Vladimir I Bogachev. Measure theory, volume 1. Springer Science & Business Media, 2007.
Leo Breiman. Population theory for boosting ensembles. The Annals of Statistics, 32(1):1–11, 2004.
Jingyi Cui, Hanyuan Hang, Yisen Wang, and Zhouchen Lin. Gbht: Gradient boosting histogramtransform for density estimation. In International Conference on Machine Learning, pages 2233–
2243. PMLR, 2021.
Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio. Density estimation using Real NVP. Proceedings of the 5th International Conference on Learning Representations, 2017.
Dheeru Dua and Casey Graff. UCI machine learning repository, 2017. URL http://archive.ics.
uci.edu/ml.
Jerome H Friedman. Greedy function approximation: a gradient boosting machine. The Annals of
Statistics, pages 1189–1232, 2001.
Jerome H Friedman. Stochastic gradient boosting. Computational Statistics & Data Analysis, 38(4):367–378, 2002.
Mathieu Germain, Karol Gregor, Iain Murray, and Hugo Larochelle. MADE: Masked autoencoder for
distribution estimation. In International Conference on Machine Learning, pages 881–889. PMLR,
2015.
Subhashis Ghosal and Aad Van der Vaart. Fundamentals of nonparametric Bayesian inference, volume 44. Cambridge University Press, 2017.

28

Trevor Hastie, Robert Tibshirani, and Jerome Friedman. The elements of statistical learning: data
mining, inference, and prediction. Springer Science & Business Media, 2009.
Jingyu He and P Richard Hahn. Stochastic tree ensembles for regularized nonlinear regression. Journalof the American Statistical Association, (just-accepted):1–61, 2021.
Chin-Wei Huang, David Krueger, Alexandre Lacoste, and Aaron Courville. Neural autoregressive
flows. In International Conference on Machine Learning, pages 2078–2087. PMLR, 2018.
David Inouye and Pradeep Ravikumar. Deep density destructors. In International Conference on
Machine Learning, pages 2167–2175. PMLR, 2018.
Rafael A. Irizarry and Amy Gill. dslabs: Data Science Labs, 2021. URL https://CRAN.R-project.
org/package=dslabs. R package version 0.7.4.
Priyank Jaini, Kira A Selby, and Yaoliang Yu. Sum-of-squares polynomial flow. In International
Conference on Machine Learning, pages 3009–3018. PMLR, 2019.
Ivan Kobyzev, Simon Prince, and Marcus Brubaker. Normalizing flows: An introduction and review
of current methods. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2020.
Zhifeng Kong and Kamalika Chaudhuri. The expressive power of a class of normalizing flow models.
In International Conference on Artificial Intelligence and Statistics, pages 3599–3609. PMLR, 2020.
Michael Lavine. Some aspects of Pólya tree distributions for statistical modelling. The Annals of
Statistics, 20(3):1222–1235, 1992.
Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner. Gradient-based learning applied to
document recognition. Proceedings of the IEEE, 86(11):2278–2324, 1998.
Qiao Liu, Jiaze Xu, Rui Jiang, and Wing Hung Wong. Density estimation using deep generative
neural networks. Proceedings of the National Academy of Sciences, 118(15), 2021.
Luo Lu, Hui Jiang, and Wing H Wong. Multivariate density estimation by bayesian sequential partitioning. Journal of the American Statistical Association, 108(504):1402–1410, 2013.
Li Ma. Adaptive shrinkage in Pólya tree type models. Bayesian Analysis, 12(3):779–805, 2017.
Llew Mason, Jonathan Baxter, Peter Bartlett, and Marcus Frean. Boosting algorithms as gradientdescent. Advances in neural information processing systems, 12, 1999.
George Papamakarios. Neural density estimation and likelihood-free inference.
arXiv:1910.13233, 2019.

arXiv preprint

George Papamakarios, Theo Pavlakou, and Iain Murray. Masked autoregressive flow for density
estimation. Advances in neural information processing systems, 30, 2017.
29

George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan. Normalizing flows for probabilistic modeling and inference. Journal of Machine
Learning Research, 22(57):1–64, 2021.
Parikshit Ram and Alexander G Gray. Density estimation trees. In Proceedings of the 17th ACM
SIGKDD international conference on Knowledge discovery and data mining, pages 627–635, 2011.
Greg Ridgeway. Looking for lumps: Boosting and bagging for density estimation. Computational
Statistics & Data Analysis, 38(4):379–392, 2002.
Saharon Rosset and Eran Segal. Boosting density estimation. Advances in Neural Information Processing Systems, 15, 2002.
Wing H Wong and Li Ma. Optional Pólya tree and Bayesian inference. The Annals of Statistics, 38(3):1433–1459, 2010.

30

A

ProofsIn the following proofs, for a tree CDF G : (0, 1]d 7→ (0, 1]d and B ∈ B(Ω), the image is denoted byG(B) = {G(x) : x ∈ B},
and the same notation rule is applied for the inverse G−1 . The Lebesgue measure is denoted by µ.

A.1

Proof of Proposition 1

By the assumption on the full support, the CDFs are invertible.
Suppose X ∼ G1 ⊕ · · · ⊕ Gk . For x ∈ (0, 1], we have
−1
P (Gi ◦ · · · ◦ G1 (X) ≤ x) = P (X ≤ G−1
1 ◦ · · · ◦ Gi (x))−1
= Gk ◦ · · · ◦ G1 (G−1
1 ◦ · · · ◦ Gi (x))= Gk ◦ · · · ◦ Gi+1 (x),so Gi ◦ · · · ◦ G1 (X) ∼ Gi+1 ⊕ · · · ⊕ Gk . The converse can be shown by transforming P (X < x), whereGi ◦ · · · ◦ G1 (X), in the same way.

A.2

Proof of Proposition 2

The tree CDFs are already shown to be bijective in Section 2.2.1, so we only show the measurability
here.
Let E be a set of hyper-rectangles that are written in the form of(a1 , b1 ] × · · · × (ad , bd ]including the null set ∅. Since B((0, 1]d ) is the Borel σ-field, E ∈ B((0, 1]d ) holds for every E ∈ E. To
show the measurability of the tree CDF G, which is defined by the measure G ∈ PT , it suffices toshow that G−1 (E) ∈ B((0, 1]d ) for every E ∈ E since E generates B((0, 1]d ).
By the definition of G, the image G(A) (A ∈ L(T )) is also a hyper-rectangle included in E, andtheir collection {G(A) : A ∈ L(T )} forms a partition of (0, 1]d . Hence E is written as a union of
disjoint sets:[
E=(E ∩ G(A)),
A∈L(T )

where each E ∩ G(A) is a hyper-rectangle that belongs to E and a subset of G(A). Hence G−1 (E ∩G(A)) also belongs to E ⊂ B((0, 1]d ). Therefore, their finite union[
G−1 (E) =G−1 (E ∩ G(A))
A∈L(T )is also an element of B((0, 1]d ).
31

A.3

Proof of Theorem 2

To prove the first assertion, define G[r] for r = 1, . . . , R − 1 as
G[r] (x) =X

GA (x)1A (x),A∈Ar

. Then the first assertion is equivalent to that if X ∼ G, thenG[1] ◦ · · · ◦ G[R−1] (X) ∼ U nif ((0, 1]d ),(13)

which we prove here.
In the proof, we let X [R] = X and for r = 1, . . . , R − 1
X [r] = G[r] ◦ · · · ◦ G[R−1] (X).
Let G[r] denote a probability measure for X [r] . With these notations, we prove (13) by induction: We
show that, if for A ∈ N (T ), andG[r+1] (Al ) = G(Al ), G[r+1] (Ar ) = G(Ar ),G[r+1] (· | Al ) = µ(· | Al ), G[r+1] (· | Ar ) = µ(· | Ar ),(14)

G[r] (A) = G(A), G[r] (· | A) = µ(· | A).

(15)

then

The conditions in Equation 14 holds if r = R − 1 because G[R] = G and G ∈ PT , and the statementin Equation 15 being true for r = 1 implies that G[1] ◦ · · · ◦ G[R] (·) = µ(·), which is equivalent toEquation 13.
Assume Equation 14 holds for some r. By the definition, G[r] is bijective, and G[r],−1 (A) = A for
every A ∈ Ar . Then for X ∼ G[r] and A ∈ Ar , we have,G[r] (A) = P (X [r] ∈ A)
= P (G[r],−1 (X [r] ) ∈ G[r],−1 (A)) = P (X [r+1] ∈ A)
= G[r+1] (Al ) + G[r+1] (Ar ) = G(A).
[r]

[r]

Hence the first equation in Equation 15 holds. To prove the second equation, let X [r] = (X1 , . . . , Xd ).
Then, we show that for zj ∈ (aj , bj ],
[r]

[r]

P (X1 ∈ (a1 , z1 ], . . . , Xd ∈ (ad , zd ] | X [r] ∈ A) =32dY
zj − ajb − ajj=1 j(16)

holds. The probability in the left hand side can be written as follows:
[r]

[r]

P (X1 ∈ (a1 , z1 ], . . . , Xd ∈ (ad , zd ] | X [r] ∈ A)
[r]

[r]

= P (X1 ∈ (a1 , z1 ], . . . , Xd ∈ (ad , zd ] | X [r+1] ∈ A)=

P (X [r+1] ∈ Al )[r][r]
P (X1 ∈ (a1 , z1 ], . . . , Xd ∈ (ad , zd ] | X [r+1] ∈ Al )P (X [r+1] ∈ A)+

P (X [r+1] ∈ Ar )[r][r]
P (X1 ∈ (a1 , z1 ], . . . , Xd ∈ (ad , zd ] | X [r+1] ∈ Ar )P (X [r+1] ∈ A)
[r]

[r]

[r]

[r]

= G(Al | A)P (X1 ∈ (a1 , z1 ], . . . , Xd ∈ (ad , zd ] | X [r+1] ∈ Al )+ G(Ar | A)P (X1 ∈ (a1 , z1 ], . . . , Xd ∈ (ad , zd ] | X [r+1] ∈ Ar ).

(17)
[r]

Let A be divided in the j ∗ th dimension. By the definition of GA , for j ̸= j ∗ , Xj

∈ (aj , zj ] ⇐⇒

[r+1]
Xj
∈ (aj , zj ]. For j ∗ , because GA,j ∗ (·) is strictly increasing,

i

 
[r]
[r]
−1
−1
∗)
∗ ), G(z(aX∈
G
Xj ∗ ∈ (aj ∗ , zj ∗ ] ⇐⇒ G−1
∗
∗
∗
∗
j
j
j
A,j
A,j
A,j
[r+1]

⇐⇒ Xj ∗∈ (aj ∗ , yj ∗ ] ,∗∗∗∗where yj ∗ = G−1
A,j ∗ (zj ). The expression of yj changes depending on whether yj ≤ cj or not, wherecj ∗ is a partition point at which A is divided. We first assume that yj ∗ ≤ cj ∗ . In this case, the second[r+1]
term in Equation 17 is 0 because Xj ∗∈ (aj ∗ , yj ∗ ] does not happen if X [r+1] ∈ Ar . Also, by thedefinition of GA,j ∗ ,bj ∗ − aj ∗zj ∗ − aj ∗G(Al | A))== G(Al | A))yj ∗ − a j ∗µ(Al | A)cj ∗ − aj ∗∗∗∗∗yj − a jzj − aj1⇐⇒=.
cj ∗ − aj ∗G(Al | A)) bj ∗ − aj ∗Therefore, it follows that
[r]

[r]

P (X1 ∈ (a1 , z1 ], . . . , Xd ∈ (ad , zd ] | X [r] ∈ A)
Y z −a y ∗ −a ∗jjjj
= G(Ar | A)
 ∗ bj − aj  cj ∗ − aj ∗j̸=jY z −a zj ∗ − aj ∗1jj
= G(Ar | A)
 ∗ bj − aj  G(Al | A)) bj ∗ − aj ∗
j̸=j

=

d
Y
j=1

zj − a j.
bj − aj

We can prove (16) for the case of yj ∗ > cj ∗ in the same way.

33

To prove the second result, let U ∼ Unif((0, 1]d ). The multi-scale CDF G is bijective (Proposition2), so we obtain for B ∈ B(Ω)
P (G−1 (U ) ∈ B) = P (U ∈ G(B))= µ(G(B))= G(B),
where the last line follows Theorem 1. (Note that the proof of Theorem 1 only uses the first result ofTheorem 2). Therefore G−1 (U ) ∼ G.

A.4

Proof of Theorem 1

Let X ∼ G. By the first result of Theorem 2, we obtain
G(B) = P (X ∈ B)
= P (G(X) ∈ {G(x) : x ∈ B})= µ({G(x) : x ∈ B}).

A.5

Proof of Lemma 1

We only need to check the countable additivity. By Proposition 2, Gk ◦ · · · ◦ G1 is bijective. Hence,for disjoint sets Al ∈ B(Ω) (l ∈ N), it follows that!
[[Gk ◦ · · · ◦ G1 (Al ) .
Gk ◦ · · · ◦ G1
Al =
l

l

Because {Gk ◦ · · · ◦ G1 (Al )}l=1,2,... are disjoint, this result implies that!
!
[[
Al = µGk ◦ · · · ◦ G1 (Al )
Fk
l

l

=

X

=

X

µ(Gk ◦ · · · ◦ G1 (Al ))l

F (Al ) .

l

A.6

Proof of Proposition 3

First we suppose X ∼ G1 ⊕· · ·⊕Gk . By Proposition 2, G1 , . . . , Gk are all bijective, so for i = 1, . . . , k,
the composition
Fi = Gi ◦ · · · ◦ G1

34

is also bijective. For B ∈ B(Ω), we obtain
P (Fi (X) ∈ B) = P (X ∈ F−1
i (B))= Fk ((F−1i (B)))= µ(Fk (F−1i (B)))= µ(Gk ◦ · · · ◦ Gi+1 (B)).
Hence Fi (X) ∼ Gi+1 ⊕ · · · ⊕ Gk . Showing the converse is now straightforward.

A.7

Proof of Proposition 4

We first show the following lemma, which implies that the conditional distributions Fk−1 and Fk(k = 2, . . . , K) are the same on subsets in a partition defined by Tk and Fk−1 = Gk−1 ◦ · · · ◦ G1 .
′Lemma 4. Let L(Tk ) be a set of the terminal nodes in Tk . Also, we let A′ ∈ L(Tk ) and A = F−1
k−1 (A ).
Then for any B ⊂ A, Fk−1 (B | A) = Fk (B | A).

(Proof) For B ⊂ A, by Theorem 1, we have
Fk (B) = µ(Fk (A))= µ(Gk (Fk−1 (B)))= Gk (Fk−1 (B)).

(18)

Since B ⊂ A, Fk−1 (B) ⊂ Fk−1 (A) = A′ . Hence Fk (B) is further rewritten as follows
Fk (B) = Gk (A′ ∩ Fk−1 (B))= Gk (A′ )Gk (Fk−1 (B) | A′ )= Gk (A′ )µ(Fk−1 (B) | A′ )µ(Fk−1 (B))µ(A′ )
Fk−1 (B)= Gk (A′ ).
µ(A′ )= Gk (A′ )By the definition of A, Fk−1 (A) = µ(Fk−1 (A)) = µ(A′ ), and by replacing B with A in (18), we have
Fk (A) = Gk (Fk−1 (A)) = Gk (A′ ).
Therefore, we obtain
Fk (B | A) =Fk (B)Fk (A)

Fk−1 (B) 1µ(A′ ) Fk (A)Fk−1 (B)=Fk−1 (A)= Gk (A′ )= Fk−1 (B | A).
35

′′(Proof of Proposition 4) Let Lk = {F−1k−1 (A ) : A ∈ L(Tk )}. By Lemma 4, the conditional distributions
Fk−1 (· | A) and Fk (· | A) are the same for A ∈ Lk . Hence the density functions of Fk−1 and Fk denotedby fk−1 and fk are expressed asX
fk−1 (x) =Fk−1 (A)fk−1 (x | A)1A (x),A∈Lkfk (x) =X

Fk (A)fk−1 (x | A)1A (x),A∈Lkwhere 1 is the indicator function. Fix x ∈ (0, 1]d and let Ak ∈ Lk such that x ∈ Ak and A′k = Fk−1 (Ak ).
By Theorem 1, we haveFk−1 (Ak ) = µ(Fk−1 (Ak )) = µ(A′k ),Fk (Ak ) = µ(Fk (Ak ))= µ(Gk ◦ Fk−1 (Ak ))= Gk (Fk−1 (Ak ))= Gk (A′k ).
Hence, we have
Gk (A′k )fk (x)Fk (Ak )==.
fk−1 (x)Fk−1 (Ak )µ(A′k )

(19)

Since x ∈ Ak ,Fk−1 (x) = Gk−1 ◦ · · · ◦ G1 (x) ∈ Fk−1 (Ak ) = A′k .
Thus the density ratio in (19) is rewritten as
fk (x)= Gk (A′k )µ(Gk−1 ◦ · · · ◦ G1 (x) | A′k )fk−1 (x)= Gk (A′k )gk (Gk−1 ◦ · · · ◦ G1 (x) | A′k )= gk (Gk−1 ◦ · · · ◦ G1 (x)),
where the second equation follows that A′k ∈ L(Tk ). Because the discussion above holds for k =2, . . . , K, we obtain the following expression
KY
fk (x)fK (x) = f1 (x)fk−1 (x)= f1 (x)k=2
KY

gk (Gk−1 ◦ · · · ◦ G1 (x)).

k=2

36

A.8

Proof of Proposition 6

The associativity clearly holds, and showing the existence of the identity element is also straightforward
because the identity transformation, which is a tree-CDF of the uniform distribution, is included inG.
Proving the existence of an inverse element for every element in G is done by showing that an inverse
function of a local-move function, which is an essential component of a tree-CDF, is also a local-movefunction. To this end we use the local-move function GA : A 7→ A defined in Section 2.2.1 and thesame notations.
Let c̃j ∗ = aj ∗ + G(Al | A)(bj ∗ − aj ∗ ) and Ãl and Ãr be a pair of nodes that we obtain by dividing
A at c̃j ∗ in the j ∗ th dimension. Then we define a function G̃A : A 7→ A such that for any x ∈ A,
G̃A (x) = (G̃A,1 (x), . . . , G̃A,d (x)) where G̃A,j (x) = x for all j ̸= j ∗ , andG̃A,j ∗ (x) − aj ∗µ(Al |A)=xj ∗ − aj ∗µ(Ãl |A)for x ∈ Ãlandbj ∗ − G̃A,j ∗ (x)µ(Ar |A)=bj ∗ − xj ∗µ(Ãr |A)for x ∈ Ãr .

From this definition G̃A works as a local-move function for the conditional distribution that assigns
the probability µ(Al |A) to Ãl . It is straightforward to see that this transformation G̃A is identical toG−1
A (see the expression of the inverse function provided in Section 2.3.2).

A.9

Proof of Lemma 2

By the definition of the KL divergence, we have
Z
Z
∗
∗
∗
KL(F ||F ) = log f dF − log f dF ∗ .
By Proposition 4, the second term
ZR

log f F ∗ is decomposed aslog f dF ∗ =K Z
Xlog gk (Gk−1 ◦ · · · ◦ G1 (x))dF ∗ (x).

k=1

By the change-of-variable formula (e.g., Theorem 3.6.1 in Bogachev (2007)), the right hand side canbe written in a form of integration with respect to F̃k ,ZZlog gk (Gk−1 ◦ · · · ◦ G1 (x))dF ∗ (x) = log gk (x)dF̃k (x).
Since the measure F ∗ is absolutely continuous with respect to the Lebesgue measure µ, by the definitionof F̃k , that is,−1dF̃k (B) = F ∗ (G−1
1 ◦ · · · ◦ Gk−1 (B)) for all B ∈ B((0, 1] ),37

F̃k is also absolutely continuous with respect to µ. Hence F̃k admits the density function denoted byf˜k (x), and the right hand side is further rewritten as follows
ZZf˜k (x)f˜k (x)dF̃k (x) − log
dF̃k (x)logµ(x)
gk (x)

Z
log gk (x)dF̃k (x) == KL(F̃k ||µ) − KL(F̃k ||Gk ).

A.10

Proof of Proposition 5

Since Tk is a finite tree, the log of gk , which is piece-wise constant on TK , is written asX

log gk (x) =log

A∈L(Tk )Gk (A)
1A (x) for x ∈ (0, 1]d .
µ(A)

(n)

Hence the improvement Dk (Gk ) is rewritten as follows:
(n)

(n)

X

Dk (Gk ) =F̃k (A) log

A∈L(Tk )Gk (A)
µ(A)
(n)

(n)

X

=

F̃k (A) log

A∈L(Tk )F̃k (A)
−
µ(A)

(n)

X
A∈L(Tk )(n)F̃k (A) log

F̃k (A).
Gk (A)

Because the second term in the bottom line takes a form of KL divergence defined for the two discrete(n)distributions, it is minimized if Gk (A) = F̃k (A) for all A ∈ L(Tk ). Under this Gk , since the second
term is 0, the improvement is maximized if Tk satisfies the condition provided in Proposition 5.

A.11

Proof of Proposition 7

In this proof, we suppose the learning rate c(A) is independent to a node A for simplicity. For every
leaf node A ∈ L(Tk ), there is a sequence of nodes {BA,r }R
r=1 such that BA,r belongs to the rth levelof Tk , and(0, 1]d = BA,1 ⊃ BA,2 ⊃ · · · ⊃ BA,r = A.
(n)With such sequences, based on the discussion in Appendix A.10 , the improvement Dk (Gk ) is
decomposed as(n)

Dk (Gk ) =X(n)F̃k (A) log

A∈L(Tk )Gk (A)µ(A)
Gk (BA,2 |BA,1 )Gk (BA,R |BA,R−1 )=log+ · · · + logµ(BA,2 |BA,1 )µ(BA,R |BA,R−1 )A∈L(Tk )X
Gk (Al | A)
Gk (Ar | A)
(n)
(n)
=
F̃k (Al ) log+ F̃k (Ar ) log.
µ(Al | A)µ(Ar | A)X(n)F̃k (A)A∈N (Tk )

38

(n)

For the bottom line the summand is 0 if F̃k (A) = 0. Otherwise, the conditional probabilities(n)(n)
F̃k (Al | A) and F̃k (Ar | A) are defined. In such a case, by the definition of Gk (Al | A),#"(n)(1 − c)µ(Al | A) + cF̃k (Al | A)
Gk (Al | A)
log
= log
µ(Al | A)µ(Al | A)(n)≥ (1 − c) log 1 + c log

F̃k (Al | A)µ(Al | A)(n)= c logF̃k (Al | A),µ(Al | A)

where the second line follows the Jensen’s inequality. The same result holds for Ar . Hence,
Gk (Al | A)
Gk (Ar | A)(n)
+ F̃k (Ar ) logµ(Al | A)µ(Ar | A)
#
"
(n)
(n)
F̃k (Ar | A)F̃k (Al | A)
(n)
(n)
(n)
+ F̃k (Ar | A) log,≥ cF̃k (A) F̃k (Al | A) logµ(Al | A)µ(Ar | A)(n)F̃k (Al ) log

where the sum inside of the brackets is the KL divergence for the two Bernoulli distributions and thus(n)
non-negative. Therefore, the improvement Dk (Gk ) is non-negative. Additionally, the last inequality(n)(n)
is strict if and only if F̃k (Al | A) = µ(Ar | A) and so Dk (Gk ) is positive.

B

Expressive Power of the Tree Ensemble

In this section, we provide theoretical results on the expressive power of the tree ensemble with thefinal goal of proving Theorem 3.

B.1

Preparations

We introduce the following notations:
1. Let T L be a collection of dyadic trees with axis-aligned boundaries with at most L maximum
resolution. When L = d, T d is a set of trees that can be formed under Assumption 1. We notethat as implied in the following proofs, T L can a set of trees that have at least one node reachthe Lth while the other leaf nodes belong to the shallower levels.
2. For a tree T ∈ T L , a set PT denotes a collection of probability measures conditionally uniform
on T such that
G(· | A) = µ(· | A) and G(A) > 0(20)

for every terminal node A ∈ T . A collection of such tree measures are denoted by G0L , that is,
G0L = {G : G ∈ PT for some T ∈ T L }.
39

For a measure G ∈ G0L defined on a tree T ∈ T L , we can define a tree-CDF as in Section 2,
which is denoted by G. We define a set G L
0 as a collection of such tree CDFs, namely,LGL
0 = {G : G is a tree CDF of G ∈ G0 }.

3. Let G L denote a set of finite composition of tree CDFs, that is,
G L = {GK ◦ · · · ◦ G1 : K ∈ N and for k = 1, . . . , K, Gk ∈ G L
0 },and define G L as a collection of probability measures defined by such finite compositions, thatis,
G L = {µ(G(·)) : G ∈ G L }.
Hence G L includes all measures that can be expressed in the form of ensemble G1 ⊕ · · · ⊕ GK .
We also need to review the definition of push-forward measures because this notation is closely
related to the operation of residualization. Let φ be a mapping Ω 7→ Ω and H be a probability
measure. Then the push-forward of H is defined in the following form:
φ#H(B) = H(φ−1 (B)) for B ∈ B(Ω).
The following lemma establishes a connection between the ensemble measure and the push-forward
measures.
Lemma 5. For a probability measure F , F ∈ G L holds if and only if there exists a mapping G ∈ G L
such that G#F = µ.
(Proof) Suppose F ∈ G L . Then there exists a mapping G ∈ G L such that
F (B) = µ(G(B)) for B ∈ B(Ω).
From Proposition 2, G is bijective. Hence for B ∈ B(Ω), we haveµ(B) = µ(G ◦ G−1 (B)) = F (G−1 (B)),
so G#F = µ. The necessity can be shown in the same way.
In the rest of the section, we first discuss the expressive power of the tree ensemble for the univariate cases and next generalize the result for the multi-variate cases. After that, this result is used
to prove Theorem 3.

40

B.2

Uni-variate Cases

The following proposition shows that any distribution with piece-wise constant and positive densities
can be represented in the form of tree ensemble.
Proposition 8. Let F be a probability measure that admits the piece-wise constant density f with thefollowing formIX
f (x) =βi 1(ci−1 ,ci ] ,i=1

where βi > 0 for i = 1, . . . , I and
0 = c0 < c1 < · · · < cI = 1.
Then, if L ≥ 2, F ∈ G L holds.
(Proof) We first show the existence of a tree CDF G1 ∈ G 20 such that the push-forward measure
G1 #F has a density f1 with the following form
f1 (x) =I−1
Xβ̃i 1(c̃i−1 ,c̃i ] ,

(21)

i=1

where β̃i > 0 for i = 1, . . . , I − 1 and 0 = c̃0 < c̃1 < · · · < c̃I−1 = 1.
Let α ∈ (0, 1) be a constant that satisfies
β2 1 − c1
1−α
=
.
α
β1 c1
Then define a measure G1 ∈ G02 such that
G1 ((0, c1 ]) = α, G1 ((c1 , 1]) = 1 − αand G1 is conditionally uniform on (0, c1 ] and (c1 , 1]. Let G1 be G1 ’s tree CDF and F1 = G1 #F bea probability measure with the density f1 . For x ∈ (0, α], we have−1
F1 ((0, x]) = F (G−1
1 ((0, x])) = F ((0, G1 (x)]) =Hence, by the chain rule, the density at this x is written as
f1 (x) =

 c1
c1
f G−1
β1 .
1 (x) =
αα

Similarly, the density at x ∈ (α, 1] is written as
f1 (x) =
1 − c1f G−1
1 (x) .
1−α41

Z G−1 (x)
1

f dµ.
0

Let c̃i = G1 (ci+1 ) for i = 1, . . . , I − 1. By this definition, α < c̃1 , and the density of f1 at x ∈ (α, c̃1 ]satisfies
f1 (x) =1 − c1c1
β2 = β1 ,1−ααwhere the second equation follows the definition of α. Hence f1 is constant on (0, c̃1 ]. Moreover, thedensity on (c̃i−1 , c̃i ] for i = 2, . . . , I − 1 is (1 − c1 )/(1 − α)βi−1 so constant. Therefore the density f1
is written in the form of Equation 21.
By using the same logic for the rest of the I − 2 discontinuous points, we can define tree CDFs
G2 , · · · , GI−1 that connect the densities at these points one by one. Hence the measure (GI−1 ◦ · · · ◦G1 )#F has a constant density and thus is the uniform measure µ.

B.3

Multi-variate Cases

In this section, we prove the following proposition that is a generalization of Proposition 8.
I

Proposition 9. For j = 1, . . . , d, let {cj,ij }ijj=1 be a sequence such that0 = cj,0 < cj,1 < · · · < cj,Ij = 1,
and L = {Ai1 ,...,id }i1 ,...,id be a partition of the sample space (0, 1]d that consists of rectangles written as
Ai1 ,...,id = (c1,i1 −1 , c1,i1 ] × · · · × (cd,id −1 , cd,id ].
If a probability measure F is piecewise uniform on L and written as
F (B) =

X
i1 ,...,ij

ai1 ,...,idµ(B ∩ Ai1 ,...,id ), for B ∈ B((0, 1]d ),µ(Ai1 ,...,id )where ai1 ,...,id > 0, then for L ≥ d + 1, there is a mapping G ∈ G L such that G#F = µ andthus F ∈ G L . In addition, we can choose G so that for every pair of indices (i1 , . . . , id ), the image
G(Ai1 ,...,id ) is a rectangle written as(G(c1,i1 −1 ), G(c1,i1 )] × · · · × (G(cd,id −1 ), G(cd,id )].
(Proof) We use induction: We assume that the statement of Proposition 9 is valid for the 1,2,. . . ,(d-1)dimensional cases. Because in this proof we handle measures and transformation defined in different
dimensional spaces, the sets G L and G L defined for the j-dimensional space are denoted by G L,d and
G L,d , respectively.
Inside of the induction, we also assume that for some l ∈ {1, . . . , Id − 1}, there are mappings
G1 , . . . , Gl ∈ G L,d such that a probability measure Fl := (Gl ◦ · · · ◦ G1 )#F is a piecewise uniform
probability measure written as, for B ∈ B((0, 1]d ),Fl (B) =

l
X
i=1

Ci

Id
X
µ(B ∩ (0, 1]d−1 × (cd,i−1 , cd,i ])+µ((0, 1]d−1 × (cd,i−1 , cd,i ])(l)X

id =l+1 i1 ,...,id−1

42

(l)

ai1 ,...,idµ(B ∩ Ai1 ,...,id )(l)µ(Ai1 ,...,id )

,

(l)

(l)

where Ci > 0 and ai1 ,...,id > 0 for all indices. Also, for the second term, Ai1 ,...,id is a rectangular
written as

i

i
(l)
(l)
(l)
(l)
(l)
Ai1 ,...,id = c1,i1 −1 , c1,i1 × · · · × cd−1,id−1 −1 , cd−1,id−1
× (cd,Id −1 , cd,Id ],(l) Ijwhere for j = 1, . . . , d − 1, {cj,i }i=1
is a sequence such that

(l)

(l)

(l)

0 = cj,1 < cj,2 < · · · < cj,Ij = 1.
(We note that this sequence’s length can be different from “Ij ” provided in Proposition 9 but to avoid
an excessive number of indices, we use Ij here because its size does not affect the logic provided inthis proof.) Under this assumption, we show that there is a measure (“Fl+1 ”) that has the same form
for l + 1.
Define a d − 1-dimensional probability measure F̂l+1
F̂l+1 =(l)
X a(l)
i1 ,...,id−1 ,l+1 µd−1 (B ∩ Âi1 ,...,id−1 )i1 ,...,id

Cl+1

(l)

for B ∈ B((0, 1]d−1 ),µd−1 (Âi1 ,...,id−1 )where Cl+1 is the normalizing constant, µd−1 is the Lebesgue measure defined for the d−1-dimensional(l)sample space, and Âi1 ,...,id−1 is a set written as

i

i
(l)
(l)
(l)
(l)
(l)
Âi1 ,...,id−1 = c1,i1 −1 , c1,i1 × · · · × cd−1,id−1 −1 , cd−1,id−1 .

(22)Because F̂ is a piecewise uniform measure defined on the partition that consists of hyper-rectangles, by
the assumption we set for the induction, there is a mapping Ĝl+1 ∈ G L−1,d−1 such that Ĝl+1 #F̂l+1 =µd−1 . With this mapping, we define a mapping Gl+1 : (0, 1]d 7→ (0, 1]d such that for x = (x1 , . . . , xd ) ∈(0, 1]d ,
Gl+1 (x) = Ĝl+1 (x1 , . . . , xd−1 ), xd
if xd ∈ (cd,l , 1] and otherwise Gl+1 (x) = x. The mapping Gl+1 moves points only in (0, 1]d−1 × (ci,l , 1],
which is a node one can obtain by dividing the sample space only once, according to Ĝl+1 , which isa mapping that is a composition of tree CDFs based on trees with L − 1 leaf nodes. Hence Gl+1 is acomposition of tree CDFs defined on trees with (L − 1) + 1 = L leaf nodes, so we have Gl+1 ∈ G L,d .
With this mapping, we define a measure Fl+1 = Gl+1 #Fl .
Fix a pair of indices (i1 , . . . , id ) and let Bd−1 ∈ B((0, 1]d−1 ) and B1 ∈ B((0, 1]) be measurable setssuch that(l)
Bd−1 × B1 ∈ Ĝl+1 (Âi1 ,...,id−1 ) × (cd,id −1 , cd,id ] .
If id ≤ l, by the definition of Fl+1 and Gl+1 ,Fl+1 (Bd−1 × B1 ) = Fl (G−1l+1 (Bd−1 × B1 )) = Fl (Bd−1 × B1 ).
43

(l)

On the other hand, if id ≥ l + 1, since Fl is conditionally uniform on Ai1 ,...,id ,−1
Fl+1 (Bd−1 × B1 ) = Fl (G−1l+1 (Bd−1 × B1 )) = Fl (Ĝl+1 (Bd−1 ) × B1 )µ(Ĝ−1(l)l+1 (Bd−1 ) × B1 )= ai1 ,...,id(l)µ(Ai1 ,...,id )µd−1 (Ĝ−1(l)
l+1 (Bd−1 ))µ1 (B1 )= ai1 ,...,id,(l)µd−1 (Âi1 ,...,id−1 )µ1 ((cd,id −1 , cd,id ])where µ1 is the Lebesgue measure defined for the 1-dimensional sample space. For such id , by thedefinition of F̂l+1 and Ĝl+1 ,(l)µd−1 (Bd−1 ) = F̂l+1 (Ĝ−1
l+1 (Bd−1 )) =ai1 ,...,id−1 ,l+1 µd−1 (Ĝ−1
l+1 (Bd−1 ))
Cl+1

(l)

,

µd−1 (Âi1 ,...,id−1 )from which we obtain(l)ai ,...,i ,i µd−1 (Bd−1 )µ1 (B1 )Fl+1 (Bd−1 × B1 ) = Cl+1 (l)1 d−1 dai1 ,...,id−1 ,l+1 µ1 ((cd,id −1 , cd,id ])µ(Bd−1 ×B1 )
(id = l + 1),
Cl+1 µ (0,1]d−1 ×(c(d,id −1 ,cd,id ])=(l+1) µ(Bd−1 ×B1 )
(id > l + 1),
ai1 ,...,id
(l+1)

(23)

µ(Ãi ,...,i )1dwhere, for id > l + 1,
(l+1)

(l)

Ãi1 ,...,id = Ĝl+1 (Âi1 ,...,id ) × (cd,id −1 , cd,id ].
and
(l)


ai1 ,...,id−1 ,id
(l+1)
(l)
ai1 ,...,id = Cl+1 (l)µd−1 Ĝl+1 (Âi1 ,...,id ) .
ai1 ,...,id−1 ,l+1

Because B((0, 1]d−1 ) × B((0, 1]) generates B((0, 1]d−1 ), from the discussion provided above, Fl+1 ispiecewise-uniform on a partition that consists of(0, 1]d−1 × (cd,id −1 , cd,id ] (id ≤ l + 1)
(l+1)

and Ãi1 ,...,id (id > l + 1), and this partition is denoted by P (l+1) . Note that by the definition of Ĝl+1 ,and Proposition 9, which we assume holds for (d − 1)-dimensional cases, the image of the hyper(l)(l+1)
rectangle of Âi1 ,...,id−1 under Ĝl+1 and Ãi1 ,...,id are a hyper-rectangle in the (d − 1)-dimensional space
and the d-dimensional space, respectively.
The following lemma states that the partition structure P (l+1) has a finer partition that has the“checker-board” form, as shown in the next lemma.
44Lemma 6. Let {Di }Ii=1 is a partition of the sample space (0, 1]d such that every Di is a hyper-rectangle.
Then, there are sequences {ej,ij }ij (j = 1, 2, . . . , d) such that0 = ej,0 < ej,1 < · · · < ej,Ij = 1
and a partition {Ei1 ,...,id }i1 ,...,id defined as
Ei1 ,...,id = (e1,i1 −1 , e1,i1 ] × · · · × (ed,id −1 , ed,id ]such that every Di is a finite union of elements of {Ei1 ,...,id }i1 ,...,id .
Its proof is straightforward because we only need to “extend” the boundaries between the rectangles{Di }Ii=1 . By applying this lemma to the partition P (l+1) , it follows that there are finite sequences(l+1){cij }ij (j = 1, . . . , d − 1) such that a checkerboard-like partition consisting of the following type ofrectangles
(l+1)

(l+1)

(l+1)

(l+1)

(l+1)

Ai1 ,...,id :=(c1,i1 −1 , c1,i1 ] × · · · × (cd−1,id−1 −1 , cd−1,id−1 ] × (cd,id −1 , cd,id ]is finer than P (l+1) . With this partition, the measure Fl+1 is written as for B ∈ B((0, 1]d )Fl+1 (B) =

l+1
X
i=1

(l)
µ(B ∩ Ai1 ,...,id )(l)ai1 ,...,id.
(l)µ(Ai1 ,...,id )
i=l+2 i1 ,...,id−1IdXµ(B ∩ (0, 1]d−1 × (cd,i−1 , cd,i ])Ci+µ((0, 1]d−1 × (cd,i−1 , cd,i ])X

Because this result holds for l = 1, . . . , Id − 1, there exists a sequence of mappings
G1 , . . . , GId ∈ G L,d such that a push-forward measure H := (GId ◦ · · · ◦ G1 )#F has a form
H(B) =

Id
X
i=1

Ci

µ(B ∩ (0, 1]d−1 × (cd,i−1 , cd,i ]).
µ((0, 1]d−1 × (cd,i−1 , cd,i ])Define an one-dimensional probability measure Ĥ as follows
Ĥ(B1 ) =

Id
X
i=1

Ci

µ1 (B1 ∩ (cd,i−1 , cd,i ])for B1 ∈ B((0, 1]).
µ1 ((cd,i−1 , cd,i ])Then, by the assumption (or Proposition 8), there exists a mapping Ĝ0 ∈ G 2,1 such that Ĝ0 #Ĥ = µ1 .
With this mapping, we define another mapping G0 : (0, 1] 7→ (0, 1] such that for x = (x1 , . . . , xd ),G0 (x) = (x1 , . . . , xd−1 , Ĝ0 (xd )).
This mapping moves input points only in the dth dimension according to Ĝ0 so it is written as acomposition of tree CDFs defined on trees with 2 terminal node and thus an element of G L,d . Hence
G0 ∈ G L,d . Fix i ∈ {1, . . . , Id }. For a measurable set Bd−1 × B1 such that
Bd−1 × B1 ∈ B((0, 1]d−1 ) × Ĝ0 ((cd,i−1 , cd,i ]),45

because H is piecewise uniform, we have
−1
G0 #H(Bd−1 × B1 ) = H(G−1
0 (Bd−1 × B1 )) = H(Bd−1 × Ĝ0 (B1 ))= Ciµ(Bd−1 × Ĝ−1
0 (B1 ))d−1µ((0, 1]× (cd,i−1 , cd,i ])= Ciµd−1 (Bd−1 )µ1 (Ĝ−1
0 (B1 )).
µ1 ((cd,i−1 , cd,i ])On the other hand,µ1 (B1 ) = Ĥ(Ĝ−1
0 (B1 )) = C1

µ1 (Ĝ−1
0 (B1 )).
µ1 ((cd,i−1 , cd,i ])Hence, we obtain
G0 #H(Bd−1 × B1 ) = µd−1 (Bd−1 )µ1 (B1 ) = µ(Bd−1 × B1 ).
Therefore, we conclude thatG0 #H = (G0 ◦ GId ◦ · · · ◦ G1 )#H = µ.
The result of Proposition 9 can be described in a simplified form as in the next corollary. This
proof immediately follows Proposition 9 and Lemma 6.
Corollary 1. Let {Ei }Ii=1 is a partition of the sample space (0, 1]d such that Ei is a rectangle with aform
Ei = (ai,1 , bi,1 ] × · · · × (ai,d , bi,d ],and F be a piecewise uniform probability measure defined on the partition:
F (B) =

I
X

βi

i=1

where βi

B.4µ(B ∩ Ei )for B ∈ B((0, 1]d ),µ(Ei )> 0. Then F ∈ G L for L ≥ d + 1.

Proof of Theorem 3

We finally provide the proof of Theorem 3, which can be obtained by adding minor modifications to
the proof of Theorem 4 in Wong and Ma (2010).
Let f ∗ denote F ∗ ’s density function, and we first assume that f ∗ is uniformly continuous. For
ϵ > 0, there exists ϵ′ > 0 such that log(1 + ϵ′ ) < ϵ. Since the function f ∗ is uniformly continuous,
there exists δ > 0 such that|x − y| < δ ⇒ |f ∗ (x) − f ∗ (y)| < ϵ′ .
Let {Ei }Ii=1 is a partition of the sample space (0, 1]d such that Ei has a rectangle shape and diam(Ei ) <δ. Define a function g̃ as

I 
X
g̃ =
sup f ∗ (x) 1Ei (x) for x ∈ (0, 1]d .
i=1

x∈Ei

46

Let C =Rg̃dµ. Because g̃(x) ≥ f ∗ (x) for x ∈ (0, 1]d , we have C ≥ 1 andZ
0≤C −1=(g̃ − f ∗ )dµ =

I Z
X
i=1

≤

I Z
X
i=1

(g̃(x) − f ∗ (x))dµ

Eiϵ′ dµ = ϵ′ .

EiDefine a density function g := g̃/C. The corresponding probability measure G is an element of G L by
Corollary 1. Hence, for the two measures F ∗ and G, we can bound the KL divergence as follows
Z
Z
Z
f∗
f∗
∗
∗
KL(F ||G) = f log dµ = f log dµ + f ∗ log Cdµgg̃≤ log C ≤ log(1 + ϵ′ ) < ϵ.
We next consider the general case, where we assume f ∗ ≤ M for some M > 0. By Lusin’s theorem,
for any ϵ̃ > 0, there exits a closed set B such that µ(B c ) < ϵ̃ and f ∗ is uniformly continuous on B.
Using this fact, we modify the first discussion as follows. The definition of g̃ is modified as follows: If
Ei ∩ B ̸= ∅, for x ∈ Ei , we letg̃(x) = sup f ∗ (x).
x∈Ei ∩B

Otherwise, g(x) = M . With this modification, we obtain
Z
Z
Z
∗
∗
0 ≤ C − 1 = (g̃ − f )dµ =(g̃ − f )dµ +(g̃ − f ∗ )dµBc

B≤ ϵ′ + M ϵ̃,which can be arbitrarily small, so the same result follows.

C

Details on Learning Probability Measures with the Pólya Tree
Process

This section provides details on the weak learner we use to fit tree measures to the residuals in theestimation. The algorithm is based on the PT-based method proposed in Awaya and Ma (2022), andinterested readers may refer to this paper.

C.1

Theoretical Justification of Using the PT-based ModelAs shown in Section Proposition 5, the improvement in the entropy loss is maximized when a fitted
tree is a solution of the problem(n)

arg max
T ∈T

X

(n)

F̃k (A) log

A∈L(T )

47

F̃k (A)
,
µ(A)

(n)

(k−1)

where F̃k is the empirical measure defined by the residuals r (k−1) = {ri}ni=1 . As n → ∞, the(n)empirical measure F̃k (B) converges to F̃k (B) for B ⊂ (0, 1]d , where F̃k is the true distribution of theresiduals defined by the previous tree-CDFs G1 , . . . , Gk−1 . At this population level, the maximization
problem is written asXarg max
T ∈TF̃k (A) log

A∈L(T )F̃k (A),µ(A)and we can show that this maximization is equivalent to minimizing the KL divergence KL(F̃k ||F̃k |T ),where F̃k |T is “a tree-approximation of F̃k under T ”, namely,F̃k |T (B) =XF̃k (A)

A∈L(T )µ(B ∩ A)for B ∈ (0, 1]d .
µ(A)Theorem 4.1 in Awaya and Ma (2022) shows that the posterior of trees also concentrates on theminimizer of KL(F̃k ||F̃k |T ), and this result implies that at the population level, or when n is large,
we can find the tree that maximizes the improvement in the entropy loss or similar ones by checking
the posterior of trees.

C.2

Details on the sampling algorithm

Suppose we have obtained the residuals at the beginning of the boosting algorithm. Since the task offitting a new measure to the residuals is essentially the same for all steps, we drop the k, the index of
the trees and measures consisting of the ensemble, from the notations for simplicity. Then the residuals
are denoted by r = (r1 , . . . , rn ), and our task at each step is to capture their distributional structure
by fitting a dyadic tree. In the section, we provide details on the prior distributions introduced forthe tree T and the stochastic top-down algorithm we use to find a tree with good fitting.
C.2.1

Prior distribution of T

As in Awaya and Ma (2022), we construct a prior of T by introducing the random splitting rule foreach node A. First, we introduce the stopping variable S(A) that takes 0 or 1, and if S(A) = 1, westop splitting A and otherwise split A. Here we set P (S(A) = 1) to 0.5. In the latter case, we next
define the dimension variable D(A) and the location variable L(A). If D(A) = j (j = 1, . . . , d), thenode A is split in the jth dimension, and the location of the boundary is determined by L(A) ∈ (0, 1),
in which 0 (or 1) corresponds to the left (or right) end point. Their prior distributions are as follows:
P (D(A) = j) = 1/d, (j = 1, . . . , d),1
P (L(A) = l/NL ) =(l = 1, . . . , NL − 1),
NL − 1
where NL − 1 is the number of grid points, which is 127 in the estimation.

48

On the tree T , we also define a random measure G̃, with which we can define the likelihood of theresiduals r. The prior of the measure G̃ is defined by introducing the parameters θ(A) = G̃(Al | A),
where Al is the left children node, for every non-terminal node A. They follow the prior distribution
specified asµ(Al )
θ(A) ∼ Beta(θ0 (A), 1 − θ0 (A)), θ0 (A) =.
µ(A)We note that this random measure G̃ is introduced just to define the marginal posterior of T , namely,
P (T | r), since our main goal is to find a tree that fits the distribution of the residuals. In theestimation, we first choose one tree according to this posterior and construct the measure to output,
which is denoted by Gk in the paper, as in Equation 12. The method to select one tree is describedin more detail in the next section.
The joint model of the tree T and the measure G̃ can be seen as a special case of the density
estimation model that is referred to as the adaptive Pólya tree (Ma, 2017) model in Awaya and Ma(2022) with the number of the latent states being 2.
C.2.2

Top-down stochastic algorithm

The particle filter proposed in Awaya and Ma (2022) is shown to be effective to sample from theposterior of trees. This original algorithm, however, has drawbacks when seen as a component of theboosting from a viewpoint of computational cost:
1. In the particle filter, we construct thousands of candidate trees, but this strategy may make thewhole boosting algorithm too time-consuming since in the boosting algorithm we need to repeat
fitting trees to the residuals many times.
2. In the original algorithm, we do not stop splitting nodes until we reach the bottom nodes unlessthe number of included observations is too small. (Technically speaking, this is because thestopping variables, or the latent variables in general, are integrated out in the sampling.) The
number of nodes generated in a tree, however, tends to be large especially when the sample sizeis large, and constructing such large trees repeatedly in the boosting algorithm is also too-time
consuming. The computation cost would become reasonable if we “give up non-promising nodes”,
that is, stop dividing nodes if no interesting structures are found there.
From these reasons, we modify the original algorithm as follows: (i) Instead of generating many
candidate trees, we set the number of particles to one, that is to say, construct a tree by randomly
splitting nodes on the tree in a top-down manner. Hence the algorithm is similar to the top-down
greedy method, but in our algorithm one selects splitting rules stochastically. (ii) For each active
node, we compare possible splitting rules and the decision of stopping the splitting, where the latter
option is added to the algorithm. This comparison is based on their posterior probabilities, and thesplitting tends to be stopped if the conditional distribution is close to uniform.

49For an active node A, the possible decisions are compared based on the following quantities thatare seen as “prior × marginal likelihood”. A conceptually very similar algorithm for supervised learning
is proposed in He and Hahn (2021). For the decision of stopping, we computeL∅ = P (S(A) = 1)µ(A)−n(A) ,where n(A) is the number of residuals included in A. On the other hand, for the splitting rule D(A) = j
and L(A) = l/NL that divides A into Al and Ar , we compute
Lj,l = P (S(A) = 0, D(A) = j, L(A) = l/NL )Z× Beta(θ | θ0 (A), 1 − θ0 (A))θn(Al ) (1 − θ)n(Ar ) dθ× µ(Al )−n(Al ) µ(Ar )−n(Ar )= P (S(A) = 0, D(A) = j, L(A) = l/NL )×

Be(θ0 (A) + n(Al ), 1 − θ0 (A) + n(Ar ))Be(θ0 (A), 1 − θ0 (A))× µ(Al )−n(Al ) µ(Ar )−n(Ar ) ,where Be(·) is the beta function. Based on these quantities, we choose to stop the splitting withprobability
L∅
Pd PNL −1
L∅ + j=1 l=1 Lj,l .
Otherwise, we choose the splitting rule D(A) = j and L(A) = l/NL with probabilityPdj ′ =1D

Lj,lPNL −1l′ =1

Lj ′ ,l′.

Details of the 48-dimensional Experiments

In the experiment, we used the following three scenarios.
Scenario A: n = 1, 000, andNormal(µ, Σ),
where
"
µ=

1/2
1/2

#

"
, Σ=

1/82
0.95/82
0.95/82
1/82#.

Scenario B: n = 5, 000, and1
3
Beta(x1 | 1, 1) × Beta(x2 | 1, 1) + Beta(x1 | 15, 45) × Beta(x2 | 15, 45)10103+ Beta(x1 | 45, 15) × Beta(x2 | 22.5, 37.5)
103+ Beta(x1 | 37.5, 22.5) × Beta(x2 | 45, 15).
1050

Scenario C: n = 2, 000, and1
1
1[0.1,0.45]×[0.35,0.9] (x1 , x2 ) + 1[0.2,0.8]×[0.45,0.5] (x1 , x2 )3
3
1
+ 1[0.7,0.9]×[0.05,0.6] (x1 , x2 ).
3

E

Additional Tables and Figures

Table 4:seeds.

(c0 , γ)

POWER

GAS

HEPMASSMINIBOONE(0.1, 0.0)
(0.1, 0.5)
(c0 , γ)

0.004
0.004
AReM

0.060
0.023
CASP

0.028
0.032
BANK

0.075
0.064

(0.1, 0.0)
(0.1, 0.5)

0.022
0.015

0.038
0.030

0.078
0.023

The standard deviations of the average predictive scores based on 30 different random

POWER

GAS

HEPMASSMINIBOONE

4.5
AReM

9.6
CASP

24.3
BANK

15.3

1.4

1.6

9.0

Table 5:The average computation time (seconds) for simulating 10,000 observations based on 30
different random seeds. The tuning parameters c0 , γ are set to 0.1 and 0.5, respectively.

51

Figure 8: The training set of the POWER data (a subset of size 10,000 is visualized) and 10,000
observations simulated from the learned probability measure.

Figure 9:
The training set of the GAS data (a subset of size 10,000 is visualized) and 10,000
observations simulated from the learned probability measure.

52

Figure 10: The training set of the AReM data (a subset of size 10,000 is visualized) and 10,000
observations simulated from the learned probability measure.

53