# Medoid splits for efficient random forests in metric spaces11footnotemark: 1

Matthieu Bulté22footnotemark: 2<EMAIL_ADDRESS>Helle Sørensen33footnotemark: 3
<EMAIL_ADDRESS>
###### Abstract

This paper revisits an adaptation of the random forest algorithm for Fréchet
regression, addressing the challenge of regression in the context of random
objects in metric spaces. Recognizing the limitations of previous approaches,
we introduce a new splitting rule that circumvents the computationally
expensive operation of Fréchet means by substituting with a medoid-based
approach. We validate this approach by demonstrating its asymptotic
equivalence to Fréchet mean-based procedures and establish the consistency of
the associated regression estimator. The paper provides a sound theoretical
framework and a more efficient computational approach to Fréchet regression,
broadening its application to non-standard data types and complex use cases.

###### keywords:

Least squares regression , Medoid , Metric spaces , Random forest , Random
objects

[label1]organization=Department of Mathematical Sciences, University of
Copenhagen

##  1 Introduction

We study the extension of random forest to regression situations where the
response takes values in a metric space. The usual expectation is replaced by
the Fréchet mean, which can be be very computationally intensive with
restrictions on possible applications to follow. Emphasis is therefore on a
new medoid-based splitting rule used in the individual trees. Application of
this new splitting rule speeds up computations without compromising
consistency or finite-sample prediction quality and makes it possible to
analyzer larger datasets in terms of sample size and dimension.

Random forest is a statistical learning method introduced by Breiman (2001)
for classification and regression. It constructs predictions by averaging the
predictions from a randomized ensemble of regression trees fitted on a
training dataset. Random forests have proven to perform well in practice on
complex and high-dimensional datasets with little parameter tuning, making
them an ubiquitous method in the machine learning toolset.

With the success of statistical analysis, more domains are being analyzed,
leading to the study of more complex use cases and non-standard data types. A
first generalization that has emerged from this need is the study of random
variables taking values in manifolds. More recently, a further generalization
has started to be studied, in which the random variables, referred to as
random objects, take values in a metric space. This general setting does not
require anything more than a distance function, and does not require common
algebraic structures such as those found in vector spaces or on manifolds. A
core quantity in the study of random objects is the Fréchet mean, as defined
by Fréchet (1948), which generalizes the concept of the expected value to
metric spaces by defining it as the minimizer of the expected squared loss.
This notion allows to generalize the definition of the variance and further
construct tools to study the variability of random objects.

The lack of algebraic structure in the study of random objects makes the
development of regression algorithms particularly challenging. Efforts in
adapting regression algorithms to non-Euclidean objects have mainly been done
for the case where the response takes values in finite-dimensional
differentiable Riemannian manifolds. This has led to the development and
analysis of parametric models exploiting the local structure of the space,
such as geodesic regression approaches (Thomas Fletcher, 2013), but also fully
non-parametric approaches via adaptations of the Nadaraya-Watson estimator
(Davis et al., 2010; Hinkle et al., 2012; Pelletier, 2006; Yuan et al., 2012;
Hein, 2009). More recently, Petersen and Müller (2019) proposed a general
framework for regression in metric spaces with Euclidean covariates in which
the regression method specifies a weighting scheme. For a given input, the
training samples are reweighted following that scheme and the prediction is
formed by computing the weighted Fréchet mean. They then define a global
weighting scheme generalizing the classical linear regression and a local
scheme which generalizes local linear regression. This original work paved the
way to the construction and study of more flexible regression techniques
adapted to Fréchet regression.

Random forest is a well-suited candidate for adaptation to Fréchet regression.
Several authors have already used random forests outside the classical setup
of Euclidean response, for example for function spaces (Nerini and Ghattas,
2007; Fu et al., 2021) by either working directly in the Hilbert space $L_{2}$
or via a transformation to $L_{2}$, but these methods lack a theoretical
analysis and their adaptation to metric space data is unclear. Another attempt
was made by Capitaine et al. (2020) by direct translation of the random forest
algorithm via a two-stage approach in which the forest prediction is
constructed as the Fréchet mean of the prediction of the trees. Their approach
has been successfully applied to general settings in which both the covariates
and response are allowed to be in a metric space and proposes a novel
splitting rule that reduces the computational cost of the tree-fitting
procedure. Unfortunately, this two-stage procedure complicates theoretical
analysis, and to our knowledge, no consistency result for their method has
been developed.

Another route for adapting random forest to Fréchet regression is to fit a
random forest and aggregate the trees into weights for the observations, which
can be used to construct a prediction. This approach was developed by
Meinshausen (2006) for quantile regression, and Athey et al. (2019) for
conditional Z-estimation under the name of Generalized Random Forest (GRF). It
has been exploited by Qiu et al. (2022) to construct a random forest adapted
to Fréchet regression. The authors combined this with the splitting rule
proposed by Capitaine et al. (2020) to alleviate some of the computational
burden of tree fitting. The paper presents results about asymptotic properties
of the estimator, but proofs are not provided. Moreover, the results regarding
rates of convergence and a central limit theorem depend on strong assumptions
regarding the random forest algorithm itself.

One essential ingredient in the adaptation of a random forest algorithm is the
construction of a suitable splitting rule used to recursively partition the
covariate space in the fitting of individual trees. In their work, Athey et
al. (2019) propose a general splitting rule for Z-estimation relying on
gradients of the objective function to build splits capturing the
heterogeneity of the response. This approach cannot be used for Fréchet
regression since it requires a differentiability structure which is not
available. A possible approach in Fréchet regression is to use the standard
splitting rule when the response is in $\mathbb{R}$, namely by finding the
splitting position that results in the highest reduction in variance in the
subgroups. However, for each split, this requires computing a large number of
Fréchet variances and hence Fréchet means. In many cases, computing the
Fréchet mean is expensive, making this approach unusable. The approach
proposed by Capitaine et al. (2020) and used in Qiu et al. (2022) based on
2-means clustering for efficient splitting reduces the number of Fréchet means
required to a great extent but, as highlighted by the authors, is still
computationally intensive.

In this article, we revisit the algorithm proposed by Qiu et al. (2022), which
implements GRF for Fréchet regression. Our main contribution is a new
splitting rule which replaces the computation of Fréchet means with medoids,
i.e., sample members. This, after the computation of a distances matrix,
allows to find splits as efficiently as in the Euclidean case. We denote the
adapted procedure Metric Random Forest (MRF). We prove that the new splitting
rule is asymptotically equivalent to using the Fréchet mean-based procedure.
Furthermore, we prove the consistency of the MRF based on classical
M-estimation results from van der Vaart (1998) and following the proof
strategy presented in Wager and Athey (2018), and illustrate the benefit of
the new splitting rule with numerical experiments.

The paper is organized as follows: Section 2 includes background material on
Fréchet regression and existing work on random forest in metric spaces. The
new medoid-based splitting rule is introduced in Section 3, and the
consistency result for the MRF estimator is stated in Section 4. We present
simulations from three scenarios in Section 5 and conclude in Section 6.
Proofs are presented in the appendix.

##  2 Background

###  2.1 Fréchet regression

Let $(\Omega,d)$ be a compact metric space. Consider a pair of random
variables $(X,Y)$, where $X$ is a vector of predictors taking value in
$[0,1]^{d}$ with density bounded away from $0$ and $\infty$, and $Y$ is a
response taking value in $\Omega$. Following the idea of Fréchet (1948), the
notions of mean and variance can be generalized to random objects to obtain
the Fréchet mean and Fréchet variance,

$\omega_{\oplus}=\operatorname*{\arg\\!\min}_{\omega\in\Omega}\mathbb{E}\left[d(Y,\omega)^{2}\right],\qquad
V_{\oplus}=\mathbb{E}\left[d(Y,\omega_{\oplus})^{2}\right].$ (1)

Building upon these generalized concepts of mean and variance, Müller (2016)
define the Fréchet regression function, denoted $m$, in terms of the
conditional distribution of $Y$ given $X=x$,

$m(x)=\operatorname*{\arg\\!\min}_{\omega\in\Omega}M(\omega;x),\qquad
M(\omega;x)=\mathbb{E}\left[d(Y,\omega)^{2}\mid X=x\right],$ (2)

where $M$ is called the conditional Fréchet function.

We are interested in fitting a regression model with response $Y$ and
predictors $X$. Given a dataset
${\mathcal{D}}={\left\\{(X_{i},Y_{i})\right\\}}_{i=1}^{n}$ of independent
random pairs following the same distribution as the prototypical pair $(X,Y)$
and a point $x\in[0,1]^{d}$ for which a prediction is to be obtained, the goal
is to construct an estimator $m_{n}(x)$ of the Fréchet regression function
$m(x)$. One generally applicable solution to this problem is to estimate the
conditional Fréchet function $M(\cdot;x)$ with a data-dependent objective
function $M_{n}(\cdot;x)$ converging to $M(\cdot;x)$ in such a way that the
maximizer $m_{n}(x)$ of $M_{n}(\cdot;x)$ also converges to $m(x)$. A natural
approach consists in using a re-weighting strategy to approximate the
conditional expectation in (2), yielding an estimator of the form

$M_{n}(\omega;x)=\sum_{i=1}^{n}w_{i}(x)d(\omega,Y_{i})^{2},$ (3)

