Fréchet Geodesic Boosting
Authors: Yidong Zhou
Su I Iao
Hans-Georg Müller
Venue: NeurIPS 2025 poster ()
Keywords: geodesic metric space
gradient boosting
non-Euclidean output
Wasserstein space
Primary area: general_machine_learning
TL;DR: We propose FGBoost, a gradient boosting framework designed to intrinsically model complex regression relationships with non-Euclidean outputs in geodesic metric spaces.

Abstract: Gradient boosting has become a cornerstone of machine learning, enabling base learners such as decision trees to achieve exceptional predictive performance. While existing algorithms primarily handle scalar or Euclidean outputs, increasingly prevalent complex-structured data, such as distributions, networks, and manifold-valued outputs, present challenges for traditional methods. Such non-Euclidean data lack algebraic structures such as addition, subtraction, or scalar multiplication required by standard gradient boosting frameworks. To address these challenges, we introduce Fréchet geodesic boosting (FGBoost), a novel approach tailored for outputs residing in geodesic metric spaces. FGBoost leverages geodesics as proxies for residuals and constructs ensembles in a way that respects the intrinsic geometry of the output space. Through theoretical analysis, extensive simulations, and real-world applications, we demonstrate the strong performance and adaptability of FGBoost, showcasing its potential for modeling complex data.

Fréchet Geodesic Boosting
Yidong Zhou"
Su I Iao' Department of Statistics University of California, Davis Davis, CA 95616 siao@ucdavis. edu
Department of Statistics University of California, Davis
Davis, CA 95616 ydzhou@ucdavis.com edu
Hans-Georg Müller Department of Statistics University of California, Davis
Davis, CA 95616
hgmueller@uccaviis.ede
Abstract
Gradient boosting has become a cornerstone of machine learning, enabling base learners such decision trees achieve exceptional predictive performance. While existing algorithms primarily handle scalar or Fuclidean outputs, increasingly prevalent complex-structuree data, such as distributions, networks. and manifold- valued outputs, present challenges for traditional methods. Such non non-Euche Euclidean data lack algebraic structures such addition, subtraction, or scalar multiplication required by standard gradient boosting frameworks. To address these challenges, we introduce Fréchet geodesic boosting (FGBoost), novel approach tailored for outputs residing in geodesic metric spaces. FGBoost leverages geodesics as proxies for residuals and constructs ensembles way that respects the intrinsic geometry of the output space. Through theoretical analysis, extensive simulations, and real- world applications, we demonstrate the strong performance and adaptability of FGBoost, showcasing its potential for modeling complex data.
Introduction
Boosting has emerged as one of the most influential learning paradigms, enabling base learners, such as decision trees, to achieve superior predictive performance. The foundational idea of boosting can understood as a functional gradient descent method applied to a cost function in function space ( Over two decades ago, explicit gradient boosting algorithms were introduced [3 19). laying the groundwork for their widespread success. Modern iterations of gradient boosting, such as XGBoost [13 and LightGBM (27), have further advanced the field, providing highly efficient and scalable solutions for numerous machine learning tasks.
The increasing availability of complex-structure data in modern science has introduced significant challenges for conventional learning methods 99. Examples of such data include functional data [59], networks [64], trees [44, distributions [48) and data residing on manifolds such as symmetric
positive-definite matrices [46]. These types of data inherently -uclldedd and can be viewed as random objects located in metric spaces equipped with suitable metrics. A notable example the Wasserstein space, where elements are probability distributions and distances are measured using the Wasserstein metric [5]. Despite the success of gradient boosting, existing algorithms are
"The first two authors contributed equally to this work.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

predominantly designed for scalar or Euclidean outputs and cannot handle outputs in general metric spaces due the absence of algebraic operations such as addition or scalar multiplication.
To address this gap, we propose Fréchet geodesic boosting (FGBoost), a novel framework designed to adapt gradient boosting for non-Fuclidean outputs, specifically for data residing in geodesic metric spaces. FGBoost enables modeling complex regression relationships between Euclidean predictors and non-Euclidean outputs by leveraging the intrinsic geometry of the output space.
1.1 Contributions
The primary contributions of this work are follows:
Methodology. We address the challenge of working in geodesic metric spaces, which lack the linear structure required for standard gradient boosting. FGBoost introduces geodesics proxies for residuals and iteratively constructs an ensemble by adding geodesies while preserving the geometric properties of the output. To achieve this, we develop novel geometric definitions that ensure FGBoost operates intrinsically and adheres to the underlying geometry of the output space. Furthermore, develop
version of Shapley Additive Explanations (SHAP) values 37 to enhance the interpretability of FGBoost. To the best of our knowledge, FGBoost represents the first boosting framework designed to effectively accommodate general non-Euclidean outputs.
Theoretical analysis. We introduce a general framework study the theoretical properties of FGBoost, requiring only that the output space is Hadamard space [56) and demonstrate that the loss function strongly convex and Lipschitz continuous, guaranteeing the existence and uniqueness of the solution. Using empirical process theory (58), we show that the empirical risk functional converges uniformly to its population counterpart, and the corresponding minimizer is consistent. Detailed proofs of theoretical results are provided in the appendix.
Simulation studies. Through extensive numerical experiments, we evaluate the performance of FGBoost across various types of non-Euclidean outputs, including distributions, networks, and compositional data. The results reveal the superiority of FGBoost over existing regression methods designed for non Euclidean outputs and demonstrate its adaptability to diverse data structures.
Experiments on real-world data. We validate the practical utility of FGBoost using real-world datasets from multiple domains. These include distributional data from human mortality studies, networks derived from New York City yellow taxi trip records, and compositional data from survey of unemployed workers in New Jersey. These applications highlight the ability of FGBoost to effectively model complex data and its relevance across many fields.
1.2 Related work
Gradient boosting. Recent ad advancements in gradient boosting have focused on extending its applicability to handle complex-structue
data. Examples include algorithms for censored survival data [23 32], functional data [8] 57 8), and online learning scenarios (3) These methods rely on well-defined loss functions (e.g., Cox partial likelihood) or a linear space structure (c.g., Hilbert spaces), enabling adaptation within the gradient boosting framework. Efforts have also been made to incorporate predictive uncertainty, e.g., within the framework of evidential learning [40] In particular, 40 leverages Wasserstein geometry to construct posterior distributions as intermediate targets for scular regression. By contrast, FGBoost directly learns regression maps into general geodesic spaces, enabling prediction for wide range of non-Buclidean outputs such as distributions, networks, SPD matrices, and compositional data.
Regression models for non-Euclidean outputs. Recent years have seen a surge regression methods for non-Euclidean outputs. Early approaches include Euclidean embeddings with distance matrices [17 and Madaraya-Watson kernel regression [22]. More recently, Fréchet regression [7] extended linear and nonparametric regression to metric space-valued outputs. To handle high- dimensional predictors, extensions have incorporated sufficient dimension reduction [60 62), single
index models [4] (0), principal component regression (54), and deep neural networks 25 However, these methods often depend on restrictive assumptions, such as linear or single-index structures, or low-dimensional manifold constraints. Random forest algorithms have also been adapted for non-Fuclidean outputs 10, 49]. In simulations and real-world applications, FGBoost demonstrates superior performance against these alternatives.
2

---

2 Preliminaries on metric geometry
Let (M,d) d) be a bounded metric space. A curve in M continuous map 7 [a, + M with sup Eizo
length L(7) suu
((t+1)), where the supremum is taken over all possible partitions of the interval [a, b] with arbitrary breakpoints
to t < .... S ti b. Two curves "/1 and ? are considered equivalent if there exist non-decreasing. continuous reparametrizations Φ1 and d2 such that 7/1 0h 7/2 02. In this case, 71 is said be reparametrisation of 7 and one has that I(n) = L(22)- A curve y [a, b] - M is said to have constant speed if for all a S 8 S t S b, [(((s.)) b-o I(7), where "[x.d] denotes the restriction of y to (s, t]. By construction, the metric d(a, 8) is always less than or equal the length of any curve connecting a and B. A metric space M is called length space if for all a .M:
d(a, ()) inf L(7);
where the infimum is taken over all curves y connecting a to 3. A length space is a geodesic space for all a, 3 .M the infimum on the right-hand side of (I) is attained.
In a geodesic space, a geodesic between two points a and 3 is defined as any constant speed curve 7: [0, 1] - .M that achieves the infimum in ① This geodesic is denoted as 7a.8- If there exists only one such geodesic for all a, 3 E M, the space M is a unique geodesic space 7
Definition 1. For a,r, M and [0, 1), define the following simple operations on geodesics,
Yo, 2 TC,8 == To,,8, SYa.B == "i V @ Ya.g = (ta,p(t): t E v0, id, == Ya.a.
These operations generalize the notions of addition, reversal, scalar multiplication, and zero from vectors to geodesics. The following example spaces frequently arise in real-world applications and will feature in simulations and real-world data applications. Additionally, the space of compositional data is discussed in Appendix A
Example 1 (Univariate probability distributions). Consider the Wasserstein space (W, dw), which consists of probability distributions on R. with finite second moments, equipped with the Wasserstein metric dw. This space both complete and separable [45]. The 2- 2-Wasserstein metric between two distributions μr and μtz is dry = (F=, (F) F
{p)} dp, where F,--! and Fyzl are the quantile functions of μr and U2,, respectively. This space offers a natural framework for analyzing distributions as geometric objects, with geodesics explicitly characterized through optimal transport maps. Denote by TH the pushforward measure of p. by the transport т. The geodesic connecting two distributions pl, N1,322 E W is given by McCann's 's interpolant (41]:
We (t) {id + t(F 1
Fm ffm--id idd] #MI E [0,1].
where id denotes the identity map and Fus is the cumulative distribution function of pt.
Example (Networks). Consider the space of simple, undirected, weighted networks with a fixed number of nodes and bounded edge weights. Each network can be represented uniquely by its graph Laplacian. The space of graph Laplacians equipped with the Frobenius metric can thus be used to characterize the space of networks [28 53 64). For any two graph Laplacians a, B, the geodesic connecting them the line segment, i.e., (t, = a++ + (B 188-ay. at)..
Example (Symmetric positive-definite matrices). Consider the space of x symmetric positive- definite matrices Sym," Common examples include covariance and correlation matrices, which play a crucial role many statistical and data analysis tasks. Depending on the application, different metrics have been proposed to equip Sym," with geometric structure. including the basic Frobenius affine-invariant (15] and the Pog-Cholesky metric (33). Under any of these metrics, Sym forms unique geodesic space.
metric as well as more advanced metrics such as the affine-inveer invariant metric (46), the power metric
Methodology
3.1 Problem formulation
Consider a unique geodesic space (M, d). Let (X,Y) be a random pair in RP x M. Suppose {(X,, ) form sample that consists of n independent realizations of ().)). FGBoost
3

---

presents novel approach to model the relationship between the non-Euclidean output Y E M and multivariate predictor X IR..
For a random object E M,, the Fréchet mean of Y, extending the usual notion of mean, is
Eg (Y)
arg min ( (Y,w)):
WE.
where the existence and uniqueness of the minimizer are guaranteed for Hadamard spaces [56] and the example spaces described in Examples
3.2 Fréchet geodesic boosting
iteratively constructs ensemble Fx of K base learners f1...., fk, starting with a constant y, that best fits the data. At each step k the ensemble is updated
For Euclidean outputs, canonical gradient boosting [19]
(z) Yo +o+uf(()))
Fa(x)) Fi-(()) vff(x)), k =
where v E (0. of each base learner to the overall model. Typically, Yo is chosen as the sample mean of (Y,)
1) is a shrinkage parameter, known as the learning rate, which controls the contribution
The central idea gradient boosting is to improve the model by iteratively reducing the remaining error of the current ensemble. At each iteration kt the base learner f& is trained to approximate the negative gradient of the loss function at the current prediction Fk-1- For squared error loss, the negative gradient corresponds to the residuals Y; (X,) Therefore, the base learner fx fitted to the data set {(X,, Y; FE- (X.)) This iterative procedure refines the ensemble, progressively reducing the overall prediction error. Gradient boosting is highly flexible and can employ various base learners, where tree-based models such as decision trees are most common 22
To explore the regression relationship between -un-Eucideean outputs Y E M and Euclidean predictors X R", we propose Fréchet geodesic boosting (FGBoost). In the Euclidean setting, the ensemble model is constructed by sequentially adding base learners, each approximating the residual, which corresponds to the negative gradient of the loss function. For unique geodesic spaces, the concept of residuals can be naturally extended to geodesics. Specifically, the residual is replaced by the geodesic connecting the current prediction and the actual observation. Consequently, the ensemble model in FGBoost is defined the addition of sequence of geodesics. However, geodesic addition is not inherently well-defined unless the geodesics are connected end to end, preserving continuity. To address this challenge, we introduce the following assumption to ensure well-defined operations within the geodesic framework.
Assumption 1. Let (M,d) be unique geodesic space. For any two points a., M, there exists geodesic transport map TT., : :M +> M with the following property: Ty.., (a) = 3 and for any w E M, there exists a unique point Ç E M such that Ty., (w) Ç.
This assumption ensures that any geodesic Ta,9 can be naturally extended from any starting point E M to new endpoint Ç E .M. In the Euclidean space TR. this map straightforward and expressed as T Tr. (w) = + (33 a). This construction extends to Hilbert spaces (e.g., T.2 I)), where geodesics are straight lines connecting points. An analogous principle can be applied for Riemannian manifolds through parallel transport 35) Specific definitions of geodesic transport maps for Examples] 3 are provided in Appendix B Using the geodesic transport map, we now extend the notion of addition to geodesics that are not connected end to end.
Definition 2. For any points a, B. Ç E M with 3 + w, define the addition between two geodesics Ya.8 and Ti 3 No.5" 8 ① "/a,2" -a,'' where Ç¹ T.w.c (3).
The above operation is intuitive in Euclidean space, where Ç' <== = (5 -ww.
For outputs in unique geodesic space M, the base learner fx+1 is trained approximate the geodesic connecting the current prediction to the actual observation, generalizing the concept of residuals the geodesic setting. The initial ensemble for FGBoost is defined as Fo(X,) idy, which corresponds to the geodesic from a fixed reference point Yo E M to itself. In practice, Yo is chosen the sample Fréchet mean of (Y;). (== The ensemble model constructed as the addition of a sequence of geodesics, giving rise to the iteration
Fx(x)) Fx-1(x) @ {u O fo(x)), k k=,,..,,KK

