Diffusion and Flow-based Copulas: Forgetting and Remembering Dependencies
Authors: David Huk
Theodoros Damoulas
Venue: ICLR 2026 Poster (2026)
Keywords: Copula estimation
dependence modelling
diffusion
non-parametric copula
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)
TL;DR: We model copulas based on the principles of diffusions and flows with marginal preserving processes.

Abstract: Copulas are a fundamental tool for modelling multivariate dependencies in data, forming the method of choice in diverse fields and applications. However, the adoption of existing models for multimodal and high-dimensional dependencies is hindered by restrictive assumptions and poor scaling. In this work, we present methods for modelling copulas based on the principles of diffusions and flows. We design two processes that progressively forget inter-variable dependencies while leaving dimension-wise distributions unaffected, provably defining valid copulas at all times. We show how to obtain copula models by learning to remember the forgotten dependencies from each process, theoretically recovering the true copula at optimality. The first instantiation of our framework focuses on direct density estimation, while the second specialises in expedient sampling. Empirically, we demonstrate the superior performance of our proposed methods over state-of-the-art copula approaches in modelling complex and high-dimensional dependencies from scientific datasets and images. Our work enhances the representational power of copula models, empowering applications and paving the way for their adoption on larger scales and more challenging domains.

Published as a conference paper at ICLR 2026
DIFFUSION AND FLOW-BASED COPULAS: FORGETTING AND REMEMBERING DEPENDENCIES
David Huk
Theodoros Damoulas
Department of Statistics
Department of Statistics, Department of Computer Science
University of Warwick
University of Warwick
United Kingdom
United Kingdom
David. Hukewarwick. .u
T. D..Damolla as @aaeaaaasssaaaaaaa. @wwrwick.
ABSTRACT
Copulas are a fundamental tool for modelling multivariate dependencies data, forming the method of choice in diverse fields and applications. However, the adoption of existing models for multimodal and high-dimensional dependencies is hindered by restrictive assumptions and poor scaling. In this work, we present methods for modelling copulas based on the principles of diffusions and flows. We design two processes that progressively forger inter-variable dependencies while leaving dimension-wise distributions unaffected, provably defining valid copulas all times. We show how obtain copula models by learning remember the forgotten dependencies from each process, theoretically recovering the true copula at optimality. The first instantiation of our framework focuses on direct density estimation, while the second specialises in expedient sampling. Empir- ically, we demonstrate the superior performance of our proposed methods over tfate-of-the-an art copula approaches in modelling complex and high-dimensional dependencies from scientific datasets and images. Our work enhances the repre- sentational power of copula models, empowering applications and paving the way for their adoption on larger scales and more challenging domains.
INTRODUCTION
Given a collection of d continuous random variables, a simple model for their joint probability density function is the product of the corresponding d univariate densities (Peterson, 1987). Of course, this omits vater-variable dependence, which given by d-dimensional density supported on the unit hypercube with uniform marginal distributions' called a copula (Sklar, 1959). Indeed, the copula uniquely and exactly represents the inter-variable dependence, unlike correlation or mutual information (Geenens, 2023), fully disentangling the marginal behaviour from the joint.
This disentanglement enables a modular approach for multivariate modelling: first, model the uni- variate variables independently, and second, model their dependence with copula. This is pre- dominant in applications where the marginal behaviour is known or requires specific properties, as copulas respect these by construction. Copulas produce well-specified probabilistically calibrated marginals to represent extreme events accurately; properties that are crucial for state-of-the- art mod- els across fields such weather forecasting (Cong & Brady, 2012; Huk al., 2023), hydrology (Salvadori & De Michele, 2007; Van de Velde et al., 2023), risk management (Kole et al., 2007; Dewick & Liu, 2022), retail (Salinas al., 2019), and causal inference (Evans & Didelez, 2024). Another popular application is to induce dependence to correct pre-existing independent models, as used in generative modelling (Tagasovska et al., 2019; Liu et al., 2024), multi-agent imitation learn- ing (Wang et al., 2021), uncertainty quantification (Chilinski & Silva, 2020; Tagasovska et al., 2023), Bayesian non- con-aaraeetries (Huk al., 2024), and classification calibration (Panda & Garain, 2025). Copulas are also invariant to monotonic transformations of variables, benefiting multi-objective Bayesian optimisation (Park et al., 2024). Finally, they isolate inter-variable dependencies in an
Work partly done during an internship within team at Unilink Software Ltd.
In this work, marginal distributions refer to univariate dimension-win marginals, unlike (multivariate) time-marginals common in the literature diffusions.

---

