Functional data analysis for multivariate distributions through Wasserstein slicing
Authors: Han Chen
Hans-Georg Müller
Venue: NeurIPS 2025 poster ()
Keywords: Distributional data analysis
multivariate distributional data
Radon transform
Sliced Wasserstein distance
Primary area: theory
TL;DR: We introduce an invertible Wasserstein slicing map that embeds multivariate distributions into a Hilbert space, enabling functional data analysis with theoretical guarantees and application to real-world data.

Abstract: The modeling of samples of distributions is a major challenge since distributions do not form a vector space. While various approaches exist for univariate distributions, including transformations to a Hilbert space, far less is known about the multivariate case. We utilize a transformation approach to map multivariate distributions to a Hilbert space via a Wasserstein slicing method that is invertible. This approach combines functional data analysis tools, such as functional principal component analysis and modes of variation, with the facility to map back to interpretable distributions. We also provide convergence guarantees for the Hilbert space representations under a broad class of such transforms. The method is illustrated using joint systolic and diastolic blood pressure data.

Functional data analysis for multivariate distributions through Wasserstein slicing
Han Chen Department of Statistics University California, Davis
Hans-Georg Müller
Department of Statistics University of California, Davis
Davis, CA 95616
Davis, CA 95616
nahchen&ucdavii.. edu
hgmueller@ucdavis.. edu
Abstract
The modeling of samples of distributions is major challenge since distributions do not form vector space. While various approaches exist for univariate distributions, including transformations to Hilbert space, far less is known about the multivari- ate case. We utilize transformation approach to map multivariate distributions a Hilbert space via a Wasserstein slicing method that is invertible. This approach combines functional data analysis tools, such as functional principal component analysis and modes of variation, with the facility to map back to interpretable distributions. We also provide convergence guarantees for the Hilbert space repre- sentations under broad class of such transforms. The method is illustrated using joint systolic and diastolic blood pressure data.
Introduction
Data that consist of multivariate distributions are becoming increasingly available. Instances of such data encompass anthropometric data distributions [24], stock price data distributions [19] and distributions of temperatures in climate modeling. Analyzing distributional data presents distinctive challenges, as they do not form vector space, rendering conventional methods designed for vector and functional data inapplicable.
For data that correspond to samples of univariate distributions, invertible global transformations such as the log quantile density transform and log hazard transform have been employed to map Hilbert space L [35]. While various maps were proposed to embed metric spaces into a Hilbert space, starting with Schoenberg's 's kernel embedding 43 4 and also including the centered log-ratio transformation [23] 30 24 these maps generally not 1:1 and therefore not invertible. However, invertibility is crucial as functional principal component analysis or regression in Hilbert space involves elements that can be located anywhere in the Hilbert space. To get back to the original distribution space, which is a must for interpretability and visualization, therefore requires invertibility of the map.
distributions
Another example such non-invertible embedding for the Wasserstein space of univariate distribu- tions is to use Riemannian type log maps map densities to tangent bundles. The resulting tangent spaces, are Hilbert spaces, where the development of principal component analysis [5) and regression models [2] is straightforward, however the reverse exp map is not defined on the entire Hilbert space but only on a convex cone and this poses a major problem, as only ad hoc solutions such as projecting elements outside the invertibility cone the cone available. Statistical regression models targeting univariate distributions as responses and Euclidean predictors can be implemented with Fréchet regression (36), where also other distribution-to doitttribuitttiiionn regression models with viable implementations have been proposed [6]
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