---

Algorithm Fréchet Geodesic Boosting
Input: data {(X,, Y:))==> new predictor level X and learning rate v E (0, Initialize the model with the estimated Fréchet mean of (Yi)==1: Fo(x:) idy, where Yg
arg wEM EE= (Y,,,. for to K do
1. Fit base learner (e.g. tree) fax to approximate the geodesic from the current prediction the actual observation using data {(x;, 7p.- "x x))== current prediction.
where Y Tr, ((x)) denotes the
2. Update the ensemble model: (xx = (x) @ {u & fz())
end for Output: prediction y =F(()) (Y6) where ((X) = (X).
P,(K)) Fo(K)
87 Tr
vof
Fx(X,) =x-1(1)0100))
Yi
Pf
87 Opo
Fo(X)) idy, vof (00)
Yx,y.
Yo
Figure I: Illustration of the general framework for Fréchet geodesic boosting. The algorithm starts with fixed reference point Yo .,, serving the initial estimate. The initial ensemble of Fréchet geodesic boosting the geodesic from Y, to itself, idy. At each step, the base learner fa+1 is trained to approximate the geodesic connecting the current prediction Y:* and actual observation Yi After K iterations, the ensemble model terminates at Fx (xX)),, and the final prediction for Y; is Yx = Tí Txxxxx)
where scalar multiplication and addition of geodesics defined in Definition aand Definition 2 respectively. The ensemble F-) ()) is the geodesic from Yo to the current prediction, which can be expressed as (F.=.(X.) (Y,) using the geodesic transport map. The complete algorithm for FGBoost detailed in Algorithm with schematic illustration provided in Figure
To define suitable loss function for FGBoost, we need well-defined metric for comparing geodesics. geodesic T,.38 in a unique geodesic space is uniquely characterized by its endpoints a and 33. This observation allows to represent the space of geodesics follows:
Definition 3. The space of geodesics on a unique geodesic space (M, d) is
G(M) == ((a,3) ar, M).
Each geodesic 7a,8 is uniquely represented in ((M) as the pair (a, B). To quantify the distance between two geodesics Tor.81 Ta2.22 we define the following metric:
dg( To13)) = / d'(a), + d2 82).
Proposition dg a valid metric on the space of geodesics (MM).
(G(M), dg) is thus metric space, which allows to formally define the loss function for FGBoost.
FGBoost seeks to find an approximation that minimizes the average loss over the training set,
f
arg min
n
5

---

where F maps the predictor X; to geodesic connecting Yo and the prediction. FGBoost starts with an initial constant model Fo and incrementally adds base learners in greedy manner. At iteration k the base learner trained
f arg min Jt
Σ AA (Y)),
(3)
where Fx,) is the current prediction and T P.. ( (vos(())) (Y6) represents the ending point of the updated prediction after incorporating the new base learner.
Although the expression for f* in (3) may appear to involve a nested optimization, in practice FGBoost adopts the same greedy approximation used in classical gradient boosting. This reduces the task to fitting decision trees to pseudo-residuatt At iteration k, we compute the geodesics 7p,-11xx .x,
connecting the current predictions Y the observed responses Y, and treat these geodesics as rscudo-residuals. decision tree is then trained on the pairs (Xt: (y.-- providing a tractable approximation to the idealized optimization problem.
Tree construction follows the standard greedy procedure of decision trees. For any candidate split defined by a feature and threshold, the data are divided into two child regions Rí and Rg Each region assigned representative geodesic obtained the Fréchet mean of the pseudo-residuals contains:
Yj arg min
Σ de (3 Tww) j jE [1,2].
Tho EG(M) S:EEER,
Under the metric dg, this optimization decouples into two simpler Fréchet mean problems: one over the starting points { ky= and one over the endpoints {Yi]. The resulting "3 is therefore the geodesic connecting the two means. This decoupling is purely a computational device for efficient leaf estimation and does not compromise the geometric integrity of the method.
The quality of candidate split is measured by the resulting mean squared error, i.c.. the sum of squared distances between the observed responses Y, and their updated predictions after applying the representative geodesics. The split that yields the greatest reduction in this loss is selected, and the process is repeated recursively until a stopping criterion (e.g., maximum depth) is met.
4 Theoretical analysis
We introduce general framework study the theoretical properties of FGBoost, Let F represent class of base learners f IRP - G(M), where G(M) the space of geodesics as per (2) Define sF as the set of all linear combinations of base learners in F. For any F E span(J) one has F(x) = idy, @ {v1 @ fi())} @ ... {vK @ fK(x)), where Yo is a fixed reference [0, oo) denote the loss func- tion, defined as vi F) = (G(r,F), where dg is the metric on the space of geodesics ((M). FGBoost aims to construct function F Rp G(M) that minimizes the empirical risk func- tional A, (F) = En (/Y%,, F(X.)). The population counterpart of this risk functional is A(F) = F(x(Tx.,,F(X))
point and fū E F for k L...., K. Let v G(M) x G(M)
Remark 1. Throughout this section, we assume that (M, d) is bounded metric space, introduced in Section 2. For spaces that are potentially unbounded, in practice the data distribution is typically supported on stochastically bounded subset, the diameter can be taken as high-probability bound. This boundedness assumption is standard in the analysis of metric space-valued data and is reasonable for practical applications [16 30).
To ensure that the optimization problem is well- polelposed it is crucial that the loss function v satisfies certain desirable properties. These properties are guaranteed when M a Hadamard space.
Definition (Hadamard space). metric space (M, d) is Hadamard space if it is complete and if for each pair of points W1, W2 M there exists point E - with the property that for all points 3 E M:
1 (8) s ) (B.ww) (B( d2 dwl) 2

---

Hadamard spaces, also known as global NPC -Non-Positivy Curvature) spaces [56], are unique geodesic spaces. The spaces discussed in Examples 3 fall within this category. The following proposition establishes key properties of the loss function vp in Hadamard space. Proposition 2. If (M,d) d) Hadamard space, then for any geodesic 7 E g(M), the function (),,)) is strongly convex over (MM) and Lipschitz continuous with respect to dg-
The strong convexity and continuity of de enable us to establish the existence and uniqueness of solution for the risk minimization problem. These properties are crucial for analyzing the asymptotic behavior of gradient boosting algorithms; see for example [63].
Theorem If (M, d) is a Hadamard space, then the optimization problems arg min F'spun()) A(F) and arg min (F)) A,,(F) each admit a unique solution.
To address the challenge posed by the absence of linear operations, we employ tools from empirical process theory [58] to study the asymptotic behavior of the empirical risk functional A, (F) and establish that An (F) converges uniformly to the population risk functional A(F) over F E F as the sample size grows, which then guarantees the convergence of the corresponding minimizer.
Theorem 2. Suppose (M, d) is a Hadamard space. Then sup FEspo((())) |A, (F) A(F) = (1). Furthermore, SuP=ER= dg [F ( F+(x)) = =,((), where F; = arg min respo(() AF(F) and F arg Fispant F) A(F).
5 Numerical experiments
We assess the performance of FGBoost through comprehensive numerical simulations involving non-Fuclidean outputs, specifically distributional data modeled in the Wasserstein space equipped with the Wasserstein metric, and network data represented by graph Laplacians using the Frobe- nius metric, as detailed in Examples [ and 2 Simulations are conducted with sample sizes of 100, 200, 500, 1000, 2000, and each scenario replicated across 500 runs. FGBoost is bench- marked against state-of ootiee--f--h--aaaeeeee the-art art regression models for non- -uc-ldcid outputs, including global Fréchet regression (GFR) 47), sufficient dimension reduction (SDR) [62], single index Fréchet regression (IFR) 4. Fréchet random forest (FRF) [10], and random forest weighted local linear Fréchet re- gression (RFWLLFR) [49. Due to their high computational cost, SDR and IFR are not evaluated at sample size = 2000. A detailed comparison of training times across all models is provided in Appendix Additional simulations for compositional data are presented in Appendix A Code for implementing FGBoost is available at httpps:///itthi ctccauu..comm// wwt.ccoo/Suu... Cttp:/SSIAOOFFFT FFBosst
Common hyperparameters. In all simulations, the learning rate is set 0.05, and the number of iterations K is fixed at 100. The depth of the tree is fixed at 3, with each leaf requiring a minimum of 10 samples. Tuning these parameters can be accomplished through grid search, assessing empirical risk with cross-validation. Additionally, 10% of the training set is reserved as the validation set in each run. The training process halts when the empirical risk on the validation set no longer shows consistent improvement.
Performance evaluation. The out-of-sample performance of FGBoost is assessed using the mean squared prediction eгror (MSPE). Write m(.) for the true regression function and m))).) for the fitted regression function for the şth Monte Carlo run. The MSPE is computed as MSPE, 100 Ei d2 [ [m, (X(st), m(x ( where [x(es)] 100 denote out-of-sample predictors and d is the metric for the corresponding metric space. The average performance over 500 Monte Carlo runs is g= MSPE3-
quantified by AMSPE-
Σ 500
500
Distributions. We consider truncated one-dimensional Gaussian distributions with random pa- rameters 77 and a on 2-2, 2] as distributional outputs, characterized by quantile functions Q(p) X
E( ) - ((((a) + + - +(a))), where Φ(-) represents the cumulative distribution
function of the standard Gaussian distribution, a -2-E(n)X) and b 2 2-E(r) To generate E(e X)
B(o
distributional outputs, the predictor X E R9 is sampled as follows:
X, ~ U0, 1), X2 ~ U(-1,1), Xa U(-2, 2), X4 ~ N(0, 1), X5 ~ 1),,1), X6 ~ (O, 1), X, ~ Ber(O.1), Xg ~ Ber(0.2), Xg Ber(0.5).
7

---

Mean 17 and standard deviation a of the truncated Gaussian distribution that serves as the distributional output are generated conditional on predictor vectors X, where
x| NN,,0.55) a X Gamma , 8-1) and μ = sin( -
1 + 2 cos(( X2/2) X
To mimic real-world scenarios where direct access probability distributions is unavailable, simulate independent data samples for each distributional output. Specifically, 100 observations {Mr) 1100 are sampled independently from each distribution Y,. Consequently, one must initially estimate the distributional output Y; from the random sample {yt )=ر introducing a bias in the regression model. This setup reflects practical challenges and aligns with prior approaches [65]. where the empirical measure is adopted as a proxy for the latent distribution Yi.
Networks. Consider simple, undirected, weighted networks with fixed number of nodes and bounded edge weights. Such networks can be uniquely characterized using their graph Laplacians 64.. which are symmetric matrices satísfying specific constraints. Formally, the space of graph Laplacians is
M (Y (vy): = y Y11, O₃; there exists W >> 0 such that W Vij for j),
where 1, and O, are l-vectors of ones and zeroes, respectively.
To construct generative model for network outputs, we draw predictors E B from the following distributions:
X1 U(-1, 1), X2 ~ U(-1,1), X3 ~ (L, X1 Garmaa33 1), X5 Gamma (4, 1), X6 Gamma(5, 1), X7 Ber(0.2), Xs Ber(0.3), X9 Ber(0.5).
For the corresponding network output, the weights of the edges are modeled using beta distribution with shape parameters 2X7 sin ( + -1-44 (1- cos"(T and B = XīXs + X5 (1- Xs). The generated edge weights are then used construct the graph Laplacian, which serves as the output. As an alternative to the proposed geometry- aware approach, one could apply XGBoost [13] to vectorized representations of the graph Laplacians. We evaluate this baseline approach in Appendix and find that it underperforms FGBoost, highlighting the importance of respecting the geometric structure of the output space.
Discussion on the simulation results. Table presents the AMSPE for FGBoost and the competing models. As the sample size increases, FGBoost exhibits a clear trend of decreasing prediction error, indicating its convergence to the target. Across all data types and sample sizes considered, FGBoost consistently outperforms the competing methods. This performance gap becomes increasingly pronounced with larger sample sizes, underscoring the scalability and accuracy of FGBoost in capturing complex relationships between multivariate predictors and non-Fuclidean outputs.
6 Real-world data applications
We evaluate the performance of the FGBoost algorithm using real-world datasets, including human mortality data with distributional outputs and New York City yellow taxi data with network outputs. An additional application involving compositional data is presented in the appendix. To assess predictor importance, developed adapted version of Shapley Additive Explanations (SHAP) values 37), with further details provided in the appendix.
6.1 Human mortality data
This analysis uses age-at- at-death distributions from 162 countries in 2015 as distributional outputs.
The life tables, sourced from the United Nations World Population Prospects 2024 (https: /// population ww..uun.rrr..................... wtg/oorgpppp hum.org/sppp//dmloaae)).. downloade)) provide death counts grouped into five-year age intervals. Using the frechet package [14), these aggregated death counts were smoothed via a local linear smoother and then normalized through trapezoidal integration to generate density estimates. Age-at- death distributions are influenced by variety of factors, and consider nine-dimensional predictor

---

Table Average mean squared prediction errois and standard deviations (in parentheses) of FGBoost, global Fréchet regression (GFR) [ sufficient dimension reduction (SDR) [62.. single index Fréchet regression (IFR) 4 Fréchet random forest (FRF) [00] and random forest weighted local
linear Fréchet regression (RFWILFR) 49 for simulated distributional and network outputs.
Output
FGBoost
GFR
SDR
IFR
FRF
RFWLLFR
100
0.034
0,053
0.098
0.041
0.048
0.097
(0.011)
(0.014)
(0.043)
(0.012)
(0.013)
(0.042)
200
0.028
0.043
0.059
0.038
0.041
0.075
(0.010)
(0.009)
(0.022)
(0.008)
(0.010)
(0.021)
Distribution
500
0.023
0.038
0.040
0.037
0.034
0.057
(0.009)
(0.008)
(0.013)
(0.008)
(0.008)
(0.013)
1000
0.019
0.037
0.035
0.037
0.029
0.046
(0.007)
(0.007)
(0.012)
(0.008)
(0.006)
(0.009)
2000
0.015
0.036
0.026
0.038
(0.005)
(0.006)
(0.005)
(0.006)
100
13.644
15.326
19.448
17.768
13.820
19.321
(3.140)
(2.570)
(11.786)
(5.529)
(2.743)
(4.157)
200
10.531
14.309
14.391
16.619
12.190
16.528
(3.371)
(2.474)
(2.862)
(3.325)
(2.372)
(3.113)
Network
500
6,912
13.831
12,473
16.382
10.659
14.572
(1.950)
(2,591)
(2.351)
(3,343)
(2.205)
(2.769)
1000
5,471
13.481
11.769
16.100
9.376
13.066
(1.481)
(2.383)
(1.865)
(2.979)
(1.957)
(2.673)
2000
4.996
13.459
8.308
11.891
(1.325)
(2.329)
(1.765)
(2.450)
Table 2: Average mean squared prediction errors and standard deviations (in parentheses) of FGBoost,
global Fréchet regression (GFR) [47), sufficient dimension reduction (SDR) (62), single index Fréchet regression (IFR) 4), Fréchet random forest (FRF) [10] and random forest weighted local
linear Fréchet regression (RFWLLFR) [49 for human mortality and taxi network data.
Data
FGBoost
GFR
SDR
IFR
FRF
REWLLFR
Human Mortality
20.35
31.41
27.60
58.36
22.02
22.55
(36.63)
(58.83)
(44.40)
(107.54)
(38.08)
(35.49)
Taxi Network
8.30
11.80
13.15
32.16
9.50
10.17
(0.18)
(0.20)
(0.38)
(2.14)
(0.15)
(0.68)
set encompassing demographic, economic, and environmental variables. Detailed deseriptions of these predictors are provided in Table in the appendix.
To evaluate model performance, leave-one-out cross-validation is employed compute the MSPE, with results summarized in Table [ FGBoost achieves the lowest MSPE, demonstrating its effec- tiveness in capturing complex distributional relationships, even with relatively small sample sizes. Figure
2 presents the SHAP summary plot, ranking predictors by their importance in FGBoost. The analysis reveals that GDP has the most substantial impact on age-at-death distributions, corroborating prior studies emphasizing the critical role of socioeconomic factors in health outcomes 43 42 Other influential factors include mean childbearing age, health expenditure, and population density.
6.2
New York City yellow taxi data
We analyze transport dynamics in Manhattan, New York, using yellow taxi trip records obtained from hattps://wwt nyc. www.gggo////sseee
reppr-daorddaaa page Manhattan is partitioned into 13 regions, for which daily transport networks are constructed. The nodes of the networks represent the 13 regions and their edge weights indicate the number of passengers traveling between them. We characterize these networks by 13 13 graph Laplacian matrices and associate these with 1 12-dimensional predictor vectors, with components including weather attributes and weekday/holiday indicators, as detailed Table [T0 the appendix. We examine how these factors impact the transport networks, based on data spanning January 1, 2017, December 31, 2018.

---

GDP
Hear Shearing Age-
Fdssy and Soturdoy
Torperaure
Feature
Talls Amoun
CPI-
Sox Rado
GVA
Pressure
Widd
Predolitóon
5053
1000
1909)
meants HAP HP value)
meants HAP value)
Figure 2: Summary plot of SHAP values for FGBoost applied to human mortality data (left) and taxi network data (right). Features are sorted by their impact in descending order.
Model performance is evaluated through five-fold cross-validation. with the MSPE averaged over 100 runs, as reported in Table FGBoost achieves the lowest MSPE, consistently outperforming all competing models in predicting transport networks. Figure Sesure displays the SHAP summary plot, ranking predictors by their influence on the predictions obtained by FGBoost. Key factors include passenger count, tip amount, and the indicator for whether is Friday/Saturday, with temperature and toll amounts also contributing, though to a lesser extent.
Conclusion
We propose FGBoost, innovative intrinsic regression method designed for geodesic metric spaces that successfully addresses the challenge of the absence of linear operations, which are essential for traditional gradient boosting. By leveraging geodesics proxies for residuals, FGBoost iteratively builds ensembles while preserving the geometric properties the output space, ensuring intrinsic compatibility with the underlying metric space structure. FGBoost is supported by theoretical results and performance is highly competitive in numerical experiments and real data analysis.
Future research could focus extending the theory of FGBoost beyond Hadamard spaces. Another extension of interest for future research will be to extend FGBoost to handle scenarios where both predictors and outputs reside in general metric spaces. One potential solution involves modifying the splitting rules for tree-based FGBoost to operate intrinsically within metric spaces [10]. This would enable FGBoost to address regression tasks with non-Euclidean predictors.
Acknowledgments and Disclosure of Funding
would like to thank the reviewers their constructive feedback. This research was partially supported by NSF grant DMS-2310450.
References
[1] Miroslav Bacák. Computing medians and means in Hadamard spaces. SIAM Journal on Optimization, 24(3):1542-1566 2014.
[2] Alexis Bellot and Mihacla van der Schaar. Multitask boosting for survival analysis with competing risks. In Advances Neural Information Processing Systems, volume 31, 2018.
[3] Alina Beygelzimer, Elad Hazan, Satyen Kale, and Haipeng Luo. Online gradient boosting. In Advances in Neural Information Processing Systems, volume 28, 2015.
[4] Satarupa Bhattacharjee and Hans-Georg Müller. Single index Fréchet regression. Annals of Statistics, 51(4):1770-1798 2023.
[5] Louis J. Billera, Susan P. Holmes, and Karen Vogtmann. Geometry of the Space of Phylogenetic Trees. Advances in Applied Mathematics, 27:733-767. 2001.
[6] Leo Breiman. Arcing classifier (with discussion and rejoinder by the author). Annals of Statistics, 26(3):801-849, 1998.
[7] Martin R Bridson and André Haefliger Metric Spaces of Non-Positive Curvature. Grundlehren der mathematischen Wissenschaften. Springer Berlin, Heidelberg, 1999.
10

---

[8] Sarah Brockhaus, Michael Melcher, Friedrich Leisch, and Sonja Greven. Boosting flexible functional regression models with high number of functional historical effects. Statistics and Computing, 27:913-926, 2017.
[9] Michael M Bronstein, Joan Bruna, Yann LeCun, Arthur Szlam, and Pierre Vandergheynst. Geometric deep learning: going beyond Euclidean data. IEEE Signal Processing Magazine, 34(4):18-44 2017.
[10] Louis Capitaine, Jérémie Bigot, Rodolphe Thiébaut, and Robin Genuer. Fréchet random forests for metric space valued regression with non Euclidean predictors. Journal of Machine Learning Research, 25(355): 1-11, 2024.
[11] Hsin-wen Chang and Ian McKeague. Empirical likelihood-based inference for functional means with application to wearable device data. Journal of the Royal Statistical Society Series B: Statistical Methodology, 84(5):1947-1963 2022.
[12] Shang-Tse Chen, Hsuan-Tien Lin, and Chi-Jen Lu. An online boosting algorithm with theoretical justifications. In International Conference on Machine Learning. 2012.
[13] Tianqi Chen and Carlos Guestrin. XGBoost: scalable tree boosting system. ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 785-794, 2016.
[14] Yaqing Chen, Yidong Zhou, Han Chen, Alvaro Gajardo, Jianing Fan, Q Zhong, P Dubey, Kyunghee Han, S Bhattacharjee, Changbo Zhu, Su I lao, Poorbita Kundu, Alexander Petersen, and Hans-Georg Müller. frechet: Statistical Analysis for Random Objects and Non-Euclidean Data. R package version 0.3.0, 2023.
[15] lan L. Dryden, Alexey Koloydenko, and Diwei Zhou. Non-Euclidean statistics for covariance matrices, with applications diffusion tensor imaging. Annals of Applied Statistics, 102- 1 1123, 2009.
[16] Paromita Dubey, Yaqing Chen, and Hans-Georg Müller. Metric statistics: Exploration and inference for random objects with distance profiles. Annals of Statistics, 52(2):757-792, 2024.
[17] Julian J Faraway. Regression for non-Euclidean data using distance matrices. Journal of Applied Statistics, 41:2342-2357, 2014.
[18] Frédéric Ferraty and Philippe Vieu. Additive prediction and boosting for functional data. Computational Statistics & Data Analysis, 5304):1400-1413, 2009.
[19] Jerome H. Friedman. Greedy function approximation: A gradient boosting machine. Annals of Statistics, 1899-122,, 2001.
[20] Aritra Ghosal, Wendy Meiring, and Alexander Petersen. Fréchet single index models for object response regression. Electronic Journal of Statistics, (7(1):1074-111) 2023.
[21] Trevor Hastie, Robert Tibshirani, and Jerome Friedman. The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer New York, NY, 2nd edition, 2009.
[22] Matthias Hein. Robust nonparametric regression with metric-space valued output. In Advances in Neural Information Processing Systems, pages 718-726, 2009.
[23] Torsten Hothorn, Peter Bühlmann, Sandrine Dudoit, Annette Molinaro, and Mark J Van Der Laan. Survival ensembles. Biostatistics, 7(3):355-37,, 2006.
[24] Su I Iao and Hans-Georg Müller. Measure selection for functional linear model. Computational Statistics & Data Analysis, page 108270, 2025.
[25] Su lao, Yidong Zhou, and Hans-Georg Müller. Deep fréchet regression. Journal of the American Statistical Association, pages 1-30, 2025. in press.
[26] Leonid losipoi and Anton Vakhrushev. Sketchboost: Fast gradient boosted decision tree for multioutput problems, Advances Neural Information Processing Systems, 35:25422-25435. 2022.

---

[27] Guolin Ke, Qi Meng, Thomas Finley, Taifeng Wang, Wei Chen, Weidong Ma, Qiwei Ye, and Tie- Yan Liu, LightGBM: A highly efficient gradient boosting decision tree. In Advances in Neural Information Processing Systems, volume 30, 2017.
[28] Eric D. Kolaczyk and Gábor Csárdi. Statistical Analysis of Network Data with R, volume 65. Springer Cham, 2nd edition, 2020.
[29] Alan B. Krueger, Andreas Mueller, Steven J. Davis, and Aysegüi Şahin. Job search, emotional well-being, and job finding in a period of mass unemployment: Evidence from high frequency longitudinal data [with comments and discussion]. Brookings Papers on Economic Activity, pages 1-81, 2011.
[30] Daisuke Kurisu, Yidong Zhou, Taisuke Otsu, and Hans-Georg Müller. Geodesic causal inference. arXiv preprint arXiv:2406. 19604, 2024.
[31] Mark K Ledbetter, Lucia Tabacu, Andrew Leroux, Ciprian M Crainiceanu, and Ekaterina Smirnova. Cardiovascular mortality risk prediction using objectively measured physical activity phenotypes in NHANES 2003-2006. Preventive Medicine, 164:107303, 2022.
[32] Donald KK Lee, Ningyuan Chen, and Hemant Ishwaran. Boosted nonparametric hazards with time-dependent covariates. Annals of Statistics, 49(4):2101-222) 2021.
[33] Zhenhua Lin. Riemannian geometry of symmetric positive definite matrices via Cholesky decomposition. SIAM Journal on Matrix Analysis and Applications, 40(4): 2019.
[34] Zhenhua Lin, Dehan Kong, and Linbo Wang. Causal inference on distribution functions. Journal of the Royal Statistical Society Series B: Statistical Methodology, 85(2):378-39) 2023.
[35] Zhenhua Lin and Fang Yao. Intrinsic Riemannian functional data analysis. Annals of Statistics, 47(6):3333-3577, 2019.
[36] Scott M. Lundberg, Gabriel Erion, Hugh Chen, Alex DeGrave, Jordan M. Prutkin, Bala Nair, Ronit Katz, Jonathan Himmelfarb, Nisha Bansal, and Su-In Lee. From local explanations to global understanding with explainable Al for trees. Nature Machine Intelligence, 2(1):56-67, 2020.
[37] Scott M. Lundberg and Su-In Lee. A unified approach to interpreting model predictions. In Advances in Neural Information Processing Systems, pages 4765 4774, 2017.
[38] Llew Mason, Jonathan Baxter, Peter Bartlett, and Marcus Frean. Boosting algorithms as gradient descent. In Advances in Neural Information Processing Systems, volume 12, 1999.
[39] Marcos Matabuena and Alexander Petersen. Distributional data analysis of accelerometer data from the NHANES database using nonparametric survey regression models. Journal of the Royal Statistical Society Series C: Applied Statistics, 72(2):294-313, 2023.
[40] Takuo Matsubara. Wasserstein gradient boosting: A framework for distribution valued super- vised learning. In Advances in Neural Information Processing Systems, 2024.
[41] Robert J. McCann. A convexity principle for interacting gases. Advances in Mathematics, 128(1)1153-179 1997.
[42] Goran Miladinov. Socioeconomic development and life expectancy relationship: Evidence from the EU accession candidate countries. Genus, 76(1):2, 2020.
[43] Abdalali Monsef and Abolfazl Shahmohammadi Mehrjardi. Determinants of life expectancy: A panel data approach. Asian Economic and Financial Review, 5(11): 2015.
[44] Tom MW Nye, Xiaoxian Tang. Grady Weyenberg, and Ruriko Yoshida. Principal component analysis and the locus of the Fréchet mean in the space of phylogenetic trees. Biometrika, 104(4)9911-22, 2017.
[45] Victor M. Panaretos and Yoav Zemel. An Invitation to Statistics in Wasserstein Space. Springer New York, 2020.
12

---

[46] Xavier Pennec, Pierre Fillard, and Nicholas Ayache. A Riemannian framework for tensor computing. International Journal of Computer Vision, 66(1):41-66, 2006.
[47] Alexander Petersen and Hans-Georg Müller. Fréchet regression for random objects with Euclidean predictors. Annals of Statistics, 47(2):691-71,, 2019.
[48] Alexander Petersen, Chao Zhang, and Piotr Kokoszka. Modeling probability density functions as data objects. Econometrics and Statistics, 21:159-178, 2022.
[49] Rui Qiu, Zhou Yu, and Ruoqing Zhu, Random forest weighted local Fréchet regression with random objects. Journal of Machine Learning Research, 25(107):1-69, 2024.
[50] Janice L. Scealy and Alan H. Welsh. Regression for compositional data by using distributions defined on the hypersphere. Journal of the Royal Statistical Society Series B: Statistical Methodology, 73(3):351-333 2011.
[51] Janice L. Scealy and Alan H. Welsh. Colours and cocktails: Compositional data analysis 2013 lancaster lecture. Australian & New Zealand Journal of Statistics, 56(2):145-169, 2014.
[52] Robert E Schapire and Yoav Freund. Boosting: Foundations and Algorithms. MIT Press, 2012.
[53] Katie E. Severn, lan I.. Dryden, and Simon P. Preston. Manifold valued data analysis of samples of networks, with applications in corpus linguistics. Annals of Applied Statistics, 16(1):368-30) 2022.
[54] Dogyoon Song and Kyunghee Han. covariate approximation. In Thirty-seventh Conference on Neural Information Processing Systems, 2023.
Fréchet regression with low-rank
[55] Erik Strumbelj and Igor Kononenko. Explaining prediction models and individual predictions with feature contributions. Knowledge and Information Systems, 6166476665, 2014.
[56) Karl-Theodor Sturm. Probability measures on metric spaces of nonpositive curvature. Heat Kernels and Analysis Manifolds, Graphs, and Metric Spaces (Paris, 2002). Contemp. Math., 338. Amer. Math. Soc., Providence, RI, 338:357-390, 2003.
[57] Gerhard Tutz and Jan Gertheiss. Feature extraction in signal regression: a boosting technique for functional data regression. Journal of Computational and Graphical Statistics, 19(1):154-174 2010.
[58] Aad Van der Vaart and John Wellner. Weak Convergence and Empirical Processes: With Applications Statistics. Springer, New York, 2nd edition, 2023.
[59] Jane-Ling Wang, Jeng-Min Chiou, and Hans-Georg Müller. Functional Data Analysis. Annual Review of Statistics and its Application, 3:257-295, 2016.
[60] Chao Ying and Zhou Yu. Fréchet sufficient dimension reduction for random objects. Biometrika. 109(4):975-992, 2022.
[61] Ying Yuan, Hongtu Zhu, Weili Lin, and JS Marron. Local polynomial regression for sym- metric positive definite matrices. Journal of the Royal Statistical Society Series B: Statistical Methodology, 74:697-719, 2012.
[62] Qi Zhang, Lingzhou Xue, and Bing Li. Dimension reduction for Fréchet regression. Journal of the American Statistical Association, 119(548):2733 2747, 2024.
[63] Tong Zhang and Bin Yu. Boosting with early stopping: Convergence and consistency. Annals of Statistics, 33(4):1538-117) 1538-1579, 2005.
[64] Yidong Zhou and Hans-Georg Müller. Network regression with graph Laplacians. Journal of Machine Learning Research, 23:1-41, 2022.
[65] Yidong Zhou and Hans-Georg Müller. Wasserstein regression with empirical measures and density estimation for sparse data. Biometries, 80(4) 2022 2024.
13

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The contributions and scope the paper are outlined in the abstract and further detailed in the introduction. Specifically, the scope is described in the main body of the introduction, and the key contributions are listed Section 1.1.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results be expected to generalize other settings.
is fine include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: limitation section included Appendix K due to space constraints the main text.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and prool's
Question: For each theoretical result, does the paper provide the full of assumptions and complete (and correct) proof?
14