Published as a conference paper at ICLR 2026
SDE-driven forward process
Classification-iifusiii Copula
Density evaluation in one function call via class probabilities dependence
DO
ts0.16 LIITE
t40.24 PERTIEE
E20
P(time 0lu)
c( (
P(time 2c
Scores-based sample generation
log Pttim = t
V Vuju llo Pitime
Dependent copula
Only dependence forgotten
Independent copula
Reflaction based
Reflection Copula
random Initial velocity
Predict given location and
t=0.00 10000
t=0.16
tw0.24
t=0.32
v uu, 1) = E |1 = u]
Follow reversed valocity prodictor from uniform generate dependent samples
Un/f (t, - c(u)
Figure 1: Overview of proposed copula models. We design forward processes to forget inter- variable dependencies but preserve dimension- -ise marginals. Our cessificaatonndiffssi copula and reflection copula learn by remembering the forgotten dependencies of these processes.
interpretable way, facilitating the study of neuron spikes in the brain (Berkes et al., 2008; Verzelli & Sacerdote, 2019), increasing the flexibility of Bayesian networks (Elidan, 2010), enabling synthetic data generation (Patki et al., 2016; Sun et al., 2019) with preacccpeeereiingggg applications (Griesbauer et al., 2025), and empowering conformal prediction (Sun & Yu, 2023; Park & Cho, 2025).
high- Yet, current copula models impose restrictive assumptions limiting their effectiveness
dimensional and non-trivial data. For instance, Gaussian copulas (see Ch. in Hofert et al. (2018)) only capture diagonal symmetric dependencies, while vine copulas (Nagler & Czado, 2016) omit parts of the dependence and explore an exponential model space with the data's dimension. Existing deep copula models suffer from mode collapse issues (Hofert et al., 2021; Janke et al., 2021) and struggle when sampling from multi-modal and high-dimensional data (Huk et al., 2025).
To address these challenges, we design copulas representing dependence purely through the for- malism of diffusions and flows (Song & Ermon, 2019; Lipman et al., 2023). The first step of our analysis is to design appropriate forward processes, motivating our first research question:
Q1: How to design a stochastic process that only forgets dependence?
We answer it by introducing two processes that maintain their marginal distributions but progres- sively forget the inter-dimension dependence. The first is grounded in stochastic differential equa- tions with time-invariant marginal distributions. The second consists of endowing samples with random velocities and reflecting them in hypercube through time. We show that, when applied to copula data, these processes provably preserve uniform marginal distributions, i.e., they define valid copulas at all times with vanishing dependence. This insight motivates our second question:
Q2: Can we remember the forgotten dependencies obtain copula densities and samples?
Motivated by the main uses copulas, for each process, we propose a model with a separate goal: in Sec. 3, we present an effective density estimator we term the classifiation-ddifision copula, and in Sec. 4, we introduce the reflection copula as a generative dependence model suited for expedient sampling (See Figure 1). Theoretically, show both models recover the ground truth copula at optimality, enabling density estimation and sampling. We summarise our contributions as follows:
1. We introduce processes to purely forget the dependence of data with time (Proposition 2, Proposition 7), defining a spectrum from a dependent to an independent copula.
2. We propose glassificatoon-ddiffssinn copula which provably learns the dependence in data (Theorem 5), outputs copula densities in a single model evaluation (Proposition 3), and leverages diffusion algorithms for sampling (Proposition 4).
2

---

Published as a conference paper at ICLR 2026
3. We introduce the reflection copula to learn the probability path of forgotten dependencies and generate samples by following this path in reverse (Proposition 8).
4. We achieve state-of of--th--rrttt results dependence modelling. surpassing existing classical and deep copulas on variety of complex real-world data. Our copulas are the first instance able scale to such high-dimensional and multimodal dependencies (see Fig. 5).
BACKGROUND ON COPULA MODELS
Notation and setup. We denote observed data X (x''...,,,") X C Rd coming from
cumulative distribution function (CDF) P(-) with probability density function (PDF) p(.). and refer X' the data scale. We assume that dimension- wise marginal CDFs {P)())
( i= and PDFs (u' ..... E [0. 17 with
(O) 1 are given. Next, we denote copula observations as u' (zx),, calling this the copula scale. Finally, for 4 the univariate standard Gaussian CDF, we also work with variables z = (z',.... ; =)) IRd for 2 (w'( which we say are on the Gaussian scale. We always denote time in subscripts, keeping superscripts for the dimension.
The task. We are given data x and know marginal distributions either through assumptions on the model or by pre-emptively estimating them. These marginal distributions are used to transform the data x to copula scale observations u encapsulating dependence. The copula modelling task then consists of providing a valid copula model for produce samples or likelihood evaluations.
An example of this multivariate time series forecasting, where one strives to describe the future joint distribution over many variables, such the state of multiple agents in imitation learning, Apdx B. Such uses of copulas are a useful way of dealing with the exceedingly large dimension- alities of multivariate timeseries by decomposing them into marginal forecasts which are stitched together through copula, decomposing the original problem into manageable subtasks. First, one estimates the individual distributions for each agent separately, possibly utilising expert knowledge to ensure calibration, and secondly, unites all univariate distributions with a copula to prescribe given dependence structure.
Copula-based modelling. A copula is joint probability distribution over univariate marginal distributions, owing its name to the Latin term for a link or bond. Indeed, the copula links together independent marginals, infusing them with dependence to form valid joint distribution:
Theorem 1 (Sklar (1959)). Let P be a d-dimensional distribution on X C Rd with continuous marginal cumulative distribution functions p .2..... pd. Then there exists unique copula dis- tribution C: [0, 1]d + [0, 1] such that for all x x2 ,..,, ....,,))= Rd:
P(x) P(x'....,,z" x²)) C(P.
And if a probability density function p X + R is available:
pla',...,, (x²) ..... (x*) c(P1 (x*),..., pd(x)))
(I)
where p' (x*)) are the marginal PDFs, and [0, + (0, oc) is the copula PDF.
This theorem allows for two-step estimation process of joint density: First, model marginal distributions, and second, model the dependence between them with a copula. This grants great flexibility, as marginals can be tailored to the specific problem in mind, while the choice of the copula model can prescribe particular joint behaviours. Eq. (1) used obtain density evaluations, and sampling is done by applying the inverse marginal distributions p pr (w') ) to copula samples u. Correctly modelling the copula therefore of utmost importance to faithfully represent data.
Popular copulas used in machine learning Most works adopt one of two models due to their ability to quickly produce densities and samples, even in high dimensions. Firstly, the Gaussian copula is generally the first copula to be implemented in applications, and is used in Elidan (2010); Patki et al. (2016); Salinas et al. (2019); Wang et al. (2021); 1); Huk et al. (2023); Panda & Garain (2025). Second, the vine copula (see Aas al. (2009): Nagler & Czado (2016)) is the de facto model for copula density evaluation and sampling, being used in the majority of cases not using a Gaussian, in Tagasovska al. (2019; 2023); Huk et al. (2024); Park et al. (2024); Park & Cho
3

---

Published as a conference paper at ICLR 2026
(2025); Griesbauer et al. (2025). Other model classes, such Archimedean copulas, have seen some applications (Liu et al., 2025) but remain limited to diagonal and symmetric dependencies and scale poorly to high dimensions. Their deep variants have overcome the large computational complexity of early works (Ling et al., 2020), and have been shown to effectively scale to dimensions as high as 20 (Ng al., 2021; 2022). While the work of Ng al. (2021) remains bound to diagonal relationships of an Archimedean copula, Ng et al. (2022) gains flexibility via the Archimax class. Fully general deep copulas have been explored in Hofert et al. (2021); Janke et al. (2021) generative metric- matching networks and in Kamthe et al. (2021) as normalising flows.
Current state of the art in copula models. Recently, Huk et al. (2025) showed that any copula density is equivalent to a binary Bayes optimal classifier between dependent and independent data with the same marginal distributions. Their framework yields flexible models, achieving state-of- the-art performance in copula modelling for complex dependencies. Unfortunately. it only estimates densities, meaning sampling requires expensive Markov Chain Monte Carlo scaling as (dd1/3) with the data dimension d, making this method inappropriate for multimodal high-dimensional settings.
In the next Section, we propose our first novel framework for copula modelling, by relying on diffusion processes specifically designed to forget dependence. In doing so, we depart from the ratio copula's 's dichotomy between dependence and independence to study a continuous range of decaying dependencies. This framework shift lets us leverage diffusion-baser sampling, making our method scalable to high-dimenniioall and multimodal dependencies (Azangulov al., 2024).
CLASSIFICATION DIFFUSIIO COPULA
Class probabilities based on dependence
Copula density evaluation
5.3
0.30
dd
s1
6.00
P(t
T 2
5.0
P(t
7
Pt
T 2
2)
Til
(w))
to Bausslan
P(
T. 1z
Copula score-based samplins lt T,
31
V. (o)
Pog
-100
CLTE)
Figure 2: Classifir-DDffuusii Copula. We map copula data to Gaussian scale data z, which we apply an Ubieebbee process up to a time t. We train a multinomial classifier cae to identify the diffusion time based z's dependence. This classifier recovers the copula density.
In this Section, we develop dlassificationddiffiiiii copulas whose strength lies in performing accu- rate density evaluation in a single function call, while also generating samples through diffusion- based sampling methods. A diagram of the approach is shown Fig.
3.1 FORWARD PROCESS TO FORGET DEPENDENCE
Since copulas density ratios (Huk et al., 2025), their densities preserved under diffcomor- phisms (Choi et al., 2021). Thus, noticing that x, u and z stem from invertible mappings of cach other, model the copula of x through the copula of given they share the same dependence. This lets us define a process on Rd instead of [0. 1]d. We augment our variables z with a time dimension denote 2, as the variables following a density p. for [0, oo), with the original data (on the Gaussian scale) at & = O. In the following, we formalise and answer our first research question:
Q1 (a): How to define a process which preserves p = p for any dimension at all times oo?
t, 2 0, while diffusing the copula e, of 2, independence
4

---

Published as a conference paper at ICLR 2026
OrasteinUUhlenbeck the Gaussian scale. As us marginally uniform, is marginally standard Gaussian with density (2;; 0, 1) by the probability integral transform. This insight motivates the following process (2) . 24) for t 0:
dz -2 dt + 22d,,
(2)
where B, denotes Brownian motion. We Eq. (2) the forward process of our model as it moves the data Zo forward in time ž₁ This an Censten-UUlleneee dimensions, and has the property of having the standard d d-dimensional Gaussian as its stationary distribution in the limit of time + oo (Särkkä & Solin, 2019), fulfilling our first desiderata on sampling, the OU process enjoys faster convergence to its limiting
(OU) process that is independent across
marginals. For diffusion-baseee distribution (Brešar & Mijatović, 2024) and easier to learn from (Reu et al., 2024) compared to other stochastic processes, which motivates our choice by Theorem 2 of Chen et al. (2023). It is also easy to simulate due to being analytically tractable. To fulfil our desiderata on the copula, we inspect what happens to u, + ((z) the copula scale, yielding the following result. The proof is given in Apdx. A.I.
Proposition (Convergence to independence copula), For copula samples Zo on the Gaussian scale, consider un ((z)) with z as defined through the process in Eq. (2). Then, this process maintains uniform marginal distributions at all times:
~ U), VI 2 2,0E ()...e(
Moreover, denote by C the copula of zı under Eq. (2). Then, c converges to the independence copula as t + oc in the Kullback- Lulbback-Leibler divergence with rate O(e-2t)
For complex distributions, a process with a limiting distribution closer to the data distribution in the Kullback-Leible divergence is preferable (Rhodes et al., 2020; Chen et al., 2023). Hence, modify Eq. (2), by adding a correlation matrix Σ:
dz; = -2 dt + V2 dB,, with ((,)) = 1 for E (d.....d).
(3)
This process converges to (0, E), and we prove its convergence in Proposition 9 of the Appendix. In what follows, we will focus on Eq. (2), and defer the derivations for Eq. (3) to Apdx. A.5.
3.2 REMEMBERING DEPENDENCIES FOR COPULA DENSITIES AND SAMPLING
We present the classification citiffccation-diffsiioo copula to learn the true copula density from process (2). As standard diffusion models require expensive numerical integration (Song al., 2021) for likeli- hoods, we instead answer the following question to facilitate density-reliant copula applications:
Q2 (a): How to obtain direct copula densities without foregoing effective sampling?
Classifying diffused data based on remaining dependence. We begin by discretising the for- ward process of Eq. (2) into times? Ti, T2...., with T; = 0, Te co, each with respective copula
densities denoted pr, . ......., PT.. We define the cuasiffiaationdiiffiiii
function which, given a data point z, outputs a vector of probabilities corresponding to the likelihood of z origi- -i=
nating from each of the i times Cdc:
Rd + (h E [0. 1) Σ h² =1]. with:
(z) = (P(t = 2)) P(t Tal z)).
(4)
The obtained convergence rates of the forward process in the previous Section motivate our choice of time steps Ti. ...... T3 in practice to achieve proportional changes in the dependence in line with our convergence rate Prop. 2, see Apdx for further details. B. This vector (4) of correct class probabilities can be linked to the copula density as follows:
Proposition 3 (Cae copula density evaluation). For the classification task as defined in Eq. (4), the true copula density c(u)) at a given value [0, 1]d is equal to the following probability ratio:
c(u)= IP(L = Tilz = h- (w)) (P(L Tulz= ==12=)-)(()),
where & (u) is the inverse standard Gaussian CDF applied dimensions-wise to
We take Ti as large time when the forward process expected he stationarity, see Apdx. B.3.
5

---

Published as a conference paper at ICLR 2026
We include proof in Apdx. A.2. Next, show how express the score of copula at given dependence level in terms of the class probabilities alone, enabling score-based sampling. The proof is in Apdx. A.3 and advances the technique of Yadin et al. (2024) to copula densities.
Proposition 4 (Cae copula score-based sampling). Consider the classification task as defined in Eq. (4), and let cs be the copula of u, ((Z.) with & the standard Gaussian CDF applied element- wise, and defined by the process of Eq. (2) at time E {T....., Tx}. Then, for @ denoting element-wiss multiplication, the true copula score Vu log (o(g),() is given by:
Vu log (c(()) = =(u) @ 1ogo = T=|z= = Vx log(P(( = Talz = =1)== ===-i(w)
where for as the Gaussian pdf. we have w(u)= == o( 1(11))
1(u4))
employ powerful diffusion model algorithms to obtain copula samples through Langevin dynamics (Song & Ermon, 2019), see Alg. 1 in Apdx. B. Importantly, under the popular assumption that the true distribution is supported on a manifold with dimension < d, diffusion-basee higher dimensions (Azangulov et al., 2024). This is in stark contrast current deep copula methods (Huk et al., 2025) requiring Hamiltonian Monte Carlo in the ambient dimension d.
Consequently, given accurate model for the class probabilities,
sampling scales more efficiently to datasets of
Training the cde model. Having demonstrated how to leverage class probabilities perform cop- ula density estimation and copula sampling, our next contribution is to prove how to estimate such a model correctly. To obtain both accurate densities and high-fidelity samples, we adopt a mixture loss between cross-entropy for the dependence level and mean squared error for the score, which has been empirically successful in previous work (Yadin al., 2024). We provide novel result that theoretically shows its correctness with proof in Apdx. A.4.
Theorem (Loss function for Cac). For the classification task defined in Eq. (4), and for a family of multinomial models Cac(2; 0) RR + {h E [0, 1]41 Σ h² = 1} indexed by parameters 8, assume there exists set of parameters O" = {0 Cat 2; 0) = (P(L Tilz =2),...,Pi P(I = Txlz = =)) },
for all z E IRd. Then, for any weight a > 0, {0 A = arg ming Le,, (6) = O" for the loss:
(0) =a
Ezppr. log de (z;
1=1
+
E Eer, ~pr e>A1(0.)
[EEa(Caoe To To .2, + V VI- -2)))))))- - =E|?
where de (z: (2: 0) denotes the s-th component of the vector caol (2; 0), and with
(k) és (s(cac(z;0)) = e (V log Cd (z: f) - Vz log cdc (z; 0) +))
We choose a such that both loss terms are of the same magnitude, following Yadin et al. (2024), with ablation study on in Apdx. B.4. As this result only holds optimality, we validate sample properties in practice in Apdx. B.5 via statistical tests, rank diagnostics and calibration metrics, with visualisations in Apdx. B. We give full overview for training the cde in Alg. of Apdx. B.
4 REFLECTION COPULA
As many copula applications only require sampling, we design a generative copula model through processes that forget dependence, inspired by flow methods. We provide a diagram in Fig. 3.
4.1 FORGETTING DEPENDENCE ON THE COPULA HYPERCUBE
In distinction from Sec. 3, here directly design our method on the [0, 1] copula scale. To obtain a stochastic path from a dependent to an independent copula, we answer:
Q1 (b): How to forget dependence directly on J0, 1)dt and maintain uniform marginals?

---

Published as a conference paper at ICLR 2026
Velocity space
Reflected to copula [0, 17
Velocity prediction
0.8
0.05
0.4
0.2
0.0
Copula observations
Copula space [0, 1) Velocity of sample u time
Predicted Velocity
Figure 3:
Reflection copula design. (Left panel) Copula data in black given red velocities to
move with time. (Middle panel) Trajectories are reflected from IR to [0, 1] following mirrored blue outlines. (Third panel) Reflected trajectories diffuse the copula with time. The reflection copula then learns a velocity predictor vuu, t) == E[ve/u u = ul for the average velocity, needed for sampling.
We start with copula observations u E [0, 1] and again augment our space with a time dimension t [0, oo). We let ug and endow samples with velocities Vo NN Ia). We define process for the evolution of sample-velocity pairs (ut, ve) through time via dimension-wis reflections.
Definition 6 (Reflection process). For initial samples ug E [0, 1]d and velocities Vo ~ (O),,
define the univariate hypercube reflection operator R. IR
R + [0, 1] x IR
(x (x]. (11 - x1 [x], -),
if | is even,
R(x, y)
(5)
if [x] is odd
The reflection process for sample-velocity pairs (u,, vz) time t is defined dimension-wise
(u vi) = R(ub vo. vo)..
(6)
The process starts from copula distribution, selects a direction according to Gaussian, and bounces around the [0, 1] hypercube through time. The only randomness in this process comes from the velocity, so the trajectory deterministic once Vo is sampled, see the diagram in Fig. 3. We prove that this process defines valid copula at all times as it converges to independence.
Proposition 7 (Reflection converges to independence). For any initial point uo [0, 1] with vo (0, Ia) as in Definition 6, the reflected point u, converges in distribution to the independent copula on the hypercube as t - co. Further, if ug follows copula distribution, then u, also follows some copula distribution e with uniform marginals for any time 0.
We give proof in Apdx. A.6. Note that the copula O generally different from the true copula and will have less dependence with time. As such, this process creates stochastic interpolation between dependent copula and independent one with time, gradually reducing dependence.
4.2 REMEMBERING DEPENDENCE BY PREDICTING VELOCITIES
To obtain a generative copula model from this process, answer the following:
Q2 (b): How to generate samples by learning from the reflection process in Definition 6?
Intuitively, as velocities are the only randomness in the system, if we learn the expected velocity, we should learn the system's 's average behaviour. We leverage connection between velocities and the probability path governing the change in copula distributions through time, as is shown in the following result. We include proof of its applicability the reflection process in Apdx. A.7.
Proposition 8 (Holderrieth al. (2024)). Consider the reflection process introduced in Eq. (6) applied to copula samples un ~ c(u) with velocities Vo ~ NN(,, Ia). Let u (u, = E[v, = u] be the expected velocity time f 2 and location [0, 1] and denote by er the copula of

---

Published as a conference paper at ICLR 2026
samples uT at time T. Then, starting from value ur cr(uu) and following the probability path given by the ordinary differential equation
di Ur - u² (u,()
(7)
backwards from T 0, the destination point comes from the initial copula distribution Ug c(u).
In other words, the expected velocity of this process is enough generate samples. However, we first require a sample Uy at time T from our process. As we show in Proposition 7, uy ~ U10 1]a for T - co. Therefore, we initialise Eq. (7) with uniform samples for suitably large time T (see Apdx. B.3) and numerically solve it to generate copula samples. We illustrate this procedure in Fig. and Alg. 3.
Optimal velocity predictor. As t) not available in closed form, our reflec- tion copula model learns to approximate the average velocity with a velocity pre- dictor ue (u) [0, 1 x [0, oc) + Rd. Analogous to flow models, for suffi- ciently rich model class, minimising the mean square error predictions against sample velocities of our systems recovers the expected velocity:
t=1.00
ts0.56
tw0.25
t-0,06
tw0 00.00
-v(,.t) t)
Figure 4: Reflection copula sampling. Initialised from uniform samples at t T, following - (u, preserves marginals and generates c(u) at t = 0.
0" = =ra min Eug Eu (ll vl
+ up- (u, 1) -[vr/u+ = =u]
e where the pair (ue, vx) follows the reflection process initialised uo Vo- We verify our models learn valid copulas in practice in Apdx. B.5 via statistical tests, rank diagnostics and calibration metries. Training details are given in Alg. 4 of Apdx.
5 RELATED WORK
Deep copulas. In Hofert et al. (2021), the dependence is modelled the Gaussian scale with moment matching network. This model is improved with a generative adversarial network design by Janke al. (2021). However, both methods are only able to output samples, precluding their use in density-reliant copula applications, and suffer from common mode collapse issues. In Kamthe al. (2021), a normalising flow approximates the copula obtain densities and samples. Theo- retically, their model also obtains valid copulas at optimality, but unlike our work, it lacks similar results to our Props. 2, 7, and on "forgetting and remembering" mechanisms to ensure the model evolves through valid copula states. Specific copula classes, such as Archimax (Ng et al., 2022) and Archimedean copulas (Ng et al., 2021), which focus on extremes, have also been estimated with deep networks. These works have model class which exactly matches given parametric copula family, which lends important robustness, especially for extrapolating to tail behaviours, but limits applications for the more general dependencies considered in work. Closest to Cde copula is the ratio copula of Huk et al. (2025), which we extend to the multinomial setting with diffusion- based dependence levels, motivated by the connection of Yadin al. (2024). Our reflection copula unique in leveraging principles of Holderrieth et al. (2024) in the context of a generative copula on [0, 1]d with preservation of marginal distributions. Our designs are the first instances of copulas scaling to complex and high-dimenssioal dependencies, such as in image data.
Copulas in diffusion works. To the best of our knowledge, our work is the first to utilise diffusion and flow principles for copula modelling. In Bibbona et al. (2016), serial dependence univariate diffusions is connected to bivariate copulas, which differs from our goal of flexible dependence models. In Liu et al. (2024), autoregressive language models are used as discrete copulas to improve sampling of text diffusion models. Finally, processes with Gaussian distributions studied in diffusion works (Kingma et al., 2021; Ou et al., 2025; Pierret & Galerne, 2025; Sahoo al., 2025), but no model knowledge leverages their univariate marginal preservation.
commonly
8

---

Published as a conference paper at ICLR 2026
6 EXPERIMENTS
In Sec. 6.1, we model challenging dependencies in synthetic data generation (Patki et al., 2016; Sun et al., 2019) and aulti-agent imitation learning (Wang et al., 2021) tasks. In Sec. 6.2, we show that our models are the first instance of copulas able to capture high-dimensional structured dependencies found in images. Appendix contains full implementation details, including computational times (B.1), analytical simulation study (B.2), showcase of the forward convergence rates (B.3), an ablation of our mixture loss (B.4), sample uniformity diagnostics with statistical tests and calibration metrics (B.5), well as sample quality visualisations (B.6). Code to reproduce results is publicly available at httpss//github.co cww.uuk.oooouuuuuoa Http//Huuu-Dii////e David/Diffusion serffusiio-nnd--aa and-aaddoaaaae www-bbasd-Copplla
Benchmarks and metrics. We compare to models whose bespoke goal to learn dependencies. We use the Gaussian and vine copulas as baselines, and the Implicit Generative Copula (IGC) of Janke et al. (2021) and the Ratio copula of Huk et al. (2025) as state of the art deep copula mod- els. Following (Nagler et al., 2017; Huk et al., 2025), we use the copula log- likeliilloood (LL) when available and use the Wasserstein-2 metric (W2) to assess samples. We compute the Frobenius norm (Frob) between samples and observations' Kendall's tau matrices (Hofert et al., 2018), measuring pair-wise dependence', and use the Frechet Inception Distance (FID) (Heusel et al., 2017) to assess image samples. The LL and Frob both purely measure the fit of the dependence structure, while W2 and FID also consider marginal uniformity of samples. All metrics are shown with standard deviations across 10 independent runs.
6.1 MULTIMODAL DEPENDENCE OF SCIENTIFIC DATASETS
Table Modelling dependence scientific datasets. Compared to existing copulas, the Cde model achieves better LL evaluations, while both the Cde and reflection copula obtain better samples.
Magie 10020, d 10i
13611,
16)
Robocup
133807,
20)
Frob
W2
Prob
W2
Frob
Causssin
3.92
1.76
0.27 10.06
20.00
137
0.10
1.20
3.50
0 0.15
Vine
6.50ga07
1.44
0.30
32.7010.14
1.35
0.95
1.80ga00
3.96...00
0.60ggo4
Ratio IGC
6.761aas
2.26
1.24
18.211000
2.5410a7
2.25
2.00joas
3.93
1.59
1.69
1.24
1.66
2.31
4.13
2.85ga12
Cde (curs)
18.654485
1.38 2.330000
0.21 20.05
50.21
1.1220.05
0.352000
3.40.0.7
3.87 os
0.51 0.51
Reflection (ours)
1.34_0.0a
0.2840.00
1.35 os
0.47
3.84
0.49
We choose datasets with complex dependencies between scientific variables for non-parametric cop- ula evaluation. Used in Nagler al. (2017): Janke al. (2021), Magic has telescope observation data, while Dry_Beat measures dry bean shapes. Following Wang et al. (2021), we model the behaviour of team of ten robots playing football by controlling their vertical and horizontal move- ments as probabilistic multivariate time-series (which is a common use of copulas, see Apdx. B for full details). In Tab. 1, our copulas achieve the best results for LL and W2, and perform comparably to the Gaussian copula for Frob. We report pair plots in Apdx. B.6 for visualisation.
6.2 HIGH-DDMENSIONAL STRUCTURED DEPENDENCIES OF IMAGES
To demonstrate methods' scalability, we require datasets with pre-estimated marginal distribu- tions. As this scale of experiments has not been fully studied in previous works, we rely on image datasets where CDFs are trivial and the complexity resides in the dependence structures; our goal is not to compete against image models. Furthermore, our choice of monochrome images provides an analogue to spatial problems, such as in climate research, where copulas are widely adopted.
Sampling and density estimation in high-dimensions. We follow Huk et al. (2025) and use the digits and MNIST datasets as high-dimensional examples with complex dependencies but simple marginals. We further use grey version of the Cifar dataset restricted to the two classes of cars and trucks. Notably, the dependence of copula data for MNIST is very specific for pixels in the middle of the image with noise the edges, while Ci Cirar requires the whole image to be dependent (see Fig. 5). In Tab. 2, the Cdc obtains the highest LL,, while for W2 our methods
GGasssian copulas advantaged on this metric due their closed-form relationship Kendall's

---

Published as a conference paper at ICLR 2026
Reflection
Cdc
IGC
Ratio
Vine
Gauss
Real
Data
3
5
5
9
Copula image samples: Data and copula scale MNIST samples (top rows), copula scale Cifar samples (bottom rows). Only our designs accurately represent the complex dependencies.
Figure 5:
Table 2: Copulas of high-dimensional data: Our methods outperform existing copulas capturing multimodal and high-dimensional dependencies according to LLL and sample-based metrics.
digita = 1797,
dd=6))
64)
GIMMA,
7867
ClEaz
I0244
Medell
11.
W2
FID
12.
W2
FID
1I.
W2
FID
daa.s
8.13
5.7
115.8ga3.4
35.50ams
102. 022,50-224
310.02200.
180.122244
Vine
11,20 10
8:20
6.06
10
36.30.006
48.2
NaM
84
100.04
Ratio
13.29
8.42
6.0420m
224.42
35.98_.05s
1348.18
49.91gm0.9
134.41
IGC
9.52
25. as
36.37 30.37
128.87
$2.09
209.08
Ca.
16.24
346
38
7.3810.19
1470.76.2490
Rellection
5.50g1.9
35 :0.000
13:0.00
2..00
42. 14.44
outperform benchmarks with the exception of the Cac for FID on digits and the reflection copula for W2 on Cifar. The vine failed to obtain LLs on Cifar regardless of hyperparameters. From Fig. 5, our Reflection copula samples are smoother, while the ede. has slightly noisier and grainier
samples, as reflected by the FID) values on di digits. and Ci far, see also Figs. 21, 22, 23. We believe this stems from the stochasticity of the sampling procedure as detailed in Alg. 1, which is beneficial feature maintain uniform marginals (see Apdx. B.5). Existing copulas do not scale well to these high dimensions and are unable to represent their dependencies.
CONCLUSION
In this work, we theoretically show that our processes only forget dependence, and that our models remember the true copula from them. Our work suggests new types of generative copulas using stochastic processes exactly targeting dependence, making copula models scale to complex depen- dencies in 1000 for the first time. Our copulas offer powerful options for copula likelihoods and sampling, empowering applications across the varied fields utilising such models.
Limitations and future work. While studied processes define valid copulas by design, our models are only guaranteed to represent copulas at optimality. An open question remains about how design architectures to preserve properties such uniformity of marginals and normalisation of the density for any model parameterisation. Relatedly, the question of how to sample with diffusion and flow models to best represent the model's true distribution is an active area of research which is equally relevant to our work. Further, diffusion models based on different marginal-preservve processes could better suit copulas with extreme dependence. The ingredients for that are CDF mapping to latent space with a diffusion whose invariant distribution is that CDF, which is applied independently across dimensions. Similarly, different velocity distributions for the Reflection copula could lead to better inductive biases. Finally, copulas require continuous marginal densities, an avenue for future research consists of extending our approach to discrete variables following the discrete copula notion of Geenens (2020).
10

