Distributional Autoencoders Know the Score
Authors: Andrej Leban
Venue: NeurIPS 2025 poster ()
Keywords: Generative Models
Unsupervised Learning Methods
Autoencoders
Nonlinear Dimensionality Reduction and Manifold Learning
Score Matching / Score-based Models
Scientific Machine Learning / AI for Science
Primary area: deep_learning
TL;DR: For Distributional Principal Autoencoders we prove a closed-form data-score–geometry identity and that any latent coordinates beyond the data manifold are uninformative.

Abstract: The Distributional Principal Autoencoder (DPA) combines distributionally correct reconstruction with principal-component-like interpretability of the encodings. In this work, we provide exact theoretical guarantees on both fronts. First, we derive a closed-form relation linking each optimal level-set geometry to the data-distribution score. This result explains DPA's empirical ability to disentangle factors of variation of the data, as well as allows the score to be recovered directly from samples. When the data follows the Boltzmann distribution, we demonstrate that this relation yields an approximation of the minimum free-energy path for the Müller–Brown potential in a single fit. Second, we prove that if the data lies on a manifold that can be approximated by the encoder, latent components beyond the manifold dimension are conditionally independent of the data distribution - carrying no additional information - and thus reveal the intrinsic dimension. Together, these results show that a single model can learn the data distribution and its intrinsic dimension with exact guarantees simultaneously, unifying two longstanding goals of unsupervised learning.

Distributional Autoencoders Know the Score
Andrej Leban Department of Statistics, University of Michigan, Ann Arbor, MI, United States l cm edu
Abstract
The Distributional Principal Autoencoder (DPA) combines distributionally correct reconstruction with this work, we provide exact theoretical guarantees on both fronts. First, we derive a closed-form relation linking each optimal level-set geometry to the data-distribution score. This result explains DPA's empirical ability disentangle factors of variation of the data, as well as allows the score to be recovered directly from samples. When the data follows the Boltzmann distribution, we demonstrate that this relation yields an approximation of the minimum free-energy path for the Müller-Brown potential in single fit. Second, we prove that if the data lies on manifold that can be approximated by the encoder, latent components beyond the manifold dimension are conditionally independent of the data distribution - carrying no additional information and thus reveal the intrinsic dimension. Together, these results show that a single model can learn the data distribution and its intrinsic dimension with exact guarantees simultaneously, unifying two longstanding goals of unsupervised learning.
interpretability of the encodings.
Introduction
The Distributional Principal Autoencoder (DPA) is an autoencoder variant recently introduced in Shen and Meinshausen [18] and related to other approaches [19; 20] built on the energy score [10]. Because the decoder is trained to match the conditional law given code, DPA guarantees distributionally correct reconstruction for all points mapped to the same code. The encoder, in turn, minimizes the residual variability given this code. Optimizing several latent widths simultaneously induces a sinccpal-component--ieee ordering while retaining an expressive nonlinear mapping. Although the original paper demonstrated strong performance on high-dimensional data, precise theoretical explanation the behavior remained open.
We close this gap by proving strong properties for both data distribution learning and dimensionaliry reduction aspects which hold simultaneously in DPA two objectives that are almost always odds in unsupervised learning. First, we derive a closed-form score-geoment identity that links the normal components of each optimal level set to the data score V log (a))t( Second, we prove that if the data lie on a manifold that is parameterizable (or best approximable) by the encoder, then coordinates beyond the manifold dimension become conditionally independent of the data, thereby revealing the intrinsic dimension. Thus, extend the analogy to PCA made in the original DPA work by proving that, instead of finding principal linear subspaces, DPA learns nonlineur manifolds shaped by the data density, with clear, testable dimensionality criterion conditional independence. The results hold for any combination of ambient dimension p and latent dimension k; in examples, we focus on lower-dime diverddimesiioall ones in order to build geometric intuition.
"Code available at: gi thub. com cwm/oma/nooo andd eb/Di stri buti i ona
Auttoencoderaacom
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

The overall organization is follows: Section 2 introduces the score- geometry identity and its consequences, Section 3 establishes conditional independence extraneous latents and related results on intrinsic dimensionality. Section 4 illustrates both findings examples, and demonstrates that the same score alignment enables a single-fit recovery of the minimum free-energy path (MFEP) for the Müller-Brown potential-a longstanding benchmark in molecular simulations. This connection underscores the work's relevance to science: encodings that approximate the MFEP can significantly speed up expensive molecular dynamics simulations by identifying reaction pathways directly from samples. Section summarizes prior work relevant to both main results, and Section concludes with broader implications.
2 Optimal encoder level sets align with the data score
2.1 Background and preliminaries
We denote the data by X ~ Panta supported on X. DPA couples deterministic encoder e: 12 R with a stochastic decoder IR* + R Optimal encoder/decoder(() will typically be denoted by The decoder is not trained predict single point; using the energy score, it learns to produce the entire conditional distribution of the data when the encoding value (or "code") is fixed; this distribution is defined the Oracle Reconstructed Distribution (ORD). Thus, points that share the same code form a level set in the data space, and samples from the decoder are draws from the data distribution restricted to that set. This is the sense in which DPA performs distributionally correct reconstruction.
The "principal" part comes from optimizing all dimensions of the encoding (up to some K) simulta- neously. The k-th term of the objective asks how much variability remains if we restrict ourselves to the first i coordinates of e(X); minimizing these terms together orders the coordinates by how much residual variability they remove, analogous to principal components but with a flexible, nonlinear mapping and distributional reconstructions. We formalize these notions below; compact notation summary appears in App. A.
Definition 2.1 (Oracle reconstructed distribution ORD, Definition 1 in [18])
For given encoder e(-) Rp - R2 and given sample x E R ~ PAota the oracle reconstructed distribution, denoted by Pr.., defined as the conditional distribution of X, given that its embedding e(X) matches the embedding e(x) of
(X|e(X) (((X))e(x))) e(xx)) P
In other words, the ORD is the distribution of the data on the level set of the encoder.
Definition 2.2 (DPA encoder optimization objective, Eq. 4 in [18])
The optimal DPA encoder seeks to minimize the expected variability in the induced orucle recon- structed distribution:
argmin Ex~ Paata
B Byy. x,x
[Ir-r10
(2)
with 3 a hyperparameter, and the norm taken to be the Euclidean norm in TR.
We will denote the level set of an optimal encoder
Te-(() A {u: c = eX){}},
(3)
2.2 Results
Assumption 2.3 (Global Assumptions)
We assume the following throughout this section:
1) The data density Paota is C² and integrable.
2) The encoder is C. and Lipschitz.
The following assumption a local one: it only required where noted, and in case does not hold, the statements that require it are silent for the violating set.

---

