Training Normalizing Flows from Dependent Data
Authors: Matthias Kirchler
Christoph Lippert
Marius Kloft
Venue: Submitted to ICLR 2023 (2023)
Keywords: Normalizing Flows

Abstract: Normalizing flows are powerful non-parametric statistical models that function as a hybrid between density estimators and generative models. Current learning algorithms for normalizing flows assume that data points are sampled independently, an assumption that is frequently violated in practice, which may lead to erroneous density estimation and data generation. We propose a likelihood objective of normalizing flows incorporating dependencies between the data points, for which we derive a flexible and efficient learning algorithm suitable for different dependency structures. We show that respecting dependencies between observations can improve empirical results on both synthetic and real-world data.



Under review as a conference paper at ICLR 2023
TRAINING NORMALIZING FLOWS FROM DEPENDENT DATA
Anonymous authors Paper under double-blind review
ABSTRACT
Normalizing flows are powerful non-parametric statistical models that function hybrid between density estimators and generative models. Current learning algorithms for normalizing flows assume that data points are sampled indepen- dently, an assumption that is frequently violated in practice, which may lead erroncous density estimation and data generation. We propose likelihood objec- tive of normalizing flows incorporating dependencies between the data points, for which we derive flexible and efficient learning algorithm suitable for different dependency structures. We show that respecting dependencies between observa- tions can improve empirical results on both synthetic and real-world data.
INTRODUCTION
Density estimation and generative modeling of complex distributions are fundamental problems in statisties and machine learning and significant in various application domains. Remarkably, normal- izing flows [Rezende & Mohamed, 2015) Papamakarios et al.) 2021 can solve both of these tasks the same time. Furthermore, their neural architecture allows them capture even very high- dimensional and complex structured data (such as images and time series). In contrast to other deep generative models such as variational autoencoders (VAEs), which only optimize lower bound on the likelihood objective, normalizing flows optimize the likelihood directly.
Previous work on both generative models and density estimation with deep learning assumes that data points are sampled independently from the underlying distribution. However, this modelling as- sumption is oftentimes heavily violated in practice. Figure Iillustaates why this can be problematic. A standard normalizing flow trained dependent data will misinterpret the sampling distortions in the training data as true signal (Figure Ic Our proposed method, on the other hand, can correct for the data dependencies and reconstruct the original density more faithfully (Figure ld]
The problem of correlated data is very common and occurs in many applications. Consider the ubiquitous task of image modeling. The Labeled Faces in the Wild (LFW, (Huang al.] (008 ) data set consists of facial images of celebrities, but some individuals in the data set are grossly overrepresented. For example, George W. Bush is depicted individuals in the data set only appear once. A generative model trained naively on these data will put considerably more probability mass on images similar to George W. Bush, compared to the less represented individuals. Arguably, most downstream tasks, such as image generation and outlier detection, would benefit from a model that is less biased towards these overrepresented individuals.
530 images, while around 70% of the
In the biomedical domain, large cohort studies involve participants that oftentimes are directly related (such as parents and children) or indirectly related (by sharing genetic material due shared ancestry)-a phenomenon called population stratification (Cardon & Palmer] 2003). These dependencies between individuals play major role the traditional analyses these data and re- quire sophisticated statistical treatment (Lippert al. 2011), but current deep-lcarning based non- parametric models lack the required methodology to do so. This can have considerable negative impact downstream tasks, we will show in experiments.
In finance, accurate density estimation and modeling of assets (e.g., stock market data) is essential for risk management and modern trading strategies. Data points often heavily correlated with one another, due to time, sector, or other relations. Traditionally, financial analysts often use copulas for the modeling of non-paramtrric data, which themselves can be interpreted as simplified version

---

Under review as a conference paper at ICLR 2023
(a) True distribution, sam- (b) True distribution, sam- pled independently
(c) Distribution learned by standard normalizing flow, trained on dependent data
(d) Distribution learned by normalizing flow adjusted for dependencies, trained on dependent data
pled with dependencies
Figure 1: Example setting with synthetic data sampled with inter-instance dependencies. Training standard normalizing flow this data biases the model. Adjusting for the dependencies during training recovers the true underlying distribution.
of normalizing flows (Papamakarios al., 2021).
Copulas commonly in use, however, are limited in their expressivity, which has led some authors even to blame the 2007-2008 global financial crisis on the use of inadequate copulas (Salmon, 2009). Many more examples appear in other settings, such data with geospatial dependencies, well in time series and video data.
In certain settings from classical parametric statistics, direct modeling of the dependencies max- imum likelihood models is analytically feasible. For linear and generalized linear models, depen- dencies are usually addressed either with random effects in linear mixed models (Jiang & Nguyen 2007) sometimes only by the inclusion of fixed-effects covariates (Price et al. (006). Recent work in deep learning introduced concepts from random effects linear models into deep learning for prediction tasks such regression and classification (Simehoni & Rosset 2021 Xiong al. 2019
Tran et al. 2020). In federated learning of generative models, researchers usually deal with the break of the non-ii.d. assumptions with ad hoc methods and without consideration of the statistical impli- cations (Augenstein et al. 2020 Rasouli et al. (020). These methods also only consider block-type, dependencies for multi- source integration. To the best of our knowledge, both deep generative models and deep density estimation so far lack the tools to address violations against the independence assumption in the general setting and well-founded statistical framework.
In this work we show that the likelihood objective of normalizing flows naturally allows for the explicit incorporation of data dependencies. We investigate several modes of modeling the depen- dency between data points, appropriate in different settings. We also propose efficient optimization procedures for this objective. We then apply our proposed method to three high-impas real-world settings. First, we model a set of complex biomedical phenotypes and show that adjusting for the genetic relatedness of individuals leads to a considerable increase in statistical testing power in a downstream genome- geneme-wide association analysis. Next, we consider two image data sets, one with facial images, the other from the biomedical domain, leading to less biased generative image mod- els. In the last application, we normalizing flows to better model the return correlations between financial assets. all experiments, we find that adjustment for dependencies can significantly im- prove the model fit of normalizing flows.
METHODS
In this section we describe our methodology for training normalizing flows from dependent data. First, we will derive a general formulation of the likelihood under weak assumptions on the depen- dencies among observations. Afterwards, we will investigate two common settings in more detail.
2.1 BACKGROUND: LIKELIHOOD WITH INDEPENDENT OBSERVATIONS
A normalizing flow is an invertible function t: R RD that maps p-dimensional noise variable u to a p-dimensional synthetic data variable x. The noise variable u is usually distributed following simple distribution (such as a (0, I,) Uj,,1r). ). for which the density is explicitly known and efficiently computable. By using the change of variable formula, the log-density can be explicitly computed as
log(p. lo)(pa(u) log(1 det (u)))
2

---

Under review as a conference paper at ICLR 2023
where u : t==t=(a0 (x) and Ju( is the Jacobian matrix of in
Given a data set 21..... Zn. if the observations are independent and identically distributed, the full log-likelihood function readily factorizes into its respective marginal densities:
(21, kn))
Σ (og(pa(z))) = E log(p.(u:)) log( det (u.)))
The function is usually chosen in such way that both the inverse t- and the determinant of the Jacobian J can be efficiently evaluated, e.g. using coupling layers (Dinh et al. 2017). Therefore, all of the terms in the likelihood can be explicitly and efficiently computed and the likelihood serves as the direct objective for optimization.
2.2 LIKELIHOOD WITH DEPENDENCIES
Assuming the data points are identically distributed. but not independently distributed. the joint density does not factorize anymore. A model trained on non-independent data but under indepen- dence assumptions will hence yield biased results for both density estimation and data generation.
We can derive the -oneinddeend setting as follows. Let T Ruxp - Rexp be the normalizing flow applied on all data points together, i.e.,
(
(x))
- =
Fn)
(2n)
X,U E RaXP are now matrix-variate random variables. We can still apply the change of variable transformation T, instead of the p p transformation t: log(px (X)) )) = log(pu (U)) log() det Jr (U))).
formula, but the
If ' is understood on np instead of space, it becomes clear that the Jacobian Jr is block-
diagonal matrix,
Je(w1)
0
Je(w2)
JT(U)=
Je Ju(n
for which the determinant is readily available: det (T) II Je(us). In other words, the log- abs-det term in the normalizing flow objective remains unchanged even under arbitrary dependence structure.
The density (U), however, is challenging and generally not tractable, and we will consider dif- ferent assumptions on the joint distribution of U.
In the most general case, we could assume that each is marginally distributed as Mp. i vari- able, with arbitrary dependence structure across observations. This direct extension of standard copulas to matrix-variate variables. As learning general copulas is extremely challenging even in rel- atively low dimensional settings (Jaworski al. 2010), we focus this work only on equivalent of Gaussian copula:
Assumption 2.1. We assume that the dependency within can be modeled by a matrix normal distribution MN with independent columns (within observations), but correlated rows (between observations):
U TMM, (0, C, Nmp (0, @ C).
Here, denotes the Kronecker product.
We can model the columns of with O-mean vector and Ip -poovaaia,, as the normalizing flow is usually chosen to be expressive enough to transform a aN, I,) into the desired data distribution. We note that this assumption means that we cannot model all forms latent dependencies so constitutes trade-off between expressivity and tractability.
Now we can state the full likelihood in the non-i.i.d. setting:
log(px (X)) =
log( det J+(ux))
np log(2.)) log(dct(C)) - 1r(U (C-)). (1)
2
2