---

Published as a conference paper at ICLR 2026
AUTHOR CONTRIBUTIONS
David Huk originated the idea, designed the project, wrote the code, ran the experiments, wrote the paper, and derived the proofs. Theodoros Damoulas supervised the project, provided valuable feedback, and helped edit the paper.
ACKNOWLEDGMENTS
We grateful for the helpful discussion and comments of the reviewers and Area Chair, well as those of Peter F. Holderrieth, Miha Brešar, Aleksandar Mijatovic, Thomas Nagler, and Ji Won Park. Furthermore, David Huk would like thank Unilink Software Ltd for supporting his attendance at the conference. David Huk funded by the Center for Doctoral Training in Mathematical Sci- ences Warwick. Theodoros Damoulas acknowledges support from UKRI Turing Al acceleration Fellowship [EP/V02678X/11.
REPRODUCIBILITY STATEMENT.
We take the following measures to ensure full reproducibility our work. We give all details of proofs and theoretical results in Appendix A with statements presented in the main text. We give full details of data preprocessing and model implementation in Apdx. B. We provide Algs. 1, 2, 3, and 4 with instructions to train and sample from our proposed models. Code reproduce experiments publicly available at wttpgggiwub.ooouuuuaaaaaai/a///e Ftw.D//iiiio------ooo---eeeeeeee
REFERENCES
Kjersti Aas, Claudia Czado, Arnoldo Frigessi, and Henrik Bakken. Pair-copula constructions of multiple dependence. Insurance: Mathematics and economics, 44(2):182-199 2009.
Iskander Azangulov, George Deligiannidis. and Judith Rousseau. Convergence of diffusion models under the manifold hypothesis in high-dimensions. arXiv preprint arXiv:2409. /8804, 2024.
Pietro Berkes, Frank Wood, and Jonathan Pillow. Characterizing neural dependencies with copula models. Advances in neural information processing systems, 21, 2008.
Enrico Bibbona, Laura Sacerdote, and Emiliano Torre. A copula-based method to build diffusion models with prescribed marginal and serial dependence. Methodology and computing in applied probability, 18:765-783, 2016.
Miha Brešar and Aleksandar Mijatović. Non Non-asy an-aayymptoicc bounds for forward processes in denoising diffusions: Ornstein-uhlenbeck hard to beat, arXiv preprint arXiv:2408. 13799, 2024.
Sitan Chen, Sinho Chewi, Jerry Li, Yuanzhi Li, Adil Salim, and Anru R Zhang. Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions. In International Conference on Learning Representations, 2023.
Pawel Chilinski and Ricardo Silva. Neural likelihoods via cumulative distribution functions. In Conference uncertainty in artificial intelligence, pp. 429. PMLR, 2020.
Kristy Choi, Madeline Liao, and Stefano Ermon. Featurized density ratio estimation. In Uncertainty in Artificial Intelligence, pp. 172-182. PMLR, 2021.
Adam Cobb and Brian Jalaian. Scaling hamiltonian monte carlo inference for bayesian neural networks with symmetric splitting. Uncertainty in Artificial Intelligence, 2021.
Gong-Gang Cong and Mark Brady. The interdependence between rainfall and temperature: copula analyses. The Scientific World Journal, 2012(1):405675, 2012.
Claudia Czado, Karoline Bax, Özge Sahin, Thomas Nagler, Aleksey Min, and Sandra Paterlini. Vine copula based dependence modeling in sustainable finance. The Journal of Finance and Data Science, 8:309-330, 2022.
RALPH D' "gggsstioo and Egon S Pearson. Tests for departure from normality. Biometrika, 60(3): 613-622, 1973.
11

---

Published as a conference paper at ICLR 2026
Philip Davies and Nicholas J Higham. Numerically stable generation of correlation matrices and their factors. BIT Numerical Mathematics, 40(4):640-6511 2000.
A Philip Dawid, Monica Musio, and Laura Ventura. Minimum scoring rule inference. Scandinavian Journal of Statistics, 43(1):123-138, 2016.
Paul R Dewick and Shuangzhe Liu. Copula modelling to analyse financial data. Journal of Risk and Financial Management, 15(3):104, 2022.
Gal Elidan. Copula Bayesian networks. Advances in neural information processing systems, 23, 2010.
Robin Evans and Vanessa Didelez. Parameterizing and simulating from causal models. Journal of the Royal Statistical Society Series B: Statistical Methodology, 86(3):535-568, 2024.
Leopold Fejér. Lebesguessche konstanten und divergente Fourierreihen. Journal für die reine und angewandte Muthematik, 1910.
Gery Geenens. Copula modeling for discrete random vectors. Dependence Modeling, 8(1):417-440, 2020.
Gery Geenens. Towards a universal representation statistical dependence. arXiv:2302.08151, 2023.
arXiv preprint
Tilmann Gneiting and Adrian E Raftery. Strictly proper scoring rules, prediction, and estimation. Journal of the American statistical Association, 102(477):359-378 2007.
Elisabeth Griesbauer, Claudia Czado, Arnoldo Frigessi, and Ingrid Hobæk Haff. Tvinesynth: A truncated c-vine copula generator of synthetic tabular data balance privacy and utility. In Proceedings The 28th International Conference on Artificial Intelligence and Statistics, volume 258 of Proceedings of Machine Learning Research, pp. 3511-3519. PMLR, 03-05 May 2025. URL https: http://pooceeeennnncce mlr press / v258 / grweerauere5a..... httm.....
Thomas M Hamill. Interpretation of rank histograms for verifying ensemble forecasts. Monthly Weather Review, 129(3)5550-560, 2001.
Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter. Gans trained by two time-scale update rule converge to local Nash equilibrium. Advances in neural information processing systems, 30, 2017.
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. Advances in neural information processing systems, 33:6840-6855 2020.
JL Hodges Jr. The significance probability of the smirnov two-sample test. Arkiv för matematik, (5):469-486, 1958.
Marius Hofert, Ivan Kojadinovie, Martin Mächler, and Jun Yan. Elements of copula modeling with R. Springer, 2018
Marius Hofert, Avinash Prasad, and Mu Zhu. Quasi-random sampling for multivariate distributions via generative neural networks. Journal of Computational and Graphical Statistics, 30(3):647- 670, 2021.
Peter Holderrieth, Yilun Xu, and Tommi Jaakkola. Hamiltonian score matching and generative flows. Advances in Neural Information Processing Systems, 37:1 10464-110493, 2024.
David Huk, Rilwan A Adewoyin, and Ritabrata Dutta. Probabilistic rainfall downscaling: Joint gen- eralized neural models with censored spatial gaussian copula. arXiv preprint arXiv:2308.0982. 2023.
David Huk, Yuanhe Zhang, Ritabrata Dutta, and Mark Steel. Quasi-Bayy Bayes meets vines. Advances in Neural Information Processing Systems, 37:40359-40333 2024.
12

---

Published as a conference paper at ICLR 2026
David Huk, Mark Steel, and Ritabrata Dutta. Your copula is classifier in disguise: classification- based copula density estimation. In Proceedings of The 28th International Conference on Arti- ficial Intelligence and Statistics, volume 258 of Proceedings of Machine Learning Research, pp. 3790-3798. PMLR, 03-05 May 2025. URL https: / httpes://rroceee /proceed ngg...g.. ..s.mma. press//2se/// huk25a html.
Tim Janke, Mohamed Ghanmi, and Florian Steinke. Implicit generative copulas. Advances in Neural Information Processing Systems, 34:26028-26039) 2021.
Sanket Kamthe. Samuel Assefa, and Mare Deisenroth. Copula flows for synthetic data generation. arXiv e-prints, pp. ar1iv-2101, 2021.
Diederik Kingma, Tim Salimans, Ben Poole, and Jonathan Ho. Variational diffusion models. Ad- vances neural information processing systems, 34:21696-21707, 2021.
Diederik P Kingma. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
Erik Kole, Kees Koedijk, and Marno Verbeek. Selecting copulas for risk management. Journal of Banking Finance, 31(8):2405-244 2423, 2007.
Chun Kai Ling, Fei Fang, and J Zico Kolter. Deep archimedean copulas. Advances in Neural Information Processing Systems, 33:1535-1545, 2020.
Yaron Lipman, Ricky T. Q. Chen, Heli Ben-Hamu, Maximilian Nickel, and Matthew Le. Flow matching for generative modeling. In The Eleventh International Conference on Learning Repre- sentations, 2023. URL httph::/ttpeeeeeeeeeee http///eepeeeeeeeee//eeeeeeeeeee wtwp:eiwwwweeeeetoooneeeeeee.dd
Anji Liu, Oliver Broadrick, Mathias Niepert, and Guy Van den Broeck. Discrete copula diffusion. arXiv preprint arXiv:2410.01949, 2024.
Xin Liu, Weijia Zhang, and Min-Ling Zhang. Hacsurv: hicrarchical copula-based approach for survival analysis with dependent competing risks. In International Conference Artificial Intelligence and Statistics, pp. 3079-3087. PMLR, 2025.
Olivia Michael, Oliver Obst, Falk Schmidsberger, and Frieder Stolzenburg. Robocupsimdata: A robocup soccer research dataset. arXiv preprint arXiv:1711.017000 2017.
Koichi Miyasawa al. An empirical Bayes estimator of the mean of normal population. Bull. 1961.
Inst. Internat. Statist,
Thomas Nagler and Claudia Czado. Evading the curse of dimensionality in nonparametric density estimation with simplified vine copulas. Journal of Multivariate Analysis, 151:69-89, 2016.
Thomas Nagler, Christian Schellhase, and Claudia Czado. Nonparametric estimation of simplified vine copula models: comparison of methods. Dependence Modeling, (11):99-120, 2017.
Yuting Ng. Ali Hasan, Khalil Elkhalil, and Vahid Tarokh. Generative archimedean copulas. In
Uncertainty in Artificial Intelligence, pp. 643 653. PMLR, 2021.
Yuting Ng, Ali Hasan, and Vahid Tarokh. Inference and sampling for archimax copulas. Advances in Neural Information Processing Systems, 35:1 10:117999-111 191191771,,,,,116111 2022.
Zijing Ou, Mingtian Zhang, Andi Zhang, Tim Z Xiao, Yingzhen Li, and David Barber. Improving probabilistic diffusion models with optimal diagonal covariance matching. In The Thirteenth International Conference on Learning Representations, 2025.
Arkapal Panda and Utpal Garain. Copula based trainable calibration error estimator of multi-label classification with label interdependencies. In Proceedings of The 28th International Confer- Artificial Intelligence and Statistics, volume 258 of Proceedings of Machine Learning Research, pp. 3745-3753. PMLR, 03-05 May 2025. URL https: : http///pcceee/geeeeeee/ee mlr.. press/v258 2a155eppana22555 html.
13