Assumption 2.4 (Local Rank)
The (optimal) encoder's Jacobian matrix evaluated at y D (y) has full row rank for almost every y on (X).
The first lemma introduces general relationship governing the expected (3thh power of the distance on a level set of an optimal encoder:
Lemma 2.5 (General integral balance for an optimal encoder)
For any B > O, assume Assumptions 2.3 and define:
Ja (U:X,,) 7) = / V.- [D. (v)) lly - w'1l Pa - -xxi) - dy',
S(X,n) fv.
V= [D] Paata where yj is a perturbation (function), and ó is the Dirac delta distribution. Next, define the level-set mass:
n(xi) (=) ex)) dz,
(4)
Z(X) - Paata(2) ((aea(2))(((()))))) e(x)) dz.
(5)
Then, for almost every sample X Pasta whose level set (X) satisfies Assumption 2.4, and any n in the same function class as e, we have;
S(X.n) E
Ey-p:- [J (x:
Z(X) BymWpp ||X 111
(6)
provided the quantities in Eqs. 4 are finite and Z(X) > O.
We can interpret this relation as a balancing act: the terms J3 and S represent "push" from the data density to deform the level set, while the "variance" term on the RHS "pulls" closer as it needs to be minimized on the level set. While the above relationship holds for any 8 > 0 and is ultimately
about the expectations on a level set, a more striking result occurs when we set , = 2, which yields pointwise balance equation exactly relating the encoder geometry to the data score:
Theorem 2.6 (When 3 = 2, the optimal encoder's level sets align with the data score) Fix 3 2 and assume Assumptions 2.3. Then, for almost every sample X ~ Paata whose level set Le (X) satisfies Assumption 2.4, the following balance equation holds for almost every E (x);
2yy c(()) DT. (v) = Sdata(y) (T)
(7)
V(X) Z(X)
- -1ly cc(X)12
where Sdata = V, y log (y is the Stein score, whenever the following quantities: the level-set center-of-mass:
c(X) 2 / Paora (1) -((X) dy,
(8)
Z(X)
and the level-set variance:
v(X) / ((x111
cx)]l2 Paata 8(e(y) - e(X)) dy.
(9)
are finite and Z(X)= > O.
The proofs of the lemma and the theorem are based on deriving balance equations from the first variation of the encoder's optimization objective (Eq. 2) and are presented in Appendix A.2.
Eq. expresses trade-off: the variance minimization objective (Eq. 2) pulls the level set toward its center of mass c(X), whereas the local data geometry, through the score, pushes back. The balancing factor compares global term, V(X)/2
local term, |yy c( 12. For the unprojected radial vector - c(X) there are three regimes: for small radii it points outward; beyond the critical radius /V X) the sign flips and it points inward; the critical shell itself has measure zero level set and is excluded from the theorem. After projection to the normal space via D.T. the level sets align these "outward/inward" directions with the data score a manner that satisfies both the variance minimization and the local distribution constraints.
3

---

The projection to the normal space D. is natural: the encoding value changes only in directions normal to the level set, so optimal level sets must align those normal directions with the change in density (the score). Furthermore, the normal space has dimension k, making explicit the trade-off between dimensionality reduction (from p to k) and how strictly alignment with the score can be enforced (i.e., in l of the p dimensions). As shown in Sec. 3, if the data lie on a k-dimensional manifold, DPA can, in principle, align with the score in all relevant directions while rendering the remaining (p k) coordinates extraneous.
As in Shen and Meinshausen [18], an optimal encoder well defined for 5 = 2, which suffices for Theorem 2.6. The only subtlety concerns the optimal decoder wniqueness: the energy score is proper but not strictly proper
so distinct conditional laws can attain the same risk even if only one equals the ORD. In our experiments we did not observe such degeneracy; the learned encodings approximate the data score (Sec. 4).
The requirement that the Jacobian has full row rank almost everywhere on the level set can be violated in practice if, e.g., the encoder is suffering from "mode collapse" and mapping large regions of the input to exactly the same constant. Another cause might be that the neural network is not sufficiently be the case throughout this work. As
expressive approximate the manifold, which is assumed mentioned, Theorem 2.6 is simply silent on those level sets; in our examples we observe well-behaved level sets.
direct consequence for extrema follows by setting the right-hand side of Eq. 7 to zero:
Corollary 2.7 (Consequence of Theorem 2.6 for extrema of the data distribution)
Adopt all the assumptions of Theorem 2.6. Then, excluding minima where 1ly c(Xil approaches infinity, an optimal encoder's level sets at extrema will either:
1) have the center-of-mass c(X) at the extremum: y" = c(X), or
2) have the vector from the extremum the center-of-mass (y" c(X)) tangent to Le (X).
The proof is located in Appendix A.2.3. Due to the variance minimization objective, it is typically suboptimal for the same level set to cross multiple maxima of the data distribution. If does occur, it will often mean that the level set between them is approximately line segment, c(X) will lie close to the connecting line (due to the density weighting) and (y - c(X)) must be tangent both.
3 On (approximately) parameterizable manifolds, extraneous latents are uninformative
3.1 Background and preliminaries
In this section we relate the intrinsic dimension K of the data manifold to encoder coordinates beyond K. The decoder is taken be stochastic network that takes noise ~ AN ((p-k)) (by default) as input, with L. X, and, by Theorem 1 in Shen and Meinshausen [18]: d" (e"
P
Given this, one can consider multiple possible encoder output dimensions k, with the remaining input dimensions of d being padded by expanding e: d ([e1a(x), ((*+1)))) Pa.er (X). Then, the joint optimization objective across components (Eq. 12 in [18]):
(e",d") arggmin a,d
1 Ex TrEy-Paacpee Paeeaa)))) [[X - Y11
E y. Pa
[IX r1[]
2
(10)
where Wk [0, are (optional) weights. In this section, we will assume uniform weights: p+1 and B E (0, 2), which makes the objective a strictly proper scoring rule [10], and thus the global optimum unique. This is technical necessity for the proofs; however, using 3 = 2 yields comparable empirical results, as demonstrated in Sec. 4.2.
1
We denote the k-th term of the optimization objective 10 as:
Lg(e, d = ExlEx ~ - Pa [[|X Y11
- Ex v.y
r11]. (I)

---

3.2 Results
We will assume the data X Paata lie on K K-dimeessiona manifold, with K p. We first focus on case where the manifold cannot be exactly parameterized via an encoder, but can be approximated in the energy score sense for instance, union of manifolds. Throughout, we assume sufficient network expressivity.
Definition 3.1 (K' parameterrizable manifold)
A K diddessioonll manifold is K' --araaeterriable if can be approximated in the minimum energy score sense in K' dimeinnsss, that is. for an optimal encoder/decoder pair, the K' -t-term in the loss Eq. 10 is globally the smallest among all terms and among all encoder/decoder pairs:
Lx [(ee, ; d''] = min Lyle, d] c.d.k
(12)
Naturally, we have K" 2 K. As K dimensions cannot be fully "captured" by lower-dimensional mapping. the first - reconstruction term the loss could always be reduced by considering K dimensions; at optimum, the second term in Eq. 11 equals the reconstruction term (without the
factor, cf. Prop. in [18]). Thus we have contradiction as that cannot minimal, either, implying
that K' 2 K.
Next, we show that this definition compatible with the usual notion of exact parameterizability:
Proposition 3.2 (An exactly parameterizable manifold is K' xl-parameterizaabl)
Consider the case where the data is located on K X-dimensional manifold C. IRP, which is exactly parameterizable by some function I2 + R* and our encoder function class is expressive enough to realize such a parameterization in its first K components.
functinnnnttion
If an encoder realizing the manifold parameterization e is optimal (together with a suitable optimal decoder), then the manifold is K' -parameterizable with K' K, and
LKcl(e" dd =
The proof can be found in Appendix A.3.2.
Definition 3.3 (K'-best wwr-bb--pppooimmaain encoder)
Suppose the data is supported on K -ddimensional manifold, which is K' parrameterizablee If solution (e", d²) satisfying Eq. 12 is also optimal among all dimension- K' encoders:
K'
(e", E argmin
Σ Lule,d, d].
e,d
k=U
denote it as the K' d'-best-approoximatin encoder (with an accompanying optimal decoder).
Before stating the main result of this Section, we summarize how the pieces fit. Definition 3.1 is notion about the data manifold: "K' parameterizable" means that the energy score term Ly(e, d] attains its global minimum the model: a "r--bes--approximaa globally optimal among all K'- latent models. Theorem 3.4 below then describes the properties of such encoders, assuming they exist.
K' for some encoder. Def. 3.3 conversely, a statement about encoder (i) achieves that globally minimal K" -term and (ii) is
K' si-bbestapprroximatin encoder are uninformative) Suppose the data is supported on a K K-dimensional manifold, which is K' -prareetteeiiallee Then the K" aprooiappo encoder is also the optimal solution when optimizing across all p dimensions, with the dimensions (K' + 11...., ,P) obeying:
Theorem 3.4 (Extraneous latents
Pa oT. (X) = Pa (X), Wk E [K' + 1...,,
(13)
Furthermore, the dimensions (K" + 1, 1.... ......,))) p) of the encoder will be conditionally independent of X, given the relevant components (e ... )
X L (X) (X), Vi E [l...., K'). In other words, they will carry no additional information about the data distribution:
(14)
I (X: eis (X)) = 0, vi E [l,...,P K'), where [(:-)) is the mutual information.
5

---

The proof is located in Appendix A.3.3 and is based on the global optimality of the K'-best- approximating encoder and the fact that the random variables ei (X) form filtration with increasing k.
The K"- best-approximmtting encoder simply repeats the K' K'-dimensional distributional approximation in the extraneous dimensions, The "extra" coordinates may be deterministic functions of the preceding ones (or the data) or stochastic, yet conditionally independent of X. This result naturally extends the property of PCA, which finds the principal linear subspace of the data, to a nonlinear encoding.
The following Corollary 3.5 is special case for when encoder exactly parameterizes K K-manifold and is globally optimal: it is sufficient condition for Thm. 3.4, and reduces the existence question to the encoder being globally optimal.
Corollary 3.5 (Optimal encoders that exactly parameterize the manifold output the data distri- bution)
Consider the case of an exactly parameterizable manifold with parameterization ei:K- If this encoder optimal among K K-dimenssonal encoders, then it is a K beest-approiimatin encoder and Theo- rem 3.4 holds with K' = K. Furthermore, together with an accompanying optimal decoder, it will output exactly the data distribution using the first K dimensions:
d' (e1:
R.K.
X Panta-
(15)
The proof is combined with that of Prop. 3.2 and can be found in Appendix A.3.2. The capacity of an autoencoder architecture to parameterize sufficiently "nice" manifolds was recently demonstrated in [5]; while the same remains as future work for the DPA, is certainly plausible.
Finally, the following remark addresses the attainability of the global optimality condition for the parameterizing encoder in Cor. 3.5:
Remark 3.6
Consider the case of exactly parameterizable K-dimensional manifold with parameterization eiKK. Then, typically when p >> K, such an encoder is an optimal encoder (for K and dimensions), with the results 3.4 and 3.5 holding.
more precise statement is given in Appendix A.3.4; a short sketch that since the parameterizing encoder achieves zero loss on the extraneous dimensions and cannot do "too bad" the relevant ones, as p > this naturally leads to it being the globally optimal encoder.
Experiments
So far, the results presented are at the population level. In all the experiments below, the decoder is an Engression network, for which non-asymptotic simite-sample error bounds are available- see Thm. in [19], so deviations from the population predictions shrink with the sample size. The encoder is standard MLP throughout, so usual generalization arguments apply.
4.1 Level Set Score Alignment
We present score alignment results for examples where the data density Pdata is known. In order of increasing complexity, we consider a standard multivariate Normal, a Gaussian mixture, and the Müller-Brown potential, which is standard benchmark in molecular simulations.
The first two examples, shown in Fig. 1, train DPA 10,000 samples from either standard Normal (a) or a three-component Gaussian mixture (b). Since the intrinsic dimension K = 2 = we visualize both encoder coordinates. Each panel shows heatmap of the latent together with selected encoder level sets as contours. To accompany Fig. 1, Tab. reports the absolute cosine similarity between the two cormal-space vectors in Eq. 7 the level-set vector and the data score after projection evaluated on a 100 x 100 grid and restricted to points with density > 0.5% of the maximum (further details can be found in App. B.2). The alignment essentially perfect in both datasets. For the standard normal (Fig. 1 a)), the level-set centers c(X) should roughly coincide with the mean of the data distribution due to their geometry. Because the standard Gaussian is rotationally symmetric, DPA can choose any orthogonal pair of encoder directions without changing the loss. As the two components are orthogonal, we can examine each component's 's alignment with the score independently, each should approximately satisfy Theorem 2.6. The first component thus roughly

---

(a)
(b)
Figure 1:
Gaussian examples. a) standard Normal: b) Gaussian mixture, Red contours: data density; black arrows: score. Left: first latent; right: second.
Table 1: Absolute cosine similarity between the normal-space vectors in Eq. 7 (points with density > 0.5% of its maximum kept).
Latent comp.
Mean abs. cosine
Points kept
Dataset
Std
-ile
Standard Normal
0
1.00
0.00
1.00
5088
Standard Normal
1
1.00
0.00
1.00
5088
Gaussian Mixture
0
1.00
3.1
10
1.00
4
Gaussian Mixture
1
1.00
3.0 10
1.00
4729
encodes the polar angle and the second the radius, both parameterized with increasing value the latent
The first component attempts to minimize both sides of Eq. 7, that is, have both y c(X) and Vy log Paota(U) (y) lie the tangent space of the level set. The second component, on the other hand, attempts maximize both sides Eq. by having y - c(X) and V log (y) almost entirely normal to the level set. Thus, we obtain (approximately) level sets that are orthogonal to the score gradient. The encoding also exhibits other desirable properties, such all the level being connected and smooth, directed variation of the latent.
For the Gaussian mixture with non-symmetric centers (Fig. 1 b). the polar picture persists but without rotational symmetry. The region of high density of the level sets coincides with the near-zero gradient region of the data distribution where the contributions of the three components cancel out. As the distribution is no longer symmetric, the score alignment must be understood jointly across latents rather than per-coordinate; nevertheless, the numeric alignment remains near-perfect.
The Müller-Brown potential [16] is a standard ww--dimensional benchmark in computational chem- istry. It features three minima and several saddle points; the minima are marked by red dots in Fig. 2, with potential contours overlaid in red. Physically. the minima represent metastable states, and the potential is designed simulate chemical processes where molecule undergoes (potentially rare) transitions between different configurations (states). The data distribution is the Boltzmann distribution for the Müller -Boown Potential U (x)) at (known) temperature T:
P Panta (z; T) exp( (S((-U(z)/nn)),
7

---

where kg is the Boltzmann constant. We generate samples by running Brownian dynamics (typically initialized at the minima). After re-arranging Eq. we obtain for each level set defined by X:
y - -
F(y) D = - Uw)) (2) TT?
De- (v),
(16)
V(X) 2(X) - 1ly c(x)11
where Fly) the force (field) at position y. This implies that the level set geometry is determined by the normal components of the force field and could, in principle, be recovered from the latter. For training DPA, we discard trajectory information and consider the data as i.id. samples from an unknown distribution. This contrasts with time-series approaches such VAMPnets [15], which explicitly assume a Markovian structure the process. Conversely, in existing i.i.d. autoencoder approaches, such as [7; 4J. the authors commonly use iterative procedure to find good encoding: first, encode the data produced by an unbiased simulation (such as the data used here), then use the resulting encoding to add bias to the potential and run another simulation, then encode the data again, and so on until convergence. Both steps of this procedure are potentially computationally expensive for larger problems.
2.0
2.
15
372
-032
-26
-11
13
(a)
con
18
155
Figure 2: Müller-Brown potential: encoder level sets and comparisons. Red contours: potential; black arrows: potential gradient; purple: the MFEP. (a) First (left) and second (right) DPA components. (b) First component for an Autoencoder (left) and VAE (right).
A quantity of considerable interest is the minimum free-energy path (MFEP), shown in purple in Fig. 2 and computed by the string method [8]. The MFEP connects minima and is everywhere tangent to VU, representing the fleast-eneryy-one ccsstly" transition between the former. In contrast to iterative pipelines, DPA approximates the MFEP from single fit of unbiased data even though samples between minima very scarce (cf. App. B. By Eq. 16, the normal components of the encoder's level sets align with the force field; thus, the first DPA component provides a monotone (approximate) parameterization of the MFEP in its latent z. Along the inter-minimun corridor, its level sets are (approximately) tangential to VU and intersect the high-density regions orthogonally, placing c(X) near the MFEP; following V, exz) from the top-left minimum with increasing thus traces the MFEP. For the second component, the level sets are approximately orthogonal to the gradient in the regions between the minima, and the level set (e(x) R 0} approximately traces out the MFEP. In Fig. 2 (b), the Autoencoder captures the general direction of the transition but does not parameterize the MFEP, whereas the VAE fails encode the transition pathway entirely.
We quantify the MFEP parameterization ability in Tab. while extending the comparison to 3- 3-VAE [11] and 3-TCVAE [6]. DPA achieves the best MFEP agreement across all distance metrics, and its first component is always the best parameterizing one (ž ), in line with the principul-components nature of the method (details and further results in App. B.3).

---

Table Evaluation of distances between the closest parameterization and the true MFEP: each entry is the mean + s.d. over the 24 seeds kept after discarding the worst- Chamfer run for every model. da represent average closest distances for both directions, while 95th -erc. error the larger the two 95th-percentiles for the two closest distance sets. Lower numbers are better for every distance
metric.
Param. comp.
Chamfer dc
Hausdorff du
Mean 95 -perc. error
Model
dhaFEr-ppth
OPUhhMMFEP
AE
0.54 + 0.51
0.387 = 0.113
0.804 + 0.142
0.760 0.110
0.467 + 0.128
0.306 = 0.100
DPA
0.00 + 0.00
0.262 = 0.053
0.730 + 0.317
0.212 0.567
0.289 + 0.066
0.236 + 0.047
VAE
0.62 0.49
0,515 = 0.469
1.461 + 0.973
1.311
0,980
0.541 0.217
0.490 = 0.809
8-VAE
+ 0.51
0.450 + 0.288
1.172
+ 0.512
1.051
0.477
0.539 + 0.180
0.360 + 0.462
3-TCVAE
0.375 + 0.49
0.377 = 0.077
1.378 + 0.501
1,228 0.433
0.591 + 0.180
0.164 = 0.101
Table 3: Determinism diagnostics for extraneous latents U. "ID-drop" reports the 2.5%, 50%, 97.5% bootstrap quantiles (200 resamples). Dataset details in App. B.4.
Dataset
R²
ID-drop
H(U Z) [nats]
Gaussian line
0.9997
0,0112 / 0.0122 / 0.0132
-7.259
Parabola
0,9997
0.0046 / 0.0048 0,0051
-9.190
Exponential
0,9996
0.0045 0.0047 0.0050
-9.162
Helix slice
0.9995
0.0041 / 0.0043 0.0045
-6.090
Grid sum
0.9986
0.0023 / 0.0029 / 0.0034
-2.759
S-curve
0.9996
-0.0031
-..000000///440 0.0000
-1.762
S-curve (3 2)
0,9990
0.0030
-2600
Swiss-roll (3D)
0,9872
0..00.00/ 0.00055 0.0000
-0.042
4.2 Independence of Extraneous Latents
Theorem 3.4 states that, the optimum, latent coordinates beyond the manifold are uninformative given the informative block Z. This can be satisfied by two regimes: (i) on many manifolds, the "extra" latents U become (nearly) deterministic functions of Z, and (ii) they are stochastic, yet conditionally independent of the data X given Z.
Deterministic regime. We test whether U g(Z) using three diagnostics (details in App. B.4): (i) R2 of a nonlinear regressor U on Z (near-l indicates determinism), (ii) the intriniieddmmnsin
drop when using Z vs. (Z,U) via the Levina-Bickel MLE: drop 0 implies no additional degrees of freedom in U, and (ii) an estimate of the conditional entropy HU Z): very negative values consistent with near-determinism. Across datasets, we observe R² Rs 1,, near-zero ID-drop, and strongly negative H(U Z); see Table 3. Results include B 2 row, which still exhibits uninformative extraneous latents (despite the theoretical requirement / E (0,2)).
Stochastic but uninformative regime. In this case, test for the null hypothesis U 1 X Z directly. For the Gaussian line dataset, a conditional randomization test yields p-values consistent with the null hypothesis distribution Uniform (0, 1]; Malnogorov-Smmmoo test gives D 0.061, a p-value p = 0.822 with 4% of 100 replications below 0.05, confirming conditional independence (details and further results in App. B.4).
5 Related work
For denoising (and related contractive) autoencoders, Alain and Bengio [1] derive formula showing that, for each fixed data point, asymptotically the noise approaches zero, the difference between the reconstructed data vector and the original will tend to the score the data. Using our notation:
d" (z))
0² V log ( (z) o (o²)
+ 0, Vr fixed.
(17)
In Bengio et al. [3], the authors further relate the score to the first two local moments of the data distribution (moments restricted to a b-ball around given point) and propose asymptotically valid

---

estimators as both the noise and 6 ++ O. In contrast, the first two moments appear directly in Eq. 7 as level-set moments, yielding a global geometric constraint that holds almost surely. It is also interesting to note that the results presented here are also not connected to any sort of regularization (whether denoising, contractive, or other), as they arise directly from the distributional regression using the energy score.
concurrent line of work underscores the advantage of energy-score training over (denoising) score matching. Shen et al. [20] replace diffusion's 's matre-matching objectíve across noise levels [21] with a sequence of energy-score regressions (as in the DPA), achieving comparable distributional fidelity in far fewer steps. This mirrors the contrast between the asymptotic nature of Alain and Bengio [1] (linked diffusion via Vincent [22]) and our exact, con-asymtootic identities. Mirroring the (implicit) recovery of the Müller- Brown force field presented here, the score matching property of diffusion models has recently been used in [2] to obtain force fields.
On the informativeness of extraneous dimensions, Liu et al. [14] analyze Chart Autoencoders [17] that build on denoising autoencoders and show that, for both exactly and approximately parameterizable manifolds, the optimal reconstruction error decays exponentially in the manifold's intrinsic (not ambient) dimension, provided the corruption acts normal the manifold thus requiring geometric knowledge which DPA does not. For VAžs with learnable decoder variance, Zheng et al. [23] prove that, on simple Riemannian manifolds, the optimal decoder variance collapses to zero and the VAE loss scales with (p - k) log o², making the intrinsic dimension identifiable. Prior work on disentangled representation learning has also investigated the phenomenon of extraneous dimensions (depending on the formulation) being "fturnedoff", typically utilizing extra regularization terms in the loss [11: 12]. In contrast, for the DPA we prove that extra latents are exactly (i.e., not uninformative across both manifold scenarios - without manifold knowledge, additional regularization, or specialized architectures.
Discussion
This paper establishes two exact properties of the Distributional Principal Autoencoder (DPA). First: a closed-form relation between the optimal level sets of the encoder and the score of the data distribution. Second: if the data lie manifold that can be approximated by the encoder, the encoder's components beyond the dimension of the manifold (or its best approximation) are conditionally independent of the data and therefore carry no additional information.
The first identity global and non-asymptotic, significantly strengthening related results for denois- ing/contractive autocncoders [1]. Likewise, the second main result unifies and expands on results from manifold learning and disentangled representation learning literature [14; 23]. Crucially, they hold simultaneously: single DPA fit can be successful in both reconstructionn/aattdiiiit
learning dimension determination via conditional
through score alignment, and disentanglemenn/inriinn independence. This contrasts with almost all unsupervised learning methods where the two properties are a trade-off: e.g., this trade-off is precisely what the titular 3 in B- VAES [11] controls. Taken together, the two properties also yield a PCA-like picture: nested ordering of coordinates and a rigorous cutoff' via conditional independence, with the "principal" directions determined by the data-distribution geometry.
Our statements assume sufficient "niceness" of the data distribution (allowing, c.g. for full-rank Jacobians enough level sets for the Theorem to be relevant) and sufficient encoder expressive- ness to approximate the manifolds. While we derive expressions for general B, the most explicit level-set-score identity uses B = 2, where the optimal decoder need not be unique; empirically did not observe degeneracy. All results are at the population optimum; finite samples and optimization errors can introduce deviations, as discussed with each result.
As demonstrated, score alignment in practice enables (approximate) recovery of force fields from samples drawn from the Boltzmann distribution, yielding a single-fit approximation to minimum free-energy paths and suggesting methods that could significantly speed up molecular simulations. also suggests recovering densities from learned encoder level sets and motivates energy-score- based generative modeling. On the representation side, the conditional independence of extrancous latents also presents a concrete approach to intrinsic dimension determination beyond heuristics such as scree plots: one can directly test for conditional independence in the learned encoding.
10

---

Acknowledgments
We wish thank Nicolai Meinshausen and Felipe Maia Polo for helpful discussions.
References
[1] G. Alain and Bengio. What Regularized Auto-Encoders Learn from the Data-Generating Distribution. Journal of Machine Learning Research, 15(110):34333377, 2014. ISSN 1533- 7928. URL http:htt///cc jmlr .rgg http:/aappssss///eeeeeeeee httpppes//v//aaaaaaa......
[2] M. Arts, V. Satorras, C.-W. Huang, D. Zuegner, M. Federici, C. Clementi, F. Noé, R. Pinsler, and R. v. d. Berg. Two for One: Diffusion Models and Force Fields for Coarse- Grained Molecular Dynamics, Sept. 2023. URL httpt://aaaii wwtpziiioorrg/aaaa HiXiv:2302.00611 [cs].
http://aaxwvv..rraaa///22//2//2/ wttg/aa///////////////// 00600.
[3] Y. Bengio, G. Alain, and S. Rifai. Implicit Density Estimation by Local Moment Match- ing to Sample from Auto-Encoders, 2012. http:/t:://axxi..oooo abs///bs////////// http//ab///110055555/
urXiv:1207.000 [es].
[4] L. Bonati, E. Trizio, A. Rizzi, and M. Parrinello. A unified framework for machine learning collective variables for enhanced sampling simulations: mlcolvar. The Journal of Chemical Physics, (59(1)01,,,,, July 2023. ISSN 0021-9606, 1089-7690. doi: 10. 10.1063/5.0156343
URL https: : /oi. /rggr//// wt..o//.////.................... 1063/5...... 0156343.
[5] J. Braunsmann, M. Rajković, M. Rumpf, and B. Wirth. Convergent autoencoder approximation
of low bending and low distortion manifold embeddings. ESAIM: Mathematical Modelling and Numerical Analysis, 58(1):335-36 Jan. 2024. ISSN 2822-7840, 2804-7214. doi: 10/ 79 URL http:: wwwp///www//.e//////////////t/// esaaeesiinnn m2an.
/mann/ abs /2024/01/ m2an220261 wtw.oo22mm2m2222222222 htmll. Number: Publisher: EDP Sciences.
[6] R. T. Q. Chen, X. Li, R. Grosse, and D. Duvenaud. Isolating Sources of Disentangle- ment in Variational Autoencoders, 2018. URL [es].
04942.
[7] W. Chen and A. L. Ferguson. Molecular enhanced sampling with autoencoders: On-the-
fly collective variable discovery and accelerated free energy landscape exploration. Journal of Computational Chemistry, 39(25):009--102,, Sept. 2018. ISSN 0192-8651, 1096-987X. doi: 10.1002/jcc.25520 URL ht http:///nnnnnnnneeeeee
wwwwyyywww.. com www/dooi.ccccoooooooooooo.o..... com/doi/abs/110 /bb//00
1002/ joc. 25520. _eprint:
[8] W. F, W. Ren, and E. Vanden-Fijnneen. String method for the study of rare events. Phys. Rev. B, 66:052301, Aug 2002. doi: 10.1 10.11033/hyseevv..6..00 or hrgg/doi/11
URL httts:://t//mm
/ wttp///yysssaeeeeeeeeeee 66. 052301.
[9] M. Giaquinta and S. Hildebrandt. Calculus of Variations 1, volume 310 of Grundlehren der mathematischen Wissenschaften. Springer Berlin Heidelberg, Berlin, Heidelberg, 2004. ISBN doi: 10. 0223-662-0222-
978-3-642-080744- 978-3-662-03278--7
URL http:
//link. www.:iprriiggeee. com 1ttp//11..110///0
[10] T. Gneiting and A. E. Raftery. Strictly Proper Scoring Rules, Prediction, and Estimation. Journal of the American Statistical Association, 102(477):359-377 -78, Mar. 2007. ISSN 0162-1459, 1537- 274X. doi: 10/c6758w. URL http::/ttt// wwwp///wwww/com//////////////w wandfannnnnnnneeee wwt.cooooooooooo abs 10. 1. /118///198/ 01621 02244440000001444
[11] I. Higgins, L. Matthey, A. Pal, C. P. Burgess, X. Glorot, M. Botvinick, S. Mohamed, and A. Ler- chner. beta- VAE: Learning Basic Visual Concepts with a Constrained Variational Framework. Nov. 2016. URL https: http://www.eaaattccccclcclcceeee
wwwpp//wpppp..aaaappp/eooo/pa 34 Saciearrinngsseeeeaaaaaaoooooeee www.:eitttaaatiiieaeea-aaeeeee-e Matthhy// a90226c41579 08020220110000007755555577777773
[12] H. Kim and A. Mnih. Disentangling by Factorising, 2018. Utthhttp://arrxwwooorr
1802. 05983. arXiv: 1802.09888 [cs, stat].

---

[13] J. M. Lee. Introduction to Smooth Manifolds, volume 218 of Graduate Texts in Mathematics. Springer New York, New York, NY. 2012. ISBN 978-1-441--9981-8 978-1-4419--9922555 doi: 10.10077978-1-444 1-4419-994--99-9- URL https:://nn....rrngggrr...e
978-1-4419 99822-5.
[14] H. Liu, A. Havrilla, R. Lai, and W. Liao. Deep Nonparametric Estimation of Intrinsic Data Structures by Chart Autoencoders: Generalization Error and Robustness, 2023. URL http:
/ /arxivv.. www.v.org//.o..ee.........oooott abs 2303. 09863. arXiv:2303.0968 [stat].
[15] A. Mardt, L. Pasquali, H. Wu, and Noé. VAMPnets for deep learning of molecular kinetics. Na- fure Communications, 9(1):5, Jan. 2018. ISSN 2041-1723, doi: 10. 10110388
URL https::/// /tt///wwwww nature........ httpe//cco/e.eeeceeesssseeeeeeee 0110000-411166000--0----3---8888 Publisher:
Nature Publishing Group.
[16] K. Müller and L. Brown. Location of saddle points and minimum energy paths by constrained simplex optimization procedure. Theoretica Chimica Acta, 53(1):75-925 1979. ISSN 0040-5744, 1432-2234. doi: 10/bkwf52. URL https: http:///oooooo///./////o/ BF00547608.
[17] S. Schonsheck, J. Chen, and R. Lai. Chart Auto-Ens for Manifold Structured Data, 2019. URL httttttp//aaree org/ abs/1912 10094. arXiv:1912.1009: [cs].
[18] X. Shen and N. Meinshausen. Distributional Principal Autoencoders, Apr. 2024. URL http: aariiv org/ abs 2404 13649. urXiv:2404.13649 [cs, stat].
[19] X. Shen and Meinshausen. Engression: extrapolation through the lens of distributional regression. Journal of the Royal Statistical Society Series B: Statistical Methodology, page qkae108, Nov. 2024. ISSN 1369-7412, 1467-9868. doi: 10/g9hp74. URL https: //doi.
org 10. 1093/19//0000 1 wwwpsssss//kaa0a qkael08..
[20] X. Shen, N. Meinshausen, and T. Zhang. Reverse Markov Learning: Multi-Step Generative Models for Complex Distributions, 2025. URL Utl.ntt:://arrdvvorr/abbbb/ arXiv:250..37477 [es].
[21] Y. Song, J. Generative Modeling through Stochastic Differential Equations, 2020. URL http://arxiv; org/abs 2011. 13456. arXiv:201 13456 [cs, stat].
D. P. Kingma, A. Kumar, S. Ermon, and B. Poole. Score-Based
[22] Vincent. A Connection Between Score Matching and Denoising Autoencoders. Neural
Computation, 23(7):1661-1674 July 2011. ISSN 0899-7667, 1530-888X. doi: 10/d7h7bn.
URL https: www.i.oor///ooooooooooooooooooo. 10. 11 1111:122/2___00112222..1 C0142.
[23] Y. Zheng, T. He, Y. Qiu, and Wipf. Learning Manifold Dimensions with Conditional Varia- tional Autoencoders. Oct 2022. httpt http:::ss////peeeeeeeeeeee net forum? id=Lvlxg_ ottpo./
H961 tett:tttaaaaiissso1eeeeeeeeeeeea
wor 2dd.ddaaaddds..
12

---

NeurIPS Paper Checklist
1) Claims
Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Question:
Answer: [Yes]
Justification: The abstract summarizes the main contributions of the paper in the same order as they are presented.
2) Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: The limitations are addressed when each result is presented, as well summa- rized in the discussion.
3) Theory assumptions and proof's
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: All theorems, formulas, and lemmas are numbered and cross-referenced. All proofs appear in the Appendix with sketches accompanying the major results in the main text. All assumptions are stated up-front.
4) Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
The experiment details are provided in the Appendix, the code is provided in public GitHub repository specified therein.
Justification:
5) Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described supplemental material?
Answer: [Yes]
The experiment code is provided in public GitHub repository specified in the Appendix and is self-contained. allowing one to replicate the experiments.
Justification:
6) Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: All experimental details are provided in the Appendix and the README, mc file at the root of the accompanying GitHub repository.
7) Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
The theoretical results are accompanied by experiments which show the statistical significance (or standard deviation) across multiple random seeds.
Justification:
8) Experiments compute resources
13

---

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: The computer resources used are described in the Appendix.
9) Code ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics https: /tt:///nuuue cc httpp//ccc/tticcssss.eeeeeeeeeee
Answer: [Yes]
Justification: As a primarily theoretical investigation working with synthetic data, this work does not touch upon any areas that could consist of a deviation from the stated Code of Ethics.
10) Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: Potential positive societal impact of the results, namely in the realm of molecular simulations, are discussed in the Discussion.
11) Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: Not applicable this work.
12) Licenses for existing assets
Are the creators or original owners of assets (e.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 existing assets are credited either in the paper, and in more detail, in the Appendix and the accompanying GitHub repository.
13) New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification:
New assets consist of the experiment code provided in the GitHub repository, and is documented there.
14) Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: The work does involve human subjects.
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?
14