While these approaches have been developed for univariate distributions, multivariate distributional data have been less investigated; while kernel embeddings are still available, the invertibility problem remains. For regression with multivariate distributions, 34 introduce regression model but its scope limited to the case of multivariate Gaussian distributions, while [17] provides some theoretical results for more general settings for an optimal transport map regression approach in analogy to the case of univariate distributions [16.. Another approach utilizes a Bayes Hilbert space over bivariate domains to model time series of copula probability density functions [19], further developed for bivariate distributional data through spline representations for bivariate densities in the Bayes Hilbert space 24]. The Bayes Hilbert space proposals do not offer theoretical support and also do not address the estimation of densities from available random samples that they generate. Wasserstein geodesic principal component analysis represents variability among probability measures along geodesics in Wasserstein space, offering closed-form solutions for the case of Gaussian distributions and neural-network approximations for more general scenarios. However, it lacks explicit inverse map and does not guarantee invertibility for non-Gaussian distributions
45
Our goal in this paper is to extend the univariate log quantile density and log hazard transformation method [35] to the case of multivariate distributions. Then functional data analysis methodology. in- cluding functional principal component analysis (FPCA), functional modes of variation and functional regression can be applied in the Hilbert space into which the distributions are mapped. The densities corresponding to the Hilbert space representations are then obtained through a regularized inverse map [9] Theoretical results regarding the convergence of these representations in density space are derived under regularity conditions, drawing strength from known results for FPCA and addressing the complexities introduced by the transformations. Given that multivariate distributions are often not fully observed and need to be recovered from random samples generated by the underlying distributions, we also consider the challenges introduced by the estimation step.
2 Preliminaries
2.1 Univariate Density Transformation
Let G = (g E CZ) g > 0, So g(2)dz = =} be the space of strictly positive densities on a
common compact interval Z. Without loss of generality, assume the support for the univariate
density function is Z = [0. 1]. (55 propose a family satisfy assumptions (I.3) that are listed in Appendix A These transformations map continuous
transformations v G + L'(S)) that
univariate densities to a new space I.) where S is a compact interval.
Here we focus on the log quantile density (LQD) transformation, while the log hazard transformation another prominent example for this class transformations and is reviewed in Appendix D Choosing S = J0, IJ. the LQD transformation th 9 + L2( is given by
(s) = log (s)) s S
where G(z) = Jo g(u)du is the cumulative distribution function and G-1 is the corresponding quantile function supported on S. The inverse of a continuous function & on S can be written as exp(-v(G, (2))), where (s) S (u Since ( is not fixed, a boundary adjustment guarantees that the inverse will be density on Z [0, 1] 337.
(s) a [- cpfu) du,
(+) a w
where 0gp- = Se ep(u) du. The LQD transformation is an element of family of a larger class of univariate density transformations; however, similar transformations are currently not available for multivariate density functions. We fill this gap by introducing generalized multivariate transforma- tion in Section 3.2 which will be illustrated with examples based on multivariate extensions of the LQD transformation.
2.2 Density Slicing Transformation
To connect multivariate and univariate distributions, we use the Radon transform R 40) which an integral transform that takes an integrable p-dimensional function and maps it to infinite set of

---

integrals over hyperplanes in R'. Let 6 be unit vector in è - {2 R², ||>||2 = 1} and an element in R. We denote affine hyperplanes by 12.0 {ir E R2 (z.6) = and select an orthonormal basis {0, el: Cp-1) ) in RP such that (0, = 0 and (ej,et ررر, for j. .,.... 1. We then define the p-dimensional Radon transform R. through the integral over 12.0 as follows,
R()(0,2) RF-1
20 +
Σ
8yej
ds d5p-1 for 0 O and E R,
where is p-dimensional multivariate density function supported on D.
The inverse Radon transform has undergone thorough investigation, both theoretically and numerically [2 B1]. The commonly used method to reconstruct the original multivariate distribution from its Radon transform @ is through the filtered back projecrooottion given by
R (6) (x) 2(2 So [ (0(0,)) (v)e /ul=-' dod0,
(3)
where the Fourier transform is [w) Jg)ee -wwe du for all R. There are cases where Radon transforms are bona fide probability densities, while the inverse is not, such as the Wigner quasi-probability distribution in quantum state tomography [2]. Moreover, the inverse map R R-1 is not continuous, and small reconstruction errors in R(f) can be magnified [IS]. To address this, regularized inverse
[5
25
using filter function the Fourier domain is commonly applied,
()) 2( / (6(0,))) (v)e (0)) 1pi
Note that the regularization parameter balances approximation accuracy and the continuity of the inverse map (9. As this regularized inverse is not guaranteed to be a multivariate density function, normalization is applied to map the resulting L function into the multivariate density space F via
()))) Sp RR (x) de |R de > O,
R₂ (x) (z) 1 [
(5)
otherwise.
Here [D] represents the Lebesgue measure of D which is the domain of the multivariate density function f. Note that the Radon transform and its inverse part of generalized slicing transform family, for which an explicit invertible embedding between multivariate and sliced univariate densities was introduced in [9] where further details can be found.
Multivariate Density Transformation
3.1 Assumptions for Multivariate Densities
By convention, ( denotes the set encompassing all continuous functions while C2 denotes the consisting of all continuously differentiable functions. Let || ll be the L metric and || ll be the uniform metric for measurable functions. Throughout, Co. Ci.... represent various constants, and their dependency on pertinent quantities will be denoted by writing Cole), Ci (g).
We assume that the multivariate distributions that we consider have density functions and that their common support is C RP, where
(DI) The support set D C Rº is compact and convex. Denote the space of strictly positive multivariate density functions R with compact support D by
F { ((x) E L' LR(PP) ((x) > 0, [ [(x)do 1, support(f)) = D. satisfies (FI)
(FI) For all f F, there exists constant Mn > 0 such that max (1/ |11//1l=> < Mo and for 2 p+ 1, is continuously differentiable of order on C RR and has uniformly bounded partial derivatives.

---

v
R())
((R(()00..)
A2
R 2)
R
X1 02
((((F((22))
100)
Figure I: Scheme for the Radon Log Quantile Density (RLQD) transformation for the bivariate case: Slicing the bivariate density function along each projection (represented by an angle) through the Radon transform, followed the slice-wise log quantile density transform.
3.2 Multivariate Density Transformation
We introduce multivariate density transformations V F (e x S) from multivariate density functions to Hilbert space as follows. Let w. ab be the generalized univariate density transforma- tion and its inverse reviewed in Section 2.1 then
(8) s) (R(f) (6, }) (s). for (0.s) e f
(6)
By using the regularized inverse of the Radon transform, the regularized inverse can be defined for any continuous function Ç as
(() R, ((V- (=)) (z).
(7)
A concrete example is the Radon log quantile density (RLQD) transformation Vo. Vo(s)(0,s) s) +/Q (R() (0, (s) log (R( R (s)) for (0, E e FE,, (8) where R2(2) = (e, u)du is the cumulative distribution function along when fixing the first argument 0 and R, (s) the corresponding sliced quantile function supported on [0, 1]. The inverse map can be defined for any continuous function as
V (S)x() R.² (vO² (S(0,-)) (2) (z).
The steps defining the RLQD transformation are visualized in Fig.
For the theoretical analysis of these transformations require certain smoothness and regularity conditions. We show that when the univariate density transformation vp satisfies assumptions (LI)- (L3) in the Appendix the corresponding multivariate density transformation Yi according (6)
satisfies the four smoothness properties (T1)-(T4) in Appendix
Representations for Transformed Distributions
4.1 Metric for Multivariate Distributions
Let B(D) be the Borel o-algebra on D. We focus on random probability distributions on (D, B(D)) and each probability distribution is characterized by density function with respect Lebesgue measure. While classical metrics such as L3, L2_ L 00 or total variation metric [14] are commonly
used, recent work has emphasized the Fisher Rao 13 50) and Wasserstein metrics [18 32 26 for
samples of distributional data. The 2- Wasseerstein distance between two densities defined as
dy (f1. f2)
inf E Zzll3) for all fi, E F.
22)
(2, 1222)/1

---

Here Z1 and Z2 are random variables on R", and II (f1. f2) is the space of joint probability measures on D D with marginal densities fi and fz. For mne-iimnnsional distributions, the 2- 2-Wasserstein metric corresponds to the 1,2 distance between their quantile functions,
= / Go) (5))² ds, for all 91.92 9,
(01
where GT GT G are the corresponding quantile functions of 91. 92. Obtaining the Wassestein distance for the multivariate distribution involves intricate algorithms and does not have analytic solution 33 38.. An alternative approach the sliced Wasserstein distance [7]. It offers a computationally more efficient solution compared to the conventional Wasserstein distance and is defined
(S,,f, = = [f
div
d0, for all f1. F,
(10)
Jee
where the Radon transform R is introduced in 2
4.2 Transformed Distribution Representations
Suppose the random densities generated as f 3 where is probability measure F. Applying the transformation in Section
Sititonnyyiildd yjields a random element in the Hilbert space X == (f) E ( x S), which corresponds to the realization of stochastic process. This enables functional data analysis 4 and linear representations the Hilbert space L with results mapped back to the density space via the regularized inverse 1 that was introduced in (9).
Denote the mean function the process X by v(0, 8) = (0, x)) and the covariance function by H (01, $1 82 $2) Cov (X (01, $1), (02, s2)), where the eigenfunctions of the covariance operator associated with the kernel H are denoted by {ox (0, } and their corresponding eigenvalues by { The transformed process can be represented through the Karhunen-Loève expansion 27
X 8) = v(),s) s)+ +
Σ
XA 7A (0,s), (0,s)
(11)
k=
with the principal components X = (X((s - v(0 (0, (0, s)dOds. This expansion exhibits an optimality characteristic where the initial K terms constitute the K K-dimensional representation of (s,,t) t)) that minimizes the unexplained variance. Empirical estimators é, H, 2kk A Y& can be obtained based on the sample of transformed processes [4] 20 3 47). For extensions to multivariate functional data, including product and marginal FPCA, see
Truncating the Karhunen-Loevv expansion to the first K components yields an I.2-optimal approxi- mation; with judicious choice of K it captures pre-specified fraction of the variation,
(x( = ((. + Σ (0,) s), (0,s) e
(12)
To interpret the contribution of each principal component in the original density space, we extend the notion of univariate transformation modes of variation [35] to the multivariate setting. Specifically, we compute modes of variation in the transformed I space and map them back the density space via the regularized inverse. The k-th transformation mode of variation is defined as family of functions indexed by a E R.
my (I, = -; (u+ (
(xr).
(13)
The data application in Section 8 demonstrates that transformation modes of variation play an important role in understanding the effect of each principal component.
5 Estimation
In practice, random densities {f:}__ are often not fully observed and one has only have observations generated by each of these densities. In this more realistic case, preliminary density estimation step
5

---

0 is constant and p denotes the is needed, which is described in Appendix [F) If the random densities fi are estimated by kernel density estimator f. through (31) the corresponding sample mean and variance are based on the new transformed sample {X = (V(())]
7
1
D(0, s) (0,ss,
Σ
s1) 82).
(14)
Another application of Mercer's theorem implies
20
H0 $1, 02 82)
01 82),
Σ
(15)
where dx and 9/: are estimated eigenvalues and eigenfunctions of the covariance operator associated with the covariance function H. The principal components Xik are estimated as
Ass - . (X. 8) E(0, s)) (0, s)d0ds
(16)
and the plugin estimator for the transformation modes of variation 13
my (z, a, T) = -,,² 1 ( Aańi (I).
(17)
Selection of tuning parameter K and are further discussed in Appendix D
6 Asymptotic Convergence
In this section we establish asymptotic properties of the empirical estimators proposed in Section A key feature in the convergence rate is the spacing between eigenvalues, where we define
O
min (A) )j+1). 1SJSK
(18)
for any integer k. For the 2th distribution, we assume there are N, independent observations available to estimate the density function fi in (31. with details in the Appendix. The lower bound (n) minigign N; is relevant for the overall asymptotic convergence rate. For the following result, we require that the sample size per distribution is large enough so that the effect of density estimation can be neglected when deriving the asymptotic rates. Specifically, we require
(E1) lim N(n)/n" 2 e for some r > + p/4, where dimension of the domain of the densities.
Convergence analysis of the representations is presented in Appendix G For the following, we require an additional assumption (SI) which imposes restrictions on the eigencomponents. It reflects typical assumptions of exponentially or polynomially declining eigenvalues.
(S1) The sequence K = ((n) of approximating K re-rpresentations is such that O 2==" 11122 O(1). With C2 as in (T5) in Appendix(A) [A], we use the following bound to derive the convergence of transformation modes of variation, where for fixed constant og 0.
K) max sup Ca (T. 11 ||x (0, s) as) doo (Xe,or ISkSK la|Sa0
(19)
=0 T max ISKCK
Tx (0
20 (Xearde.a)))
with X4,a = = + Ar7k and Xy.a = 0 av An Under (SI), we have doo (Xr.rr Xxa) = (1). When is fixed, a| and ||0Xx,a (0, ()0.s) as can be uniformly bounded, implying S(T. K) = Op())) When K = K(n) + co, the convergence of S(T, K) requires bounded- ness of || lo and |0x,, 8)/a=lloo. which is implied by the boundedness of |/ and | (0, s)/Oal l The following result provides the rate of convergence of the estimated transfor- mation modes of variation
to the true modes in 13 uniformly over the parameter (a] S 0g
and for any constant ag > 0.

---

Theorem 1. Under assumptions (DI), (FI), (TI)-(T4), (EI), (KI)-(K),, for fixed K and on. my (z, a, T) and m a,T) r) as per 13 and (17).
sup dsw (ma(x, a, (, = 00((Tn IgkgR |a|saa
If we assume K = (n) 7 3o satisfies (SI), then
max (gkgk(m)
sup dsw (ma (2, a, ix, a3,
|a|Sav
max = OpT
{e [e 10 8) /00 Jo 1/1),
where T, Yk. σK are as in equations
and (88)
Corollary 1. When assuming exponential decay for the eigenvalues, i.e., A cc ce-kk for some constants 0 > O. with K (n) = -o log(n)] + cc, Theorem provides that
max sup dsw (m=(a, a, (x, (n) SalSaa
max Op(T"
(6lhall 10,x(0, s)/0s|l<0) 1/4).
1gRSk
We also present the convergence rate of the truncated density estimator in 25 See G for additional asymptotic convergence results in the Appendix.
Theorem 2. Under assumptions (DI), (FI), (TI)-(5)) (EI), (KI)-(K)), assume exponentially decaying eigenvalues 2 cc 0k for some constants e and O. With K ((n)
L log(n)] oc,
max dsw (AO),.(KKTT) =0p 18i5
10007-10 (T-(M-P) + T nn
and choosing T 2/6k leads to max daw ()))) ssih
= =0p (n- (n-(p=))/o))
7 Experiments
We illustrate our methods by comparing the performance of three approaches: (1) Directly applying FPCA and modes of variation to estimated densities viewed as elements of L ignoring density constraints 48). The resulting functions are truncated to their positive part and renormalized to ensure valid densities; (2) Performing PCA in the tangent space (TPCA) of the Wasss-Wassersteein manifold of covariance operators at a reference density and maps the resulting components back to the density space via the exponential map [29] and (3) the proposed RLQD method, which takes the constraints into account and where K-component approximations as well as modes of variation are guaranteed to reside in the density space and additional adjustments are not necessary.
We conducted Monte Carlo experiments under two simulation settings: (A) random bivariate normals and (B) mixtures of two bivariate normals. Details of the data generation process appear in Table In both settings, we generate 50 random densities, each from latent parameters, and sample 200 observations per density. As evaluation criterion we use the fraction of variance explained Vk / /V as defined in Section E. where a higher fraction of variance explained for the same choice of K reflects better representations. Boxplots depicting the fraction of variance explained with respect to the sliced Wasserstein distance for FPCA and RLQQ approaches are in Figure 4 It is clear that in comparison with both
these settings the proposed RLQD method performs the best throughout the standard FPCA approach and the tangent FPCA approach. The advantage of RLQD becomes particularly evident beyond the first two modes, where FPCA's explained variance plateaus and TPCA yields diminishing returns. In contrast, RLQD continues capture meaningful higher-order variation, demonstrating stronger structural representation and efficiency with fewer retained components.
We further visualize the modes of variation for the first two principal components under Setting A in Figures 2 and 3 Under the Gaussian assumption, both the proposed RLQD transformation and the classic FPCA exhibit similar patterns of variation. Specifically, the first mode primarily captures the overall intensity or concentration the density around the mean, whereas the second mode reflects a rotational change in the orientation of the joint distribution.
7

---

Table Simulation settings for bivariate normal and mixture distributions.
Setting
Random components
Resulting density
14 ~ -
( 200 ~ N (1i, 2,) j=1
P: ~ (004400.5
truncated on [-2, 2]
log $1: ~ -0.3
= 50
log S2i) ~ U[3. 0.3
8i 81:82104
8168280:
$21
Ml, M2,i 0.22,
2000
3 (ui E1,i)
P1 22, [40.4,0088
22N (u2, Σ2,5)
11,181
log 82,16)
~ -.33,
truncated on 2- 2)
log (51,21)))
log 22,21)
U[ 0.3
l..... 50
$1.1
11,16$1,2:P3
Σ1,
$1
$
2,11
2,16
82,1882,287
$2.21
gammal,-2
gammal,-1 -1
gammall 1,0
gammal. 1,1
gamma 1,2
5.0
2.5
density
0,0
9amma2.7
gamma2,-1 -1
g2mmeeee
gamma2,1
gamma2.2
5.0
2.5
0.05
0.0
-25
0.00
-5.0
-5.0
-1.
0.0
2.5
5.0-5.0
-2.5
0.0
255
5.0-500-5
0o
255
5.0-5.0-11
0.0
2.5
5.0-5.0-25
0.0
2.5
5.0
Figure Proposed RLQD transformation modes of variation as per 13 for Setting
gammal,-2
gammal,- -1
gammall
gamma1.3 1,1
gamma1.2
5.0
2.5
density
0.0
-255
-5.0
0.10
gamma2,-2
gammaz,-1 -1
gamma2,0
gamma2,1
gamma2,2
0.05
25-
-255
0.00
-500
-5.0
-1.5
0.0
2.5
5.0
5.0
-.5
0.0
2.5
5.0-5.0
-15 ol
2s
5.0-5.0
-1.5
o.o
2.5
5.0-
0.0
2.5
5.0
Figure 3: Classic FPCA Modes of variation for Setting A
Data Illustrations
We analyze blood pressure data from the Baltimore Longitudinal Study of Aging (BLSA) https://wwwbblsa.n gov. The dataset comprises 16,715 bivariate measurements (SBP, DBP) for approximately 2,800 individuals, with the number of visits per individual ranging from to 26 and ages spanning from 17 75. For each age group, we estimate joint two-dimensional densities using kernel density estimator. The estimator applied over 51 cquidistant grid points in each direction,

---

0.90
11:00
0.05
0.60
0.80
0.55
0.30
0.75
0.50
0.70
RLOD
(a) Setting A,
(b) Setting A, K
(c) Setting A,
0.40
0.00
0.50
0.65
0.55
0.70
(d) Setting B.
(e) Setting B,
(i) Setting B, K
Figure 4: Boxplots show the fraction of variance explained (via sliced Wasserstein distance) across 100 simulations. Each row represents a different setting, with varying component numbers K. Larger values indicate that higher fraction of variation explained and are preferred. The proposed RLQD approach is compared with unconstrained FPCA and tangent PCA for densities.
covering the domains [50, 205] for SBP and [40, 125] for DBP with Gaussian kernel and bandwidths hsap 24 and IDPP = 15. The fitted densities (Figure in the Appendix B) show a clear positive correlation between SBP and DBP, along with increasing variability older ages, highlighting the need for low-dimensional representations that capture these distributional shifts.
Using the RLQD approach with K 2 components, achieve strong dimension reduction while explaining 87% of the variance in sliced Wasserstein distance as shown in Figure 8 of the Appendix. The corresponding transformation modes of variation in Figure 5 provide insight into principal patterns of change. The predominant mode of variation in the top row reflects a horizontal mean shift coupled with increasing variance, primarily in systolic blood pressure. This suggests that the primary source of variation is shift in the mean and increased variance of systolic blood pressure. The second mode of variation highlights changes in the variation of diastolic blood pressure and in the covariance structure between diastolic and systolic blood pressures. The first two FPCA components explain 70% and 11 11%. Under standard threshold of 85% cumulative variance, RLQD requires only components for effective dimension reduction, whereas FPCA requires at least 3. The cigenfunctions and modes of variation produced by FPCA typically resemble cauussian-saapoe perturbations around the mean as presented in Figure of the Appendix, which can be restrictive when modeling data with inherent skewness more flexible, non non-Gaun Gaussian modes that better reflect the underlying distributions. For example, the blood pressure data exhibit clear asymmetry, with noticeable skew toward lower pressure levels. This is statistically confirmed by Mardia's 's test (28 for multivariate normality: both skewness (3063.39) and kurtosis (42.75) statistics are highly significant (p-values < 0.001), indicating a strong departure from multivariate normality. These findings demonstrate that RLQD not only provides more informative components with fewer dimensions but also produces representations that are more faithful to the structural properties of complex real-world distributions.
Gaussian features. In contrast, RLQD captures

---

2- 2-TEuuLs S
2
$27
000
density
0.00100
T- crunna
commani
0.00075
125
dascolicc
100
0.00025
75-
0.00000
50
100
150
200
50
100
150
200
50
100
150
200
100
150
200
50
100
150
200
systolic blood pressure
Figure 5: Transformation modes variation as per (33) for bivariate blood pressure distributions. The first row illustrates the action of the first mode of variation m and the second row that of the second mode variation mg as the values are varied from -2 to 2 as per 17
One can also the K K-component representations- with K 2, yielding bivariate score vectors-to characterize evolving blood pressure patterns. Local linear smoothing of the scores against age is shown in Figure using a Gaussian kernel with bandwidth selected via leave-one- out cross-validation. The first score, X1. associated with ?1, shows an increasing trend with age, suggesting rising blood pressure variation and increased diastolic pressure in older individuals. The second score, X2, linked to 72, increases until midlife then declines, indicating that diastolic pressure peaks in midlife. additional case study on temperature data in the United States provided in Appendix H
Scorel
Scora
(Ycar)
(Year)
(a) Scores associated with
(b) Scores associated with
Figure 6: Blood pressure data score smoothing. Local linear smoothing score representations
along age. Left: Smoothed first scores XI versus age. Right: Smoothed second scores X2 against age.
Conclusion and Limitations
This work presents transformation-based framework for the functional analysis of multivariate distributions. By combining Radon slicing with univariate log quantile density transforms, we map multivariate densities into Hilbert space where classical FDA tools such FPCA and modes of variation become applicable. The regularized inverse transform ensures that the structure of the original density space is preserved. While the focus of this paper, we that this framework facilitates distributional regression once vectorized representations are obtained. Limitations include the dependence accurate density estimation from finite samples that is subject the curse of dimensionality and the need for regularization parameter.
Acknowledgement and Disclosure of Funding
We would like to thank the reviewers for their constructive feedback. This research was partially supported NSF grant DMS-21 23155000
10

---

References
[1] Habti Abeida, Qilin Zhang. Jian Li, and Nadjim Merabtine. Iterative sparse asymptotic minimum variance based approaches for array processing. IEEE Transactions Signal Processing, 61(4):933-444, 2012.
[2] Luis Artiles, Richard D Gill, and MI Gut a. An invitation to quantum tomography. Journal of the Royal Statistical Society Series B: Statistical Methodology, 67(1):109-134, 2005.
[3] Juan Lucas Bali, Graciela Boente, David E Tyler, and Jane-Ling Wang. Robust functional principal components: A coojection-pursuit approach. The Annals of Statisties, 39:2852-2882 2011.
[4] Philippe Besse and James O. Ramsay. Principal components analysis of sampled functions. Psychometrika, 51(():285-311 1986.
[5] Jérémie Bigot, Raúl Gouet, Thierry Klein, and Alfredo López. Geodesic PCA in the Wasserstein space by convex PCA. Annales de I'Institut Henri Poincaré B: Probability and Statistics, 53:1-
26, 2017.
[6] Nicolas Bonneel and Hanspeter Pfister. Sliced Wasserstein barycenter of multiple densities. Technical report, Harvard Technical Report. TR-02-13. ftp://ftp. hwtpttpppeas.aaaeeee deeas.araaar.eed.. waaarard..duueechrppotttttt 02- pd....pdf. 2013.
[7] Nicolas Bonneel, Julien Rabin, Gabriel Peyré, and Hanspeter Pfister. Sliced and Radon Wasser- stein barycenters of measures. Journal of Mathematical Imaging and Vision, 51(1):22-45, 2015.
[8] Denis Bosq. Linear Processes in Function Spaces: Theory and Applications. Springer- Verlag, New York, 2000.
[9] Han Chen and Hans-Georg Müller. arXiv:2306. 10601, 2023.
Sliced Wasserstein regression.
arXiv preprint
[10] Han Chen, Hans-Georg Müller, Vasilis G Rodovitis, Nikos T Papadopoulos, and James R Carey, Daily activity profiles over the lifespan of female medflies as biomarkers of aging and longevity. Aging Cell, page e14080, 2024.
[11] Kehui Chen, Pedro Delicado, and Hans-Georg Müller. Modeling canntiio-valucee stochastic processes, with applications to fertility dynamics, Journal of the Royal Statistical Society, Series B (Theory and Methodology), 79:177-196, 2017.
[12) Yaqing Chen, Zhenhua Lin, and Hans-Georg Müller. Wasserstein regression. Journal of the American Statistical Association, 118:869-882, 2023.
[13] Xiongtao Dai. Statistical inference on the Hilbert sphere with application to random densities. Electronic Journal of Statistics, 16(1):700-736, 2022.
[14] L. Devroye and L. Györfi. Nonparametric Density Estimation: the LI view. Wiley, 1985.
[15] Charles L Epstein. Introduction to the Mathematics of Medical Imaging. SIAM, 2007.
[16] Laya Ghodrati and Victor M Panaretos. Distrbbution-nn-dissribttuoo regression via optimal transport maps. Biometrika, 109(4):957-974, 2022.
[17] Laya Ghodrati and Victor M Panaretos. Transportation of measure regression in higher dimen- sions. arXiv preprint arXiv:2305. /7503, 2023.
[18] Alison Gibbs and Francis Edward Su. On choosing and bounding probability metrics. International Statistical Review, 70(3):419 2002.
[19] Dominique Guégan and Matteo lacopini. Nonparametric forecasting of multivariate probability density functions. arXiv preprint 18X3:v1803.066 2018
[20] Peter Hall and Mohammad Hosseini-Nasah. On properties of functional principal components analysis. Journal of the Royal Statistical Society: Series B, 68(1):109-126, 2006.
11

---

[21] Sigurdur Helgason and S Helgason. The Radon Transform, volume 2. Springer, 1999.
[22] Gabor T Herman. Fundamentals of Computerized Tomography: Image Reconstruction from Projections. Springer Science & Business Media, 2009.
[23] K Hron, A Menafoglio, M TTepll K Hruzova, and P Filzmoser. Simplicial principal component analysis for density functions in Bayes spaces. Computational Statistics and Data Analysis, 94:330-350, 2016.
[24] Karel Hron, Jitka Machalová, and Alessandra Menafoglio. Bivariate densities in Bayes spaces: Orthogonal decomposition and spline representation. Statistical Papers, 64:1629-1667, 2023.
[25] Avinash Kak and Malcolm Slaney. Principles of Computerized Tomographic Imaging. SIAM, 2001.
[26] Leonid V Kantorovich. On the translocation of masses. Journal of Mathematical Sciences, 133(4):1311-1382 2006.
[27]
Kari Karhunen. Zur Spektraltheorie stochastischer Prozesse. Ann. Acad. Sci. Finnicae, Ser. A, 34, 1946.
[28] Kanti V Mardia. 9 tests of unvariate and multivariate normality. Handbook of Statistics, 1:279-320 -200, 1980.
[29] Valentina Masarotto, Guido Masarotto, Maintainer Valentina Masarotto, and Panaretos Masarotto. Package 'fdwasserstein', R package version, 2024.
[30] Alessandra Menafoglio, Marco Grasso, Piercesare Secchi, and Bianca Maria Colosimo. Profile monitoring of probability density functions via simplicial functional pea with application to image data. Technometrics, 60(4):497-510, 2018.
[31] Russell Manning Mersereau and Alan Oppenheim. Digital reconstruction of multidimensional signals from their projections. Proceedings of the IEEE, 62(10))1319-1333 1974.
[32] Gaspard Monge. Mémoire sur la théorie des déblais et des remblais. Histoire de l'Académie Royale des Sciences Paris, 1781.
[33] Gans-Georg Müller and Fang Yao. Functional additive models. Journal of the American Statistical Association, 103((84): 1534-1544, 2008.
[34] Ryo Okano and Masaaki Imaizumi. Distribttion-onddiiirbbuuui regression with Wasserstein metrie: Multivariate Gaussian case. arXiv preprint vrXiv:23377.06.77 2023.
[35] Alexander Petersen and Hans-Geon Georg Müller. Functional data analysis for density functions by transformation to a Hilbert space. The Annals of Statistics, 44(1):183-218 2016.
[36] Alexander Petersen and Hans-Georg Müller. Fréchet regression for random objects with Euclidean predictors. The Annals of Statistics, 47:691-719, 2019.
[37] Alexander Petersen, Chao Zhang, and Piotr Kokoszka. Modeling probability density functions as data objects. Econometrics and Statistics, 21:159-178, 2022.
[38] Gabriel Peyré and Marco Cuturi. Computational optimal transport: With applications to data science. Foundations and Trends in Machine Learning, 1(5-6):355-600 2019.
[39] Julien Rabin, Gabriel Peyré, Julie Delon, and Mare Bernot Wasserstein barycenter and its application to texture mixing. In Scale Space and Variational Methods in Computer Vision: Third International Conference, SSVM 2011, Ein-Gedi, Israel, May 29 2, 2011, Revised Selected Papers 3, pages 435 446. Springer, 2012.
[40] J. Radon. Über die Bestimmung von Funktionen durch ihre Integralwerte längs gewisser Mannigfaltigkeiten. Akad. Wiss., 69:262-277, 1917.
[41] GR Ramesh, N Srinivasa, and K Rajgopal. An algorithm for computing the discrete radon trans- form with some applications. In Fourth IEEE Region 10 International Conference TENCON, pages 78 1. IEEE, 1989.
12

---

[42] Alessandro Rinaldo and Larry Wasserman. Generalized density clustering. The Annals of Statistics, 38(5):2678-222 2722,, 2010.
[43] Isaac Schoenberg. On certain metric spaces arising from Euclidean spaces by a change of metric and their imbedding in Hilbert space. Annals of Mathematics, pages 787-793, 1937.
[44] Isaac J Schoenberg. Metric spaces and positive definite functions. Transactions of the American Mathematical Society, 44(3):522-553 1938.
[45] Nina Vesseron, Elsa Cazelles, Alice Le Brigant, and Thierry Klein. On the Wasserstein geodesic principal component analysis of probability measures. arXiv preprint arXiv:2506. 04440000004 2025.
[46] Martin J Wainwright, High-dimensional Statistics: A Non-asymptotic Viewpoint, volume 48. Cambridge University Press, 2019.
[47] 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.
[48] Yidong Zhou, Han Chen, Su lao, Poorbita Kundu, Hang Zhou, Satarupa Bhattacharjee, Cody Carroll, Yaqing Chen, Xiongtao Dai, Jianing Fan, Alvaro Gajardo, Pantelis Z Hadji- pantelis, Kyunghee Han, Hao Ji, Changbo Zhu, Paromita Dubey, Shu-Chin Lin, Hans- Georg Müller, and Jane-Ling Wang. fdapace: Functional data analysis and empirical dynamics. URL Bet.aaaae--ppeeeeooo///a/ccooooe
R package version 0.6.0, 2024.
[49] Changbo Zhu and Hans-Georg Müller. Autoregressive optimal transport models. Journal of the Royal Statistical Society Series Statistical Methodology, 85(3):1012-1033, 2023.
[50] Changbo Zhu and Hans-Georg Müller. Spherical autoregressive models, with application to distributional and compositional time series. Journal of Econometrics, 2023.
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 abstract and introduction accurately summarize the paper's key contribu- tions and are consistent with the scope and content presented the main text.
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 can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: Please refer Section for limitations
Guidelines:
The answer NA means that the paper has no 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 a separate "Limitations" section 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 how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low 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 of 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 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.
the paper. The authors should use their best
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
14

---

Justification:
We present three main theorems: one in Section of the main paper and two in Appendix All are accompanied by detailed proofs in Section
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 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 short proof sketch 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 it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: Please refer to Section 5f the main paper and Section the Appendix for the choice of tuning parameters.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that appropriate to the research performed.
While NeurlPS does not require releasing code, the conference does require all 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 new 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 instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome deseribe 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 or verifying the results.
5. Open access to data and code
15

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification:
Code is provided along with the submission.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.co cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
//nips //iiss.cc/pill wttp:ccpppppccceeeeeeeeee..c.ett hotaa//gggudee//eeeeabiiiisssss /Co///ooddiisssseeeeeeee ssiiiiiisiyyyyyyy 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. If only subset of experiments are reproducible, they should state which ones are 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 understand the results?
Answer: [Yes]
Justification: Please refer to Section of the main paper and Section E of the Supplement for the estimation details. Code is also provided along with the submission.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper 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 error bars provided in Section of the main paper.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi-
dence intervals, 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 with given experimental conditions).
16

---

The method for calculating the error bars should be explained (closed form formula, call 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 enror 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 range (e.g. negative enror rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification: Please refer to Section we use and execution time.
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, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs 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 or 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 NeurIPS Code of Ethics htttp::///eurpp::::: http://aaaarrs...cppo..c.ccctttt
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 not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. if there is special consid- eration due to laws or regulations 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]
Justification: The paper presents foundational research that is not directly linked to specific applications, and therefore a discussion of societal impacts not applicable.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
17

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g., deployment 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 partieular applications, let alone deployments. However, if there is direct path any negative applications, the authors should point it out. For example, it is legitimate to point that 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 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 of ML).
11. Safeguards
Does the paper describe 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: The paper does not involve the release of models or datasets with high risk misuse. It focuses on theoretical and methodological contributions rather than deployable systems.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have 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 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 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: All external assets used in the paper, including datasets and code, are properly cited with references to the original sources. The licenses and terms of use for these assets are explicitly stated, and their usage complies with the specified terms. Version information and URLs are included where applicable.
Guidelines:
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 used and, possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
18

---

For scraped data from a particular source (e.g., website). the copyright and terms of service of that souree should be provided.
assets are released, the license, copyright information, and terms use in the datasets.paperssiii paperswi thcode. www.com//aaaa datasets package should be provided. For popular datasets,
has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are 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 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 new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/mmodel as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to 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 participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
This paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor 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 as 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 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 an 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.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
19

---

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 we 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
Does the paper describe the usage of LLMs if is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: The core methods and contributions of this research do not involve the use of large language models (LLMs) in any important, original, non-standard way.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to our LLM policy (https://neurd) /ttpr/ippp..c//oooeeeee/e////22/ for what should or should not be described.
20

---

A Additional Assumptions
The univariate transformation & is required to satisfy the following assumptions (L.1)-1L3), where we denote the L and uniform metrics as d2 and doc, respectively.
(LI) Let 01.92 with gi differentiable and | lo < oc. Set
Do 2 max (llgllse: || / 20 || |/ ll
There exists const Co depending only Do such that d2(v(g)) ((g))) & Cod2(g1. 92), dx (v(())) v022)) & Codoo (91.92).
(1.2) Assume E 9 differentiable with I1''l < oo and D1 is a constant bounded below by max (119|]oo, |1// lg/'ll.)) Then (g)) is differentiable and there exists constant C₁ > 0 depending only Dj such that
1+(g) C1, 11 SO..
(L3) Assume 91, 92 continuous functions, where y is differentiable S with 1| co. There exist constants c, C22 O3 > where C2 o1 11p11 S ce ll 110pp (0,8))0s) and O3 (doo (P1, (2)) S (p1:p))) such that
d2 (4)) x (p2)) s C2C,d2(911
where C2 and C3 are increasing their respective arguments.
The multivariate transformation V is required to satisfy the following assumptions (TI)-(T5).
(TI) There exists a constant Ci such that
d2 (f)),, (f2)) S f2), deo (v(()), (f2)) S Cidso(fi for all fi, F.
(T2) For all f E F, (0,)(0,,)) is continuous and has a continuous partial derivative with respect to with constant C2 such that
(f)l C2, [00 s)/ C2. [V) s) (f))(02, C2 101 02112
(T3) An inverse transformation - exists such that - o V()) for all f F.F.
(T4) There exists sequence of approximating inverses 1 T + 0o, and a constant C3 = C3(T) S
CT (x-p) for a constant e, such that
doo ( o (f),, S C3, for all f F,
where C3 is decreasing 0 as 7 co.
(T5) Let (1, Co O + be continuous functions, where Ci has a continuous partial derivative
0ç1 (0, /0s There exist constants c, C2 0 with
C₁ = C. (7, |K// ||0 (0, $) 0s| xo, (G1, (G1, s0) o flkil (2)
such that
do (V, V, (G2) < Cada (51, C2),
where C2 is increasing in each of its four arguments.
Assumptions (TI) and (T2) ensure the continuity and the boundedness of the forward transformation. while assumption (T3) guarantees invertibility on the image space (FF) Note that & (FF is not guaranteed to cover the entire C space and the inverse transformation - is only defined on the image space V (F). A sequence of approximating inverses V is provided in assumption (T4) and (T5) and is used to map 1,2 elements back to F. Note that the regularization parameter r controls the approximation accuracy and the roles of (1 and in assumption (T5) could be reversed.
21

---

B Additional Figures
(29.311
131,331
(37.39)
(39.411
143,451
(47,40)
125
100
demsily
75
0.00100
50
diastolic bload
(51.53)
(55.57]
(39.61)
(63.65)
169.71]
(73,75)
0.00050
125
100-
0.00025
75-
0.00000
50-
50
100
150
20050
iôn
150
2boso
150
20050
1ào 150
zñosó
iôn
:50
zboso
100
150
200
systolic blood pressure
Figure 7: Fitted kernel density surfaces for bivariate blood pressure data for selected age groups.
(29.311
(31.33]
(37.39)
(39.411
[43,45]
(47.49)
125
100
density
75
0.00100
50
0.00075
169.71]
(73,75)
0.00050
diastolk
100
75
0.00000
50
50
100
150
zóos0
16o
150
zbos0
100 150
20050
100 150 záasa iô0
150
zboso
100
150
2ó0
systolic pressure
Figure 8: Reconstructed two-dimensional blood pressure distributions using the proposed RLQD approach, with sliced Wasserstein variance explained of 0.87-00 00.3 from the first component and 0. 0.14 from the second. The wwo-dimensional representations capture most of the variation in the original density space and enable substantial dimension reduction.
gammall-22
gammall.-1
gammal100
gammal.11
gamma1,2
125
100
densily
75
0.00100
50
0.00075
125
gemma2.-2
gamma2.-1
gamma2.0
gemmaa.1
gemma2.2
0.00050
diastolk
100-
0.00025
75-
0.00000
50
100
150
20
100
150
200 50
100
150
200
so
15o
150
200
50
100
150
200
systolic blood pressure
Figure 9: Modes of variation of the first two components by using standard FPCA approach
C Reconstruction of the L2 process
Proposition 1. Assume (DI), (FI) and (T2)-(T5). Let f E F and be continuous function I x S). Then
dsw (J, v o c) = o (Cs Cad2 (V()),)), as oo.
Here, Cs. Co, and W are defined per assumptions (T4)-(T5) and in equation 6
22

---

This bound indicates that by approximating the transformed density within the transformed L space, one can obtain consistent reconstruction of the multivariate density through the regularized inverse map V introduced in (9). where the rate of convergence can be further specified as in the Proposition.
D Radon Log Hazard and Log Quantile Density Transformations
We review here an additional specific univariate density transformation, the log hazard transformation 35].
We adopt the assumption that has been made in Section 2 namely that the support for univariate density function is = [0, 1]. As hazard functions diverge the right endpoint of the support of the distribution,
consider quotient spaces by equating densities that equal
subdomain S [0, Lul. where 1,
for some 0 < u < The log hazard transformation
(S) is expressed follows:
g(s) S (g) log II-
SES,
(20)
where G(s)) So 9(2)d2 a is the cumulative distribution function. The inverse of a continuous function on S is
= exp (p(2)= - So es(w) du) for 2 E [0, 1e]. exp { -So ex(w) du for E [L., 1].
())
Employing the Radon transform R as defined in (2, the Radon log hazard transformation (RLH) L # is
V (f(()) s) S (R(f)(B.))) llog ((,)(O,s)) - Re(s)
(21)
where = R( u)du is the cumulative distribution function along (O))( fixing the first argument at A.
Corollary 2. For C3, C4 as in assumption (T4)-(T5), when using the RLH transform V as per (I it holds that Cs(T) (-(-p) and C4 (T, 11allo, (0, s) doc (2))
T Galle ee as T + DC.
Similarly, for Radon log quantile density (RLQD) transformation, we have
Proposition 2. If (DI), (FI) and (LI)-(L3) hold, assumptions (TI)-(T5) are satisfied for the
mullivariate density transformation V T6) When using the RLQD transformation Vo
as per we have Cz(T) ~ (x-pp and C4 (T, (I,llls 110c, (0, s)/Osll (So (2))
T a2| (6, /)sls e bl (Ci
The following corollary follows directly from Corollary 2 and Theorem
Corollary 3. When assuming exponentially decaying eigenvalues, A ce for some constants c, 0 > O, with K(n) = Lzg log(n)) | - X., holds that
max sup dsw (m a, r), (x =Op =0p (gk5K(n) lai Sag
max ISKSK ( n -1/4
Furthermore, the eigenfinctions correspond to the trigonometric basis, then
sup dgw (ma a, ( (1P (28n-1/4)
(gkgK(n) Saa
E Estimation and Selection of Tuning Parameter
We discuss the empirical estimation in Section now present the population version to provide complete formulation and to facilitate the proofs.
23

---

Suppose we have a sample of independent random densities (fil ~ 8. Writing X; = (fs). empirical sample mean and variance are
s)0,s) = Σ x),,,s), 91,,, 82) = x((0,,)))X((11.>)))) (00h:s))X{(0000-0-0000:000000000
(22)
Mercer's theorem implies that the spectral decomposition of H leads to
H $1, 02, 82) (1) $2),
(23)
where A and ii are empirical eigenvalues and cigenfunctions of the covariance operator with the covariance kernel H.
In terms of the selection of tuning parameter, By selecting the first K components in the expansion (IT) and applying the regularized inverse transformation, the truncated representations can be obtained through
K
f1 K,T)
+
Xu k=l
E
(x),
(24)
with estimated version
K xurin
(x, V.- = 0
(
+*EExsse k=1
(o)).
(25)
simple and often well-working approach for the selection of the tuning parameter K can be based on the fraction of variance explained in the density space. To implement this method, we need to first define the sliced Wasserstein Fréchet mean with respect to the metric dsw, which is
fe arginf E (dsw (dsw(f.q))) (f.
(26)
gEF
with empirical estimates fa = arrgminngg 9E. n- - & sw (fi, q). We note that Sa and f might be
sets with more than one element. Further results regarding the sliced Wasserstein distance and sliced Wasserstein barycenter can be found in [38 71. To quantify the variance explained in the density
space, we use the Fréchet variance
=
Edsw (dsw (Z)))))
(27)
which is unique irrespective of whether the Fréchet mean is unique not. Its empirical estimate is V. Σi dsw(fr. fo) The variance explained by the K component representation is
VK Voo E (dsw(f. )²)
(28)
where fk.- = (v Ek= X47x) is the truncated representation of the density process. This variance explained estimated by Vx V. E. dsw (f.)). (:. K,r))* Control of the tuning parameter is implemented by adjusting the sampling rate in the density domain, which determines the Nyguist frequency- -th highest frequency representable without aliasing. This Nyquist frequency directly sets the cutoff for the regularized inverse Radon transform, thus controlling
the level smoothing in the reconstructed density (55 41
The fraction of variance explained is the ratio VK /Vo.. with estimated ratio V / The choice of K is then the smallest K that explains a fraction of variance above given threshold & with 1, formally
VK K
K argmin
0
(29)
Vo
K
with plug-in estimator
VK > .).
)
R² argmin
K
(30)
Vo
Here the default choices of 9 are 1 .95 or 1 = .9.
24

---

F Random Distribution Estimation
We assume that the multivariate distributions we are dealing with are well-defined and have density functions fixed domains, which might be unknown but encompass fixed known domain D. The domain is as in assumption (D1) and is determined either by prior knowledge or selected by the user in practical scenarios. Our focus is distributions truncated the common domain D and we presume that these truncated distributions are well-defined with density functions on D. Similar to the preliminary density estimation step as in Section S.3 of [9], assume that there exists such that D. = USED (2, <) belongs to the domain of all underlying multivariate distributions, where B2, c) is a ball with radius c centered at 2. This assumption guarantees that boundary effects can be ignored in the density estimation step. The continuous differentiability assumption (F1) is extended to (F1') below. Additional assumptions (P1), (K1), (K2) are also listed below.
(F1') There exists constant Mo 0 and integer 2 p + 1 such that for all I E F. the density function f is the density of a truncated distribution on the domain D, where the original distribution has a domain D, with D, D. and density [D, that is defined on Dy. It holds that Mo on D. and that JD, is continuously differentiable of order k on De. with uniformly bounded partial derivatives. The target distributions with densities are the truncated versions of the fD, on D, so that for all 2 E D one has fzz = fD, (2)/ JD JDD (u) (u)du.
max max(I]f (Ilf,, 11= f11 )
(P1) ((n) min N; 2 en', where N; the number of random observations for the i-th distribution, and the constants satisfy r + p/4 and c>> 0.
(KI) The kernel function K per 31 is probability density function that has compact support and a symmetric, bounded and & times continuously differentiable function (without loss of generality, the support is assumed to be contained in the unit cube of R?).
(K2) For some A and 0, the class of functions Is [K (1K())) RP, 0} satisfies
A
sup M (Io, L2(P), ((IF)L p
where MT, d, E) denotes the e-covering number of the metric space (T, d), F is the envelope function of Ia and the supremum is taken over the set of all probability measures on R The quantities A and w correspond to the
(VC) characteristic K.
Note that assumption (P1) ensures that the L convergence of the multivariate density is faster than the parametric rate n -1/2 while assumptions (K1)- are standard to derive the T.2 and uniform convergence of the multivariate density estimator. Let jAD, be random probability distribution with density function SD, that satisfies (F1'), from which random observations Wi..., WN are independently sampled. A standard kernel estimator f is then applied to estimate fD,, and its truncated version ff is employed to estimate the density f, which is the distribution μDd truncated to the domain D.
1 N
W
f(x)= I(x) / [(w)dd DEDC RP.
f(x)
Σ
KC
(31)
h
The following proposition is direct consequence of Propositions S2 and S3 in
Proposition 3. Assume (DI), (F1') and (K1). Choosing h N pi. the kernel density estimator
in (11) satisfies
d2 (F.1)² =0 -ONNS/0444
sup E /EF
(32)
If we assume further (K2), any 8 > 0.
sup P doc (f.f) N sea sha o(1) as N xo.
N
(33)
fEF
25

---

Proposition iiniccates the convergence rate of the proposed density estimator based on the observation sample size N from the random distribution. The following proposition provides an enhanced form of uniform convergence for the multivariate kernel density estimator under assumption (P1), which is essential for establishing uniform convergence the representation the I,² space.
Proposition 4. Assume (DI), (F1'), ) (KI), (K2) and (P1). For any R > 0,8 > 0,
n. sup FEF
P doc (f.5) RN शंं -
o(1) as - c0.
(34)
G Additional Asymptotic Convergence Analysis
Convergence analysis of the representations of the process is presented as follows.
Theorem 3. Under assumptions (F1), (TI)-(12), (EI) and (KI)-(K2), and with D, H, de, 7 as per and (7). for the k-th component of the eigen-expansion,
d2 (v, D) = Op (n-1/2), do (v, D) = =Op (
(log(n)
1/2
(log() 1/2
d2 (. = Op (n (1-1:2). dx
H, = =0p
n
|>x du = Op (12
d2 (7k,, 916) = xx 0 Op (1-1
doo (7k, 91k) = Op (>=ton n 1-1))
Furthermore, the eigenfunctions correspond to the trigonometric basis, then |/ = 0(1) and [[0Tx(0, (0, /)/0slloc Ok)), so that,
max sup dsw (m=(I, at T), (z, a, ()) = log(Poog -1/4).
(m)) jalgo
The following additional assumptions impose additional restrictions on the eigencomponents. They are satisfied for eigenvalue sequences with exponential decay and are needed to show convergence of the K-representations the target densities in Theorem 4 Both conditions involve sequence of included components K = ((n) with K + oo as n + cc.
(S2)
K
X
n- 1/2 0(1).
(S3)
Σx
K (En>k Ab) = 0(1).
Ax
The following bound is essential,
R(T) max max C4 (T. 1/ 110X, s)/08|loo doo 81t) 1sks 1515
(35)
2d 2oo (X:.Xtt)
0 PP KEkSK.IsiSn
(36)
Here ||x and (0, s)/Os} | are bounded, as per (T2). The convergence of doc Xx,, Xis
can be decomposed into two parts, the convergence of doc (Xits X and that of doo (X,, Xsk)-
Both are (1), using assumptions (S2), respectively, (S3).
The convergence rate of the truncated density estimator in (5 is as follows.
Theorem 4. Under assumptions (DI), (FI), (TI)-(T5), (EI), (KI) if K(n) + satisfies conditions (S2)-(S3), then
max dsw ((,),,((,KTT) = =0, (1( 7 n -1/2 1
K
8
1/2
Σ
A
k=
K>K
26

---

Under exponentially decaying eigenvalues, one can derive convergence rates that depend on the dimensionality the domain and the growth behavior of the regularization parameter т. These rates are obtained for specific choices of K and 7, as stated in Theorem which follows as a corollary of Theorem
H Temperature Data Case Study
Our second data illustration involves analyzing maximum and minimum temperature data obtained from the National Centers for Environmental Information, with raw data conveniently accessible at https://www. http:w/wwwc..... nede.noaa. noaa.gov/ This extensive dataset comprehensively covers weather observations spanning the entirety of the United States. We focus on climate stations positioned at major inter- national airports within the 44 largest U.S. cities, since the quality the data from airport weather stations is superior compared to other general weather stations due to the obvious relevance of weather conditions for the operation aireraft. The temporal scope of investigation spans from January 2020, to December 31, 2020, capturing a comprehensive spectrum of climatic conditions throughout the year 2020. In preliminary data processing, we converted the min/max daily temperature into bivariate vectors (min,range), where range=max-min temperature for each day. This simple linear transformation yields bivariate data that are equivalent the original data but avoids the somewhat awkward triangular support set of the original (min/max) data.
Utilizing kernel density estimation, illustrated in Figure TO with Gaussian density kernel and bandwidth set at 20 for the maximum temperature and 10 for the temperature range, we obtain estimated densities Dimension reduction is achieved by employing only five components to represent the infinite- dimensional space while explaining fraction of variance with respect the sliced Wasserstein distance that exceeds 0.85.
depicted in Figure [1 Here utilize five-dimensional representations.
ATLANTA
BOSTON
CHICAGO
HOUSTON
LAS VEGAS
75
50
temperature range
density
0-
-25
75-
NEW YORK
PORTLAND
RALEIGH
SAN JOSE
SEATTLE
6 -44
50
3e-04
25-
0e 400
0-
-25
100
150
50
100
150
50
100
150
50
100
150
0
50
100
150
maximum temperature
Figure 10: Observed smoothed densities for few randomly selected airport weather stations.
27

---

ATLANTA
BOSTON
CHICAGO
HOUSTON
LAS VEGAS
temperature range
25
density
0-
9e-04
75-
NEW YORK
PORTLAND
RALEIGH
SAN OSE
SEATTLE
6e-04
3e-04
25
0e+00
0-
-25-
50
100
50
100
150
50
150
50
100
150
0
50
100
150
maximum temperature
Figure 11: Fitted density surfaces for a few airports with sliced Wasserstein fraction of variance explained at level 0.85 with five components, where the first component contributes 0,41 and the second 0.18.
The first two transformation modes variation are shown Figure 12 while scores related to eigenfunctions are presented in Figure 13 and 14 The first model of variation reflects both the intensity of variation in the underlying distribution and the average level of temperature range. A lower score of the first eigenfunction indicates reduced variation and lower temperature range, while a higher score corresponds to amplified variation, particularly in maximum temperature, and an overall higher temperature range. The second mode of variation reflects the shape of the joint distribution. For smaller scores it reflects flatter joint distribution, where the temperature range exhibits lower overall variation and a weaker correlation with the maximum temperature, while higher scores are associated with a more positive correlation between the temperature range and maximum temperature and is accompanied by increased overall variation in the temperature range.
gammal, 1,
gammal,-1 -1
gammal 1,0
gammal,1
gamma1.2
75
50
temperature range
density
0-
-25
0.0012
gamma2,-2
gamma2,-1
gammaz,0
gammaz,1
gammaz.2
75-
0.0008
50
0.0004
0.0000
0-
-25
100
150
100 150
100
100
50
100
150
maximum temperature
Figure 12: Transformation modes of variation as The first row corresponds to the first mode of variation and the second row the second mode of variation, for a varying from -2 to 2 as per (17)
13 for bivariate temperature distributions.
28

---

45
PORTLAND
MINNEAPOLIS
PONTIAC
BOSTON
CHICAGO
NEW YORK
Scorel
SACRAMENTO
DENVER
KANSAS
INDIANAPOLS OLIS
PHILADELPHIA
10
Latitud
FRANCISCO
TULSA
NASHVILLE
JOSE
LAS VEGAS
LOS ANGELES
LÖNG SEACH
PHOENIX
HEBER
ALBUQUERQUE
OKLAHONA
DALLAS
MEMPHIS
-10
-
SAN DIEGO
TUCSON
FORT
ARLINGTON
-15
AUSTIN
30
JACKSONVILE
HOUSTON
SAN ANTOVIO
TAMPA
MIANI
-120
-110
100
-70
Longitude
Figure 13: Visualization of the scores for the first eigenfunction.
SEATTLE
45
PORTLAND
KINKEAPOLIS
PONTIAC
BOSTON
CHICAGO
DETRDIT
NEW
score2
40
KANSAS
CBLI
SACRAMENTO
INDIANAPOLIS
LHIAA
Labtude
SAN FRANCISCO
LOUISVILLE
BALTIMORE
FRESNO
TULSA
NASHVILLE
SANJOSE
LAS VEGAS
RALEIGH
LOS ANGELES
PHOENIK
ALBUQUERQUE
OKLAHOMA
WICHITA
DALLAS
KEMPHIS
ATLANTA
4
SAN DIEGO
FORT WORTH
ARLINGTON
TUCSON
EL PASO
AUSTIN
JACKSNNVLLE
30
HOUSTON
SAN ANTONIO
TAMPA
NIAMI
25
-120
-110
-100
-80
-70
Longitude
Figure 14: Visualization of the scores corresponding to the second eigenfunction.
Take San Jose and Portland. Both have positive scores for the second eigenfunction, indicating a higher correlation between temperature range and maximum temperature. Nevertheless, San Jose's variance is smaller, attributable its negative score for the first eigenfunction, in contrast to Portland's positive score On the other side of the spectrum, New York exhibits a negative score for the second eigenfunction, resulting weak association between the temperature range and the maximum temperature, accompanied by lower overall variation in the temperature range. San Jose and Portland's 's increased correlation between temperature range and maximum temperature likely is due to regional factors, such as maritime influence due their proximity to the Pacific Ocean, resulting in temperature amplitudes that are relatively smaller in winter compared to summer. New York's weakened association and flatter distribution also likely linked its specific geography with relatively little variation in the temperature range.
flattened shape of the joint distribution. This pattern reflects
29

---

I
Computation Resources
All data analyses were performed on MacBook Pro with an Apple MI chip and 16GB RAM, running R version 4.3.1 (2023-06-16) under macOS Big Sur. The full analysis for the BLSA dataset and the temperature dataset (presented in the Appendix) each took approximately minutes complete.
Proofs
J.1 Proof of Proposition 3
Proof of Proposition B follows directly from the proofs of Propositions $2 and in 99.
J.2 Proof of Proposition 4
Denote the density SD, on the domain D f f= fD, D. We first show that for any R >> 0,8 0,
n sup P (doo (f.i) RN pår 1 {)) o(1) as N 00. co. fEF
By the triangle inequality,
d (1.1) do (f, do0 (F())
From equation (8.19) in 9 we find
sup sup [(x)) - E Ef(x)) - =ED fEF
By Proposition 9 in [42], there exist constants I, and L2, which depend only on the VC characteristics of KC in (K2), such that for any R > 0 and sequences aN. Li- O(h2).
P (doc Ran 1) s P Pd (f. E. Rax (f) + P do (F Rax(s) L11exp(-L2NN1 R2 2 +1Lh RaN
where / is the indicator function. Let aN N 2/2((p+4)) for any ō 0, and h N sta, Thus
sup P doc (F.i) Rax 1) L1 expl -sp 12 fEF
Under assumption (P1), N(n) n², for constant >> + 4/4, therefore
sup P dox (f. Rax(s) I nLi JEF
28r3 =(1) as cc.
Next, we show the convergence of the truncated density f. where the target is fxx ((x)/S E D. Note that when N is large enough we have
J, flu)du
|D| Mo'
H.Hujduu 2 2M
DD]
whence
doo (f.5) s dxc
f(x)
f(x)
f(x) (e
f(x))
doc fuuddrr So (()ydu
Jo fu))dui
Jo f(u)du
2Mo
2M3
de (f, +
doo (f. (f.f).
D
[t]
follows that for any 0and8 0,
sup P doo (f.f) > RN 220 o(1) as co. feF
30

---

J.3 Proof of Proposition 2
Lemma 1. Under (DI) and (F1), there exists a constant such that,
sup
((ff)(1,s)) 8) ((f)(22, s)1 S A| 02 l2, for all 81, B2 E O
leF
Proof of Lemma ! We first show the result for = 2. Write ((a) = (cos sin a)* E S where [0, 2x] and x x2)" R². Let IIa be the two-dimensional rotation operator that employs the angle A and fuA = f(lls o x). Note that for all 2,
((.fn))(0(a))))
(()))((a+)).2))
Calculating the derivative with respect to ((f) A), - R(f)(((a)) = R(fua ()(a),
(37)
Note that the two-dimensional rotation operator can be represented as rotation matrix,
[cos A (sin A
sin A COS
II
222
22 sin A 2 sin A + I2 COs A
21 COS A
(38)
Then.
fiLA (I1 22) f(x1 cos A 22 sin A Ii sin A + x2 cos A) f(11, 22)
Of (2 222 ) (<) (cos A 1) - x2 sin A)
dzi
05 (s 22) sin A x2 A 1),
Orz
where 2 (1 co) 111 + col211 cos
x2 sin A) and 23 = - (1-co)a3 + col211 sin A 22 cos A)
for cp E (0, Thus,
SU,
-I2 0x1 (111,02)
0s (I1, + σι 0r,
0s
lim A00
(2),
(39)
A
Combining (37 and 33) we have
(21. (f(a), 2) lim
R(f(x1, -2)
da
A
R (f11 (x1, - f(11, (0(a),
lim
A-0
Of
Df
R
32
(2)) (0(a),
Our
Under (DI) and (F1), there exists a constant Ao such that
@R(()) (0(a), 2)
< Ao.
da
Let 0, 0(a1) and 02 ((22) S2, an æ1,02
E [0,2 2T... where we note that |00 - 02112
[2sin((() az) a2)/2) Z laj - azl when |a] 0g
+ 0. Thus, for any 0,, 02 E S , E R, f E F,
- - 2) < Agla; 0gl S 2Aol 02112. as 1 - 02112 + 0.
Next, we provide the proof for p =
Let [0, 2T], (0. m), f(a, 3) (sin cos a, sin 2 sin a, B)T e g²
and x = (21, 22, 23)
E R3 Noting that
OR (f 22, x2a) da
(0(a, 2)
0f an 3YR
(1), 22 22, u3)
(0(a, 8), + (cosarsin 3yR
0f
(01, x2, ma) (((a,)).2)
Ox
31

---

0R 0R (f 22, xa))
(0(a, 3),
03
-sin cos 3)R
0f
(21, 22 m3) (0(a,)),2) 2), + (cos a cos R
0f
(œ1, x2 2s)) (o(a,
Ozi
0f
(sin B)R
(I), (2,33) (),a,,),2))
0xz
under (D1) and (F1), there exists constant Ag such that
@R(f)
@R(f) 03
(0(a, 3), < 2Apl sin B1,
(0(a, B), 2) 33Ag.
da
Let 01 ((1,, 31) and 02 0(a2, 32) g 01,02 [0, 2ml and 31, 32 E 10, Since ||0
021l2 /2- (sin 31 sin 82 cos(a2 01) + cos 31 cos B2) 2 sin 81 8i||a2 (000--a)- ail [32 B1| when
||0 0il12 O, for any 01, 02 8 R. F, as |0 02112 0.
(R(f)(01, ((f)(O2, =)1 - 2A,1 sin 31 || an + 3.A 33 811
< |0 01112-
Following these arguments, it is clear how the proof can extended to the multivariate case when p > 3, and we omit the details.
The analytical tools in Section 2.2 provide the following property for the sliced density function, which implies that the sliced univariate density function adheres to the domain assumption (A1) in
35 This observation is basic for our subsequent exploration of multivariate density transformation.
Lemma 2. Under assumptions (DI) and (F1), (, J E F, is n-differentiable, k. 2 p + 1, and
there exists constant Co such that
@nR(F) (0,
sup
[[(]]l Co. fEF
sup
|1/R(f)llc
Co.
sup JEF
Co.
JEF
Oz
Proof of Lemma From assumption (F1) which states that the density f uniformly bounded above and below, it follows that
R & Mo |1/R(()| S MM
where [D] represents the Lebesgue measure the domain set D. From (S1) in 9 and assumption (F1), there exists constant Co 2 Mo DI such that
0 R(()) O= (0, 2)
S Co..
Note that the constant Co chosen independently of the function f, which completes the proof of the proposition.
Proof of Proposition
Denote the slice-wise density transformation き as
(0,) s) (h(((0)))) h E R(F)
and the slice-wise inverse transformation 5 as
(0), = (C
(EV(F..
(40)
The multivariate density transformation in 66) is v (0,,s) s) To ((,)(O,s). From assumption (D1) and Proposition 2 in [ for any f1. E F, there exists a constant A₁ such that
d2 (R(Ji), R(fx)) A,
d (.ff))R((f)) < A (2)
32

---

Furthermore, from Proposition 2 there exists a constant Ag such that
(R(
| S A2, |1/R(/)| S A2,
S A2,
for f E EF.
Oz
From (LI), for any S1, E there exists a constant A3 such that
(oS JncS (JoR(fi((. u) -8oR(()(0. 1/2
d2 (w R(fi), 30 (
( A3d2 (R(T)))))))
1/2
()d))R((f),R((a)))
Similarly,
dx ( R( Agd ((1)))
which proves the existence of C₁ A1A3 in (TI). Furthermore from (L2), there exists a constant A such that
a 2 R (f) (0, 2)
sup fEF 15 ())R)(()) (1)
S A4,
sup fEF
S
A4-
Oz
From Lemma D and (L1), it follows that there also exists constant A5 such that for any f E F,
sup ( (f)(0..)) (NoR) (f) s)
- Co sup (R(f)( R(S((2, z)1
< Co As| 02|/2-
This shows the existence of C2 in (T2). Assumption (T3) is supported by the inverse transformation formula (40) and the injectivity of the Radon transform on rapidly decreasing functions, as established in Theorem 2.4 of (21]. The regularized inverse 7 is V= R- vh. From Corollary 1 of (9]
there exists constant depending on Cs(r) = O(T p)) as r such that
d 7 (f),)f d (R," : T o - R(f).f)
doo (R=1 R(f),f) S Cs(T),
which shows the existence of C3 in (T4). Furthermore, there exists constant Ao(T) O(7P) as T -> 2o, such that
dos (G)) (Ga)) = do (R. (G), o (
5 Ao(T)d2 (()) -0,,,60-0 (G))
(41)
From assumption (L3), there exist constants
Ar(1IGi 11 |10 s) (0) a0) and As (do (S1,
( / (G(0))) ((= ) 2 1/2
(G1,, (G2)
( JESSP-1
121Rdz (GI(0,) 1/2
Agdd
(42)
Combining 41 and 42 we have
do (V, V (G2) < G4C (2).
33

---

which proves the existence of C4 = As in (T5), From Proposition 2 of [35], when using the LQD transformation t'Qq.
A7 (1001(0..)//85|lceelll)
Ag = O (0a==022a0laa (c2dm(Gtnss).
Thus, when using the RLQD transformation QQ. we have
Ca(T, 11 | (0, s)/0s>///0ssssss (21, -PI s)00sss e (61 For the LH transformation N/r.
A7 = 0 (11l5 As 0 (edau((ris))))
Thus, when using the RLH transformation Nr, we have
CT (0t (x01// (2)) ~ ? (e
J.4
Proof of Proposition| 1
Proof. From assumption (T2),
[f S C2, (2)1 CC2.
From the fact that f and ( is continuous function on (2(() S),
dx (*))) 11 + +| oc.
Thus
C4 (411:///00llle000000000000000110 < oo.
From assumption (T5),
dna ((((()), < Cada (V))),
and from assumption (T4),
dxc (f)) SC3
Then
d (f) C Cada (V))).))
From Assumption (DI), we have
d2 (f- (C)) I C3 + Cide (f), From Lemma S3 in and assumption (D1), follows that
dsw (f, / (4)) O (Ca Cid2 (f).
J.5 Proof Theorem
Lemma 3. Assume (12) and let & per 22 then
d (ū, Op (11)
(43)
Proof of Lemma 3 As S C R is compact interval, one can find grid Ts = {si, .... sa) with [Tsl O()) such that for any s S, there exists a point 8m Ts satisfying |s 5ml 1 /n. Similarly, one can find set Te (0,,.., 0z) such that for any 0 E θ, there exists point 0m. E Te satisfying | 0, < 1/n. The size of the set [Tel is determined by the covering number of O From Lemma 5.7 and Example 5.8 in [46), can choose Te such that its size satisfies
34

---

Tel Onn)). Forming the product of Ts and To leads to the grid T = Te x Ts on Q S with IT| (n
From the smoothness assumption in (12),
sup
[0(0,5) 8) - 1/0, s)
sup 9,
[000, s)
v(0, +
Under (T2), [X: I C2. Applying Bernstein inequality for fixed (0,s) E T, for any
R > O,
log(n) 1/2
n
log(n) 1/2
P
[0(0,)) - - 8) R
P
X((0,8) - -(00,))((0ss) R
71
n
3(n R?
S 2cxp
GnC2 + 2C2nR( log 1/22
Then
P sup (0(0,)) v(0,8) 2 R log(n) 1/2)
s 2n exp (-
3(n n R)
6nC 2C2l2(n log 1/2
B,
2n
(44)
where R (l)g///)()/// For fixed dimension p. choose a large enough such that R
3R
1. Then
[5(0,)) - v(0 = Opp ((1os
|log
1/2
(45)
8EB,ES
Proof of Theorem
Assumption (T2) implies that X 2 Denote by i, G the population mean function and covariance function as per (22) Using Theorems 3.9 and 4.2 8. it follows that
d2(v, D) Op(n (12), 1/2), d2(G, G) (p(n)1/2)
From Proposition 3
sup F dz (J. [| 0 (tr
JEF
The convergence rate of d2(2, is
P((d2(D,)) > RN गंअ }} < P Σ de (X, X,) > RN नM
1 Liz1 ((d Xi))
RN मं
(d2(X,, X.)2)
RN
1 C E (d2(se +
RN
O
which indicates that d2( ONNN then d2(v, D) On 112 - Similarly, have d1, .I) = (n 1/2 + N -f²). In terms of uniform convergence, Lemma provides the uniform convergence for & and similarly for A.
loggn n2
1/2
d:o (H, H) Op ((10) n
1/2
doe (b,v)) =0p Op
35

---

From Proposition 3 for any R > 0,8 > 0,
P ((dxo(2 > RN p+a +6)
max Pf > CT RN 2 - 18i5n
S
(1))
Thus (0) = =0p (
1/2
log
+ N पीス
and similarly,
doa (H, H) = Op (( 1og
=0p
1/2
+ N =1x+
+6
Based on Lemma (4.2) and (4.3) of [8], it follows that
AA - Axl = N
N
(2 7x) = o On (2 1/2 N -
In terms of the uniform convergence rate of 7kk employing the techniques in Lemma of 33 and the proof of Corollary [47].
doo (N )x S doo (AeTiss SxFin) doe (NeTTO du
S doo (11 + d2 (7k 9N) 1111 + Ivs - Su 11
Thus,
doo (7k, 91k) =0p (2
J.6 Proof Theorem
Lemma 4, Under (D1), (FI) and (T2), there exists constant C such that
S SC lnlll I CX (0, C²
||||00 S C,
(46)
(47)
Proof of Lemma
The derivation of 446 follows directly from (T2). Note that
[ jexs
(0, Xx
H s, (a,
(T2) and the fact that |////2 then imply (47
Proof of Theorem 7 Assume K is fixed. Then from Lemma 4
Xx S |1ullo + an Aillaall S sc(1 ao (Vo+-iz))
(48)
Similarly,
(0.)aa (c (00
(49)
Employing Theorem
(Virdo (Yks, 9x) l1 Vxx v Sra)
doc (Xi.as xk(a) & doo (v, + 000
= =ppo -1 \ - =).
Here we use the fact that
D& Sel
VTe Vie
-dol n co.
Vds Vxx
36

---

For the L convergence,
(Vāida(Ta, (w) + VAo - V
d2(X Xkoo) S d2(v, + ao
= =Oplox'n n 11/2).
From assumption (T4),
do (mx U, m) 1)
I Ca(T, aa ||0X (8, doo Xxaa))d((Xk. (k,)).
Case 1: When K is fixed,
max ||0X o (0, s) s)/0s|l do (X (xal =pp().
Writing K)) as per (9). by (T5) follows that S(T, K)) (7T)). Thus
doo (m=(a, T), ma M,T)) = Op +).
From Lemma S3 in 9 and assumption (D1),
dsw (mila,, v, ma v)) Op(T"n b).
Case 2: When K [(n) + co.
doo (Ki ==aa==0pp)
From assumption (T5), 48) and 49
do (m V, r), my (a, T))
S C4 T: oko, 110 (0, s)/0s|l doo (NR,as {k.a) 0x nn
Op T" max 1sksk
(5 110 8) /)/ |o0 ok n 1/2
From Lemma S3 in 9 and assumption (DI),
(7"
} o o n1122)
dsw (my(a, m. v,,)) =Op
max ISkSK
J.7 Proof of Theorem
Lemma 5. Note that under assumption (S3),
d (x,xx))=0p00
(0, s) and let Zx - x -xxx- For the k-th principal component X = Sexs X X(O, 8) (0, s)d0ds < C2 from the assumption (T2). We first show that the process Z is Lipschitz continuous,
Proof of Lemma 5 Recall that (x, =((+)))=v((s)++++ v)
Z (01, $1) Zx (02, s2) S Lx (110 0212 1s2 s2l).
where L2 = =OEE 1) From assumption (T2).
(8(0 - ((22, 82) < (0,, $1) X(2,, s2)l + IX(01, 81) X (02, 82)
Σ C22 (1(0 - 02112 + |s1 - -22)).
Similarly,
(0(01,81) - v(02, s2)1 5 C2 (110, - 02112 + |8| - -2
H (1,,st,, H(0, 82, a. t) < 2C, (110 - 02112 + |s| =22
37

---

Expanding the eigenfunction (x,0,)) = >x Texs H(0, s, a,,,t, t)dadt,
(0x(01)))) ->:: (62, s2)l S xxx / (H(0,, a, H(2,, $2, a, t)) (a t)dadt
888
s 2CC2 (110) 02 + +181- - -22)).
Then
K
(Xx(01,81) - (x2, 82) 5 (1) - -v(02, 82)1 + Σ
1x 17x (01,81) - nn (02, 82)
k=1
K
S2
1+ 102 Σ
(110 - 02112 + +1s1-52)).
1=1
whence
[Zx(O1,81) - Zx(02, 82)1 S C2 2 + 2C2
Σ xxx
(100 021l2 + 1ss 82))
k=1
Lx (110) 02112 |s1 s2
where LK = 0 ). Similar to Lemma 1 in 35 it follows that
12xl S 2 max (10) (191-12|s-~"/ Zx LR 1/3 1| 2/3
(50)
Regarding the 1 norm of Zx.
E (1]x = 1 E(X - -Xx)² d0ds
Se E Σ
E XOTK
d0ds
KSK
=Σu Ak.
k>k
Then
7 11
||X X Op
Σ
Nk
(51)
WSK
Combining (50) and assumption (S3),
1/3
doc (X, =0p
x
Σ
A
-Op().
k
Proof of Theorem 4 Writing X = Σ Xıa 74 and Xix 6 + Ek= & Sr note that
Xik = Jaxs (, (0, (0, (0, sddOds and λia = Sexs (1, s)in s)ddds. By the triangle inequality
and Theorem 3
1/8 exs
Ii Xacl
(X,(0, ())(0))) X,, s) s)) dAds
1/1 Xi(0, 8) (7) (0, 8)) dOds 11800 (X:(0, (x(0, (0,
S dolTk, 11x) + dg (Xi, X,)
-,,00 n- -1/2).
38

---

Here the last inequality follows from Proposition B and assumption (EI). Assumption (S3) and Op(1). The triangle inequality and results Theorem and
Lemma 5 imply dx (X,,,
Lemma imply
K
do Xix S (u, +
do (xa + doo (Xix: X;)
s (v, +
(xxx= Rul + Raddo (nkr (n)) + doo (XiK, Xi)
k=1
=0
Σ k 1
n n-1/2
From (S2) and assumption (T2),
max { X (8)) doo (x,, Xxx] Op(1). 15i5n
From assumption (T2) and (T5),
R(T) = max max C1 T, (1x, do Xix) igksk Isisa
Op(7").
Regarding the 2 convergence.
K
d2 XIx) s d))
+
Σ
d2 (x) + d2(X
S d2(1, +
Σ (Ixik Rial + dan 9x0) d2 (XiK,
K
1/2
=0p
Σ
n n1/22
Σ A
Here use (51) from Lemma and d2 (x(K Xi) = Op (EEnzRde) (T4) and (T5). do (G.,A.KKKW doc (V) (V-((x)) 17' (xix)
Nx) l2) Using assumptions
s do (4 (X.), v (X.)) +
( ( U. ((xxx)
1/2
O,
C3(7) + R((T)
o
1
-1/2 1/2
Σ
NR
k=1
1
1/2
Op (k-p) + 7
Σ
Ok
11/2 +
Σ KSS
d
k==
From Lemma $3 in 99 and assumption (D1),
1/2
dsw v.T) =Op
Op
(-- 7D 7D
ox
n
11/1 +
Σ
AR
KSS
With eigenvalues d ce- for some constants and 0, it follows that
K
K
9k
e 8K
e
km
km
l
39

---

Similarly,
Σ
A
Σ NSc
e e 0k E OK
Choosing ((n) -1
-8 log(n)),
= 0p (T(TT + + 100n
dsw
and choosing 166
dsw (f. (.))) =0p (T-
40


Reviews:
[{"review_id": "3iC0KEVpzV", "reviewer": "Reviewer_X7xr", "summary": "The submission presents an invertible map from the space of multivariate probability distributions to a Hilbert space,  a method to determine modes of variation, and a means of computing these empirically from data. The result builds primarily on the univariate  Hilbert space embedding [33] by Petersen and M'\"uller, and draws significantly on [9] for the multivariate case.\n\nSection 2.2 reviews the Radon transform. Since the Radon transform is unstable, in eq. 4 a low pass filter is applied, and the resulting function renormalized to be a valid density function (eq. 5).  The review is very similar to that in [9].\n\nSection 3.3 introduces the Hilbert Space embedding, which is the univariate embedding of the  Radon transform. Section 4.1 recalls the sliced Wasserstein distance, as similarly reviewed in [9, Section 3], and relates it to the Hilbert space embedding of Section 3.3.\n\nSection 4.2 recalls from [33] the method of obtaining \"modes of variation\" of the Hilbert space embeddings.\n\nIn practice, it is usual to have access to data, and not closed-form densities. Thus in Section 5, approximate densities are obtained via a kernel density estimate on the data, and the modes of variation are computed on the basis of these estimates. \n\nSection 6 describes convergence of the modes of variation, under the estimation strategies and with the assumptions made in the paper.  I did not verify the claims. \n\nSection 7 describes synthetic experiments against functional PCA directly on kernel density estimates, treating them as elements of L2, and comparing with RLQD. The fraction of variance explained is higher for the new method.\n\nSection 8 shows a 2-D blood pressure distribution experiment, where the first two components are sufficient to reconstruct most of the variation, and reveal the nature of changes in blood pressure distribution with age.", "questions": "l. 79: should be l_[z,\\theta} ?  It seems this is an error copied over from [9 p. 5].   In any case not sure why notation l_{z,\\theta} is actually needed in the neurips submission.\n\nThe notation J_1 in eq. 4 is undefined (we are told J in eq. 3 denotes the Fourier transform). I assume this notation is taken from [9 p. 5], i.e. the 1-D Fourier transform?\n\nReferences [36] and [37] are the same.", "limitations": "not applicable", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths\n----\n\nThe paper clearly presents the sliced Wasserstein distance, and demonstrates how to use the Hilbert Space distribution representation implicit in this distance in order to discover modes of variation. The procedure is shown to give a small gain in explaining variance on synthetic data, compared with functional PCA on the density functions. On real data experiments, the resulting modes of variation are interpretable. \n\nThe main contribution of this submission appears to be primarily the calculation of the modes of variation in Section (4.2) their  corresponding empirical estimates in Section 5, the convergence results for these modes in Section 6, and the experimental evaluation of the proposed method. \n\nThe generalization of the modes of variation in [33] to the multivariate representation in [9], and the corresponding empirical estimates, are a sufficient contribution for NeurIPS. \n\nWeaknesses (theory)\n-----\n\nThe abstract makes stronger claims than the paper warrants. First, the claim that the multivariate transformation is \"introduced\" in this submission is overly strong in light of the results of [9] (which the authors cite) - see  the paper summary above for details. Second, \"functional regression\" is claimed as a contribution, however this does not appear in the paper. If it did, then a NeurIPS submission should demonstrate a thorough experimental comparison against relevant ML methods that solve the same problem, such as\n\nZoltán Szabó, Bharath K. Sriperumbudur, Barnabás Póczos, Arthur Gretton\nLearning Theory for Distribution Regression\nJMLR 17(152):1−40, 2016\n\nwhich the authors of [33] are aware of.\n\nThe paper suffers from lossy compression of some aspects of the manuscripts [33,9].   The paper also reads like a longer journal paper rearranged for conference submission. This was not always done with sufficient care, for instance with an entire Appendix H Temperature case study never mentioned in the main document.\n\nWeaknesses (experiment)\n----\n\nFigure 2: performance for RLQD and FPCA seems very similar for cases A,B and K=1,2. Only for K=3 is there a larger difference. Thus the advantage on synthetic data seems marginal.\n\nFigure 3: the argument is made that RLQD is useful as it provides an interpretable functional PCA decomposition of the densities. However, for the blood pressure dataset, would FPCA provide an equally interpretable visualization? On real data is there any significant difference in variance explained between RLQD and FPCA?  It is notable that this is omitted, which weakens the claims made for the method.", "was_revised": true, "final_justification": "Thank you for the reply. I appreciate the additional experimental comparisons with FPCA, which strengthen the claims made. The updated abstract is also now better aligned with the paper.\n\nI agree with reviewer qbBC on the importance of including tangent PCA experiments, and experiments on Gaussian distributions. The authors have committed to conducting these experiments in a final version of their submission.\n\nI will increase my Rating to 5.", "rebuttal": "## **Response to Reviewer X7xr**\nWe appreciate the reviewer’s constructive feedback and have taken several steps to address the concerns raised.\n\n### **Theory-Related Concerns**\n\n- **Abstract revision:** We have rewritten the abstract to more accurately reflect our contributions. Specifically, we no longer claim to “introduce” the multivariate transformation, and instead clarify that the novelty lies in combining Wasserstein slicing with log quantile density transforms in a framework suitable for multivariate distributional data analysis.\n\n- **Functional regression claim:** We agree that the original manuscript may have overstated this. Functional regression is not the focus of this paper, and we have now moved this discussion to the limitations section. The relevant literature (e.g., Szabó et al., 2016) is now appropriately cited.\n\n- **Appendix reference:** We have reviewed the manuscript to ensure that all appendices, including Appendix H (Temperature case study), are explicitly referenced in the main text to avoid confusion.\n\n---\n\n### **Revised Abstract**\n\n*The modeling of samples of distributions is a major challenge since distributions do not form a vector space. While various approaches exist for univariate distributions, including transformations to a Hilbert space, far less is known about the multivariate case. We utilize a transformation approach to map multivariate distributions to a Hilbert space via a Wasserstein slicing method that is invertible. This enables the use of functional data analysis tools—such as functional principal component analysis and modes of variation—with the ability to map back to interpretable distributions. We also provide convergence guarantees for the Hilbert space representations under a broad class of such transforms. The method is illustrated using joint systolic and diastolic blood pressure data.*\n\n---\n\n### **Revised Discussion and Limitations**\n\nThis work presents a transformation-based framework for the functional analysis of multivariate distributions. By combining Radon slicing with univariate log quantile density transforms, we map multivariate densities into a Hilbert space where classical FDA tools such as FPCA and modes of variation become applicable. The regularized inverse transform ensures that the structure of the original density space is preserved. While not the focus of this paper, we note that this framework facilitates distributional regression once vectorized representations are obtained. Limitations include dependence on accurate density estimation from finite samples and the need for regularization due to the fact that the transformed objects lie in a subset of $L^2$.\n\n---\n\n### **Experimental Concerns**\n\nWe thank the reviewer for their constructive comments regarding the experiments. In the updated manuscript, we have expanded the simulation section to include a more comprehensive comparison between RLQD and FPCA across multiple modes.\n\n- **Variance explained across modes:** While RLQD and FPCA perform similarly in lower modes ($K = 1, 2$), the advantage of RLQD becomes more pronounced in higher modes. This is shown in the updated simulations:\n\n  - **Setting 1:**\n    - \\( K = 4 \\): RLQD median = 91.2%, FPCA = 81.7%\n    - \\( K = 5 \\): RLQD median = 95.6%, FPCA = 84.3%\n  - **Setting 2:**\n    - \\( K = 4 \\): RLQD median = 86.6%, FPCA = 75.4%\n    - \\( K = 5 \\): RLQD median = 91.3%, FPCA = 82.7%\n\n- **Practical implications:** In our experiments, the advantage of RLQD becomes particularly evident beyond the first two modes, where FPCA’s explained variance tends to plateau earlier. This highlights RLQD’s strength in modeling higher-order structure more effectively, allowing it to achieve robust performance with fewer retained components. This improvement stems not only from RLQD's better capture of signal in higher-order modes but also from its intrinsic modeling approach that incorporates the structural constraints of distributional data directly into the analysis. We have included this part in the revised manuscript.\n\n  When higher-order eigenfunctions are used, the estimation becomes increasingly unstable. As shown in Appendix G (Theorem 3), the convergence rate deteriorates for higher-order components. Specifically, for the $k$-th component, the estimation error scales as $d_\\infty(\\gamma_k, \\hat{\\gamma}_k) = \\mathcal{O}_p(\\lambda_k^{-1} \\sigma_k^{-1} n^{-1/2})$ and $d_2(\\gamma_k, \\hat{\\gamma}_k) = \\sigma_k^{-1} \\mathcal{O}_p(n^{-1/2})$. \n\n  Importantly, a similar convergence rate also holds for traditional FPCA, which means that the instability in estimating high-order components is not unique to RLQD but is a general statistical limitation of functional eigenanalysis. Therefore, reducing the number of retained components is not only computationally advantageous but also statistically critical for ensuring the stability and interpretability of the estimated eigenfunctions and the final representation. RLQD's ability to achieve high explained variance with fewer components is thus a practical strength, especially when modeling complex or noisy data.\n\n- **Blood pressure dataset:** We have also applied FPCA to the blood pressure dataset and observed two key advantages of RLQD: (i) higher variance explained in fewer components and (ii) more interpretable and data-consistent modes of variation, both of which align with our findings from synthetic experiments.\n\n  In terms of variance explained, the first two RLQD components account for 73% and 14% of the total variance, respectively, while the first two FPCA components explain 70% and 11%. Under a standard threshold of 85% cumulative variance, RLQD requires only 2 components for effective dimension reduction, whereas FPCA requires at least 3. This reflects RLQD's superior efficiency and parsimony in real-world settings.\n\n  In terms of interpretability, the eigenfunctions and modes of variation produced by FPCA typically resemble Gaussian-shaped perturbations around the mean, which can be restrictive when modeling data with inherent skewness or non-Gaussian features. In contrast, RLQD captures more flexible, non-Gaussian modes that better reflect the underlying structure of real-world distributions.\n\n  For example, the blood pressure data exhibit clear asymmetry, with a noticeable skew toward lower pressure levels. This is statistically confirmed by Mardia’s test for multivariate normality: both skewness (3063.39) and kurtosis (42.75) statistics are highly significant (p-values < 0.001), indicating a strong departure from multivariate normality.\n\n  These findings demonstrate that RLQD not only provides more informative components with fewer dimensions but also produces representations that are more faithful to the structural properties of complex real-world distributions. We have produced an additional figure with the modes of variation derived from FPCA, which show systematically less variation in the direction of diastolic blood pressure (to be included in the next version).\n\n---\n\n### **Clarification of Notation and References**\n\nWe thank the reviewer for pointing out these issues:\n\n- **Line 79:** The reviewer is correct — the notation should be \\( l_{z,\\theta} \\) rather than \\( l_{u,\\theta} \\). We have corrected this in the revised version.\n\n- **Equation 4:** We acknowledge that the notation \\( J_1 \\) was undefined. We now explicitly define \\( J_1 \\) in the revised text.\n\n- **Duplicate references [36] and [37]:** Thank you for catching this. We have removed the duplicate and updated the numbering accordingly.\n\n### **References**\n\n- Zoltán Szabó, Bharath K. Sriperumbudur, Barnabás Póczos, and Arthur Gretton. *Learning Theory for Distribution Regression*. Journal of Machine Learning Research, 17(152):1–40, 2016."}, {"review_id": "8carVqx2wH", "reviewer": "Reviewer_moDX", "summary": "This paper proposes a new method to analyze multivariate probability distributions using functional analysis. The key idea is to slice a distribution into many directions using Radon transform and then apply log quantile/hazard transformation to each slice. Such resulting transformation is invertible so it could be mapped back to the original space for interpretations. The authors also prove theoretical guarantees on asy. convergence and did experiments on simulation data to justify the effectiveness of the proposed method.", "questions": "1.  How can be the method applied to high dimensional space if it is difficult to get a good slice from Radon transform? How many samples of slice do you need to make sure the method is working?\n\n2. Can the author explain more on the implications of assumptions made for convergence? I'm not really familiar with this proof and it could be helpful to explain more get me have a better understanding.\n\n3. The asymptoic convergence is based on inf sample size. How it applies to finite sample size? I think this is important as in practice, you have limited data usually.\n\n4. Line 160 mentioned \"data application play an important role\". This looks like disconnected from Section 4.2. Can you explain the important role specifically?\n\n5. I'm not familiar with functional approaches on PCA, but I think the experiment comparison on the standard FPCA and the proposed method seems to be not enough. I asked LLM if there are other methods, it gives me \"Wasserstein geodesic PCA or tangent space PCA\". Can the author explain it to me if it could be possible to extend the experiment sections and if not why the above methods are not available for comparison?", "limitations": "see above", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strength: The proposed method is novel and it allows for invertible mapping that allows classical approaches like pca, regression to be applied meaningfully. This does address the limit of current approach.\n\nWeakness: One weakness is that the proposed method by authors might suffer from high dimensional space, as Radon transform requires sample from an informative direction to get a meaningful slice. The experiments are also too simple. I also feel that the paper mainly focuses on combining two existing approaches for application, although it provides some theoretical analysis on asymptotic convergence, the depth of theory may not fully meet the expectation and emphasis typically preferred by NeurIPS.", "was_revised": true, "final_justification": "Concerns have been addressed in the rebuttal.", "rebuttal": "## **Response to Reviewer moDX**\n\nWe thank the reviewer for the constructive feedback. Below we address each concern in detail:\n\n---\n\n### **Applicability in High-Dimensional Spaces**\n\nOur method is primarily designed for multivariate but *low- to moderate-dimensional* settings (e.g., 2D or 3D distributions), which are common in practical applications such as joint modeling of anthropometric, physiological, or environmental variables. While the Radon transform may face limitations in very high dimensions due to the curse of dimensionality, our convergence analysis (Section 6) explicitly shows how the estimation error scales with the dimension $p$.\n\n---\n\n### **New Approach**\n\nIndeed, the method draws on two existing ideas (LQD and sliced Wasserstein) but introduces new tools—specifically, modes of variation and K-projected representations for distributional data. For both, we develop non-trivial theoretical guarantees with convergence rates.\n\n---\n\n### **Number of Slices**\n\nWe found the method performs well with 50–75 Radon slices per distribution, sampled uniformly over the unit sphere. This balances computational efficiency with numerical stability. Empirically, we observe that beyond this range, the performance gains saturate, while using too few slices can lead to poor approximation and reconstruction quality.\n\n---\n\n### **Implications of Convergence Assumptions**\n\nOur convergence theorem involves multiple assumptions, each serving a distinct purpose:\n\n- **(D1)** and **(F1)** regulate the properties of the input distributions:  \n  - (D1) controls the domain of the distributions  \n  - (F1) imposes smoothness conditions on the underlying density functions.\n\n- **(T1)–(T4)** provide a generalized framework for the class of transformation-based methods considered in our theorem. Our proposed method (RLQD) satisfies all of these assumptions, as shown through a concrete verification. We present the theorem in general form to allow applicability beyond RLQD.\n\n- **(K1)** and **(K2)** are standard regularity conditions on the kernel function used in multivariate kernel density estimation, necessary to establish convergence guarantees.\n\nTwo additional assumptions introduced in Section 6 include:\n\n- **(E1)**  $\n\\lim_{n \\to \\infty} {N(n)}/{n^r} \\geq c \\quad \\text{for} \\quad r > 1 + {p}/{4}\n$  \nThis ensures that the number of within-distribution observations \\\\( N_i \\\\) grows sufficiently fast relative to the number of distributions \\\\( n \\\\) and the domain dimension \\\\( p \\\\).  \nAs shown in Proposition 3 (Appendix F), \\\\( N_i \\\\) directly impacts the accuracy of density estimation. (E1) is critical for uniform convergence when estimating multiple distributions jointly (Proposition 4, Appendix F).\n\n- **(S1)**  $\\sigma_K^{-1} \\lambda_K^{-1} n^{-1/2} = \\mathcal{O}(1)$  \nThe condition governs the estimation stability of higher-order eigenfunctions. It is only required when the number of components \\\\( K \\\\) grows with \\\\( n \\\\). Estimating higher-order modes is increasingly sensitive due to the decay of eigenvalues \\\\( \\lambda_k \\\\) and narrowing spacing \\\\( \\sigma_k \\\\).  \nAs shown in Theorem 3 (Appendix G), we have the convergence bound:\n$\nd_\\infty(\\gamma_k, \\hat{\\gamma}_k) = \\mathcal{O}_p\\left( \\lambda_k^{-1} \\sigma_k^{-1} n^{-1/2} \\right)\n$.\n(S1) ensures this remains bounded as \\\\( n \\to \\infty \\\\). These assumptions are standard in FDA and are shown to be achievable. For instance, under exponential eigenvalue decay (e.g., \\\\( \\lambda_k = c e^{-\\theta k} \\\\)), these conditions are easily satisfied (Corollary 1 and Theorem 2).\n\n---\n\n### **Finite-Sample Interpretation**\n\nOur theoretical results are asymptotic and provide basic guarantees about convergence as the sample size increases. They do not guarantee finite-sample error bounds, which we will acknowledge explicitly in the limitations section.\n\n---\n\n### **Line 160 – Clarification on the Role of Data Application**\n\nWe will omit this sentence in the revised version to avoid confusion.\n\n---\n\n### **Comparison to Wasserstein Geodesic PCA and Tangent Space PCA**\n\nWe appreciate the reviewer’s interest in broader method comparisons.\n\n- **Tangent Space PCA:**  \nThis approach has been investigated for univariate distributions (e.g., Chen 2020), but its extension to multivariate settings is nontrivial due to the lack of a globally invertible mapping from the Wasserstein manifold to its tangent space. In particular, the exponential map is not defined on the entire Hilbert space but only on a convex cone within it. This poses a significant limitation for practical use, as elements lying outside the invertibility cone cannot be mapped back to the Wasserstein space without relying on ad hoc strategies—such as projecting them onto the cone—which can distort the recovered distributions.\n\n- **Wasserstein Geodesic PCA (WGPCA):**  \nWhile theoretically elegant, WGPCA generally lacks a closed-form inverse map and does not guarantee invertibility for general non-Gaussian multivariate distributions (e.g., Vesseron 2025). This limits its utility for reconstructing densities or visualizing interpretable modes of variation. In contrast, our method offers an explicit and invertible transformation to a Hilbert space, enabling both functional analysis and exact mapping back to the distributional domain.\n\nAdditionally, our data application in Section 8 falls in a non-Gaussian regime where we do Mardia's test for multivariate normality; see the response to Reviewer X7xr. This reinforces the advantage of RLQD, which is designed to capture non-Gaussian modes of variation. In contrast, methods like Wasserstein geodesic PCA approaches can only achieve invertibility under a Gaussian assumption. For this reason, we believe it would be inappropriate to emphasize comparisons in purely Gaussian simulation settings when the real-world application of interest is demonstrably non-Gaussian. As such, we focus our comparisons on methods that are both practically applicable and invertible.\n\n---\n\n### **References**\n\n- Yaqing Chen, Zhenhua Lin, and Hans-Georg Müller. *Wasserstein regression*. Journal of the American Statistical Association, 118:869–882, 2023.\n\n- Nina Vesseron, Elsa Cazelles, Alice Le Brigant, and Thierry Klein. *On the Wasserstein geodesic principal component analysis of probability measures*. arXiv preprint arXiv:2506.04480, 2025."}, {"review_id": "tJCdLR8FLI", "reviewer": "Reviewer_kYBw", "summary": "This paper addresses the challenge of analyzing samples of multivariate distributions by proposing a novel, invertible transformation that maps them into a Hilbert space using a Wasserstein slicing approach. This enables the application of functional data analysis tools—such as functional PCA, modes of variation, and regression—while preserving the ability to recover the original distributions. Theoretical convergence results are established, and the method is demonstrated using joint systolic and diastolic blood pressure data.", "questions": "no questions.", "limitations": "yes.", "rating": 5, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 3, "originality": 4, "strengths_and_weaknesses": "This paper extends the univariate log quantile density and log hazard transformation approach to multivariate distributions, enabling the use of functional data analysis techniques—such as FPCA, functional modes of variation, and regression—within a Hilbert space framework. The original densities are recovered via a regularized inverse map, and convergence results are established under regularity conditions. The paper also addresses practical challenges related to estimating multivariate distributions from sampled data. \n\nHere are some comments based on my reading. \n\nFirst, how to handle density functions with multiple isolated support sets? \n\nSecond, another limitation of the Radon transformation is the curse of dimensionality.", "was_revised": false, "final_justification": "", "rebuttal": "## **Response to Reviewer kYBw**\nWe thank the reviewer for the insightful comments. Below we provide detailed responses to the main points raised.\n\n\n### **Handling Densities with Multiple Isolated Support Sets**\n\nWhile the Radon transform remains well-defined for general multivariate densities—including those with non-connected support sets—certain univariate transformations (e.g., log quantile density transformation or log hazard transformation) require that the densities are positive on the domain. This implies that the corresponding probability measures must also be positive on the entire domain, and thus separated subdomains are excluded for the theoretical guarantees to hold.\n\nHowever, in practical applications, when dealing with densities that exhibit multiple isolated support sets, we adopt a regularization strategy by modeling a perturbed version of the true density:\n$$\ng = \\frac{f + \\varepsilon}{1 + \\varepsilon \\cdot |\\mathcal{D}|},\n$$\nwhere \\\\( \\varepsilon > 0 \\\\) is a small constant and \\\\( \\mathcal{D} \\\\) denotes the domain that covers the full support of \\\\( f \\\\). This modification ensures that \\\\( g \\\\) is strictly positive and smooth over \\\\( \\mathcal{D} \\\\), which enables stable transformation and reconstruction. The original density \\\\( f \\\\) can then be approximated either by taking the limit \\\\( \\varepsilon \\to 0 \\\\), or by applying the inverse adjustment:\n$$\nf = g + \\varepsilon \\cdot g \\cdot |\\mathcal{D}| - \\varepsilon.\n$$\nAn alternative practical strategy is to model each disconnected component of the density separately, if prior knowledge of the support structure is available. The trade-off, however, is that this approach requires performing separate PCAs for each component, rather than learning a unified set of principal components across the entire distribution.\n\n---\n\n### **On the Curse of Dimensionality and Practical Focus**\n\nIndeed, the curse of dimensionality presents a challenge for multivariate distributional data analysis with the Radon transform. Our method is theoretically applicable in arbitrary dimensions \\\\( p \\\\), where assumption (E1) clarifies how the number of within-distribution samples must grow with \\\\( p \\\\) and \\\\( n \\\\).\n\nThe convergence rate explicitly reflects the dependence on \\\\( p \\\\), and deteriorates quickly as \\\\( p \\\\) increases. Our practical focus is on low-dimensional settings (\\\\( p = 2 \\\\) or \\\\( 3 \\\\)), which are most commonly encountered in real-world applications (e.g., physiological measurements, environmental variables).\n\nOur theoretical framework still provides general convergence guarantees in the multivariate setting. We will include further discussion of this limitation in the revised Limitations section."}, {"review_id": "M5hFaM6rrM", "reviewer": "Reviewer_qbBC", "summary": "The authors propose an embeddings from the space of multivariate probability distributions into an Hilbert space to enable principal component analysis in functional spaces. While several embeddings of this type exist, the authors emphasize that their proposed embedding is invertible, and therefore allows for the interpretation of modes of variations directly in the space of probability distributions. They also provide the asymptotic convergence of the plug-in empirical estimators of the principal components (i.e. modes of variations), and of the reconstructed density estimators. Their method is illustrated with toy examples and experiments on real-world datasets.", "questions": "- Why is it necessary to assume that the support domain $D$ is convex?\n\n- In standard PCA, the explained variance curve typically exhibits an elbow. Did you observe a similar behavior in your experiments?\n\n- I may be missing something, but it is unclear to me whether the embedding proposed in this work is truly original. In what way does it differ from the embedding introduced in [9] for regression purposes?\n\n- Why use the log quantile density transformation instead of the quantile function directly?\n\n- In the Baltimore Longitudinal Study of Aging (BLSA) experiment, one could have approached each point cloud with a Gaussian, and then apply Wasserstein tangent PCA on the estimated covariance instead. How different would have been the results?\n\n- I understand that the family of modes of variation is parametrized by a variable $\\alpha$ within a compact interval $[-\\alpha_0,\\alpha_0]$, which is needed for the proofs. However, what happens as $\\alpha_0$ goes to infinity?\n\n- Typos. Line 183 \"convergene\". Line 503 : Theorem 6 and Corollary 7 do not exist.", "limitations": "Yes.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "The main advantage on using the proposed embedding for principal component analysis of probability distributions lies in its invertibility, which is particularly valuable for representing the modes of variation of PCA in the space of probability distributions. Additionally, the rate of convergence of the estimated truncated density is particularly relevant for assessing a reconstruction error of the method.\n\nMy major concern with this paper lies in its comparison to the existing literature. First, in my opinion, the comparison setting is not fair. Indeed, it is well known that performing PCA on densities without accounting for the density constraints leads to poor results. A more relevant baseline would be the classical tangent PCA approach, which consists in mapping the probability measures to the tangent space at the Wasserstein Fréchet mean denoted by $\\bar{f}$ (the reference distribution) using the Monge maps (i.e., optimal transport maps from $\\bar{f}$ to $f_i$), and then performing functional PCA in the space $L^2(\\bar{f})$. This is a standard technique in the analysis of the probability density spade endowed with the Wasserstein distance (see, e.g., [a, b]). Another fair baseline would be to apply a classical kernel mean embedding followed by PCA. Even without the invertibility property, I believe it remains worthwhile to compare the proposed method with these (linear Wasserstein PCA and kernel mean embedding PCA) more robust approaches.\n\nSecondly, a more in-depth analysis of the Gaussian case would have been valuable, especially given the existing literature on this setting. Although the proposed method assumes compact support of the measures, it is still possible to consider truncated Gaussians for comparison. For example, performing PCA directly on the parameters space (mean and covariance) could offer a fair baseline (e.g. with tangent PCA in the Bures Wasserstein space). In the experiments of Setting A, it could have been interesting to visualize the reconstructed $k$-th transformation mode of variation to assess whether the method recovers Gaussian-like behavior. In particular, in my opinion, the importance of embedding invertibility is not adequately emphasized in the experiments.\n\n[a] Wei Wang, Dejan Slepčev, Saurav Basu, John A. Ozolek, and Gustavo K. Rohde. A linear optimal transportation framework for quantifying and visualizing variations in sets of images. International Journal of Computer Vision, 101:254–269, 2013.\n\n[b] Seguy, Vivien, and Marco Cuturi. \"Principal geodesic analysis for probability measures under the optimal transport metric.\" Advances in Neural Information Processing Systems 28 (2015).\n\n\n\n------\n\nEDIT : I have increased my rating from 3 to 4 in light of the authors' response to their contributions, which I had certainly underestimated. I also appreciate that the authors plan to add a comparison with tangent PCA and extend the experiment with Gaussian distributions so that readers can better understand the geometry  of the proposed embedding.", "was_revised": true, "final_justification": "The authors propose using an embedding of multivariate a.c. probability measures in a Hilbert space to perform PCA on the measures, by taking advantage of the invertibility of the embedding. Although the use of this invertible mapping is very interesting, I feel that it has not been properly justified by experiments and compared to existing methods. Indeed, it is known that functional PCA will fail when considering probability distributions. In practice, tangent PCA generally gives a \"good\" approximation, even in \"ill-posed cases\" (e.g. discrete measures). Thus, in the setting of this paper, namely a.c. distributions supported on a compact convex set, I expect the results to be satisfactory. In my opinion, the paper would strongly benefit from being compared to tangent PCA, as it is unclear if the proposed method is actually \"better\" than existing methods. Including experiments on the geometry of the proposed embeddded PCA for Gaussian distributions is also important in my view. Subject to these two points, I recommend acceptance of the article.", "rebuttal": "## **Response to Reviewer qbBC**\nWe thank the reviewer for the insightful comments. We address each concern below:\n\n---\n\n### **Fairness of Baseline Comparisons**\n\nWe agree that accounting for density constraints is crucial in distributional PCA and appreciate the suggestion to include Tangent PCA and kernel mean embedding PCA. Our current comparison focused on invertible methods for which reconstruction and visualization of modes of variation are well defined. \n\nIn contrast, Tangent PCA and kernel mean embedding PCA are known to be non-invertible, and Tangent PCA with general Wasserstein transports is computationally very expensive in dimensions \\\\( p > 1 \\\\). It is also heavily dependent on the choice of tuning parameters for regularization (e.g., entropic regularization), where the regularized version can introduce distortion. For example, the modes of variation for large factors may easily venture into the non-invertible area of the linearized space, making invertibility a serious issue.\n\nDue to time constraints in the rebuttal period, we were unable to integrate these additional methods. However, we will include them in the revised version to the extent that well-developed code is available. We will also expand our discussion of these alternative approaches and add references [a] and [b].\n\n---\n\n### **On the Convexity Assumption for the Support Domain**\n\nConvexity of the support domain \\\\( \\mathcal{D} \\\\) is not strictly necessary for the Radon transform itself, which remains well-defined for general multivariate distributions. However, convexity is required due to our use of univariate linearizing (and invertible) transformations for each slice—such as the log quantile density transform (LQD)—which require the density to be positive on the entire support. \n\nThis condition is mainly relevant for the theoretical guarantees. In practical implementations, regularization (as discussed in our response to Reviewer kYBw) can be used to address such challenges.\n\n---\n\n### **Explained Variance and Elbow Behavior**\n\nWe do observe an “elbow” in the cumulative variance explained curve, similar to standard PCA. This pattern is evident in both our synthetic and real-world experiments. For instance, in the blood pressure dataset, the first two components explain 87% of the total variance, after which the contribution of additional components diminishes. We will make this observation more explicit in the revised version.\n\n---\n\n### **Originality of the Embedding**\n\nWhile our embedding builds on the slicing approach from [9], there are important distinctions. \n\n- Our framework generalizes to multivariate distributions using Radon transforms combined with univariate LQD, while [9] focused on univariate regression.\n- Our goal is to develop interpretable modes of variation and convergence guarantees for FPCA in this transformed Hilbert space, supported by theoretical analysis.\n\nThe embedding, theoretical development, and application scope in our work are thus quite distinct from those in [9].\n\n---\n\n### **Why Use Log Quantile Density (LQD) Instead of the Quantile Function Directly?**\n\nThe set of quantile functions is not closed under algebraic operations such as addition or subtraction, whereas log quantile density functions (LQD) are. These operations are crucial for our framework to define principal components and variation modes; see [33] for details.\n\n---\n\n### **Alternative Approach for BLSA Using Gaussian Approximation**\n\nWe appreciate this suggestion. However, we show that the blood pressure data are not Gaussian, as expected. See our response to Reviewer X7xr for statistical evidence supporting this.\n\n---\n\n### **Interpretation as the Parameter \\\\( \\alpha \\to \\infty \\\\)**\n\nIn our implementation, the modes of variation are defined over a compact interval, typically \\\\( \\alpha \\in [-2, 2] \\\\)(or another compact interval), to ensure both interpretability and boundedness. \n\nWhile pointwise convergence holds for any \\\\( \\alpha \\\\), uniform convergence can only be established over compact sets—specifically intervals \\\\( [-\\alpha_0, \\alpha_0] \\\\) for arbitrary \\\\( \\alpha_0 > 0 \\\\). This is consistent with standard practices in functional data analysis.\n\n---\n\n### **Typos and Reference Errors**\n\nThank you for catching these:\n\n- Line 183: Corrected “convergene” to “convergence”.\n- Line 503: Removed erroneous references to “Theorem 6” and “Corollary 7”.\n\n---\n\n### **References**\n\n[a] Wei Wang, Dejan Slepčev, Saurav Basu, John A. Ozolek, and Gustavo K. Rohde.  \nA linear optimal transportation framework for quantifying and visualizing variations in sets of images.  \n*International Journal of Computer Vision*, 101:254–269, 2013.\n\n[b] Seguy, Vivien, and Marco Cuturi.  \nPrincipal geodesic analysis for probability measures under the optimal transport metric.  \n*Advances in Neural Information Processing Systems*, 28, 2015."}]

Author rebuttal:

Meta-review: This paper presents a Hilbert space embedding approach to uncover modes of variations in distributions of multivariate probability distributions. The considered multivariate probability distributions are assumed to be compactly supported and possessing positive densities, and the embedding relies on applying apply log quantile/hazard transformations to slices resulting from a Radon transform. While other embeddings of this type exist, the authors emphasize that the proposed embedding is invertible. Asymptotic convergence of the plug-in empirical estimators of the principal components (i.e. modes of variations), and of the reconstructed density estimators are provided. Their method is illustrated with toy examples and experiments on real-world datasets. The paper was referred by four reviewers who lean towards acceptance (two weak accept and two accept gradings). The discussion allowed to improve some points (e.g., X7xr noted before discussion that “the abstract makes stronger claims than the paper warrants” and progresses in that respect were achieved and acknowledged).  Yet qbBC stressed that the paper would strongly benefit from being compared to tangent PCA, as it is unclear if the proposed method is actually "better" than existing methods. Including experiments on the geometry of the proposed embedded PCA for Gaussian distributions is also important in my view. Subject to these two points, I recommend acceptance of the article”, which I concur with.
Decision: Accept (poster)