---

Published as a conference paper at ICLR 2026
Ji Won Park and Kyunghyun Cho. Semiparametric conformal predictioo. In Proceedings of The 28th International Conference on Artificial Intelligence and Statistics, volume 258 of Proceedings of Machine Learning Research, pp. 3880-3880 PMLR, 03-05 May 2025. URL https:// proceedi hgs mml. www.srresssssv. http///2288p2akkk2/ park25c..... html.
Ji Won Park, Natasa Tagasovska, Michael Maser, Stephen Ra, and Kyunghyun Cho. Botied: Multi- objective Bayesian optimization with tied multivariate ranks. International Conference on Machine Learning. pp. 39813-39833. PMLR, 2024.
Neha Patki, Roy Wedge, and Kalyan Veeramachaneni. The synthetic data vault. In 2016 IEEE international conference on data science and advanced analytics (DSAA), pp. 399-410, IEEE, 2016.
Carsten Peterson. A mean field theory learning algorithm for neural network. Complex systems, 1: 995-1019, 1987.
Emile Pierret and Brumo Galemn. Diffusion models for gaussian distributions: Exact solutions and wasserstein enrors. In Forty-second International Conference Machine Learning, 2025.
Prajit Ramachandran, Barret Zoph, and Quoc V Le. Searching for activation functions. arXiv preprint arXiv: arXiv:1710.0594) 2017.
Teodora Reu, Francisco Vargas, Anna Kerekes, and Michael M Bronstein. To smooth cloud or to pin it down: Expressiveness guarantees and insights on score matching in denoising diffusion models. In The 40th Conference on Uncertainty Artificial Intelligence, 2024.
Benjamin Rhodes, Kai Xu, and Michael U Gutmann. Telescoping density-ratio estimation. Ad- vances in neural information processing systems, 33:4905 4916, 2020.
Subham Sekhar Sahoo, Justin Deschenaux, Aaron Gokaslan, Guanghan Wang. Justin T Chiu, and Volodymyr Kuleshov. The diffusion duality. In Forty-second International Conference Ma- chine Learning, 2025. URL httph: opwnrevieeeeeeeeeeeeeeeeee .nt/ http://oommmiiix9ex0000
David Salinas, Michael Bohlke-Schneider, Laurent Callot, Roberto Medico, and Jan Gasthaus. High-dimensional multivariate forecasting with low-rank gaussian copula processes. Advances in neural information processing systems, 32, 2019.
Gianfausto Salvadori and Carlo De Michele. On the use of copulas in hydrology: theory and prac- tice. Journal of Hydrologic Engineering, 12(4):369-380, 2007.
Simo Särkkä and Arno Solin. Applied stochastic differential equations, volume 10. Cambridge University Press, 2019.
M Sklar. Fonctions de répartition à n dimensions et leurs marges. In Annales de I'ISUP, volume 8, pp. 229-231, 1959.
Jiaming Song, Chenlin Meng, and Stefano Ermon. Denoising diffusion implicit models. In Interna- fional Conference on Learning Representations, 2020.
Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. Advances in neural information processing systems, 32, 2019.
Yang Song. Jascha Sohl-Dickstein Diederik Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. In Interna- tional Conference on Learning Representations, 2021. URL
PorumiidPXIGGII wwt.dpxx2x2RRSS...... 2RRHS.
Akash Srivastava, Seungwook Han, Kai Xu, Benjamin Rhodes, and Michael U Gutmann. Esti- mating the density ratio between distributions with high discrepancy using multinomial logistic regression. Transactions on Machine Learning Research, 2023.
Sophia Huiwen Sun and Rose Yu. Copula conformal prediction for multi-step stepi time series prediction. In The Twelfth International Conference on Learning Representations, 2023.
14

---

Published as a conference paper at ICLR 2026
Yi Sun, Alfredo Cuesta-Infante, and Kalyan Veeramachaneni. Learning vine copula models for syn- thetic data generation. In Proceedings of the aaai conference on artificial intelligence, volume 33, pp. 5049-5057, 2019.
Natasa Tagasovska, Damien Ackerer, and Thibault Vatter. Copulas high-dimensional generative models: Vine copula autoencoders. Advances in neural information processing systems, 32, 2019.
Natasa Tagasovska, Firat Ozdemir, and Axel Brando. Retrospective uncertainties for deep models using vine copulas. In International Conference on Artificial Intelligence and Statistics, pp. 7528- 7539. PMLR, 2023.
Jom Van de Velde, Matthias Demuzere, Bernard De Bacts, and Niko Verhoest. Future multivariate weather generation by combining bartlett-lewis and vine copula models. Hydrological Sciences Journal, 68(1):1-15, 2023.
Vatter and T Pyvinecopulib Nagler. 0.6. 1, 2022.
Pietro Verzelli and Laura Sacerdote. study dependency features of spike trains through copulas. Biosystems, 184:104014, 2019.
Hongwei Wang, Lantao Yu, Zhangjie Cao, and Stefano Ermon. Multi-agent imitation learning with copulas. In Machine Learning and Knowledge Discovery in Databases. Research Track: European Conference, ECML PKDD 2021, Bilbao, Spain, September 13-17, 2021, Proceedings, Part 1 21, pp. 139-156. Springer, 2021.
Yunshu Wu, Yingtao Luo, Xianghao Kong, Vagelis Papalexakis, and Greg Ver Steeg. Your diffusion model is secretly noise classifier and benefits from contrastive training. Advances in Neural Information Processing Systems, 177:3230--33999 2024.
Shahar Yadin, Noam Elata, and Tomer Michaeli. Classification diffusion models: Re- vitalizing density ratio estimation.
In Advances in Neural Information Process-
ing Systems, volume 37, pp. 9837-9863. Curran Associates, Inc., 2024. https: / httpps://proceedingg neuri ca cccpaper_fiiiiccceceeee /ttp://appeeepe/22//////// 131888222200771aafcccbbllllle. Pappaee6eee---eer--eeeeeeeee-eee Con CeapperreCofeeeneeeecc
URL
A PROOFS
We restate the propositions and theorems for completeness, and provide their proofs afterwards.
A.1 PROPOSITION 2
For copula samples Zn on the Gaussian scale, consider uf ((t)) with zf defined through the process in Equation (2). Then, this process maintains uniform marginal distributions at times:
Proposition
-~~ U,,, W20,6C (1.....d).
Moreover, denote by at the copula of zt. Then, ct converges to the independence copula co in the KL divergence with rate O -.)..
The proof consists of converting results from stochastic differential equations on from the Gaussian scale to u on the copula scale.
Proof. To begin, rewrite Equation (2) in the following identity:
Z, =
+ 21 e~ (0,1).
For single dimension, note that z is standard Gaussian by assumption and cl is standard Gaussian by construction, and so z is also Gaussian. Its mean is thus zero and its variance is e 2t 1. Therefore, with O the univariate standard Gaussian CDF, u, (22)) is uniform.
+1-
For the second part, with the OU process as defines through Equation (2), under suitable assump- tions the data distribution p. the convergence the stationary (O,0 T) is well known with rate
15

---