---

Answer: [NA]
Justification: The work does not involve crowdsourcing nor research with human subjects.
16) Declaration LLM usage
Does the paper describe the usage of LLMs if it 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, formatting purposes and does not impact the core methodology. scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [Yes]
Justification: LLM usage scope is described the Appendix.
15

---

A Theoretical Appendix
A.1 Notation summary
Table 4: Notation used throughout the paper.
Symbol
Meaning
Paata
mata-generating density
X E R
Data vector drawn from Pdata
I2P 1R
Encoder
d: R X RD
Stochastic decoder
P
Oracle reconstructed distribution (ORD) conditioned on e(X)
e(x)
(())
Level-set mass at
c(X)
Level-set center-of-mass at X
V(X)
Level-set variance at X
De(v)
Jacobian of the encoder at point y
A.2
Proofs Section 2
First, let's refresh the assumptions used in this section:
Global assumptions 2.3: Pdata C²,, integrable; E C²,, Lipschitz.
Local rank 2.4 rank (y = a.e. on the level set (x(X for the X under consideration.
The local rank assumption will be invoked when needed, while the global assumptions be taken to hold throughout.
A.2.1 Proof of Lemma 2.5
Writing out the first variation condition, have:
8
E
ExxPaato
E y x
|||Y - 1
= 0.
(18)
de
Equivalently, Vn, we have for Y ~ p Peton;x:
(e* (e*+ & (e* + (X),
(19)
by the definition of the ORD. We adopt the same assumptions as we did on as e needs be valid encoder. Furthermore, we will be dropping the in e* we will be assuming be optimal in the arguments to follow to clear up the notation.
Now, is the norm in calculus of variations, assume the following form of the variation for small
e(Y) + en(Y) e(X) + ex(X),
(20)
for Y P Peten.x.ic i.c. ty (e + (n) (E + en) (X)} ((=--en)(X) that is, from the perturbed
level set.
More precisely, to the first order in e, we have:
e(Y) + +vn = + (x()
Which leads the following property of the perturbation 7y:
e(Y) e(X) -=((()) - -(x))))((x)))) O
(21)
16

---

for Y ~ Pe e
Furthermore, due to e) Cº and both being Lipschitz, we have (Y) nx(x 0(1)
Vn, the stationarity condition 18 thus becomes:
Ex Pauta
E Y Prienxx
[IY- =0
(22)
de
As
Using the Dirac delta functions, we can express the ORD as:
Patta (y) 8(e(y) e e(X))
P2 (v)
(23)
Pdata 5(c(2) c(X)) de²
and, equivalently, the "perturbed" ORD as: Using the definition of the level-set mass 5:
Peten.X (v) Pdota (1) + + (e + sX)) Pdata (z) + - (e + en)( de
(24)
/
((X)
- dz,
can define the "perturbed" mass:
Ze(X) = / Pánta()) 8((e+ (e + sn)( dz
Writing out Eq. 22 explicitly:
[f
d 1/11
[ (Pttepx(V)) Poten.((())) dy drf de ==0
= Ex-Paat
where we used Leibnitz' rule move the derivative inside the integral, which is justified since the ORD terms P P (v) vanish exponentially (as they inherit the tail behavior from Paota. while lyy - y'l grows at most polynomially, hence a dominated-convveggncc requirement is satisfied.
We note:
d P Po+em.x(y)
P2.x(0) d log(Piten,x(1))
(25)
de ==0
de ==0
Using the product rule the product of the two ORDs:
d
(Pttenx(())
Pe Peten_x(y)) (v'))
da
=== Potenx(()) P (1) + Pix(y) d
d
Potenx((1))
de
de lelo
and substituting the above, we obtain
d [Poten.(()) P Pete (V)) PP Pix())
d x d
log(P* (og(Pttmnxx())))
de ===
de delezo
de ==
Thus, Eq. 22 is equivalent to:
[f 1/
d
log log(Pitem.x())) x(/)) dy dey'
0 Ex~Paate
log(P* (og(Petenx(())
de c=
==0
(26)
17

---

Now, let's examine the derivatives of the perturbed log-ORDs explicitly using 24 and the linear
variation approximation 20:
log(Pe+en.x (y)) log(Pdata (v)) + log()(e()) =nly) e(X) =n(())))
log (/ Paasaa ((e(() + en(=) - e(X) en(())dd
After taking the derivative evaluated at O, we get for the second term:
(n(w) (x)) 8' (e(v) e(X))
8(e(g) e(X))
and equivalently for the second log-ORD for y'.
We have used d delta undefined, for now.
de - (()(z) + =g(x)) = g(x) 8'(f(z)) where we leave the derivative of the Dirac
For the third term we get:
d de ==D
log Zx(X)= Z(x) / Panto(z) d
-1
de e= 8(e(2) + En(2) e(X) =n())) dz - (m(=)-n(X)) 8'(e(z) e(X)) dz.
/
Paatal))
Z(X)
where we used our definition of the level-set mass 5. The first-order variation condition Eq. 22 thus becomes:
= ExaPa // 1lv - 1'1l3 Poata (1)) Panta (1') 5(e(v) 5(e(v)- (e((v))e(()))))) ((x)))((())) e(())
Z(X)2
(nrlv) - v(x)) (e(eV) e(X))
8(e(y) - e(x)) + (wn(2) - (x)) (e(l)) - e(X))
5(e(x') e(X))
Z(X) / Pantaa() (ı](z) - n(()) o' d'((() -((x)) dz dy dy'
(27)
We want use integration by parts to get rid of the pesky derivatives of the Dirac delta functions.
Now 5(e(2) - e(X)) maps R + IR, and e: TRP + Rt By the higher-dimensional chain rule: = (u
(x()))((x)=) v.ō(u)) (====))((( Dz(2), where D is the total derivative, and D. the Jacobian. We are treating the gradient V.ō(u) row-vector preserve the overall chain rule form.
D (6(e(2)
Above, we have used 8 (e(2) ((ss))e(x))
e(X)) as shorthand for the k-dimensional gradient V..(u) evaluated e(2) e(X), and the product of the delta function with the other terms is dot product from the higher-dimennionna chain rule:
d 8(u - (cl)) + =(1())- (x)))
(a (n(2) nxX)).
==
Thus, have terms of the following form:
1
/
- nx)] Vuō(u)
dz
((X)
u=e(=)-e(X))
Next, the gradients of the 5 function are be interpreted distributionally: for any (smooth, compactly supported) test function 1b(2): IRP + R that vanishes at the boundary, we have:
f ((2) e(X)) dz = - 8(e(z) v.x)))
18

---

or, for the dot product with vector-valued function + that, likewise vanishes:
/ (()) V=δ(...) d= =
/ ((=))((-.) dz
by a "distributional" integration by parts (via the divergence theorem).
Thus, the above chain rule is to be interpreted as the following distributional equality:
I ((e)V e(()) d= = f xb())
De(z) dz.
(28)
Now, let's apply this to Eq. 27, starting with the third term inside the square brackets:
f Paatale (z) (n1() - xx()) d'(e(z) - =(X)) dz
Denote (=() = Pauta (z) (17(2) (X)) E RR
Introduce scalar test function f. i.c., examine the following integral:
/ ((()) [V l f()) dz
Now define: G(z) == (z) Rt Note that G vanishes the boundary due to the Paato factor, hence we don't t need that requirement on f. Thus, using u = =e(z) - e(X):
/ p=)) [V)ō(u)] [()) dz
G(z) [Vo(()) dz
Now examine the integral where we multiply the integrand with D (z) from the left, and (z(=) from the right:
f (D.( G()) ([Váó(u)) De(2)) dz,
assuming the Jacobian (=) has full row rank & a.e. the level set Le(x)- that is the integration domain induced by the function. In other words, this is the point at which we invoke Assumption 2.4. Note that, as stated in the discussion around this assumption, Lemma 2.5 and Theorem 2.6 are simply silent on level sets that do not satisfy this condition."
We will not be simplifying the linear algebra, but instead note that by the distributional equality of the chain rule - Eq. 28, this must equal:
J [II (Vuō(u) De(2)) de / (DI G()) -())))(V - e(X)) dz
Since this holds for any f. we get for / = 1:
/ w(t)) [V)) dz (D. p(=)) [V,ō(e(z) e(X)) de
(29)
We now integrate by parts and obtain:
f ID p(=)) (V=5(e(z) e(x)() dz / V= /V=------- (D. (2) p(z)) 8(e(z) - e(x)) dz (30)
We have discarded the boundary terms since &p + 0 as |||||| + oc due to Pdata being a valid (integrable) density, and 17 being Lipschitz.
Thus, we arrive at the following expression for the third term:
Paota (=)((()--n(X)) (d(=)))))())
V ID. (e) (n(2) ()a))))(X))) ((e())-e(X)) dz
"We believe that adopting stricter requirement the encoder's smoothness and continuity would allow us invoke Sard's theorem and drop this assumption; however, consider that scenario less realistic.
19

---

Let's now apply the same to the other two, "symmetric" terms in Eq. 27: Focus on the first one which depends on y. Ignoring the outer expectation over X, we have:
|| 11u v11 Panta Patta (v') ((x))) 2(X)²
(n(u))n((3)) 8'(e(y) ((X)) 8(e(y) - -(xx)) dy dy'
1
/ Pauta (v - e(X)) dy' / Ilv y/1" Panaa (w)((()) 8' (e(y))
dy.
Z(X)²
Note that the potentially problematic ratio o/4 was merely to clean up the notation.
Denoting now so()) lyy -'111 Panaa (v) (m(() n ()),,, we get by the same integration by parts
(cf. Eq. 30):
-1 =(X)²
/ Pauaa (v') (()))() dul / Vy. (v) lly - w'll Pasa (((l) dy
And, due to symmetry, for the other term:
((x) / Panna (v) ((e(v) dy / VD. (v)) lly 1 1'1l Pana (2')(n(')) ((e(n')-e())) du'
Putting all together, we have the following equivalent of Eq. 27:
1
J dy dy' (((v(e(y)) e(X)) 8(e(y') e(X))
0 Ex ~ Paaca (X
Pauna (5') Vv [D_ (v) lly 1/1 Pasaa (v)(1(()) (X))]
Pauta (v) V [ (1') |1l - 2'1l Paita (1/ (T((()()) mx))
11 1/1 Pdota (v) (d)sa)))))))))) / V. (x) Paona (2) (n(z) - n(X)]] 8(e(2) e(X)) d2
+ Z(X)
(31)
We can now spot that the first two terms inside the curly brackets are identical since Y, iid P hence we can simplify. Thus, Vny
2
/ Pana (w') 8(e(1') - e(X)) dy' / dy b(e(()) - e(X))
0 = ExaPasta
y(X)²
Vy . ID. (v) 11y - 111 Paana (v)(1(()) xxx))
+ Z(X) Ily - 111 Pdoca /aaaaai//V Vv= 1D (z) Paata (z) (17(2) - (x))] 8(e(z) - e(x)) dz)
(32)
The optimality condition will be zero in full generality - An - if the terms inside the brackets are zero, implying VX almost surely:
(v) / Pana (V') 5(e((')) e(x)) dy' /v. DD. (v) 11y I'11 Pamn (V)(n()) v(X))) ō(e(s) e(X)) dy
z(x) / Pana (v') 8(e((")) e(x)) dy' / 11s - Paata()))) 6(e(v) - e(x)) dy
1
V. (z) Paoca (2) (17(2) n(x)) 8(e(2) - e(X)) dz
(33)
20