where the weights ${\left\\{w_{i}(x)\right\\}}_{i=1}^{n}$ depend on the data
and the point $x$ for which the prediction is being made. The construction of
such weights has been the subject of a large corpus of literature. A common
and well studied method is the Nadaraya-Watson estimator weights which defines
weights through a kernel function $K$ and a bandwidth parameter
$h\in(0,\infty)$, giving $w_{i}(x)=K_{h}(X_{i}-x)$ with
$K_{h}(\cdot)=h^{-1}K(\cdot/h)$. This method is well understood and easy to
implement, but it requires a crucial tuning of the bandwidth parameter $h$,
suffers from the curse of dimensionality and is not adaptive, in the sense
that the weights are invariant to the part of the covariate space where the
prediction is made. This lack of adaptivity makes the Nadaraya-Watson
estimator practically ill-suited in case of heteroskedasticity of the response
for different values of the predictors.

###  2.2 Random forest in metric spaces

The original and commonly used random forest algorithm for regression proposed
by Breiman (2001) relies on constructing an ensemble of randomized regression
trees. A prediction at a given point is then formed by averaging the
predictions across the ensemble of trees. The fitting of the trees is
randomized by bootstrap aggregation, meaning that each regression tree is
fitted on a random subset of size $s$ of the training data. Each tree is built
greedily from a recursion of splits, each attempting to partition the training
sample available at the node into two subsets of minimal prediction error
(Breiman et al., 1984). The procedure results in a partitioning of the space
of predictors into a collection of rectangles called leaves. The resulting
random forest estimator then takes the form

