COMET Flows: Towards Generative Modeling of Multivariate Extremes and Tail Dependence

1 Introduction

Modeling multivariate extremes involves two main challenges at odds with the usual assumptions of independence and Gaussianity in learning theory: (i) the challenge of modeling heavy-tailed marginal distributions [Resnick, 2007], and (ii) the challenge of modeling tail dependence within the joint distribution [Joe, 2014]. In this work, we propose COMET (COpula Multivariate ExTreme) Flows-a novel deep generative modeling architecture combining Extreme Value Theory and Copula Theory with normalizing flowsin an effort to more accurately model multivariate extremes exhibiting tail dependence.

Deep generative models estimate an unknown probability distribution from an input data for the purposes of density estimation and sampling. Popularized by advances in generative adversarial networks [Goodfellow et al., 2014], variational autoencoders [Kingma and Welling, 2014], autoregressive networks and normalizing flows, deep generative models offer a promising path forward as limits to labeled data constrains progress in supervised learning. Among deep generative models, normalizing flows set the gold standard of tractability, allowing exact sampling and exact computation of log-density by modeling a transformation between data vectors x ∼ p(x) and latent vectors z ∼ p(z). Given sufficient data and architectural complexity, normalizing flows are universal density approximators [Kobyzev et al., 2020].

Yet, normalizing flows still struggle to address challenges (i) and (ii) when modeling multivariate extremes (Figure 1). As discussed at length in [Wiese et al., 2019; Jaini et al., 2020], vanilla normalizing flows mapping heavy-tailed data distributions to light-tailed (e.g., Gaussian) latent distributions cannot be Lipschitz-bounded and hence become numerically unstable in practice. Likewise, viewing tail dependence among two or more dimensions as inducing a lowdimensional manifold structure in a region of feature space, we may reasonably expect vanilla normalizing flows to perform poorly in that region of space [Cunningham et al., 2020; Brehmer and Cranmer, 2020; Kim et al., 2020; Horvat and Pfister, 2021]. In this work, we address challenge (i) by employing a copula transformation to decouple heavy-tailed marginal distributions pi(xi) from the joint dependence structure in p(x), and address challenge (ii) by using a conditional perturbation of inputs during training to capture tail dependence in a stable manner.

The value of addressing challenges (i) and (ii) is apparent in a number of real-world applications, including climate science and meteorology (Figure 2). Extreme weather has become more frequent and severe in recent years due to climate change. Of particular concern is the projected increase in compound extreme events [Zscheischler et al., 2018], in which simultaneous multivariate extremes pose a greater risk to life and property than comparable univariate extremes spread over time. Accurate deep generative models of such compound extremes allow one to ask how likely an event is to occur, and may even be used to simulate extremes and their effects under various climate change scenarios [Puchko et al., 2019; Ayala et al., 2020; Klemmer et al., 2021].

Our main contributions are as follows:

1. We demonstrate that vanilla normalizing flows parameterized by Lipschitz-bounded neural networks are unable to capture heavy-tailedness and tail dependence in multivariate data;
2. We present COMET Flows, a novel architecture designed to capture asymmetric heavy-tailedness and tail dependence in multivariate data; and,
3. We empirically show that COMET Flows present a promising approach for modeling multivariate extremes using synthetic, climate, and benchmark data.

2 Preliminaries

Normalizing Flows Normalizing flows [Kobyzev et al., 2020; Papamakarios et al., 2021] offer a systematic approach to modeling multivariate probability distributions with deep learning. In particular, normalizing flows facilitate tractable sampling and density estimation under a d-dimensional probability distribution p(x) by parameterizing the transformation of each data vector x ∈ Rd sampled from p(x) to a latent representation z ∈ Rd following some simple probability distribution p(z) using an invertible neural network fθ : Rd → Rd . The network is trained to maximize the log-likelihood of observed data {xi}Ni=1 , which may be computed in closed form using the change-in-variables formula from calculus as a sum of the log-likelihood of the latent vector z = fθ(x) and the log-Jacobian determinant log |dfθ(x)/dx| of fθ:

Following training, the normalizing flow fθ may be used for density estimation by evaluating p(x) under (1). Assuming fθ is invertible, we may sample from p(x) by sampling z ∼ p(z) and considering f -1 θ (z).

Architecturally, normalizing flows present two challenges: (i) they must be invertible in order to allow sampling in the generative (z → x) direction, and (ii) they must allow for efficient computation of the Jacobian determinant |df(x)/dx| in order to be trained and to allow density estimation in the normalizing (x → z) direction. The model we propose in this work satisfies both (i) and (ii) with a coupling-layer based flow backbone [Dinh et al., 2017].