---

This is the general balance equation that must hold on the level sets of the encoder (as these are induced by the ō-functions in the integral). Defining:
Jo X,n) / V (D. BIv 1/11 Paata (v)(n(v) - ,x(X)] ō(e(u) e(X)) dy
S(X,n) = / /V.- V. [D] P (=)((())))))))))))))))))))))))) - 1f (x))] d(e(z) e(X)) dz,
and WLOG switching around and v', we arrive at the desired expression:
Ey~pi- | (Y: S(X,1) E y. - 8110
Z(X)
A.2.2 Proof of Theorem 2.6
While Lemma 2.5 holds for any 3, we now focus on the case of the squared Euclidean norm: , = 2. As we're continuing from the proof of the Lemma, we adopt the same assumptions and notation as above.
First, let's restate the definition of the level-set center of mass (Eq. 8):
1
/ Paara (v) 8(e(y) - e(X)) dy
c(X) Z(X)
Likewise, for the level-set variance (Eq. 9).
v(x) = f [1y c(x)]l Patcaa (1) ō(e(y) - e(X)) dy
Expanding the norm, get: Ily v/ 2| = 1 c(X) +(X) - =1V c(x))
1/ - c(x)11 22 c(), cxx)
On both sides of Eq. 33, the C term vanishes: First, the RHS:
-2
/ Pasm (v') d(e(s/)) - e(X)) duf / (w ((), - cc()) Pyotaly) 5(e(y) e(X)) dy f J......... !..da dz
7(X)
/ ((x,(x), - c(()) (5') 5(e(1') e(x)) dey) Patta (v) 5(e(v) el (x)) dy f...de ...d =Z(X)
-2
=0
where we used the bilinearity property of inner product (and Fubini's theorem), and:
[w - c(X)) (2/) 8(e(y') - e(X)) dy' = -xx))(())c(xx((xx) =
by the definition of the center-of mass.
On the LHS, get term like:
{ Panaa (g') (()o((e)))-e()
f [(l)))
(v) D.. (x)((())-)(x))) dy
Since the divergence operator is again linear and acts on y-only, we can bring the integral over y/' inside the inner product, getting this term to vanish, again.
The B term gives us on the RHS:
((X) / llv' -((x) Pama (v') 8(e(y') - e(X)) dy' [(eel))e(()) - e(())) (1) dy
/ V (D. (2) (1((2) - n(x)) 8(e(x) - e(X)) dz = Z(X) v(x)z(x) f...
-v(x) /V.
V. V.-ID. (2) Pauta (2) (n(=) n(X)) 8(e(z) e(X)) dz
21

---

On the LHS:
/ Pasaa (v/)
du' / b(e(v) ( (D] (v) 11 - cxx)]| Pasaa (v)(r(v)) n(())) dy
= =v(x) v. 1D (v) (v)((())) - m(x))] o(e(v) - e(x)) dy
Thus, we have for the optimality condition 33 in terms of c()) and V(X) on the RHS:
1 / Paam (w') 8(e(1/) e(X)) dy/ ʃ 8(e()) c(x)))ly 1/1 " Panta(l) dy
Z(X)
f v (D (2) (1(2) (xx)) 8(e(z) - e(x)) dz
z(X) / Pana (g') 8((((')) - e(x) dy / 11y c(X) Pa Pautaly) (v) 8(e(y) e(X)) dy / dz
1
+ +V(X) VV-1ID? (2) (n(=) n d(e(z) - -(x)) dz
= -2 / V= - (x) Paota (n(=) (x))) B(e(z) - e(X)) dz
On the LHS:
/ Paata (v') -(x)))-((()))))) dol I V. 1D. (v) 1evllly/ y'1l2 Paaa (V)(T(v)) r(x)) d(e(y) - e(x)) dy
= -f Panm (((())))))) - e(x)) del [ V 11y- - c(x)11 Pasn (y) D.) - x(x))) dy
+(X)) / vy DD Pana (v)(n(v) - (xx)) 6(c(y) e(X)) dy
= -x) / Vy, llly c(x)1 Paun D.) (n(() - (xx))] 6(e(v) - e(X)) dy
+ +V(X) / V, PPana (.) xx)) 6(e(v) e(X)) dy
Thus:
z(x) / V. -1ly---- 11ly-cl -11y-((x)1 Pasaa (v) (v)((())) - ((x)) 8(e(v) - e(x)) dy
+V(X) / -
[Pama (v) (T)(1())))))))))))))))()))))(() (xx) ō(e(v) - e(X)) dy
-2((x) / V= Pa (z) DJd (=)()()) - ((X)) 8(e(=) e(X)) dz
Since y and z are just integration variables, we can rename them and we bring over the second term on the LHS to the RHS. Thus we get, Vn, VX almost surely:
/ vv, llllyyd 11y-o(X8 2(x)11 Pat (w)()) m(X))) 8(e(y) - e(X)) dy v(X) z(x) / (V,- (y) DJ. (n)(() - (x)))) 5(e(y) e(X)) dy
(34)
In the above, we are dealing with divergences of the following form: Vy [f(y where fi(w) = |yy c(X))1 Panta (v) is scalar (with fz(y) V(X) Pasta (v) on the other side), M(y) = z(X)
22

---

De (v) a x & matrix, and v(v) = 7()) -x( vevtor We will only expand the first term, namely:
V. [f()) v - f()) ( [M(())(())]
Eq. 34 thus decomposes into:
/ Vyfi(y) [M())(()) B(e(y) e(X)) dy / fi(VVV, [M()) 8(e(y) e(X)) dy
/ V fa(()) (M())(()) 8(e(v) e(X)) dy / fa(u)V (M()))(V) 8(e(y) e(X)) dy
(35)
We wish to go from integral equality to statement about the integrands. The argument will be made as follows:
1) We will argue that the second, divergence terms become negligible in the first-order station- arity conditions (Eq. 22).
2) At this order (i.e., at order E), we will show that the integrands of the first terms must coincide almost surely.
From matching the second terms (with the perturbation inside the divergence), one might expect that on the level sets, would have
f1(1) a.c. f2(1),
which would lead to spherical level sets:
V(X) 1l (1--(()) xX]]] Z(X)*
which cannot be justified.
Another option is for the divergence V, [y) to vanish. The "trivial" solution My))u(y) == D.E (v(()) - (x))
0 would imply that all perturbations are tangential the level set; this would fly in the face of the variational argument where we allow v) to vary freely under mild assumption (i.e., that it is differentiable and smooth).
Thus, we aim to show that the divergence V, (y))())) vanishes when integrated when considering the first-order optimality conditions. As refresher, have denoted the unperturbed level set (manifold) as:
Le(x) ty R" e(v) = e(X)}
We have already assumed that D. has full row rank h a.e. on Le(x)- Thus, dim (La(x)) k. The normal space at E Le(x) is spanned by the rows of De(y) (or, equivalently, the columns of D.T (v)). Thus, any small displacement öy E R of can be uniquely decomposed into the normal and tangential component:
8y dy Syl, with De(y) dy 0, De(y) dyg # 0
Let's consider now the perturbed level set for a small E:
(()te)(x)) = {u (e + (n) = (e + {X)
As we've already shown: to the first order in e, y is on L [(e+)n)( then we have Eq. 21:
e(y) e(X) = -E(n(y) - (X)]
Set y yo ōy for some reference point yo Le(x)) and expand around yo:
e(y) e (yn) R D. (yu) (y yu)
Again, for y the (close-by) perturbed manifold, that difference must equal -5 (n(v) 17 (30)). So
23

---

De (30) ly yo] ~ - (m(2) 1 (yo))
We have assumed that is Lipschitz, meaning that () // (vn) is O (lly yo 1)-
We now wish to show that the normal componens the displacement öy = y yo is forced to be (E)).
Denote the projection onto the row space of De (yo) as IT (i.e., the projection into the normal space of the level set at yo). We have:
Da (vo) by De (vo) (ōyl + by.) ) = Dr (vo) IT, by
Thus:
[e öy = 11, 1-8 (n(w) - 7 (vo)) 0(=)) O(s)
Since (y0) has full row rank by assumption, its pseudo-inverse 3, hence state:
I By = IT [De(wo) (-E(n(v) - v(v6)))] + O(e²)
Taking norm both sides:
|لا = 11l1 [Do)t (-e(n(v) v(uo)))) lel De (1(r r73)))
$ |<| (bo)* (1(((y) n(6)))l,
where repeatedly applied the Cauchy- Gauchy-Schwarz inequality and used the fact that the norm of projection operator is one.
We will use the fact that both and v are smooth and Lipschitz with constants L, and La Since De is C and has full row rank k at yo, its rank remains & in small open neighborhood of yo rank can only change in discrete steps, and continuity prevents drop for small enough neigborhood. Consequently, the norm of the pseudo-inverse is also (locally) bounded by finite constant C. Hence, we have
by || < |e| C T,,l|6v|| + Os²)) -
i.., any normal displacement from the original level set is O(s).
Moving on to the divergence: we note that the divergence operator is linear operator and can be decomposed into the tangential and normal component (relative the current level set):
(V 9.11)) + (Vy.1)
We have, Voy: V ky. = 0 and Vy tyg = 0
Since the terms we are taking the divergence over are of the form:
Vy [D] (v)(n(v) (x))),
we thus have
Vy - = V..1 [D. (v) (n(v) - n(x))]
(since the column space of DTT spans exactly the normal space). Thus, only the normal components of the divergence will play part.
Denote the vector field that the divergence is taken over as
u = [M(v)((w)]= = [D_ (v)l(()) n xX]]]
Namely, since | || 0ss) (as shown above), we have again by Lipschitz-ness of the perturbation:
(1(7(1) - v(())) = O(ly voll) O11 II) = O(e). So [ O(e), since the norm of the Jacobian is bounded, as well, due to its regularity.

---

Next, let's define a "cylindrical" region R₂ around the old level set L((x):
R = (y d (v, (x(x)) S ce)
for some small constant e that the perturbed level set contained within R. Finally, let's apply the divergence theorem:
/ R.
V. u(y)dy S (w)-n(y)dS - 11u//oo Area( (
[ (s) for y in the region, and the measure ("area") of RR, is most (e)) times the measure of the level set, as the "outer faces" of the "cylindrical surface" are (ss²). Hence the above integral is at most ((=²)x x the p- k-dimensional measure of Le(x). After accounting for the fi and f2 factors in Eq. 35, we note:
1) JR. f12(1VV u(y)dy I SuPgER, f12 JR. vy. dy. As we'ne we're assuming V(X) (Eq. 9) finite, this means that f1,2 is bounded. Hence the inclusion of the /1,2 term doesn't change the (s=2) scaling of the flux.
2) If Le(x) were extend to infinity (the "times the measure" part), the Pdata factor in fi. would kill this contribution, as Pdota vanishes quickly enough for the level set variance
((X) to be finite (by assumption).
Thus the (flux) integral is at most O (s*).
This means, that when considering the first-order optimality condition 22, for which we have obtained and expression of the form:
0 0= lim 1 (F(()) F(O)) - -1 Vyf(v) ... dy like terms) + / f(y) .... =+0
...ddd.. dy like terms),
where we combine the terms from both sides Eq. 35. We have shown that the second group terms is (s))), thus it vanishes in the limit and cannot play role in the first-order optimality conditions.
Now the second point. That is, assume that the second terms, where the variation 7) appears inside the divergence, vanish in the first-order stationarity condition.
Thus have the following integral equality, VX, Vny.
/ V. [1ly c(X)) Panta (v) D. (v) ō(e(v) - e(X)) (n(v) - mx()) dy / (X)
V Paata (1)) D. (v) d(e(y) e(X)) (1(()) mxx)) dy
((x)
Now, denote:
P ()) = Vy [1l e(X))) Panta (w) D. (1) d(e(v) - e(X))
and
X) V(X)
V, PPanaa (v) D. 6(e(v) e(X))
Z(X)
So, the identity becomes:
f [i((.x) (n(v)) - n(X)) dy = / X2(1, (n()) - -(x)) dy
While we have put the 5 functions inside the integrands to be compared as to leave the (what are to be) test functions 7) clearly separated, one needs keep mind that they will induce the integrals to be over the level set manifold surface measure. Additionally, Fi and F2 now distributions. Also
25