Published as a conference paper at ICLR 2026
O 2) for some constant C. which depends the furthest mode of the initial data distribution Po (see e.g. Brešar & Mijatovié (2024) for the (DATA) assumption and Proposition 1.2). To convert this convergence from Z₁ the Gaussian scale u, the copula scale, notice that the mapping from one to the other is the standard Gaussian CDF applied dimension-wiss which is an invertible reparametrisation. The result is obtained via the KL's invariance under change of variable.
A.2 PROPOSITION 3
Proposition 3. For the classification task as defined in Equation (4), the true copula density c( (()) given value u E [0, 1] is equal to the following probability ratio:
Pt == 5-
c(u) = Pit Ta p (u))
where (1) is the inverse standard Gaussian CDF applied dimensions-wise to
Proof. We consider the classification problem of Equation (4), and in particular the class probabili- ties given by (P(1 = /2 2)) 2))................. P(I Ta =::=]]]] 2))). We define the mixture of data likelihoods between classes by p+
Σ -== T. Pr. (z). where the prior probabilities л, assumed iden- tical for all classes 8 E {1..... ......... k). By Bayes's rule, we show the following equivalence, similar to
identities used in Density Ratio Estimation works (Srivastava al., 2023):
pn (z)
pr. (z 2|6= Ti)
Pr,(2 = 2, = Ti)/
(z)
PT, (z 2/1 =T)
pT. (z 2,1 Ts)/ TT.
PT. (1 = Ti 2) (2) pr, (t = Ti z 2)
Pr, (1 = Ts/2 2)
pr, Th =
where the prior odds cancel out in the first line, as they are assumed to be identical. To conclude, note that copula is a ratio of densities shown in Section 3 of Huk et al. (2025), and such ratios are preserved under invertible functions by Lemma in Choi et (2021), meaning
( (w)) IT ( 10,1) Pr. (τ-1
(wC)
Pr. (w))
Pt Ti z +b (u))
cuu)
P(t = Talz - Φ-
where T-(()) the standard Gaussian CDF applied dimensions-we wise u.
A.3 PROPOSITION
Proposition 4. Consider the classification task as defined in Eq. (4), and let c, be the copula of u, = (22)) with & the standard Gaussian CDF applied element-wise, and 2, defined by the process of Eq. (2) at time E (Ti...., Tx). Then, the true copula score Va log c, (2) is given by:
Vn loglog. = w(u) 3 V= (P = Talz= =))===)))>> V log(P(( Talz =2)) (1))
where O is element-wise multiplication and w(u)= = ( (188((1)))"" (
Proof. We begin by showing a link between the class probabilities and the gradient of the joint density, which has been identified in existing works (Yadin et al., 2024; Wu et al., 2024), and then utilise this link within our diffusion copula context to derive an expression for the copula score. start, consider the classification task as defined in Equation (4) with class probabilities
(P(L = Ti = =))...., ...., ,PLL Ta/2 z))), and mixture of likelihood p+ (z) with weights s E {1...., k} as in the proof of Proposition 3. Rewrite the following two conditional probabilities:
pz = z|t T.) Pt
Tg zz
2) D+
(8)
T
and
P(t = Th Iz 2) p+
pT. (z= z|t Ti.
(9)
TTr
16

---

Published as a conference paper at ICLR 2026
Next, replace p+ (2) in Equation (8) with the appropriate quantity from Equation (9):
TT P( = = TT P(
p= T.)
PT, (z z|t TT)
Ws
where the class probabilities cancel by assuming equal prior odds. We take the log and then the gradient with respect to on both sides (where densities are positive and differentiable), yielding:
v Log(p, (z = z|t = T.))
(10)
v lo T2Iz =)) v log(I(( Tilz 2)) + V log(pn, (z 2|t Tx))
(11)
Here, we note that copula decomposition following Sklar's theorem leads to the following expres- sion of the log scores of a copula and joint density p with marginals p', E (l.....d}:
d
& log el
(²(24) V logp(z)
log (z')
In our case, for the joint density on the Gaussian scale (z 2|t T.). its marginal densities are standard Gaussian by construction as a consequence of the process in Equation (2). Thus, their marginal scores are =2i resulting in the following expression for the score of the log copula e, corresponding to the inter-variable dependence at time 8:
V x log c. (9(2)),
)x=)))
v log = =|t= TT) +
(12)
We can now replace the score of the joint log density in Equation (12) with the expression in Equa- tion (11), obtaining:
V log c. ( ((.,)) ()²n)) V log(P(( - T2|z - =)) V log(P(( = Ta z = =))
where v, log(Pn, (z = z|t =|t==100 Tk) cancels out with x as pT. is a standard Gaussian by construc- tion. Finally, we use change of variable from dimensions-wiss to u), obtaining the final expression:
a (2) (the standard Gaussian CDF applied
V log (u) w(u)=
V IP TzIz = 2)) V log(P(t Tal x))
2-1
where @ is element-wise multiplication and w(u) is defined
1
1
w(u)= = ( o ( )(b-L(ua))
A.4 THEOREM
Theorem 5. For the classification task defined in Equation (4), and a sufficiently broad family of multinomial models cac(z;0) =dd + {h E 10, 1 h' = 1} indexed by parameters 0, assume there exists set of parameters @ = {0 cac(2;0) = (IP( = Tilz = z)....,1PP IP(t Talz = for which the model is equal to the vector of class probabilities for all E Rd. Then,
for any weight a > 0, {0: 0 = arg ming Lea. (0)} = 0 for the loss:
Loge (8) =a Egnpr, [- -llog (a) (2:
+
E Ezn, upr, c>N(1..)
[lea(Caole To, + cT..-- - cc
where Cac (.) (z: 0) denotes the s-th component of the vector cac(2; (z; 0), and with
(cac((2;))) V1 ---.. (V log de (2; 0) Vz log ((s) (2; =).
(k)
17

---

Published as a conference paper at ICLR 2026
Proof. We first describe the set of minimisers for each component of the loss and draw conclu- sion afterwards.
The first term corresponds to the cross-entropy loss with prior probabilities per class:
(O) = Σ E. Bz-pr, [-log log Ce (2: (2:
Since cross-entropy is strictly proper scoring rule (Dawid et al., 2016), under the assumption that the model class contains the true distribution, minimising the loss with respect to the model parameters is equivalent to recovering the true data-generating process. Thus, we have OCr ==
arg ming (0) if and only if:
Cao(z; (GE) = = Ti|z - P((
= =2)
for all z over the support R4. Denote the set of all such minimisers as OcB.
The second term is the mean squared error between the predicted noise and the true added noise to a sample:
(MSE( = Σ Bzn ~rr; [elC_ Tx + V1 - (0, cll].
s=1 (0N().1)
This expression is minimised when, for all E (Ti...., Tk), and all Bd (except possibly regions of 0 mass):
Es(cao(2; z; 0)) Elelz x]:
for c and 2g such that 23 = c TT 2n + c 2.T. c. Then, by Tweedie's identity (Miyasawa
al., 1961), get
VI 1 (cde 0)) Vz log (z).
By using the definition of é (cde(z; 8)) and Equation (11) from the proof of Proposition 4, we obtain:
V log (5) (z; 0) V2 log Cd (z;0)
(k)
= log(P(( = Tg z = 2)) V log(I(( = Tklz = 2)) + Vz log(pn, (z = z|| Tk). =b))
Since PT, is standard Gaussian in d dimensions, the gradient of its log density is -2, leading
Cdc (z; 0) (k)
(x)
Pt = Tg zz = 2)
# log
Vx log
(13)
dc (z; 8)
p(t Ti iz 2) C.
de (2; 0)
(e)
P(t = T3 |z = z)
(14)
de (2:6)
P( Tkle z 2)
P(t = T,Iz = x) C. Cdd (z;0).
(8) de (z:0)
(k)
(15)
PP( Tyl
where C. R a constant, possibly different for each class s. Note that Cs by necessity. Furthermore, since cac 0) is normalised, we have
(2; da
(8) Cd (z;
+ 1 Co
(2;0)
sf1
(15)
do (2:0)
(k)
P(I TzIz = 2)
C. -. (2:
P(t ==
s
Cdc (2;0)
1 Σ
Ptt T,
C.}
*ith
x=)
18

---

Published as a conference paper at ICLR 2026
By substituting the last equality into Equation (15), obtain for any E fl,...,k)
P(t T,I2 = 2). C.
Cdc (z:
Pt = Ta +)+ Eith (P(t Tilz = 2) -c.
3 (5)
ddc (2;0)
(P(t = T, 2) C.
(P(t Ei=1
T = 2) Ci
Thus the set of minimisers for (M) can be expressed
PLL Tg 1z =) C.
{0= : de (z; 0)
= MSE.
Σ [P(1 = T: = 2)- CC]
To conclude, note that
(a) the set OCE is subset of OMSE. setting C₁ Vi in OME recovers OCE.
(b) parameters 8" for which cae (2;(2;;") (P(( Ti z ....... P(t Ta = 2)) the definition of both θ* and OcE.
By (a), and since a > 0, only solution of OcE can minimise Lew. (f). By (b), this equal 0r.
A.5 PROPOSITION 9, PROPOSITION 10, PROPOSITION 11
Here, we show the analogous properties from above for the case of the correlated OU process defined in Eq. 3.
Proposition 9 (Convergence to Gaussian copula). For copula samples 2o on the Gaussian scale, consider uf 2zz2) with zi as defined through the process in Equation (3). Then, this process maintains uniform marginal distributions at all times:
u ~ (0,1), 2 0,00 (I...a
Moreover, denote by c the copula of under Equation (3). Then, c converges to Gaussian copula with correlation matrix E as t + oo in the Kullback-Leible divergence with rate O 21).
The proof is similar to that of Proposition 2, and follows the same steps.
Proof. To show uniformity of u, at all times, decompose follows:
2, zo + VI- E~ N~~
Since the distribution of a single dimension 2 does not depend on the correlation matrix Σ, we can identify it standard Gaussian, meaning u, is marginally uniform.
For the convergence rate, similar existing results establish the limiting distribution of Equation (3) to be correlated Gaussian AO0, Σ) with rate OCCe- -2). The same reasoning in the previous proof applies, meaning we can derive the appropriate rate in KL divergence due to reparametrisation of of ((=))-
Proposition 10 (Cae copula density evaluation under Eq. (3)). For the classification task defined in Eq. (4) following the diffusion process of Eq. (3) with correlation E, the true copula density c(u) given value u E [0, 1]4 is equal to the following probability ratio:
c(uu (tt = T T (- N( (()-1(u)
P
N(6-1(u),Ia)
where th (1)( is the inverse standard Gaussian CDF applied dimensions-wise to
19

---

Published as a conference paper at ICLR 2026
Proof. The proof is identical to the steps in Proposition 3 with the exception that the terminal density N ④ (u), Σ). Thus the optimal classifier recovers the ratio:
P(4 Ti 2 l (w)) ýT, ( (u))
P( - Tx z = (u) N (u): Σ)
We obtain the copula density by multiplying the ratio by correction term
pr. + (u)) M( (w) du), N (2)
N (w)
which concludes the proof.
Proposition 11 (cde copula score-based sampling under Eq. (3)), Consider the classification task as defined in Eq. (4) for the process of Eq. (3), and let ey be the copula of us = ((2,)) with & the standard Gaussian CDF applied element-wise, and by defined by the process of Eq. (3) time
8 E (T.....,,Ty, T. Then, the true copula score V log c, (u) is given by:
V. log (u) = w(u) @ V log(((t T, 2)) V log(IP(t Talz 2)) + (Ia - (In-E=))
==
where @ element-wise multiplication and w(u) = ((
of (
The proof mirrors that of Proposition 4, with the exception that we need to account for the correlation matrix E and the different ratio target, in Proposition 10. In fact, when E Ig, exactly recover Proposition 4.
Proof. We define the mixture of likelihood (z) = Σx= T. p(2) with priori equal class
weights ns 8 E {1,...,kk kk. We can again rewrite the following two conditional probabilities as:
p. = 2/1 = =T,)=,)- PLL
Talz = 2) (z)
(16)
and
P(t = T z = 2) (+
PT. ( I=)=
(17)
TTr
Next, replace p+ (2) in Equation (16) with the appropriate quantity from Equation (17):
Th.
P(t=T.|z==)
p=
Pr. (z = x|t Tk)
Pi-=TiZ==
where the class probabilitics cancel due to equal prior odds. We take the log and then the gradient with respect to z both sides (where densities positive and differentiable), yielding:
V_ log(p, (z 2|t T.))
(18)
V= log(P(( = T2 = =)) - v= log(IP(t = Talz = 2)) + V log(pn, (z = 2|t = Tx))
(19)
Here, we note that a copula decomposition following Sklar's theorem leads to the following expres- sion of the log scores of copula and joint density p with marginals p', i E fl.....d].
d
Vz logg((=) Vx Σ log (x')
Vx log
In case, for the joint density the Gaussian scale pz =/t T,). its marginal densities are standard Gaussian by construction as a consequence of the process in Equation (3). Thus, their marginal scores are -2, resulting in the following expression for the score of the log copula c, corresponding to the inter-variable dependence at time
V. log c, (0:=)),, (())),((:))))
((=9))
V D2 2/1 T.) +
(20)
20

---

Published as a conference paper at ICLR 2026
We replace the score of the joint log density in Equation (20) with the expression in Equa- tion (19), obtaining:
V = log Cs ($(= ...;
-. {)=
Vx log(PP( T2 2)) = log(P(( T3 2)) + - -
Where we used that PT. - A0 Σ) by construction, so its gradient is -E-1 2. Finally, we use a change of variable from z to u ((2) (the standard Gaussian CDF applied dimensions-wise to u), obtaining the final expression:
Vn log Ga (u() = w v log(P(( = T, = -))- log(P(1 Talz - = ==0==))))+[4 +))+ 3-
where ( element-wise multiplication and vo(2) () defined as:
1
w(u)= = ((
o( ((b-((u))
A.6 PROPOSITION
Proposition 7. For any initial point uo E 10, 11t with vo (o, the reflected point Ur R (un +L-.v) vo) converges distribution to the independent copula on the hypercube as I - co. Further, if ug follows copula distribution, then u, also follows some copula distribution with uniform marginals for any time t O.
Proof. We first prove convergence the independence copula and then show that u, follows some copula ot at any time t > 0.
Convergence to the independence copula. As vo is independent across dimensions, we initially only focus the one one-dimnnii case of ul R₂ (uj + t vo), for {1...., d), with the reflection operator in one dimension on the u² component only being
x x] R1(x) = 1x (- [x]),
if [x] is even,
R.
if
| is odd
Notice that R1 is 2-periodic function. For bounded and continuous function f. have that g(x) == ((R1 is bounded, measurable and 2-periodic. Then, with ((2) as the one-dimensional standard Gaussian density, from definitions
[ (u/) Eif (Ri +t b))] BBg (né + + ++)]
glub +t 2) 6(2)d2.
Next, we use Fejér's Theorem:
Theorem 12 (Fejér (1910)). Let g: R be bounded, measurable, and k-periodic, and let ip be Lebesgue integrable. Then for any sequence on + to0 and any real constants an.
k (1." (() (())ddy)
(fC pl())z))
lim g( + an) (pxx) dr
In our case, since g is 2-periodic, and & is integrable, we have:
[ g(v) dy
glv) dy 2 L g(v) dy
1
lim slub + (2) ((=) dz
where we split the right-han hand side based on Ly] being even or odd, over [o, 1] and [1, 2] respectively. For the first part, using the definition of R1, identify the expectation of uniform variable
2 1. g(v) dy S/' f(v) f(y)dy dy 2 [B((U)) with U ~ 410
21

---

Published as a conference paper at ICLR 2026
For the second part,
achieve the same with substitution 2 y-
2 [ glu) dy / 2 [ f(2 dy = S. f(z)dz= -[[((U) with U ~ U[0, 1].
Thus,
lim Bl - with U u[0 1].
which shows that uf converges in distribution to the uniform density. Finally, to conclude, we note that since velocities vo are independent across dimensions, is the distribution of u,- The limiting distribution of u, is therefore the product of the limiting distributions of uz- This is the uniform distribution on [0, 1]" the independent copula.
Marginal uniformity at all times. Consider the marginal variable uf obtained from the reflection R1 (u: + 12.-v v2, vi) for any E {1....,)). d). Note that R2 is measurable as it is piecewise continu- operations, namely translations, and reflections
ous. It also combination of measure-preserve across the vertical axis, meaning RR is also measure-pressrring This is equivalent to having, for any bounded measurable f:
( + L. vo)) = = (()())))
By this reasoning, assuming Ug follows a copula distribution, every marginal remains uniform for t > 0, and so u, also defines a valid copula distribution at any
A.7 PROPOSITION
Proposition 8 Consider the reflection process introduced in Eq. (6) applied to copula samples ug c(u) with velocities vo ~ (O, La). Let u² (u,t) = B[vi = ul be the expected velocity at time O and location u [0, 11 and denote by or the copula of samples ur at time T. Then, starting from a value ur ~ Cr (ur) and following the probability path given by the ordinary differential equation
di =
backwards from T 0, the destination point comes from the initial copula distribution Un c(u).
Our proof largely follows Proposition in Holderrieth al. (2024), and include below deriva- tion for our reflection setup for completeness.
Proof. We first describe our reflection process ordinary differential equation (ODE) with initial distribution of samples Ug ~ c(u), and initial velocity distribution as Vo ~ (O Ia). As both components independent at time 0, our initial joint probability of being at location-veloce pair (uo: Uo) ~ II (u,, vo) = c(up) (v O, Ia). Therefore, the ODE governing the distribution (u,, v.) the interior (0, 1)d through time is defined as:
dt u v and d
d
Vt =,=0, 0,
di
with a Neumann boundary condition for reflection, where for öΩ the boundary of the hypercube, and nu)) the outward unit normal vector at a point u E An, we have for all u E 02 and 0:
foITA [T, v)do nuu) =0
(21)
We will work with the distribution of being location time t, given by c.(u) ;= l (u, v )dv which has been shown to correspond to a copula for all 2 0 Proposition 7. For our ODE, the change of the density [u, v) with t is described by the fokker-Plande equation:
Om, (u, v) = - u ((v, II, [1, v)).
Øt
22

---

Published as a conference paper at ICLR 2026
By integrating both sides in v, obtain
dc,(u)
/ Vu,o ((v,0) v) d
(22)
0t
de((u) vu (u)
V,0- 7700.1, (u, d)
(23)
Öt
(c((u) - v (c(()) (c(())-B((vlle = u|)
(24)
When U1 E 0S, boundary condition in Eq. (21) guarantees there is loss of mass at the boundary as [ = u| n(u) = O. Therefore, Eq. (24) proves that the marginal continuity equation for our process satisfied by o (u,t) =: E[ ur u].
B EXPERIMENTAL DETAILS
Pre-processing. As copulas have rarely been applied to the settings aim to explore, select datasets with simple marginal distributions but where the complexity resides in the dependence, the only exception being Robocup. For all datasets other than Robocup, we pre-process the data in Ied by applying the empirical CDF to each marginal, obtaining copula observations in [o, 1)d. For image datasets, 10 repeated runs, we split the copula observations into a train and test split according to Table 3, where splits for digits s and MNIST were chosen by following Huk et al. (2025). We also give the terminal times T at which the stochastic process has practically reached its stationary distribution, which we found be quickly attained experiments. We report more details on the time choices in Apdx. B.3.
add small Gaussian noise to dequantise it and obtain copula ranks. For
Table 3: Experimental setups:
For experiments from the main paper, we show the train and test percentages of the full dataset, as well the terminal times used for convergence of the forward process.
magic
Dry.Bean
Robocup
digits
MNIST
Cifar
Train/Tes: Split
80/20%
80/20%
80/20%
50/50%
50/50%
80/20%
Cde Terminal 7k 3
3
3
3
3
Ref. Terminal T 1.5
1.5
1.5
1.5
2
2.5
Pre-processing for Robocup. The Robocup dataset, first studied with copulas by Wang et al. (2021), consists of a 20-dimensional time series containing the vertical and horizontal positions of a team of 10 robots. The dataset is obtained from simulations of robot football matches Michael et al. (2017), specifically from the 25 games between the two teams cyrus2017 and helios2017, where we model the cyrus2017 team. To model the movements of the robot team, we adopt the following decomposition of the multivariate time series. For (xl, yy) R x R the horizontal and vertical movement of robot at time f, and (X,, Yd) E R. 10 x TR. 10 the vector of aggregated horizontal and vertical movements of all 10 robots, we model the next position
x; (X= V,-1, Xr-2, 22- + e"", yi (X,'' V X,-2, V2-2) + edd
where the mean functions f(x'), f(r) predict the expected next position based on the past two po- sitions of the team, with noise coming from s* s²',e"", e"", independent from the movements of the team. We model the mean functions network fo 40 + E. taking input the 40-dimensional position of the team over the last two steps, with one hidden layer of dimension 32, using the ReLU activation function. We use the Adam optimiser (Kingma, 2014) with 0.0001 learning rate, and optimise for 250000 steps, which takes 5 hours on CPU. We train separate model for each player's movement directions, resulting in 40 marginal mean models.
which we assume not predictable it is
fully connected
However, we assume that the noise vector ex with copula-based decomposition. We first estimate the marginal CDFs (using the empirical CDF)
is dependent, as such, modelling
23

---

Published as a conference paper at ICLR 2026
p=',p", of ex ey respectively, for i 1,.... 10, and use them to map the noise vector to the copula scale with u == (PM (s* Py (5t ).... (pr (E*)), py (e)"). We then model u with our copulas.
This approach is popular application of copulas dependent time series analysis, such as in finance and climate science (Czado et al., 2022). The mean functions serve to absorb the non-stationarity of the time-series, leaving i.i.d. noise which can be modelled with copula.
Further application of copulas. Another application could be synthetic data generation or data imputation. For instance, consider sensitive data such as criminal records about recidivism. The variables interact in complex ways with non-trivial dependencies, and obtaining more data is not feasible. Our copula models can then be used model the data first step, and in a second step, impute missing values of incomplete records based on this dependence. First, samples are generated from the fitted copula model (condional on the observed variables), and second, these generated copula samples are transformed via marginal distributions to the data scale, filling in the missing data.
Time discretisation. As in our theoretical derivations, for the Cae we pick equal priori class odds, meaning that during training we sample classes for times uniformly among the time discretisation. However, this discretisation does not need to be uniform in the time interval [7i, T&). For the Cdts on preliminary results from the training set, we found that Kullback-Leible discretisation worked best for the scientific datasets, while linear optimal transport discretisation (Lipman et al., 2023) worked best for image data. Following our result on the convergence rate Proposition 2, we define the KL discretisation as the time steps required for constant change in the KL. That is, for the total number of timesteps, the 8hh time step is given by:
t. - in(1- e-2x) T s 0, ......kk
The linear optimal transport discretisation is simply a uniform time grid from 0 to Ta with N points,
tg Tk. 0, k
For the Reflection copula, the following power-law discretisation across all experiments:
4,3 = Tk (
0.125
We provide an ablation for the choice of scheduler. For the Robocup experiment, we train our Cde model on 90% of the train set and use the last 10%for validation. In Fig. 6, we report the LL for three different choices schedulers, namely a linear optimal transport scheduler, a power-law scheduler and KL scheduler. While the scheduler does not affect the LL of the train set, we show that the KL scheduler obtains better generalisations and more stable performance on the test set, while the other two schedulers quickly overfit.
Log-Likellhood (Train)
Log-Likelihood (Test)
L1 PERT
KL Scheouler Mean = Std
KL Scheduler Mean = Std
Power law Scheduler Mean
Powerlaw Scheduler Mean
Linear Schoduler Mean + Std
Lincar Scheduler Mcan + Shd
Epach
Epoch
Figure 6: Scheduler ablation: For three different schedulers, we show the training and validation LL values through epochs. The KL scheduler is more stable and generalises better from the train to the test set.
24

---

Published as a conference paper at ICLR 2026
Practical implementation details. For all experiments, we first tune hyperparameters the train set, finding that the Cde converges quickly, while the reflection copula benefits from more training.
scientific datasets, we use a simple time-conditioned multi-layer perceptron architecture for the reflection copula, with 6 hidden layers of size 512 each. For the Cde on scientific data, we use a ResNet backbone with 6 hidden layers of size 512 each, followed by classifier head to output 50 class probabilities with a softmax activation to normalise them. We use the Swish activation function (Ramachandran al., 2017) for both models.
For experiments
For image datasets, for both models, we use the DDPM architecture from Ho et al. (2020), again with extra classifier head for the Cac- We use the Adam optimiser for training (Kingma, 2014) with learning rate 0.00005 for the cde and 0.0001 for the reflection copula. We train for 1K, 6K, 250, 250, 75K, 50K epochs for the cde- on datasets ordered as in Tab. The batch size is 1024 for scientific datasets and 128 for image data with the cdc. For the reflection copula, we train for 100K epochs on the scientific datasets, and for 500 epochs on image data, with a batch size of 512 for scientific data and 128 for images, iterating over the whole datasets at each epoch for images. We use mixture loss weights a 0.05, 0.05, 0.15 for the Cdte on scientific datasets, and use classes with a 0.005 for digits, and 256 with o = 0.005 for the other image datasets.
For the <de on digi ss and OU process on all other datasets. This is because their dependencies are concentrated on the diagonal when inspecting their pair plots.
we use correlated QU process, and we use an uncorrelated
For the implementation of benchmark methods, the ratio copula uses the same models as the Cdc except for the MNIST dataset, where better LL was obtained with a convolutional network model following Huk et al. (2025). We train the ratio copula for 3K, 25K, 3K, 1K, 2K, 5K epochs for datasets as ordered in Tab. 3. We sample the ratio copula using Hamiltonian Monte Carlo imple- mented in ham 1 torch by Cobb & Jalaian (2021), and tune the hyperparameters to have least 60% acceptance rates. For the IGC, in scientific datasets, we found that the same multi-layer perceptron architecture with 6 hidden layers of size 512 each led to frequent mode collapse. We therefore use a smaller network size with only 2 hidden layers of size 512 for scientific datasets, and use the full layers for image data. We train the IGC for 100K epochs. The Gaussian copula is fitted using the empirical covariance of the training data on the Gaussian scale. The vine cop- ula is fitted using the pyvinecopulin ib package of Vatter Nagler (2022) using non-parametric bivariate copulas.
Algorithms for the Cdce Here we give Alg. 1 for sampling, and Alg. for optimising the Cac copula. Note that like most copula methods (vines, Gaussian copulas, and IGC), sample on the Gaussian scale and transform samples back to the copula scale at the end. This slightly alters the form of the score from Prop. 4. We employ DDPM-style sampler for scientific datasets, and use a DDIM-style sampler (see Song et al. (2020)) for image data, as it resulted in less noisy outputs. Below, we present the DDPM-style sampler for the Cde. copula model, followed by the training algorithm.
Algorithmic changes for cae with correlated OU. Note that for the correlated OU process from Eq. 3, one has to compensate for the correlation matrix Σ. Concretely in Alg. 1, Line 2 becomes
ZTx ~ EO,E),
line 6 becomes
V ST, ((((Zz))) = Ven, log = Ty Tz(zn.) V xn, P(t = Tal [2))
line 7 becomes
OrT, . 2T, + (1 ar.) E V 2m G(
ZT
Var.
and finally line 8 becomes
ZT, + 2T, H VT- az with e ~ (N Ia) for H the lower triangular matrix from Cholesky decomposition of Σ.
For Alg. 2. the only changes are line 5 becoming
2, = 2T, + v 26-266 for e N( E),
25

---

Published as a conference paper at ICLR 2026
and line becoming
V1 2x. .20 (Vx, log P(t Thl=a) vz, log Ti|z,)) + 2z).
Algorithm Sampling algorithm for the ede copula model.
1: Initialize diffusion timesteps .
with class probability model cde (d)(2) (IP(t Til=)........ ..... P(t Ta 12).
2: Sample TT₃ ~ (0, Ia)
3: for t = A to do
Pre-compute constants for OU process: ar exp(2(
5: Compute copula using Prop. 3: a (((rr)) P(t = T, 21, (P = Ta 221
Compute copula score using Prop. on the Gaussian scale:
zT, log er (#(z))) V zm, log P(t 2T,) V log P(I Th (k(2TT)
DDPM update to Gaussian scale sample:
(aT, ZT, + (1 oT.)) Ven, log eT, ( (2T.))
2T.-1
VOT.
8: Add exploration noise: 2T,-1 + 2T,-, + c. vI -1-at at with c ~ (O)
end for
10: Map final variables back to copula space: uo ((2T,)
11: Return ug c(u)
Algorithm 2 Training the cde copula model.
1: Initialize model parameters number of epochs Nepoctas weight a > 0, timestep [Ti...... Tik
2: for epoch to NNpochs do
3: Sample Gaussian scale dependent data ZT, PT.
4: Sample diffusion time uniformly ~ Ui, 1.......
5: Compute perturbed data: 2, e - 2T, + /I e- 25 for NN Ia)
6: Compute class probabilities: Cde [=: f) = (P(t Ti Tilz.)...... Pt = Th [x|z,))
Compute cross-entropy loss for class probabilities: LcE = - log (t = s/2a)
8:
Compute score-based noise estimate:
' = /1 c 2s (V log P(t = Talo.) V.. log P(t = T 22,) + z.)
Compute mean squared error loss for added noise: LMSE =1(E -E(() 112
10: Take gradient step with respect to on combined loss: Less = a Lce + LMSE
11: end for
12:
Return trained model (2)(2;
Algorithms for the reflection copula. We now present Alg. to sample, and Alg. 4 to train reflection copula. Our sampling algorithm similar to an Euler Mallr--aruyymaa scheme for numerical integration. Even though we reflect samples even during sampling, in practice, the velocity predictor generally does not point outside the hypercube. Our training algorithm is simply a mean square error minimisation for velocity predictions.
26

---

Published as a conference paper at ICLR 2026
Algorithm 3 Sampling algorithm for the reflection copula model.
1: Initialize discretisation timesteps Ti..... Tk with velocity predictor v(u, 1).
2: Sample un, ~ [[ 1]"
3: for t = to do
4: Predict velocity at time t: êr, + u(UT, T:)
5: Update copula sample: UT.-1 ur, (Tt-1 Tτ) 6r;
6: Apply reflection to stay within [0, 1]a using Eq. 6, for i l...., d:
u T--1
-Lut
if LL is even else
u T
u
+ [1/
7: end for
Return ug c(u)
Algorithm 4 Training algorithm for the reflection copula model.
1: Initialize model parameters number of epochs Nepoctas terminal time T
2: for epoch = 1 to N Nepochs do
Sample copula data ur, c(u)
4: Sample velocity e ~ ~N I2)
T. - for U10, 1] 6: Compute reflected location and velocity using Eq. 5, for = ,,
Sample forward time uniformly
(u v R(un; + s v v
7: Predict velocity:
0. +- v(us 8)
8: Take gradient step with respect to on mean squared velocity error:
MMSE = | vell"
end for
10: Return trained reflection copula model u(u,t; 0)
B. COMPUTATIONAL TIMES FOR COPULA MODELS
Our copulas have different design goals; the Cae performs density estimation with single function evaluation, while the Reflection copula specialises sampling. While both designs use iterative sampling procedure (see Algs. 1, 3), each iteration of the reflection copula is a simple function evaluation of the network, while the Cae is function evaluation followed by gradient computation of the output with respect to the input. As such, the Reflection copula is preferred for applications heavily dependent on fast sampling, while the cde shines in LL-dependent applications. In Tab. NVIDIA
provide the training times and sampling times for producing 1000 samples, all done L.40S GPU except for the Vine copulas which used an AMD EPYC 9555P processor. The Cde uses number of sampling steps equal the number of classes, while the reflection copula uses 50 steps for all experiments. We do not provide times for the Gaussian copula as training simply taking the empirical covariance of the data the Gaussian scale, and sampling is just sample from that Gaussian with probability integral transform. This estimation simplicity comes at the cost low modelling capacity.
The vine copula training time scales mostly with data size and dimension, which explains the higher time of MNIST compared to Cifar. As vines inherently sequential and cannot take advantage
27

---

Published as a conference paper at ICLR 2026
of GPU compute, methods obtain faster sampling times compared vines. The IGC has fast training and the fastest sampling times, which are partially due to the smaller network size we used consequence of the frequent model collapse with larger network sizes. The Ratio copula has a faster training time than our methods, but requires much longer to produce samples due to the reliance on Hamiltonian Monte Carlo (differences in sampling times explained by differences in the HMC sampler on different datasets to yield 60% acceptance rate).
Table 4: Training and sampling times: For all methods, we measure the time to train and sample. Compared to the Cdc: the reflection copula trains for longer but samples much faster.
magic
Pry_Bean
Robocup
digits
MNIST
Cifar
Vine training
1.58
1.1s
13s
2s
2h
21min
IGC training
2s
3s
9s
11s
12min
14min
Ratio training
7s
62s
7s
56s
23min
37min
Cde. training
6min
7min
4 mim
100s
2h
4h
Reflection training
33min
40min
46min
90s
40m
Ih
Vine sampling
0.09s
0.17s
0.23s
0.77s
78s
131s
IGC sampling
1.2 10 10-5s
1.5 * 10
1.7 * 10 2
0.0008s
0.0013s
0.0015s
Ratio sampling
165s
327s
151s
554s
2h
2
sampling
0.07s
0.08s
0.08s
2s
42s
67s
Reflection sampling
0.01s
0.02s
0.02s
0.58
8s
15s
We additionally perform a dedicated computational time study in Fig. From our computa- tional study, we conclude that the network model is the most influential, followed by the number of timesteps used simulate the SDE/ODE when sampling.
In more detail, for d the data dimension. w the network width and / the number of hidden layers, Cde training has complexity of Old w + h- w²). This is because the model performs one forward pass (for class probabilities) and one backward pass (for scores) before taking gradients with respect to the loss during training. During inference, to sample the cdo, we perform a forward and backward pass for each of the & classes, following SDE discretization scheme, resulting in complexity of Okk (d wo + h- w²)). For the reflection copula, one training epoch consists of a single forward pass followed m network evaluations where m is the number of steps to simulate the ODE for, which scales as Om (d w + h w²)).
weight update, with complexity O(d h ) Generating samples requires
28

---

Published as a conference paper at ICLR 2026
CDC trais TIME Dmmrsssonn
COC simuason TIME Omarssaee
CUC LL eva Time Demertions
0.0040
CBC SIAAETA
CDC DEaEDES-
0.044
0.042
CULL SIMPSIR- 29
LEFU UNAO 200
0.09133
CSA teredejseble
0.0400
0.30333
016000
0.0039
00
(a) Cde training time
(b) Cde sampling time
(c) Cae LL eval time
Belleciu Tirs Demrersaoess
Belleciam siraaaionnn
Demeraiores
30
Redlctiion Eistedeppoii
6.0444
2.034
0.659
0.0000
Dirwniionn
Cimenalion
(d) Reflection training time
(e) Reflection sampling time
Figure 7: Computational time study: For varied dimensions and time steps, we report our models' training, sampling and LL eval (for the cao) times in seconds, reported as the mean over 25 runs with + one standard deviation.
B.2 SIMULATION STUDY
To show that our Cdte copula obtains correct density estimates, we study an example with an analytical density. Concretely, we use the following data generating process as the ground truth copula:
Copula of a mixture of Student's distributions. Let = 4 be the number of mixture com- ponents, with probabilities [T1, n2, T3, u4] = [0.3, 0.3, 0.2, 0.2]. each multivariate Student's T distribution with 10 degrees of freedom, with densities denoted by St For a mixture component, define the mean μ4
Fa = 4x with x~ x~N(;,,L.)
with correlation matrix Σx obtained using the Davies Higham Algorithm (Davies & Higham, 2000) based eigenvalues constructed by sampling 4 uniform variables in [0, 1]. and normalising the vector to sum to d. The final copula density at a value u = (F' (xx) (x)))........ F (x")) with F the univariate mixture of Student's T distribution functions, is obtained by the ratio of the joint mixture density divided by the product of mixture marginal densities:
Ei 5 NR, Ex)
c(u)=
d IH ( 5 10
(x:: (Ex)")}
where 5x. a univariate student's density with 10 degrees of freedom, 1u is the ith entry of pa and )2 is the th diagonal term of Mi. The resulting copulas are multimodal with shape not well captured by common parametric copulas, making them fit to assess our models' flexibility while still maintaining access to the ground truth density. They also display traits of extreme dependence due the heavy tailed nature of student's distributions.
Simulation study results. We assess cae model's LL values against those of the Ratio copula, in dimensions 10, 50, 100 with 8000 samples to train on and 2000 samples in the test set. We the same model architectures for the ede and ratio copulas for the scientific datasets, with the ede having time classes with a KL-based discretisation (see discussion on time discretisations above in Apdx. B), while the ratio copula only has 2.
In Tab. 5, we report the mean absolute error (MAE) and the mean squared error (MSE) LL values of aggregated over 10 independently and identically sampled datasets from fixed mixture of
29

---

Published as a conference paper at ICLR 2026
given dimension. Standard deviations over the 10 rums shown in subscripts. Our cdc obtains more accurate density estimates than the Ratio copula, and does so with less variability between runs too. This difference is accentuated in higher dimensions.
We visualise a single run from the 10-dimensional setup in Fig. 10, where we show violin plots displaying the densities of LL errors, with dots overlayed top represent samples for which the specific LL error was obtained. We colour-code these samples based on the ground truth copula LL,, allowing us discern patterns with respect the original copula LL. For instance, the ratio copula overestimates the LL of high-density samples under the original copula (blue samples tend to have c(u) > ratio(u)), and underestimates the LL of low-density samples (red samples tend to have c(u) < ratio(ur u)). Our Cde slightly displays a similar pattern with red dots being more frequent near the bottom of the violin plot, but the pattern less pronounced. Most of our model's errors con- centrate near 0, showing it correctly represents the dependence of the ground truth copula. Bivariate visualisations are displayed Fig. as aggregate views of all but two dimensions (showing just ut ) showcasing their multimodality.
Table Simulation study: For analytically tractable mixture copula densities in increasing dimen- sions, Cdc model obtains more accurate density estimates with lesser variance.
10
d = 50
= 100
Cde MAE
3.
13.7340.56
26.5641.8
Ratio MAE
3.40+4.22
34.18
65.87
Cde MSE
18.0341.23
291.43+2121
1077.30=100
Ratio MSE
19.30
2324.58
4808.75+1239.08
1.0
1.0
0.8
0.8
0.6
0.6
5
50
100
0.4
0.4
0.4
0.2
0.2
0.2
0.0
0.0
20
0.0
(a) d = 10 mixture copula
(b) d = 50 mixture copula
(c) d = 100 mixture copula
Figure 8: Simulation study: Plots of u u2 copula samples with colours representing LLs.
B.3 CONVERGENCE OF PROCESSES TO THEIR LIMITING DISTRIBUTION
As both methods use terminal time which the copula of the process has converged to indepen- dence, we show here a practical method for selecting such times. For observed copula data, we apply the forward process with respect to uniform distribution, which is the stationary distribution here. We also measure the W2 uniform distribution to itself using two sets of uniformly sampled data, which are shown in orange as a baseline. The reflection process is shown in blue while the OU process is shown in purple. In Fig. 9, we see that already for t 1, both processes obtain similar W2 values to uniform samples, indicating they suitably close to stationarity. This exemplifies the fast convergence rates derived in Section for the OU, and provides support for the reflection process as good choice of forward process.
and measure the Wasserstein-2 distance of the samples at different times
30

---

Published as a conference paper at ICLR 2026
Ref. =1 500, dev. (10 runsi
dev. nunst
Ret +1 std. dev. 110 runs:
sid, dee (10 nard
snd (10 nunal
1233
OU m1 sid thes ram))
Uriform
Liiffrn samples
semples
21 sta dev. (20 runst
Redlaction
Ref Retlection o' data
1315
12000
data
11.000
140
1250
24
25
2.4
(a) W2 magic 19020, d=10
(b) W2 Bean 1= 1=13611, d=16
(c) W2 robocup
II 135607, d=20
Ref. #1 std. dev. 110 runsi
Ref. 1 sti dev. (30 nns)
UU des
OU al atd dev. (20 rurs)
=1 dev. runs)
102 00/700400
0.82
Unifeem
oction
(1000 carmetee)
10.00
Unform samples
Roflect daea
=1 ctd. dev. (10 runst
12:02
12:8
Rufluct on of cata
=3 ctd dcv. (20 runsi
form samples
OU data
(il cate
12:88
CtI of cata
12:4
1240
250
25
10:18
Time
Time
Time
(d) W2 digits
(c) W2 MNI 6=0000, d=784
( W2 Ci far n= 10000, d=1024
n= 1797, d=64
Figure 9: Convergence to uniformity of forward copula process. Initialised at 1000 copula obser- vations, we run the processes forward time and measure the Wasserstein-2 distance with respect to the Uniform distribution. We show the result with one standard deviation in blue across 10 mea- surements, depicting the fast convergence to uniformity.
31

---

Published as a conference paper at ICLR 2026
True LL
LL true LL model
5
4
0
12
-5
-10
0
Left: Errors CDC, Right: Errors Ratio
(a) LL Error
True LL
12
10
0
4
MAE
8
6
12
0
0
1
Left: MAE CDC, Right: MAE Ratio
(b) Absolute Error
True LL
140
120
100
MSE
80
8
12
60
40
20
i
Left: MSE CDC, Right: MSE Ratio
(c) Squared LII. Error
Figure 10: Simulation study: We show the distribution of errors for our Cdo and the ratio copula. Overall, our model achieves more accurate density estimates, with emrors of smaller magnitudes.
32

---

Published as a conference paper at ICLR 2026
B.4 ABLATIONS ON THE MIXTURE LOSS
For our experiments, we use a mixture weight a to equalise both loss terms, following previous work Yadin et al. (2024). Here, we provide an ablation investigate the sensitivity of the Cdc to this choice. We perform a sweep over values = (0.0, 0.01, 0.25, 0.5, 0.75, 1.0, 5.0, 20.0, 500.0) for the experiments on scientific datasets. We report the resulting LL values on held-out validation set (10% of the train set) in Fig. 11. Our model is shown to be robust to values of a >> 0.01 with stable performance across different values. This indicates that incorporating the CE term is important to achieve optimal performance, which echoes our theoretical analysis in Thm. 5.
Magi
Dry Bean LL
Robocup
180
25-
65 80
X033
Epozht
-100
24
600
60%
0076-300
Figure 11: Ablation on mixture loss: We train with varied values of the mixture hyperparameter a and report LLs (# std. as dotted lines) on the validation set. Our method is robust to all values 0.01 with stable performance over hyperparameters.
B.5 UNIFORMITY OF GENERATED SAMPLES
For copula to correctly sample a joint density, its samples must be uniformly distributed for each dimension. This property is not obvious to maintain with deep learning models. As diagnostic tool following probabilistic calibration literature (Hamill, 2001), in Figs. 12, 13, 14, 15, 16, and 17, we show aggregate rank histograms over dimensions, split over 10 intervals of length 0.1 over the [0, 1] support'. For each subplot, the histogram comes from 10 independent runs on that dataset. As the Gaussian, Vine, and IGC copulas enforce uniform marginals samples by construction, we only show plots for our two models and the Ratio copula. In grey, we show truly uniform samples averaged over the same quantities. Therefore, deviations from the grey histogram indicate non-uniformity of samples. Deviations to the right compared to the uniform histogram show an over-representation of samples in that percentile range, while deviation the left shows an under-representation of samples in that percentile range. See Hamill (2001) for more details on rank histograms. The aggregate rank for run over bin by in (b., (b,,...,bio) {[0, 0.1]. ..., [0.9, 1]} is computed as:
model
N
Ranker
10
Iop (ulter))
N d
where (i...)) denotes the n-th sample in dimension from run r. d is the dimensionality. and N = 1000 is the sample size per run, with the indicator function [x)x) 1 if A taking the value ( else.
Overall, our two methods produce histograms similar to uniform samples, while the Ratio copula obtains more spread-out histograms. This means that the Ratio copula obtains biased samples each run, but the type of bias is different across runs. This issue is indicative of the difficulty of accurately sampling complex high-dimensional distributions with MCMC methods. All methods generally per- form worse on digits and Ci Cifar, possibly indicating the difficulty that those datasets represent for copula sampling task.
*Here, ranks refers the samples on copula space.
33

---

Published as a conference paper at ICLR 2026
samples weraget
Feeuery OfC sarplesu Serric aEES Wetegramn dreession
Requetcy tftr unbom samples
Frqqann Aatosunplo -vouge aotou fologiero ctriraate
hrsqurcy unlarr sumplex
Figure 12: Uniformity of samples for Magic:
rineodans
Ernpery oEE orrpleci sonage Hitsy wos per drwwdin
Requatry the arfim Omrol
Ants
cirendont
Figure 13: Uniformity of samples for Dry Bean.
Tezsshny 1et Satocaan autplo unsreges ouna talograto dicutto
hesusnoy clrs urbbr parpao
Daquundy of Auts euopleti vienged Nit.e Netecrone sec dirandont
Treos of unloen eueples.
Figure 14: Uniformity of samples for Robocup.
34

---

Published as a conference paper at ICLR 2026
swtrajel Neting
Reggeny OfC sarples aveage SIRSS Ietegnns diressions
Requercy ofiee unter sarples
Figure 15: Uniformity of samples for digits.
Seletioo diaesios
Pegpery Of carelss sorrged Nomgens direndins
Requerc owweerteeeeeeeee www.emmmmmmmmmmmm.mmmmm.mmm.mm.m CHISIEE
knds
Figure 16: Uniformity of samples for MNIST.
hezsenay of fatecian arrplo unaraged ucana faisg wna memkro
Teccanor toe uniber parpão
of
Deqquncy Bat eumples seeraged wwc.rnn cbreecont
Erqueccy off undeon curcles
Figure 17: Uniformity of samples for Cifar.
35

---

Published as a conference paper at ICLR 2026
We further investigate the uniformity of our samples with quantitative metrics in Tab. 6. For models to produce accurate copula samples, it is required that generated samples follow the respec- tive marginal distributions. For the Cdce the model operates on the Gaussian scale, we include an additional check for the normality of the generated samples via the statistical test of D d'agossino & Pearson (1973). For the reflection copula which operates the hypercube, mity of samples with a two-sample Kolmogorov-Sming -orinooro--mmooo test for goodness of fit (Hodges Jr, 1958). Along with the Ratio copula (operating on the Gaussian scale), we report the rejection rate of the hypothesis that marginals are preserved, averaged over dimensions and model runs, at the 5% alpha level with Bonferroni correction due multiple testing over runs and dimensions.
verify the unifor-
The results indicate that for most experiments, our models produce samples that do not significantly deviate from prescribed marginals and do so more effectively than the existing Ratio copula model. In particular, our Cat achieves a small rejection rate, which we attribute to the sampling procedure using Gaussian noise during sampling (see step 8 in Alg. 1). Note that this test only indicates departures from uniformity, but does not inform about their severity.
To strengthen these findings, we additionally report the Continuous Ranked Probability Score (CRPS) as univariate probabilistic calibration metric to compare our generated copula samples to uniform data, dimension-we wise. As a strictly proper scoring rule (Gneiting & Raftery, 2007), is minimised if and only if both samples come from the same distribution (Dawid al., 2016), here being the uniform on [0. 1]. In Tab. 6, we report the mean CRPS averaged across dimensions and runs, from which subtract the CRPS obtained with pairs of uniform samples (computed to be 0.1666 over 100 repetitions of 1000 vs 1000 draws). Thus, the ideal value for this metric would be uniquely achieved by uniform sample.
Our results corroborate the statistical tests, showing that our models achieve samples close to uni- formity, with a noticeable improvement over the Ratio copula model.
Table 6: Marginal uniformity metrics: For the Cat, Reflection, and Ratio copulas, for each dimen- sion, we report the rejection rate of statistical tests for the respective marginal distributions. We also report CRPS values to quantify deviations from uniformity. Our models generally achieve samples close to uniformity, with a noticeable improvement over the Ratio copula model.
Magic
Robocup
di digits
MNIST
Cifar
Cdee test rejection Reflection test rejection
0.01
0.01
0.0
0.03
0.00
0.00
0.02
0.09
0.00
0.16
0.01
0.70
Ratio test rejection
0.47
0.93
0.16
0.30
0.09
0.50
CRPS 0.0029
0.0021
0.0021
0.0068
0.0022
0.0041
Reflection CRPS
0.0024
0.0030
0.0022
0.0036
0.0024
0.0142
Ratio CRPS
0.0169
0.0128
0.0028
0.0150
0.0070
0.0680
B.6 VISUALISING SAMPLES
Visualising scientific datasets. We show pair plots for each dataset below, obtained from 1000 samples from our model, shown in the lower left, compared to 1000 randomly drawn observations from the test set (only one our 10 runs is shown). We note that pair plots are unable to convey higher-order dependencies; they only provide an aggregate view of the relationship between two variables. However, this lets us analyse the accuracy of our models in sampling from the correct regions. As can be seen in Fig. 20, both methods closely match the pair plots of observed data, displaying the same patterns across dimension pairs. We also display the accuracy of different methods targeting specific regions of the support in 18..
36

---

Published as a conference paper at ICLR 2026
1.0
0.3
0.6
0.4
0.2
0.0
0.0
0.0
(a) Observations of Magic.
(b) Reflection copula Magic.
Cae copula Magic.
Figure 18: Bivariate plot copula samples vs observations. We show the ability of our copula models' samples to match the target distribution. The reflection copula is particularly apt at capturing details.
Visualising high-dimensional image samples. We further show 25 samples per method for image datasets Figs. 21, 22, 23. From the samples shown on the copula scale Figs. and 23, note how the MNIST dataset consists of variation mostly in the centre of the image, while the edges mostly noise. Thus, models need to target specific dependencies in parts of the [0, 1| hypercube, while ignoring most of the remaining dimensions. In contrast, for Ci far, the whole image mean- ingfully contributes to the dependence, meaning models need to make the full 1024 dimensions fit specific dependence pattern generate image. We further inspect the dependence of images by downscaling the MNIST ST images to lower resolutions until equaling that of the digits data in Fig. 19.
8
16xi8
24 24
28
/
P.
2
2
3
3
3
5
7
7
8'
8
8
8
q
9
Figure 19: Dependence downscaled images. We downscale the MNTST data to the resolution of digits and observe the effects this has on the dependence via the copula scale. Dependence is fundamentally different, as the digits are no longer discernible on the copula scale and appear more like random noise.
37

---

Published as a conference paper at ICLR 2026
Dm.
Sin.
Dim.
Din.
Din, 6
mim. 10
Dm. 10
(a) Reflection copula Magic c.
(b) Cdc copula on Magí
(c) Reflection copula Dry_Er Dry_Bean.
(d) cde copula on Dry_Bean.
(c) Reflection copula on Robocup.
(f) Cdc copula on Robocup.
Figure 20: Pair plots of copula samples vs observations. In the lower left of each plot, we show our copula samples u with the reflection copula in blue the left column and our Cat copula in green in the right column. Observed data is shown in the upper right of each plot in grey. Note that instead of being mirrored, the dimension pairs are reversed for the bottom-left parts of the plots to match the orientation of the upper-right pairs. Best viewed digitally.
38

---

Published as a conference paper at ICLR 2026
Reflection
Cdc
1223
8
5
8
9
6
IGC
Ratio
Vine
Gauss
Figure 21: Samples for MNIST: Random samples from all methods, shown on the copula [0, 784 scale with lighter values being closer to 1. Notice how most of the variation takes place in the centre, where the shape of the number appears as high u values. However, image edges are not very dependent and take random values.
39

---

Published as a conference paper at ICLR 2026
Reflection
Cdc
6
5
4
1
3
1
0
7
8
55
3
9
5
7
5
O
6
8
9
2
I
8
6
8
9
0
6
8
2
7
8
7
9
6
9
9
8
9
3
IGC
Ratio
Vine
Gauss
Figure 22: Samples for MNIST: Random samples from all methods, shown on the data 784 scale with lighter values being higher. Our reflection and eae copulas correctly produce samples resem- bling digits, while competing copulas struggle to produce coherent samples.
40

---

Published as a conference paper at ICLR 2026
Reflection
Cdc
IGC
Ratio
Vine
Gauss
Figure 23: Samples for Ci far: Random samples from all methods, shown on the copula [0, 1] 1024 scale with lighter values being closer to 1. Notice how here, all dimensions of the image depen- dent and matter in producing coherent samples. As a consequence, the images are already visible on the copula scale without needing transform them the data scale.
41

---

Published as a conference paper at ICLR 2026
Sanity checks generated image samples. For both of our copulas, we verify that our models did not simply remember the training data but truly learned the data distribution. We qualitatively inspect the k 10 nearest neighbour (KNN) from the training set our samples, measured with the Fuclidean distance in IRd, reported in Fig 24. Our samples differ from the nearest observations, showing that the model learned the data distribution beyond pure memorisation of training examples.
Sm
COO
DEG2
DIS
CER 5
OBS
ONL
Vn
CBC
CBC
DISA
OH
COC.
0N6
OBS
0IS
DRG 10
5
S
S
5
5
S
5
5
5
S
5
9
4
9
9
9
4
7
9
9
7
7
Q
O
3
O
6
0
6
0
O
O
8
6
5
5
5
8
8
6
8
6
6
6
6
6
6
6
6
6
6
6
6
6 6
6
6
6
6
6
6 8
6
9
4
4
5
5
9
9
9
3
8
8
8
8
8
8
8
8
8
6
6
6
6
6
6
6
6
6
6
6
2
2
2
2
2
2
2
2
2
2
2
8
8
9
8
8
8
8
8
8
7
7
7
7
7
7
Reflection copula KNN MNIST.
(b) Cde copula KNN MNIST.
Cis
Olo
Ott
(c) Reflection copula KNN Cifar.
(d) Cdc copula KNN Cifar.
Figure 24: 10-nearest neighbours of generated samples: For random samples from our models, we show the 10 nearest neighbours with respect to the Euclidean distance. is apparent that the generated samples differ from the observations, meaning the models did not simply remember the training dataset. Best viewed digitally.
C THE USE OF LARGE LANGUAGE MODELS (LLMs)
In this paper, we only used LLMs to prototype data visualisations, to check parts of the code for er- rors, and to supplement the search for related works. All outputs of LLMs were verified by humans. All other aspects of this research were done by humans, specifically, originating the research idea, deriving proofs, producing diagrams and final figures for the paper, writing code for the models, conducting experiments, reviewing related works, and writing the paper.
42


Reviews:
[{"review_id": "ckrcfonRlK", "reviewer": "Reviewer_JPHs", "summary": "The paper proposes a principled framework that unifies copula theory with continuous-time generative modeling. It defines marginal-preserving forward processes (diffusions or flows) that progressively remove dependence while keeping all marginals uniform, so every intermediate distribution remains a valid copula. The core mechanism uses an Ornstein–Uhlenbeck diffusion in Gaussian space followed by the Gaussian CDF, which preserves uniform marginals and converges to the independence copula. Reverse dynamics are learned either as a score-based diffusion (likelihood) or a deterministic flow (fast sampling). Experiments report strong performance on high-dimensional data versus vine and neural-flow baselines. Conceptually, the work reframes dependence modeling as a diffusion/flow on the copula manifold, yielding a mathematically grounded and scalable approach to complex dependencies.", "questions": "Your Proposition 2 shows marginal preservation for the OU + CDF forward process. Does this extend to other forward processes (e.g., non-OU or state-dependent drifts/β-schedules)? If not, is OU intended as the recommended path?\n\nYour results report overall log-likelihoods and FID-type metrics, but these conflate marginal and dependence contributions. Since the defining advantage of a copula-based model lies in capturing joint dependence given fixed marginals, how do you disentangle and assess the quality of dependence fit separately from the quality of marginal fit?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Originality: Defining stochastic/deterministic dynamics directly on the copula manifold (uniform marginals at all times) is novel and technically \n\nQuality: The marginal-preserving construction (OU + CDF map) is well-motivated and, as stated, supported by a formal proposition; the dual reverse instantiations (score vs. flow) are complementary and well engineered. Empirically, results indicate improved likelihoods and scalability where vines or generic neural copulas struggle.\n\nClarity: Exposition and notation are clean.\n\nSignificance: Establishes a geometric generative foundation for high-dimensional copulas and clarifies a long-standing issue of marginal inconsistency in prior deep copula models.", "weaknesses": "In practice, sampling or training in copula space requires numerical integration on $[0,1]^d$ or its Gaussian transform.  The paper does not describe how boundaries $u_j \\in \\{0,1\\}$ are handled under discretization, nor how marginal uniformity is maintained when the OU process is approximated with finite steps.  Since the theoretical result assumes continuous time and ideal mapping, discretization could introduce marginal leakage.\n\nThe paper benchmarks primarily against classical vine and neural-flow baselines but omits newer methos such as [1].\nAlthough the paper claims linear scaling in dimension, no runtime, or memory/parameter analysis are presented compared to other methods and as the dimension increases (in addition to Table 4).\n\nThe weighting function governs how quickly dependence is forgotten. This choice may strongly influences the likelihood–sample-quality trade-off.  The paper fixes a single schedule without ablation, leaving unclear how robust the results are to this design.\n\n \nRefs:\n[1] Kamthe, Sanket, Samuel Assefa, and Marc Deisenroth. \"Copula flows for synthetic data generation.\" arXiv preprint arXiv:2101.00598 (2021).", "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": "jLo7z0b88X", "reviewer": "Reviewer_if98", "summary": "This paper introduces a novel framework for modeling multivariate dependencies using copulas, leveraging principles from diffusion and flow-based models. The core idea is to design forward processes that progressively \"forget\" inter-variable dependencies while provably preserving the marginal distributions. The authors propose two such processes: one based on an Ornstein-Uhlenbeck process on a Gaussian-transformed space, and another based on a reflection process on the copula hypercube.\n\nBased on these forward processes, the paper develops two distinct copula models. The first, the \"classification-diffusion copula\" (cdc), is designed for accurate density estimation. It learns the copula density by framing the problem as classifying the \"time\" or level of dependence in the forward process. The second, the \"reflection copula,\" is a generative model optimized for efficient sampling by learning to reverse the reflection process. A key achievement of this work is demonstrating that these models can effectively capture complex, multi-modal dependencies and scale to very high dimensions (d > 1000), a significant advancement over existing copula methods. Empirical results on scientific datasets and images show good performance.", "questions": "The reflection copula initializes velocities from an isotropic Gaussian distribution. Have you considered alternative or learnable initial velocity distributions? Could this provide a better inductive bias for certain types of dependencies and improve performance or training efficiency?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "1. The core idea of modeling dependence by learning to reverse a dependence-forgetting process is both elegant and powerful. It provides a new perspective on generative modeling that explicitly disentangles marginal behavior from joint dependence.\n2. The paper is built on solid theoretical ground. The authors provide proofs for the validity of their forward processes and demonstrate that their proposed models can recover the true copula at optimality. This rigor adds significant credibility to the methods.\n3. The models achieve impressive results, outperforming existing copula methods on complex scientific datasets.\n4. The most significant strength is the demonstrated ability to scale copula models to thousands of dimensions, as shown with the image datasets. This was previously considered infeasible for flexible copula models and represents a major leap forward for the field.\n3. The paper is clear, well-structured, and easy to read, which is commendable given the technical depth of the topic.", "weaknesses": "1. The loss function for the cdc model (Theorem 5) involves a hyperparameter α that balances the cross-entropy and score-matching terms. The paper states this is chosen to balance the magnitude of the terms, which is a common but heuristic approach. The paper would be stronger with a more detailed analysis of the sensitivity to this parameter or a more principled selection method.\n2. The paper correctly notes that the models are only guaranteed to be valid copulas (i.e., have uniform marginals) at optimality. While the appendix includes rank histograms as a diagnostic, a more quantitative analysis (e.g., using statistical tests for uniformity) of how close the generated samples are to being marginally uniform in practice would be beneficial.\n3. The FID score for the cdc model on the digits dataset is notably worse than other methods, including the simpler reflection copula. While the authors suggest the samples are \"noisier,\" a deeper explanation for this specific performance drop would be helpful. It raises questions about the robustness of the cdc approach on certain data types.", "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": "9sTm0Hw4ml", "reviewer": "Reviewer_dHoi", "summary": "The authors propose a type of flow mapping that maps a copula to the independence copula, i.e. the authors define a mapping that forgets dependencies. The model borrows techniques inspired by diffusion models to define a mapping that forgets dependencies such that an appropriate independent copula can be reversed and transformed to new samples that map to the observed data. The authors propose two variants: one that provides likelihoods and another that only provides samples. The authors illustrate the method on a number of high dimensional datasets to showcase the extent to which the method can be applied.", "questions": "While the experiments on image data are quite interesting to showcase the capabilities of the copula, can the authors describe some additional sets of data that would more amenable to this method? Most image data would be well suited for diffusion models etc, but the interesting components of the copula are specifying possibly different parametric families for the marginals and it would be nice to describe that within the paper. \n\nWhat is the computational cost of this method? \n\nHow often are the conditions of the copula empirically violated in practice?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The authors propose an interesting direction to enable high dimensional sapling from copulas, a difficult task that extends the applicability of copulas.\n\nThe method scales reliably to high dimensions which has been somewhat of an elusive property of copulas up to now.", "weaknesses": "The architectures prescribed do not preserve the necessary components to ensure a valid copula is being recovered, which can make some of these tools a bit difficult to employ in practice (e.g. if one needs specific families of marginals to be used).\n\nThe authors also miss some relevant works on deep learning and scaling copulas which should be taken into account, especially since these consider high dimensional questions using the stochastic representations [1,2,3].\n\n[1] Inference and Sampling for Archimax Copulas, NeurIPS 2022\n\n[2] Generative Archimedean Copulas, UAI 2022\n\n[3] Copula Flows for Synthetic Data Generation, arxiv:2101.00598", "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": "CoojtwIrk8", "reviewer": "Reviewer_5pKF", "summary": "This paper introduces a novel framework for copula modeling using diffusion and flow-based methods, focusing on the processes of \"forgetting\" and \"remembering\" dependencies within data. The authors propose two new models: the classification-diffusion copula and the reflection copula. The CDC is designed as an effective density estimator by framing the problem as a classification task to estimate the ratio of copula densities at different stages of a diffusion process. The Reflection Copula is a generative model built on a flow architecture that enables efficient sampling. Theoretically, both models are shown to recover the true underlying copula at optimality.", "questions": "No other issues; please refer to the Weaknesses section for specific points.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "**S1**. The core idea of explicitly modeling the \"forgetting\" (diffusion to independence) and \"remembering\" (density estimation/generation) of dependencies provides a fresh and theoretically grounded perspective on copula modeling.\n\n**S2**. The introduction of the classification-diffusion copula (leveraging classification for density estimation) and the reflection copula (a dedicated flow for generation) are sound and well-motivated solutions.\n\n**S3**. The most significant strength is the demonstrated ability to achieve SOTA results on complex, high-dimensional data, including image datasets. The models clearly outperform classical (e.g., Vine) and other deep copulas, which struggle with scalability and complexity. The paper also provides thorough empirical validation across diverse datasets, using multiple metrics and visualizations to support its claims.", "weaknesses": "**W1**. While the reflection copula samples quickly, the training times for both models, especially the CDC, are reported to be very long (e.g., up to 4 hours for CDC on CIFAR). This high computational cost could hinder practical application.\n\n**W2**. For image generation, the CDC samples are noted to be \"slightly noisier and grainier\" compared to the smoother samples from the reflection copula, indicating a potential trade-off between accurate density estimation and sample fidelity.\n\n**W3**. The paper could benefit from more detailed ablation studies to isolate the contribution of specific architectural choices or loss components to the overall performance.", "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)