Extreme Value Theory Extreme Value Theory [Coles, 2001; Beirlant et al., 2005] characterizes the behavior of extremes in precise statistical terms, and may be used to model the realizations of a random variable greater than some threshold. Indeed, under relatively mild assumptions [Coles 2001, Theorem 4.1] one may show that the excesses of a univariate random variable x ∼ p(x) over some threshold µ follows the generalized Pareto (GP) distribution with probability density function:

Here, σ > 0 is a scale parameter, ξ ∈ R is a tail heaviness parameter, and p(x) has support on x ∈ [µ,∞) for ξ ≥ 0; alternatively, on x ∈ [µ, µ - σ/ξ] for ξ < 0.

Note that the GP distribution may be used to model left excesses (i.e., values under a threshold µ) as right excesses by transforming x → -x and µ → -µ. In this work, we utilize the GP distribution to model heavy-tailed marginals of the random vector x ∼ p(x) in both left and right tails.

Copulas Copulas [Nelsen, 2006; Joe, 2014] describe the dependence structure between uniform random variables, and in particular may be used to model the d-dimensional dependence between the components of a multivariate probability distribution after transforming the marginals to be uniformly distributed. Through the lens of Sklar's Theorem [Nelsen 2006, Theorem 2.3.3], a copula C(u) : [0, 1]d → [0, 1] is uniquely defined by a multivariate cumulative distribution function (CDF) F(x) : Rd → [0, 1] having univariate marginal CDFs Fi(xi) : R → [0, 1] through the relationship

If F is continuous with well-defined marginal inverse CDFs F -1 i (xi), then (3) may be equivalently expressed as

Taking partial derivatives of (3) with respect to each component gives us a similar expression relating p(x) : Rd → R to the copula density function c(u) : [0, 1]d → R which is often more convenient in practice when we are interested in density estimation:

In this work, we leverage copulas to model heavy-tailed marginals separately from the dependence structure of the components of a multivariate probability distribution.

Tail Dependence Tail dependence [Joe, 2014] characterizes the degree to which the behavior of two components of a multivariate probability distribution p(x) coincide at high (u → 1-) and low (u → 0+) quantiles. More formally, let p(x) be a multivariate probability distribution with joint CDF F andmarginal CDFs Fi, i = 1, . . . , d. Then the coefficients of upper and lower tail dependence λUij , λLij describing the tail relationship of components i and j are given by

As conditional probabilities, both λUij , λLij take values in [0, 1]. If λUij > 0, we say that components i and j of p(x) exhibit upper tail dependence; if λLij > 0, we say that components i and j of p(x) exhibit lower tail dependence.

In the extreme case of tail dependence between components xi and xj, it is possible that F -1 i (u) = xi = xj = F -1 j (u) for u < a ≈ 0 or u > b ≈ 1, such that there exists a local manifold structure of dimension Rd-1 embedded in Rd in the ij-th tail of p(x). More generally, if there exists such co-linearity between k components of x in the upper or lower tail of p(x), the support of p(x) may be a manifold of dimension Rd-k with measure zero in Rd . Such manifolds are difficult to model using normalizing flows given the requirement that the learned mapping fθ : Rd → Rd be invertible-in fact, if fθ successfully maps a manifold of dimension Rd-k ⊂ Rd to a latent space of Rd , it necessarily has an exploding Jacobian determinant. As a result, the problem of manifold learning with normalizing flows remains an active area of research [Cunningham et al., 2020; Brehmer and Cranmer, 2020; Kim et al., 2020; Horvat and Pfister, 2021]-yet to the best of the authors' knowledge, this work is the first to consider the connection between manifold learning, normalizing flows, and tail dependencies in multivariate probability distributions.

3 RelatedWork

Extremes and Normalizing Flows [Wiese et al., 2019; Jaini et al., 2020] prove the inability of normalizing flows to capture heavy-tailed marginal distributions, showing that any mapping of heavy-tailed distributions to light-tailed (e.g., Gaussian) distributions cannot be Lipschitz-bounded. To address this issue, [Wiese et al., 2019] propose Copula & Marginal Flows (CM Flows) which disentangle heavy-tailed marginals from the joint distribution by employing a hybrid deep-parametric marginal flow mapping inputs xi ∈ R to uniform space ui ∈ [0, 1], followed by a RealNVP-inspired [Dinh et al., 2017] copula flow, mapping the (non-uniform) copula distribution to a latent uniform distribution. However, the authors omit the implementation of marginal flows in experiments, consider only synthetic bivariate data on [0, 1]d , and suggest the use of Deep Dense Sigmoidal Flows [Huang et al., 2018] for constructing marginal flows despite the fact these are not analytically invertible and hence do not permit easily tractable sampling.

Alternatively, [Jaini et al., 2020] propose Tail Adaptive Flows (TAFs) which replace the latent isotropic Gaussian p(z) with a latent isotropic Student's t distribution having learnable degrees-of-freedom parameter ν > 1, arguing that mapping heavy tails of the input distribution p(x) to heavy tails of the latent distribution will maintain Lipschitzboundedness. However, the degrees of freedom parameter ν and the log-likelihood of the data under the model (which depend on one another) must be optimized simultaneously. Furthermore, by using an isotropic Student's t distribution in latent space, Tail Adaptive Flows implicitly assume symmetry in tail heaviness and hence are unable to capture asymmetric tail heaviness across dimensions.

Copulas and Normalizing Flows [Bossemeyer, 2020] develops a generalization of CM Flows to higher dimensions using vine copulas [Czado, 2019]. Orthogonally, [Kamthe et al., 2021] present Copula Flows motivated by the task of synthetic data generation, and utilize Neural Spine Flows [Durkan et al., 2019] to design a generative model able to handle both discrete and continuous data. In [Laszkiewicz et al., 2021], the authors propose to optimize a normalizing flow against a copula distribution in latent space instead of a standard isotropic Gaussian. Other relevant works at the intersection of generative modeling and copulas include [Ling et al., 2020], which considers deep estimation of Archimedean copulas, and [Tagasovska et al., 2019], in which vine copulas are fit upon the latent space of an autoencoder to construct a generative model.

Manifolds and Normalizing Flows A novel contribution of our work is to view tail dependence in multivariate probability distributions as inducing a local manifold structure in feature space. Recent works on manifold learning in normalizing flows entails density estimation on a d-dimensional manifold embedded in D > d dimensional space. Noisy Injective Flows presented in [Cunningham et al., 2020] utilize stochastic inversion for sampling. M-Flows are developed in [Brehmer and Cranmer, 2020] to carry out simultaneous manifold learning and density estimation. SoftFlow [Kim et al., 2020] estimates a conditional distribution of data given a perturbation noise level as input to the flow. This idea is advanced by [Horvat and Pfister, 2021] who derive theoretical conditions under which it is possible to recover the underlying manifold using a noising inflation-deflation approach.

4 Proposed Framework: COMET Flows

Vanilla normalizing flows struggle to capture (i) the heavytailed marginals and (ii) asymmetric tail dependencies in multivariate probability distributions. To address these weaknesses, we propose COMET (COpula Multivariate ExTreme) Flows, a normalizing flow architecture that decouples marginal distributions from the joint dependence structure within a multivariate probability distribution using a copula, models heavy-tailed marginals using a combination of kernel density and generalized Pareto (GP) distribution functions, and captures tail dependence by considering cases of λUij > 0, λLij > 0 as inducing a low-dimensional manifold in feature space.

COMET Flows has two main architectural components: a marginal transform fm : Rd → [0, 1]d mapping each marginal distribution pi(xi) to a uniform distribution such that the dependence structure among dimensions is expressed as a copula, and a copula transform fc : [0, 1]d → Rd mapping this copula distribution to a tractable latent distribution. Thus, COMET Flows define the normalizing flow (x → z) by the composition fc ◦ fm, and the generative flow (z → x) by the composition f -1 m ◦ f -1 c . We describe the details of each component in the following subsections, and present a schematic diagram of COMET Flows in Figure 3.

Marginal Transform fm In order to capture heavy-tailed marginal distributions pi(xi), we construct the marginal transform fm as a semi-parametric mixture in which the center of each marginal distribution is modeled using a 1dimensional Gaussian kernel density function, and the tails are modeled using a parametrically-fitted GP distribution motivated by the threshold exceedance paradigm of EVT.

More formally, let ai, bi ∈ [0, 1] be the lower and upper quantiles beyond which to apply parametric GP tails in dimension i, and let αi, βi ∈ R be the values in data space corresponding to ai, bi obtained through an inverse empirical CDF transform. Then fm : Rd → [0, 1]d is defined to be a decoupled mapping acting on each dimension i = 1, . . . , d of its input independently of other dimensions j ̸= i, with i-th component fm,i given by

Here, Ai and Bi are the rescaled CDFs of parametric GP distributions fit to the data in the left and right tails of pi(xi), respectively, and Ki is the rescaled CDF of a 1D kernel density estimate fit to points in the center of pi(xi):

Here, ˜ Ai, ˜Ki, and ˜Bi refer to the original CDFs fit to points xi < αi, xi ∈ [αi, βi], and xi > βi, respectively. By construction, the inverse marginal transform f -1 m,i : [0, 1] → R is well-defined by the rescaled inverse CDFs (A-1 i , K-1 i , B-1 i ) on (ui < ai, ui ∈ [ai, bi], ui > bi).

Note that the only learnable portion of fm,i are the parameters of the GP distributions used to model each tail, which can be fit using maximum likelihood estimation (e.g., using scipy.stats.genpareto.fit). This is by design, and vastly simplifies the learning dynamics of the end-to-end fc ◦ fm model: the GP parameters of each fm,i may be fit and fixed before fc is trained. This comes in contrast to the approaches of [Wiese et al., 2019; Bossemeyer, 2020], who allow each fm,i to be learnable flows but struggle to successfully train the coupled fc ◦ fm model.

Copula Transform fc Given a copula distribution in the unit hypercube [0, 1]d , we construct fc to be a normalizing flow with codomain Rd . To guarantee that the generative direction (z → u) maps vectors into the unit hypercube, we apply a logit transform as the first layer of fc such that the last layer in the generative direction comprises a sigmoid transform. As the logit transform has no learnable parameters, we use fc to refer only to the learnable transformation applied after the logit transform in what follows.

The learnable transformation fc is implemented as a coupling layer-based RealNVP [Dinh et al., 2017], repeatedly alternating dimensions 1:k and k + 1:d and applying

in each layer, where s, t : Rd → Rk , 1 < k < d, and ⊙ denotes elementwise multiplication. With Ik representing the k × k identity matrix and ∗ being any arbitrary (d - k) × k matrix, the Jacobian of each coupling layer takes the form

and has a readily-computable log determinant ∑di=1 s(x1:d)i. Because this layer is invertible with tractable log-Jacobian determinant regardless of the form of functions s and t, we implement s and t as concatsquash [Onken and Ruthotto, 2020] networks conditioned on a contextual input σ:

To justify this construction, we recall that modeling multivariate extremes presents two main challenges: capturing (i) heavy-tailed marginal distribution and (ii) asymmetric tail dependencies. The marginal transform fm presented previously addresses (i), leaving the copula transform fc to address (ii). We accomplish this by viewing tail dependencies as inducing a manifold structure within subspaces of the unit hypercube, and draw inspiration from SoftFlow [Kim et al., 2020].

More precisely, we interpret the conditioning input σ to the concatsquash networks s, t to be a noise level, and perturb inputs x to s and t by x + ϵ where ϵ ∼ Nd(0, σId) for a noise level chosen uniformly at random σ ∼ U[0, σmax]. This noise level σ is passed to the flow as a conditioning contextual input through networks s and t throughout training, allowing the network to learn the conditional distribution p(x|σ) while sidestepping the exploding gradients introduced by the manifold structure of tail dependencies. Yet, this approach still allows the network to recover such tail dependencies at inference time by setting σ = 0.

5 Experimental Evaluation

5.1 Baselines

RealNVP A coupling layer-based RealNVP [Dinh et al., 2017] implementation optimized against an isotropic Gaussian latent distribution.

SoftFlow A modified SoftFlow [Kim et al., 2020] implementation in which the FFJORD [Grathwohl et al., 2019] backbone is replaced with a RealNVP backbone.

Tail Adaptive Flow (TAF) A modified TAF [Jaini et al., 2020] implementation composed of a RealNVP flow backbone with an isotropic Student's t latent distribution, parameterized by degrees-of-freedom parameter ν. We consider ν to be a fixed hyperparameter as no algorithm is presented to learn ν in the original work.

Copula Marginal Flow (CMF) A modified CMF [Wiese et al., 2019; Bossemeyer, 2020] implementation in which the unspecified tail distribution and non-analytically-invertible Deep Dense Sigmoidal Flows [Huang et al., 2018] used for fm are replaced with an analytically-invertible Gaussian KDE and generalized Pareto tail mixture model for fm.

5.2 Data

Synthetic Data We design a synthetic dataset exhibiting heavy tails and asymmetric, heterogeneous tail dependence by sampling from a uniform distribution U on [0, 1]8 with probability 0.95, and from independent univariate GP marginal distributions X parameterized by (µ, σ, ξ) = (0, 1, 1) with probability 0.05. These samples are combined such that x1, x2 exhibit a heavy upper tail with tail dependence; x3, x4 exhibit a heavy lower tail with tail dependence; x5, x6 exhibit heavy upper and lower tails, each with tail dependence; and x7, x8 are perfectly collinear with heavy upper and lower tails. We generate 200,000 vectors in x ∈ R8 for training, 25,000 vectors for validation, and 25,000 for a heldout test set. Key bivariate marginals of this distribution are shown in the rightmost column of Figure 4.

Climate Data To evaluate the practical applicability of COMET Flows, we apply them to a subset of CLIMDEX 1 extreme climate indices [Sillmann et al., 2013]. Motivated by the projected increase in extreme temperatures and rainfall under climate change, we extract data from a model run under the RCP 8.5 scenario of continued high greenhouse gas emissions, and aim to model the joint distribution of such extremes over North America. We consider the region encompassing 25◦N-75◦N and 170◦W-50◦W and simulation data from the years 2091-2100, in which the effects of climate change on extremes will be most pronounced. We choose four indices of extreme temperature (TXx, TNx, TXn, TNn), along with four indices of extreme precipitation (prcptot, r1mm, r10mm, and r20mm), then add latitude and longitude to form a 10dimensional dataset characterizing joint extremes and their spatial distribution. Indices are reported at every 0.94◦ latitude and 1.25◦ longitude, giving a total of 97 × 53 = 5,141 vectors x ∈ R10 per year. We use data from 2091-2098 for training, 2099 for validation, and 2100 for testing. A subset of this data was presented in Figure 2.

1 https://www.climdex.org/

Benchmark Data To contextualize the performance of COMET Flows, we consider the POWER and GAS benchmark datasets originally presented in [Papamakarios et al., 2017]. We follow the standard preprocessing steps provided in the MAF GitHub 2 repository. Note that we do not expect our baseline results to match state-of-the-art negative loglikelihood values reported in [Kobyzev et al., 2020] due to the simplified and adjusted architectures we implement.

5.3 Results

Synthetic Data As shown in Figure 4 and Table 1, the copulabased CM and COMET models dominate in both the normalizing (x → z) and generative (z → x) directions on synthetic data due to their functional form: by construction, they are naturally able to account for the asymetrically heavy-tailed margins of training data where vanilla normalizing flows cannot, and where TAF are forced to assume symmetry in tail heaviness. In particular, COMET's ability to capture the lowdimensional manifold structure in the upper tail region of feature space is highlighted by its qualitative performance in the generative direction (Figure 4), and its ability to avoid the instabilities and exploding gradients associated with manifold structures in feature space by learning perturbed manifolds conditional on noise is apparent by its quantitative performance in the normalizing direction (Figure 1). COMET's edge over CMF is most evident in its ability to capture tail dependence, though both significantly outperform TAF, SoftFlow, and RealNVP.

Climate Data COMET's advantage in the normalizing direction is more evident in this experiment. Key bivariate marginals in Figure 2 appear to exhibit local manifold structures (e.g., between tails of TXn and TNn), explaining COMET's superior performance. This observation suggests that explicitly accounting for such manifold structures in the joint tails of multivariate probability distributions is profitable in practice, particularly when redundancy across dimensions is likely to exist as an artifact of the latent physical processes generating the observed data.

2 https://github.com/gpapamak/maf

Benchmark Data COMET is not as performant as CMF on the benchmark data in the normalizing direction, likely due to the absence of heavy-tails and tail dependence in POWER and GAS. Nevertheless, COMET outperforms TAF, SoftFlow, and RealNVP, suggesting that it remains competitive even on datasets it is not explicitly designed to handle.

6 Conclusion

Modeling multivariate extremes involves two fundamental challenges: (i) modeling heavy-tailed marginal distributions, and (ii) modeling tail dependence within the joint distribution. To address these challenges, we propose COMET Flows, a novel normalizing flow architecture which draws inspiration from Extreme Value Theory and Copula Theory, and which views tail dependence as inducing a local manifold structure in feature space. We demonstrate the efficacy of COMET Flows through experiments on synthetic, climate, and benchmark datasets, observing competitive performance in sampling and density estimation.

Acknowledgements

This research is partially supported by the National Science Foundation under grant IIS-2006633. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government. The authors thank Asadullah Hill Galib and Tyler Wilson for their input in discussions regarding the paper.