---

note that we have assumed no specific constraints on 7) besides them being smooth and Lipschitz they are in the same function class as (cf. Assumption 2.3).
We will proceed via a proof by contradiction. Assume AF S {(y, X)) with nonzero measure w.r.t. dy data (2)) such that:
Fi ())x(((((()))) E A} + F2 X) (V, E A}
Furthermore, assume that Ty (A) C {y c(v) = c(X)) where Ry the projection to y- In other words, assume that the integrands differ on "subsection" of the (unperturbed) level set (the delta functions inside F would make the opposite - A & Le(x)) impossible, anyway).
We need show that 17 s.t. the two integrals differ on a subset S R' with Pdota (S) as the integral identity is a.s. in X. To that end: S x (A) = {X Byst s.t. (V, E A} and Paata (5)
A has non-zero mass under the joint measure.
Let's pick AX E S an 7/x to be a smooth function such that: [0,
outside a small open neighborhood of Ux inside Ux
(x)) = nonzero and positive in all components
where Ux A Lx(x) is open set. That is, we're using the hump function and partition of unity approach common in calculus of variations [13; 9].
In this approach, we wish to build a global that belongs to the function class of the encoder. Pick The compactness guarantees a finite sub-cover a compact subset Sx S S s.t. Paata perturbation 17 by "gluing" together the bumps 7/x for (Ux, ) Then, we can construct global that sub-cover using a smooth partition of unity p,(((()). As we can select the cover to exclude the non-full-rank points the encoder's s Jacobian, this guarantees that the partition of unity (y(e(y) is smooth and exists. Then, we define our "global" perturbation as:
m
(g))
Σ ps (e(()))n, (v)
j=1
As 1) is then a finite sum of smooth bumps multiplied by the smooth partition-of-unity functions, both C² and Lipschitz.
Furthermore, have for such n at least of the nx, terms active in the integral AX E Sk. leading to:
f (R(, (nll)) - n(X) dy #
- dy
Since the integral equality must hold for any 17. 17-s form a rich function class, and our picked satisfies the requirements the class, we arrive contradiction on non-measure-zens set as Paota (Sx) > 0, thus must hold (y, X) X),
Thus we have, in the first order of e:
V(X) Z(X)
Vy [llv 2 (v)] D. (1)
V, [Pasta (v)) D. (v), (v),
almost surely in X on the level set of e. Taking the gradients, we obtain:
as in V(X)
2 ((axa) + -(x)]| Pana (v) D.T (v)
V, Panaa (v) D.. (v),
Z(X)
and finally:
2(y - c(X))
(T) in V y
Pa Paaca (v)
D.T (v).
v(X) x |I c(()]l2
Patta (3)
For the existence of the bump functions and the smooth partition of unity, e.g. Lee [13], $2.
26

---

A.2.3 Proof of Corollary 2.7
We adopt all the assumptions of Theorem 2.6. By the fact that we are at extremum y", we have Vy Patta (y²) = This gives us the following relation for Eq. 7:
2(y" - c(())
(V") 0.
((X)
11y" - c(x)]1
2(X)
V(X) finite (by assumption) and minimized on the level set due to the encoder's optimization objective 2, is Z(X), with additionally ((X) > O, since smooth level that contains extremum must have some mass. Thus, the only way the denominator could go to (-) infinity
Ily" dx)il was itself approaching infinity, which is not considered by assumption: this scenario
represents "trivial" minima at distance approaching infinity.
As the denominator on the LHS cannot go to infinity, must have:
cx()) D (y))
This satisfied either by:
1) (y' c(X)) = O, that is, the encoder aligns the level set so its center of mass coincides with the local extremum (most likely maximum),
2) -y"- c()) D (g²) = that is, the normal projection of (v" c(X) is exactly zero, meaning that (y* c(X)) lies in the tangent space of the level set at the extremal point.
A.3 Proofs Section 3
A.3.1
Comment: DPA recovers the data distribution
To refresh: for optimal encoder/decoder pair, we have by definition:
d" (X), x
That is:
d² d (x)e (X)= = =); z); Vz
Treating Z (X) as a random variable, we have:
d'(Z= = 2, (xl =)
Thus, for any measurable set A:
P(X A)= = P(XE 2 P2(2) dz
and
IP(d (2) P(d - Z 2) Pz(z) dz
d (xl = 2). the two integrals match, and by the law total probability:
Since (Z= =
d' (c(()))))
(36)
A.3.2 Proof of Propositions 3.5 and 3.2
We assume that the manifold M can be exactly parameterized by smooth, injective in the first K components e1:k MM- + R and that our encoder function class expressive enough achieve this.
27

---

If we assume this 61:K is optimal among K --dimensional encoders, this must necessarily imply by the definition of the ORD:
d' (lexx (x); ((++))p)) P Ki
where d' is the accompanying optimal decoder.
Now, due the injectivity of the parameterization, we must have for the RHS:
P Ei.k (x) Law(X)e(X) - e(x)) 8(x), Va.
(37)
In other words, we get an almost-sur sure equality:
d' (e)xx((X), a X
This implies for the two terms in Lx [x(en: d'] (Eq. 11):
Ex By ~ (x.e1:
[||X --1 = Ex [[XX - Y|1
and
Ex E Y.Y P.
|1 - r'1 = = -1-1111 Ex x
Since Lile, df 2 O, this is indeed the global minimum. Thus, the encoder is an K -bestt-aprrxiimatig encoder, and the manifold is K' srraparmetteiiaalleee by (fI.x, with K' K.
A.3.3 Proof of Theorem 3.4
We consider the setting where the data X are still supported on a K -ddimensional manifold, which, however, might not be globally parameterizable by single, smooth encoder e with K output dimensions.
We have defined the K" o'-heet-apppoxxmmaing encoder as the encoder that minimizes the loss the when considering terms up the K'-th term only:
(e", d") arrmin Σ [x d]
again taking the weights to be uniform, and which achieves the globally best energy score in its K'-th term:
L [(e", dd)] min Lufe, d].
e.d.k
This also automatically makes the manifold K' --parameterizable, with any remaining manifold variance being non-explainable with DPA.
By the fact that for E (0, 2), the energy score is strictly proper, this encoder/decoder pair is indeed the unique global optimum for -i-dimensional encoders, with Pa+ Pi (x) being the unique distribution minimizing [(e,(("",")))
Now consider the overall p-dimensional problem given by Eq. 10. The terms for K' + 1... cannot do better than L je", d"] by the definition of our K" d'-best-approxiiaain encoder, thus the best an encoder can do is to output the same distribution Pa-.ej.e.
Next we observe that the K' dimeinnnllo optimization is nested subproblem of the p-dimensional one. Thus, when optimizing over p dimensions, the optimal encoder must coincide with the above K" K''bbest-approxee aperxxaproxi encoder in terms up to K", otherwise it would contradict the latter being the global optimum among K' K'-dimensional encoders: one would obtain lower loss by simply using first K' components of the p-dimensional encoder.
Hence, the K'-best-any approxippriiiinnnggg encoder is the global optimum for p dimensions and we obtain for the optimal solution:
Pa- (X) (X))) for k K"
*This has been given example on page 6 of Shen and Meinshausen [18] for general invertible e".
28

---

Next, we wish to show that the "extra" dimensions are independent of the data conditioned on the relevant components K'. Consider the (K' + 1)-th component. Note that (X) and +i,k+1(X) when viewed as joint distributions (over dimensions of e"), form a filtration; denote FK- =
S FK++1 = (X(e)(X)))....... cx (X))).
By Eq. 13, we have
X| X|FK',
hence
B[((X) = Fx)
for any (Borel-measurable) function /.
Let Z be FK - but not Fg  measurable. The claim from the theorem then equivalent to the following conditional independence:
Z LL X|FK
or
flg(z) = [((
for any pair of integrable functions f.g.
We can prove this claim in the following way:
Elg(2)
E[o(Z)
Fic4
Zig(Z) ((
Blo(Z)|F<d) ([)
[((())F F(
[F FR! meas.
The other dimensions K' 2.... follow by the fact that the filtrations are nested, meaning that the above derivation holds when F is replaced by Fx etc. Thus, we have shown that the "extra" dimensions are independent of the data, given the relevant first K' components.
A.3.4 Discussion on Remark 3.6
We are considering the "exactly-parrae "aanetttypprrmmeerizaale"""" manifold scenario. We have shown in Proof A.3.2 that the K-th terms L [e, d. are identically zero and thus the global minimum.
As in Proof A.3.3, we observe for the higher terms (e.g., K 1-th):
Exhev~ (x) (8) [||
Yl
>
and
Ex E
[||'' -Y113] >
Y (1
unless
Pareixx (x) = 8(X) Pa (X)
In other words, the only zero variance distribution is the delta distribution, which is (by assumption) the distribution induced by Par.ei,x- Thus, the encoder that parameterizes the manifold in the ei:K dimensions and outputs the same distribution for terms K + p is the optimal encoder for the terms K p. We will keep e* to denote this encoder in the following argument.
Next, we argue that typically for p > K, this also the optimal encoder for the first K dimensions, thus the K f--bes-approximmatng one, or equivalently, the global optimum.
Again, assume that all the weights Wa [0, 1] are uniform, i.e. p+I. As discussed in Shen and Meinshausen [18], it remains open question whether optimal encoder is necessarily the one that minimizes all the terms in the loss simultaneously (which the case for the terms K p when the encoder is the K- happen for parameterizable manifolds as >> K.
one), the following argument will examine what is likely to
Suppose that there is another, different (in the sense that it reconstructs different ORDs) globally optimal (ē, d) pair that, by "sacrificing" perfect reconstruction at dimensions K improves
29

---

for terms k = I.....K K to such an extent as to strietly beat the manifold parameteriiin encoder e* in the aggregate loss. Note that for the latter, the "partial" manifold parameterization should typically (i.e., for non-pathological manifolds) already be reasonably good encoding for dimensions < K (in the energy score sense), so T& (e", d²) is unlikely to be very far from the minimum for each of these k.
We will be using Proposition (together with Theorem 1) of [18], which states that the optimum. the two terms in the loss are equal and thus focusing on the reconstruction one; as we are conjecturing the existence of another global optimum, this is valid.
Since the data manifold is K dimensional, the encodings for k < K cannot describe the manifold, leading to imperfect reconstruction: ||X Y1 > O. This is true for any encoder/decoder pair, and we can denote the global minimal loss when considering l dimensions (i.e.. single term in the optimization objective) as:
Re,min ==
min Lx d] eask.d
Furthermore, typically for k < K each encoding has to "describe" K - k additional "directions" of the manifold, meaning that P P (X e1a(X) 2) becomes more spread out on the manifold as k decreases; hence, typically:"
Ri.min Ry min for ls < K,
which implies
Lule, d] 2 L di for k k'
for any reasonable (assumed-to-be) optimal ē as the manifold's dimension is K k' k.
Suppose now that for some S C {l...., K 1} this pair reduces the loss Lu(ē, d) below La (e*,d") which implies (by "breaking" the manifold parameterization):
L [2, d] [x dd] = 0, vk: 2
Denote the hypothesized improvements (in absolute value) on the k K terms as AS. and the latter costs for k 2 K as
Ak Lale",d dd] - Lu(ē, d] k<K, K, (> Of there is improvement)
Eb == Lu(e,d d) - dule",d") k2 K
0).
Thus we have the following best-case trade-off that needs to hold in order to (ē, d) to be optimal: (K 1) max A > (p K) mines
(38)
Eg are nonzero 5 is the only zero variance distribution:
E = Lx(ē, La le", d'|20>0,
and the improvements Ag are bounded above by:
A < Lx(e", ,d' Run
k << K.
20
As discussed, typically, La(e", d=] are unlikely to catastrophically large as the partial manifold parameterization should do reasonably well on the energy score loss. Denoting:
M max La [e", d'), kcK
min Rk.min
k<R
we can express Eq. 38 as:
(K 1)(M m) (p K)
As the LHS is constant for non- pathpoloogiaaa manifolds and the RHS is linear function of p. Eq. 38 cannot hold as increases. Stated more formally:
'While there might be other pathological scenarios. in most cases this would imply the true manifold dimension to be K k', which contradicts the assumption.
30

---

Lemma A.1 (The parameterizing encoder is typically globally optimal) Fix M, m,c 0 as noted above. The manifold manifold-parameettt parameterizing encoder e* is the unique global minimizer whenever:
K K-13
p
(M m) + K.
(39)
Thus, typically, an encoder that parameterizes the manifold in the first K dimensions is an optimal encoder. Then, due to Prop. 3.5, the results in Theorem 3.4 follow.
31

---

B Experimental Appendix
B.1 Data distributions for the experiments
(a) The data for the standard Nor- (b) The data for the Gaussian mix- mal. Figure 3: The data for the Gaussian examples.
ture.
Figure 4: The data for the Müller-Brown potential example. Note the dearth of samples in between the potential minima.
B.2 Score alignment
Our method of extracting level-sets is relatively rough and uses matplot contour function. Thus, there are significant inaccuracies in estimating the level set statistics (e.g., c(X)), prompting us to discard scalar factors on both sides of Eq. 7, and evaluate alignment instead. Such inaccuracies are especially impactful around the critical shell, where slight numeric error can erroneously flip the direction of the vector on the LHS of Eq. 7, as illustrated in Fig. 5. As evident in the figure, the directional alignment remains very good. To account for these drawbacks, opted for absolute cosine similarity as a measure the match between the level-set geometry and the score as predicted by Theorem 2.6.
LHS (radial term DT)
RHS (score D_e
Figure 5: Signed score alignment: sign flips due to the inaccuracy of estimating the level-set statistics.
32

---

B.3 Minimum Free Energy Path further results and discussion
In the use of unsupervised learning methods in computational chemistry, one typically wishes to find a good collective variable [4]: representation that can approximate the important transitions between physical states. For the Buller-BBownn potential, the states are represented as minima of the potential energy. The ideal collective variable is, naturally, the minimum free-energy path (MFEP), as traveling along the latter is the least-costly way for the system to undergo transitions. Thus, the better method is able to parameterize the latter in a single component, the better this encoding can serve in downstream tasks, such as enhanced sampling [4], where the encoding is used modify the potential to "guide" the system towards undergoing transitions between states.
As discussed in the main sections, the fact that the DPA aligns the level sets with the force field naturally induces an approximate parameterization of the MFEP, as changing the encoding value represents moving to the closest level set, and this move aligned with the score the force field. In this section, we present further results of the automated experiment used to generate Table 2.
In the experiment, first train the encoders the data (ef. Fig. 4) across 25 different random seeds. Next, we encode the MFEP, obtained by the string method [8]. We find the best encoding component by regressing the latter against the cumulative component with the largest Re coefficient. If the encoding is the exact parameterizing parameter of the MFEP, the latter should be almost 1. Next, we select the range of the encodings travel over to obtain our path approximation by finding the closest encoding points to the start and the end of the MFEP by using the same regression. With this, we travel over the range of the encodings and, at each step, predict the next position z by finding the point corresponding to the next encoding in the range via a noot-finding algorithm; that is. by inverting the encoder: Inext e (Znezt). Additionally, we observe that for the DPA. the best encoding is always a monotonically increasing (approximate) parameterization of the MFEP, and the level sets are connected. Hence, we also present an alternate method, where we generate the MFEP approximation by simply taking step in the direction of the gradient of the encoder, obtained by automatic differentiation. The results essentially coincide with the more involved root-finding results, and are used in Fig, 6 together with the results obtained from root-finding for the other model architectures.
length of the MFEP and picking the
Figure Best MFEP parameterizations for teft-to-right: DPA, Autoencoder, VAE, 3-VAE, 8- TCVAE for two arbitrarily selected random seeds: top: 43, bottom: 63.
The results consistently demonstrate that the DPA performs much better scalar parameterizer of the MFEP; much tighter approximations can be obtained (without retraining the model) by decreasing the step size (in the latent 2) and manually adjusting the range of the encoding, as illustrated in Fig. That is, the method of obtaining the approximate paths is relatively crude, but unbiased with regard the model architectures. As noted in the main work, the Autoencoder typically obtains the correct direction between the first and last minimum, while the VAE does not and performs worse, with the level sets often being disconnected, as can be observed in the bottom row of the figure.
Below, we summarize the metrics used in Table 2, where we present the results obtained from 25 runs with different random seeds. We retrain the models on each run and find the best parameterization of the MFEP possible by the encoder.
Average closest distance from the parameterizing path to the MFEP dputt dpaih--MFEP. -MEEP.
Average closest distance from the MFEP the parameterizing path dMFER-ppath
33

---

14
0.
0.8
62
1.0
0.5
0.0
0.8
0.
90
02
0.4
Figure 7: Tighter approximation of the MFEP for the DPA.
The Chamfer distance de: the average of the above two distances.
The Hausdorff distance du: the maximum distance between either paths.
95th -prr-enteetlee error: similar the Hausdorff distance, we present the larger of the two 95th -percentiles for the two distance sets.
Additionally, we examine which encoder component was the closest parameterization - denoted as I|| and present the average over the random seeds. Of note here is the prineipal-componets aspect of the DPA, as the first component always represents the closest parameterization, as the MFEP is the single most important feature for describing the Müller-Brown potential. This does not hold for other models.
B.4 Independence of extraneous latents
Table 5: Datasets used for the extraneous-latene independence experiments (Sec. 4.2). is ambient dimension; K intrinsic dimension. "Generating map source" summarizes the construction.
Dataset
P
K
Generating map source (brief)
Gaussian line
2
Line embedded RR² with orthogonal Gaussian noise (I-D signal small normal perturbation).
Parabola
2
1
Curve (t, (²) sampled over bounded interval.
Exponential
2
1
Curve (t, exp(t)) sampled over bounded interval.
Helix slice
3
2
Helical sheet (cos t, sin f, t) (thin slice of a 3D helix).
Grid sum
3
2
(a, 2) with constraint = x ae-ty (two free coordinates; one dependent).
S-curve
3
2
Standard sklearn S-curve surface RR viv-ddimesioonl manifold).
S-curve (,3
2)
3
2
Same S-curve as above; DPA trained with B = 2 to show empirical
robustness.
Swiss-roll (3D)
3
2
Standard sklearn Swiss-roll surface in 3 (two-dimensionall manifold).
Details about determinism checks:
For the regression part, we fit the following nonparametric regressors: cubic B-splines with up to 64 knots, random forests with 400-1000 trees, high- degree polynomials. For the conditional entropy H(U Z) H(U, Z) H(2) we use k-NN Kozachenko-Leonenle estimation with 5.
We adopt double Conditional Randomization Test (CRT). We (i) fit Gaussian process to model the conditional distribution of U given Z obtaining the mean μ Z) and variance dZ(:): (ii) compute the observed test statistic the HSIC (Hilbert-Schmidt Independence Criterion) between U and X; (iii) generate B bootstrap null samples by drawing U, ~ (N( o² d2()) and decoding the latent representation [Z, U.) through a trained decoder to obtain synthetic data Xx; (iv) compute HSIC statistics for each bootstrap sample; and (v) calculate the p-value as the proportion of bootstrap statisties that exceed the observed statistic. The test is repeated N times with fresh decoder samples to assess the validity of the test under the null hypothesis that U 1 X Z.
Details about CIT:
34

---

Additional diagnostics on the S-CURVE
We train DPA with port unconditional and conditional (given Z1:2) distance correlation and mutual information be- tween each latent Z; and the data X, plus conditional linear R summary. For distance correlation we use the dcor package; its partial estimator is used as proxy for conditional independence (near-zero values are consistent with 74 IL X 71:2- For mutual informa- tion we use sk1 earn. fww.trreeeeeeeeeeeeeeeeeeeeeeeee ion. mutua 1 iino_regreese on; we report the "per-coordinate maximum" MImax = max, / (2; Xj), and its conditional analogue after residual- izing both Zi and X; Z12 with kernel-ride regression. This yields a conservative check for residual dependence. All metrics are computed on held-out test set after standardizing and Z.
on the S-CURVE dataset (K intrinsic dimension). We re-
Table 6: S-curve diagnostics (N = 10,000, A 2). "cond." conditions on 1:2
Distance corr.
Mutual info. (nats)
R² (linear, cond.)
Latent
uncond.
cond.
uncond.
cond.
value
Z1 (relevant) 0.71
1.49
Zz (relevant)
0.57
0.67
Za (extra) 0.27
0.07
0.29
0.13
0.0006
As sanity check on the intrinsic dimension, conditioning on only one latent (K 1 hypothesis) leaves substantial linear dependence: [²[Xx; Zz Z1) = 0.94; under the correct K = hypothesis, Rx(;; Zz Z122 0.0006. Together with the main tests in Sec. 4.2, these complementary diagnostics support the uninformativeness of the extraneous latents result.
B.5
Experimental Details
B.5.1
Experiment code
The code required to reproduce the results is provided at github. com/ aadd eblleb/Diii Disstribuiion Aut ont README md document therein for full details on the code structure and details.
Please refer to
The hyperparameter choices, the optimizer types, used for generating the results are provided in Sec. B.5.2 below.
B.5.2 Experimental details
Figure 1 (Gaussian Score): Uses 2D standard Gaussian data with analytical function.
Random seed: 42. DPA models trained with 3 = 2, deterministic encoder, stochastic decoder, latent dimensions (A= 3), 4-layer networks with 256 hidden units per layer, residual blocks enabled. Training used standardized inputs.
Table 1 (Score Alignment): Evaluates cosine alignment between V Vxia and Vx log p(x) on a 100 x 100 grid over [-4, 4/2. Tests two distributions: (i) 2D standard Gaussian and (ii) trimodal Gaussian mixture with means μn (-1.1, -111,, #2 (1.1, 00.9.. μ3 (-0.33, 1.0), equal weights wj = 1/3, and shared variance o² = 0.66². Density cutoff:: 1% of maximum density.
Reports mean cos f), standard deviation, and 95th percentile across grid points with sufficient density.
Figure 2 Müiller-Browr Potential): Uses Müller Broler potential energy surface with standard parameters. Collective variables learned via DPA with 3 = 2, 2D latent space (k: = 2), 4-layer net- works with 100 hidden units. Training trajectories generated via molecular dynamics at temperature T = 300K. Random seed: 42.
Table Figures 6-7 (MFEP Comparisons): Compares five autoencoder variants on minimum
free-energy path (MFEP) parameterization for the Müller-Brown potential: (i) DPA (3 = 2, deter- ministic encoder, stochastic decoder), (ii) standard autoencoder (AE), (iii) vanilla VAE, (iv) 3-VAE
(3 4, Higgins loss), and (v) 8-TC-VAL VAB (a 1, B 6, = 1). Each method trained across
35

---

N = 25 random seeds (42-66). Training data: 10,000 Müller-Brown MD samples at 7 300K. All models use 2D latent space, encoder/decoder with two hidden layers 100 units, ELU activations. Training: 1200 epochs, Adam optimizer with Ir 10-3 3 (VAR variants) or 5 x 10-2 (DPA), batch size 5000 (DPA/AE/VAE) or 256 (TC-VAE to ensure batch 2 for log-density-ratio estimation). Reference MFEP computed via string method with 32 nodes, tolerance 10-7. For each model, summary statistics computed by dropping the worst-performing seed (by Chamfer distance) and reporting mean + standard deviation over remaining 24 seeds.
Table (Conditional Independence Deterministic Tests): DPA models: 3 E (1, 23. informative extrancous latent, 4-layer encoder/decoder with 100 hidden units, batch size 128, trained for 1000 epochs. Diagnostics: (i) R² via cubic B-splines (up to 64 knots), random forests (400-1000 trees), high-degree polynomials; (ii) intrinsic dimension via Levina Levina-Bickel MLE with 200 bootstrap samples; (iii) conditional entropy HU Z) via k-NN Sozachenko-Leonenie estimator (k: 5). Random seed: 42.
Section 4.2 (Conditional Randomization Test): Double CRT with = 500 bootstrap samples, Neep 50 repetitions. Gaussian process fits conditional U2 Z1 using RBF kernel. Test statistic: HSIC with Gaussian kernels (median heuristic for bandwidth). Null samples drawn as Up ~ (N)(Z)), o² (Z))) and decoded via trained decoder. Significance level 0.05.
B.5.3 Existing assets
The experimental code makes use of modified routines from the micolvar project [4]: https: //gith siihhub... com/ haamm//lgloooaatteoooo /ml molcolaas available under the MIT license. Additionally, the data for the Müller-Brown potential examples is bundled with the examples provided in the mlcolvar repository, and reproduced in this supplement for convenience.
The original license texts are kept verbatim thi higa aeseenseespoojecttt/CCCSSS
All other dependencies are installed via pip (cf. the requirements txt file provided); third parrt_lliennsttteeaaaaaaaeeeeeee wwceenseeeeeiR THIRD PARTY_ LICENSES md.
their names, versions and licences appear in
B.5.4
Compute resources
The results presented the main paper were run on a laptop and required about 10 minutes for each model training and about 20 seconds for the detailed plots.
The MFEP parameterization experiment was run on single Nvidia V100 GPU, taking hour and 38 seconds of wall time. The Independence experiments took roughly hours on single Nvidia V100 GPU.
B.6 LLM usage
We used GPT-4 interactively as a search tool and sanity checker for content created by the authors. All LLM-generated suggestions were independently verified against the original literature.
36


Reviews:
[{"review_id": "kAUO7wlKIR", "reviewer": "Reviewer_Dpaa", "summary": "The paper analyzes the model called Distributional Principal Autoencoder (DPA), a recently proposed method which performs distribution matching using the energy score. The main contributions are showing that optimal encoders exhibit certain structure depending on the data distribution (Lemma 2.5), and the extra dimensions of an optimal encoder for DPA beyond the dimensionality of the manifold do not carry extra information (they are noise, Theorem 3.4) for \\beta < 2 in the DPA. In particular, it is shown that for \\beta = 2, the optimal encoder for DPA has level sets that align with the Stein score of the data distribution (Theorem 2.6), and numerical examples are given for this \\beta = 2 regime.", "questions": "1. How do we interpret Equation (6) for values of \\beta other than 2?\n2. Are there properties similar to those established in Section 3 when \\beta = 2?\n3. The connections to dimensionality estimation are interesting, do you have any empirical validation for it?", "limitations": "The main limitation is that the theoretical results rely on finding the optimal encoder, which is not well-defined for \\beta = 2, and is likely not computationally tractable for most problems. This is a common limitation in machine learning theory papers. The numerical examples support the overall hypothesis of the paper, but there are no numerical examples for the manifold results.", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths:\n1) The theoretical analysis is novel and seems correct (though I did not have time to check the details of all proofs).\n2) The result that DPAs can reveal the intrinsic dimension is fascinating. \n3) Connection to the Stein score function can be useful and connects this paper to other research on distribution estimation.\n\nWeaknesses:\n1) The is heavy mathematical notation that I found difficult to follow at first, especially since I was not familiar with DPA. This is a lesser known model in this community at the moment, and readers may not be familiar with the topic and appreciate the analysis. This is the main reason for the lower score.\n2) The numerical examples provided were low dimensional and simplistic. and the connection with score alignment was hard to see. Only comparisons to simple autoencoders and VAEs were given. Details on the setup for the experiments were not included in the paper (though code was included in supplementary material).", "was_revised": true, "final_justification": "The author clarified my main question in their response, and the authors responses to other comments helped clarify a lot of the confusion I had with the work.", "rebuttal": "We thank the reviewer for their careful reading and questions; we have grouped and will address each concern in order below. \n**In short:** we clarify notation, $\\beta$‑regimes, and provide new quantitative score‑alignment and intrinsic‑dimension numbers; and add stronger baselines on the MFEP task.\n\n### Clarity and significance\n\nWe understand that the notation can be a bit heavy and would like to point the reviewer to Appendix A.1, which contains a summary of the notation. We will also expand the introduction to DPA (Sec. 2.1) in the final work. The gist of it is that the DPA is an autoencoder trained with an **energy score** loss (Eq. (10)) that ensures *distributionally correct reconstructions*, so all data points that were encoded with the same encoding value are *guaranteed* to be distributed under this conditional of the *original data distribution* (the *ORD*) on reconstruction. \n\nWhile the use of the Energy score in deep learning is somewhat novel, similar models have been picking up considerable interest; for examples, see reference [20], De Bortoli et al. 2025 - “Distributional Diffusion Models with Scoring Rules”, and Shen and Meinshausen, 2024 -  “Engression\".  From a more general perspective, most unsupervised learning techniques **trade-off data distribution approximation with dimensionality reduction** (cf. lines 321-326). Here, we show that a single method *provably* achieves **both** *without* any targeted regularization, approaching PCA in the second aspect *while* being a flexible *non-linear* method.\n\nAs a practical example: because MFEP is defined by the **Stein score of the Boltzmann density**, a single latent that aligns with the score lets DPA recover the reaction pathway **in one pass**, whereas current approaches need expensive iterative refinement (lines 244-252).\n\n\n### Encoder optimality & $\\beta$ Regimes\n\nWe'd like to clarify that the optimal **encoder** is well-defined for $\\beta=2$ (lines `102-107`). $\\beta=2$ might induce the **only the decoder** to have a *non-unique* optimum, which would mean that the ORD is not a conditional of the original data distribution. This possible degeneracy has not been observed in any of the experiments performed. \n\n\nAs for the general case of non-optimality of neural networks due to training, we'd again like to point to the experiments as even approximately optimal (by necessity) encoder/decoder pairs used seem to validate the theoretical results. Furthermore, we believe that understanding the mechanism behind a method is valid even if the exact optimum might not be achievable with finite data/training time.\n\nFor a concrete illustration, we provide  **numeric results of the *cosine similarity*** between both sides of Eq. (7) - that is, the level set vector and the true data score - for the two examples in Fig. 1. The results were evaluated  on a 100 × 100 grid with points where the pdf exceeds 0.5 % of its maximum kept:\n\n| data set         | latent component | mean  cosine | std        | 95 %‑ile | points kept |\n| ---------------- | ---------------- | ------------ | ---------- | -------- | ----------- |\n| Standard Normal  | 0                | 1.000000     | 0.000000   | 1.000000 | 5 088       |\n| Standard Normal  | 1                | 1.000000     | 0.000000   | 1.000000 | 5 088       |\n| Gaussian Mixture | 0                | 1.000000     | 3.1 × 10⁻⁸ | 1.000000 | 4 729       |\n| Gaussian Mixture | 1                | 1.000000     | 3.0 × 10⁻⁸ | 1.000000 | 4 729       |\n\n\nFor the questions concerning $\\beta$: \n1. A general way to summarize Eq. (6) would be that a general push/pull relation holds (in expectation) that matches an expected \"flux\" (as it is an integral of a divergence) $J_\\beta$  on the LHS to the scaled variance of the level set on the RHS. The terms $J_\\beta$ and $S$ represent a \"push\" from the data density to deform the level set, while the \"variance\" term on the RHS \"pulls\" it closer as to minimize the level-set variance.\n   \n2. We have performed additional experiments for the conditional independence result using $\\beta=2$ (next section) and find that the Theorem *seems to hold in practice*. The requirement for $\\beta<2$ is a technical one for the proof, as the issue is the same discussed above with regards to the possible degeneracy of the optimal *decoder*. \n\n\n### Intrinsic Dimension Results \n\nWe'd like to point to **Sec. 1.2. of the supplement** which empirically demonstrates the validity of the intrinsic dimension result. The Supplement also explains the provided code. *Suppl. §1.2, Tab. 2‑3* thus show for the *noisy* S-curve dataset that the *conditional mutual information* `I(X;Z3|Z12)=0.13 bits` and the *conditional distance correlation* `dCorr=0.07`, meaning that the extraneous latent is effectively independent.\n\nFor many manifolds, the uninformative latents (denoted as $U$) become a *deterministic* function of the informative ones ($Z$), and satisfy Thm 3.4 directly. We present the following estimators:\n\n| Estimator           | Idea                                                                                                                       |\n| ------------------- | -------------------------------------------------------------------------------------------------------------------------- |\n| **R²**              | Fit a non-linear regressor $U\\approx g(Z)$ ; if $U=g(Z)$ deterministically the fit should explain ≈ 100 % of the variance. |\n| **ID‑drop**         | Compare intrinsic dimension of $(Z,U)$ vs. $Z$ via the Levina-Bickel MLE; $\\approx$ 0 implies no extra d.o.f. in $U$.      |\n| **$H(U \\mid Z)$** | Conditional entropy; large negative implies that $U$ tightly concentrated given $Z$.                                       |\n\nThe results for various datasets are presented below:\n\n| Dataset                   | Description                                               | $R^2$  | ID‑drop \\[2.5 %, 50 %, 97.5 %\\] | $H(U \\mid Z)$ \\[nats\\] |\n| ------------------------- | --------------------------------------------------------- | :----: | :-----------------------------: | :----------------------: |\n| Gaussian line             | 1-D line in $\\mathbb{R}^2$ with orthogonal Gaussian noise | 0.9997 |  0.0112 / **0.0122** / 0.0132   |          −7.259          |\n| Parabola                  | $(t,t^{2})$                                               | 0.9997 |  0.0046 / **0.0048** / 0.0051   |          −9.190          |\n| Exponential               | $(t, exp(t))$                                             | 0.9996 |  0.0045 / **0.0047** / 0.0050   |          −9.162          |\n| Helix slice               | $(\\cos t,\\sin ⁡t, t)$                                     | 0.9995 |  0.0041 / **0.0043** / 0.0045   |          −6.090          |\n| Grid sum                  | $z = x + y$                                               | 0.9986 |  0.0023 / **0.0029** / 0.0034   |          −2.759          |\n| S-curve                   | from `sklearn`                                            | 0.9996 | −0.0031 / **−0.0014** / 0.0000  |          −1.762          |\n| Swiss-roll (3-D)          | from `sklearn`                                            | 0.9872 | −0.0048 / **−0.0025** / −0.0003 |          −0.042          |\n| S‑curve **($\\beta$ = 2)** | ibid                                                      | 0.999  |     0.0030/ 0.0034/ 0.0037      |          ‑2.600          |\n\nFurthermore, for the Gaussian line we test Thm. 3.4 directly with a **Conditional randomization test**; under the null, we expect the p-values (100 replications) to be uniform. We confirm this with a K-S test:\n\n| K-S statistic $D$ | K-S p-value | p < 0.05 |\n| ----------------- | ----------- | -------- |\n| **0.061**         | **0.822**   | **4 %**  |\n\n\n### Further Experimental validation\n\nFirst, we'd like to address the relatively simple examples used for the experiments. We wish to provide visual intuition and verification of the theoretical results, which, for score alignment, is only possible for 2D examples with a known data density. The original DPA paper already includes high-dimensional and complex data examples that demonstrate DPA's *empirical* performance; here, we provide the *why*.\n\nThe score alignment is additionally illustrated in Sec. 1.1 of the Supplement. We  **add  more advanced autoencoder types** - $\\beta-VAE$ (Higgins et al. 2017) and  $\\beta-TCVAE$ (Chen et al. 2018) to *Tab 1.* below:\n\n| Model         | Param. comp.$\\bar z_{\\parallel}$ | Chamfer $\\bar{d}_{\\mathrm C}$ | Hausdorff $\\bar{d}_{\\mathrm H}$ | Mean 95th-perc. error | $\\bar{d}_{\\text{MFEP}\\rightarrow\\text{path}}$ | $\\bar{d}_{\\text{path}\\rightarrow\\text{MFEP}}$ |\n| ------------- | -------------------------------- | ----------------------------- | ------------------------------- | --------------------- | --------------------------------------------- | --------------------------------------------- |\n| $\\beta$-VAE   | 0.5 ± 0.51                       | 0.450 ± 0.288                 | 1.172 ± 0.512                   | 1.051 ± 0.477         | 0.539 ± 0.180                                 | 0.360 ± 0.462                                 |\n| $\\beta$-TCVAE | 0.375 ± 0.49                     | 0.377 ± 0.077                 | 1.378 ± 0.501                   | 1.228 ± 0.433         | 0.591 ± 0.180                                 | **0.164 ± 0.101**                             |"}, {"review_id": "5kfhhyeF3H", "reviewer": "Reviewer_kCFt", "summary": "The paper presented two theoretical results regarding the Distributional Principal Autoencoder (DPA). First, a closed-form relation linking each optimal level-set geometry to the data-distribution score explains, to some extent,  DPA’s empirical ability to disentangle factors of variation of the data. Second, the paper proved that if the data lies on a manifold that can be approximated by the encoder, latent components beyond the manifold dimension carries no additional information. The above two theoretical results are helpful in deepening the understanding of DPA.", "questions": "Can the authors provide a proof for the second result based on a more general definition of K'-parameterizability? As an example, is it possible to prove that data with an intrinsic dimension of K' (and satisfies a set of loose conditions) can always be effectively learned by DPA using K' latent parameters? This kind of results will significantly improve the applicability of Result Two.", "limitations": "N/A", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 4, "strengths_and_weaknesses": "Strengths: The results of the paper may arouse theoretical interests and are mathematics-inclined. \n\nWeakness: The generalization of the second major result in the paper relies on the definition of K'-parameterizability. However, the definition of  K'-parameterizability depends on DPA itself. That is to say, the K'-parameterizability here refers to the K'-parameterizability of DPA. This definition method obviously weakens the significance of this theoretical result.", "was_revised": true, "final_justification": "I maintain my previous assessment.", "rebuttal": "We thank the reviewer for the attentive review and for pointing out a possible generalization of our results.\nBelow we clarify what is already covered and what remains open.\n\n1. **K′‑parameterizable manifold always exists**: For every finite-dimensional dataset there is some smallest index $K′$ where $L^*_{k}$ stops decreasing; that K′ is “parameterizable” in Def. 3.1 with the trivial case  $K′=p$ . This is not necessarily connected to DPA per-se.\n   \n2. **K′‑best‑approximating encoder (Def. 3.3).** An encoder is **K′‑best-approximating** if  \n\t1. it achieves the global optimum among all K′‑latent models *and*  \n\t2. its *last* reconstruction term $L_{K′}[e,d]$ equals the optimal $L_{K′}^*$ .    \n\t   \n       When those two conditions hold, adding any surplus latent cannot lower the joint objective; strict propriety of the energy score then forces the surplus latents’ ORD to *repeat*  the optimal K′ ORDs, yielding the conditional‑independence conclusion in Theorem 3.4. Thus the applicability of Thm. 3.4 is connected to the **existence** of a **K′‑best-approximating** encoder.\n\n3.  **Exactly parameterizable manifolds (Corollary 3.5).** If the data lie on a K‑manifold and an encoder $e_{1:K}$ *exactly* parameterizes it, then $L_K[e,d]=0$ . Therefore any global K‑latent optimal encoder/decoder is automatically **K‑best-approximating** and the existence of a K′‑best-approximating encoder  is reduced to the optimality of the parameterizing encoder in this case. The existence of the parameterizing encoder is not too stringent of a requirement (cf. lines 191-193), and its optimality is explored in Remark 3.6.\n   \n4. A **further generalization** could thus require that for a certain class of $K$-manifolds and a certain class of neural networks used for DPA, there exists a a network achieving $L_{K}[e,d] \\approx 0$ . This would extend Cor. 3.5 and get around the existence issue, and is an open but feasible direction.\n\n\n *In short,* Theorem 3.4 already applies whenever a K′‑best-approximating encoder exists; Corol. 3.5 simplifies the latter under the clean “exact manifold‑mapping” scenario, and we think universal‑approximation approaches can be used to tackle the more general case."}, {"review_id": "63M2SFdknd", "reviewer": "Reviewer_NE6s", "summary": "In this submission, the author(s) investigate the theoretical properties of Distributional Principal Autoencoder (DPA): (1) They establish a closed-form relationship connecting optimal level-set geometry to the data-distribution score; (2) For data residing on a manifold approximable by the encoder, encoder components exceeding the manifold's dimensionality (or its optimal approximation) exhibit conditional independence from the data, thus conveying no supplementary information.", "questions": "1. Can the author(s) discuss the potential impact of the theoretical results?\n2. Can the author(s) provide some additional numerical results on the intrinsic dimension estimation? This is an important contribution of the work but there are no numerical demonstration on this point.\n3. Can the author(s) provide some discussions on their assumptions and how they are related to the existing literature? How strong they are?\n4. Can the author(s) discuss how the dimensionailty impacts the theoretical results? This is related to Q1. Because to my best, I note that many other generative learning techniques e.g. diffusion models performs quite well in the high-dimensional cases, and currently it is still unclear how the DPA performs in high-dimensional data setup.", "limitations": "See my comments on the weaknesses and questions.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths: The theoretical analysis is thorough and solid. The connection between DPA and the data-distribution score is interesting. At least for me, I haven't seen the similar result somewhere else.\n\nWeaknesses: It is unclear how the theory impacts the community. Although the author(s) provide an application of the link between DPA and the data-distribution score in the case of Boltzmann distribution, it is unclear why this application matters. In addition, it is unclear the performance of the intrinsic dimension estimation, which is related to the second contribution of the work.", "was_revised": true, "final_justification": "The author(s) have addressed all my concerns and I raise my score accordingly.", "rebuttal": "We thank the reviewer for their interest and informative feedback. We have grouped and will address each aspect of it in order below. \n### Impact\n\nMost approaches to unsupervised learning **trade-off data distribution approximation with dimensionality reduction** (cf. lines 321-326). Here, we show that a single method *provably* achieves **both** *without* any targeted regularization, approaching PCA in the second aspect *while* being a flexible generative model. The empirical findings (some additional ones presented here) further demonstrate that the theoretical findings can be observed even in practical experiments.\n\nWhile the use of the Energy score in deep learning is somewhat novel, similar models have been picking up considerable interest; for examples, see reference [20], De Bortoli et al. 2025 - “Distributional Diffusion Models with Scoring Rules”, and Shen and Meinshausen, 2024 -  “Engression\".\n\nAs a practical example: in molecular simulations, one typically wishes to find a compressed representation that can approximate the important, but *rare* transitions between physical states. Ideally, such a representation would parameterize the minimum free-energy path (MFEP), as traveling along the latter is the least-costly way for the system to undergo transitions. This representation can then be used to guide simulations. Because the MFEP is defined by the **Stein score of the Boltzmann density**, the fact that the latents align with the score lets DPA recover the reaction pathway **in one pass**, whereas current approaches need expensive iterative refinement (lines 244-252). This is further illustrated  in Sec. 1.1 of the Supplement (`supplement.pdf` in the provided archive). We furthermore **add more advanced autoencoder types** - $\\beta-VAE$ (Higgins et al. 2017) and  $\\beta-TCVAE$ (Chen et al. 2018) to *Tab 1.* below and demonstrate that DPA's score alignment is *provably* uniquely powerful for this scientifically important application:\n\n| Model         | Param. comp.$\\bar z_{\\parallel}$ | Chamfer $\\bar{d}_{\\mathrm C}$ | Hausdorff $\\bar{d}_{\\mathrm H}$ | Mean 95th-perc. error | $\\bar{d}_{\\text{MFEP}\\rightarrow\\text{path}}$ | $\\bar{d}_{\\text{path}\\rightarrow\\text{MFEP}}$ |\n| ------------- | -------------------------------- | ----------------------------- | ------------------------------- | --------------------- | --------------------------------------------- | --------------------------------------------- |\n| $\\beta$-VAE   | 0.5 ± 0.51                       | 0.450 ± 0.288                 | 1.172 ± 0.512                   | 1.051 ± 0.477         | 0.539 ± 0.180                                 | 0.360 ± 0.462                                 |\n| $\\beta$-TCVAE | 0.375 ± 0.49                     | 0.377 ± 0.077                 | 1.378 ± 0.501                   | 1.228 ± 0.433         | 0.591 ± 0.180                                 | **0.164 ± 0.101**                             |\n\n### Empirical Results on Intrinsic Dimension\n\nIn case this was missed, we'd like to point the reviewer to  **Sec. 1.2. of the supplement**  (`supplement.pdf` located at the root of the provided archive), which empirically evaluates the intrinsic dimension results. *Suppl. §1.2, Tab. 2‑3* thus show for the **noisy** S-curve dataset that the *conditional mutual information* `I(X;Z3|Z12)=0.13 bits` and the *conditional distance correlation* `dCorr=0.07`, meaning that the extraneous latent is effectively independent.\n\nWith DPA, for many manifolds the uninformative latents (denoted as $U$) become (in practice) a *deterministic* function of the informative ones ($Z$), and satisfy Thm 3.4 directly. Then the results can be verified using the following estimators:\n\n| Estimator        | Idea                                                                                                                       |\n| ---------------- | -------------------------------------------------------------------------------------------------------------------------- |\n| **R²**           | Fit a non-linear regressor $U\\approx g(Z)$ ; if $U=g(Z)$ deterministically the fit should explain ≈ 100 % of the variance. |\n| **ID‑drop**      | Compare intrinsic dimension of $(Z,U)$ vs. $Z$ via the Levina-Bickel MLE; $\\approx$ 0 implies no extra d.o.f. in $U$.      |\n| **$H(U\\mid Z)$** | Conditional entropy; large negative implies that $U$ tightly concentrated given $Z$.                                       |\n\nThe results for various datasets are presented below:\n\n| Dataset                   | Description                                               | $R^2$  | ID‑drop \\[2.5 %, 50 %, 97.5 %\\] | $H(U\\mid Z)$ \\[nats\\] |\n| ------------------------- | --------------------------------------------------------- | :----: | :-----------------------------: | :-------------------: |\n| Gaussian line             | 1-D line in $\\mathbb{R}^2$ with orthogonal Gaussian noise | 0.9997 |  0.0112 / **0.0122** / 0.0132   |        −7.259         |\n| Parabola                  | $(t,t^{2})$                                               | 0.9997 |  0.0046 / **0.0048** / 0.0051   |        −9.190         |\n| Exponential               | $(t, exp(t))$                                             | 0.9996 |  0.0045 / **0.0047** / 0.0050   |        −9.162         |\n| Helix slice               | $(\\cos t,\\sin ⁡t, t)$                                     | 0.9995 |  0.0041 / **0.0043** / 0.0045   |        −6.090         |\n| Grid sum                  | $z = x + y$                                               | 0.9986 |  0.0023 / **0.0029** / 0.0034   |        −2.759         |\n| S-curve                   | from `sklearn`                                            | 0.9996 | −0.0031 / **−0.0014** / 0.0000  |        −1.762         |\n| Swiss-roll (3-D)          | from `sklearn`                                            | 0.9872 | −0.0048 / **−0.0025** / −0.0003 |        −0.042         |\n| S‑curve **($\\beta$ = 2)** | ibid                                                      | 0.999  |     0.0030/ 0.0034/ 0.0037      |        ‑2.600         |\n\nOne can also test Thm. 3.4 directly. For the Gaussian line example, we use a **Conditional randomization test**; under the null, we expect the p-values (100 replications) to be uniform. We confirm this with a K-S test:\n\n| K-S statistic $D$ | K-S p-value | p < 0.05 |\n| ----------------- | ----------- | -------- |\n| **0.061**         | **0.822**   | **4 %**  |\n\nWe will consider moving some of these results to the main paper in the final version, and are performing further experiments.\n\n### Assumptions\n\nThe assumptions we make can be summarized as follows:\n\n- **Global optimality of encoder/decoder:** while this is a relatively strong assumption, it is not uncommon in theoretic treatments of  deep learning models.  Unlike, e.g., ref. [1], our results are **non-asymptotic and exact** given the optimum (lines 272-285). The analysis of optimal solutions can reveal important insights, even if in practice one can only approximate them. Furthermore, our empirical experiments corroborate the theoretical findings despite only approximate optimality of the models. \n  \n  For example, we provide **numeric results of the *cosine similarity*** between both sides of Eq. (7) - that is, the level set vector and the true data score - for the two examples in Fig. 1. The results were evaluated  on a 100 × 100 grid with points where the pdf exceeds 0.5 % of its maximum kept:\n\n| data set         | latent component | mean  cosine | std        | 95 %‑ile | points kept |\n| ---------------- | ---------------- | ------------ | ---------- | -------- | ----------- |\n| Standard Normal  | 0                | 1.000000     | 0.000000   | 1.000000 | 5 088       |\n| Standard Normal  | 1                | 1.000000     | 0.000000   | 1.000000 | 5 088       |\n| Gaussian Mixture | 0                | 1.000000     | 3.1 × 10⁻⁸ | 1.000000 | 4 729       |\n| Gaussian Mixture | 1                | 1.000000     | 3.0 × 10⁻⁸ | 1.000000 | 4 729       |\n- **Encoder’s function class is sufficiently expressive:** We assume the encoder can approximate the data manifold. This is akin to a **universal approximation** assumption restricted to the manifold structure. This is not unheard of, and has been empirically evaluated (lines 191-193, [5]).\n  \n- **Data lies on (or near) a low-dimensional manifold:** This is a common assumption to all representation learning. If this is not the case, Theorem 3.4 does not apply, but then again such a problem is not suitable for dimensionality reduction techniques. This does not impact Sec. 2.\n  \n- **Jacobian full row rank almost everywhere on the level set:** We would like to point out that:\n  1. This only applies to Thm. 2.6\n  2. If it's violated for a specific level set, the Theorem is simply **silent** on that specific level set; that is, this is not a pre-condition of the theorem.\n   Further discussion is provided in lines 108-112; in our experiments, we observe level sets that are well-behaved for models whose training has converged (cf. the above alignment result).\n\n\n\n### Data Dimensionality and High-Dimensional Performance\n\n\nAll theoretical results presented hold for **any** intrinsic/ambient dimensionality combination. In order to provide visual intuition, we opted for simple examples with a known data density (in case of the score alignment).\n\nThe original DPA paper already includes high-dimensional and complex data examples (e.g., MNIST, climate data) that already demonstrate DPA's strong *empirical* performance; our analysis *explains* it.\n\nIn ordinary diffusion models, the latent dimension of the process is equal to the ambient one. While they learn the score, DPA thus additionally learns a *compressed representation*, for which we prove is *maximally informative* in case the data *can* be compressed. Thus we show that the method combines two desirable properties that are not found together in other models."}, {"review_id": "Ys9aJUJtje", "reviewer": "Reviewer_YCG7", "summary": "This paper provides a rigorous theoretical analysis of the Distributional Principal Autoencoder (DPA), a model that combines principled distributional reconstruction with interpretable latent representations. The authors derive an exact closed-form relation connecting the encoder’s optimal level sets to the data distribution’s score function, offering a geometric understanding of how the model captures the data structure. They further prove that for data lying on a low-dimensional manifold, the latent dimensions beyond the manifold’s intrinsic dimension become conditionally independent and carry no additional information, enabling precise intrinsic dimension estimation without extra regularization. The theory is supported by experiments on synthetic distributions and the Müller-Brown potential.", "questions": "$\\textbf{Questions for the Authors:}$ \n\n1) The theoretical analysis assumes that the encoder and decoder achieve global optimality in minimizing the DPA loss (e.g. Definition 3.1., and the related theoretical results Proposition 3.2, Theorem 3.4). Given the non-convexity of neural network training and the common presence of local minima, how realistic is this assumption in practice?\n\n2) The theory requires the data manifold to be exactly or approximately parameterizable by the encoder’s function class. How strong is this assumption, and how does it hold up when modeling complex or noisy real-world data?\n\n3) The experiments primarily involve low-dimensional synthetic data and the Müller-Brown potential. How well do these empirical findings generalize to complex, high-dimensional datasets, such as MNIST, climate simulations, or single-cell genomics, which were explored in the original DPA paper?\n\n4) The main results rely on population-level optimality, assuming infinite data and exact minimization of the DPA loss. Could the authors clarify theoretically or empirically how finite samples and imperfect training affect the level set alignment (Theorem 2.6) and the disentanglement of extraneous dimensions (Theorem 3.4)?\n\n5) The results assume that the encoder’s Jacobian has full row rank almost everywhere on the level sets. How restrictive is this assumption in practice?", "limitations": "Yes.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "$\\textbf{Strengths}$\n\n$\\textbf{Quality:}$\nThe paper presents mathematically precise, non-asymptotic results connecting the DPA encoder’s level sets to the data score function. The derivations provide a solid framework explaining the DPA’s latent space geometry. Additionally, by proving that extraneous latent dimensions are conditionally independent, the work offers a novel formalization of intrinsic dimension estimation without extra regularization or architecture changes.\n\n$\\textbf{Clarity:}$\nThe paper clearly defines key concepts. The logical flow from assumptions to lemmas, theorems, and corollaries is coherent and supports the main narrative effectively.\n\n$\\textbf{Significance:}$\nThe connection between encoder geometry and the data score may inspire new generative model architectures and more interpretable latent variable models.\n\n$\\textbf{Originality:}$\nWhile building on the existing DPA model, the exact characterization of level sets in terms of score functions and conditional independence results for extraneous dimensions appear original.\n\n\n$\\textbf{Weaknesses}$\n\nThe theoretical guarantees provided by the DPA rest on several strong assumptions that may limit their applicability in practice. First, the results require that the encoder and decoder achieve global optimality in minimizing the DPA loss, an assumption that is difficult to satisfy due to the non-convex nature of deep network training and the potential for convergence to suboptimal local minima. Second, the encoder function class must be sufficiently expressive to realize an exact parameterization of the underlying data manifold, which may not hold for networks with limited capacity or when faced with highly complex or irregular data structures. Finally, the data distribution is assumed to lie exactly or approximately on a low-dimensional manifold, an idealization that real-world datasets often violate due to noise, outliers, or intrinsic complexity. \n\nConsequently, while the theoretical framework offers valuable insight into the capabilities of DPA under ideal conditions, translating these guarantees into practical success requires careful consideration of these assumptions and further empirical validation. The empirical validation is also limited to low-dimensional synthetic examples and a classical benchmark (Müller-Brown potential), which may not fully demonstrate the model’s robustness or scalability to more complex, high-dimensional real-world data. For example, the original DPA paper evaluated the model on diverse datasets including MNIST, synthetic disk images, regional and global climate data, and single-cell genomics, demonstrating its applicability across varied domains and data complexities.", "was_revised": true, "final_justification": "The discussion on global optimality and manifold parameterizability helped me clarify my earlier misunderstanding regarding the assumptions behind Theorem 3.4 and Corollary 3.5.", "rebuttal": "We thank the reviewer for their questions; we have grouped and will address them in order.\n\n**In short:** we clarify possible misunderstandings, provide new quantitative score‑alignment and intrinsic‑dimension results, and add stronger baselines on the MFEP task; together, these demonstrate that the theoretical findings transfer to applications.\n\n### Global optimality\n\nWhile this is a relatively strong assumption, it common to many theoretic investigations of deep learning models. The goal is to explain the performance of DPA, which was already demonstrated in the original paper on complex data; in practice, our experiments show that even with approximate optima, one can observe the presented theoretical properties *and* that they lead to very desirable outcomes.\n  \nFor example, we provide **numeric results of the *cosine similarity*** between both sides of Eq. (7) - that is, the level set vector and the true data score - for the two examples in Fig. 1. The results were evaluated  on a 100 × 100 grid with points where the pdf exceeds 0.5 % of its maximum kept:\n\n| data set         | latent component | mean  cosine | std        | 95 %‑ile | points kept |\n| ---------------- | ---------------- | ------------ | ---------- | -------- | ----------- |\n| Standard Normal  | 0                | 1.000000     | 0.000000   | 1.000000 | 5 088       |\n| Standard Normal  | 1                | 1.000000     | 0.000000   | 1.000000 | 5 088       |\n| Gaussian Mixture | 0                | 1.000000     | 3.1 × 10⁻⁸ | 1.000000 | 4 729       |\n| Gaussian Mixture | 1                | 1.000000     | 3.0 × 10⁻⁸ | 1.000000 | 4 729       |\n\nSimilarly, The conditional independence results were presented in **Sec. 1.2. of the supplement**  (`supplement.pdf` located at the root of the provided archive). *Suppl. §1.2, Tab. 2-3*  show for the **noisy** S-curve dataset that the *conditional mutual information*  and the *conditional distance correlation* indicate conditional independence even with the presence of noise and imperfect training.\n\nFurthermore, the stability of the results (across numerous random seeds) presented in Sec. 1.1 of the supplement indicate that finding a *good* optimum that displays the score alignment property is *not at all difficult* for DPA. We add further comparisons to two more advance autoencoder models: - $\\beta-VAE$ (Higgins et al. 2017) and  $\\beta-TCVAE$ (Chen et al. 2018) below as an extension to Tab 1. in the Supplement:\n\n| Model         | Param. comp.$\\bar z_{\\parallel}$ | Chamfer $\\bar{d}_{\\mathrm C}$ | Hausdorff $\\bar{d}_{\\mathrm H}$ | Mean 95th-perc. error | $\\bar{d}_{\\text{MFEP}\\rightarrow\\text{path}}$ | $\\bar{d}_{\\text{path}\\rightarrow\\text{MFEP}}$ |\n| ------------- | -------------------------------- | ----------------------------- | ------------------------------- | --------------------- | --------------------------------------------- | --------------------------------------------- |\n| $\\beta$-VAE   | 0.5 ± 0.51                       | 0.450 ± 0.288                 | 1.172 ± 0.512                   | 1.051 ± 0.477         | 0.539 ± 0.180                                 | 0.360 ± 0.462                                 |\n| $\\beta$-TCVAE | 0.375 ± 0.49                     | 0.377 ± 0.077                 | 1.378 ± 0.501                   | 1.228 ± 0.433         | 0.591 ± 0.180                                 | **0.164 ± 0.101**                             |\n\nWe will consider moving some of these results to the main paper in the final version, and are performing further experiments.\nIn short, all experiments point to the *observability of the derived theoretical properties in practice*.\n\n### Manifold Parameterizability Assumption\n\n**Encoder capacity:**\n> “…the encoder function class must be sufficiently expressive to realize an exact parameterization of the underlying data manifold, which may not hold …”\n\nWe apologize for any confusion—our Theorem 3.4 does **not** assume the encoder can already parameterize the manifold **exactly**.\n**Definition 3.1** introduces *K′‑parameterizability* **of the data manifold**, not of the network. It states only that there exists some latent width K′ at which the **best attainable (nested) loss term stops decreasing**. This is always the case for any finite-dimensional manifolds, with $K=p$ in the \"worst\" case.\n\n**Exact manifold maps are merely a special case:** *Corollary 3.5* shows that **if** an encoder exactly parameterizes a K‑manifold **and** is globally optimal among K‑latent models, then the K dimensions are *automatically* K′‑best with K′=K. This is a *sufficient* but **not necessary** scenario. *Theorem 3.4* *still applies* whenever a K′‑best-approximating encoder exists—even if $L_{K′}^*$ is positive; that is if:\n\t1. the data are *noisy* and the  cannot be learned as an exact manifold\n\t2. The manifold can only be approximately learned due to its *complexity*\nThen, the latents beyond $K'$ capture, at most, repetitions of that irreducible error.\n\nWe wish to thank the reviewer for prompting this clarification and will add a sentence further distinguishing \"*K′‑parameterizable data*\" (a definition) from *K'-best-approximating* encoder (a requirement for Thm. 3.4) and the *exactly-parameterizing* encoder (Cor. 3.5, a special case).\n\nThus, in practice for a very complex data, DPA might use slightly more latents $K' > K$ to \"explain\" the noise (or complexity); however, past $K'$ the latents will still be uniformative. It is unlikely that any other dimensionality reduction technique would perform much better in that case as the limitation is on the neural network as such.\n\n\n### Generalization to More Complex data \n\nAll theoretical results presented hold for **any** intrinsic/ambient dimensionality combination. In order to provide visual intuition, we opted for simple examples with a known data density (and hence the score). For example, the Muller-Brown example is an example with a known score that nonetheless demonstrates a large practical impact of the findings, as the method could be used to significantly speed up simulations.\n\nThe scalability question is largely resolved by **prior empirical work**:  the original DPA paper already includes results for e.g.,  MNIST, climate data, etc., demonstrating competitive reconstruction and an automatically inferred latent dimension of $\\approx$ 32 for MNIST ([17], Fig. 7). Our contribution is to **explain why** the method performs so well. For example, the global precipitation field results ([17], Fig. 9), which represent the periodic `month` dimension in a polar representation can be explained by our results (cf. lines 212-219).\n\n\n### Finite Samples & Imperfect Training\n\nWe believe that the experiments presented already demonstrate robustness to finite sample sizes and imperfect training.\n\nFor a more theoretical guarantee:\n1. The DPA decoder is an *Engression* network. The original paper provides explicit finite‑sample error bounds (Theorem 3, Shen & Meinshausen 2024).\n2. The DPA encoder is a generic neural network, so standard results apply there.\n\n### Full-Row-Rank Assumption\n\n1. This assumption is **only invoked for Theorem 2.6**\n   \n2. If it's violated for a specific level set, the Theorem is simply **silent** on that specific level set; that is, this is not a pre-condition of the theorem.\n   Further discussion is provided in lines 108-112.\n3. In our experiments, we observe level sets that are well-behaved for models whose training has converged (cf. the above alignment result)."}]

Author rebuttal:

Meta-review: The paper studied the properties of distributional principal autoencoder (DPA) and provided solid theoretical analysis. The reviewers value the insightful theoretical results, but also raised questions regarding derivation, empirical evaluation and the practicality of the results. The authors provided additional results and clarified most of the concerns. The rebuttal successfully convinced two reviewers to raise their scores. Now, all reviewers share a positive view about the paper.  The AC agrees with the reviewers’ assessments and feels enthusiastic about the potential of the paper. The final version should incorporate the reviewers’ comments and also provide more clear clarifications as provided in the rebuttal.
Decision: Accept (poster)