---

Under review as a conference paper at ICLR 2023
2.3 SPECIFIC COVARIANCE STRUCTURES
We investigate different assumptions on the covariance structure in the latent dependency model. The most general case is fully unspecified covariance matrix, e-g. parametrized the lower- triangular Cholesky decomposition of its inverse, C = 1-- with n(n (n + 12/2 parameters. In this case, the determinant can be efficiently computed, as det(C) det(L ) = (e)t(L1122
(L=1(- Matrix products with C²11 can also be evaluated reasonably fast. However, this parametrization requires optimizing On²²) additional parameters, which is unlikely to yield useful estimates and may be prone to overfitting.
Instead, consider two different assumptions C that very common in practice and give reasonable trade-off between expressivity and statistical efficiency.
2.3.1 KNOWN AND FIXED COVARIANCE MATRIX
In many settings, side information can yield relationship information, given in the form of fixed relationship matrix G. The covariance matrix then becomes C' = >I, +1-A(1 with only parameter A E [0, to be determined.
This setting is commonly assumed for confounding correction in genetic association studies, where G is a genetic relationship matrix (where the entries are pairwise genetic relationships computed from allele frequencies (Lippert et al., 201 or based on pedigree information (e.g.. .. parent-child pair receives a relationship coefficient of 0.5 and a
pair of 0.25 (Visscher ett al. 2). Similarly, for time-related data, can define relationship via, c.g., a negative expo- nential function: C.j = (tp - )), where the hyperparameter y > 0 is a dime-decay factor and f2 and t are the measurement time points of observations and j. respectively.
More generally, G itself can again be a mixture of multiple relationships G - G,,, where G, denote multiple sources of relatedness. In this work, we consider G to be fully specified and only estimate
If the sample size is moderate (say, below 50k), efficient approach to optimizing À (Lippert et al. 2011) consists of first computing the spectral decomposition of = QAQ (with diagonal A and Then, the deg-determinant
orthogonal Q) and noticing that + +1-++G and the trace are
log(det(C))
Σ log(A + (1 (1-A))A,;), and
1=1
tr(U -=)
(AI, + (1 Q)A)-
The rotation matrix Q makes mini-batch estimation of the trace term inefficient, will either mix across batches or requires a full re-evaluation of Q(A,, + (1 2)A) after each update X, i.e., in every mini-batch. Instead, we optimize the parameters of the normalizing flow and in an alternating two-step procedure, see Section 2.4.2 Note that the main additional cost of this procedure, the spectral decomposition of G, is independent of and only needs be performed once for given relationship matrix G.
For larger sample sizes, there still exist practical algorithms for estimating the variance component (Loh et al. 2015). In practice, G is also often sparse can be approximated sparsely (e.g.. by setting all elements with absolute value below fixed threshold to 0). This can greatly accelerate parameter estimation and is usually accurate enough in practice (Tiang et al. 2019). More generally, different matrix structures may allow for additional speed-ups, but we defer this investigation future work.
2.3.2 BLOCK-DIAGONAL EQUICORRELATED COVARIANCE STRUCTURE
In the setting, consider block-diagonal covariance matrix C with equicorrelated correla- tion matrices C₁ E Rm, pn, as blocks:
C1
Ps
Ps
0
C2
where C: =
Pi
1
0
CN
Pi
1
4

---

Under review as a conference paper at ICLR 2023
with P: E 10 (we ignore the case of potentially anti-correlated blocks). In other words, there is no dependence between blocks, and there a constant dependence within blocks. We assume that the block structure is known ahead and we only need to find the parameters Pą- For each block there either no (n: = 1) or only one (n; > 1) parameter to be learned.
The assumption of equicorrelated blocks is reasonable in settings with repeat measurements of iden- tical objects or individuals. E.g., in a facial image data set, certain individuals may have multiple images. This setting similar to the setting of high-cardinality categorical features in prediction models (Simchoni & Rosset (02T).
Both the determinant and the inverse of each block can be efficiently computed (Tong 2012) 2). Prop- 5.2.1 & 5.2.3):
1+ I+ - 21p (1-p))++(m= 1)p)
if = k
det = + (n, 1)p;)( -1a))m- and (C,
-P: (1- (L-
otherwise.
2.4
OPTIMIZATION
2.4.1
MINI-BATCH ESTIMATION
The full likelihood in Equation can be computed explicitly but does not lend itself easily to stochas- tie optimization with mini batche. Note that the log-abs-det term decomposes nicely into indepen- dent observations and the next two terms independent of the observations. Only the trace term is problematic for mini-batch estimation, so we propose an unbiased stochastic estimator for it.
Proposition 2.2. Given mini-bateh of size b 2 and & {0,1)* th variable indicating batch inclusion (i.e., Ti is in batch iff E 1; 1==1 & = b) and A == C--1, the stochastic trace estimator >x=
Σ au us + n(n 1)
Σ
trg =
(2)
b(b 1) i<j
is unbiased, i.e. Eg(frg] Uu AU).
The proof can be found in Appendix A The trace estimator frę only depends on observations u, = (x;) within the batch and can be efficiently computed, assuming A = C-1l can be efficiently evaluated, which is the case for the parametrizations discussed in Section
2.4.2 TRAINING SCHEDULES
From here on, we distinguish between the true parameters À and ps, and the parameters estimated by model, S, Ps.
Known & Fixed Covariance Joint optimization between Ā and the parameters of the flow is possible, but would require in each step a full re-evaluation Ir(UT C-10) across the full data set, instead of just the current mini-batch. This makes this training scheme infeasible. Instead, we propose two different methods to optimize both the flow parameters and variance component Å. First, we can use a simple grid search over different possible values for Å and choose the best according to performance
validation set.
Second, we can use an alternating descent approach. In this case, we alternate between optimiz- ing only the parameters of the flow model for a number of epochs (with version of mini-batch stochastic gradient descent) and only optimizing Ả for number of epochs (with gradient descent). At the beginning of every flow-parameter training stage, we compute the current = C-1 for the given and can then compute all mini-batch likelihood estimates using the trace estimator in Equation
2 without the need for recomputation. At the beginning of every À training stage, we only once compute the rotated noise variables QTU for the full data set and can then optimize the derivative of the full objective with respect to very efficiently. The trace can be computed as tr((QT)) + (1 cyclical trace property, but in our experiments found that this was not bottleneck computation.
or as tr(Q (QQUUQQUU) (SI, + (1 \)A) due to the
5

---

Under review as a conference paper at ICLR 2023
Table 1: Results in terms of the test-data negative log-likelihoods for synthetic data with equicorre- lated blocks (top) and fixed covariance (bottom), averaged over 10 random seeds (lower better). Significantly better results are bold (one-sided paired t-test, a 0.05). Baseline the same model without taking dependencies into consideration.
Algorithm
Abs
Crescent
CrescentCubed
Sign
SineWave
Equicorrelated
Baseline
1.513
2.021
3,010
1.519
2.070
Blocks
Grid Search
1.379
1.885
2.938
1.420
1.983
Joint
1.475
2.005
3.067
1.501
2.087
Fixed
Baseline 1.590
2.144
3,262
1.818
2.353
Covariance
Grid Search
1.376
1.866
2.944
1.443
2.040
Alternating 1.361
1.869
2.915
1.430
2.064
To yield values in the interval [0, 1], we chose to parametrize Ả as the output of sigmoid function Å of (Araw), where Xrow E the raw optimization parameter. We tried different sigmoidal parametrizations, but those had little effect on the outcome.
Equicorrelated Blocks In the case of equicorrelated blocks, we also propose two different training schemes. First, we can again use a simple grid search over a single joint parameter ý -= وI = ÅN- Alternatively, if there only very few blocks, grid search for all P. is possible, although the exploration space grows exponentially with the number of blocks N.
Second, due to the simple computations of det(C)) and C-1 in this case, can also perform joint optimization over the flow parameters and all P.. We again parametrize pas with raw parameters pushed through a sigmoid function as for
EXPERIMENTAL EVALUATION
We validate on both synthetic and real-world data that our novel training scheme can help alle- viate sampling biases when training normalizing flows. On real-world data with non-independent data, the ground-truth dependency structure is usually not known, making the evaluation inherently challenging. Therefore, we first investigate simulated settings where we can explicitly control the dependencies. Our evaluation metric in all settings is the negative bog-likelihood (NLL) on a holdout test set. For the imaging experiments, we also report bits per dimension (bpd), linear transforma- tion of the NLL. Additional details for all experiments can be found in Appendix B
3.1 SYNTHETIC DATA EXPERIMENTS
3.1.1
EQUICORRELATED DATA
In the first setting, we simulate draw with repeat measurements, inducing an equicorrelated depen- dency structure described in Section 2.3.2 For each block, draw one p: Unifjo.5.0.9) and define the full covariance matrix in Section 2.3.2 Using this covariance matrix, sample from several non-parametric
distributions provided by Durkan et al. (2019). data set can be seen in Figure For modelling the equicorrelated blocks, we choose both grid search over fixed parameters and joint gradient-based optimization of D and the flow parameters.
example for the Abs
Table (top part) shows the result. Surprisingly. while the grid search clearly outperforms the baseline, the joint optimization does not improve upon the model. For the Crescent data set we also computed the distance of learned A, to true p,s for the best models (each block is counted only once, independent size). The baseline model had an average MSE 0.57 (MAE: 0.74), while grid search and joint optimization had MSEs of only 0.023 (MAE: 0.13) and 0.08 (MAE: 0.25), respectively. In an additional experiment, this time only on the Crescent data set, we investigate how sensitive only change the sampling of true dependency parameters from a Unif, distribution, with inter- val E { 0.2), [0.2, 0.4], [0.4, 0.6), [0.6, 0.8), [0.8, 1.0]). The results are shown in Figure 2
model is the strength of dependencies. In the data creation,

---

Under review as a conference paper at ICLR 2023
Table Results on real-world data, negative log-likelihoods on test data set, averaged over 10 ran- dom seeds for UKB & Stock Pair data (lower = better). P-values for one-sided paired t-test against baseline in parentheses. Baseline is same model without taking dependencies into consideration.
Algorithm
UKB Biomarkers
Stock Pairs
ADNI (bpd)
LFW (bpd)
Baseline
24.50
-5.69
7836.5 (2.760)
6481.6 (3.044)
Grid Search
24.27 (p 0.002)
-5.72 (p 0.002)
7725.8 (2.721)
6334.7 (2.975)
Joint
-5.71 (p 0.003)
7715.8 (2.718)
6346.5 (2.980)
Alternating
24.04 (p 0.000083
At each of the five data set settings, one-sided paired t-test shows that the normalizing flow incorporating de- pendencies outperforms the baseline (a = 0.05). As ex- pected, for small dependencies in the true data, both mod- els perform similarly, but our method very robust and barely decreases in performance up until the highest range of sampling distortions.
Baseline
2,10
Grid Search
2.05
2.00
1.95
3.1.2 KNOWN COVARIANCE
1.90
We next simulate the setting of known covariance ma- trix between different samples but with unknown variance component A. We use the covariance structure AI + (1 XJG as the equicorrelated case to generate correlated bi- variate standard normal samples that again get non-linearly transformed. In Table [aaleT[borr (bottom part) we compare the re- sults. Both the simple grid search and the alternating de- scent approach perform considerably better than the naive baseline algorithm that ignores the dependencies in the data.
1.85
[0.0.2]
[0.2.04]
[0.4.0.6 [0.6,0.8 [0.8.]] p range
Performance of baseline model versus model adjusted for de- pendencies, for different strengths of dependencies (p).
Figure 2:
3.2
REAL-WORLD DATA
3.2. 1
UKB BIOMARKERS
The UK Biobank (UKB, (Bycroft al. 2018) provides rich phenotyping and genotyping for large cross-section of the UK population. We investigate number of blood biomarkers, whose distribution starkly deviates from standard parametric distribution families. Usually, the data needs to be transore-r to match normal distribution (Monti et al., 2022), which, however, can decrease power of many statistical methods (MeCaw et al.] 2020). These biomarkers are well-known to be highly heritable and subject to population stratification, type of confounding due to joint ancestry of unrelated individuals Siinott-AAmsteee et al., (021). We perform two experiments on this data set, building non-paue aaraammettric density model that can incorporate the distorting genetic correlation between individuals.
Density modeling In the first experiment, select the 3,223 individuals for whom all 30 biomarkers are available. Relatedness between two individuals is computed as the correlation co- efficients between the individuals' first 40 (unnormalized) genetic principal components (computed from SNP microarry chip data provided by the UKB resource). We use this Matrix as the fixed covariance structure and optimize over Å. This way of measurement of genetic relatedness between individuals is very common in genetic association studies and has been shown to reliably correct for population stratification (Price et al.) 20066. We investigate the density estimation on the test data. Due to the relatively small data set size, we re-run the same experiment times with differ- random splits between train, validation, and test set and also different network initializations. The results in Table 22 first column, indicate that incorporating the dependencies can significantly improve model fit, both using grid search and using the alternating optimization scheme.
7

---

Under review as a conference paper at ICLR 2023
Application in association studies A genome-wee wide association study (GWAS) is frequentist hy- pothesis testing procedure, in which a phenotype is tested for association against a large number of individual genetic mutations (typically on the order of hundreds of thousands or millions of variants). GWAS are a fundamental tool within multiple disciplines in the life sciences, such as in the medical domain, in plant and in lifestock breeding, and have considerably contributed to the understanding of the genetic architecture of complex traits (Visscher et al. (017). State-of-the-art GWAS algo- rithms model dependencies between individuals with random effects in a linear mixed model (LMM) framework and can effectively control for both population stratification and (known and cryptic) re- latedness between individuals (Yu et al. (006). In this experiment, we perform multivariate GWAS, testing for association between individual genetic variants with multiple phenotypes together. Due to the high computational cost of multivariate LMMs (mvLLMMs), we split the 30 available biomarkers into six disjoint groups of related biomark- ers and subsample 10,000 individuals per group. Rank-based normal trans- formations are insufficient to transform vector of arbitrarily distributed ran- dom variables into multivariate normal distribution, as would be necessary for mvLMMs. This is due to the fact that not all random vectors whose marginals are normally distributed are also mul- tivariate normally distributed; see Fig- frr illustration on the biomarker data. Hence, mvLMMs can not be ap- plied to muantile-transform stead, a standard method is to test for association with each biomarker in the group independently, take the minimum of the p-values over all biomarkers the group, and perform a (Boneeroni-correetinn for the number of biomarkers in the group. We propose to instead use normalizing flow to transform the biomarker group into a multivariate normal vector and then apply the mvLMM on this transformed data. We use both a baseline normalizing flow without consideration of the data dependencies, and our proposed method with the alternating optimization scheme. More details can be found Appendix B.2.
Table 3: Number of loci associated with biomarker groups at genome wide significance level, avg. over random seeds. Single: univariate, crautile-transformee LMMs; Baseline: mvLMM on flow-transformed biomarkers; Alternating: mvLMM on biomarkers trans- formed with flow correcting for dependencies.
Biomarker group (# biomarkers)
Single
Baseline
Alternating
Bone and joint (4)
18.7
12.0
16.3
Cardiovascular (8)
55.0
58.0
61.0
Diabetes (2)
5.3
5.3
6.7
Hormonal (4)
6.7
5.7
7.0
Liver (6)
29.7
32.3
35.0
data. In-
Renal (6)
18.3
19.0
18.7
All
133.7
132.3
144.7
We report the number of indepedent loci significantly associated with each group of biomarkers in Table 3 While the baseline normalizing flow performs similarly to the naive single-dimension approach, taking care of the dependencies can boost the number of found loci by more than 8%.
We believe these findings may also significantly improve the analysis of more complex intermediate phenotypes, such as in full-imaging GWA studies (Kirchler et al. (022). To the best of our knowl- edge, this is the first time that normalizing flows have been used for GWAS in this style, although Hansen et al. 1 recently used normalizing flows in different GWAS setting.
3.2.2 IMAGE MODELING
Image modeling is major research area for normalizing flows, with applications in image synthesis (Kingma & Dhariwal 2018), outlier detection (Schirmeister et al., (020)), and semi-supervised learning (Izmailov et al. 2020p. Repeat measurements are very common in image data sets, and without adjusting for dependencies, overrepresentation biases will translate into biased generative models, as well. We investigate two prominent examples.
ADNI brain imaging The Alzheimer's Disease Neuroimaging Initiative (ADNI, (Jack Jr et al.) 2008 is longitudinal study of Alzheimer's Disease (AD) progression, so many of the individuals in the study are imaged multiple times. Prior work on similar data has shown that causal effects can be modeled in generative image models using explicit confounding factors such as age and sex (Pawlowski et al., 2020). Here we show that we can also model the i...d-violations using proposed method. The data set comprises 1,820 individuals with each individual having between 1 and 35 images (mean: 7.03, median: 6) and a total of 12,799 images. We model these repeat
8

---

Under review as a conference paper at ICLR 2023
measurements with the equicorrelated model and use Glow-type image normalizing flow Kingma & Dhariwal 2018) as our base architecture.
LFW face images LFW (Huang al. 2008 consists of 13,233 facial images of 5,749 celebrities, where each individual has between and 530 images (mean: 2.3, median: 1). We again model these repeat measurements with the equicorrelated block model and the same Glow-type architecture as for the ADNI data set.
The results both data sets show that incorporating dependencies improves the likelihood fit on the holdout test data set. We note that this does not necessarily translate into a higher quality for individual images, but rather into better fit of the full distribution.
3.2.3 STOCK DATA PAIRS
A range of different stock trading and risk management strategies require accurate modeling of the behavior of different stocks (Kole al. 2007). We focus on modeling the daily returns for two pairs of correlated stocks, which is used, e.g., in pair trading strategies (Stander et al. 2013). A pairs trading strategy can utilize a probabilistic model of stock returns as follows: each day, one can assess if a given stock pair lies outside of high-confidence region given the model. If the pair behaves anomalously and one stock underperforms compared to the other stock, trader can hedge these two stocks against each other. The trader would "buy long" the underperforming stock and "sell short" the overperforming stock, with the implicit assumption that in the future the two prices will revert back to high-confidence region. Here, we use the pairs AAPL-MSFT (Apple & Microsoft) and MA-V (Mastercard & Visa), each starting from initial public offering (IPO) of the later of the pair, until late 2017, using publicly available data at close time. A single data point is the 2d daily logarithmic return of one of the two pairs of stocks, For example, MA closed on 2012/002 2012/06/21 with a price of $40.737 and on 2012/06/22 at $42.080, while V closed at $28.661 and 296.76 for those two days. The associated data point then is (ogg(4008/440) (07)). log 2969999999977 288.661) (0.0324, 0.0449), and corresponding data point for the same days for AAPL-MSFT is in the data set. We split data into train (70%), validation (15%), and test (15%) data temporally (non-randomly) to counteract information leakage. Since Apple and Microsoft had their respective IPOs in the 1980s and Visa and Mastercard theirs in the 2000s, the AAPL-MSFT pair is overrepresented in the training data, while both pairs are equally represented in the validation and test data. We use the equicorrelated dependency model with two blocks, one for LAPL-MSET and one for MA-V. The distribution fit for the equicorrelated model is slightly improved using the data dependencies, but again shows that joint optimization appears be inferior to a simple grid search.
4 CONCLUSION
We have shown that through a simple adaptation in the likelihood loss of normalizing flows, we can integrate flexible data dependencies into the training objective, which can also be trained with mini-batch SGD. Experimental evaluation of synthetic and real-world data showed that our method can significantly improve the fit of probabilistic models. In future work, we' 'e especially interested how our method can be extended to other generative models such as VAEs and if it can be combined with other debiasing methods such as causal DAGs as done by Pawlowski al. 2022)).
Our method is not without limitations. The trace estimator in Equation 2 is unbiased, but has a high variance due the overweighting of the off-diagonal terms. This can lead unstable gradient estimates, especially in the early stages of training. In addition, joint optimization of p:s with the flow parameters counterintuitively only sometimes leads to better results. We believe further improvements to the optimization schemes might alleviate these issues.
We also note that, if the goal density estimation or generative modeling, incorporating depen- dencies into the normalizing flow objective is not necessary in all cases with dependent data. For example, in the case of equicorrelated repeat measurements with identical block-sizes, no improve- ments can be expected. This is because region of the sampling space is overrepresented relative to the other regions. Only when some blocks are larger than others (or with more general, unbalanced covariance matrices), adjustment for dependencies makes sense. However, if interested in full likelihood evaluation over the whole data set instead of just density estimation in individual data points, the results will differ.

---

Under review as a conference paper at ICLR 2023
REFERENCES
Sean Augenstein, H. Brendan McMahan, Daniel Ramage, Swaroop Ramaswamy, Peter Kairouz, Mingqing Chen, Rajiv Mathews, and Blaise Aguera y Arcas. Generative models for effective ml on private, decentralized datasets. In International Conference on Learning Representations, 2020. URL httptttt//ppeeeeee..ee ott.a//www....etooooooooooooaaaa
Clare Bycroft, Colin Freeman, Desislava Petkova, Gavin Band, Lloyd T Elliott, Kevin Sharp, Allan Motyer, Damjan Vukcevic, Olivier Delaneau, Jared O' O'Connell, et al. The uk biobank resource with deep phenotyping and genomic data. Nature, 562(7726):233-099 2018.
Lon R Cardon and Lyle Palmer. Population stratification and spurious allelic association. The Lancet, (61 2003.
Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio. Density estimation using real NVP. In International Conference on Learning Representations, 2017. URL /ttpsstss////
net/ ftr.//ormmmmmm do= x
Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios. Neural spline flows. Ad- vances in neural information processing systems, 32, 2019.
Derek Hansen, Brian Manzo, and Jeffrey Regier. Normalizing flows for knockoff-free controlled feature selection. arXiv preprint S.Xiv:2106.01222 2021.
Gary B Huang, Marwan Mattar, Tamara Berg, and Eric Learned- Miller. Labeled faces in the wild: A database forstudying face recognition unconstrained environments. In Workshop on faces in Real-Life Images: detection, alignment, and recognition, 2008
Pavel Izmailov, Polina Kirichenko, Mare Finzi, and Andrew Gordon Wilson. Semi-supervised learn- ing with normalizing flows. International Conference on Machine Learning, pp. 4615-4630.
PMLR, 2020.
Clifford R Jack Jr, Matt A Bernstein, Nick C Fox, Paul Thompson, Gene Alexander, Danielle Har- vey, Bret Borowski, Paula J Britson, Jennifer L. Whitwell, Chadwick Ward, et al. The alzheimer's discase neuroimaging initiative (adni): Mri methods. Journal of Magnetic Resonance Imaging: An Official Journal of the International Society Magnetic Resonance in Medicine, 27(4))685- 691, 2008
Piotr Jaworski, Fabrizio Durante, Wolfgang Karl Hardle, and Tomasz Rychlik, Copula theory and its applications, volume 198. Springer, 2010.
Jiming Jiang and Thuan Nguyen. Linear and generalized linear mixed models and their applica- fions, volume 1. Springer, 2007.
Longda Jiang, Zhili Zheng, Ting Qi, Kathryn E Kemper, Naomi R Wray, Peter M Visscher, and Jian Yang. A resource-efficient tool for mixed model association analysis of large-scale data, Nature genetics, (1(12):144-11555 2019.
Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. CoRR, abs/ (bb//442,.980, 2015.
Durk P Kingma and Prafulla Dhariwal. Glow: Generative flow with invertible 1x1 convolutions. Advances in neural information processing systems, 31, 2018.
Matthias Kirchler, Stefan Konigorski, Matthias Norden, Christian Meltendorf, Marius Kloft, Claudia Schurmann, and Christoph Lippert. transfergwas: Gwas of images using deep transfer learning. Bioinformatics, 38(14):3621- 3628,, 2022.
Erik Kole, Kees Koedijk, and Marno Verbeek. Selecting copulas for risk management. Journal of Banking & Finance, (1):244455-444 2007.
Christoph Lippert, Jennifer Listgarten, Ying Liu, Carl M Kadie. Robert Davidson, and David Heckerman. Fast linear mixed models for gido association studies. Nature methods, (10):833-835, 2011.
10

---

Under review as a conference paper at ICLR 2023
Po-Ru Loh, George Tucker, Brendan K Bulik-Sullivan, Bjarni J Vilhjalmsson, Hilary K Finucane, Rany M Salem, Daniel I Chasman, Paul M Ridker, Benjamin Neale, Bonnie Berger, et al. Efficient bayesian mixed-model analysis increases association power in large cohorts. Nature genetics, 47(3):28--200 2015.
Zachary R McCaw, Jacqueline Lane, Richa Saxena, Susan Redline, and Xihong Lin. Operating characteristics of the based inverse normal transformation for quantitative trait analysis in genome-wide association studies. Biometrics, 76(4): 1262-1272222 2020.
Remo Monti, Pia Rautenstrauch, Mahsa Ghanbari, Alva Rani James, Matthias Kirchler, Uwe Ohler, Stefan Konigorski, and Christoph Lippert. Identifying interpretable gene-biomarker associations with functionally informed kernel-based tests in 190,000 exomes. Nature communications, 13(1): 1-16, 2022.
Didrik Nielsen, Priyank Jaini, Emiel Hoogeboom, Ole Winther, and Max Welling. Survae flows: Surjections to bridge the gap between vaes and flows, Advances in Neural Information Processing Systems, 23:1568-11266 2020.
George Papamakarios, Eric T Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan. Normalizing flows for probabílistic modeling and inference. J. Mach. Learn. Res., 22(57):1 2021.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary VeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high-performance deep learning library. In H. Wallach, H. Larochelle, A. Beygelzimer,
('Alcéé-Buc, E. Fox, and R. Garnett (eds.), Advances in Neural Information Processing Systems 32, pp. 8024-8035. Curran Associates, Inc., 2019. URL ht httpt:///pppaee..srrrspp.ooo.oot
ott.://ww........aaaie.eeeeeaeaa librrar..
Library.
pdf
Nick Pawlowski, Daniel Coelho Castro, and Ben Glocker. Deep structural causal models for tractable counterfactual inference. Advances in Neurul Information Processing Systems, 33:857- 869, 2020.
Alkes L Price, Nick J Patterson, Robert M Plenge, Michael E Weinblatt, Nancy A Shadick, and David Reich. Principal components analysis corrects for stratification in genome- wide association studies. Nature genetics, 38(8)-9049909 2006.
Shaun Purcell, Benjamin Neale, Kathe Todd-Brn Todd-Brown, Lori Thomas, Manuel AR Ferreira, David Bender, Julian Maller, Pamela Sklar, Paul IW De Bakker, Mark J Daly, et al. Plink: a tool set for whole genome association and population-based linkage analyses. The American journal of human genetics, 81(3):559 (753))5555555200000000 2007.
Mohammad Rasouli, Tao Sun, and Ram Rajagopal. Fedgan: Federated generative adversarial net- works for distributed data, arXiv preprint arXiv:2006.0722. 07228, 2020.
Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows. In Interna- tional conference on machine learning, pp. 1530-1538. PMLR, 2015.
Felix Salmon. Recipe for disaster: the formula that killed wall street. Wired Magazine, 17(3):17-03, 2009.
Robin Schirmeister, Yuxuan Zhou, Tonio Ball, and Dan Zhang. Understanding anomaly detec- tion with deep invertible networks through hierarchies of distributions and features. Advances in Neurul Information Processing Systems, 33:21038-21049, 2020.
Giora Simchoni and Saharon Rosset. Using random effects to account for high-cardinality categor- ical features and repeated measures in deep neural networks. Advances in Neural Information Processing Systems, 34, 2021.
11

---

Under review as a conference paper at ICLR 2023
Nasa Sinnott-Armstrone Yosuke Tanigawa, David Amar, Nina Mars, Christian Benner, Matthew Aguirre, Guhan Ram Venkataraman, Michael Wainberg, Hanna M Ollila, Tuomo Kiskinen, et al. Genetics of 35 blood and urine biomarkers in the uk biobank. Nature geneties, 53(2):185-99,, 2021.
Yolanda Stander, Daniël Marais, and Ilse Botha. Trading strategies with copulas. Journal of Eco- nomic and Financial Sciences, 6(1):83-107, 2013.
Yung Liang Tong. The multivariate normal distribution. Springer Science & Business Media, 2012.
M-N Tran, Nghia Nguyen, David Nott, and Robert Kohn. Bayesian deep net glm and glmm. Journal of Computational and Graphical Statistics, 29 2020.
Peter M Visscher, Matthew A Brown, Mark I McCarthy, and Jian Yang. Five years of gwas discov- ery. The American Journal of Human Genetics, 90(1):7-24, 2012.
Peter M Visscher, Naomi Wray, Qian Zhang, Pamela Sklar, Mark I McCarthy, Matthew A Brown, and Jian Yang. 10 years of gwas discovery: biology, function, and translation. The American Journal of Human Genetics, 101(1):5-22 2017.
Yunyang Xiong, Hyunwoo J Kim, and Vikas Singh. Mixed effects neural networks (menets) with applications to gaze estimation. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp- 7743 -7552, 2019.
Jianming Yu, Gael Pressoir, William H Briggs, Irie Vroh Bi, Masanori Yamasaki, John F Doebley, Michael McMullen, Brandon Gaut, Dahlia M Nielsen, James B Holland, et al. A unified mixed-model method for association mapping that accounts for multiple levels of relatedness. Nature genetics, 38(2)2203-208, 2006.
Xiang Zhou and Matthew Stephens. Genome-wide efficient mixed-modee analysis for association studies. Nature genetics, 44(7):821-824 2012.
A
PROOF OF PROPOSITION 2.2
We have
[g = n2
E Ar, 4s + n(n - 1)
Eg EelEvE)LAA, us.
b(b - 1)
For the first term, we know that Eg| = b/n, 50
n A ui-
Σ EglEd Ar. Ui
For the second term, we first note that
b EglEce)] = EglE
EglE =1|= b(b - 1) n(n 1)
n
Then we get
2 (n 1)
Σ EelE uj Σ AnjuT ug. i<j
b(b - 1)
izj
Adding the two terms back together, we get the trace term.
B EXPERIMENTAL DETAILS
All experiments were implemented in PyTorch (Paszke et al. 2019) and PyTorch Lightning, using the normalizing flow implementations provided by Nielsen et al. (2020). In all settings, Adamax optimizer (Kingma Ba, 2015) and reduce the learning rate with an exponential decay. Weight decay (chosen as described below) is always only applied to the weights of the normalizing flows, not the dependency parameters À and Pc-
the
12

---

Under review as a conference paper at ICLR 2023
B.1 SYNTHETIC DATA
B.1.1 EQUICORRELATED DATA
We sample block-sizes from Pareto II distribution with shape parameter mum value 1, rounded to integer values. We clip block-sizes to a maximum of 1,000 and draw new blocks until all blocks together sum n 10, 000 samples. For each block, we draw one Ps ~ Unif(o.5.0.99 and define the full covariance matrix as in Section 2.3.2 Using this covari- ance matrix, we sample non-independently from a bivariate standard normal distribution. We non- linearly transform these data into complex shapes (Abs, Crescents CrescentCubeds Si gn, and S Sinewave) provided by Durkan et al. (2019), for a more challenging density estimation task. We repeat all experiments 10 times with different random seeds.
0.5 and mini-
As a base flow model, we choose rational quadratic spline flows (Durkan et al. 2019), which are state-of-the-art for these challenging data sets. For modelling the equicorrelated blocks, we choose both a grid search over fixed parameters Ó E (0.01, 0.025, 0.05, 0.1, 0.175, 0.25, 0.375, 0.5, 0.6, 0.67, 0.75, 0.9} and joint gradient-based opti- mization of D; and the flow parameters, with starting values for P. E (0101,0.1, 0.25, 0.5). We train all models for 100 epochs, perform small hyperparameter sweep over learning rate (in (0.001, 0.003, 0.01, 0.03)) and weight decay (in {0.001, 0.01, 0.13), and choose the best model for each setting based on early stopping and validation set performance (which is sampled without dependencies).
B.1.2 KNOWN COVARIANCE
In this setting, simulate known covariance matrix between n 5, 000 different samples but with unknown variance component A. We first draw a lower- triangular matrix L, with diagonals all set to I and all elements below the diagonal drawn independently from Unif(o.s.0.99) We use G norm [L as our covariance structure, where norm normalizes the covariance matrix correlation matrix (with all- on the diagonal). We then use the covariance structure AI (1-A1G as in the equicorrelated case to generate correlated bivariate standard normal samples that again get non-linearly transformed. We fix À 0.5 in all experiments, and experiments are again repeated 10
times.
For the grid search, we choose E {0.99, 0.975, 0.95, 0.9, 0.825, 0.75, 0.625, 0.5, 0.4, 0.33, 0.25, 0.13 (note that À corresponds from (0.99, 0.9, 0.75, 0.5). We train for 100 epochs the baseline and in grid search; for the alter- nating optimization, gradient descent updates of Á inbetween. Remaining parameters are chosen as in the equicorrelated simulations.
and for the alternating optimization scheme, we initialize
train for 5 stages of 25 epochs for the flow optimization, with 4 stages of 100
B.2
REAL-WORLD DATA
B.2.1 UKB BIOMARKERS
We used an architecture with 16 affine coupling layers. where the fully connected networks have layers input-122 1eattt-188-12oouutpp for each block, Swish activation functions, and batch-size of 256, as well as a step-wise exponentially-decaynn learning rate schedule.
As in the synthetic experiment, for grid searches we search over E (0.99, 0.975, 0.95, 0.9, 0.825, 0.75, 0.625, 0.5, 0.4, 0.33, 0.25, 0.1} and for the alternating optimization scheme, initialize from (0.99, 0.9, 0.75, 0.5]. For both grid search and baseline, we do a hyperparameter sweep over the learning rate (in {0.001, 0.003, 0.01, 0.03)), weight decay (in {0.001, 0.01, 0.1)), and number of epochs (25 or 50; in a preliminary exploratory sweep we found that more epochs only lead to overfitting). For the alternating optimization, we chose the same learning rate weight decay grid, and additionally optimized over the learning rate for Ā (in f0.03, 0.1, 0.3)) and the number of epochs in the main stages (5 or 25). We alternated for 4 stages, and Ś-optimization stages went for 100 steps each.
13

---

Under review as a conference paper at ICLR 2023
(a) Original data.
(b) Data after marginal
quantil-trrnsfoormatonn
(c) Data after transformation with standard normaliz- ing flow.
(d) Data after transformation with normalizing flow correcting for dependencies.
Figure 3: Marginal histograms (diagonals) and 2-d histograms of pairwise joint distributions for the group of "Hormonal" biomarkers.
14

---

Under review as a conference paper at ICLR 2023
GWAS experiment For each biomarker group, we selected 10,000 individuals at random from those individuals that had values for the corresponding biomarkers. For flow training, we used the same architecture and training for the previous Biomarker experiment. Based on the results from the previous experiment, we fixed learning rates at 0.03 (learning rate for À = 0.1) and weight-decay on the weights at 0.01. To adjust for fixed covariate effects (age, sex, and genotyping batch), projected out covariates from the raw phenotypes with a standard linear regression. For the baseline model, we trained for 250 epochs (this performed considerably better than the fewer epochs in the prior experiment). For the alternating flow, we again trained for 4 alternating steps with 25 (flow- stage) and 100 (A-stage) epochs each. We performed each experiment three times with random seeds for both the selection of individuals and for flow initialization & data loading.
F shows the pairwise joint distributions for the Biomarker group "Hormonal" (after covariates were projected out). Figure 3a shows the original data; Figure 35) shows this data after marginals were transformed using quantile transformation to standard normal values- it is clearly visible that although the marginals are normally distributed, the joint distribution is far from multivariate normal. Figures 3cl & 3d show the data after being transformed with normalizing flows, without and with correcting for dependencies.
Genotype filtering was performed with Plink (Purcell et al., 2007), setting minimum minor al- lele frequency MAF 2 0.1% and Hardy-Weinberg equilibrium p-value p 0.001; and linkage- disequilibrium (LD) pruning with R² = 0.8 and 500kb window. Both univariate (""ingle") and multivariate ("Baseline", "AAternating") GWAS sion 0.98.5 (Zhou & Stephens| (012) with score tests (option - -Imm 3) and centered relatedness matrix (option -k 1). For other GEMMA options we used the defaults, hence, final results were further pruned for MAF 2 5%. This resulted in approximately 500,000 genotypes per experiment, but slightly varying between different random seeds and different biomarker groups, and resulting genome-wied significance threshold of a 0.05/num numm geno 2 10
performed using the GEMMA software ver-
Loci were identified using the Plink clumping utility, defining a locus as group of significantly associated SNPs single-nucleotide polymorphisms) that were both close spatially (within a 250kb window) and in LD with R² 2 0.1.
B.2.2 IMAGE MODELING
For both data sets we used a Glow-like architecture with 2 scales and 12 steps per scale, as imple- mented by Nielsen al. (2020). We grid-searched for A {0.01, 0.025, 0.05, 0.075, 0.1, 0.15) and for joint optimization we initialized with the same parameters. ADNI models were trained for 200 and LFW models for 400 epochs, all with a batch-size of 64. All models were trained with a batch-size of 64 and learning rate and weight decay of 0.001 single A100 GPU. Due to compute constraints, no further hyperparameter exploration was performed.
ADNI brain imaging The data are T1-weighted MRI, preprocessed and standardized with a brain atlas registration pipeline, using brain extraction, linear alignment, non-linear alignment, and debi- asing. The resulting images are more homogeneous than the raw images and thus easier to model. We select the axial-view centered slices and resize them to 64 x 64 grayscale images. The PÚ chosen by the best final model with joint optimization ranged between 0.066 and 0.081.
LFW Here, we used 32 x 32 RGB images. The A chosen by the best final model with joint optimization ranged between 0.052 and 0.15, while the best model with grid optimization with A = 0.15.
B.2.3 STOCK DATA PAIRS
For the stock data, we used an affine coupling normalizing flow with 8 layers of input-64-64-outppe
dimensions and swish activation function. Grid search and joint search were initialized with the same values in the synthetic experiment. We performed a hyperparame- ter sweep over learning rate ((0.001, 0.003, 0.01, 0.03)), weight decay (f0.001, 0.01, 0.13) and ran all models for 100 epochs and batch size of 256.
15

---

Under review as a conference paper at ICLR 2023
C COMPUTATIONAL CONSIDERATIONS
Additional compute & memory requirements for incorporating dependencies depend mostly on the lype of dependencies and on the optimization scheme. In our implementation. baseline runs were implemented as special cases of the flow with dependencies (i.e., 0 or X = 1), which makes fair empirical comparison challenging.
C.I EQUICORRELATEDD BLOCKS
Grid optimization A single with fixed dependency parameter p; O will have almost iden- tical rum times as the baseline method with = O, as the base distribution likelihood evaluation not bottleneck. Since all Pi identical, there is virtually no additional memory requirement. However, as the full network needs to be trained for each the Mgnia grid values tested, the grid evaluation scheme takes roughly Mnselle
Joint optimization In this setting, N (number of blocks) parameters Ps need to additionally be estimated and stored memory, but all cases considered this paper this was strongly dominated by the number of parameters in the model (e.g., in LFW, the normalizing flow model had ~ 90M parameters, but only few thousand extra parameters for the individuals. For very slim models and very large number of blocks, this relationship may change.
C.2 FIXED COVARIANCE
For the fixed-covariance case, a full spectral decomposition is necessary prior to training, which is (in practice) (nn²) operation. It also requires storing the full spectral decomposition in memory. Standard linear algebra libraries used in PyTorch or Numpy & SciPy only support spectral decom- positions up several 10k and oftentimes become unreliable beyond that. Therefore, using fixed covariance schemes is infeasible for larger-scale problems using out-of-the-boy software.
Grid optimization For the fixed grid schedule, mini-batch estimation requires quadratic time in the size of the mini-batch, due to the stochastic trace estimator in Equation 2 However, for batch- sizes used in our settings, this was still dominated by the neural architecture shared with the baseline flow architecture. The log-det-Jacobian can be cached and the remaining parts are identical the baseline flow, so each individual epoch has very similar time requirements to the baseline model. Analogously the equicorrelated blocks grid optimization, although the same spectral decomposition can be used for all those runs.
still need to perform Mgrid runs,
Alternating optimization The main training stage for the flow parameters identical computa- tional considerations as the grid optimization procedure. However, for optimizing every other training stage, first the full data set needs to pushed through the normalizing flow and then rotated with the orthogonal matrix Q from the spectral decomposition. Despite this, the alternating train- ing procedure dominated by the original spectral decomposition and the main training stage of the flow.
D ADNI IMAGES
16

---

Under review as a conference paper at ICLR 2023
VII
II
XEM
Train images
(b) Baseline
(c) Grid Search
(d) Joint
Figure 4: Random samples of ADNI train images and images generated by the normalizing flow models.
17


Reviews:
[{"review_id": "n1Ip_vUYE2U", "reviewer": "Reviewer_gW39", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "The paper proposes a modification to the standard normalizing flow training loss to account for certain dependencies in data. Specifically, the authors consider the special case where datapoints are dependent in latent space, but are transformed pointwise to the observed space. In addition, the authors specialize to the case where the correlation structure is either known up to a convex combination with the identity, or has a block-diagonal structure with known block assignment and constant positive correlation within blocks. The authors then propose training strategies for learning the parameters of these dependencies. The authors evaluate experimentally on simulated and real datasets, using NLL on held-out datasets as the core metric.", "summary_of_the_review": "The proposed method seems fine in terms of specifying an objective that reflects certain forms of independence. However, the actual tasks being solved require far more explicit discussion, since there are many different tasks that one might try to solve in a dependent data context, and it is not clear which of these the method is designed to address. Because of this concern, the paper feels incomplete.", "strength_and_weaknesses": "Strengths:\n + Dependence between datapoints can be an important issue, especially in cases where one cares about being able to make joint or conditional predictions where the task depends on multiple datatpoints simultaneously.\n + The derivations are clearly explained. Although they are straightforward, they are not obfuscated to appear more complicated than they are.\n + Several practical issues in training are addressed.\n + Experiments explore degrees of freedom that were raised in the development of the method.\n\nWeaknesses:\n  - The primary weakness to me is that the method seems undermotivated. Although I know that there exist tasks where dependency between datapoints matter, it is not clear how accounting for dependencies between data translates to better solutions for important tasks. The authors raise several applications, including genomics, finance, and neuroimaging, but stop short of indicating what problems could be solved in these areas if only we modeled dependencies between examples better. NLL is not itself self-motivating; it measures the model's consistency with the data, but does not tell us what this extra fit buys in terms of performance on a real task. What is the unsolved problem in genomics, finance, or neuroimaging that can be solved with a good model of between-example dependence?\n - Along the same lines, when we consider dependent data, the actual prediction task needs to be specified with more care. What exactly is being measured by performance on the test data? Do we care about our having low expected prediction risk on the next datapoint? Low joint prediction risk on a set of $m$ new datapoints? Some kind of statistical inference? In the first case, dependency structure would not matter. In the second case, it might matter, but one would need to motivate why we care about $m$ datapoints simultaneously, and the composition of this test set would need to be motivated carefully (e.g., would this test set represent measurements from a completely new individual? from an individual we have observed before?). Also, we might need the test set to be composed of $n$ replicates of $m$ datapoints. In the third case, the form of statistical inference would need to be made clear (e.g., Are we doing anomaly prediction and constructing p-values? Are we using the normalizing flow to do covariance estimation for a downstream regression / GWAS task?).\n - For example, it is not clear what the goal of the task illustrated in Figure 1 is. If the goal were density estimation, it seems like the standard normalizing flow is doing the right thing. Meanwhile, it is not clear why we want the dependency-adjusted normalizing flow to match the distribution of data sampled independently. This comes back to the issue that when one has dependent data, the goals are not self-evident in the same way that they are when predicting with independent data. Goals and notions of bias need to be made explicit.\n - Given the note in the discussion that \"equal block sizes\" would be a case where the dependent data objective would produce no changes, did the authors consider the baseline of weighting datapoints by inverse block size? It seems that for some tasks this might address many issues and not require major modifications to training.\n - Many tricks for working with covariance matrices have been discussed before in literatures that do efficient ridge regression or work with Gaussian processes. I am not sure that the linear algebra tricks need to be rehashed in the main body of the paper, and this could make room for better discussion of tasks where dependence matters.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "_xuxec_HggN", "reviewer": "Reviewer_Krj4", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "The paper tries to tackle the problem of learning a generative model for sampling identically distributed but dependent data points from a distribution. Like previous normalizing flows, they try to learn a bijective map of data to an assumed latent space but unlike previous works which assume a diagonal Gaussian, they assume a correlation structure within the latent space as well. This makes the computation of probability density estimates in the latent space intractable in general but by restricting to specific classes of correlation structures the authors design a computation scheme that works for their case of dependent data. The authors test their method on toy and real world datasets.\n", "summary_of_the_review": "I feel the authors need to strengthen the experiments section considerably in an attempt to justify that dependency modelling is actually useful [other than slight NLL benefits]. Other than that I think the paper is a good attempt at incorporating data dependency assumptions into the normalizing flow literature. I am willing to change my rating based on discussions on the forum.", "strength_and_weaknesses": "Strengths:\n\n- The authors chose an oft ignored area in generative modeling: that of dependent sampling which is interesting to analyze.\n- The paper is mostly well written and flows well.\n\nWeaknesses:\n- I am not quite sure why the specific covariance structures were chosen. As an aside do the authors think that adding a regularisation functional that minimises the norm of $\\rho$ will benefit the joint optimization strategy. In a way, would that choose the most decorrelating bijective map?\n- The authors should provide a compute budget difference (both time and memory) between the baseline and proposed methods. It seems the baseline is just a quadratic spline flow whereas the proposed method has additional parameters for the latent correlation structure. An aside: have the authors attempted to optimize the baseline latent mean and covariance?\n- The experiment section lacks many details:\n  -- In UKB: how is the relatedness used? Do you assume G is known?\n  -- In stocks: a few samples would help  understand what is log return referring to? Is it tracking per day movement of a pair? How is the temporal nature of stocks over time incorporated in the model?\n  -- ADNI: To my eyes all images in figure 3 look very similar. Does dependent modeling help in downstream tasks? If so, it would strengthen the author’s case a lot if they could show something of that sort.\n\nOverall I feel the authors would greatly benefit in relegating some exp details to the appendix (around values in grid search, extra dataset details) and adding more results from each of the datasets. For ex. in UKB they say a quantile transformation makes statistical modeling hard. It would strongly benefit the authors to show that their way of modeling makes that easier.\n\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "1Xmq7iCaOb", "reviewer": "Reviewer_xwvU", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "The authors propose an extension of normalizing flows to non-IID data. This extension can be seen as the generalization of a Gaussian copula to vector-valued variables. The authors propose parameterizations of the resulting model and associated learning algorithms that yield efficient training.", "summary_of_the_review": "I can see this paper being a useful addition to this line of literature, although I don't see it as completely ground-breaking, It could benefit from a more developed and more thoroughly studied use case.", "strength_and_weaknesses": "STRENGTHS:\n* The paper studies an important problem. Although the motivation describes it as non-IID data, this method also applies to high-dimensional correlated vectors, and may be useful for modeling IID data with interesting correlations among the features.\n* This work is novel, and to my knowledge this techniques of extending copula to higher-dimensional datasets is new.\n* The algorithmic ideas proposed in the work are novel, technically sound, and non-trivial. The level of technical depth required to derive the methods is significant.\n* The experiments are extensive and look at a wide range of interesting real-world problems.\n\nWEAKNESSES:\n* These are all relatively minor comments. I think it would have been interesting to describe the motivation and the various applications of the method early on (in the introduction). I only \"got\" the idea and why it's interesting after reading the whole paper. I didn't feel like the intro did a good job at selling the work.\n* While the experiments look at many domains, I didn't feel like the method achieved very strong improvements. Highlighting specific examples where this method really shines could make the paper stronger.\n* The paper mentions GWAS and ancestry-based confounding a lot, and it sounds like an interesting motivating application. It could help make the paper stronger to focus and commit on one application problem more deeply, use it as a clear consistent running example, show strong results on it, etc.\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject