Pareto Variational Autoencoder
Authors: Mincheol Cho
Yedarm Seong
Joong-Ho Won
Venue: ICLR 2026 Poster (2026)
Keywords: Variational Autoencoder
Symmetric Pareto distribution
Information geometry
Heavy-tail learning
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)

Abstract: This paper
introduces a new class of multivariate power-law distributions---the symmetric Pareto (symPareto) distribution---which can be viewed as an $\ell_1$-norm-based counterpart of the multivariate $t$ distribution,
with the motivation of capturing the heavy tail of the target distribution in generative modeling and bringing robustness to noise in downstream tasks such as image denoising. 
The symPareto distribution possesses many attractive information-geometric properties with respect to the $\gamma$-power divergence that %naturally %\red{characterizes the geometric structures of power-law families.}
is a natural alternative to the Kullback-Leibler divergence, the core of the conventional variational autoencoder (VAE) models, for power families.
Leveraging on the joint minimization view of variational inference, this paper proposes the ParetoVAE, a probabilistic autoencoder that minimizes the $\gamma$-power divergence between two statistical manifolds.
ParetoVAE employs the symPareto distribution for both prior and encoder, with flexible decoder options including multivariate $t$ and symPareto distributions. 
Empirical evidences demonstrate the effectiveness of ParetoVAE across multiple domains through varying the types of the decoder. 
The $t$ decoder achieves superior performance in sparse, heavy-tailed data reconstruction and word frequency analysis; the symPareto decoder enables robust high-dimensional denoising.

Published as a conference paper at ICLR 2026
PARETO VARIATIONAL AUTOENCODER
Mincheol Cho Yedarm Seong Joong-Ho Won
'Department of Statistics, Seoul National University Intterdisciplinary Program in Artificial Intelligence, Seoul National University (code 1 4 78, mybirth0407, won. j) @ ac.kr
ABSTRACT
This paper introduces new class of multivariate power-law distributions- sym- metric Pareto (symPareto) distribution which can be viewed as an (--oorm-based counterpart of the multivariate t distribution, with the motivation of capturing the heavy tail the target distribution in generative modeling and bringing robustness to noise in downstream tasks such image denoising. The symPareto distribution possesses many attractive information -ooomationggommetiicc properties with respect to the y-power divergence that is natural alternative the Kullback-Leibler divergence, the core of the conventional variational autoencoder (VAE) models, for power fam- ilies. Leveraging the joint minimization view of variational inference, this paper proposes the Pareto VAE, a probabilistic autoencoder that minimizes the "-power divergence between two statistical manifolds. Pareto employs the symPareto distribution for both prior and encoder, with flexible decoder options including multivariate t and symPareto distributions. Empirical evidences demonstrate the effectiveness of ParetoVAE across multiple domains through varying the types of the decoder. The / decoder achieves superior performance in sparse, heavy-tailed data reconstruction and word frequency analysis; the symPareto decoder enables robust high-dimensional denoising.
INTRODUCTION
For over decade since its introduction, the Variational Autoencoder (VAE, Kingma & Welling 2013)
has stood the test of time as a fundamental probabilistic generative model for scalable probabilistic inference and representation learning. Conventional VAEs typically employ exponential-fanm family distributions, most notably multivariate Gaussian, for their probabilistic model due to their math- ematical tractability. However, the exponential tail assumption often falls short when modeling real-world data that exhibit heavy tails and extreme events. For example, scale-free network degree distributions (Barabási & Albert 1999) and long-tailed class-frequency distributions (Zhang et 2023) challenge models built on light-tail assumptions. In particular, the VAE literature has pointed out that Gaussian-based VAEs tend to underestimate tail probabilities, over-regularize latent codes, and fail to capture rare but informative events (Lafon al.
2023 Kim et al. 2024).
To address these limitations, recent work Takahashi et al. 2018 Abiri & Ohlsson 2020 Kim et al. 2024) has explored VAEs based multivariate Student's 's ditrrbbtion, yet the latter constitutes only one option among many heavy-tailed distribution families. In particular, classical extreme-value theory suggests that univariate Pareto distribution-basse modeling is well suited to capturing tail and extreme behavior (Balkema & de Haan 1974 Pickands) 1975). Motivated by this, we introduce a new multivariate generalization of the Pareto family-the the symmetric Pareto (symPareto) distribu- tion-and and corresponding heavy-tail VAE framework, named ParetoVAR VAE. ParetoVAD incorporates symPareto distributions key components, overcoming the representational limitations of Gaussian VAEs in satreme-valuee scenarios and inducing sparsity in the latent embedding space through the inherent E formulation.
The use of symPareto distributions in conventional VAEs raises computational challenges, especially for the evidence lower bound (ELBO) estimation. Specifically, the Kullback-Leihle (KI.) divergence between two symPareto distributions lacks a closed-form expression, necessitating numerical integra-
Equal contribution.

---

Published as a conference paper at ICLR 2026
Table 1: Comparison of heave- VAE variants. Each tuple denotes (reconstruction loss, latent
regularization), where the components are either { or squared lg (E)
Latent distribution (Encoder Prior)
Decoder distribution (Generator)
Student's t
SymPareto
Student's SymPareto
Cg 2g)
Kim al. (C,,
2024)
(6g, 4 (L,, 81)
tion that becomes computationally prohibitive as the dimension grows, or Monte Carlo methods that cost extra variability. We address this issue by adopting an information-geomeet framework based on the y-power divergence, serving as a tractable alternative to the KL divergence. This reformulation admits closed-form expressions between power-law distributions, thereby enabling efficient optimization for modeling heavy-tailedness.
joint minimization
This work naturally extends the family of heavy-tailed variants of VAEs from Student's distri- butions to symPareto. In particular, we generalize the VAE (Kim et al. 2024) structure into a symPareto-based formulation. Table provides tion combinations and their corresponding reconstruction loss/regularization structures (either & or ()). It shows how our proposal broadens the /-based heavy-tailed VAE family a symPareto-based classes. Notably, incorporating symPareto distributions induces an & term in the objective, featuring sparsity and robustness in heavy-tailed VAE models.
provides an overview of the resulting latent/decoder distribu-
Related Work Recognizing the limitations of Gaussian distributions in standard VAEs, many researchers have proposed approaches that modify the ELBO by altering the prior, encoder, or decoder. On the prior side, spike-and-slab priors (Tonolini al. 2020), Dirichlet priors (Joo al. 2020). and tilted-Gaussian variants (Floto et al., 2023) have been explored as non-Gaussian alternatiess. distinct line of work incorporates heavy-tailed distributions directly into the VAE components. Student's 's & decoders improve robust density estimation (Takahashi et al. (2018;; heavy-tailed radii and conditional angles target multivariate extremes (Lafon et al. 2023); conditional t-based decoders enhance cass-conditional generation quality (Bouayed et al. 2025): and spherical Cauchy encoders combined with a uniform prior effectively capture directional structure (Sablica & Hornik 2025).
Beside modifying the model distributions, other VAE approaches focus replacing the divergence measure itself. For example, Rényi's a-divergence (Li & Turner 2016), geo-geeomttric Jensen- Shannon divergence (Deasy al.) 2021). and 3-divergence (Akrami et al. (0022) have been used in VAEs for specific purposes such as robust estimation. Beyond alternative divergences, some approaches directly reformulate the ELBO via joint minimization problem, such as Info VAE (Zhao et al. 2019) with mutual-information regularization and Coupled VAE (Hao & Shafto 2023) with entropic optimal transport. Notably, the VAE (Kim et al. 2024) introduces joint minimization framework that leverages the 7-power divergence and replaces the prior, encoder, and decoder distributions with multivariate distributions, based information geometric principles (Amari
2016 Eguchi 202T)
Other generative models, especially those based on the generative adversarial network (GAN) framework, have also been developed to explicitly capture heavy-tail behaviors by changing the prior distribution or divergence. Some GAN-based models use power-law priors, such as univariate & (Feder et al. 2020) and Pareto (Huster et al. (021) distributions. Others apply some exotic divergences, such as the (Laszkiewicz et al., 2022. Hickling & Prangle| 2024) and diffusion models (Pandey et al. 2024 Lian al.] 2025) have garnered attention for their ability to handle rare or extreme data points by moving beyond the Gaussianity assumption.
a-divergence (Chen al.) 2024). Furthermore, normalizing flows
2 THEORETICAL BACKGROUND
2.1 VARIATIONAL AUTOENCODER (VAE)
VAE aims to approximate the true data distribution Pamaa (x) by modeling the marginal likelihood (x) po (=( dz, where represents latent variable. Due the intractability of
2

---

Published as a conference paper at ICLR 2026
directly computing the posterior po(2|æ), VAE employs variational inference to approximate this distribution. The VAE framework consists of two primary components, encoder (inference model) gs(2|z)) FF and decoder (generative model) Pe (a|z) G, where F and & are spaces of density functions parameterized by deep neural networks. The encoder approximates the posterior po (z|æ), while the decoder generates sample x conditioned the latent variable
The primary objective is to minimize the KL divergence ( II- between the true posterior pa (=(2ar) and its approximation qo (x(22) (encoder). This is equivalent to maximizing the evidence lower bound (ELBO) of the marginal lgglikelihood, expressed as
BLBOg, (x) (- (lo DxL (1o(= (l
Here pz(=) refers the prior density, which is commonly assumed be standard m-variate Gaussian Pz(z) ~ (0, I). Typically, both the encoder and decoder are also chosen as multivariate
Gaussian
(=(=/æ) Ni (up(a) (x)x),
pa(a|=)
with diagonal encoder variance [x(x) ( (z)) simplify computation. As result of the Gaussianity assumption, the objective function of the VAE can be written as the sum of the mean squared error (MSE, Eg
20 -|>~ ||>--((==]|
and a KL regularization term.
VAE as Joint Minimization Problem VAE objective can be viewed as joint minimization problem between two statistical manifolds HHan et al.) (020)). The model distribution manifold Mmodel and the data distribution manifold Masta are defined as
Mmodel = {pe(z, (z, 2) = pe (=( θ}, Maata = {go(x, 2) = Paata (z)0(()))))) e Φ}.
The KL divergence between these manifolds can be expressed follows:
DKL (40 || pe) = -Ez-pan ELBO,, (z)] H (Pdat);
where H is the differential entropy of the data distribution. Therefore, maximizing the ELBO can be reformulated as the following joint minimization problem:
(po- {e-)
ME DKI p))
This perspective allows us to interpret VAE as minimizeing the KL. divergence between Mmodel and Másta- which can be solved by the ieformation-geomentii em-algorithm (Csiszár) 1984 Han et al. 2020). The key point is that, when constructing VAE structures, one can bypass an explicit ELBO maximization and instead directly solve joint minimization problem. Moreover, this view naturally motivates replacing KL divergence with alternative divergences to move beyond Gaussianity.
2.2 DISTRIBUTIONS IN POWER-LAW FAMILIES
random variable follows a (univariate) power-law distribution if its probability density or mass function p(x) is proportional to - 0 (Newman 2005). We introduce some multivariate extensions of paw distributions relevant to heavy-tarled VAE models.
Multivariate distribution The n-variate Student's f distribution with degrees of freedom has its density function as
tn (2
Tnw
1
r 11
T
I' uoto
1
(VT) T
where := IT ore Here y E R and σ E R denote the location and scale parameters, respectively, and the division (r p)/a is to be understood element-wise, with slight abuse of notation. It is widely regarded as a heavy- generalization of the multivariate Gaussian, since converges to Nr, (j4, r²T,))
3

---

Published as a conference paper at ICLR 2026
P.0. 1.6)
1400, 1.6)
N1 f
OCDF (lgg acie)
109
10
Radius
10
(b)
(a)
Figure 1: (a) 2D scatter plots of SymPareto (left), Student's t (middle), and Gaussian (right). (b)
Log-scale CCDF comparison of SymPareto with varying p.
Multivariate Symmetric Pareto distribution Although several formulations of the multivariate Pareto distribution have been proposed (Aroll et al. 1993 Hanagal 1996; Rootzén et al. 201 7)).
most lack explicit density representations suitable for the computation of divergence and are restricted to positive supports. Motivated by Mardia's 's multivariate Pareto of the first kind (Mardia 1962), we introduce the multivariate symmetric Pareto distribution (symPareto). Its density is defined
Pa(x (,,,N)) Cn. σ
(1 x )
1
-v++n)
I( + n) (20))
CR.PI.P
(3)
This distribution can be regarded heavy-tailed analogue of the product of univariate Laplace distributions, or an fg-norm-basee counterpart of the multivariate t distribution. Further properties of the symPareto distribution are provided in Appendix
Figure [illustrates characteristic behaviors of the proposed symPareto distribution. Figure la displays scatter plots of 5,000 samples drawn from the symPareto, Student's 4, and Gaussian distributions in two dimensions. Most Gaussian samples concentrate within radius S 3 due to the exponentially decaying tail, whereas Student's t produces noticeably more samples beyond this range. The symPareto distribution exhibits even heavier tails, with many extreme samples lying outside radius or 10. Moreover, the samples of symPareto tend to align the coordinates, producing a cross-like shape. This reflects the l fnorr structure in its density.
Figure Lb shows the complementary cumulative distribution functions (CCDFs) the symPareto distribution log scale for various values of For small v, the CCDF exhibits polynomial decay, indicating heavy tails. As v increases, the tail becomes lighter, and the limit v -> oc the symPareto converges the Laplace distribution, which decays exponentially.
2.3 INFORMATION GEOMETRY AND POWOER DIVERGENCE
In this section, we briefly introduce the 7-power divergence an alternative to the KL divergence from an information peermaation-gommerii perspective. More formal and detailed exposition these concepts can be found in Amari (016)).
Divergence in information geometry In information geometry. family of probability distributions is regarded as a statistical manifold S, and specific distribution in the family is point on S. To quantify the discrepancy between two points on S, one typically employs differentiable function D S x S [0, oc). called divergence. A divergence not only measures separation of points but also induces geometric structures on S, including an affine connection that governs how tangent vectors transported between nearby points. The connection V in turn defines V-geodesics, i.e., one-dimensional submanifolds that generalize the notion straight lines on the statistical manifold.
Formally, each V-geodesic associated with V. In this sense, the manifold S is said to V-flat. This property enables powerful tools of geometric analysis, such joint minimization algorithms (Csiszár 1984 Han et al. 2020).
be expressed linearly with respect to affine coordinate systems
4

---

Published as a conference paper at ICLR 2026
q
(0,0z)
(ps
(A.01) Pb 90
Ga, m))
Sy (μ 0)
Figure 2: Illustration two symPareto manifolds with different μ and their translation to the 7-flat manifold S, for analyzing the "y-power divergence.
KL divergence well-known example of divergence is the KI divergence. One of the corre- sponding geodesics is the e-geodesic (e for exponential), defined by
(x:) & exp((1 (og)(æ) + Llo (())), L 10.
In particular, e-geodesics characterize distributions in the exponential families, implying that these families form e-flat submanifolds. Thus, the KL divergence naturally induces geometric structure with respect exponential families. However, not all statistical manifolds exhibit such family. For instance, for distributions in the power-law family, such as the t and generalized Pareto distributions, the KL divergence does not induce dually flat geometric structure.
y-power divergence A divergence suitable for the power-law families is the y-power divergence Kim et al. 2024), generalization of the notion due to Eguchi (202T). It is defined as
Dy > 7 H4 (-1, \
Hp := (--p))==1111::: 11lli (/ p( ((o))ttdd))
+ (,(,,) = :=
p(x) |P
((x))
dx.
where (.) and C..). are called the y-power entropy and y-power cross-entropy, respectively. This divergence parallels the e-geodesic for the KI. divergence by introducing the 7-power geodesic, which induces 7-flat submanifolds as
S = {pa(x) x (1 (1T s(x)) 00
(4)
where <xx) the sufficient statistic of the given distribution. When the symPareto has / s(x) |x| E Ra is valid sufficient statistic. In this case, the symPareto density can be expressed in the f-lt form (4 with y = 24n
Upper bound of 7-power divergence for flatness However, when μ * 0, the symParcto distri- bution no longer admits an sufficient statistic s(x) of the required form even though the 7-power divergence between two symPareto distributions with different μ be written in closed form. Consequently, 7- flatnness is not preserved in the noncentral symPareto distributions.
To handle this case, we consider translation of the distributions. By shifting the location parameters Jn, 122 to the origin, the corresponding zero-centered distributions Po = Pn, (0, σ1, v) and go = P,, o0, v) lie the 7-flat manifold S., with 7 _1. Figure
2 illustrates two manifolds of the symPareto distribution with different μu, and their translation to the f-at manifold Sy, We then obtain an upper bound of D,, || q) by y-power divergence between Pn and go together with an additional term reflecting the translation cost in terms of the E -ddittncce between μir and μ2 relative the scale.
Theorem 2.1 (Upper bound for p-power divergence between two noncentral symPareto distributions). Let 2 1, v > 1, and 7 = vtn' For Po = P., v) and go = (0, a2 v), the 7-power divergence between Po and is finite when > and can be expressed in closed-form
---Ti= +î (1
1
))
- -
V 1ー (u+ ))
Dx || go) = C²
(5)
5

---

Published as a conference paper at ICLR 2026
where C == (n, - (v
- n.u-1vv Furthermore, for two noncentral symPareto distribu-
tions Pn (1 v) and P, (12, 02, 02, v), the divergence D, (p q) has an upper bound
Dy(p11g) - (20 qo) +
M1
12
(6)
02
with 3 (1 9) C n.--1, 02 - In 6 equality holds when μ1 P2-
THE PARETO VARIATIONAL AUTOENCODER
3.1 THE PARETOVAE STRUCTURE
In this section, we introduce the Pareto VAE framework and compute the corresponding 7-loss objective via the 7-power divergence. All detailed derivations and proofs are in Appendix B
Prior, decoder, and encoder Motivated by ( we set a minimization problem for VAE as mini- mizing the 7-power divergence between two joint manifolds:
(ge-, pa-) = argmingee, de Dy pe).
(7)
where (x 2) = go (x) and pa(a, 2) = pa 2 p(2). Under this problem setup, the construction of Pareto VAE starts from the heavy-tailed joint decoder model:
2(v 2irlm))
|| 1 ==(=) ||=||
(o(x, x 1
(8)
o²
From g
deduce the prior p(=) and decoder (2r =):
p(2) Pm 0, 1m, v), po(2 2) = in x uo(2, + ||
In +
2(v m)
For the encoder, we also adopt a symPareto distribution consistent with the prior, while increasing the degree of freedom by reflect the contribution of the data dimension:
4e (a x)) = (= (o(x), (x), n/2).
Reparameterization Trick for symPareto To enable gradient-based learning, the VAE framework employs the reparameterization trick (Kingma & Welling 2013). Exactly as Student's t-distribution can be represented in terms of Gaussian and chi-squared variables, symPareto admits a representation of Gamma mixture, allowing easy reparameterization.
Proposition 3.1. Let C. (O, In) and W Gamma(v, 1) be independent, where C, (0, In)) has i.i.d. univariate Laplace L, (0, 1) components. Then T == (v/W)Z satisfies T ~ P. 1,, v).
3.2 DERIVATION OF THE L-LOSS
The computation of (go |p4) with ? = 20-+2m+ results in a closed-form expression. By
simplifying the constants and utilizing the alternative prior expression (cf. Kim et (2024). we obtain the loss function for the Pareto VAE under the mild condition (m + 1)y + 0:
(0, = Exg-psn
20 Eewgs(112) ||:c - po(=))1 + aDD (1
(9)
where Pau = (m, k0,kmmm v mn22n, k v+2n " - &Cn/2p-111) wtnto-1 and a
y with
2Cs
C2 C2= 0== CmmcooTT....
2m &
2
-+m+ n/2
In other words, the y-loss can be understood as a sum of the mean squared enror (MSE) and the 7-power divergence regularizer between the encoder and the alternative prior, Patt-

---

Published as a conference paper at ICLR 2026
Moreover, incorporate y-flatness into the optimization process, apply the result of Theorem 2.1 to the regularizer:
(0, ④) S Enr-pi
E 2x
x~ | mo(=)l2+ +a. [1)( (0, +),
(10)
where 3 is the value in Theorem 2.1 and $o.0 90.0 Pm(2 0, (2))), + n/2). Our practical objective C. (0, Φ) thus consists of three parts: the fg-reconstruction loss, the 7-power divergence under y-flatness, and an & penalty term on (z(x))
Limiting behavior theoretically analyze the limiting behavior a, and 3 in (9) as v + oo.
Proposition 3.2 (Limiting behavior the y-loss). For For fixed n, m and x, and y = 22/4
1
lim P1000
lim k lim 3 = lim Dy Pati (x, Patr.so):
2'
20000
where 46.00 = Lm (Ma (x), (z) and Palnoc Lm (0, 1m) are limiting distributions of and Patt
Proposition 3.2 shows that as v - oc, the y-loss converges to the LaplaceVAE. (LVAE) objective with the regularizer weight 4; LVAE is VAE with Laplace prior and encoder distributions together with Gaussian decoder (Geadah et al. 2024). In this sense, Pareto VAE can be regarded as a heavy-tailed extension of the LVAE.
Morcover, in practice, the weight a and corresponding 3 (coupled with ai) in (10) can be fine-tuned in the same spirit as the B-VAE (Higgins et al. (017). Under such tuning, the modified y-loss smoothly converges to the that of B-LVAE in the V - oc limit.
3.3 DECODER SELECTION: SYMPARETO DECODER WITH & RECONSTRUCTION ERROR
We may also modify the joint decoder distribution by replacing the &g reconstruction error term || (2)11 in (8) with an & --noo-obbbee variant:
[1+
1 " (1==l1 + ||
| -
(++m+n)
po(x, 2)
(11)
This leads to the 7-loss function in which the MSE is replaced by the mean absolute error (MAE) when 7 = utm+n"
(, 6) Egapan
1 Ezngu(1))) ||2 (z(2)) + 2a-D
Pat) + 2a,3 [|
(12)
where 40,0 = Pm (0, 0 (00,0o((),,00 n), Palt = (0, klm, + n),, and
17
4: v+n (Cn,u-1,u0 -
a
Tv 202
C₂ =
(Cu+mee D D + m
0
Thanks to the MAE in the loss function, the model promotes robustness to extreme values. In the same way, we can apply symPareto decoder VAVAE, thereby obtaining MAE-based ($VAE variant; see Appendix B.5.2
EXPERIMENTS
In this section, we evaluate ParetoVAE on both low and high-dimensional datasets under different decoder choices. We provide comprehensive comparison against VAEs with different distributions: Gaussian VAE (VAE), Laplace (LVAE), and VVE. Experiments details are given in Appendix C We additionally report wall-clock comparison of l-loss training versus ELBO-basce training: Appendix C.4 Throughout, the hyperparameter y is held fixed for t VAE and Pareto VAE (without tuning); Appendix D for related discussion.
7

---

Published as a conference paper at ICLR 2026
70-
101
ParstoVA
10
LVAE
PVAAE
101
SWD
AE
40
100
Caunr flog
100
201
109-
0
20
30
(30
70
Epoch
(b) SWD convergence curve on tail (log-scale)
Model
SWD (Civeral)) (U)
(Tai) ())
p-valuc (Tail)
PAretoVAE
0.044 1: 0.005
0.170 1: 0.029
0.221 1 0.144
LMAE
0.355 + 0.000
0.301 0.084
0.119 0.000
40
69
50
PVAE
0.018
0.402 +00055
0:881 0.000
0.042
0.110
AE
0.000
(a) 2d hexbin of joint in/out-degree counts
(c) Overall SWD, SWD, tail MMD p-values
(mean + standard deviation over 10 runs). Overall
p-values are omitted since all models fail to reject.
Figure 3: Consolidated results for graph-degree reconstruction on the Epinions network.
4.1 SPARSE HEAVY-TAILED DATA ANALYSIS WITH D-DEODEER
We evaluated Pareto VAE with a f-decoder on two representative tasks: low-dimensional graph degree distributions and sigh-dimensional word token vectors. As a point of comparison for tail modeling, we also include a deterministic autoencoder (AE) baseline, We employed the two-sample MMD test (Gretton al. 2012) test the hypothesis
Ho Pattta
Precom
vs.
H1
Piata
Precon
(13)
and report the resulting p-values, where Precon is the distribution of the reconstructed samples.
constructed a joint in- and out-degree count distribution of directed graph using the SNAP Epinions social network (Richardson et 2003). In addition to the MMD test, we also measured the sliced I-Wasserstein distance (SWD; Bonneel et al. 2014) in both head and tail regions, defined the top 10% of samples according to their &g norm.
Graph degree reconstruction
Figure 3a shows a hexagonal binned scatterplot (hexbin) of the joint degree counts, revealing a sparse, power-law shape. Figure presents SWD convergence curves during training. We observe that the models with l f--regularization (ParetoVAE, LVAE) converge more rapidly than those with (E-regularization (VAE, VAE), indicating enhanced robustness extreme values. Figure 3c summarizes SWD metrics and tail MMD p-values. In all metrics, ParetoVAE achieved the lowest SWDs in both parts, showing its ability capture tails and sparse extremes. By contrast, VAE and AE often rejected Ho, indicating poor tail fits.
Word frequency analysis Based on the WikiText-2 dataset (Merity et al. 2017). we constructed 19, 99,622diimensoonn bag-of-words representations from all tokens occurring least 5 times. For head-tail analysis, the head was defined as the set of the 2,241 most frequent words, and the tail as the set of the 2,241 least frequent words among those with frequency 2 5. For each part, we report the overlap ratio between the sets induced by reconstructions and the ground truth and the Jaccard similarity between the two sets.
As summarized in Table PAretoVAE achieved the highest overlap and Jaccard scores for both head and tail. It also attains non-nee -eglliiiee MMD p-values on both parts, whereas baselines typically rejected Ho in at least one part, indicating that Pareto VAE better captures the plwer-law structure. In contrast, AE captured neither tail nor head, and the others were only good reconstructing head.
8

---

Published as a conference paper at ICLR 2026
Table 2: Metrics for word frequency reconstruction on Wikitext-2. Reported values are the mean + standard deviation over three runs. v denotes cases where the p-value does not reject the null hypothesis at conventional significance levels (0.05 or 0.01), whereas indicates rejection.
Head
Tail
Model
Overlap (1)
Jaccard (1)
p-value
Overlap (%)
Jaccard (1)
p-value
VarctoVAE
0.981 + 0.001
0.964 + 0.001
0.417 + 0.237 /
0.717 0.035
0.560 0.043
0.178 0.161 /
LVAE
0.772 0.008
0.629 0.010
0.233 + 0.148 /
0.230 0.003
0.130 0.002
0.001
0.000
VVE
0.739 0.002
0.586 + 0.002
0.665 + 0.116 /
0.226 = 0.001
0.127
0.001
0.001
0.000
VAE
0.775 +: 0.017
0.633 = 0.022
0.229 + 0.200 /
0.224 = 0.009
0.126
0.006
0.001
= 0.001 X
AE
0.642 0.007
0.473 + 0.008 0.001 + 0.000 X
0.197 = 0.004
0.109 0.003
0.001 0.000 X
Table 3: Quantitative results for denoising task with a noise probability of 0.5 on various datasets.
Mean + standard deviation reported for 5 runs (class-relateds metrics are not reported for CelebA).
Dataset
Model
PSNR (7)
SSIM (1)
Accuracy(1 (1)
Consistency (1)
VAE
24.185 1 0.074 0.950 +: 0.001
0.909 +: 0.005 0.983 + 0.001
LVAE
20.246 =: 0.154
0.891 0.025
0.828 + 0.011
0.954 + 0,005
MNIST
VVAE
22.985 = 0.389
0.935 1: 0.005
0.799 0.019
0.979 + 0.001
VAE
18.516 = 0.138
0.840 0.741 0.002
0.004
0.913 0.002
PAE
25.980 + 0.630
0.878 0.010
0.215 +: 0.005
0.782 + 0.021
LVAE
23.276 = 0.374
0.786
# 0.018
1:204:00000
0.528 + 0.057
SVHN
VAE
23.488 = -.864
0.797
+: 0.076
0.198 1: 0.003
0.570 + 0.208
VAE
21.028 = 0.309
0.690
1 0.013
0.197
0.276 + 0.020
VAE
20.552 1 0.119
0.023
0.368
+.000200002
0.341 + 0.006
LVAE
16.622 = 0.105
0.498 + 0.003
0.285
0.181 + 0.005
CIFARIO
VAB
19.584 = 0.064
0.664 0.004
0.001
0.259 + 0.014
VAE
15.726 = 0.068
0.457
1 0.006
0.016
0.165 1 0.005
VAE
20.783 1 0.389
0.903 1 0.007
0.257 1 0.003
0.564 0.038
LVAE
11.953 = 0.031
0.716 + 0.001
0.222
0.026
0.000 1 0.000
Omniglot
VAB
11.918 + 0.034
0.712
+ 0.001
0.032
0.005
0.000 0.000
VAE
11.917 = 0.035
0.712 + 0.001
0.033
0.002
0.000 0.000
Pareto VAE
25.125 1 0.065
0.818 1 0.002
LVAE
21.487 = 0.027
0.708 1 0.003
CelebA
VVR
22.406 + 0.452
0.741 + 0.012
VAE
18.554 = 0.017
0.598 1 0.000
4.2 IMAGE DENOISING APPLICATION WITH SYMPARETO DECODER
We conducted a denoising task using ParetoVAE with a symmsycee Pareto decoder on benchmark datasets encompassing MNIST (Deng 2012). SVHN (Netzer et al. 2011). Omniglot (Lake et al. 2015), and CelebA (Liu et al. 2015). Evaluation metrics were comprised of PSNR and SSIM to assess the quality of images reconstructed from corruption. For class-related metrics, we report: (1) accuracy, obtained by applying linear probing directly to the latent variables sampled from the corrupted images. (2) (class) consistency, quantified as the agreement between the class predictions of a clean and reconstructed image using an pretrained external classifier exclusively on clean images.
Similarly Vincent et al. (2008), we corrupted the original data obtain noisy inputs x'. which are passed through the VAE models to reconstruct 4. The main objective includes an MAE loss between x and é. This procedure has been shown to yield theoretically valid objective for VAEs (Im et al. (016). and we extend this argument to Pareto VAE in Appendix C.5 The results are summarized in Table Across all datasets, Pareto VAE consistently outperformed baseline models. We attribute this ad advantage to the MAE term in the ylooss, which improves robustness to outliers and enhances recovery of fine-graing ireeggraied structures. The V'VAE was competitive on MNIST but its performance degraded complex datasets like Omniglot, highlighting the advantage of Pareto in high-dimensional denoising.
Figure commpare denoised images from datasets corrupted with salt-and- -alt-and- noise, applied either in RGB or grayscale. On CelebA (Figure 4a), e VVAEE yielded visually sharp reconstructions than VAE, but often generated samples not too close to the input. LVAE better preserved input-specific features,

---

Published as a conference paper at ICLR 2026
Original
5
Noisy
J S
3
"
ParetoVAE
(b) Original
(c) Noisy
LVAE
~VAE
6
3
VAE
3
4
(a) CelebA
(d) ParetoVAE
(c) The others
Figure 4: Denoising results on the benchmark datasets (with noise probability of 0.5). (a) displays the results on CelebA. For Omniglot, (b) shows the original images, (c) the noisy images, and the denoised images from ParetoVAL and other VAEs are shown in (d), and (e), respectively.
though with lower fidelity. Pareto appears to combine both strengths, retaining distinctive features while achieving superior reconstruction quality, illustrating the benefit of the symPareto distribution.
Interestingly, Omniglot, all models except Pareto VAE failed to reconstruct the original images from noisy inputs, even after fine-tuning hyperparameters Frigure4d]lnn and Figure This outcome may stem from the fact that Omniglot is highly sparse in terms of classes, making difficult for light-tailed or t- b-based regularization schemes to capture its underlying structure.
5 CONCLUSION
We have introduced the multivariate symmetric Pareto distribution (symPareto) and proposed Pareto- VAE, which employs symPareto prior and encoder with flexible decoders (Student's t symPareto). By leveraging the information-gcometii 7-power divergence and the joint minimization viewpoint to variational inference, ParetoVAL VAE induces tractable y-loss as alternative the ELBO-based one. Empirical results show that Pareto captures the power-law behavior of the data and enhances robustness to outliers. Further discussion is provided in Appendix D
ACKNOWLEDGEMENTS
All authors were supported by the AI-Bio Research Grant through Seoul National University (No. 0413-20222 20230000), Artificial Intelligence Graduate School Program (Seoul National University), and by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT, No. RRS-2024-00336611,
REFERENCES
Najmeh Abiri and Mattias Ohlsson. Variational auto-encoders with student's t-prior, 2020.
Haleh Akrami, Anand A. Joshi, Jian Li, Sergül Aydöre, and Richard M. Leahy. A robust variational autoencoder using beta divergence. Knowledge-Base Based Systems, 238:107886, 2022. ISSN 0950-
7051.
Shun-ichi Amari. Information geometry and its applications, volume 194. Springer, 2016.
Barry C. Arood, Enrique Castillo, and Jose María Sarabia. Multivariate distributions with generalized pareto conditionals. Statistics and Probability Letters, 1765)3611336, 1993. ISSN 0167-7152.
August A. Balkema and Laurens de Haan. Residual life time at great age. Annals of Probability, (5)-792-804, 1974.
10

---

Published as a conference paper at ICLR 2026
Albert-László Barabási and Réka Albert. Emergence of scaling in random networks. Science, 286 (5439):509-51) oct 1999. ISSN 1095-9203.
Nicolas Bonneel, Julien Rabin, Gabriel Peyré, and Hanspeter Pfister. Sliced and radon wasserstein barycenters of measures. Journal of Mathematical Imaging and Vision, 51, 04 2014.
Aymene Mohammed Bouayed, Samuel Deslauriers-Gautinn Adrian laccovelli, and David Naceache. Conditional-t Equitable latent space allocation for fair generation, 2025.
Ziyu Chen, Hyemin Gu, Markos A. Katsoulakis, Luc Rey-Bellet, and Wei Zhu. Robust gener- ative learning with lipschitz-regularized a-divergences allows minimal assumptions on target distributions, 2024.
Imre Csiszár. Information geometry and alternating minimization procedures. Statistics and Decisions, Dedewicz, 1:205-237, 1984.
Jacob Deasy, Tom Andrew McIver, Nikola Simidjievski, and Pietro Lio. a-VAEs Optimising variational inference by learning data-dependent divergence skew. In ICML Workshop on Invertible Neural Networks, Normalizing Flows, and Explicit Likelihood Models, 2021.
Li Deng. The mnist database of handwritten digit images for machine learning research. IEEE Signal Processing Magazine, 29(6): 141-142, 2012.
Shinto Eguchi. Pythagoras theorem in information geometry and applications to generalized linear models. In Handbook of Statistics, volume 45, pp. 15-42. Elsevier, 2021.
Kai-Tai Fang, Samuel Kotz, and Kai-Wai Ng. Symmetric Multivariate and Related Distributions. Classics in Applied Mathematics. Chapman and Hall/CRC, 2017.
Richard M. Feder, Philippe Berger, and George Stein. Nonlinear 3d cosmic web simulation with heavy-tailed generative adversarial networks. Phys. Rev. D, 102: 103504, Nov 2020.
Griffin Floto, Stefan Kremer, and Mihai Nica. The Tilted variational autoencoder: improving out-
of-distribution detection. In The Eleventh International Conference on Learning Representations, 2023.
Hao Fu, Chunyuan Li, Xiaodong Liu, Jianfeng Gao, Asli Celikyilmaz, and Lawrence Carin. Cyclical annealing schedule: A simple approach to mitigating KI. vanishing. In Jill Burstein, Christy Doran, and Thamar Solorio (eds.), Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics; Human Language Technologies, Volume (Long and Short Papers), pp. 240-250, Minneapolis, Minnesota, June 2019, Association for Computational Linguistics.
Victor Geadah, Gabriel Barello, Daniel Greenidge, Adam S. Charles, and Jonathan W. Pillow. Sparse-coding variational autoencoders. Neural Computation, 36(12):2571-2601, 2024.
Arthur Gretton, Karsten M Borgwardt, Malte J Rasch, Bernhard Schölkopf, and Alexander Smola. A kernel two-sample test. The Journal Machine Learning Research, 13(1):723-773, 2012.
Tian Han, Jun Zhang, and Ying Nian Wu. From em-projections to variational auto-encoder. In
NeurlPS 2020 Workshop: Deep Learning through Information Geometry. 2020.
David D. Hanagal. A multivariate pareto distribution. Communications in Statistics Theory and Methods, (377))1471-14833 1996.
Xiaoran Hao and Patrick Shafto. Coupled variational autoencoder. In Proceedings of the 40th International Conference on Machine Learning, 2023.
Tennessee Hickling and Dennis Prangle. Flexible tails for normalizing flows, 2024.
Irina Higgins, Loic Matthey, Arka Pal, Christopher Burgess, Xavier Glorot, Matthew Botvinick, Shakir Mohamed, and Alexander Lerchner. beta- VAE: learning basic visual concepts with con- strained variational framework. In The Fifth International Conference on Learning Representations, 2017.
11

---

Published as a conference paper at ICLR 2026
Todd Huster, Jeremy Cohen, Zinan Lin, Kevin Chan, Charles Kamhoua, Nandi O. Leslie, Cho- Yu Jason Chiang, and Vyas Sekar. Pareto gan: Extending the representational power of gans to heavy-tailed distributions. In Marina Meila and Tong Zhang (eds.), Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pp. 4523-4532. PMLR, 18-24 Jul 2021.
Daniel Jiwoong Im, Sungjin Ahn, Roland Memisevie, and Yoshua Bengio. Denoising criterion for variational auto-encoding framework, 2016.
Weonyoung Joo, Wonsung Lee, Sungrae Park, and Il-Chul Moon. Dirichlet variational autoencoder. Puttern Recognition, 107:107-514 2020.
Juno Kim, Jaehyuk Kwon, Mincheol Cho, Hyunjong Lee, and Joong-Ho Won. - autoencoder: Learning heavy-tailes data with student's and power divergence. In The Twelfth International Conference Learning Representutions, 2024.
Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization, 2017.
Diederik P Kingma and Max Welling, Auto-encoding variational bayes. In The Second International Conference on Learning Representations, 2013.
Nicolas Lafon, Philippe Naveau, and Ronan Fablet. A vae approach to sample multivariate extremes, 2023.
Brenden Lake, Ruslan Salakhutdinov, and Joshua Tenenbaum. Human-level concept learning through probabilistic program induction. Science, 100:1332-1335 12 2015.
Mike Laszkiewicz, Johannes Lederer, and Asja Fischer. Marginal tail-adaptive normalizing flows. In Proceedings of the 39th International Conference on Machine Learning, volume 162. PMLR. 2022.
Yingzhen Li and Richard E Turner. Rényi divergence variational inference. Advances in Neural Information Processing Systems, 29, 2016.
Qi Lian, Yu Qi, and Yueming Wang. Cauchy difffusio:: A heavy-tailed denoising diffusion proba- bilistic model for speech synthesis. Proceedings of the AAAI Conference Artificial Intelligence, 39(23):24549-244 24557, Apr. 2025.
Ziwei Liu, Ping Luo, Xiaogang Wang, and Xiaoou Tang. Deep learning face attributes in the wild. In Proceedings of International Conference on Computer Vision (ICCV), December 2015.
Ilya Loshchilov and Frank Hutter. Decoupled weight decay regularization, 2019.
K. V. Mardia. Multivariate Pareto Distributions. The Annals of Mathematical Statistics, 33(3):1008-
1015, 1962.
Stephen Merity, Caiming Xiong, James Bradbury, and Richard Socher. Pointer sentinel mixture models. In International Conference on Learning Representations, 2017.
Yuval Netzer, Tao Wang, Adam Coates, Alessandro Bissacco, Bo Wu, and Andrew Y. Ng. Reading digits in natural images with unsupervised feature learning. In NIPS Workshop on Deep Learning and Unsupervised Feature Learning 2011, 2011.
MEJ Newman. Power laws, pareto distributions and zipf's law. Contemporary Physics, 46(5): 323-351, September 2005. ISSN 1366-5812.
Kushagra Pandey, Jaideep Pathak, Yilun Xu, Stephan Mandt, Michael Pritchard, Arash Vahdat, and Morteza Mardani. Heavy-tailed diffusion models, 2024.
James III Pickands. Statistical Inference Using Extreme Order Statistics. The Annals of Statistics, 3 (1):119 131, 1975.
Matthew Richardson, Rakesh Agrawal, and Pedro Domingos. Trust management for the semantic web, 2003.
12

---

Published as a conference paper at ICLR 2026
Holger Rootzén, Johan Segers, and Jennifer L. Wadsworth. Multivariate generalized pareto distribu- tions: parametrizations, representations, and properties, 2017.
Lukas Sablica and Kurt Hornik. Hyperspherical variational autoencoders using efficient spherical cauchy distribution, 2025.
Henry Scheffe. Useful Convergence Theorem for Probability Distributions. The Annals of Mathematical Statistics, 18(3):434 438, 1947.
Hiroshi Takahashi. Tomoharu Iwata, Yuki Yamanaka, Masanori Yamada, and Satoshi Yagi. Student-t variational autoencoder for robust density estimation. In Proceedings of the 27th International Joint Conference Artificial Intelligence, pp. 2696-2702, 2018.
Francesco Tonolini, Bjöm Sand Jensen, and Roderick Murray-Smith. Variational sparse coding. In Ryan P. Adams and Vibhav Gogate (eds.). Proceedings of The 35th Uncertainty in Artificial Intelligence Conference, volume 115 of Proceedings of Machine Learning Research, pp. 690-700. PMLR, 22-25 Jul 2020.
Pascal Vincent, Hugo Larochelle, Y. Bengio, and Pierre-Antoine Manzagol. Extracting and composing robust features with denoising autoencoders. In Proceedings of the 25th International Conference on Machine Learning, pp. 1096-1103, 01 2008
Yifan Zhang, Bingyi Kang, Bryan Hooi, Shuicheng Yan, and Jiashi Feng. Deep long-tailed learning: A survey. IEEE Transactions on Puttern Analysis and Machine Intelligence, 45(9):10755-18816 2023.
Shengjia Zhao, Jiaming Song, and Stefano Ermon. Infovae: Balancing learning and inference in variational autoencoders. Proceedings of the AAAI Conference on Artificial Intelligence, 33(01): 5885-5892, Jul. 2019.
13

---

Contents
Introduction
2
Theoretical Background
2.1
Variational Autoencoder (VAE)
2.2
Distributions Power-law Families
2.3
Information Geometry and 7-power Divergence
3
The Pareto Variational Autoencoder
The PAretoVAE, Structure
3.2
Derivation the y-loss
3.3
Decoder Selection: SymPareto Decoder with Reconstruction Error
4
Experiments
7
4.1
Sparse Heavy-tailed Data Analysis with E-decoder
8
4.2
Image Denoising Application with SymPareto Decoder
Conclusion
10
Some Properties of the Multivariate Symmetric Pareto Distributions
15
A.1
From Mardia's One-sided Pareto to the Symmetric Pareto via Reparameterizztion
15
A.2 Statistics of the SymPareto Distribution
16
A.3 Convergence in Distribution of the SymPareto
17
Proof of Theorem
18
Computation for Pareto VAE
19
B.
Proof of Proposition 3.1
19
General Decoder Construction: Student's s and SymPareto
20
B.3
Derivation the y-loss for ParetoVAE with t-decoder
21
B.3.1
Derivation of the Raw y-loss
21
B.3.2
Alternative Representation the 7-loss
23
B.4
Proof of Proposition] 3.2
25
B.5
Derivation of the Heavy-tailed VAES with Sym Pareto Decoder
26
B.5.I
ParetoVAE with SymPareto Decoder
26
B.5.2
VVE with SymPareto Decoder
28
Experimental Details
31
Details on Graph Degree Reconstruction
31
C.2
Details on Word Frequency Analysis
31
Details on Image Denoising Application
31
Comparison Computational Cost versus
32
Theoretical Justification of Denoising with Pareto VAE
32
Discussion
34
14

---

Published as a conference paper at ICLR 2026
A SOME PROPERTIES OF THE MULTIVARIATE SYMMETRIC PARETO DISTRIBUTIONS
In this section, provide some properties of the multivariate symmetric Pareto distributions.
A.1 FROM MARDIA'S ONE-SIDED PARETO THE SYMMETRIC PARETO VIA REPARAMETERIZATION
We start from Mardia's n-variate Pareto distribution of the first kind (Mardia 1962 whose joint density is given by
-(++)
Pn (21, ................................ In 0, T(v + 1 T(u)
Σ 2
> 0, リ>
where ē = IIm og, This form has a lower-bounded support anchored at the vector a.
To shift the support to the positive orthant with unit offset, reparameterize y; I; o.. The density becomes
-(++))
T(u + n) r(v)
Pr (11 Un l0,)
+
> 0.
Next, to align this form with the y-flat family in a; + voz. This yields the one-sided (positive-orthant) Pareto density
rescale the coordinate-win scales by
(v+n)
p (
a,u)
I(v+n) n) 1 V²T(v) σ
1
1
yk
> 0.
Finally, obtain symmetric distribution the whole space by extending the support from the positive orthant to IR.. Concretely, replace yi by [2; p4l (introducing location parameter (1) and adjust the normalizing constant by factor 2- since the support is extended all 2 orthants. This gives the multivariate symmetric Pareto (symPareto) density
Pa(x U, CR.P.
(1
1
1 (u+R)
I(v1 + n2)
which coincides with 33. Here Cn.
(2v²) (2 [')] and the division (x - p)/o is understood
elsment-wise. In this way, the construction proceeds seamlessly from Mardia's one-sided form, through reparameterization and vo-scaling that reveals the f-flat structure, to the desired symmetric extension on I...
Marginal distribution of the symPareto distribution For the n-variate standard symPareto distri- bution, the marginal distribution of each component X; is itself univariate symPareto distribution. Its density can be computed
pı(21) = 1 /
Cnew (1 + (11111))
(utn)
dxy. drn
I(v +n)
n- (v .)... (u 1++1) (1 b (u+1) -(u+1)
21(2)
(1 + (nnl
= 2
Thus each marginal has the same symmetric Pareto form in one dimension. By restricting to the positive half-line > and writing v 1/5, this reduces to
(xx) (1 (1/s-1
0,
15

---

Published as a conference paper at ICLR 2026
which coincides with the generalized Pareto distribution (Pickands 1975). This fact is directly connected to classical extreme-value theory. According to the Pickands Balkema- De Haan theorem (Balkema & de Haan 1974; Pickands [975), the GPD arises the universal limit law for threshold exceedances of a broad class of distributions. Therefore, the symPareto distribution inherits marginal tail behavior that is theoretically consistent with the limit structure predicted by extreme-value theory.
A.2 STATISTICS OF THE SYMPARETO DISTRIBUTION
Mean and variance Let X ~ (0 1, 1,) The expectation exists for v/ > 1, and by symmetry of the density we have F[ = 0. When > 2, the variance is finite and can be computed Var(X) = E[X?]
(v++1)
x
dx
2
1
-(++1)
00 x²
-(++1)
1 1
-I
dx
+
1 =
dr
2
V
0
2
3
00
(1
y
(u+1) dy +
(z 1)2² (u+1) ((++1) dz)
1 (v=1- -x,
1 1+ x
2
1
2
+ 1)
2v²
3
2
v
(v- 1)
For the multivariate case X ~ (., a. v) the covariance between distinct components also vanishes by spherical symmetry (Fang et al., 2017). although this does not imply independence.
The first centrallnnon-cenrrll absolute moment absolute moment of X is
Let ~ P1(0, 1). For v > 1, the first central
(1++))
||X = /
1
x 1 + -x
dx
J
v
111 + 1)
(v+1)
dt ー 1
The first non-central absolute moment about t E R is given by
V/
1
E| +
1
(14)
To compute this, split the integral into two cases depending the sign of t:
1w++1)
I + 1
E||X + t||
1+ |2l
dx
2
(x + 1)
1 1 101)
- (w+1)
x +
1 -
-(++1)
1 +
dr
1
|s||
dr
1.
12
1)t20 0
1 ((u+1)
1
(1-11²
1 =
(x + 1)
x w dx
2
0
1
(u(1 - y) + +tly (v++) vdy
2
v+t
1 t 1)
u²
(1 1)
-ut1
+ 2(v 1)
2
/
(1+ t
-0+1
2(p - 1)
16

---

Published as a conference paper at ICLR 2026
1
1
0
(x+1) 1) (1 1 - xt
1 -(++1)
dr +
0
(x (x+
- =D
-
-(u+1)
de
2
v
2 2
2 2/1
1
1+//v
(v(1 -y) + 1)y1
(1++1) vdy
2 2/1
1
(v(z - 1) 1)=- (v-1) vdz
++t
1 t
D²
1 t 11+1
v²
t - v
1
1 -
+
+
2
2(
2(v
1)
2
v
1 1
-v+1
=t+ 2(V- -yy 1)
1/
t -141
:. El + t|| 11 12
1+
.t20 0.
w-1 1
V/
2)
1 2 _000 (o (o+1) (1- 1 -(++1)
1 o
(1-
1 (x+ 1) (1+
- -(++1)
dr
dx
2 J
1 ( (v(1 11) + 1)y (w++)) vdy
1-1//v
(u(2 1) + 1) 1)==((ut)) (1+1) vdz
-2 2
2
(v+
v²
レ-v v
-1)
D²
1 11- -2+1
2 + 2(1 -1)
1
-
1 -
2
2(1 1)
-0+1
-l+ + 2(v - 1)
v
1 11_
t
v
-(++1)
1
(x + 1)
11 :
-x
dx
(u(2
1) + t)z (v+1) -(++1 vdz
1-t/2
v2
(1-4)
-++1
1
+ 2(v - 1)
2
1/
(1
(1
t
-++1
2(v -1] 1)
v
-2+1
V 1 (1-1
:. F|I + 1|| = 1 + T2 = -1
0.
1-
Combining both cases leads the unified expression in
A.3 CONVERGENCE IN DISTRIBUTION OF THE SYMPARETO
We now investigate convergence in distribution of the symPareto with respect to the parameter v. Notice that according to Scheffé's s theorem (Schetfe 1947), it's enough to check the convergence of the probability density functions as assuming regular condition for the measure space.
Theorem A.1 (Convergence in distribution of the symPareto under v goes infinity). For n 2 1, fix 1 RR, and E IR. Then, the limiting distribution Pn(x u, v) is given by:
Palx U, S. Ln(x (1,)),
(15)
where Ln(x I (),0) has i.dd univariate Laplace Li 03) components, which is then given by
r -ー 4 1)
Lx(x μ. 0) = 2n exp -
17

---

Published as a conference paper at ICLR 2026
Proof. We analyze the limit of the normalizing constant Cn,Y, - oc. Using the Stirling's approximation,
T(x + 1) ~ V2RO
I
x x> 0,
obtain
T(v 11)
lim Capu 40000
lim
2 1-000 vmI(v))
v-1
1 +
v-1
lim
2 1-+00
V
2
(16)
2n
For the unnormalized part (kernel) of the density, we observe that
((++n)
14
1
T.
lim 210000
1+
exp
3/ 14
)
exp
which shows that the symPareto converges pointwisely to the product of univariate Laplace.
A.4 PROOF OF THEOREM
Proof. Consider two n-variate symPareto distributions p (n, 11, and q = Pu 02 v). To analyze these distributions within 7-flat manifold, we set y utn' integration for 7-power entropy and y-power cross-entropy terms as follows:
We then compute the
/ p(a)) t2 = C o. / (
01 --
(v-1+n)
dx
n.v.v
= Cn.--ivv
+
= C n2u.
μ-1 1
/ ((x))+1da dx or C
V 12 1 v-1
since
[(u+n)
Onen
[(0)
v+n
(17)
Cn.--i
T(++n-1)
w-1
1)2v1)(2v)m
Also,
/ p(r)g(x)ddr = C2 " EEnp
x 12
1
O2
= -C
1
E Baup
Fj
12.
0z
i=l
1).
-G,
1
01.
E
10.11...))
2
(MA. = pi. (2,x))
|
V
0
|لممد|
-tt1
02
1 1
1
Equation 14
01,3
i=1
=C 7
1+
UA
1
+
1 |لساا
02
02
18

---

Published as a conference paper at ICLR 2026
Hence, the 7-power entropy and cross-entropy
H,,() - - -f ) itr
1
T Chiyje
1/ #+ 1
-
u v-1
(p = - ( ((o))(a))) (ʃ a(x))tt)
- 1,
(+ 1
1 02
us
91 1+
1
d2
Combining the two terms, the 7-power divergence between p and q yields closed-form expression: - (p)
(y(p11 = (D.
C² 02
1+ 1/
01
-++1
5j 177
V 1
02
1/
02
1/ σ1
v
(18)
where C' = (n(,, (v + n) n)C n, and 17 is used in the last equation. We then put Ma 0 and get the closed form of D,, (po go)
Next, let Po - (0, oi, and go = Pn 02 v). To obtain an upper bound Dz (p) || q) in terms
of DD, ll 9u), we apply the triangle inequality to [x; -22,
I [2, M1,il
+ |M1, M2.a As
result, the cross- entsopy term admits the following bound:
f p(z)g(x)) = Ch 02poo2
12
12,2
1
E Ezop
1/
2
im
12
01,1
34 - 11,
11,1 - 122,1
I Ch 02
1
Expup
+
1/ i=1
22
01
22
1
Ma,s - 112,5
Ch 02
1
+ E Exup
2
V 1
= Ch
1
01
11 12z
Using this result, we have an upper bound:
Dy(ell (0l1))) C
02 rテ
1+
1
σ1
+
11
12
|))
v + n
)
v-
02
1
02
1ー
SDy(Po (1) +
Cn.V-IN
(1
)
n
141 12
02
02
and equality holds when #1 μtz
B COMPUTATION FOR PARETOVAE
This section presents details of computation involved with the Pareto models. Appendix B. provides the proof of Proposition]] 3.1 The subsequent subsections describe the construction process for heavy-tailed VAE variants, including Pareto with t-decoder, Pareto with symPareto decoder, and VAE with symPareto decoder.
B. PROOF OF PROPOSITION
Proof. Let Z R and W RR be independent random variables such that:
7 (7,, 7n), 7: Lid. (0, W Gamma(y, 1),
19

---

Published as a conference paper at ICLR 2026
and denote their probability density function as fz) and fw (w)):
1
fw(w)
[(v) w e-
fz(2)
2 exp(-11=11)) 2 RR..
Let S W. Then the probability density function of S given by:
fs(s)
[s. w)dw
Izw w) wn dw
fz fww)) w" dw ZW)
f. w" ( 2 exp (xoxp( (xxp(11ss(ll)))
( ((v)) w--1 exp(- du
wuta-l exp(-w([|s|) + 1))
2"I(v)
[(u+n) n) (*T(v)
(1|s|1 + 1111++1--00 1)(t)))
Therefore, the probability density function of T vS is
l'l2 n) fx(t)
(20)aI(p) (1+ ||d||, -(++n)
P (t 0, 2).
Sampling algorithm From Proposition 3.1
we summarize the sampling algorithm for symPareto distribution from Laplace and Gamma distribution Algorithm
Algorithm Sampling algorithm for the symmetric multivariate Pareto distribution
Require: N: # of samples, m Latent dimension, v Hyperparameter for the Pareto prior Ensure: m-variate symPareto sample E R
Sample w; lid. Laplace(O, 1), for ,, m. Assign (#1,..., Tn).
Sample w Gamma(y, 1). Assign (v/w).
B.2 GENERAL DECODER CONSTRUCTION: STUDENT'S S / AND SYMPARETO
Following Table we provide explicit construction of the joint decoder distributions followed by heavy-tailed VAE variants.
po (x) x 1+ +
1 |==1|p (1=IB= + σa || - -mo(=)11
]
P. q E (1,2)
The prior can be generalized by
p(z) x 1
p E (1,2).
20

---

Published as a conference paper at ICLR 2026
which implies the symPareto prior (p 1) or t-prior 2). From this prior, the conditional decoder distribution can be derived as follows:
+
+
(1 + tt || (=(2))l
po(x 2)
[ + ||=11p)
+
g/p(v + m) 1 = 1=llz 09 || - ue (2)1l²
x +
q/p(v + m)
By simplifying the ratio of the joint and the prior kernels, we obtain the following four cases depending the choice of p0 2):
-+,|ا|| w+m
Pn x [=()))
oln; v + m
(p, 4) (1, 1)
ln
na(2),
++|=l1
I²I,, 2v + 2m
(p, q) (1,2)
po(x) 2)
20+2m
2(++11 01 u+m
Pn
x (2),
(p, (2, 1)
v+m
in
x (=(2))
=+12
o Ins v + m
(p. q) (2
*+m
which implies symPareto decoder (q 1) and (-decoder (q 2). Note that the symPareto decoders employ element-wise scale parameters (a11,,, whereas t-decoders use covariance scaling (a I,).
B.3
DERIVATION OF THE 7-LOSS FOR PARETOVAE WITH IDECCDEE
B.3.1
DERIVATION OF THE RAW 7-LOSS
The VAE starts from the heavy-tailed joint decoder model:
(u+m+
+ || - (1))
Po(x, 2) x 1 +
|I=||,
o E R+ V> 1.
o²
We first compute the normalizing constant C by separating the integration over x and 2:
{ gm
(o( 2) dædz
Jg"
1+ 1 |z||
-1++++
1
2(v + m)) 1
2n 2mmmm
2(v + m) + ||=111 02 1| uo()))1
drdz
mm
R"
2 2m
-3
{
1 |1=111
-(u+m)
dz
Rm
2m
= =T R
+
C-1
Thus, the normalizing constant C' is
2m
C =
m 2
Remark that Tn, 2mm is the normalizing constant of t,, (0, I, 2m). This decomposition naturally yields the prior and conditional decoder distributions:
p(z)
Pm (2 0, 1m, v)
1/ + ||=|||
(x 2) t, / (2(2),
2 2n, 2v + 2m)
2v + 2m
21

---

Published as a conference paper at ICLR 2026
We also define the encoder symPareto distribution with + n/2 to align the tail order with the joint exponent ( m + n/2):
q (=|=)) (m un/2))
Next, we compute the 7-power divergence between two joint distributions go and pv- Setting
7 2v+2m+n' the 7-entropy and y-cross entropy corresponding to the joint distributions pa and 4s can be expressed follows:
l| pa dodd
dxdz = Ezaup Bxmpo(-1=) [20 2>>]
1
EsNpe Exmpu(tla)
1
||=|||
0² ||= (=(=)
1
2v + 2m
V + 112211 12
CTEENp
1t
|0||
Tr
1
va²
22 + 2m -
2v + 2m
|+ 1 ||=11] )
CTEEup
|1
2v
2m
v
n/2
)
C
1+
1
1/ + m
V
V m n/2
C
1>4
v
and
(o 2))o(1 ddddz Erewpan Bsnga(laa) [2o =)
[ 1
1 ||=|| + न || 1e1==-1a 1|===m((=)13
1
Ce' Benga(tla) (110)
v n +
1
Ex
1 + [| + u(v + n/2 12/2-11 - [r + noz E Ba (l1) ||= no(2)113
where
ro(x)= ax(x) (o(a) (1
|mp(2))
(u4nn2211
(v n/ (2)
Thus, the y-power entropy becomes
(9s) = - (ʃ Goloo 2) doddz
itr
= - ( go (0(2)0)))7 Peina (x))+
147
City
v v + n/2
itr
(foola))
Pasna (x) +
it
(19)
+サ u+ 10/22-1 -
m:
P.T.
and the 7-power cross-entropy is
Cg(go.Po) = [l] go(x,, 2)po 2)" dadd)
(Jf
-17
po pa(a, 2)=+7
(v m n/22
-
1 1-1
C
1 1 ||m(())]l
n+n u(v + n/2 - 1)
Exzppon
||a(
vo² E Es || mo()))
22

---

Published as a conference paper at ICLR 2026
To approximate the power-form integral P.I. in 19). adopt the following proposition from VAE (Kim et al. 22024, which provides a first-order approximation of integrals involving Pdata with an ((7)) error rate.
Proposition B.1. Let a be any positive continuous function, y E (-1,0) 0), and suppose the values
Hi = Hy..((aana) = Pdana (x) 1+17
( |log Pautaa Paita( (z) dx
finite for each j {0, 11, & (1,
Then for any compact set Q S supp Paltoa
(1 " "Ddaaa(x)
149 - -A, [())) 17. Pdana (r)dr
r / Pduta (x) log Pdtata (x) (x)dr 0(
Proof. See Proposition 5 in Kim et al. (2024).
Using this result, the 7-entropy of the joint variational distribution qu (z, 2) becomes:
/ o 00(2) Hl((o) -c
it Pea (x)dx + C2 H. 0(?2),
where (pdam)) denotes the differential entropy.
Therefore, the 7-power divergence Dy pa) is:
1 Dy(9 po) Cg(,.. (90)
7
1 |mx(x))| + = n/2
x Ex-pin C2
ulv n/2 - 1) ||r((x)]1 + vo? Esngo(ile) || 140 (=)
-Cāo(T) # + C(H((a + +072)
To make the result simpler, let us define Gee as the ratio of C1 to C2. By discarding additive constants and the order y2 error, and applying an appropriate rescaling to the 7-power divergence, we obtain the raw y-loss:
Craw (0, ②) = Exoppi
E Ez-au((la) ||x ( 0²
+ n//
+ ||
+
(v n/2 - ||r((>)] - vorey 17n
(20)
=:fl (regularizer term)
This loss can decomposed into an &g reconstruction term and regularization term Rg. although the latter is less interpretable.
B.3.2 ALTERNATIVE REPRESENTATION OF THE 7-LOSS
We now present alternative form of the loss,, expressed as the sum of reconstruction term and a- power divergence regularizer. Consider the y-power divergence between two symPareto distributions Pak Pm, (0, kl., + n/2) and gø Pm (no (H)(()) (a(r), + n/2) with 2m Denote
23

---

Published as a conference paper at ICLR 2026
P/ = V n/2. Under the mild condition 1 + (m + 1)y o O, the divergence can be written as:
1
(96 || Patt) = (g. Patt)
-,, (Patt)
(2)
1
(2)
(1 1 |ms(iم1
-
C
h:
I+7
1
(21)
1/n
Un
h:
un 0 (x
-xp(x)) (ar))-i)))))) 17a
Un + m 1
1/1 -
my -1 1 =C*
||
+
Un
Op(x) (
1 1 [ma(v\)
V/n
1/12 -
in (x)
-oo(z) 175 (1n)
I/n + m
- 1
mi +1 h:
1
+ const,
1/2 -1
where C" = (v/ + m) -1,u, To align Rs with the above powerrr divergence, we determine k: such that
n/2 + m 1 V + n/2
DCne = (v+
Thus, the 7-lloss be rewritten the sum of the l2 reconstruction term and the 7-power divergence between encoder and alternative prior. Finally, multiplying by a factor of 2 ensures consistency with the Gaussian MSE scaling, yielding the alternative y-loss:
C (0, x) = Ez-pum
20 E. got- ||x u=(2))l + aDD, (1s ll par))
(22)
The simplified forms of and a are given as follows, with derivations provided below. 7/V a
2
(o (o"">~"/ n/ Cn///////// 1/(44n/2-1)
k 2+ n/v
> 0.
2C2
Derivation of ls and a (Pareto VAE with decoder) First, we rewrite the ratio
Cm.utn/2uutn/2
C
(2 +
2m
T(v)
V
[( + 2)
V
m
= c-1
n/2,0,0
1/
2
Using this identity, 11+7 +1 reduces to
V
v V n -
Creg
n/2
n/2 m -
v
m.vtn/2,utn/22
V uv+n22 + m - 11
17.
(v m +
n/2 - 1 17.
n/2
C1+t
+ n/2 -
v
+ n/2 1
147
n/2
V
V + 112 -
[145]
0" C-1
n/2
1
n/2,up
14 m2
v
[ = c-
+ n/2
n/2 1,
24

---

Published as a conference paper at ICLR 2026
Since
2
my
-2m
2v + - 2
+
+ 1
++7
2v + 2m
2
1 7
2v 2m + - 2
2v 2m
we simplify the constant follows:
( u n/2
n/2
T=(m=177
Czeg
+ n/2 + m -
2
(o (--
2 + 'n
Finally, we simplify the constant by using the above result:
C"
1I+ m. +1
v + m/22)
2
C -
2 C
(V
+ n1/2 m))
1C
v
2v + 2m + nt C2
20
B.4 PROOF OF PROPOSITION 3.2
Proof. We begin by deriving the limits of the normalizing constants Cr, and Tnp as - 0o. In 16 we already show that
lim Cr,
1-000
2n
by using the Stirling's approximation. Using the same approximation, we obtain
I( M
lim TT = π
lim
2-000
+ v $ ((K)
2-2
1
1+n-2
V #+n- 11
0-2
2e
π
lim #
リ-2 2
(2x) 2
Since THE as v and the base inside the exponent remains bounded away from zero and
infinity, it follows that
2m _
)
#
セー
)
lim C2 lim 1-000
2
v + m + n/2
and thus
2v
lim a lim
(23)
-0000
2 + 2m +
2C² 2'
Next, we derive the limit of and 8, which follows directly from the expression.
2
lim k = lim
vin/2 = -=
10000
2-00
+ vl
lim B lim
(1 + n2
C.4y
1-+00
1-000
D)
1.
25

---

Published as a conference paper at ICLR 2026
Next, we compute the limit of (96 Pai)) in 23
As oo, each term converges to:
V + n/2 2 + m
}
lim C"
k
- Tx(x)
m (x) - m
1-000
+ n/2
C
lim
(1
p 00000
+
k
C"
0x(x))
(x))
(0 (v+n/22+11
|1o(6)1 [x(x)
lim
1+ (1+ + n2/2 (x(xx)
(x(x) exp
+ -n2
k
Combining these limits, we obtain
T (2) (- + [[,,((x)) + log (x)) - 1]
lim (go Pat) = [[oas(o) exp 4-000
To verify that this limiting form coincides with the KL divergence between two product form of univariate Laplace, consider ((2) = I (z(2) and p(z) =
where qi = L1
Mr,i, b1,s) and p; [1(2;] b2.))
The KL divergence between p(2) and ((2) can be factorized as
(xl(allp) = Σ DKL (9: Ps)
m
b exp
[ + |M
M2,2l
$2,
+ log
1.
61,5
By setting M1, = (x), #2, 0, and bi. d (2). the distributions p(z) and g(2) correspond to the limiting forms of do and (z)(()) in Theorem A.1
We therefore conclude that
lim Dy ll Pan) DEL (90. Pat.
B.5 DERIVATION OF THE HEAVY-TAILED VAES WITH SYMPARETO DECODER
B.5.1 PARETOVAE WITH SYMPARETO DECODER
To match the exponent in we set 7 u+m+n Since the joint decoder follows a symPareto structure, its normalizing constant is C Cu+m, on We then compute the 7-entropy and y-cross entropy for the joint distributions po and gp.
The only difference from Appendix B.3 lies in the terms involving pe (x]2). Repeating the same steps as Appendix B.3 obtain
po 2) 1+ drd= Addd
Ezapz Bampo(lle) (x)
1
1
CC"Ezopz E Bamppo(l)
1
+
|21l
σ po(=))l1)
C'E
1
1
||
+
V - m
(++112111)
1
+ m ++m
)
1+ v ||=11| )
1
C'Ezapz
11
v m
(1 [+ m
n
C
1+
1/ + m 1
v /サ n -"+
C
vー1
26

---

Published as a conference paper at ICLR 2026
/ (ʃ (0(2)))(+(dz) Pdata (2) dx
City Coptwuta C Cootn -utn / [oo(x))~>>eer Panta (2) dx,
|| 2o (a, =)po dadz Eunppe (angg(-la) (Po (2, 2)>]
C Ceeppp (- (-1e) 1 |=||| + | | |
1 [|(( + vlv + 1 [[o(())111 + _ va
+
|| no(2)ll
where
-(vtn-1)
Fø(r) : ax(i)
(u m))()()))) (1+oon))ooot)
Hence, the y-cross entropy is given by:
C7(96.P)) (ʃʃ
(J)
96 (z, 2)po( 2)7 ddrdz)
po 2) dxdz
V + n + m
= -
1
C:C²
X E Eupppn
1
[|x(())l1 + V a
[[x(1)]
+ E2- || -
1
11v + 1)
vo
and the y-power divergence (46 Po) becomes:
Dy(eo po)
Cr(9g.
H. (90)
1 Ezngs((12) || (o(2))
*+ n [|mo(t))l + v(u + - - [
Eyrpan C2
E
vo
-Iās(x)) - + H( + 0(73)
where Cj is as defined in tion B.I
19 and the error term 0(72) follows from the application of Proposi-
Discarding additive constants and up anirrelevant positive scaling factor (cf. Appendix B.3 and use Grey := C1/C2 again, raw 7-loss can written as
(-, ()) Exropa
E Eelg | pa (=)I1
(24)
v
[lla(x))l1 + (v +
(25)
= Rg (regularizer)
Compared with (20), the only difference lies in the reconstruction term, changing square of & norm & one.
We further represent alternative form, by solving the equation with respect to k: +
μ 1) k k17+1
DCreg (v + n))
Then, the alternative form of the y-loss becomes
Call 6) = E Exrgpe Ez~gs(t<> || Ho(2)]l + +aD aDx 11 pan)
where Pak = (0, klm, v n)) and k:
v
(Cn
v4n
C
27

---

Published as a conference paper at ICLR 2026
Derivation of i and (ParetoVAB with symPareto decoder) expressions for and a: Since
We first simplify Ceg to obtain
Cmptnuun
T(v + n + m) 2
Cn+muv
2n 11) + n)
[ + n m)
(T))
T(v+ n)
v
m
=Cr.uv
V +
follows that
-+
C++t -eetnuta
C1
Creg =
C2
C
-in+m- V-1
IT
V v + m 1
Cmmtavtn
(v - 1) (n 1)
Catmppo
v m -
179
(v- 1) (v+ 145
V m
This yields that
1t
V v n - + n V v n + m
V
Creg
I+(++
V
1/ ! v-1
Cr, 1T 245
1+(m+1)y
V n
V w n2
v nn n
V w-
v
v+n
Finally, a is given by
1 U+1
C C
a
V + n + m
1/ V+N +
C F.
T
v+n-1 C₁
/ n + m
m.utnnlyttn' vtn V + m - 1-1C2
y+ + m, Cootavtn
V
0-1
V m 1C2
TV
C2
B.5.2
VVE WITH SYMPARETO DECODER
From Appendix B.2
the joint decoder distribution is
po 2) x 1 , (1ea113 || -mo(t)]l)] )1
(- +n)
The corresponding prior and decoder are
p(z)= ty (z 0, Im v), Pola 2) Pa mo(z), 2(v + [|=12)
m
o1,,
- m
2
28

---

Published as a conference paper at ICLR 2026
We also the encoder distribution
(= x) = mm (= (d) (x)) M,, + 2n), xx(x) E R4
which differs slightly from Kim et al. (2024) (diagonal covariance, no (1 [1+ scaling). This
simplification is harmless, since (x)
(x) is learnable and typically assumed diagonal in practice.
We first derive the normalizing constant of the joint decoder:
= Tmmmu C R. Dam vtm
2v
Set = ++m+2n to match the exponent in the joint decoder distribution. We then derive the y-entropy and 7-cross entropy for the joint distributions pe and qs- From Kim al. (2024), the y-entropy is
2
rir
V V+ 2n -2 2
(Galit) Poama ( tx
Hasss) = - -
m,
+ 2n
C1
For the joint distribution pe, we compute
pe (xx +) dxdz = Es-pp (-1=))))=)) [po ]
CC"Eaap Bewpolt1=)
1+
|111 +
|2 [o(2))l1)
2(v + ||=112)
v + m
|| +
1/
V v + m
μナ m
+ m - 2 )
2n
1 v 11=13 +1+
1
C
2n
1
m
v + m -
v
(v 2n + m
C
V -
// 40 (0,, dxdd = Ennpue Ez (xk)) [20 =
C C' Emnapan (-1))
|| + || - (o-=(())ll))
C gesspsn
[
w(v + 2n 2) |1a(())]
+ 2n
+ vo
1
F
(o()))) || m=(2))l
Thus,
co( -
(ʃʃ
=)po(0, 2) 2)dzd))) (J po(s, 2) 2)2+
-
(v+ 2n - 2
+
-
1 2
+ 2n
1 Ex ||> no(=)]l
X Earopan 1 [| +
||
v(v + 2n 2)
Da
29

---

Published as a conference paper at ICLR 2026
and with Proposition B.1
the y-power divergence (96 po) becomes:
Dy( ll po)
= (gp. Hg))
9
x Egpopin
C2
E B2mgs((a) 1| - (o(())1 + [1lm(())12
v(v + 2n 2) [|ao(2))l1)
+ 2n
(C) rेr + (pH(pa))))
+ +0(72).
Following the same steps as before, the raw y-loss is
L (0, 6) = Eupppa Bewgot-1a) || - me(2)]l
v 2n
+ +1lo(1))122 (v+ 2n 2)
|,ا(2))l1 - (2)) 17
=: (rogulaaizz))
Compared with the previous alternative 7-loss Pareto VAE models (22 and 24
notable that the l regularization term of (z) is changed to the E term and ||r(())] (or (z))l)) is simplified
to (z)()))
We finally obtain alternative representation of the 7-loss by solving the equation with respect to k:
Dereg (v + 2n)
1/ 2n + m - 2
1
my
+2
2n 2
Accordingly, the y-loss can be rewritten as
(0, (0, 6) = Ex-pens E (x | + (Dx(40 1l pan)
where Pak fm (0, 22 21),
20 and
$
$
vT
k=
v
(20") (20")-
2 Cn.
v
V 2n
+
Derivation of k and a (t VAE with Pareto decoder)
We first compute Creg as:
T mn
(ttm-2n-2
Grog = G
2+2n-2
C2
(Tm,CC umm zim -1 umm
2
+
2-2 v42n+m-2
147
-17
1 + m + 2n 2
1/- 13-1
m,0+2n Tm, On
1/ + m
pam pott
2n 2
1/ + m + 2n 2
2n 18
Cn.2.2 1
1
(202)
F 2n 2
v m
This yields a simpler form of A as
14
+ 2n -
21(m11)
k
Grog
V + 2n
4 n 2m
1/2
$ v12n
(20-) Cn,
V +
30

---

Published as a conference paper at ICLR 2026
Finally, a is given by
+ 2n
147
C"
V 2n 2
-
V 2n 2 C
> +2++ 2n)
- 2n + m
m,u+2n
V ++ 2n v 2n + m 2 C2
C2
C EXPERIMENTAL DETAILS
Common Environment All experiments are conducted using Python 3.12 and PyTorch 2.8 on Linux Ubuntu 24.04 LTS system equipped with Intel Xeon Platinum 8568Y+ @ 2.3GHz (48 cores) processors and an NVIDIA RTX 6000 Ada GPU (48GB), with CUDA 13.0 and cuDNN 8.9.7. Additional details for each experimental setting are provided below.
C.1 DETAILS ON GRAPH DEGREE RECONSTRUCTION
The SNAP Epinions social network (Richardson et al. 2003) consists of 508,837 directed edges and 75,879 nodes. From this dataset, we compute the in-degree (the number incoming edges) and out-degree (the number of outgoing edges) for each node. This results in two-dimensional dataset with 75,879 (in-degree, out-degree) pairs, which exhibit a power-law distribution.
The encoder and decoder are implemented 2-layer MIPs with 256 units per layer. The latent dimension is set to 4. with batch size of 512. We use the LeakyReLU activation function (negative slope le-4), followed by batch normalization. The model is trained with the Adam optimizer (Kingma & Ba 2017), using learning rate of 0.0001. The dataset are randomly split into vainn/alidation//estt sets with proportions of (0.6, 0.2, 0.2).
C.2 DETAILS ON WORD FREQUENCY ANALYSIS
We utilize the WikiText-2-raw dataset (Merity et al. 2017). All non-English characters are also filtered out. Although the full vocabulary contains 65,653 tokens. we choose to exclude words with frequency of 4 less. This is because these rare words, totaling 44.791, constitute approximately 70% of the full vocabulary and serve as significant impediment to stable model training. The remaining vocabulary of 19,962 tokens retains sufficiently heavy-tailed distribution, which is characteristic of natural language.
For the analysis, the vocabulary are represented as 19.962-dimensional bag-of-words vectors. The head-and-tail split is defined as the 2,241 most frequent words (head) and the 2,241 least frequent words (tail), corresponding to words appearing exactly five times.
The encoder and decoder are symmetrically structured, each consisting of a 3-layer MLP with 256 units per layer. The latent dimension is set to 64, and the batch size was 512. The ReLU activation function is used for all layers, and batch normalization is applied without the use of dropout. All models are trained using the Adam W optimizer (Loshchilov & Hutter (219)) with a fixed learning rate of 0.0005, without learning rate scheduling. The regularization weight was linearly scheduled.
C.3 DETAILS ON IMAGE DENOISING APPLICATION
Salt-and-pepea men-ee noise was applied to the datasets, with noise probability defined the likelihood of each pixel being corrupted. For corrupted pixels, a salt ratio 0.5 was used, indicating an equal probability of pixel being set to its maximum (salt) minimum (pepper) intensity value.
The MNIST, SVHN, and CIFAR10 datasets are resized to 32 x 32 pixels, while the Omniglot and CelebA datasets are resized to 64x64 pixels. The latent dimensions are set 16 for MNIST, 32 for SVHN, 64 for both CIFAR10 and Omniglot, and 128 for CelebA. batch size of 512 is used for all datasets.
31

---

Published as a conference paper at ICLR 2026
The encoder and decoder constructed symmetrically in all experiments, with the number of stacked residual blocks being the architectural difference between datasets. Specifically, the networks for MNIST, SVHN, and CIFARIO consist of 6 residual blocks, while the larger networks for Omniglot and CelebA utilized 8 residual blocks.
The models were trained using the AdamW optimizer with fixed learning of 0.0005. No learning rate scheduler is utilized. The regularization weight is selected from three schedules conndulle-cosssatt, linear, and logistic-tuned specifically for each model and dataset. We observe that the logistic schedule is more effective for 1-channel images, while the linear schedule yields better results for 3-channel images.
The external classifiers are trained for each of the MNIST, SVHN, CIFARIO, and Omniglot. Their classification accuracies are 0.9976, 0.9611, 0.8452, and 0.8497, respectively. Specifically for Omniglot, the 1,623 character classes are re-categorized to 50 language classes.
C.4 COMPARISON OF COMPUTATIONAL. COST VERSUS ELBO
This section quantifies the computational trade-off between ELBO-based training and our proposed yy-loss for learning the Pareto We focus the setting used in experiments: a symPareto encoder/prior with a Student's t decoder.
ELBO decomposition. For a datapoint æ, the ELBO can be decomposed as
ELBO(x) [gg(=/=) (log pa [z]=)] Egg (=1=) [log [=(=a=)] + plo))
With a t-decoder, the reconstruction term is of the form
| - (=) log 1+
) (v+ mjo²
E ( Pog oc (z D
The entropy term Eg [log 2o (2|x:) admits a closed form for our symPareto encoder; in our notation,
2(v + n + m) + n2 1)2
E (= (log go(()))))))
Cn.u-n.vn
Σ log (x,( + const.
In contrast, the prior enoss-entropy term (=]=) [(=] has no closed-form expression under the symPareto prior. Consequently, ELBO training requires Monte Carlo approximation,
E [o)
1 K
Σ log ogg((
200 go(r)))
which (i) increases the wall-clock cost roughly linearly in K and (ii) introduces additional estimator variance. The latter effect is amplified when (2(2)) is heavy-tailed.
Wall-clock time comparison We measured the wall-clock time per training step (forward + backward optimizer update) on MNIST with batch size 128, and report relative runtime normalized so that the y-loss runtime equals 1.0. For reliable timing on GPU, we performed a warm-up up phase and synchronized CUDA before and after timing.
Table summarizes results for latent dimensions m E {32, 64) and Monte Carlo sizes K E (1,4, 16). At K 1, the runtime difference is small. However, the gap grows either the Monte Carlo size K or the latent dimension increases. This is consistent with the fact that ELBO training must repeatedly evaluate log ((2 for heavy-tailed samples, increasing both compute and estimator variance. Overall, these measurements support the practical advantage of the 7losss formulation when one seeks stable training without increasing Monte Carlo budget.
C.5 THEORETICAL JUSTIFICATION OF DENOISING WITH PARETOVAE
In this section, we formalize the denoising task in Section and show that its y-loss can be expressed of an e reconstruction loss and 7-power regularizer.
32

---

Published as a conference paper at ICLR 2026
Table Relative runtime per training step (forvard-bacctwaddduptttt)) compared the proposed y-loss (normalized to 1.0) on MNIST with batch size 128.
Method
Latent dim m = 32
Latent dim m = 64
Pareto VAE (y-loss)
1.00
1.00
ParetoVAA (ELBO, K = 1)
1.02
1.04
Pareto VAE (ELBO, K 4)
1.06
1.33
Pareto VAE (ELBO, K = 16)
1.24
1.40
Problem setup We first describe our problem for denoising task. This setting is comparable to the denoising criterion Im et al. (2016), which established an amenable loss formulation for Gaussian VAEs in the denoising task.
Let ~ Plara denote data samples, and let x' ~ r(æ' x) denote corrupted samples, where r(x' a)) is the corruption distribution (e.g. x' = x + with Puoiser I L e). The corrupted input x' is fed into the VAE, consisting of a conditional encoder 96 (z x') and decoder pe(x x z), yielding reconstruction output åt
Our goal is to recover x from x', thereby achieving a denoising effect. For theoretical consistency, we further introduce auxiliary corruption process - A that both sides define compatible joint distributions. This auxiliary process is only a theoretical device and does not need to be implemented in practice. We assume that the 7-power entropy H., (r)) finite. Moreover, impose a Markovian assumption, under which the following conditional independencies hold:
x',
2 1 H' T.
r(x'la (x)
go(=|x²)
po (x]
r(x' (z)
$
3'
Figure 5: Graphical model of the denoising task
Analogous to 7
we consider two joint manifolds involving the corrupted data x':
Mmodel (po x', (<=112) (o(2))) θ),
Manta (g0( Pamm ( (=(xx) Φ).
We then formulate the joint minimization problem for the y-loss, with the goal of deriving its explicit
form:
arg min (go(2,, po(, $E0,0E9
Derivation of the rylloss for denoising task To derive the denoising y-loss, first state useful lemma.
Lemma C.1 (7-power divergence for product-fom form densities), Let I (i(az) and =
IH gi (2x)) Then the y-power divergence from g to is
Dy(p1l =
- c. (p. 9))
II
33

---

Published as a conference paper at ICLR 2026
Proof.
1
Dy(p11 = cr(2,)) (D)
- 1 / II (x3))
Ih- (x(zs)
de (Aiost))
IIIT==I 8.1l247
1=1
j=1
-
II (, II
gi(20))
de II (f p))(s))
|1 l1q:ll1t7
7
ém1
(ʃ (p((s))
12.
doc +
II ==1
112ill144
(pi. ②:)
II Hy Ho)
Applying Lemma C.1 the 7-power divergence between the joint distributions factorizes as
C. p(2)) (un.aaaaheniiiiii..io.oooooahahd ?
H. (Psan (4
y
Since (r,r))) (r), this simplifies to
(po (o) (r(z'' x)) (Poma x²), p(2) po 22) - (Panaa (x) 2'))
(Poana (x) go x4) (2)) 2)).
The resulting divergence reduces the same form in 7b. which implies that the denoising 7-loss takes the form
(0, Ex-pan Ey (2))))
Ex || [o(=)]l 2o 2. Pat) + 20:3 (x)()))
The only modification compared to the standard y-loss 12 is the additional expectation over the corrupted input x' ~ r(- x)).
D DISCUSSION
Computation cost for training In contrast to models that involve numerical integration or iterative estimation of tail parameters, Pareto leverages the closed-form expression the 7-loss and the simple Laplace-Gamma Gamma mixture reparametrization of the symPareto distribution. This construction avoids computational bottlenecks and ensures stable training dynamies. Although the reparameteri- zation involves one additional sampling step compared to light-tailed VAEs (e.g.. Gaussian VAE or Laplace VAE), the runtime of Pareto VAE remains virtually identical to these models.
Tuning hyperparameter // Recent VAE models based on Student's distributions (Kim al. 2024 Bouayed et al., 2025) highlight the challenge of tuning the hyperparameter v. Pareto VAE encounters a similar issue. Previous works have shown empirically that selecting v from a suitable range (c.g.,
[2.5, 20] in [Bouayed al. 2025) is sufficient for stable performance in high- dimensional settings, although this criterion lacks theoretical guarantees, AAbii & Ohlsson (2020) proposes learning in the latent space, but learning may require additional constraints gradient estimators; in our preliminary attempts with additional structures to learn v, we did not observe consistent improvements.
To illustrate how performance depends on v, we vary the tail parameter v over a wide range and report the resulting PSNR values Table 5 To keep this appendix concise, report PSNR only;
34

---

Published as a conference paper at ICLR 2026
Table 5: PSNR sensitivity to the tail parameter V.
Model
== 1.1
// 2.1
3.1
5.1
10.1
ParetoVAE
24.268
23.424
22.670
17.590 19.403
19.443
VAE
17.812
17.878
19.044
Model
== 15.1
V = 20.1
V 25.1
V = 50.1
100.1
ParetoVAR
17.590
12.064
11.991 18.931
12.041
12.047
VAE
19.294
20.652
17.329
16.263
SSIM and accuracy exhibit the same qualitative trend. We also note that VAE requires whereas PAretoVAE only requires > 1.
In setup, Pareto attains its best performance at moderate within relatively narrow range compared to Student's i (-based VAFs. A plausible explanation is that symPareto already induces sufficiently heavy tails for relatively small v, so the robustness benefit saturates over shorter interval. Moreover, the f bi-based structure improves stability against sparse large deviations, further reducing the need search over a wide range of v. Nevertheless, there is still no theoretical foundation for optimal y selection, and we regard this as an important open problem in the study of heavy-tailed generative models.
Training strategies for regularizer Balancing the individual terms of combined loss function (reconstruction + regularizer) is a key challenge in training VAE-style models, as exemplified by methods such as (Higgins et al.) 2017) and (Fu et al. 2019). Similarly, it is important for ParetoVAE to achieve an appropriate balance between the reconstruction and regularization terms. To this end, explore various scheduling methods and empirically found that adjusting the regularization weight with logistic or linear function is more effective than using cyclic or cosine annealing. This effect particularly evident when combining loss terms with disparate properties, such & based reconstruction term with an &g based regularization term, or vice versa.
35


Reviews:
[{"review_id": "ciYL9iliFp", "reviewer": "Reviewer_gUPd", "summary": "This paper introduces the symmetric multivariate Pareto distribution (symPareto) and leverages it to build ParetoVAE, a new class of heavy-tailed variational autoencoders. The authors motivate the need for power-law latent distributions when modeling real-world data with extreme values or rare events, addressing well-known limitations of Gaussian-based VAEs. The work is grounded in a solid information-geometric perspective: the γ-power divergence is employed to derive a tractable VAE objective, circumventing difficulties associated with ELBO estimation for power-law distributions. Empirical studies across multiple domains (including noisy image reconstruction and denoising tasks) demonstrate improved robustness to outliers, better sparse-data handling, and enhanced recovery of fine-grained features relative to baseline VAEs and heavy-tailed t-VAEs. Overall, this paper makes a meaningful and well-motivated contribution to the study of robust generative modeling.", "questions": "Please see above weakness for questions to answer.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. Strong conceptual motivation for replacing exponential-family assumptions with heavy-tailed distributions in real-world generative modeling.\n\n2. Novel modeling component: introduction of the symPareto distribution as a multivariate power-law prior/encoder, extending beyond the Student’s-t VAE line.\n\n3. Theoretically grounded via γ-power divergence and joint minimization formulation, providing a tractable objective.\n\n4. Empirical benefits validated across several benchmarks, especially in high-noise and heavy-tailed settings. and Practical insights into when symPareto vs Student’s-t decoders are beneficial.\n\n5. Clear writing and strong experimental motivation, making the work accessible despite its mathematical basis.", "weaknesses": "1. Computational trade-offs not fully quantified: while the γ-loss improves tractability, a more detailed cost vs. ELBO comparison would help practitioners.\n\n2. Ablations on γ-power choices and sensitivity analysis would improve clarity on hyperparameter robustness.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "iVqqSbpAjJ", "reviewer": "Reviewer_dmgt", "summary": "The paper introduces the Pareto variational auto encoder. The model addresses the problem of modeling heavy tailed data within the VAE framework, which most typically use exponential family distributions. The paper offers a detailed and clear mathematical exposition for the multivariate Pareto distribution, uses information geometric tools to construct gamma flat manifold to analyze gamma power divergence, and presents detailed simulation results comparing to relevant baseline models on data with heavy tailed structure including graphs and language.", "questions": "I don't have major questions. I enjoyed the paper, thanks!", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "Strengths of this paper include:\n- The mathematical exposition was very nice. \n- The core motivation of modeling heavy tailed data was clearly borne out in experiments\n- The proposed model raised interesting challenges that were well resolved mathematically. \n- The approach was compared with a set of baselines that gave a clear sense of why the Pareto model was succeeding", "weaknesses": "Weaknesses include:\n- The introduction could have been more concrete sooner. Many of the arguments in the introduction would have been stronger if grounded in specific examples of data that are heavy tailed, specifically which models are under consideration, etc. The remainder of the paper answered these questions well. \n\n\nDetailed comments:\n- \"Incorporating robustness in generative modeling has enticed many researchers\nof the field.\" This sentence is not necessary. \n- \"To this end, ...\" To the end of enticing researchers into the field? \n- \"the γ-power divergence that naturally populates power-law families\" what does populate mean here?\n- \"Conventional VAEs typically employ exponential-family distributions, most notably the Gaussian, for their probabilistic model due to their mathematical tractability.\" This is typical but certainly not universally true. Therefore the next sentence \"To address these limitations,\" doesn't follow. It is necessary to list the exceptions to the convention and argue why they don't address the problem you want to solve. \n- \"exhibit heavy tails and extreme events\" What are the events you have in mind here? Is there evidence that standard VAEs haven't succeeded in modeling these events? Have modelers just not tried to use VAEs?\n- What do you mean precisely for a numerical integration to be intractable? (Line 044)\n- Table 1 is not particularly informative. Is there something important that I am missing there? \n- typo: exponentiall\n- typo: includde\n- Two papers that are not cited but may be of interest are the InfoVAE and Coupled VAE (cites below) which connect to entropic OT (itself connected to Information Geometry). \n\nZhao, S., Song, J., & Ermon, S. (2017). Infovae: Information maximizing variational autoencoders. arXiv preprint arXiv:1706.02262.\nHao, X., & Shafto, P. (2023). Coupled variational autoencoder. arXiv preprint arXiv:2306.02565.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "sP1kgmRJGK", "reviewer": "Reviewer_gNb7", "summary": "The manuscript discusses changing enforcing the Gaussianity of the distribution in the latent space of the variational autoencoders (VAE).\nThe manuscript proposes  ParetoVAE, a probabilistic autoencoder that minimizes the γ-power, instead of the Gaussian assumption as in the classical VAE, divergence between two statistical manifolds. Also Student's t and symPareto distributions are considered. These types of VAEs show good results in sparse, heavy-tailed data reconstruction and word frequency analysis.\nThe symPareto decoder enables robust high-dimensional denoising.", "questions": "What is the performance of the ParetoVAE in the case when representing images corrupted by Gaussian noise?\n\nMore references should have been provided for the statistical derivations from the Appendices, in the supplementary material.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The manuscript aims to develop a VAE that would be more appropriate for certain types of data.\n- The corresponding loss function for the new VAE is derived.\n- An application seems to be in denoising, when images are corrupted by Salt and Pepper noise", "weaknesses": "- Lack of theoretical development to justify why the new distributions are necessary in certain data representations. While the manuscript provides statistical distributions known in statistical analysis, it does not connect these with data reconstruction.\n- Experimental results are rather limited, especially in real data.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "19xNDi49og", "reviewer": "Reviewer_iLe7", "summary": "This paper proposes a novel VAE based on the information geometry of $\\gamma$-divergence together with a multivariate symmetric Pareto (symPareto) distribution. Since the symPareto distribution is heavy-tailed and promotes sparsity, it improves both the performance and robustness of VAEs. A variety of experiments validate the proposed method.", "questions": "- In the proposed VAE, symPareto is used for the encoder and the prior, and either a t-distribution (L2) or symPareto (L1) is used for the decoder. Which component is particularly effective? For example, if the goal is to induce sparsity in the latent variables, using only a symPareto prior may be sufficient. If the KL can be computed, it even seems acceptable to keep the encoder Gaussian. On the other hand, when observations are noisy or sparsity is required, using a t-distribution or symPareto for the decoder appears effective. For the image denoising experiment, for instance, it seems a t-decoder should handle the task; yet the proposed method achieves higher performance than t-based VAEs on all datasets. Where does this performance gain come from?\n- Using toy data, can you visualize how the latent variables behave under a symPareto prior? For example, as in Figure 6 of [1], train on a four-dimensional one-hot vector dataset and visualize the two-dimensional latent variables. I believe this would make the effect of the symPareto prior easier to understand.\n- VAEs often suffer performance degradation due to over-regularization by the prior. Does employing a symPareto prior resolve this issue? Moreover, it is known that the aggregated posterior (a infinite mixture distribution of the encoder) can address this problem [2]. Could your method be combined with such approaches?\n\n[1] Mescheder, Lars, Sebastian Nowozin, and Andreas Geiger. \"Adversarial variational bayes: Unifying variational autoencoders and generative adversarial networks.\" International conference on machine learning. PMLR, 2017.\n\n[2] Tomczak, Jakub, and Max Welling. \"VAE with a VampPrior.\" International conference on artificial intelligence and statistics. PMLR, 2018.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Introducing complex distributions into VAEs is difficult because one must compute the KL divergence between the encoder and the prior and also sample from the encoder via the reparameterization trick. This paper incorporates a complex distribution like symPareto into a VAE in a theoretically sound and tractable way. The empirical results are also strong.", "weaknesses": "Please see the Questions section.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)