---

Answer: [Yes]
Justification: The assumptions for the theoretical results are stated in Sections and 4, and the complete proofs provided in the appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: Detailed experimental setups and results are provided in Section 5, The supplementary zip file contains code for data simulation with fixed random seeds, real data files with descriptive names, and README with full usage instructions.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, No answer this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data provided or not.
the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how replicate the results, access a hosted model (e.g.. in the of large language model), releasing of a model checkpoint, or other means that are appropriate the research performed.
While NeurIPS does not require releasing code, the conference does require submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing verifying the results.
15

---

5. Open access to data and code
Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in the supplemental material?
Question:
Answer: [Yes]
Justification: The supplementary materials provide open access to both data and code. Reproduction of the main results supported with detailed instructions, including fixed random seeds for training. The data application is organized in clearly labeled folders, with implementation guidance outlined in Section 6 and further detailed in the README.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (httpss:/miis.cce)
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
/ //nips. httpppsccppbblccgguiiieeeeeepppp /obblli/ggiddss//deeuubbiiii (roodeSSumsssiooniii ssiliii for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e-g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Key training settings described in Section with detailed information provided the README file.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: The standard error of the mean square prediction error estimation is provided in all tables.
Guidelines:
The answer NA means that the paper does not include experiments.
16

---

The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (c.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification:
We provide the information on compute resources in Appendix
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics httppp///.aarrss.ec/aaaaaa/s/aae
Answer: [Yes]
Justification: The research conducted in the paper conforms with the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
17

---

This paper proposes a new boosting method and focuses on theoretical, simulation results, and application results for regression with metric space-valued outputs. Given the scope and focus of the paper, the discussion of negative societal impacts is deemed irrelevant and therefore not included.
Justification:
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative socictal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deploymentt of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: There is no risk of misuse of the proposed method and the datasets used in the paper are open-sourced.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (c.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: We have cited the original paper or attached the link to the existing assets used in this paper.
Guidelines:
18

---

The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g.. CC-BY 4.0) should be included for each asset. For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: The paper does not release any new assets.
Guidelines:
The answer NA means that the paper does release assets.
Researchers should communicate the details the dataset/code/mooel as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and sereenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve experiments or research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included in the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
The paper does not involve crowdsoureing nor research with human subjects.
19

---

Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with
human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods this research? Note that the LLM used only for writing, editing, or formatting purposes and does not impact the core methodology. scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
Justification: The paper's 's core method development does not involve LLMs for any impor- tant, original, or non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, or non standardd components.
Please refer to our LLM policy http::::/eerrrppsscc/aoeeeeerree for what should or should not be described.
20

---

A Compositional data
In this section, we explore another important example: the space of compositional data. We demon- strate the applicability and effectiveness of FGBoost for this type of data through numerical simula- tions and a real-world case study.
Example 4 (Space of compositional data). Compositional data, represented proportions summing to I. reside in the simplex: A { R yj 2 0, j L....,d. d, and Σj=1 V, = 1}. Using the square-root transformation /y = /yh,,...... yya)² the simplex can be mapped the positive orthant of the unit sphere (50 5TT: S = fz E S 2j 2 0, j ......d.. Equipping S with the geodesic (Riemannian) metric the sphere, dz(21, 22) = arccos [T] x2) for 21, 22 S4
induces a unique geodesic structure. The geodesic connecting two points 21, 22 E S is explicitly defined as:
T=L.82 (t) = cos()))== + sin 22 (2 22)
t [0
22 (2|
where 0 arccos(=fz2) 22) is the angle between 2 and 22.
A.1 Experiments for compositional data
Consider three-dimensional compositional data residing on S7, the positive segment of the unit sphere in R², equipped with the geodesic metric. The Euclidean predictor X E R includes X1...., Xy, distributed identically the predictors in the network simulation, and an additional variable X ~ Ber(O.1).
The true regression function is modeled m( (1 Xs)mo((X Xgm (X), where mi (X) (cos(e), sin( v/3 (in))))) and á π (f) E [=/8, 37/ The function f(X) ) is defined fs b(X)) (la(((((X)))))) [b(x)]} where (X)) = 3.X sin s (= 3(1 Xi)) cos" (x and b(X) -X X1 XX5
mo (X) (cos(o),
To generate the random output Y on S for Xg 0, we add small perturbation the true regression function mo! X). First, we construct an orthonormal basis (e1, e2) for the tangent space at mo (X) where
e1 = (sin()), 3 (Ø)) - ccos()/2), C2 = (0, 1/2, -3/22).
Next, we consider random tangent vectors U = Z1e1 + Z2e2, where Z1, Zz2 are independent random variables uniformly distributed on [-0.1, 0.1]. The random output Y is then obtained by applying the exponential map at (X) to U,
U
Y = (x) (U) cos(|IU]))mo(X) + sin((
For X3- = 1, a similar procedure is followed with the orthonormal basis (e1, e2) for the tangent space
m1 X) defined as
c = (sin(Ø), - cos(0)/2 - - coc) e2 = (0, 33/2,-1/21
Figure illustrates randomly generated outputs using the above generation procedure for sample size n
500.
FGBoost consistently outperforms competing methods across range of sample sizes, shown in Table The sole exception occurs at a smaller sample size n = 100, where GFR and FRF exhibit slightly better performance. However, the sample size grows, FGBoost not only catches up but also showcases a marked improvement, further solidifying its superiority in handling larger datasets.
A.2 Emotional well-being for unemployed workers
A survey of unemployed workers in New Jersey 29 was conducted during the fall of 2009 and early 2010, a period when the U.S. unemployment rate peaked 10% following the 2007-2008 financial crisis. The analysis includes data for 3301 unemployed workers with complete measurements.
21

---

Table Average mean squared prediction errors and standard deviations (in parentheses) of Fréchet geodesic boosting (FGBoost), global Fréchet regression (GFR) [47]. sufficient dimension reduction (SDR) 62,, single index Fréchet regression (IFR) [4, Fréchet random forest (FRF) [00] and random
forest weighted local linear Fréchet regression (RFWLLFF)) 49 for compositional outputs.
Output
n
FGBoost
GFR
SDR
IFR
FRF
RFWLLER
100
0.0104
0.0099
0.0382
0.0800
0.0088
0.0456
(0.0041)
(0.0012)
(0.0068)
(0.0082)
(0.0012)
(0.0639)
200
0.0074
0.0089
0.0368
0.0782
0.0076
0.0243
(0.0018)
(0.0009)
(0.0058)
(0.0078)
(0.0011)
(0.0302)
Compositional
500
0.0047
0.0085
0.0364
0.0780
0.0064
0.0112
(0.0008)
(0.0008)
(0.0057)
(0.0073)
(0.0009)
(0.0017)
1000
0.0038
0.0083
0.0361
0.0770
0.0056
0.0092
(0.0007)
(0.0008)
(0.0054)
(0.0077)
(0.0009)
(0.0014)
XB X7 -0-00
XB 1,
1,
0.2
0.4
Figure 3: Visualization of simulated compositional data for 500.
The key response variable is the proportion of time spent in each of the four moods while at home: bad, low/irritable, mildly pleasant, and very good. The compositional response vector is represented as y (y1, y2, Us, y4)* where y/j denotes the proportion of time spent in the jth mood (j 1,..., 4). Applying square-root transformation, z (21, 22, 23, 24)= (VTi, V /Bs. VJ4) maps the outputs to the positive orthant of the unit sphere S The corresponding geodesic metric d(21, #2) arccos(2] z2) is used for the analysis. The predictors for this application consist of 10 bascline socio-cconomic and demographic variables collected through the questionnaire: (1) life satisfaction (2) highest education level (3) marital status (4) the number of children, (5) the number of people in the household,
(6) total annual household income, (7) hours per week working at the last job, (8) how the last job ended, (9) weeks spent looking for work, and (10) credit card balance.
Model performance is assessed using ten-fold cross-valle validatii,,,,,, with the MSPE averaged over 100 runs, as reported in Table FGBoost demonstrates substantial improvement, achieving more than a 50% reduction MSPE compared to GFR, SDR, and IFR. FRF and RFWLLFR are not included in the comparison as their implementations only work for three-dimensional compositional data. Figure | presents the SHAP summary plot, ranking predictors by their importance. The analysis identifies life satisfaction the most influential predictor, followed by credit balance, weeks spent job secking, and household size, though their impacts are considerably smaller than that life satisfaction.
B Geodesic transport maps
In the Wasserstein space from Example I Assumption is satisfied with the geodesic transport the cumulative distribution function of and map Ty,.8 F Fa., where F, and Fø
22

---

Table Average mean squared prediction errors and standard deviations (in parentheses) of Fréchet geodesic boosting (FGBoost), global Fréchet regression (GFR) [47]. sufficient dimension reduction
(SDR) 622, and single index Fréchet regression (IFR) for emotional well-being data.
FGBoost
GFR
SDR
IFR
0.2074 (0.0005) 0.4163 (0.0007) 124112 (0.0007)
0.4356 (0.0015)
Life Satisacction
Card Balance
Weeks Searching
Num of People [Househokd]
Feature
Education
Houns/wnek-
Marita slatus
hildren
Last ob E nd-
0.00
0.05
0.10
0.15
0.20
mean(S value)
Figure 4: Summary plot of SHAP values for FGBoost applied to emotional well-being data. Features are sorted by their impact descending order.
the quantile function of B, respectively. The resulting endpoint of the geodesic Ya, is given by Fr' - F3 o Fa F51, where F-1 and P-- denote the quantile functions of w and Ç, respectively,
For the space of networks or symmetric positive-definite matrices equipped with the Frobenius metric, the geodesic corresponds to the line segment connecting the starting and ending points. Assumption | satisfied for Examples and 3 with the geodesic transport map T-... - (w) = w + (B ax).
For the space of compositional data described in Example 4 the geodesic transport map can be interpreted as a rotation of the point along the geodesic determined by a and B. The tangent vector for the geodesic Too is given by ta,3 3 (a' 3)a, whose magnitude and direction encode the geodesic length and directionality needed to move from a toward along the sphere. The geodesic transport map is then defined as:
T T (w) = ExP (0
cos( + sin(()
v
|e||
|
where 0 arccos(a 33 is the angle between and Va, (w' g y is the projection of 2a onto the tangent space at w, and (Ou/,(0)/||lll) denotes the exponential map at w applied to the tangent vector 0c/l|ui/.
This map T (w) moves the point w along the geodesic connecting it to a new point determined by and B, with the direction and distance dictated by the original geodesic "a The construction ensures that T (w) lies on the sphere and preserves the geodesic structure. Assumption [ is satisfied with this geodesic transport map.
Beyond these examples, extending FGBoost to a new geodesic space requires specifying only two ingredients: the distance function and the associated transport map. For smooth Riemannian manifolds, geodesics and parallel transport are classical and well-studied, with closed-for or numerically stable algorithms widely available. For discrete structures such as trees networks, geodesics correspond to shortest paths under the chosen metric, and transport can be defined by propagating along these paths. For specialized metrics, such the BHV metric for phylogenetic trees 5. geodesics and transport maps are explicitly described in the literature.
The transport map assumption, therefore, does not pose major barrier in practice. For empirical distributions, transport maps are based on empirical quantile functions, which are simple step functions and computationally straightforward. For Riemannian manifolds, parallel transport is standard and efficiently implemented. These principles make FGBoost broadly applicable and provide practitioners with clear guidelines for adapting the method to new domains.
23

---

C Shapley Additive Explanations for Fréchet geodesic boosting
Tree-based regression models, such boosted trees and random forests, widely used for their flexibility and ability to model complex non-linear relationships. However, explaining their predic- tions often receives less attention. Shapley values (55 377 provide principled way to measure feature importance for predictive models. Shapley values require retraining the model on all subsets of features _ San, where Sa, is the of all features. They assign importance value each feature based on its effect on the model's prediction. To quantify this effect for feature j. two models considered: [su(i): trained with feature included, and fs. trained without feature j. The contribution of feature; ÿ is measured as the difference in predictions between the two models on the same input, {suns (T)))))) (2ss, where æs represents the values of the input features in the subset S. Since the effect of withholding a feature depends on interactions with other features, the differences are computed for all subsets S S Sw {i}. The Shapley values are then computed as weighted average of these differences:
S| (1
[S| 1)! 18
(,( = Σ
fs(æs)]
SCS
To extend Shapley values non-Buclidean outputs, the difference term is replaced by the metrie:
[S - -1| 1)!
05({,z)=
Σ SC5.ND
d(Jsut) (oss))
SSal
where d is the metric of the output space. Because d(..-) 0, these Shapley values quantify the magnitude of each feature's effect and are non-negative by construction.
Shapley Additive Explanations (SHAP) values, introduced by 37 generalize Shapley values to quantify feature contributions for the conditional expectation function of the model's output. Tree SHAP [3). a specialized algorithm for tree-based models, exploits the hierarchical structure of decision trees to efficiently compute SHAP values without retraining the model and evaluating all possible subsets. SHAP value enables both global and local interpretability. Globally, the mean absolute SHAP values highlight the importance of each feature across the dataset, providing insights into which predictors drive the model's 's behavior. Locally, SHAP values for individual predictions explain how specific feature values contribute to the output. By combining accuracy and transparency, SHAP enhances the interpretability of complex models and builds trust their predictions. Such an extension can also be useful for other regression models for non-Euclidean outputs.
D Proof of Proposition
1
Proof. To prove that dg is valid metric the space of geodesics (MM), verify the following four axioms:
Identity: For any Ta.3 E G(M),
do(T0 Ta.3) + ((, =0.
Positivity: For any 8, Ta G(M), + 02 or + 82, then Ja, + Taa..2 In this case, at least one of (a 02) or dB B2) is strictly greater than O, implying
dg dg(Yan Bn Tog.3,) = / (2 d3 0.
Symmetry: For any Ton To2.82 E G(M).
dg d 7o. Vdl(a,
(((11,8)) ( dg(Tan, 711,8))
24

---

Triangle inequality: For any 0: Pas.8 9, E (.M), the following holds:
de (7an,8: (a,,23) kg) = (3),
((02, (3)
= +( 0z) d
+ [F
( + de d (32, B3)]
=( (v)(ar, (11. + (a2, 0z) - dp (32. B3)
= =dd .1 Tog.9g) + dg (Tog.82: Tag.8a)) }3
Taking the square root on both sides, the triangle inequality follows:
da .1 Tog.82) I dg 81 Tos.8s) + dg ( (72 Yog.9g)
Since all four axioms are satisfied, dg is valid metric on the space of geodesics (
E Proof of Proposition 2
Since (M, d) is a Hadamard space, the space of geodesic (G(M), dg), as product metric space, also Hadamard space. From Proposition 2.3 (56), for any pair of geodesics %o. Vi there exists a unique geodesic r [0, 1] > G(M) connecting them, and the intermediate points Yt = (t),, & [0, depend continuously on the endpoints To: 71. Furthermore, according to the definition of Hadamard space, for any 7 ((M)
((M).
(1) $ (1 t) (7, 70) + to (7: 7) - t(1 - ()) 71)
This inequality demonstrates that the function ( :) is strongly convex over g(M).
Next, we prove that v'(7,) :) is Lipschitz continuous. Let 71, 72 E g(M) be two geodesics, it follows that
(1()))) - (11)))-((700000 72)1 (6(()))))66m
do(7) dg 72)
-dg(7)))))))))
22diam(9(M)) do(r) do(7)
(2diam(9(M)) dg dg(71,72)
where diam(G(M)) denotes the diameter of (MM) and is finite since M is a bounded metric space. Thus ((x, :) is Lipschitz continuous with respect dg.
F Proof of Theorem
According to Proposition 2 the risk functional A(-) is strongly convex and continuous. Furthermore, A(-) is bounded as M is bounded. Thus, in particular,
A(F) = inf A(F), FF)
inf
FEspan[.")
where span(F) is the closure of span(F). Since A(.) is strongly convex, there exists unique function F* E sFan(F) such that
F
= arg min A(F).
FE)pan(T)
The uniqueness follows directly from the strong convexity of the loss function. Similar arguments apply to the empirical risk functional A-).)
25

---

G
Proof of Theorem
Proof. Define the metric on span(F) as:
dx (Fi, P2)
sup (F F2(2:))
=eRr
is straightforward to verify that dy is valid metric. Specifically, for any F1, F2, Fs E span(F),
sup do [i (ar), F2(2:) v sup (de(Fi (s)) do (F) Fa())))
mERY
#E3*
sup dg Fa(=)) sup (2)))))) Fa())). zeRr
Thus, (Span(F): dr) forms a metric space.
Let 1 (span (F) represent the space of bounded functions on span(F), To establish that in probability, it suffices to show that A.(+) A(-) weakly (spann (J)). Once this weak convergence is shown, Theorem 1.3.6 of [58] can be applied to conclude the result. For detailed definition of weak convergence in this context, we refer readers Definition 1.3.3 in [58]. By Theorems 1.5.4 and 1.5.7 of (58), the weak convergence follows upon verifying the following two conditions:
suprespan(F) (F) An (F) A(F)
converges
(i) AF. A(F) (1) for all F E span(F) and (ii) A.((-) A()) is asymptotically equicontinuous in probability, i.e., for all e, > there exists 8 > O such that
lim sup P(
sup (a)
(An(F)) A(F)) (An(F)) - A(F2)H >{) < n)
To address ii
note that for any F E span(F), both F163(T (oyo,,F((X)) F(X))) and F F(X))) are finite since MM a bounded metric space. By law of large numbers, (F) - A(F) = (1( for any F
span
For
ii) consider any F1, F2 E span(F), then
4 (A,(F)) A(R)) (A, (F2) A(FE))
- An, + (AFF)) - A(F2)
ddl(7,,.. FI(X)) 1100--gg(ho do( Pa(xx))lld((xxx.x (X().) + dg(nY... (2(X)))
+ (T ( -4 dgy
((diam(((M))) (F
=Op F2)).
Thus.
sup (x(Fi.Pa)sc
(Fi( - A(R)) (A, A(F)) Op(0))
which implies (ii) (1).
Finally, by Corollary 3.2.3 in (58 and Theorem | follows that d (FF F*)
H Comparison to XGBoost and SketchBoost using vectorized graph Laplacians
To further evaluate the effectiveness of FGBoost, we conducted an additional simulation study com- paring against popular variants of gradient boosting methods, under the same network simulation setup described in Section Since gradient boosting algorithms are not inherently designed to handle metric space-valued valued responses, we adapted them by vectorizing the graph Laplacians. Specifically, owing to the symmetry and zero row-sum constraints, each Laplacian fully determined by its strict
26

---

Table Average mean squared prediction errors and standard deviations (in parentheses) of Fréchet geodesic boosting (FGBoost), XGBoost and SketchBoost for network outputs.
FGBoost
XGBoost
SketchBoost
100
13.644 (3.140)
15.234 (3.319)
15.391 (2.825)
200
10.531 (3.371)
11.989 (2.686)
12.162 (2,421)
500
6.912 (1.950)
9.035 (2.145)
8.887 (1.869)
1000
5.471 (1.481)
7.096 (1.703)
6.793 (1.443)
upper triangular entries, which we flattened into vector-valued output. We considered baselines: (i) coordinate-wise XGBoost [13] where an independent regressor is trained for each coordinate of the vectorized output, and (ii) SketchBoost [26). recent multi-output gradient boosting method that jointly models all coordinates. For both methods, default hyperparameters were used. In contrast, FGBoost directly operates on the graph Laplacians as objects a geodesic metric space.
Table D presents average MSPE (with standard deviations) over 500 Monte Carlo replications for varying sample sizes. FGBoost consistently outperforms both XGBoost and SketchBoost across all sample sizes. While SketchBoost improves over XGBoost by leveraging joint modeling, its advantage appears only at larger sample sizes and it remains inferior FGBoost. These findings highlight key limitation of vectorization- based approaches: although they enable the application of standard boosting algorithms, they disregard the intrinsic geometry and structural dependencies of graph Laplacians. By directly respecting the non-Euclidean nature of the output space, FGBoost achieves substantial improvements in predictive performance.
I Choice of hyperparameters
The hyperparameters for Fréchet geodesic boosting can be selected using a grid search over the candidate values listed in Table The optimal combination of hyperparameters is chosen minimize the mean squared prediction error for the validation dataset.
Table 6: Hyperparameter settings.
Learning rate
0.01
0.03
0.05
0.1
Number of iterations
50
70
90
100
Depth each tree
2
4
5
J Training time and computational complexity
Computational efficiency is key consideration in the practical deployment of regression methods, particularly in modern applications involving non-Euclidean outputs such as probability distributions and networks. In this section, we provide systematic comparison of the training times of FGBoost and several state-of-the- art baseline methods across a range of sample sizes. All experiments were conducted on local machine equipped with an Apple M3 Max chip rumning macOS Sequoia.
Table reports the training times (in minutes) for sample sizes n 100, 200, 500, 1000, 2000. Among the baseline methods, global Fréchet regression (GFR) is consistently the fastest, as it generalizes linear regression to non-Euchi Non-Euclidean settings without introducing significant algorithmic complexity. However, this computational simplicity comes at the cost of substantially reduced model flexibility, since GFR imposes linearity assumptions that may be too restrictive in practice.
Random forest-based methods, Fréchet random forest (FRF) and random forest weighted local linear Fréchet regression (RFWLLFR). are also computationally efficient. Their parallelizable tree-based architectures facilitate fast training, especially on multi-core systems. In contrast, FGBoost trains trees sequentially, leading to a moderately higher computational cost. Nonetheless, this sequential nature allows FGBoost to iteratively correct model bias and effectively capture complex nonlinear relationships, which is particularly advantageous when modeling outputs in curved or high-variance metric spaces.
27

---

Table Training time minutes across different sample sizes.
n
FGBoost
GFR
SDR
IFR
FRF
RFWLLFR
100
0.54
0.003
1.27
2.24
0.24
0.24
200
0.95
0.006
3.89
7.00
0.55
0.55
500
2.41
0.020
19.68
25.65
1.56
1.55
1000
3.93
0.080
71.45
137.82
3.41
3.37
2000
7.14
0.270
6.75
6.77
Table 8: Average mean squared prediction errors and standard deviations (in parentheses) of Fréchet geodesic boosting (FGBoost), global Fréchet regression (GFR) [47, sufficient dimension reduction (SDR) [62]. single index Fréchet regression (IFR) 4. Fréchet random forest (FRF) [T0] and random forest weighted local linear Fréchet regression (RFWLLFR) 49 for National Health and Nutrition Examination Survey data.
FGBoost
GFR
SDR
IFR
FRF
RFWLLFR
0.054 (0.001) 0.059 (0.001) 0.065 (0.006)
0.071 (0.012)
0.058 (0.001) 0.073 (0.003)
Methods based on dimension reduction, including sufficient dimension reduction (SDR) and single index Fréchet regression (IFR), are substantially more computationally intensive. These methods involve iterative estimation of latent structures and repeated geodesic evaluations, resulting in poor scalability. At n 2000, the computational cost of SDR and IFR became prohibitive, and we were unable to obtain results within a reasonable time frame.
Overall, FGBoost achieves a favorable trade-off between computational cost and modeling flexibility. While it is not the fastest method in absolute terms, its ability to scale to large datasets and to accommodate complex, non-Euclidean output structures makes it competitive and practical choice in modern regression settings.
K Additional real-world data application: National Health and Nutrition Examination Survey
To further assess the empirical performance of FGBoost, we analyzed a fourth real-world dataset from the National Health and Nutrition Examination Survey (NHANES) 2005-2006. NHANES is large-scale survey that evaluates the health and nutritional status of U.S. adults and children through interviews, physical examinations, and laboratory tests. this cycle, participants aged six years and older were asked to wear an ActiGraph 7164 accelerometer consecutive days. The device recorded physical activity intensity in counts per minute (CPM) at 1-minute resolution, beginning at 12:01 am the day following the health examination and removed only for sleep, swimming, or bathing. These accelerometer data have been widely used to study the relationship between physical activity and health outcomes 31 24
the right hip for seven
We focused on modeling the distribution of physical activity intensity using demographic and health-related variables predictors. For each participant, activity values equal to zero or exceeding 1000 CPM were excluded, since zeros may correspond various low- activity states (e.g., sleep or device non-wear) and values above 1000 CPM are typically considered measurement artifacts. The remaining activity counts over the seven days were concatenated to form the empirical distribution of each participant's 's activity intensity. Similar distributional representations have been employed in recent studies [T 34 39 The predictor set comprised 13 demographic and anthropometric variables: gender, age, race/cthnicity, veteran status, education (college or above), household income (S 35,000), marital status, weight, height, body mass index (BMI,, thigh cireunference, waist cireumference, and upper arm length. To ensure data quality and reliable coverage, we selected the 200 participants with the most valid observations and performed 10-fold cross-validation over 20 runs for model evaluation.
non-Euclidean response,
Table reports the AMSPE for FGBoost and competing regression methods. FGBoost achieves the best predictive accuracy, outperforming all alternatives. These results demonstrate FGBoost's 's ability
28

---

to capture complex regression relationships when the outcome is an empirical distribution derived from high-frequency sensor data.
L Limitations
While FGBoost provides flexible framework for regression with metric space-valued outputs, it has several limitations. First, although boosting reduces bias and can help control variance, it is still susceptible to overfitting, particularly when the base learners are overly complex or the number of boosting rounds is large. FGBoost, we address this by limiting tree depth and applying early stopping, but the risk remains, especially in small-sample
high-noise settings.
Second, our theoretical analysis depends on the assumption that the output space is a Hadamard space, a condition met by many practical metric spaces, but nonetheless restrictive. Broadening the analysis to encompass more general geodesic metric spaces would improve the generality of the theoretical guarantees. From implementation perspective, however, FGBoost can be applied in any geodesic space. For example, our experiments on compositional data (Appendix (A] involve the positive hypersphere, which is not a Hadamard space, and demonstrate strong empirical performance. This suggests that the method is practically robust beyond the confines the Hadamard assumption, even though formal guarantees do not yet extend to these cases.
Third, while Section esstablishes strong convexity of the risk functional, we do not provide formal convergence proof for FGBoost. Classical convergence analyses for boosting rely on Banach space structures, where linear operations enable the use of Taylor expansions or Gâteaux derivatives [63]. In geodesic metric spaces, such tools are unavailable, and new geometric techniques would be required to establish descent guarantees. Developing these tools is an important avenue for future research.
Fourth, FGBoost is inherently more computationally intensive than scalar-based boosting methods due to the need for metric evaluations and Fréchet mean computations at each iteration. While many commonly used spaces admit closed-form solutions (c.g.. Wasserstein distributions, SPD matrices with power metries, networks with Frobenius metric) or efficient iterative algorithms (e.g.. proximal point methods) []] these steps still add overhead compared to simple arithmetic operations. Our implementation adopts modular design that separates metric-specce metric-specific primitives from the core boosting loop, allowing extensibility across different geodesic spaces. Nevertheless, developing fully optimized and unified library that achieves the efficiency established boosting frameworks such as XGBoost or LightGBM is an ambitious but promising avenue for future work.
Finally, while we extend SHAP values to interpret FGBoost predictions, the model's ensemble structure, built from numerous weak learners, makes it inherently difficult to interpret. This limits transparency and may pose challenges in domains where understanding the model's decision process is essential.
Future work could explore additional regularization strategies, such as penalizing leaf weights incorporating dropout-like mechanisms, aiming to enhance robustness.
29

---

M Additional tables
Table 9: Predictors of human mortality data.
Category
Variables
Explanation
1. Population Density
population per square kilometer
number of males per 100 females in the popu- lation
Demography
2. Sex Ratio
average age of mothers at the birth of their children
Mean Childbearing Age
4. GDP
gross domestic product per capita
percentage of agriculture, hunting, forestry, and fishing activities of gross value added
GVA by Agricultrre
consumer price index treating 2010 as the base year
Economics
6. CPI
percentage of unemployed people in the labor force
Unemployment Rate
8. Health Expenditure
percentage of expenditure health of GDP
Environment
Arable Land
percentage of total land area
Table 10: Predictors of New York City taxi network data.
Category
Variables
Explanation
1. Temp
daily average temperature
2. Humidity
daily average humidity
Weather
3. Wind
daily average windspeed
4. Pressure
daily average barometric pressure daily total precipitation
5. Precipitation
6. Mon to Thur Friday or Saturday
indicator for Monday to Thursday indicator for Friday Saturday
Day
8. Passenger Count
daily average number of passengers daily average trip distance daily average fare amount daily average tip amount daily average tolls amount
9. Trip Distance
Trip
10. Fare Amount
11. Tip Amount
12. Tolls Amount
30


Reviews:
[{"review_id": "Hw2cUzeLsi", "reviewer": "Reviewer_HbNk", "summary": "The paper proposes a new gradient boosting framework which adapts for non-Euclidean output spaces. The new framework starts from a Frechet mean of regression labels, and adds geodesics in each boosting round to the current prediction of each sample. FGBoost uses decision trees as base learners to partition the data into regions. Optimal geodesics to add are calculated for samples in each region as the Frechet mean of the residual geodesics of these samples. When the output space is a Hadamard space, there exists a consistent esitmator to the global unique minimizer of the boosting loss.\n\nNumerical experiments are conducted in both simulation datasets and real-world applications, covering probability and graph outputs. FGBoost shows superior performance over existing regression models for non-Euclidean outputs.", "questions": "1. Why are all the SHAP values shown in Figure 2 positive? Are the absolute values plotted? Or just the SHAP values of features with highest positive values are plotted?\n2. How the trees are constructed in FGBoost?\n3. In the experiment comparing with XGBoost, XGBoost is trained to fit each value of the Laplacian matrix. It would be interesting to see how the performance changes with traditional gradient boosting when fitting multiple values together, given the recent progress of multiple output gradient boosting decision trees (e.g., https://proceedings.neurips.cc/paper_files/paper/2022/file/a36c3dbe676fa8445715a31a90c66ab3-Paper-Conference.pdf).", "limitations": "1. It would be better to further address the limitations due to geodesic transport map.\n2. It would be better to discuss the potential of efficient implementation of FGBoost, given that existing gradient boosting algorithms are know for its highly optimized implementations for fast training.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 4, "strengths_and_weaknesses": "Strenghts:\n1. As a powerful out-of-the-shelf machine learning algorithm, gradient boosting has been applied to many regression/classification/ranking problems. It is definitely quite valuable to extend gradient boosting to more diverse output spaces. In this case, a geodesic space.\n2. Gradient boosting (with decision trees) has been always taken as black-box machine learning tools without prior knowledge applied in the algorithm. However, the results in Appendix H show the potential advantage of explicit geometry awareness through regression in the graph space, over simple vectorized numerial representations of graph with existing gradient boosting algorithm. I personally think this is interesting if the advantage can be more widely and thoroughly verified.\n3. The approach lays on solid theoretical foundations. The adaption from traditional gradient boosting is natural and clear.\n4. The paper is clearly written.\n\nWeknesses:\n1. Some details are missing, i.e., how the trees are constructed in FGBoost.\n2. It would be better if FGBoost could be evaluated over more real-world applications, as well as comparing with traditional gradient boosting approaches.\n3. The requirement of geodesic transport map may limit the application of the algorithm. For example, the geodesic transport map of distributions requires the calculation of quantiles, which is not trivial for most of the probability distributions.\n4. Calculation of Frechet mean is more expensive than calculation of mean of scalar values in traditional gradient boosting. Thus FGBoost  could be more expensive. It would be good to think how to implement FGBoost efficiently. Especially when facing different metric spaces, the definition of distance can be very different. It can be hard to propose a unified library or framework which can handle all metric spaces efficiently.", "was_revised": true, "final_justification": "Thanks the authors for the response. I will keep my score as \"accept\".", "rebuttal": "We are very grateful for your thorough review and positive assessment. We are particularly encouraged that you found the extension of gradient boosting to geodesic spaces to be \"**quite valuable**,\" our theoretical foundations \"**solid**,\" and the paper \"**clearly written**.\" Your suggestion about highlighting the advantages of geometry awareness is a key insight that we appreciate.\n\nWe thank you for the opportunity to clarify the missing details and discuss the limitations you raised. Below, we address each of your points.\n\n---\n\n> Tree construction\n\nYou raised an important point about missing details on how the trees are constructed. We will add a dedicated paragraph to the main text to make this procedure explicit.\n\nThe tree is constructed greedily, similar to a standard CART algorithm. At each node, the algorithm searches for the optimal split (i.e., the feature and threshold) that maximizes the reduction in the overall loss function. Specifically:\n1. For every potential split, the data is partitioned into two child nodes.\n2. For each proposed child node, we compute a single representative geodesic, which is the Fréchet mean of the pseudo-residuals ($\\gamma_{\\hat{Y}_i^{k-1}, Y_i}$) of all samples falling into that node.\n3. We then calculate the total Mean Squared Error (MSE) that would result from this split. The MSE is the sum of squared distances between the true observations $Y_i$ and the newly updated predictions. A sample's newly updated prediction is its current prediction updated with the representative geodesic from the child node it belongs to.\n4. The algorithm selects the split that results in the lowest total MSE. This process is repeated recursively.\n\nThis procedure is conceptually identical to standard regression tree construction, which chooses splits to minimize MSE. The key difference is that our \"MSE\" is defined by the metric of the geodesic output space, as per the loss function in lines 184-186.\n\n---\n\n> Expanding real-world evaluation\n\nThank you for your suggestion. We agree that including more real-world datasets would strengthen the empirical evaluation of the proposed method. In addition to the two real-world datasets presented in the main text (**human mortality distributional data** and **New York taxi traffic network data**), we present a third dataset in Appendix A.2, where we demonstrate the application of FGBoost for **emotional well-being compositional data**. To further enhance the empirical evaluation, we analyzed a fourth real-world dataset from the **National Health and Nutrition Examination Survey (NHANES)** 2005–2006, which we will include in the revision. NHANES combines interviews and physical exams to assess health and nutrition. During this cycle, participants aged six and older wore an ActiGraph 7164 accelerometer for seven days, recording physical activity intensity in 1-minute epochs, with values ranging from 0 to 32,767 CPM.\n\nWe modeled the distribution of physical activity intensity using demographic and health-related predictors. For each participant, we exclude physical activity intensity values above 1000 CPM or equal to zero, as zeros may correspond to various low-activity states (e.g., sleep or swimming). The remaining physical activity intensity values over the seven days are concatenated to form the empirical distribution of each participant's physical activity intensity. Similar analyses of physical activity intensity have been conducted in previous studies (Chang & McKeague, 2022; Lin et al., 2023).\n\nPredictors include 13 variables: gender, age, race/ethnicity, veteran status, education (college or above), income (≥$35,000), marital status, weight, height, BMI, thigh circumference, waist circumference, and upper arm length. We selected 200 participants with the most valid data to ensure reliable coverage and evaluated performance using 10-fold cross-validation.\n\nResults show FGBoost outperforms competitors in prediction accuracy, achieving the lowest error (0.0544) and consistently low standard error; see lines 232-236 for the comparison methods.\n\n_Table: Average mean squared prediction errors and standard deviations (in parentheses)._\n|FGBoost|GFR|SDR|IFR|FRF|RFWLLFR|\n|-|-|-|-|-|-|\n|0.0544 (0.0010)|0.0590 (0.0011)|0.0649 (0.0058)|0.0711 (0.0121)|0.0575 (0.0008)|0.0729 (0.0028)|\n\n---\n\n> Requirement for geodesic transport maps\n\nWe acknowledge that the need for a transport map is a key assumption. However, for many practical cases, this is not a major barrier. Regarding your specific example of distributional data, when working with empirical distributions constructed from observations, the transport map relies on **empirical quantile functions**. These are **simple step functions** and are computationally straightforward. For the broader class of Riemannian manifolds, the transport map is well-defined as parallel transport. We will expand the discussion in the appendix to cover these practical considerations.\n\n---\n\n> Computational efficiency and unified framework\n\nWe agree that computational efficiency is important. Our implementation is designed to exploit closed-form solutions where possible and use efficient iterative methods otherwise, while maintaining a modular structure. In FGBoost, the metric-specific operations (distance evaluation and geodesic calculation) are separated from the core boosting procedure, which remains unchanged.\n\nWhile computing Fréchet means introduces additional cost compared to scalar averaging in standard boosting, this step is not a major bottleneck in many common settings:\n- Closed-form or low-cost computation:\n  - **Univariate distributions with the 2-Wasserstein metric**: The Fréchet mean is the $L^2$ average of quantile functions.\n  - **SPD matrices with power metrics**: Computed using the Fréchet mean in the power-transformed space, i.e., by mapping matrices through the power transform, averaging in that space, and applying the inverse transform.\n  - **Networks with the Frobenius metric**: Computed as the Fréchet mean under the Frobenius norm, which reduces to the arithmetic mean of graph Laplacians.\n- Iterative methods for complex spaces: When a closed form is unavailable, established algorithms such as the **proximal point method** (Bacák, 2014) are efficient and converge quickly under convexity.\n\nWe acknowledge that FGBoost is inherently more computationally intensive than scalar-based boosting because of metric and Fréchet mean calculations. To ensure transparency, we report training time comparisons with other non-Euclidean methods in Appendix J.\n\nDeveloping a fully optimized, unified library for all metric spaces is an ambitious engineering task, similar to the evolution of frameworks like XGBoost and LightGBM from the original gradient boosting concept. Our contribution here is to provide a general algorithmic framework and demonstrate its feasibility across multiple spaces. Creating highly optimized implementations for specific domains is a promising direction for future work, and we will note this in the Limitations section.\n\n---\n\n> SHAP values (Fig. 2)\n\nThis is a sharp observation. In our adaptation of SHAP values for a geodesic space, we replace the standard (signed) difference in predictions with the **distance** between predictions, i.e., $d\\big(f_{S \\cup \\{j\\}}(\\mathbf{x}_{S \\cup \\{j\\}}), f_S(\\mathbf{x}_S)\\big)$ in lines 547-548. Since distance is always a non-negative quantity, the resulting feature contributions (our SHAP values) are also always non-negative. Figure 2 thus plots the raw SHAP values, which represent the magnitude of each feature's impact. We will add a sentence to Appendix C to explicitly clarify this.\n\n---\n\n> Comparison with multi-output gradient boosting\n\nThank you for this excellent suggestion and for pointing us to the relevant reference (Iosipoi & Vakhrushev, 2022). We implemented SketchBoost, a recent multi-output gradient boosting method, and compared it with FGBoost and the coordinate-wise XGBoost baseline under the same network simulation setup (see Appendix H). Since SketchBoost supports multi-output prediction, we vectorized the strict upper triangular entries of each graph Laplacian into a single vector-valued response and trained the model using default settings.\n\nThe table below reports average mean squared prediction errors (with standard deviations) over 500 Monte Carlo replications. FGBoost consistently outperforms both XGBoost and SketchBoost across all sample sizes. While SketchBoost improves over XGBoost by leveraging joint modeling, its advantage appears only at larger sample sizes and it remains inferior to FGBoost. These results suggest that preserving the intrinsic non-Euclidean geometry of the output space, as FGBoost does, is critical for accurately capturing network structure.\n\n_Table: Average mean squared prediction errors and standard deviations (in parentheses) for network outputs._\n|$n$|FGBoost|XGBoost|SketchBoost|\n|-|-|-|-|\n|100|13.644 (3.140)|15.234 (3.319)|15.391 (2.825)|\n|200|10.531 (3.371)|11.989 (2.686)|12.162 (2.421)|\n|500|6.912 (1.950)|9.035 (2.145)|8.887 (1.869)|\n|1000|5.471 (1.481)|7.096 (1.703)|6.793 (1.443)|\n\n---\n\n**References**\n- Bacák, M., 2014. Computing medians and means in Hadamard spaces. *SIAM Journal on Optimization*, 24(3), pp.1542-1566.\n- Iosipoi, L. and Vakhrushev, A., 2022. SketchBoost: Fast gradient boosted decision tree for multioutput problems. *Advances in Neural Information Processing Systems*, 35, pp.25422-25435.\n- Lin, Z., Kong, D. and Wang, L., 2023. Causal inference on distribution functions. *Journal of the Royal Statistical Society, Series B (Statistical Methodology)*, 85(2), pp.378-398.\n- Chang, H.W. and McKeague, I.W., 2022. Empirical likelihood-based inference for functional means with application to wearable device data. *Journal of the Royal Statistical Society, Series B (Statistical Methodology)*, 84(5), pp.1947-1968."}, {"review_id": "wuBQB8euOt", "reviewer": "Reviewer_aXWy", "summary": "This paper introduces Fréchet Geodesic Boosting (FGBoost), a novel method that extends the gradient boosting framework to regression settings with Euclidean predictors and non-Euclidean outputs residing in geodesic metric spaces. The key idea is to use geodesics—the shortest paths between points in these spaces—as a proxy for residuals, which are not well-defined in non-Euclidean domains. This viewpoint allows a smooth extension of the gradient boosting framework to non-Euclidean settings. The author(s) establish theoretical guarantees for FGBoost, including the existence and uniqueness of the solution and the consistency of the risk minimizer, under the assumption that the output space is a Hadamard space. The author(s) also provide SHAP values for interpretability of FGBoost. The method's improved performance against competitors is demonstrated through three simulation studies and three real-world datasets.", "questions": "I alluded to some of these questions in my strengths and weaknesses discussion:\n\n- The theoretical justification of the methodology seems grounded in Hadamard spaces. The paper does well in explaining the three settings presented. However, I wonder how robust the results are to metrics that are not those listed. Perhaps an experiment with a metric that does not lead to a Hadamard space could be helpful to see how robust this approach could be. A response here would clarify the practical boundaries of the method.\n\n- The general framework for geodesic metric spaces is great, but it seems a practitioner would have to derive the geodesic transport map themselves if their problem isn't one of the three examples in the paper. This might limit how widely the method is adopted. Perhaps more guidance on the principles for constructing these maps could be helpful for new users. Addressing this could make the work much more accessible.\n\n- The implementation of tree-based FGBoost relies on finding the Fréchet mean of residual geodesics in each leaf. In the analysis, the chosen metric for the space of geodesics, $d_\\mathcal{G}$ conveniently decouples the metric to align with the empirical loss. Does this decoupling on the geodesic space have any tradeoffs among computability/convenience and geometry of geodesics space? May not be necessary, but a comment or two could help clarify here.", "limitations": "Yes.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Quality:** \n\nI believe this paper is of high quality. The presentation of the method and its theoretical and experimental justifications are strong. The evaluation across three different non-Euclidean output types, plus several real-world datasets, makes the demonstrated performance advantages convincing.\n\n**Clarity:** \n\nThe paper is very clear and well-written. The main contributions are highlighted and previous work in regression in non-Euclidean settings is presented. The development of FGBoost is smooth and well-motivated, drawing a clear link that FGBoost is a generalization of standard gradient boosting methods for geodesic metric spaces. Both theoretical and experimental justifications are convincing, and the three clear examples are helpful. I also appreciate the discussions on scalability and limitations in the Appendix, though it would help with transparency if a reference to this was added in the main text.\n\n**Significance + Originality:**\n\nI’m not very familiar with the related literature, but it seems that the work is novel in tackling general regression problems with Euclidean predictors and non-Euclidean outputs using a boosting framework. However, the paper references Matsubara (2024), which has a very similar motivation and core insight, leveraging Wasserstein gradients in the space of distributions. This could diminish the novelty a bit, as a key idea of using a geometric proxy for residuals is captured in that prior work. Still, this paper's contribution is distinct in its generality beyond distributions and its theoretical guarantees.\n\nI think a couple points may limit the method's scope of impact. \n\n- While placing FGBoost in the general context of geodesic metric spaces is a strength, for actual usage, it seems a practitioner would have to derive the geodesic transport map if their problem is not one of the three examples provided. This may limit the adoptability of the method, and more guidance here could be helpful.\n- It is also not clear how other metrics might impact FGBoost, especially those that result in non-Hadamard spaces where the theoretical guarantees do not hold. An experiment exploring this could possibly strengthen the paper by defining the practical boundaries of the algorithm.\n\n[1] Takuo Matsubara. Wasserstein gradient boosting: A framework for distribution-valued supervised learning. In *Advances in Neural Information Processing Systems*, 2024.", "was_revised": true, "final_justification": "I thought the initial paper was strong, and the authors were able to address my questions. Their overall discussions with other reviewers were also very thorough.", "rebuttal": "We are very grateful for your thoughtful and encouraging review. We are delighted that you found the paper to be of \"**high quality**\" and \"**very clear and well-written**,\" with a \"**smooth and well-motivated**\" development of the FGBoost framework. We especially appreciate you noting the strength of our evaluation across multiple data types and the helpfulness of our examples.\n\nYour feedback is invaluable, and we will use it to further improve the paper. Below, we address your comments on novelty and scope, as well as your specific questions.\n\n---\n\n> Relation to Matsubara (2024) and novelty\n\nWe appreciate your observation regarding Matsubara (2024) and welcome the opportunity to clarify the key differences. Although both works share a geometric motivation, the problem formulations and objectives are fundamentally different.\n\nLearning target:\n- FGBoost generalizes regression by learning a map from predictors to **outputs in a general geodesic space**. The prediction target is the observed non-Euclidean object itself.\n- WGBoost (Matsubara, 2024) focuses on evidential learning and uncertainty quantification for **scalar-output problems**. It constructs an intermediate target for each observation, the posterior distribution, using Wasserstein geometry (see Example 1, Sections 4.1 and 4.2 in Matsubara, 2024).\n\nAlgorithmic and theoretical differences:\n- **Generality**: FGBoost applies broadly to geodesic spaces (e.g., univariate distributions, networks, SPD matrices, compositional data), whereas WGBoost is tailored to uncertainty quantification in scalar regression.\n- **Theory**: We establish risk convexity and the convergence of the empirical risk functional under Hadamard space assumptions, while WGBoost primarily offers methodological and empirical contributions.\n\nWe will update the Related Work section to make these distinctions explicit. Thank you for prompting this clarification.\n\n---\n\n> Practicality of geodesic transport maps\n\nThank you for highlighting this important point on accessibility and ease of adoption. While our paper illustrates FGBoost on four representative settings including univariate distributions, networks, SPD matrices, and compositional data (Appendix A), we agree that offering general guidance for practitioners will improve clarity.\n\nOur framework is inherently modular and unified: extending FGBoost to a new geodesic space primarily requires specifying (i) the distance function $d$ and (ii) the geodesic transport map. Both components are standard and easily accessible for most geodesic spaces of interest.\n- **Riemannian manifolds**: For smooth Riemannian manifolds, geodesics and their associated transport maps (e.g., parallel transport) are classical and generally well-studied, with closed-form or numerically stable solutions.\n- **Other spaces**: For discrete structures such as trees, geodesics are typically defined via shortest paths, which can be computed efficiently using standard algorithms (e.g., Dijkstra), or are explicitly available for specified metrics such as the BHV metric (Billera et al., 2001). \n\nWe will make these principles explicit in the revision and provide additional guidance in the Appendix to support practical adoption.\n\n---\n\n> Robustness beyond Hadamard spaces\n\nYou correctly note that our theoretical guarantees require the assumption that the underlying geodesic space is a Hadamard space. We would like to clarify that this assumption is primarily for the **theoretical analysis** (e.g., proving strong convexity of the loss and uniqueness of the solution). From an **implementation perspective**, FGBoost is applicable for any geodesic space. Our experiments on **compositional data**, presented in Appendix A, provide a practical demonstration of FGBoost in a **non-Hadamard space** (the positive hypersphere) and demonstrate good performance. This showcases the practical robustness of our method beyond the confines of the Hadamard space assumption. We will add a sentence to the main text to highlight this further. \n\n---\n\n> Metric on the space of geodesics\n\nThe choice of $d_\\mathcal{G}$ was indeed motivated by its convenient property of decoupling the Fréchet mean calculation into two independent, tractable problems. This makes the tree-based implementation efficient.\n\nThe important thing to note is that this choice does not introduce trade-offs or limitations related to the geometry of the underlying output space $\\mathcal{M}$. The metric $d_\\mathcal{G}$ is on the space of geodesics $\\mathcal{G}(\\mathcal{M})$, not on $\\mathcal{M}$ itself. The core of our algorithm, the pseudo-residuals and aggregation, operates directly on $\\mathcal{M}$ using its intrinsic metric $d$. The decoupling is purely a computational convenience for the leaf node estimation that does not compromise the geometric integrity of the approach. We will add a brief comment to this effect to clarify that this choice is for computability without geometric sacrifice.\n\n---\n\nFinally, we thank you for the formatting suggestion. We will correct the section heading capitalization in the final version. Your detailed feedback will significantly strengthen our paper, and we are very grateful. Please let us know if our response has sufficiently resolved your concerns, and we are more than happy to take any additional questions.\n\n---\n\n**References**\n- Billera, L.J., Holmes, S.P. and Vogtmann, K., 2001. Geometry of the space of phylogenetic trees. *Advances in Applied Mathematics*, 27(4), pp.733-767.\n- Matsubara, T., 2024. Wasserstein gradient boosting: A framework for distribution-valued supervised learning. *Advances in Neural Information Processing Systems*, 37, pp.118994-119024."}, {"review_id": "v6hb1OY94P", "reviewer": "Reviewer_XGWn", "summary": "The authors propose Fréchet Geodesic Boosting (FGBoost), a boosting method that works with outputs residing in unique geodesic spaces.\nSome brief theoretical analysis is conducted, establishing, Under the additional assumption of the geodesic space being a Hadamard space, that the underlying optimization is well-posed and a uniform convergence bound (showing that the minimizer of the empirical risk functional converges to that of the minimizer of the populational risk functional).\nA significant amount of numerical experiments are also presented, showcasing the superiority of the proposed FGBoost algorithm throughout.\nThe authors also generalize SHAP explanations to their models (a standard tool in practical applications of existing boosting models) in the supplementary material.", "questions": "I'd appreciate some clarifications on the points raised in the 'Strengths&Weaknesses' text:\n- How is each iteration's argmin solved?\n- Any blockers on a convergence proof?\n- Is working with the MSE-like objective fundamental to your approach?\n- And the clarifications on the proof of Proposition 2.", "limitations": "Limitations were appropriately discussed. (Except perhaps for a minor boundedness assumption that could be more explicitly stated, but which seems fairly mild.)", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 4, "originality": 4, "strengths_and_weaknesses": "The problem of boosting with non-Euclidean outputs is a very relevant one, and the focus on geodesic spaces is quite sensible.\nI find the work to be fairly easy to read, interesting and with convincing results.\nI am leaning towards acceptance, and just have a couple of points in which I would appreciate some clarification on.\n\n- On the step-wise fitting procedure: in lines 182-186, it is said that $\\hat{f}_k$ is chosen as the minimizer\n  $$ \\hat{f}\\_k = \\arg \\min\\_f \\frac{1}{n} \\sum\\_{i=1}^n d^2 ( Y\\_i, T\\_{\\hat{F}\\_{k-1}(\\mathbf{X}_i) \\oplus \\\\{\\nu \\odot f(\\mathbf{X}\\_i)\\\\}}(Y\\_0) ). $$\n  However, it is not 100% clear how this argmin is solved in practice; this is not a simple MSE minimization problem, which e.g. CART solves for trees; the outputs of $f$ gets added to the current model and then passes through the transport map $T$, to only then be piped through a (possibly non-MSE?) squared metric.\n  Overall it looks like quite a fairly tricky optimization problem, somewhat similar to those posed in the combinatorial boosting literature (e.g. [ExactBoost](https://proceedings.mlr.press/v151/csillag22a.html) and [IPBoost](https://proceedings.mlr.press/v119/pfetsch20a.html); in fact, it seems particularly close to Equation (4) in ExactBoost, which requires meta-heuristic approaches to be solved).\n  On the other hand, the authors say that they use tree fitting for FGBoost, so it feels as though I am missing something.\n\n- In the theoretical analysis, I felt the lack of a proof of convergence: FGBoost is trained in a greedy fashion, and so is (in principle) a descent algorithm. In Section 4 it is proven that the risk functional is strongly convex (and thus has a unique minimizer), but this need not imply that FGBoost will converge to it. While not critical, this would be a nice result to have. Is there some particular difficulty in proving such a result for FGBoost?\n\n- Other loss functions: one of the great boons of gradient boosting is that one can easily change the loss function being fitted without having to tweak the base learners. FGBoost, in comparison, seems hard-wired to an MSE-like objective based on the space's metric. Is this fundamental to the approach?\n\n- Minor clarifications on the proof of Proposition 2:\n    - It's not immediately clear why the first display equation (line 573-574 upper bound on $d_{\\mathcal{G}}^2 (\\gamma, \\gamma_t)$) holds; I presume it is due to the Hadamard space assumption?\n    - Why must the intermediate points $\\gamma_t$ (line 573) depend continuously on the endpoints $\\gamma_0, \\gamma_1$?\n    - Also, this proof invokes the diameter of $\\mathcal{G}(\\mathcal{M})$ (as does Theorem 2), which is stated to be finite because $\\mathcal{M}$ is bounded; however, I could not find this assumption stated elsewhere. (The assumption is fine, just needs to be stated more explicitly in the main text.)", "was_revised": true, "final_justification": "My main concerns (explaining how the tree fit is done, convergence guarantees and limited flexibility wrt loss) were appropriately clarified in the rebuttal phase, and corresponding improvements to the writing of the paper were promised. I believe both the problem and proposed solution of the paper are of interest to the NeurIPS community, and thus recommend acceptance.", "rebuttal": "We sincerely thank you for your positive assessment and detailed, constructive feedback. We are encouraged that you found the problem \"**very relevant**,\" the work \"**fairly easy to read, interesting and with convincing results**,\" and the focus on geodesic spaces \"**quite sensible**.\"\n\nWe appreciate the opportunity to clarify the points you raised regarding the practical optimization, the convergence proof, and the choice of loss function. Below, we address each of your weaknesses and questions.\n\n---\n\n> Practical implementation of the `argmin` step for $\\hat{f}_k$\n\nThank you for this insightful comment. We will use your excellent feedback to improve the manuscript's clarity. You are correct in observing that the objective function for the base learner $\\hat{f}_k$, as written in the paper,\n\n$\\hat{f}_{k}=\\underset{f}{\\operatorname{argmin}}\\frac{1}{n}\\sum\\_{i=1}^nd^2(Y_i, T\\_{\\hat{F}\\_{k-1}(\\mathbf{X}_i)\\oplus\\{\\nu\\odot f(\\mathbf{X}_i)\\}}(Y_0))$\n\nis a complex, nested optimization that is not directly solved by standard tree-fitting algorithms.\n\nThe key detail, which we will make more explicit in our revision, is that we employ a practical, greedy approximation common to the gradient boosting framework. The procedure is a two-stage process:\n- **Compute pseudo-residuals**: At each iteration $k$, we first compute the target geodesics connecting the current prediction to the actual observation, $\\gamma_{\\hat{Y}\\_{i}^{k-1}, Y_{i}}$. These geodesics act as pseudo-residuals for the current fitting stage.\n- **Fit a tree to pseudo-residuals**: We then train a standard decision tree on the data pairs $\\{(X_{i}, \\gamma_{\\hat{Y}\\_{i}^{k-1}, Y_{i}})\\}_{i=1}^{n}$. For each terminal leaf, the prediction is the Fréchet mean of the pseudo-residuals within that leaf. As we note (lines 192-194), this calculation is computationally feasible because the optimization decouples into two independent components.\n\nThis also clarifies the distinction of our proposed method from the combinatorial boosting methods like `ExactBoost` and `IPBoost` that you mentioned. In our FGBoost framework, the base learner is trained on a **pre-computed target**, the pseudo-residual geodesics. This differs from combinatorial methods where the base learner itself must often solve a more complex optimization problem (e.g., finding an optimal subset or partition) rather than fitting to a simple, pre-defined target.\n\n---\n\n> Convergence guarantee\n\nThis is an excellent point. The primary difficulty lies in the geometric setting. Classical convergence proofs for boosting (e.g., Zhang & Yu, 2005) rely on tools like **Taylor expansions** or **Gâteaux derivatives** to analyze the **functional descent** at each step. These techniques require at least a Banach space structure, where linear operations like addition, subtraction, and scalar multiplication are well-defined. In contrast, our setting involves general geodesic metric spaces, where such operations are not available. As a result, these analytical tools cannot be directly applied without imposing much stronger assumptions. We believe that developing the geometric tools necessary for such a proof is a significant undertaking and an exciting direction for future research. We will add a discussion of this limitation and future direction to the paper.\n\n---\n\n> Flexibility of the loss function\n\nYou raise an excellent and fundamental point about the flexibility of the loss function, which is indeed a hallmark of the standard gradient boosting framework. Our approach is best understood as a direct geometric generalization of the most common boosting setup: using squared error (MSE) loss for regression.\n\nIn that classic Euclidean setting, the negative gradient of the loss function simplifies to the residual ($Y_i - F_{k-1}(X_i)$). The base learner, typically a tree, is then tasked with fitting these residuals.\n\nFor outputs in a geodesic metric space, we propose that the most natural and principled generalization of the \"**residual**\" is the **geodesic** path connecting the current prediction to the true observation ($\\gamma_{\\hat{Y}\\_{i}^{k-1}, Y_{i}}$). Our use of the squared-metric loss is directly motivated by this analogy; it allows the target of the base learner to be this generalized residual, which is the direct counterpart to the negative gradient in the MSE case.\n\nConversely, other common regression losses, such as Huber or log-loss, are **not well-defined** for outputs in a general metric space as they rely on algebraic structures (e.g., subtraction, addition) not present in these spaces. Defining meaningful and well-behaved counterparts for these other losses in a geometric context is a significant research challenge in its own right. Therefore, our choice is a consequence of providing a direct and geometrically sound generalization of the canonical gradient boosting algorithm. We will revise the manuscript to make this motivation clearer.\n\n---\n\n> Minor clarifications on the proof of Proposition 2\n\nWe thank you for carefully reading the proofs and for these helpful questions.\n- **First display equation in Prop 2 (line 573-574).** Your presumption is exactly right. This inequality is a defining property of convexity in Hadamard spaces and follows directly from the assumption that the space is Hadamard. We will add a more explicit reference to make this clearer.\n- **Continuous dependence of intermediate points (line 573).** This is also a standard property of geodesics in Hadamard spaces; see Proposition 2.3 of Sturm (2003). The map that defines the geodesic path is continuous with respect to its endpoints. We will clarify this in the proof.\n- **The boundedness assumption.** Thank you for noting this. The assumption is stated in the main text in Section 2, line 78: \"Let $(\\mathcal{M}, d)$ be a bounded metric space.\". However, we agree that its importance for the proofs in the appendix could be made more prominent. We will reiterate this assumption in Section 4 before the theoretical results to improve clarity.\n\n---\n\nWe sincerely value the time you've dedicated to reviewing our work and deeply appreciate your insightful questions, which have enhanced the quality of our paper. Please let us know if our response has sufficiently resolved your concerns, and we are more than happy to take any additional questions.\n\n---\n\n**References**\n- Sturm, K.T., 2003. Probability measures on metric spaces of nonpositive curvature. *Heat Kernels and Analysis on Manifolds, Graphs, and Metric Spaces (Paris, 2002). Contemp. Math., 338. Amer. Math. Soc., Providence, RI*, 338, pp. 357–390.\n- Zhang, T. and Yu, B., 2005. Boosting with early stopping: Convergence and consistency. *Annals of Statistics*, 119(548), pp. 1538-1579."}, {"review_id": "KT0ZOsTyRk", "reviewer": "Reviewer_rLbt", "summary": "The paper Fréchet Geodesic Boosting generalizes gradient boosting to Riemannian manifolds by introducing a novel algorithm Fréchet Geodesic Boosting (FGBoost). Traditional boosting methods operate in Euclidean spaces, where residuals are computed as vector differences. However, in non-Euclidean settings, subtraction is not well-defined. FGBoosting replaces residuals with tangent vectors derived from the Fréchet mean, allowing the model to minimize a loss defined over geodesic distances.", "questions": "1. What assumptions are required for the manifold to ensure the algorithm is well defined? Like, completeness, injectivity radius, etc.\n2. How does the algorithm ensure consistency across tangent spaces at different base points during updates?", "limitations": "Yes", "rating": 4, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n1. The paper clearly identifies the limitations of classical boosting in curve spaces and provides a nice solution.\n2. This is a novel step in extending boosting to settings with non-Euclidean geometry.\n\nWeaknesses:\n1. It would be better a further theoretical discussion. The paper doesn't provide formal convergence guarantees.\n2. The algorithm relies on having nice exponential and logs maps, which can be a limitation on certain manifolds.\n3. Learning tangent vectors using regressors may introduce approximation errors, especially in curved regions where tangent space linearization is less accurate.\nFGBoost is a powerful method for learning on Riemannian manifolds and significantly extends boosting to curved domains but the results present computational challenges, limited theoretical analysis, and dependence on manifold-specific operations.", "was_revised": true, "final_justification": "The rebuttal is well explained and the authors solved the most of my concerns.", "rebuttal": "We thank you for your time and for providing detailed feedback on our work. We are encouraged that you recognized the novelty of extending boosting to settings with non-Euclidean geometry and appreciated our solution to the limitations of classical boosting.\n\nWe would like to begin by clarifying a central point about our paper's scope, as we believe it may resolve many of the concerns you raised. Our framework, FGBoost, is developed for the broad class of **geodesic metric spaces** (Bridson & Haefliger, 2013), which is a much more general setting than smooth Riemannian manifolds.\n\nWhile smooth Riemannian manifolds are an important example of geodesic spaces, our methodology does not require the additional structure of a manifold, and specifically does not in any way use tangent spaces or Riemannian exponential and log maps. Our algorithm operates intrinsically in a geodesic metric space using only the metric and the geodesics (shortest paths) it induces. This is why our \"pseudo-residuals\" are defined as geodesics themselves. We specifically note that they are not tangent vectors. This general formulation is a key strength, allowing FGBoost to be applicable for a wide range of data.\n\nWe only mention tangent planes and tangent vectors in Appendix A since in this special case the data are located on a sphere and we use Riemannian maps and tangent spaces in this special case to generate random data for the simulation but not for our FGBoost approach, which does not utilize these tools as they are generally not available. Indeed, for our primary experiments in Section 5 of the main paper, we use only the geodesic structure of the spaces and no tools that are specific to Riemannian manifolds.\n\nWe hope this clarification helps address your specific points below.\n\n---\n\n> Reliance on exponential/log maps and tangent spaces\n\nYour concerns about requiring \"nice exponential and log maps\" and errors from \"tangent space linearization\" would be valid points if our approach were using Riemannian tools. However, as explained above, this is not the case. As mentioned, since our method is formulated for general geodesic spaces, it does not use tangent vectors, exponential maps, or log maps. All operations are defined directly in the original metric space using geodesics.\n\n---\n\n> Required assumptions (completeness, injectivity radius, etc.)\n\nBecause we do not assume a manifold structure, our assumptions are on the metric space itself. Our primary theoretical assumption, stated in Section 4, is that the output space is a Hadamard space (a complete geodesic space with non-positive curvature, which guarantees the existence and uniqueness of the geodesics and Fréchet means). This assumption is only required for the supporting theory but not for the implementation of the method. Since there is no assumed manifold structure, there are no assumptions on completeness or injectivity radius in the context of a manifold. \n\n---\n\n> Consistency across tangent spaces\n\nAs our method does not operate in tangent spaces, the issue of ensuring consistency between them does not arise. All updates are performed intrinsically within the original metric space via the defined geodesic operations.\n\n---\n\n> Convergence proof\n\nThis is a good point. As we mentioned in the response to **Reviewer XGWn**, classical convergence proofs for boosting rely on tools like Taylor expansions or functional derivatives, which are not available in our general geodesic setting. Developing the geometric tools for the proof in general geodesic metric spaces (especially spaces that are not Hadamard spaces) is a challenging but important direction for future work. We will add a discussion of this limitation to the paper.\n\n---\n\nWe thank you again for your time and feedback. We hope our response has clarified the central scope of our work, focusing on **general geodesic metric spaces** rather than the more restricted case of smooth Riemannian manifolds, and in doing so, addressed the technical concerns you raised. We believe this **generality** is a key strength of our contribution. In the revised manuscript, we will make this distinction even clearer in the introduction and methodology to prevent future misinterpretation. We would be grateful for your feedback on whether these clarifications have provided a clearer assessment of our contributions.\n\n---\n\n**References**\n- Bridson, M.R. and Haefliger, A., 2013. *Metric Spaces of Non-Positive Curvature* (Vol. 319). Springer Science & Business Media."}]

Author rebuttal:

Meta-review: The paper proposed a gradient boosting procedure with is adapted to non Euclidean output spaces (geodesic metric spaces). Experimental validation indicates that the method performs favorably compared to competitors. The reviewers were initially enthusiastic, and the rebuttal allowed providing some clarifications. Please integrate all the suggested changes (especially regarding the tree fitting procedure) to the revised version.
Decision: Accept (poster)