${\textrm{RF}}(x)=\frac{1}{B}\sum_{b=1}^{B}{\textrm{T}}_{b}(x)\qquad{\textrm{with}}\qquad{\textrm{T}}_{b}(x)=\sum_{i=1}^{n}\frac{\mathbbm{1}{\left\\{X_{i}\in
L_{b}(x)\right\\}}}{\\#{\left\\{j:X_{j}\in L_{b}(x)\right\\}}}Y_{i},$ (4)

where $B$ is the size of the ensemble, ${\textrm{T}}_{b}$ is the $b$-th
regression tree and $L_{b}(x)$ is the subset of training points used for
fitting ${\textrm{T}}_{b}$ falling in the same leaf as $x$.

In order to generalize this idea to the context of Fréchet regression,
Capitaine et al. (2020) directly translate the average-of-averages
construction of the random forest by replacing averages with Fréchet means.
The resulting estimator, called the Fréchet Random Forest (FRF), is
constructed as the Fréchet mean of an ensemble of randomized Fréchet trees
(FTs) in which each tree is built to minimize the prediction error, as done in
the Euclidean regression setup.

In this manuscript, a different approach, first proposed in Qiu et al. (2022),
is followed to constructing a random forest prediction. It does not aggregate
the prediction of each individual tree; instead the results of the fitted
ensemble are aggregated following the idea presented in Meinshausen (2006) and
Athey et al. (2019). For a given $x$, each tree contributes with a weighting
of the observations reflecting the partitioning of the predictor space fitted
by the tree. Consider a source of auxiliary randomness $\xi$ and a subset of
the data $D\subset{\mathcal{D}}$, we denote by
${\left\\{w_{i}(\cdot;\xi,D)\right\\}}_{i=1}^{n}$ the set of weights resulting
from fitting a single tree on $D$. These weights take the form

$w_{i}(x;\xi,D)=\frac{\mathbbm{1}{\left\\{X_{i}\in
L(x;\xi,D)\right\\}}}{N(L(x;\xi,D))},$ (5)

where $N(C)=\\#{\left\\{j:X_{j}\in C\right\\}}$ for any $C\subset[0,1]^{d}$,
and $L(x;\xi,D)$ represents the observations in $D$ falling in the same leaf
as $x$. Note that for observations $(X_{i},Y_{i})\in{\mathcal{D}}\backslash
D$, the associated weight function $w_{i}(\cdot;\xi,D)$ is zero. Given $B$
random subsets ${\left\\{D_{b}\right\\}}_{b=1}^{B}$ of ${\mathcal{D}}$ and
auxiliary randomness ${\left\\{\xi_{b}\right\\}}_{b=1}^{B}$, the weights
constructed by the individual trees are averaged within the ensemble to
construct an adaptive set of weights ${\left\\{w_{i}\right\\}}_{i=1}^{n}$ for
a given prediction point $x$ given by

$w_{i}(x)=\frac{1}{B}\sum_{b=1}^{B}w_{i}(x;\xi_{b},D_{b}).$ (6)

These weights are then used to form an estimator of the conditional Fréchet
function defined in (2),

$M_{n}(\omega;x)=\sum_{i=1}^{n}w_{i}(x)d(\omega,Y_{i})^{2}.$ (7)

Minimizing this adaptively weighted sum yields our Metric Random Forest (MRF)
estimator of the Fréchet regression function in (2),

$m_{n}(x)=\operatorname*{\arg\\!\min}_{\omega\in\Omega}M_{n}(\omega;x).$ (8)

It is important to note that, in general, this optimization problem needs not
have a unique solution. In practice, a user of the MRF must provide a solver
for (8) that can return a unique solution based on any preferred heuristic
implemented by the user (e.g., minimal norm, first in lexicographic order).
Since this technicality does not influence the implementation or analysis of
the MRF algorithm, we can safely assume that $m_{n}(x)$ exists and is unique.
Note that we later assume in Assumption 1 that the population version of this
quantity exists and is unique.

##  3 Splitting criterion in metric spaces

The recursive partitioning used to construct each individual tree is driven by
the choice of a splitting criterion. Considering an arbitrary cell
$C\subset[0,1]^{d}$, the goal at each partitioning step is to find the feature
index $j\in{\left\\{1,\ldots,d\right\\}}$ and cut position $z\in[0,1]$
defining an axis-aligned plane along which to split the cell into two
subcells, $C_{l}$ and $C_{r}$, given by

$C_{l}={\left\\{x\in A:x^{(j)}\leq z\right\\}}\qquad{\textrm{and}}\qquad
C_{r}={\left\\{x\in A:x^{(j)}>z\right\\}},$

where $x^{(j)}$ is the $j$th coordinate of $x$. Intuitively, each splitting
pair $(j,z)$ should be chosen as to improve the estimation of the Fréchet
regression function $m$. This is done by minimizing the CART splitting
criterion of Breiman et al. (1984). Let us denote by
$\hat{\omega}_{\oplus}(C)$ and $\hat{V}_{\oplus}(C)$ the empirical Fréchet
mean and variance computed on the observations with $X$ contained in $C$,

$\displaystyle\hat{\omega}_{\oplus}(C)$
$\displaystyle=\operatorname*{\arg\\!\min}_{\omega\in\Omega}\sum_{i=1}^{n}\mathbbm{1}{\left\\{X_{i}\in
C\right\\}}d(Y_{i},\omega)^{2},$ $\displaystyle\hat{V}_{\oplus}(C)$
$\displaystyle=\sum_{i=1}^{n}\frac{\mathbbm{1}{\left\\{X_{i}\in
C\right\\}}}{N(C)}d(Y_{i},\hat{\omega}_{\oplus}(C))^{2}.$

The CART splitting criterion is then given by

${\textrm{err}}(j,z)=\frac{N(C_{l})}{N(C)}\hat{V}_{\oplus}(C_{l})+\frac{N(C_{r})}{N(C)}\hat{V}_{\oplus}(C_{r}).$
(9)

Note that minimizing this criterion corresponds to finding a split that
maximizes the decrease in the weighted average of Fréchet variance accross the
splitted subcells. A splitting pair $(j^{\star},z^{\star})$ can be found by
iterating over the feature index $j\in\\{1,\ldots,d\\}$ and possible threshold
values $z\in\\{X^{(j)}_{i}:X_{i}\in C\\}$ determined from values of coordinate
$j$ taken by the samples present in the cell $C$.

As highlighted by Athey et al. (2019) and Capitaine et al. (2020), this greedy
search algorithm may be computationally prohibitively expensive since it
requires for each split and candidate threshold to compute Fréchet means
$\hat{m}_{l}$ and $\hat{m}_{r}$, which overall requires $O(dn^{2})$
computations of Fréchet means. In many cases, computing a Fréchet mean can
only be done approximately through computationally expensive algorithms. This
is the case for computation of Fréchet means in Wasserstein spaces (Panaretos
and Zemel, 2020), in function spaces equipped with the amplitude and phase
variation distances (Srivastava and Klassen, 2016) or in Riemannian manifolds.

To overcome this computational limitation, Athey et al. (2019) introduce the
gradient tree algorithm. This algorithm utilizes an alternative splitting
criterion based on a linearization of the target optimization problem in
regression and enables faster calculation of the relevant quantities in each
candidate subcell. However, this approach cannot be applied to general metric
spaces where no vector space structure is available for differentiability. In
Capitaine et al. (2020), an alternative splitting criterion is proposed which
offers computational advantages over the CART splitting criterion and is also
usable in metric spaces. This approach, which is also used in Qiu et al.
(2022), reduces the number of Fréchet mean computations by testing a single
possible partitioning for each coordinate, where the partitioning is found by
performing a $2$-means clustering of the observations in the cell based on the
tested coordinate. This allows to compute only $2d$ Fréchet means per split.
However, this still results in $O(dn)$ computations of Fréchet means which
renders the algorithm unusable in situations where the Fréchet mean is
expensive to compute. Notice that trees are fitted to subsamples of
$\mathcal{D}$ of size $s$, but the abovementioned computational complexities
are still valid with $n$ replaced by the size of the subsamples.

We propose a solution that is usable in any metric space and completely avoids
the need to compute Fréchet means during the fitting of the tree. To do that,
we replace the minimization of the Fréchet function over the entire metric
space by a minimization over the available sample. Denoting by $\Omega_{n}$
the set of observed responses,
$\Omega_{n}={\left\\{Y_{1},\ldots,Y_{n}\right\\}}$, we define an approximation
to the empirical Fréchet mean $\hat{\omega}_{\oplus}(C)$, the Fréchet medoid
estimator $\tilde{\omega}_{\oplus}(C)\in\Omega_{n}$, by

$\tilde{\omega}_{\oplus}(C)=\operatorname*{\arg\\!\min}_{\omega\in\Omega_{n}}\sum_{i=1}^{n}d(\omega,Y_{i})^{2}\mathbbm{1}{\left\\{X_{i}\in
C\right\\}}.$

Computation of distances between all pairs of elements in the training data is
requires, but once these distances have been computed, the Fréchet medoid is
straightforward to find. Using the Fréchet medoid hence provides an efficient
alternative to the empirical Fréchet mean in the computation of the CART
splitting criterion (9).

The algorithm we propose to fit the random forest is based on replacing the
computation of the empirical Fréchet mean by the Fréchet medoid in the split-
finding procedure. As a first step, we pre-compute the pairwise distances
matrix $\Delta\in\mathbb{R}^{n\times n}$ with $\Delta_{ij}=d(Y_{i},Y_{j})$.
Then, at each step, we find the split $(j,z)$ minimizing the approximate
splitting criterion

$\tilde{{\textrm{err}}}(j,z)=\frac{1}{N_{n}(C)}\sum_{i=1}^{n}\Delta_{i,i^{\star}(C_{l})}^{2}\mathbbm{1}{\left\\{X_{i}\in
C_{l}\right\\}}+\Delta_{i,i^{\star}(C_{r})}^{2}\mathbbm{1}{\left\\{X_{i}\in
C_{r}\right\\}},$

where $i^{\star}$ is the index of the Fréchet medoid,
$\tilde{\omega}_{\oplus}(C)=Y_{i^{\star}(C)}$.

Provided that the discrete set $\Omega_{n}$ grows dense in $\Omega$ and
further technical conditions, the Fréchet medoid $\tilde{\omega}_{\oplus}(C)$
consistently estimates the true Fréchet mean (Proposition 1), and furthermore
the associated splitting criterion $\tilde{{\textrm{err}}}(j,z)$ inherits this
consistency. Notice that Assumption 1 is stated in Section 4, and that a proof
for the proposition is given in the appendix.

######  Proposition 1.

Let $Y_{1},\ldots,Y_{n}$ be independent copies of a prototypical random
variable $Y$ taking values in a compact metric space $(\Omega,d)$ with Fréchet
mean $\omega_{\oplus}$. Assume that the Fréchet mean is well-separated as
defined in Assumption 1, and that the Fréchet mean is a possible value for
$Y$, in the sense that for every $\varepsilon>0$,

$\mathbb{P}\left[d(Y,\omega_{\oplus})<\varepsilon\right]>0.$

Then, the Fréchet medoid estimator is a consistent estimator of the Fréchet
mean, that is, $d(\tilde{\omega}_{\oplus},\omega_{\oplus})=o_{P}(1)$.

Together with the Lipschitz continuity of the squared distance function (Lemma
1 in the appendix), one can see that for some constant $K>0$ independent of
$n$,

$\left|\tilde{{\textrm{err}}}(j,z)-{\textrm{err}}(j,z)\right|\leq
K{\left\\{d(\hat{\omega}_{\oplus}(C_{l}),\tilde{\omega}_{\oplus}(C_{l}))+d(\hat{\omega}_{\oplus}(C_{r}),\tilde{\omega}_{\oplus}(C_{r}))\right\\}}.$

Since $\tilde{\omega}_{\oplus}$ and $\hat{\omega}_{\oplus}$ are both
consistent estimators of $\omega_{\oplus}$, this implies that the proposed
splitting criterion converges to the CART splitting criterion.

##  4 Asymptotic Theory

We now consider the pointwise consistency of the MRF. We show that for a fixed
vector of covariates $x\in[0,1]^{d}$, the MRF (8) converges in probability to
the output of the regression function (2). Consistency of the random forest
estimators have been established in various settings: for simplified forest
models Breiman (2004); Biau et al. (2008); Scornet (2016), forest with real
number responses (Wager and Athey, 2018; Scornet, 2016) and also in the more
general Z-estimator setting (Athey et al., 2019). However, due to the reliance
on the outcome taking values in an Euclidean space, both in the assumptions
and in the proofs of their theoretical results, existing asymptotic analysis,
in particular that of the GRF, cannot be directly applied to the MRF setting.

Instead, we develop a proof combining results from standard M-estimation
theory and existing results of pointwise consistency of the classical
Euclidean random forest algorithm. Based on standard results in M-estimation
theory (see van der Vaart (1998)), the first step to establishing pointwise
consistency of the MRF is showing the convergence of the approximate objective
function $M_{n}$ given in (7) to the conditional Fréchet function $M$ from (2)
uniformly in $\omega\in\Omega$. We show this based on results from Newey
(1991) and the pointwise convergence of the MRF objective function based on
the proof presented in Wager and Athey (2018). Their proof combines the
Lipschitz continuity of the objective function with an assumption on the
construction of the trees which implies that the diameter of the leaves of the
tree shrinks as more datapoints become available (Specification 1 from Athey
et al. (2019), see below).

Similarly to Wager and Athey (2018), we study in our analysis a theoretical
construction of the random forest in which the bagging procedure is replaced
by fitting trees on each possible subsets $D$ of size $s$ of the dataset
${\mathcal{D}}$. The effect of using an approximation of (11) based on a
finite number of subsamples of ${\mathcal{D}}$ is ignored here, but a detailed
investigation by Mentch and Hooker (2016) and Wager et al. (2014) suggest
using $B=O(n)$.

The set of weights ${\left\\{w_{i}(\cdot;\xi,D)\right\\}}_{i=1}^{n}$ given in
(5) can be used to define the contribution of a single tree to the forest
objective function

$T(\omega,x;\xi,D)=\sum_{i=1}^{n}w_{i}(x;\xi,D)d(\omega,Y_{i})^{2}.$ (10)

The tree-level weights are then combined to form the theoretical random forest
weights

$w_{i}^{\star}(x)=\binom{n}{s}^{-1}\sum_{\begin{subarray}{c}D\subset{\mathcal{D}}\\\
\left|D\right|=s\end{subarray}}\mathbb{E}_{\xi}\left[w_{i}(x;\xi,D)\right].$
(11)

Based on these weights, we can define the theoretical objective function and
regression function studied in this section

$M_{n}^{\star}(\omega;x)=\sum_{i=1}^{n}w_{i}^{\star}(x)d(\omega,Y_{i})^{2}\qquad{\textrm{and}}\qquad
m_{n}^{\star}(x)=\operatorname*{\arg\\!\min}_{\omega\in\Omega}M_{n}^{\star}(\omega;x).$
(12)

We further consider that the forest is implemented in a way that satisfies
Specification 1 of Athey et al. (2019): All trees in the forest are symmetric,
meaning that the order of the input is not relevant to the fitting of each
tree. The splits are balanced in that each split separates the observations on
two subsets, each with a proportion of at least $\alpha>0$ of the parent.
Additionally, each tree is grown to depth $k$, for some $k\in\mathbb{N}$ and
each leaf contains between $k$ and $2k-1$ observations. The forest is honest
as described in Wager and Athey (2018), and each tree is built on a subsample
of the data of size $s<n$ satisfying $s/n\rightarrow 0$ and
$s\rightarrow\infty$.

The following assumption is common in M-estimation, see for instance
Assumption (P0) in Müller (2016), and is used in most proofs of consistency of
Fréchet regression methods.

######  Assumption 1 (Existence of a population solution).

For every $x\in[0,1]^{d}$, the conditional Fréchet regression function $m(x)$
given in (2) exists and is a well-separated solution to the conditional
Fréchet function $M(\cdot;x)$, that is for every $\varepsilon>0$,

$\inf_{\omega:d(\omega,m(x))>\varepsilon}M(\omega;x)>M(m(x);x).$

We further require the conditional Fréchet function to be Lipschitz continuous
in $x$. This assumption is also common in Fréchet regression, and is used here
to exploit the shrinkage of the tree leaves to show the convergence of the
approximate Fréchet function defined by the random forest.

######  Assumption 2 ($x$-Lipschitz Fréchet function).

For any fixed $\omega\in\Omega$, the function $x\mapsto M(\omega;x)$ is
Lipschitz continuous with Lipschitz contant $C(\omega)$ possibly depending on
$\omega$.

As indicated, we first state a theorem saying that the objective function
$M_{n}^{\star}(\cdot;x)$ defined by the random forest algorithm converges
uniformly to the true Fréchet conditional Fréchet function $M$. This is a
necessary element in consistency proofs of M-estimators and thus used to prove
consistency of the random forest algorithm (Theorem 2). Proofs are given in
the appendix.

######  Theorem 1.

Under Assumption 2, the random forest objective function
$M_{n}^{\star}(\cdot;x)$ converges uniformly in $\omega$ to $M(\cdot;x)$, that
is, for all $x\in[0,1]^{d}$

$\sup_{\omega\in\Omega}\left|M_{n}^{\star}(\omega;x)-M(\omega;x)\right|=o_{P}(1).$

######  Theorem 2.

Under Assumptions 1 and 2, the random forest regressor defined in (12) is a
pointwise consistent estimator of the conditional Fréchet regression function,
that is, for each $x\in[0,1]^{d}$,

$d(m_{n}^{\star}(x),m(x))=o_{P}(1).$

##  5 Numerical experiments

In the following, we perform a simulation study in order to evaluate the
benefit of our medoid approach (MRF) as compared to the random forest weighted
local constant Fréchet regression (RFWLCFR) algorithm presented in Qiu et al.
(2022). The main difference between the two algorithms is the choice of
splitting criterion, cf. Section 3, and we will compare them in terms of
computational efficiency and prediction quality.

We study three scenarios with similar data generation processes in three
different metric spaces with different computational costs associated with the
computation of empirical Fréchet means. The first example is that of one-
dimensional density functions over the real line $\mathbb{R}$ in the
Wasserstein space. This space is isomorphic to a convex subspace of $L_{2}$ in
which computations of distances and Fréchet means are straightforward. In a
second example, we study data lying on a Riemannian manifold, the sphere
$S^{2}\subset\mathbb{R}^{3}$. In this example distances are straightforward to
compute but Fréchet means can only be computed through a more expensive
gradient descent algorithm. In a last example, we study functional data
presenting phase variability through the study of their warping functions as
described in Srivastava and Klassen (2016), in which both the evaluation of
distances and of the Fréchet mean are computationally demanding.

In each of the three metric spaces, we generate datasets of different training
sample sizes $N\in{\left\\{100,200,400\right\\}}$ and covariate dimension
$d\in{\left\\{2,5,10,20\right\\}}$. For each combination of a metric space, a
training sample size and a covariate dimension, we generate $100$ random
datasets with the desired properties. In order to generate a dataset with
covariate dimension $d$, two parameters $\alpha\in\mathbb{R}$ and
$\beta\in\mathbb{R}^{d}$ are sampled with independent standard normal
components. The Fréchet regression function is then generated from a single
index regression model

$m(x)=g(\eta)\qquad{\textrm{with}}\qquad\eta=\alpha+(x-0.5)^{\top}\frac{\beta}{\sqrt{d}}$
(13)

for some function $g:\mathbb{R}\rightarrow\Omega$, which is fixed and specific
to each metric space. The scaling $1/\sqrt{d}$ has been chosen to ensure that
the distribution of $\eta$ does not change too much as the number of
parameters increases, thus keeping the different simulation settings
comparable. Each observation $(X_{i},Y_{i})$ is finally constructed by first
sampling $X_{i}\sim{\textrm{Unif}}[0,1]^{d}$, computing the conditional
expectation $m(X_{i})$, and applying an independently sampled noise function
$T_{i}:\Omega\rightarrow\Omega$, $Y_{i}=T_{i}(m(Y_{i}))$. Note that $m$ varies
is fixed within each dataset but varies between datasets.

In a given experiment, we fit the MFR and RFWLCFR algorithms to the same
training dataset and evaluate them on an independently generated test dataset
of size $N_{\textrm{test}}=100$ with the same parameters $\alpha,\beta$ as in
the training set. We measure the training time of each algorithm and evaluate
the performance of the algorithm using the mean squared error (MSE). For a
test set ${\left\\{(X_{i},Y_{i})\right\\}}_{i=1}^{N_{\textrm{test}}}$ and a
fitted regressor $\hat{m}$, the MSE of the algorithm is given by

${\textrm{MSE}}(\hat{m})=\frac{1}{N_{\textrm{test}}}\sum_{i=1}^{N_{\textrm{test}}}d(\hat{m}(X_{i}),m(X_{i}))^{2}.$

The focus of our comparison is the runtime of the algorithms, while mean
squared errors are used to assess whether the approximation in MRF has a
negative impact on the quality of the predictions. The simulations show the
MFR provides significant performance gains over RFWLCFR, yet without
compromising the quality of the predictions, and therefore enables us to
analyze larger datasets in terms of sample size and dimension.

All simulations and analysis are done in Python. The implementation of both
methods, together with the code to reproduce the results is available at
https://github.com/matthieubulte/MetricRandomForest.

###  5.1 Univariate distributions

We consider the 2-Wasserstein metric space $\Omega=W_{2}=W_{2}(\mathbb{R})$,
the space of probability distributions over the real line with finite second
moment endowed with the 2-Wasserstein metric $d_{W_{2}}$. Specifically, for
$\mathbb{P},\mathbb{Q}\in W_{2}$ with cumulative distribution functions
$F_{\mathbb{P}}$ and $F_{\mathbb{Q}}$, the 2-Wasserstein distance between $P$
and $Q$ is given by

$d_{W_{2}}(\mathbb{P},\mathbb{Q})=d_{L_{2}[0,1]}(F_{\mathbb{P}}^{-1},F_{\mathbb{Q}}^{-1})=\sqrt{\int_{0}^{1}\left|F_{\mathbb{P}}^{-1}(u)-F_{\mathbb{Q}}^{-1}(u)\right|^{2}\text{d}u}.$

abc Figure 1: (a) shows the mean function $m(X)$ for 10 sampled values of $X$.
(b) and (c) compare the runtime and MSE, respectively, of the MRF and RFWLCFR
algorithms in Wasserstein regression for $N\in{\left\\{100,200,400\right\\}}$
and $d\in{\left\\{2,5,10,20\right\\}}$.

As shown in existing literature on Wasserstein regression (Petersen and
Müller, 2019; Ghodrati and Panaretos, 2022), the solution to the weighted
Fréchet problem given via the sample
${\left\\{\mathbb{P}_{i}\right\\}}_{i=1}^{n}\subset W_{2}$ and weights
${\left\\{w_{i}\right\\}}_{i=1}^{n}$ can be found in a two-steps procedure: a
pointwise weighted average of the sample inverse cumulative distribution
functions is first computed and then projected onto the space of non-
decreasing functions. This results in the following optimization problem

$\operatorname*{\arg\\!\min}_{\mathbb{Q}\in
W_{2}}\sum_{i=1}^{n}w_{i}d_{W_{2}}(\mathbb{Q},\mathbb{P}_{i})^{2}=\operatorname*{\arg\\!\min}_{\mathbb{Q}\in
W_{2}}\left\lVert
F_{\mathbb{Q}}^{-1}(\cdot)-\sum_{i=1}^{n}w_{i}F_{\mathbb{P}_{i}}^{-1}(\cdot)\right\rVert_{L_{2}[0,1]}^{2}$

In practice, each distribution is represented by a vector
$y_{i}=(y_{i1},\ldots,y_{iM})$ of evaluations of its inverse cumulative
distribution function on an equispaced grid
${\left\\{u_{m}\right\\}}_{m=1}^{M}$ on $[0,1]$. The weighted Fréchet mean
$\bar{y}$ is then given by

$F^{-1}_{\bar{y}}=\operatorname*{\arg\\!\min}_{\begin{subarray}{c}v\in\mathbb{R}^{M}\\\
v_{1}\leq\cdots\leq
v_{M}\end{subarray}}\sum_{i=1}^{n}w_{i}\sum_{m=1}^{M}\left|v_{m}-y_{im}\right|^{2}.$

In our application, we chose $M=100$. As highlighted in Ghodrati and Panaretos
(2022), this is an isotonic regression problem which can be efficiently solved
using scikit-learn’s implementation based on Pool Adjacent Violators Algorithm
(PAVA), see Pedregosa et al. (2011); Best and Chakravarti (1990).

Following the example in Petersen and Müller (2019) for data generation, the
conditional mean in this experiment maps a covariate vector to a normal
distribution with variable mean and variance parameters,

$\displaystyle F^{-1}_{m(x)}(u)=\mu(\eta)+\sigma(\eta)\Phi^{-1}(u).$

Specifically, we use the following functions for the mean and variance
components,

$\mu(\eta)=\eta\qquad{\textrm{and}}\qquad\sigma(\eta)=\sigma_{0}+\gamma{\textrm{logit}}^{-1}(\eta),$

with $\sigma_{0}=1$ and $\gamma=2.5$. The random responses are generated by
adding random noise to the normal quantiles as described in Panaretos and
Zemel (2016). Given a continuous and non-decreasing random noise function
$\varepsilon:\mathbb{R}\rightarrow\mathbb{R}$, the noise map
$T:\Omega\rightarrow\Omega$ is defined via the composition of $\varepsilon$
with the quantile function of $\mathbb{P}$, giving
$F^{-1}_{T(\mathbb{P})}=\varepsilon\circ F^{-1}_{\mathbb{P}}$. The noise
functions are given by $\varepsilon(x)=x-\sin\left(\pi kx\right)/\left|\pi
k\right|$, where $k$ is an integer frequency uniformly sampled from
${\left\\{-4,\ldots,4\right\\}}\backslash{\left\\{0\right\\}}$. The
observation $Y$ given $X$ are thus given by $F^{-1}_{Y}=\varepsilon\circ
F^{-1}_{m(X)}$. Note that adding randomness via composition with a non-linear
function results in random distributions outside the class of normal
distributions.

Figure 1 compares the time performance and MSE of the methods. Due to the
isomorphism of $W_{2}$ to the subspace of $L_{2}[0,1]$ of quantile functions,
the computation of Fréchet means is rather straightforward. However, the
$2$-means procedure adds an overhead to the RFWLCFR methodology which makes
the overall computation slower than our MRF. The methods have similar MSE
distributions in all setups showing that the medoid approximation in the MRF
does not affect the quality of the predictions.

###  5.2 Spherical data

abc Figure 2: (a) shows the mean curve on the sphere (solid line) as a
function of $\nu$ in (14) together with random observations (crosses). (b) and
(c) compare the runtime and MSE, respectively, of the MRF and RFWLCFR
algorithms for regression on the sphere $\mathbb{S}^{2}$ for
$N\in{\left\\{100,200,400\right\\}}$ and $d\in{\left\\{2,5,10,20\right\\}}$.

We consider the example where $\Omega=\mathbb{S}^{q}\subset\mathbb{R}^{q+1}$
is the $q$-dimensional sphere equipped with the geodesic distance

$d_{\mathbb{S}^{q}}(x,y)=\arccos\langle x,y\rangle.$

In this case, computing a weighted Fréchet mean is less straightforward and
one must resort to optimization method on manifolds. We do this using the
trust region algorithm implemented in Pymanopt (Townsend et al., 2016).

We focus on data lying on the sphere $\mathbb{S}^{2}\subset\mathbb{R}^{3}$
generated as follows. First, the mean function is given by the single index
model described earlier, where we first transform $\eta$ to $\nu\in(0,1)$ via
$\nu={\textrm{logit}}^{-1}(\eta)$, and define

$g(\eta)=\left(\sqrt{1-\nu^{2}}\cos(\pi\nu),\sqrt{1-\nu^{2}}\sin(\pi\nu),\nu\right).$
(14)

Given a random $X\sim{\textrm{Unif}}[0,1]^{d}$ and its associated mean
function $m(X)$, the response $Y\in\mathbb{S}^{2}$ is generated by
transforming a bivariate random vector $U$ from the tangent space
$T_{m(X)}\mathbb{S}^{2}$ at $m(X)$ through the corresponding exponential map,

$Y={\textrm{Exp}}_{m(X)}(U)=\cos(\left\lVert
U\right\rVert)m(X)+\sin(\left\lVert U\right\rVert)\frac{U}{\left\lVert
U\right\rVert}.$

The random noise $U$ is sampled from a bivariate normal distribution with
independent components and variance $\sigma^{2}=0.1$, that is
$U\sim\mathcal{N}(0,\sigma^{2}\mathbbm{1}_{2\times 2})$.

The results of the experiment are displayed in Figure 2. The time performance
profile is similar to the one in the Wasserstein experiment with MRF still
providing a clear advantage over the RFWLCFR in terms of runtime. In this
example the higher cost of computing Fréchet means does not seem to translate
into further runtime improvements of MFR over RFWLCFR. The two methods again
yield estimators with similar error distribution.

###  5.3 Space of warping functions

We now consider the study of functional data displaying phase variability,
with observations lying in a non-linear subset of $L_{2}[0,1]$. As described
in Tucker et al. (2013), and more generally Srivastava and Klassen (2016), the
variability of functional data can be decomposed in phase and amplitude
components. Intuitively, phase variation corresponds to variation on the $x$
axis and amplitude variation corresponds to variation along the $y$ axis.

Let $\Gamma$ be the set of boundary preserving diffeomorphisms on $[0,1]$,
that is,
$\Gamma={\left\\{\gamma:[0,1]\rightarrow[0,1]\mid\gamma(0)=0,\gamma(1)=1,\gamma{\textrm{
is a diffeomorphism}}\right\\}}$. For a given function $y_{0}\in L^{2}[0,1]$,
we define the orbit of $y_{0}$ as the set of functions that only differ from
$y_{0}$ in their phase component,
$[y_{0}]={\left\\{y_{0}\circ\gamma\mid\gamma\in\Gamma\right\\}}$. While
$y_{0}$ is not unique in this notation, whenever the starting point $y_{0}$ of
an orbit $[y_{0}]$ is fixed in a context, we call it template function.

We consider a distribution for $Y$ in which the observations are all part of
the same orbit $\Omega=[y_{0}]$. Such spaces of equivalence classes, which are
all equivalent to the space of warping functions $\Gamma$, form a non-linear
infinite-dimensional manifold which is difficult to analyze and work with when
considering the $L_{2}$ distance.

Instead, a warping function can be represented by its square-root velocity
function $\psi_{\gamma}=\sqrt{\dot{\gamma}}$. Since the warpings $\gamma$ are
boundary preserving diffeomorphisms on $[0,1]$ we have that $\dot{\gamma}$
exists and is positive, making $\psi_{\gamma}$ a well defined quantity.
Furthermore $\gamma(0)=0$ and $\gamma(1)=1$, hence
$\left\lVert\psi\right\rVert_{L_{2}}^{2}=\int_{0}^{1}\psi(t)^{2}\text{d}t=\int_{0}^{1}\dot{\gamma}(t)\text{d}t=1$
making $\psi$ an element of the $L_{2}[0,1]$ sphere,
$\mathbb{S}^{\infty}={\left\\{f\in L_{2}[0,1]\mid\left\lVert
f\right\rVert_{L_{2}[0,1]}=1\right\\}}$, and more specifically on the positive
hemisphere of $\mathbb{S}^{\infty}$ since $\psi>0$. This allows to define a
distance on $\Gamma$ via the intrinsic distance on $\mathbb{S}^{\infty}$. Let
$\gamma_{1},\gamma_{2}$ be two warpings with square-root velocity functions
$\psi_{1},\psi_{2}$, then the distance between $\gamma_{1}$ and $\gamma_{2}$
is given by

$d_{\Gamma}(\gamma_{1},\gamma_{2})=d_{\mathbb{S}^{\infty}}(\psi_{1},\psi_{2})=\arccos\langle\psi_{1},\psi_{2}\rangle_{L_{2}}.$

See Chapter 4 of Srivastava and Klassen (2016) for more details on this metric
space. Let now $y_{1},y_{2}\in[y_{0}]$ be two functions within a same orbit,
then these functions can be represented by their warpings from the template
$y_{0}$, denoted $\gamma_{1},\gamma_{2}$, with $y_{i}=y_{0}\circ\gamma_{i}$
and we can define

$d_{[y_{0}]}(y_{1},y_{2})=d_{\Gamma}(\gamma_{1},\gamma_{2}).$

The computation of quantities required to evaluate distances and the Fréchet
mean are based on implementations provided in scikit-fda (Ramos-Carreño et
al., 2019) and fdasrsf (Tucker et al., 2013).

abc Figure 3: (a) shows the template curve $y_{0}$ (solid) and a randomly
sampled curve (dashed). (b) and (c) compare the runtime and MSE, respectively,
of the MRF and RFWLCFR algorithms for regression in $\Gamma$ for
$N\in{\left\\{100,200,400\right\\}}$ and $d\in{\left\\{2,5,10,20\right\\}}$.

Next, we turn to the data generating process for our experiment. The template
function is chosen as $y_{0}(u)=(1-(u-0.5)^{2})\sin(9\pi u)$, see the solid
curve in panel (a) of Figure 3. For a given random vector of covariates $X$,
one obtains the linear component $\eta$ as defined in (13) and the mean
warping response is $\gamma$ is given for all $u\in[0,1]$ by

$\gamma(u)=(\exp(4au)-1)/(\exp(4a)-1),$

where $a=3({\textrm{logit}}^{-1}(\eta)-0.5)$. The warping is perturbed by
adding noise to its square-root velocity function $\psi\in\mathbb{S}^{\infty}$
by sampling a random Gaussian process $V$ with exponential covariance kernel
from the tangent space $T_{\psi}\mathbb{S}^{\infty}$, and applying the
exponential map at $\psi$,

${\textrm{Exp}}_{\psi}V=\cos(\left\lVert V\right\rVert)\psi+\sin(\left\lVert
V\right\rVert)\frac{V}{\left\lVert V\right\rVert}.$

The response $Y$ is then obtained by warping the template $y_{0}$ with the
warping function obtained by taking the inverse square-root velocity transform
of ${\textrm{Exp}}_{\psi}V$. An example of a sample is displayed in panel (a)
of Figure 3 (dashed curve).

The results of the experiment are displayed in Figure 3. Computing a single
empirical Fréchet mean can take up to several seconds and computing Fréchet
variance is thus a very expensive operation. The benefit of the 2-means
approach implemented in RFWLCFR is dominated by the high cost of computing
Fréchet variances. On the other hand, the runtime of MFR remains low in all
scenarios. Runtime ratios are 1 to 50 in the $N=100$ and $d=20$ scenario and 1
to 30 for $N=200$ and $d=20$. Similarly to the previous two experiments, the
runtime gain does not impact the fitting of the forest and the MSE
distributions stay similar for the two methods.

##  6 Discussion

Defining non-parametric regression method for general metric spaces that are
performant and work well in complex setups is a challenge. We have proposed a
new version of the random forest Fréchet regression that allows practitioners
to perform a random forest regression in high-dimensional problems and in
metric spaces where computation of Fréchet means is expensive. We achieved
this by replacing the classic computation of the Fréchet means in random
forest splits with a computation based on the Fréchet medoid. We proved both
the consistency of the random forest Fréchet regression as well as the
consistency of our medoid-based approximation. Furthermore, a set of numerical
experiments in metric spaces of different characteristics demonstrated
important reduction in the time required to fit a random forest without
showing any impact on the mean squared error of the method.

One possible extension of this work could be to study rates of convergence, at
best based on assumptions on the data generating process only. This, combined
with a further investigation of the medoid approximation, could refine the
understanding of the asymptotic behaviour of the metric random forest. Our
work also paves the way to exploring further approximations to random forest
regression in metric spaces. This could be through the use of algorithmic
improvements generally applicable to random forests, such as the improvement
proposed in Tiwari et al. (2022), or new approximations specific to the use of
metric spaces. One could for instance study the use of approximate distance
computations or other approximations to the Fréchet variance, and incorporate
such approximations in the theoretical analysis as well as in the
implementation.

##  7 Acknowledgement

This work has received funding from the European Union’s Horizon 2020 research
and innovation program under the Marie Skłodowska-Curie grant agreement No
956107, ”Economic Policy in Complex Environments (EPOC)”.

##  Appendix A Proofs

###  A.1 Statement of useful theorems

We start by stating two theorems that will be used in the proof of
consistency. First, Corollary 2.2 of Newey (1991) provides the conditions
under which pointwise convergence of a sequence of random functions can be
translated to uniform convergence.

######  Theorem 3.

Let $(\Omega,d)$ be a compact metric space,
$M_{n}:\Omega\rightarrow\mathbb{R}$ be a sequence of random functions and
$M:\Omega\rightarrow\mathbb{R}$ be a fixed function such that

  1. 1.

(Continuity): $M$ is continuous.

  2. 2.

(Stochastic equicontinuity): There exist some $B_{n}=O_{P}(1)$ such that for
all
$\omega,\omega^{\prime}\in\Omega,\left|M_{n}(\omega)-M_{n}(\omega^{\prime})\right|\leq
B_{n}d(\omega,\omega^{\prime})$.

  3. 3.

(Pointwise convergence): For every $\omega\in\Omega$,
$M_{n}(\omega)-M(\omega)=o_{P}(1)$.

Then

$\sup_{\omega\in\Omega}\left|M_{n}(\omega)-M(\omega)\right|=o_{P}(1).$

Further, Theorem 5.7 of van der Vaart (1998) can be used to translate uniform
convergence of a sequence of objective functions to convergence of the
minimizers.

######  Theorem 4.

Let $(\Omega,d)$ be a metric space, $M_{n}:\Omega\rightarrow\mathbb{R}$ be a
sequence of random functions and $M:\Omega\rightarrow\mathbb{R}$ be a fixed
function such that for every $\varepsilon>0$,

$\displaystyle\sup_{\omega\in\Omega}\left|M_{n}(\omega)-M(\omega)\right|=o_{P}(1)$
$\displaystyle\inf_{\omega:d(\omega_{0},\omega)>\varepsilon}M(\omega)>M(\omega_{0}).$

Then any sequence of estimators $\hat{\omega}_{n}$ with
$M_{n}(\hat{\omega}_{n})\leq M_{n}(\omega_{0})+o_{P}(1)$ converges in
probability to $\omega_{0}$.

The first condition of this theorem states that $M_{n}$ converges uniformly to
$M$ while the second condition states that there exists a well-separated
minimizer $\omega_{0}$ of $M$.

###  A.2 Consistency of the MRF estimator

We start with a lemma regarding Lipschitz continuity of the distance and
squared distance functions.

######  Lemma 1.

Let $(\Omega,d)$ be a compact metric space and $\omega_{0}\in\Omega$, then the
maps $\omega\mapsto d(\omega_{0},\omega)$ and $\omega\mapsto
d(\omega_{0},\omega)^{2}$ are Lipshitz continuous.

###### Proof.

Starting with the first map, let $\omega,\omega^{\prime}\in\Omega$, then we
have by the triangle inequality that

$\displaystyle d(\omega_{0},\omega)\leq
d(\omega_{0},\omega^{\prime})+d(\omega^{\prime},\omega)$
$\displaystyle\Rightarrow
d(\omega_{0},\omega)-d(\omega_{0},\omega^{\prime})\leq
d(\omega^{\prime},\omega)$ $\displaystyle d(\omega_{0},\omega^{\prime})\leq
d(\omega_{0},\omega)+d(\omega^{\prime},\omega)$
$\displaystyle\Rightarrow-\left(d(\omega_{0},\omega)-d(\omega_{0},\omega^{\prime})\right)\leq
d(\omega^{\prime},\omega),$

and hence

$\left|d(\omega_{0},\omega)-d(\omega_{0},\omega^{\prime})\right|\leq
d(\omega^{\prime},\omega).$

For the second map, we have

$\displaystyle\left|d(\omega_{0},\omega)^{2}-d(\omega_{0},\omega^{\prime})^{2}\right|$
$\displaystyle=\left|d(\omega_{0},\omega)+d(\omega_{0},\omega^{\prime})\right|\left|d(\omega_{0},\omega)-d(\omega_{0},\omega^{\prime})\right|$
$\displaystyle\leq 2{\textrm{diam}}\,\Omega\,d(\omega,\omega^{\prime}),$

where the first term was bounded by the diameter of $\Omega$ and we used the
Lipshitz continuity of $\omega\mapsto d(\omega_{0},\omega)$ for the second
term. ∎

Next, we prove that the random forest objective function converges uniformly
in $\Omega$ to the conditional variance. Recall definitions (12) of
$M_{n}^{\star}(\omega;x)$ and (2) of $M(\omega;x)$.

###### Proof of Theorem 1.

This proof is done by verifying that for a fixed $x\in[0,1]^{d}$, the maps
$\omega\mapsto M_{n}^{\star}(\omega;x)$ and $\omega\mapsto M(\omega;x)$
respect the conditions of Theorem 3.

#### Continuity

This follows directly from Lemma 1 since for
$\omega,\omega^{\prime}\in\Omega$,

$\displaystyle\left|M(\omega;x)-M(\omega^{\prime};x)\right|$
$\displaystyle=\left|\mathbb{E}\left[d(\omega,Y)^{2}\mid
X=x\right]-\mathbb{E}\left[d(\omega^{\prime},Y)^{2}\mid X=x\right]\right|$
$\displaystyle\leq\mathbb{E}\left[\left|d(\omega,Y)^{2}-d(\omega^{\prime},Y)^{2}\right|\mid
X=x\right]$ $\displaystyle\leq
2{\textrm{diam}}\,\Omega\,d(\omega,\omega^{\prime}),$

showing that $M$ is Lipschitz continuous and hence continuous.

#### Stochastic Equicontinuity

Similarly to the continuity of $M$, the stochastic equicontinuity of
$M_{n}^{\star}$ is inherited. Let $\omega,\omega^{\prime}\in\Omega$,

$\displaystyle\left|M_{n}^{\star}(\omega;x)-M_{n}^{\star}(\omega^{\prime};x)\right|$
$\displaystyle=\left|\sum_{i=1}^{n}w_{i}^{\star}(x)d(\omega,Y_{i})^{2}-\sum_{i=1}^{n}w_{i}^{\star}(x)d(\omega^{\prime},Y_{i})^{2}\right|$
$\displaystyle\leq\sum_{i=1}^{n}w_{i}^{\star}(x)\left|d(\omega,Y_{i})^{2}-d(\omega^{\prime},Y_{i})^{2}\right|$
$\displaystyle\leq\sum_{i=1}^{n}w_{i}^{\star}(x)2{\textrm{diam}}\,\Omega\,d(\omega,\omega^{\prime})$
$\displaystyle=2{\textrm{diam}}\,\Omega\,d(\omega,\omega^{\prime}).$

#### Pointwise convergence

Let $\omega\in\Omega$ be fixed. We need to show that
$M_{n}^{\star}-M=o_{P}(1)$. Recall the contribution of a single fitted tree
(10) to the theoretical random forest objective function $M_{n}^{\star}$,

$T(\omega,x;\xi,D)=\sum_{i=1}^{n}w_{i}(x;\xi,D)d(\omega;Y_{i})^{2}.$

We simplify the notation and remove references to $x,\omega,\xi$ and $D$ in
the remaining of this sub-proof. We now define the Hájek projection
$\accentset{\circ}{M_{n}^{\star}}$ of $M_{n}^{\star}$,

$\accentset{\circ}{M_{n}^{\star}}=\mathbb{E}\left[M_{n}^{\star}\right]+\sum_{i=1}^{n}\mathbb{E}\left[M_{n}^{\star}\mid(X_{i},Y_{i})\right]-\mathbb{E}\left[M_{n}^{\star}\right].$

This can be rewritten in terms of expectations of the individual trees as

$\accentset{\circ}{M_{n}^{\star}}=\mathbb{E}\left[T\right]+\frac{s}{n}\sum_{i=1}^{n}\mathbb{E}\left[T\mid(X_{i},Y_{i})\right]-\mathbb{E}\left[T\right].$

We proceed with the proof of pointwise convergence by showing the stronger
result of $L_{2}$ convergence of $M_{n}^{\star}$, that is
$\mathbb{E}\left[\left(M_{n}^{\star}-M\right)^{2}\right]=o(1)$. We do this via
the convergence of the Hájek projection of the forest weights $M_{n}$, in a
similar fashion to the analysis in Wager and Athey (2018). The mean squared
error $\mathbb{E}\left[\left(M_{n}^{\star}-M\right)^{2}\right]$ can be upper
bounded by the following sum in which each term is more amenable to analysis:

$\mathbb{E}\left[\left(M_{n}^{\star}-\accentset{\circ}{M_{n}^{\star}}\right)^{2}\right]+\mathbb{E}\left[\left(\accentset{\circ}{M_{n}^{\star}}-\mathbb{E}\left[M_{n}^{\star}\right]\right)^{2}\right]+\left(\mathbb{E}\left[M_{n}^{\star}\right]-M\right)^{2}.$

By the specification of the construction of the random forest, we have that
Lemma 7 in Wager and Athey (2018) applies and we can bound the first summand
as

$\mathbb{E}\left[\left(M_{n}^{\star}-\accentset{\circ}{M_{n}^{\star}}\right)^{2}\right]\leq\frac{s^{2}}{n^{2}}\textrm{Var}\left[T\right].$

As shown in the proof of Theorem 5 from Wager and Athey (2018), the variance
of a single tree is asymptotically bounded

$\displaystyle\textrm{Var}\left[T\right]\lesssim
k\textrm{Var}\left[d(\omega,Y)^{2}\mid X=x\right]<\infty,$

where $k$ is the depth at which the trees are grown, as specified in Section
4. Hence

$\mathbb{E}\left[\left(M_{n}^{\star}-\accentset{\circ}{M_{n}^{\star}}\right)^{2}\right]\lesssim\frac{s^{2}}{n^{2}}=o(1)$

For the second summand, we note that
$\mathbb{E}\left[M_{n}^{\star}\right]=\mathbb{E}\left[\accentset{\circ}{M_{n}^{\star}}\right]$
and hence the second term is the variance of the Hájek projection. We can
again use the fact that the variance of a single tree is asymptotically
bounded to obtain

$\textrm{Var}\left[\accentset{\circ}{M_{n}^{\star}}\right]=\frac{s^{2}}{n}\textrm{Var}\left[\mathbb{E}\left[T\mid
Z_{1}\right]\right]\leq\frac{s}{n}\textrm{Var}\left[T\right]\lesssim\frac{s}{n}=o(1)$

Consider at last the third summand. Since $M_{n}^{\star}$ is constructed as
the mean of identically distributed trees, we have that
$\mathbb{E}\left[M_{n}^{\star}\right]=\mathbb{E}\left[T\right]$. We thus
analyze instead the quantity

$\mathbb{E}\left[T\right]-\mathbb{E}\left[d(\omega,Y)^{2}\mid X=x\right].$

As done in the proof of Theorem 3 in Wager and Athey (2018), we have by
honesty that

$\displaystyle\mathbb{E}\left[T\right]-\mathbb{E}\left[d(\omega,Y)^{2}\mid
X=x\right]$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[\mathbb{E}\left[d(\omega,Y)^{2}\mid X\in
L(x)\right]-\mathbb{E}\left[d(\omega,Y)^{2}\mid X=x\right]\right],$

where $L(x)=L(x;\xi,D)$. Using Assumption 2, $M$ is Lipschitz continuous with
constant $C(\omega)$ and we get

$\left|\mathbb{E}\left[d(\omega,Y)^{2}\mid X\in
L(x)\right]-\mathbb{E}\left[d(\omega,Y)^{2}\mid X=x\right]\right|\leq
C(\omega){\textrm{diam}}\,L(x).$

Using now that the trees are constructed following Specification 1 of Athey et
al. (2019), Lemma 2 of Wager and Athey (2018) applies, implying that the size
of the leaves are $o_{P}(1)$. Hence

$\left|\mathbb{E}\left[T-\mathbb{E}\left[d(\omega,Y)^{2}\mid
X=x\right]\right]\right|\leq\mathbb{E}\left[C(\omega)L(x)\right]=\mathbb{E}\left[o_{P}(1)\right]=o(1)$

Altogether, we thus have that
$\mathbb{E}\left[\left(M_{n}^{\star}-M\right)^{2}\right]=o(1)$ and hence we
have for any fixed $(\omega,x)$ the converge
$M_{n}^{\star}(\omega;x)-M(\omega;x)=o_{P}(1)$. ∎

We can now prove that the random forest estimator is a pointwise constistent
estimator of the Fréchet mean.

###### Proof of Theorem 2.

Since $\Omega$ is compact and $m_{n}^{\star}(x)$ minimizes
$M_{n}^{\star}(\cdot;x)$, we have by definition that
$M_{n}^{\star}(m_{n}^{\star}(x);x)\leq M_{n}^{\star}(m(x);x)$. Furthermore, by
Theorem 1, the random forest objective function $M_{n}^{\star}$ converges
uniformly to $M$. Together with the assumption that the Fréchet mean is well-
separated, we have by Theorem 4 that $m_{n}^{\star}(x)$ is a consistent
estimator of $m(x)$. ∎

###  A.3 Consistency of the splitting rule

We finally prove consistency of our medoid-based splitting rule. The proof
follows a similar structure as the proof of Theorem 2.

###### Proof of Proposition 1.

Let $\Omega_{n}={\left\\{Y_{1},\ldots,Y_{n}\right\\}}$ and $\varepsilon>0$. By
assumption, there exists an $\alpha>0$ such that
$\mathbb{P}\left[d(Y,\omega_{\oplus})<\varepsilon\right]>\alpha$. We then have
by independence that

$\mathbb{P}\left[d(\Omega_{n},\omega_{\oplus})>\varepsilon\right]=\mathbb{P}\left[d(Y,\omega_{\oplus})>\varepsilon\right]^{n}<(1-\alpha)^{n}\rightarrow
0,$

showing that $d(\Omega_{n},\omega_{\oplus})\rightarrow 0$ in probability.

Let $M(\omega)=\mathbb{E}\left[d(\omega,Y)^{2}\right]$ be the unweighted
population Fréchet function, minimized by $\omega_{\oplus}$, and
$M_{n}(\omega)=\frac{1}{n}\sum_{i=1}^{n}d(\omega,Y_{i})^{2}$ be the unweighted
empirical Fréchet functions with a sequence of minimizers
${\left\\{\hat{\omega}_{n}\right\\}}$. We start by showing that $M$ is
continuous, which follows directly from the Lipschitz continuity of the
squared distance since for $\omega,\omega^{\prime}\in\Omega$,

$\displaystyle\left|M(\omega)-M(\omega^{\prime})\right|$
$\displaystyle=\left|\mathbb{E}\left[d(\omega,Y)^{2}\right]-\mathbb{E}\left[d(\omega^{\prime},Y)^{2}\right]\right|$
$\displaystyle\leq\mathbb{E}\left[\left|d(\omega,Y)^{2}-d(\omega^{\prime},Y)^{2}\right|\right]$
$\displaystyle\leq 2{\textrm{diam}}\,\Omega\,d(\omega,\omega^{\prime}).$

Furthermore, we can also use the Lipschitz continuity of the squared distance
to show that $M_{n}$ is stochastic equicontinuous since for each
$\omega,\omega^{\prime}\in\Omega$,

$\displaystyle\left|M_{n}(\omega)-M_{n}(\omega^{\prime})\right|$
$\displaystyle=\left|\sum_{i=1}^{n}d(\omega,Y_{i})^{2}-\sum_{i=1}^{n}d(\omega^{\prime},Y_{i})^{2}\right|$
$\displaystyle\leq\sum_{i=1}^{n}\left|d(\omega,Y_{i})^{2}-d(\omega^{\prime},Y_{i})^{2}\right|$
$\displaystyle\leq\sum_{i=1}^{n}2{\textrm{diam}}\,\Omega\,d(\omega,\omega^{\prime})$
$\displaystyle=2{\textrm{diam}}\,\Omega\,d(\omega,\omega^{\prime}).$

Finally, by the weak law of large numbers, $M_{n}$ converges pointwise to $M$,
which verifies the conditions of Theorem 3. This shows that $\left\lVert
M_{n}-M\right\rVert_{\Omega}=o_{P}(1)$. Together with the assumption of well-
separatedness of $\omega_{\oplus}$, this shows that $M_{n}$ fulfills the
conditions of Theorem 4, giving that
$d(\hat{\omega}_{n},\omega_{\oplus})=o_{P}(1)$.

Let $\varepsilon>0$ and $K=2{\textrm{diam}}\,\Omega$, under the event
$E_{n}={\left\\{d(\Omega_{n},\omega_{\oplus})\leq\varepsilon/2K\right\\}}\cap{\left\\{d(\hat{\omega}_{n},\omega_{\oplus})\leq\varepsilon/2K\right\\}}$,
there exists for each $n\in\mathbb{N}$ an $\omega^{\dagger}_{n}\in\Omega_{n}$
such that $d(\omega^{\dagger}_{n},\omega_{\oplus})\leq\varepsilon/2K$. Hence,

$\displaystyle\min_{\omega\in\Omega_{n}}M_{n}(\omega)-\min_{\omega\in\Omega}M_{n}(\omega)$
$\displaystyle=M_{n}(\tilde{\omega}_{n})-M_{n}(\hat{\omega}_{n})\leq
M_{n}(\omega^{\dagger}_{n})-M_{n}(\hat{\omega}_{n})$ $\displaystyle\leq
Kd(\omega^{\dagger}_{n},\hat{\omega}_{n})\leq
K(d(\omega^{\dagger}_{n},\omega_{\oplus})+d(\omega_{\oplus},\hat{\omega}_{n}))\leq\varepsilon.$

By the two first results in this proof,
$\mathbb{P}\left[E_{n}\right]\rightarrow 1$ and we obtain that
$M_{n}(\tilde{\omega}_{n})-M_{n}(\hat{\omega}_{n})=o_{P}(1)$, and therefore

$M_{n}(\tilde{\omega}_{n})=M_{n}(\hat{\omega}_{n})+o_{P}(1)\leq
M_{n}(\omega_{\oplus})+o_{P}(1).$

This shows that the sequence of minimizers
${\left\\{\tilde{\omega}_{n}\right\\}}$ fulfills the third condition of
Theorem 4, hence $d(\tilde{\omega}_{n},\omega_{\oplus})=o_{P}(1)$. ∎

## References

  * Athey et al. (2019) Athey, S., Tibshirani, J., Wager, S., 2019.  Generalized random forests.  The Annals of Statistics 47, 1148 – 1178.  URL: https://doi.org/10.1214/18-AOS1709, doi:10.1214/18-AOS1709.
  * Best and Chakravarti (1990) Best, M.J., Chakravarti, N., 1990\.  Active set algorithms for isotonic regression; A unifying framework.  Mathematical Programming 47, 425–439.  URL: https://doi.org/10.1007/BF01580873, doi:10.1007/BF01580873.
  * Biau et al. (2008) Biau, G., Devroye, L., Lugosi, G., 2008.  On the Performance of Clustering in Hilbert Spaces.  IEEE Transactions on Information Theory 54, 781–790.  URL: http://ieeexplore.ieee.org/document/4439834/, doi:10.1109/TIT.2007.913516.
  * Breiman (2001) Breiman, L., 2001.  Random Forests.  Machine Learning 45, 5–32.  URL: https://doi.org/10.1023/A:1010933404324, doi:10.1023/A:1010933404324.
  * Breiman (2004) Breiman, L., 2004.  Consistency for a simple model of random forests.  University of California at Berkeley. Technical Report 670.
  * Breiman et al. (1984) Breiman, L., Friedman, J., Stone, C.J., Olshen, R.A., 1984\.  Classification and Regression Trees.  Taylor & Francis.
  * Capitaine et al. (2020) Capitaine, L., Bigot, J., Thiébaut, R., Genuer, R., 2020\.  Fréchet random forests for metric space valued regression with non euclidean predictors URL: http://arxiv.org/abs/1906.01741.
  * Davis et al. (2010) Davis, B.C., Fletcher, P.T., Bullitt, E., Joshi, S., 2010\.  Population Shape Regression from Random Design Data.  International Journal of Computer Vision 90, 255–266.  URL: https://doi.org/10.1007/s11263-010-0367-1, doi:10.1007/s11263-010-0367-1.
  * Fréchet (1948) Fréchet, M., 1948.  Les éléments aléatoires de nature quelconque dans un espace distancié, in: Annales de l’institut Henri Poincaré, pp. 215–310.  Volume 10, issue 4.
  * Fu et al. (2021) Fu, G., Dai, X., Liang, Y., 2021.  Functional random forests for curve response.  Scientific Reports 11, 24159\.  URL: https://www.nature.com/articles/s41598-021-02265-4, doi:10.1038/s41598-021-02265-4.
  * Ghodrati and Panaretos (2022) Ghodrati, L., Panaretos, V.M., 2022\.  Distribution-on-distribution regression via optimal transport maps.  Biometrika 109, 957–974.  URL: https://doi.org/10.1093/biomet/asac005, doi:10.1093/biomet/asac005.
  * Hein (2009) Hein, M., 2009.  Robust Nonparametric Regression with Metric-Space Valued Output, in: Advances in Neural Information Processing Systems, Curran Associates, Inc.  URL: https://proceedings.neurips.cc/paper_files/paper/2009/hash/92977ae4d2ba21425a59afb269c2a14e-Abstract.html.
  * Hinkle et al. (2012) Hinkle, J., Muralidharan, P., Fletcher, P.T., Joshi, S., 2012\.  Polynomial Regression on Riemannian Manifolds, in: Fitzgibbon, A., Lazebnik, S., Perona, P., Sato, Y., Schmid, C. (Eds.), Computer Vision – ECCV 2012, Springer, Berlin, Heidelberg. pp. 1–14.  doi:10.1007/978-3-642-33712-3_1.
  * Meinshausen (2006) Meinshausen, N., 2006.  Quantile Regression Forests.  Journal of Machine Learning Research 7, 983–999.  URL: http://jmlr.org/papers/v7/meinshausen06a.html.
  * Mentch and Hooker (2016) Mentch, L., Hooker, G., 2016\.  Quantifying uncertainty in random forests via confidence intervals and hypothesis tests.  Journal of Machine Learning Research 17, 1–41.  URL: http://jmlr.org/papers/v17/14-168.html.
  * Müller (2016) Müller, H.G., 2016.  Peter Hall, functional data analysis and random objects.  The Annals of Statistics 44, 1867 – 1887.  URL: https://doi.org/10.1214/16-AOS1492, doi:10.1214/16-AOS1492.
  * Nerini and Ghattas (2007) Nerini, D., Ghattas, B., 2007\.  Classifying densities using functional regression trees: Applications in oceanology.  Computational Statistics & Data Analysis 51, 4984–4993.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0167947306003550, doi:10.1016/j.csda.2006.09.028.
  * Newey (1991) Newey, W.K., 1991.  Uniform Convergence in Probability and Stochastic Equicontinuity.  Econometrica 59, 1161–1167.  URL: https://www.jstor.org/stable/2938179?origin=crossref, doi:10.2307/2938179.
  * Panaretos and Zemel (2016) Panaretos, V.M., Zemel, Y., 2016\.  Amplitude and phase variation of point processes.  The Annals of Statistics 44, 771 – 812.  URL: https://doi.org/10.1214/15-AOS1387, doi:10.1214/15-AOS1387.
  * Panaretos and Zemel (2020) Panaretos, V.M., Zemel, Y., 2020\.  An Invitation to Statistics in Wasserstein Space.  SpringerBriefs in Probability and Mathematical Statistics, Springer International Publishing, Cham.  URL: http://link.springer.com/10.1007/978-3-030-38438-8, doi:10.1007/978-3-030-38438-8.
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., Duchesnay, E., 2011\.  Scikit-learn: Machine Learning in Python.  Journal of Machine Learning Research 12, 2825–2830.  URL: http://jmlr.org/papers/v12/pedregosa11a.html.
  * Pelletier (2006) Pelletier, B., 2006.  Non-parametric regression estimation on closed Riemannian manifolds.  Journal of Nonparametric Statistics 18, 57–67.  doi:10.1080/10485250500504828. publisher: Taylor & Francis.
  * Petersen and Müller (2019) Petersen, A., Müller, H.G., 2019\.  Fréchet regression for random objects with Euclidean predictors.  The Annals of Statistics 47, 691–719.  URL: https://doi.org/10.1214/17-AOS1624, doi:10.1214/17-AOS1624.
  * Qiu et al. (2022) Qiu, R., Yu, Z., Zhu, R., 2022\.  Random Forests Weighted Local Fréchet Regression with Theoretical Guarantee URL: http://arxiv.org/abs/2202.04912.
  * Ramos-Carreño et al. (2019) Ramos-Carreño, C., Suárez, A., Torrecilla, J.L., Carbajo Berrocal, M., Marcos Manchón, P., Pérez Manso, P., Hernando Bernabé, A., García Fernández, D., Hong, Y., Rodríguez-Ponga Eyriès, P.M., Sánchez Romero, A., Petrunina, E., Castillo, A., Serna, D., Hidalgo, R., 2019.  GAA-UAM/scikit-fda: Functional Data Analysis in Python.  URL: https://doi.org/10.5281/zenodo.7573921, doi:10.5281/zenodo.7573921.
  * Scornet (2016) Scornet, E., 2016.  On the asymptotics of random forests.  Journal of Multivariate Analysis 146, 72–83.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0047259X15001542, doi:10.1016/j.jmva.2015.06.009.
  * Srivastava and Klassen (2016) Srivastava, A., Klassen, E.P., 2016\.  Functional and Shape Data Analysis.  Springer Series in Statistics, Springer New York, New York, NY.  URL: http://link.springer.com/10.1007/978-1-4939-4020-2, doi:10.1007/978-1-4939-4020-2.
  * Thomas Fletcher (2013) Thomas Fletcher, P., 2013.  Geodesic Regression and the Theory of Least Squares on Riemannian Manifolds.  International Journal of Computer Vision 105, 171–185.  URL: http://link.springer.com/10.1007/s11263-012-0591-y, doi:10.1007/s11263-012-0591-y.
  * Tiwari et al. (2022) Tiwari, M., Kang, R., Lee, J., Piech, C., Shomorony, I., Thrun, S., Zhang, M.J., 2022\.  Mabsplit: Faster forest training using multi-armed bandits, in: Koyejo, S., Mohamed, S., Agarwal, A., Belgrave, D., Cho, K., Oh, A. (Eds.), Advances in Neural Information Processing Systems, Curran Associates, Inc.. pp. 1223–1237.  URL: https://proceedings.neurips.cc/paper_files/paper/2022/file/08857467641ad82f635023d530605b4c-Paper-Conference.pdf.
  * Townsend et al. (2016) Townsend, J., Koep, N., Weichwald, S., 2016.  Pymanopt: A Python Toolbox for Optimization on Manifolds using Automatic Differentiation.  Journal of Machine Learning Research 17, 1–5.  URL: http://jmlr.org/papers/v17/16-177.html.
  * Tucker et al. (2013) Tucker, J.D., Wu, W., Srivastava, A., 2013.  Generative models for functional data using phase and amplitude separation.  Computational Statistics & Data Analysis 61, 50–66.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0167947312004227, doi:10.1016/j.csda.2012.12.001.
  * van der Vaart (1998) van der Vaart, A.W., 1998.  Asymptotic statistics.  Cambridge series in statistical and probabilistic mathematics, Cambridge University Press, Cambridge, UK ; New York, NY, USA.  doi:10.1017/CBO9780511802256.
  * Wager and Athey (2018) Wager, S., Athey, S., 2018.  Estimation and inference of heterogeneous treatment effects using random forests.  Journal of the American Statistical Association 113, 1228–1242.  URL: https://doi.org/10.1080/01621459.2017.1319839, doi:10.1080/01621459.2017.1319839.
  * Wager et al. (2014) Wager, S., Hastie, T., Efron, B., 2014.  Confidence intervals for random forests: The jackknife and the infinitesimal jackknife.  Journal of Machine Learning Research 15, 1625–1651.  URL: http://jmlr.org/papers/v15/wager14a.html.
  * Yuan et al. (2012) Yuan, Y., Zhu, H., Lin, W., Marron, J.S., 2012.  Local polynomial regression for symmetric positive definite matrices.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 74, 697–719.  URL: https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1467-9868.2011.01022.x, doi:10.1111/j.1467-9868.2011.01022.x.