Training Energy-Based Normalizing Flow with Score-Matching Objectives
Authors: Chen-Hao Chao
Wei-Fang Sun
Yen-Chang Hsu
Zsolt Kira
Chun-Yi Lee
Venue: NeurIPS 2023 poster (2023)
Keywords: flow-based models
score-matching methods
TL;DR: In this paper, we establish a connection between the parameterization of flow-based and energy-based models, and present a new flow-based modeling approach called energy-based normalizing flow (EBFlow).

Abstract: In this paper, we establish a connection between the parameterization of flow-based and energy-based generative models, and present a new flow-based modeling approach called energy-based normalizing flow (EBFlow). We demonstrate that by optimizing EBFlow with score-matching objectives, the computation of Jacobian determinants for linear transformations can be entirely bypassed. This feature enables the use of arbitrary linear layers in the construction of flow-based models without increasing the computational time complexity of each training iteration from $\mathcal{O}(D^2L)$ to $\mathcal{O}(D^3L)$ for an $L$-layered model that accepts $D$-dimensional inputs. This makes the training of EBFlow more efficient than the commonly-adopted maximum likelihood training method. In addition to the reduction in runtime, we enhance the training stability and empirical performance of EBFlow through a number of techniques developed based on our analysis of the score-matching methods. The experimental results demonstrate that our approach achieves a significant speedup compared to maximum likelihood estimation while outperforming prior methods with a noticeable margin in terms of negative log-likelihood (NLL).

Training Energy-Based Normalizing Flow with Score-Matching Objectives
Chen-Hao Chao', Wei-Fang Sun!?, Yen-Chang Hsu Zsolt Kira*, and Chun-Y Yi Lee!
Elsa Lab, National Tsing Hua University, Hsinchu City, Taiwan NVIDIA AI Technology Center, NVIDIA Corporation. Santa Clara, CA, USA
3 Samsung Research America, Mountain View, CA, USA Georgia Institute of Technology, Atlanta, GA, USA
Abstract
In this paper, we establish a connection between the parameterization of flow-based and energy-based generative models, and present a new flow-based modeling ap- proach called energy-based normalizing flow (EBFlow). We demonstrate that by optimizing EBFlow with score score-matching objectives, the computation of Jacobian determinants for linear transformations can be entirely bypassed. This feature enables the use of arbitrary linear layers in the construction of flow-based models without increasing the computational time complexity of each training iteration from (D) to ( L) for an L-layered model that accepts D D-dimensional in- puts. This makes the training EBFlow more efficient than the commonly-adope maximum likelihood training method. In addition to the reduction in runtime, we enhance the training stability and empirical performance of EBFlow through a number of techniques developed based on our analysis of the score- mooee-mtchhing meth- ods. The experimental results demonstrate that our approach achieves a significant speedup compared to maximum likelihood estimation while outperforming prior methods with a noticeable margin terms of negative log-likelihood (NLL).
Introduction
Parameter estimation for probability density functions (pdf) has been a major interest in the research fields of machine learning and statistics. Given D-dimensional random data vector TR. the goal of such a task is to estimate the true pdf (x(+) of x with a function p(- 0) parameterized by 0. In the studies of unsupervised learning, flow-based modeling methods (e.g.. [I-T] are commonly-adopted for estimating Px due to their expressiveness and broad applicability in generative tasks.
Flow-based models represent p(:; 0) using sequence of invertible transformations based on the change variable theorem, through which the intermediate unnormalized densities are ce-noraaliedd by multiplying the Jacobian determinant associated with each transformation. In maximum likelihood estimation, however, the explicit computation of the normalizing term may pose computational challenges for model architectures that use linear transformations, such as convolutions [4] 5 and fully-connected layers 6 77]. To address this issue, several methods have been proposed in the recent literature, which includes constructing linear transformations with special structures [8] 12 and exploiting special optimization processes [7]. Despite their success in reducing the training complexity, these methods either require additional constraints the linear transformations biased estimation on the gradients of the objective.
Motivated by the limitations of the previous studies, this paper introduces approach that reinterprets flow-bbsed models as energy-basee based models [13), and leverages score-matching methods [14 to
*Corresponding author. Email: cyleeêcs nthu, edu.tu tw
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

optimize p(:;0) according to the Fisher divergence [14] [8] between Px(-) and p(+;0), 0). The proposed method avoids the computation of the Jacobian determinants of linear layers during training, and reduces the asymptotic computational complexity of each training iteration from O(D I.) to O(D2L)
for an I--layered model. Our experimental results demonstrate that this approach significantly improves the training efficiency as compared to maximum likelihood training. In addition, we investigate a theoretical property of Fisher divergence with respect to latent variables, and propose Buuch-after-Prepooeesst
(MaP) technique to enhance the training stability of score-matching methods. Finally, our comparison on the MNIST dataset [T9] reveals that the proposed method exhibit significant improvements in comparison to baseline methods presented in [17] and [7] in terms of negative log likelihood (NLL).
Background
In this section, discuss the parameterization of probability density functions in flow-based and energy-based modeling methods, and offer a number of commonly-used training methods for them.
2.1 Flow-based Models
Flow-based models describe (xx using a prior distribution P...) of a latent variable u E RD and an invertible function g = 9L o .... o 9t, where gi(-;0) : IB RD,, Vi E {I.... ;..,,,,, ز and is usually modeled neural network with L layers. Based on the change of variable theorem and the distributive property of the determinant operation det p):p(.;)) can be described as follows:
L
p(x; 0) = Pu (g(a;; 0)) [det(J, 0)) = Pu (9(x; 0)) II |ett(J, -1)-
(1)
s=1
where x; - o 91 (z; 0), 2o J, (0; respect to æ. and g. (z)-1;0) with respect to 1. This work concentrates on model architectures employing linear flows (0) to design the function These model architectures primarily utilize linear transformations to extract crucial feature representations, while also accommodating non-linear transformations that enable efficient Jacobian determinant computation. Specifically, let S) be the set of linear transformations in g. and Sn {9. i {1.... }}} \ S be the set of non-linear transformations. The general assumption of these model architectures is that I |det (Jg.) ) in Eq. (1) can be decomposed as lIg /det (J gt Il,es, |et (J ), where II |det 9.)) and Ig |det (J,,)) )l can be calculated within the complexity of ODD* L) and O respectively. Previous implementations of such model architectures include Generative Flows (Glow) 4 Neural Spline Flows (NSF) (5), and the independent component analysis (ICA) models presented 71.
Ba (0; 0) represents the Jacobian of g with 9. represents the Jacobian of the i-th layer of g
Given the parameterization of p(- ;)). commonly used approach for optimizing 0 is maxi- mum likelihood (ML) estimation, which involves minimizing the Kullback-Leibler (KL) diver-
gence DKL [px(=)))p(<;: 0)] (p_(=) loo Px(x)) between the true density Px(x=) and the param- p(a;;0]
eterized density p(z; 0). The ML objective (f) is derived by removing the constant term (px(#) [logp>(>)) Px with respect to from DRL. [px(<)) 0)), and can be expressed as follows:
CML (0) = (m)) [- p(x; 0)]
The ML objective explicitly evaluates p(:; f),, which involves the calculation of the Jacobian determi- nant of the layers in St. This indicates that certain model architectures containing convolutional [4] 5 fully-connected layers (6] 7 may encounter training inefficiency due to the O(D*L) L) cost of evaluating I9 jdet (Jg,)l- Although a number of alternative methods discussed in Section 3 can be adopted to reduce their computational cost, they either require additional constraints on the linear transformation biased estimation the gradients of the ML objective.
2.2 Ensrgy-based Models
Energy-based models are formulated based on a Boltzmann distribution, which expressed as the ratio of an unnormalized density function to an input-independent normalizing constant. Specif-

---

ically, given a scalar-valued energy function E(.;0) 0) RD - R, the unnormalized density func- tion defined as exp (=E(x; f)). and the normalizing constant Z(0) is defined as the integration ImeR" exp (- Fa;;; 0)) dx. The parameterization of p(:;0) is presented in the following equation:
p(x;0) = exp (- (=E(x;;)) 0)) Z (0).
Optimizing p(:; in Eq. 3 through directly evaluating LMI. in Eq. (2) is computationally infeasible, since the computation requires explicitly calculating the intractable normalizing constant Z(0). To address this issue, a widely- technique [13] is to reformulate de (ML( as its compling-based
variant 0 LsMu(()) which is expressed as follows:
88
LSMI. (0) (x)) [a;;; 0)] = [F(a; 0)),
where sg(-) indicates the stop-gree cuop-gradeent operator, Despite the fact that Eq. (4 prevents the calculation of (0), sampling from p(-;0) typically requires running Markov Chain Monte Carlo (MCMC)
process (c.g., [1 (2)) until convergence, which can still be computationally expensive as it involves
evaluating the gradients of the energy function numerous times. Although several approaches 23 24
were proposed to mitigate the high computational costs involved in performing an MCMC process, these approaches make use of approximations, which often cause training instabilities high- dimensional contexts (25).
Another line of researches proposed optimize p(-;0) through minimizing the Fisher divergence
Dy [px (1)))pp((<;)) 0)] (x)) b11 ē Dm
log
(x) p(m;9)
1|2) between Px(x) and p(æ; 9) using the score-
matching (SM) objective LSM E (x)
! Da
Fla;; T
0a
Ea:; 0))
14 avoid the
explicit calculation of (0) well as the sampling process required in Eq. (4). Several computation-
ally efficient variants of LsM. including sliced score matching (SSM) [16], finite difference sliced
score matching (FDSSM) [17, and denoising score matching (DSM) [1), have been proposed.
SSM is derived directly based on LsM with an unbiased Hutchinson's trace estimator (26). Given random projection vector E IR' drawn from Py and satisfying (n, [u) = I, the objective function denoted as LSSM: is defined as follows:
( 0)
| - (px(()))
x 82 Fla;; (rFFF((;)))) 0)
LSSM = 2
(5)
dx
dx
FDSSM is a parallelizable variant of LSSM that adopts the finite difference method 27 to approximate the gradient operations in the objective. Given uniformly distributed random vector e, it accelerates
the calculation by simultaneously forward passing F((;;) 0), F(xx + E; 0). and E( E; 0) as follows:
LFDSSM( (x))
E(x; f)] (x) [E(x 0) + E(x 0)]
1 + B (m))(()) 8 B
(6)
(E(x + e; E(x -(x==;))) 0))
where Pę(E) = U(E E ||||| = E), and is hyper-paramette that usually assumes a small value.
DSM approximates the true pdf through surrogate that is constructed using the Parzen density estimator Po (σ) [28). The approximated target Pa (æ) = SaeRe Po ( (x))da da defined based on an isotropic Gaussian kernel p, (#|x)= = A o² I) with a variance o², The objective LDSM. which excludes the Hessian term LSSM, is written as follows:
0E 0)
x T
LDSM (0) = Epp
To conclude, CSSM is an unbiased objective that satisfies bu CSsM - (x) 0) 16 while CRDSSM and LOSM require careful selection of hyper-parameter & and a, since 8 (Dr [px (xx()))l((() 0)] + o (E)) [17 contains an approximation error o (E), and
8 CFDSSM(O) 0e
Po in a 56 @ DD [po (5) (( 0)] may bear resemblance Px only for small 15]
of (D)M

---

3 Related Works
3.1 Accelerating Maximum Likelihood Training of Flow-based Models
A key focus in the field of flow-based modeling is to reduce the computational expense associated with evaluating the ML objective [7-12 29]. These acceleration methods can be classified into two categories based on their underlying mechanisms,
29 have attempted to accelerate the computation of Jacobian determinants in the ML objective by exploiting linear transformations with special structures. For example, the authors in [8] proposed constrain the weights in linear layers as lower triangular matrices to speed up training. The authors in 10 proposed to adopt convolutional layers with masked kernels to accelerate the computation of Jacobian determinants. The authors in [29] leveraged orthogonal transformations to bypass the direct computation of Jacobian determinants. More recently, the authors in [12] proposed to utilize linear operations with special butterffy structures [30] to reduce the cost of calculating the determinants. Although these techniques avoid the O(D* computation, they impose restrictions on the learnable transformations. which potentially limits their capacity capture complex feature representations, as discussed in 7 31 32). Our experimental findings presented in Appendix [A.5] support this concept, demonstrating that flow-based models with unconstrained linear layers outperform those with linear layers restricted by lower upper triangular weight matrices g or those using lower-upper (LU)
Specially Designed Linear Transformations. majority of the existing works
I2
decomposition
4
Specially Designed Optimization Process. To address the aforementioned restrictions, recent study 7 proposed the relative gradient method for optimizing flow-based models with arbitrary linear transformations. In this method, the gradients of the ML objective are converted into their relative gradients by multiplying themselves with W"W, where W E RDDD represents the = W, evaluating relative
weight matrix in a linear transformation. Since 3 aw log /det gradients more computationally efficient than calculating the standard gradients according to (W")-11 While this method reduces the training time complexity from O(D*L) to O(D² a significant downside to this approach is that it introduces approximation errors with a magnitude of o(W). which can escalate relative to the weight matrix values.
aW log dee (W)l
3.2 Training Flow-based Models with Score-Matching Objectives
The pioneering study [4] is the earliest attempt to train flow-based models by minimizing the SM objective. Their results demonstrate that models trained using the SM loss are able to achieve comparable or even better performance to those trained with the ML objective in a low- dimensioonn experimental setup. More recently. the authors in [16] and [17] proposed two efficient variants of the S.M loss, ie.. the SSM and FDSSM objectives, respectively. They demonstrated that these loss functions can be used to train non-linear independent component estimation (NICE) [ model on high-dimensional tasks. While the training approaches of these works bear resemblance to ours, our proposed method places greater emphasis training efficiency. Specifically, they directly implemented the energy function E(x; 0) in the score-matching objectives as - log p(æ; 6). resulting in significantly higher computational cost compared to our method introduced in Section In Section 5 we further demonstrate that the models trained with the methods in [6 17 yield less satisfactory results in comparison our approach.
Methodology
In this section, we introduce a new framework for reducing the training cost of flow-based models with linear transformations, and discuss number of training techniques for enhancing its performance.
4.1 Energy-Based Normalizing Flow
Instead of applying architectural constraints to reduce computational time complexity, we achieve the same goal through adopting the training objectives of energy-based models. We name this approach as Energy- Normalizing Flow (EBFlow). key observation is that the parametric density function of a flow-based model can be reinterpreted that of an energy- based model through identifying

---

the input- independent multipliers in p(-; 0). Specifically, p(=;0) () can be explicitly factorized into an unnormalized density and a corresponding normalizing term as follows:
p(:; 0) Pu (g(:; 0))
II det (Jg, (a:;
Pu (0(2; 0))
II 9.E8,
deet (J,, (2; 0))1
II 9.E3:
[det(J (0)) - exp(<E(ax;) 2 2
Unnormalized Density Norm. Const.
Unmomaalized Density
Norm. Const.
(8)
where the energy function E(:; 0) and the normalizing constant (0) are selected as follows:
E(x;; 0) _ log Pu (g(x;; 0)) II (Jtt(Jg, (x:-1;;)) 2 (0) II S.
|det (J,, gJ,,
(9)
S.ES.
The detailed derivations of Eqs. (8) and (9) elaborated in Lemma Section A.1.2 By iso-
lating the computationally expensive term in p(:;0) as the normalizing constant 20),, the parametric pdf defined in Eqs. (8) and ( becomes suitable for the training methods of energy-based models. In the subsequent paragraphs, we discuss the training, inference, and convergence property of EBFlow.
Training Cost. Based on the definition in Eqs. ( and (9). the score-matching objectives specified in Eqs. ( can be adopted to prevent the Jacobian determinant calculation for the elements in S, As result, the training complexity can be significantly reduced to (()² the O L) calculation of (0) is completely avoided. Such a design allows the use of arbitrary linear transformations in the construction of flow-based model without posing computational challenge during the training process. This feature is crucial to the architectural flexibility of a flow-based model. For example, fully-connected layers and convolutional layers with arbitrary padding and striding strategies be employed in EBFlow without increasing the training complexity. EBFlow thus exhibits an enhanced flexibility in comparison to the related works that exploit specially designed linear transformations,
Inference Cost. Although the computational cost of evaluating the exact Jacobian determinants of the elements in S still requires O L) time, these operations can be computed only once after training and reused for subsequent inferences, since 2()) is a constant as long as is fixed. cases where D is extremely large and (() cannot be explicitly calculated, stochastic estimators such as the importance sampling techniques (e.g., B3 34) can be used as an alternative to approximate Z(A).
We provide brief diseussion of such scenario in Appendix A.33
Asymptotic Convergence Property. Similar to maximum likelihood training, score-matching meth- ods that minimize Fisher divergence have theoretical guarantees on their convistency [14] [6] This property essential in ensuring the convergence accuracy of the parameters. Let N the num- ber of independent and identically distributed (i.id.)) samples drawn from Px to approximate the expectation in the SM objective. In addition, assume that there exists set of optimal parameters 0* such that p(æ; 0*) = Px(x). Under the regularity conditions (i.c., Assumptions A. A.7 shown in Appendix I.I.I). consistency guarantees that the parameters O, minimizing the SM loss converges
(in probability) to its optimal value 8* when N - oc, i.c., 0 2 0" as N + oc. In Appendix 1.1
we provide a formal description of this property based on [16 and derive the sufficient condition for g and Pa to satisfy the regularity conditions (i.e., Proposition 1.
4.2 Techniques for Enhancing the Training of EBFlow
As revealed in the recent studies [6 17]. training flow-based models with score-matching objectives is challenging as the training process is numerically unstable and usually exhibits significant variances. techniques: match after preprocessing (MaP) and exponential moving average (EMA), which are particularly effective in dealing with the above issues according to our ablation analysis in Section 5.3
To address these issues, we propose to adopt
MaP. Score-matching methods rely on the score function - E(x; f) to match & log Px Px(x), which requires backward propagation through each layer in g. This indicates that the training process could be numerically sensitive to the derivatives of For instance, logit pre-processing layers commonly used flow-based models (e.g., 4 7 g 351) exhibit extremely large derivatives near 0 and 1, which might exacerbate the above issue. To address this problem, propose to
de
5

---

exclude the numerically sensitive layer(s) from the model and match the pdf of the pre processed variable during training. Specifically, let Xx gt. o.... 91 (x) be the pre-processed variable, where ls represents the index of the numerically sensitive layer. This method aims to optimize a parameterized
pdf p): & Pu(gL o ..... o 9k+11:; (::0) II
|det (Jgz)l that excludes (....... .9) through minimizing the Fisher divergence between the pdf Pxa (-) of X and ()) by considering the (local)
behavior of Dr, as presented in Proposition 4.1
Proposition 4.1. Let px, be the pdf of the latent variable of x, - 9, (1(x) indexed by
j. In addition, let P,(+) be a paf modeled as Pulgr, o 1)-1(+)) Ih=j+t1 det (Jg.), where
13. j E {0, .. L 1). It follows that:
Dp [Px, [pp,]] = + Dr lx [pol] = O, Vj (1,
13.
(10)
The derivation is presented in Appendix A.1.3] In Section 5.3 we validate the effectiveness of
the MaP technique on the moore-mat methods formulated in Eqs. (5)-7) through an ablation analysis. Please note that MaP does not affect maximum likelihood training, since it always satisfies
DRL [xx, [p]] = D2 [px [Ipo].
Vj {I, L 1} as revealed in Lemma A.12
EMA. In addition to the MaP technique, we have also found that the exponential moving average (EMA) technique introduced in [3 is effective in improving the training stability. EMA enhances the stability through smoothly updating the parameters based on 0 md (1 m m))0 each training iteration, where Ô is set of shadow parameters [36], is the model's parameters iteration i, and m is the momentum parameter. In our experiments presented in Section S we adopt m 0.999 for
both EBFlow and the baselines.
5 Experiments
In the following experiments, we first compare the training efficiency of the baselines trained with LMI and EBFlow trained with LSMI., LSSM. LFDSSM. and LDSM to validate the effectiveness of the proposed method in Sections 5. and 5.2 Then, in Section [.3 we provide an ablation analysis
of the techniques introduced in Section 4.2 and a performance comparison between EBFlow and number of related studies [6 [7]. Finally, in Section we discuss how EBFlow can be applied to generation tasks. Please note that the performance comparison with 12 29 is omitted, since their methods only support specialized linear layers and are not applicable to the employed model architecture (7 that involves fully connected layers. The differences between EBFlow, the baseline, and the related studies are summarized in Table A4 in the appendix. The sampling process involved in the calculation of LSML, is implemented by g "(u 0), where Pu. The transformation g) for each task is designed such that S + á and Sn of é. For more details about the experimental setups, please refer to Appendix A.2
5.1 Density Estimation on Two-Dimensional Synthetic Examples
In this experiment, we examine the performance of EBFlow and its baseline on three two-dimensional synthetic datasets. These data distributions are formed using Gaussian smoothing kernels to en- sure Px (xr is continuous and the true score function Da log ps (x) well defined. The model g(: constructed using the Glow model architecture (4. which consists of actnorm layers, affine coupling lay- ers, and fully-connected layers. The performance are evaluated in terms of the KL divergence and the Fisher divergence between (xxxx) and p(æ; 0) using independent and identically distributed (i.i.d.) testing sample points.
M
M
Reseline (ML) (SML)
FBFFlow
FRFIow (SSM)
FBFlow (DSM)
EBFlow (FDSSM)
Trie
Figure 1: The visualized density functions on the Sine, Swirl, and Checkerboard datasets. The column "True" illustrates the visualization of the true density functions.
Table I and Fig. demonstrate the results of the above setting. The results show that the performance of EBFlow trained with LsSM, LEDSSM, and LDSM in terms of KL divergence is on par with those trained using LSML. as well as the baselines trained using LMIL. These results validate the efficacy of training EBFlow with score matching.

---

Table 1: The evaluation results in terms of KL-divergence and Fisher-divergenee of the Row-based models trained with LML, LSML. LSSM, LDMM, and LEDSSM on the Sine, Swirl, and Checkerboard datasets. The results are reported as the mean and 95% confidence interval of three independent runs.
Dataset
Metric
Baseline (MI.
EBFlow (SML)
EBFlow (SSM)
EBFlow (DSM)
EBFlow (FDSSM)
Fisher Divergence (4)
6.86 = 0.73 e-]
6.65 = 1.05 e-1
6.25 0.54 e-1
6,66 = e-M
6.66 + 1-33 4.57 ::
KL Divergence (4)
0.00
0.00
4.56
Swirl
Fisher Divergence (1)
K1. Divergence (1)
4.21 0.00 e10
0.48 c+0
4.21 + 0.01 ev
0.53
4.25 0.04 exd
1.35
4.22 : 0.02 eid
e+0
4,25 :: 0.08 ei0
1.37
Checkerhoard
Fisher Divergenee ($)
7:34 + 11. c1.50
1.23 + 0.75 e40
7.07 2: 1.93
7,03 = 1.99
7:08 + 1.62
KL Divergence (1)
4,89 0.02
4.81 # 0.02
4,85 0.05 c+0
4,82 = 0.05 c+O
4.83 : 0.03 c=0
Table 2: The evaluation results in terms of the performance (i.e., NLL and Bits/Dim) and the throughput (i.e., Batch/Sec.) of the FC-based and CNN-based models trained with the baseline and the proposed method MNIST and CIFAR-10. 10. Each result is reported in terms of the mean and 95% confidence interval of three independent runs after A is converged. The throughput is measured on NVIDIA Tesla V100 GPUs
MN - 7840
Midell
IC-boso
CNN-tessd
Punar.
0.027
Metbod
Bescline (ML)
EBRow (SMLL)
EBFox (SSM)
EBRlww (DSM)
EBElcw (FDSSM)
Buselie (MIL)
EBPlww (SMIL)
EBFlww (SSM) EBHlw (DSM:
BBFlow (FDSSM)
NLL (1)
100244:001 0.1
1092.3 a
10
1099.2 =
=104.:
1101.3
1048,3 6.6
1107.5
11095 2.4
=12
201 :: 0.00
2.00 0.00
201 0.60
2.02 0.00
2.63 = 0.00
208 a 0.00
2 0.01
2.08 à 0.00
2.04 0.00
2,06 2 0.01
Besa/Se. (T)
8.00
12.27
33.81
06.67
176.31
0.21
0.29
7.0M
18.32
36.76
CIPAR 10 3,002)
Model
Resed
CNN Iusel
Nom. Ponun
18.831 M
00241 M
Mathal
Buncline (ML)
EBR (SMIL)
EBFam (SSM)
EiRRow (DSM)
EBFlow (FUSSM)
Saddii: (MI))
EBFlow (SMIL)
EBFlww (SSM) ESPPor (DSM) EBBlow (FDSSM)
(i)
11912.91 1: 10.5
11915.6 56
11917.7 I8S
11940.0
12547. 1.
11488.7 26.7
1188388 151.7
114855 a 12.0
11462.5 1.9
11766.0 36.8
BraDii (.)
5.00
5.00
0.1
566
0.00
=:
536 =: 0,01
0.07
5.38
0.00
BascrSec: (T)
5.05
7.35
57.14
62.50
0.2
7.35
39.84
5.2 Efficiency Evaluation on the MNIST and CIFAR-I0 Datasets
In this section, we inspect the influence of data dimension D on the training efficiency of flow-based models. To provide a thorough comparison, we employ two types of model architectures and train them on two datasets with different data dimensions: the MNIST [9] (D 1 28 x 28) and
CIFAR-10 37 (D = 3 x 32 x 32) datasets.
The first model architecture is exactly the same that adopted by 77. It is architecture consisting of two fully-connected layers and a smoothed leaky ReLU non-linear layer in between. The second model is parametrically efficient variant of the first model. It replaces the fully-connected layers with convolutional layers and increases the depth of the model to six convolutional blocks. Between every convolutional blocks, squeeze operation (2) is inserted to enlarge the receptive field. In the following paragraphs, we refer to these models as 'FCebased' and 'CNN-based' models, respectively.
The performance of the FC-based and CNN-based models is measured using the negative log likeli- hood (NLL) metric (i.e., -- log p(:; log 0)), which differs from the intractable KL divergence by constant. In addition, its normalized variant, the Bits/Dim metric (38), is also measured and reported. The algorithms are implemented using PyTorch [39] with automatic differentiation [40], and the runtime is measured on NVIDIA Tesla V V100 GPUs. In the subsequent paragraphs, we assess the models through scalability analysis, performance evaluation, and training efficiency examination.
Scalability, To demonstrate the scalability of KL-divergence--aat based (i.e., LML and LsMl) and Fisher- divergence- based (i.e., LSSM. LDSM. and LFDSSM) objectives used in EBFlow and the baseline method, we first present a runtime comparison for different choices of the input data size D. The results presented in Fig. 2 (a) reveal that Fisher-divergence--asss efficiently than KL-divergence-basee objectives. Moreover, the compling-based objective LSMI. used in EBFlow, which excludes the caleulation of 7()) in the computational graph, can be computed slightly faster than LML adopted by the bascline.
objectives can be computed more
7

---

Bootin (ila
(MI1)
(SMM)
EIiikon (BSA))
Rarninn
6:22
MNNSS
Blts: Dim
6.0
03-
24
22
Dimemion (1,
22
00
26
10
60
25
73
10.
Suntion
(8cc (8a22
6:30
53-
CIPAR 10
11mm
12
0.0
0.3 1.31
2.3
15.8 28.2
320
300
Dimersion D...0) FC-hased
Dimassion (L CNN-based
Time (lour)))
(heupo)
PC
CNN--bused
Figure 2: (a) A runtime comparison of calculating the gradients of different objectives for different input sizes (D). The input sizes are (1, n. n and (3, n,, n)), with the x-axis the figures representing In the format (c, h, u),, the first value indicates the number of channels, while the remaining values correspond to the height and width of the input data. The curves depict the evaluation results in terms of the mean of three independent runs. (b) A comparison of the training efficiency of the FC-based and CNN-based models evaluated on the validation set of MNIST and CIFAR-10. Each curve and the corresponding shaded area depict the mean and confidence interval of three independent
Wio MaP Teclnique
WT MaP Toclnique
Table 3:
The results in terms of NLL of the FC- based and CNN-based models trained using SSM, DSM, and FDSSM losses on MNIST. The perfor- mance is reported in terms of the means and 95% confidence intervals of three independent runs.
130
Grallant Nom Nomn
50
FC-hosdd
EMA
MaP
EDFlow(DMM
EDFlow(DDSM)
18.0
1720.5 088
44550
3166
300
1000
10928 + 0.3
1099.2 0.2
1104.1 = 0.5
Training
CNN-bbssed
80 CSSM (0) of an FC- based model trained on the MNIST dataset. The curves and shaded area depict the mean and 95% confidence interval of three independent runs.
Figure 3: The norm of
EMA
HBFioWSSM)
LBFHowDDM)
BBFFOWFDDSMM
35180 4. 33,9
3170.0 7.2
3593.3 12.5
3504.5
11075 + 1.4
1109.5 2.6
11221 =
Performance. Table demonstrates the performance of the FC-based and CNN-based models in terms of NLL the MNIST and CIFAR-IO datasets. The results show that the models trained with Fisher-divergencc--asedd objectives are able to achieve similar performance as those trained with KL-divergenice-based objectives. Among the fisher-divergence--alee
objectives, the models trained using LSSM and LDSM are able to achieve better performance in comparison to those trained using LFDSSM- The runtime and performance comparisons above suggest that LSSM and LDSM can deliver better training efficiency than LML and LSML. since the objectives can be calculated faster while maintaining the models' performance on the NLL metric.
Training Efficiency. Fig. 2 (b) presents the trends of NL. versus training wall time when LML. LSML. LSSM. LDSM. and LFDSSM are adopted the objectives. It is observed that EBFlow trained with SSM and DSM consistently attain better NLL in the early stages of the training. The improvement is especially notable when both D and L large, revealed for the scenario of training CNN-based models on the CIFAR-10 dataset. These experimental results provide evidence to support the use of score-matching methods for optimizing EBFlow.
5.3 Analyses and Comparisons
Ablation Study. Table 3 presents the ablation results that demonstrate the effectiveness of EMA and MaP techniques. It is observed that EMA is effective in reducing the variances. In addition, MaP significantly improves the overall performance. To further illustrate the influence of the proposed MaP technique on the score-matching methods, compare the optimization pro-

---

D
4
4
8
D
9
0
r
3
6
3
2
8
7
3
8
3
5
7
2
3
7
9
5
9
2
4
9
7
2
7
7
3
7
8
0
3
1
4
0
9
9
1
4
9
6
(a)
(b)
Mask: MNIST
Mask: KMNIST
Figure 4: A qualitative comparison between (a) Figure 5: A qualitative demonstration of the FC- based model trained using LSM on the data im- putation task.
our model (NIL=728) and (b) the model in [1 (NII==16337) on the inverse generation task.
cesses with 30 DF [xx [Ipu] and
2. Dr [px [po] = 20 Ep,, (=x) 12 8 log( Pxy, (mx)
J 1x.17 ||
(mxg
(i.e., Lemma A.13) by depicting the norm of their unbiased estimators 80 CssM (0) calculated with
and without applying the MaP technique in Fig. 3 is observed that the magnitude of ae (SsM(
significantly decreases when MaP is incorporated into the training process. This could be attributed to
the fact that the calculation of ao Dr Pxx, [Ip]] excludes the calculation of II Jg. in an Dr [Px IPo). which involves computing the derivatives of the numerically sensitive logit pre-processing layer.
Comparison with Related Works. Table 4 compares the performance our method with number of related works the MNIST dataset. Our models trained with score-matching objectives using the same model architecture exhibit improved perfor- mance comparison to the relative gradi- ent method (7]. In addition, when com- pared the results in [16] and [17), our models deliver significantly improved perfor- mance over them. Please note that the results of 7 16 17 presented in Table 4 are ob- tained from their original papers.
Table 4: A comparison of performance and training complexity between EBFlow and number of related works 16 17 on the MNIST dataset.
Method
Complexity
NLL ())
Baseline (ML)
O(D" L)
1092.4 + 0.1
DKL. Based
EBFlow (SML)
L(D*))
1092.3 + 0.6
Relative Grad. t
((DFL)
1375.2 + 1.4
EBFlow (SSM)
O(D)))
1092.8 + 0.3
EBFlow (DSM)
L(D*L)
1099.2 + 0.2
Dy-Based
EBFlow (FDSSM)
((DFL)
1104.1 + 0.5
SSM
3355
DSM 17
:
3398 + 1343
FDSSM 17
1647 +: 306
5.4 Application to Generation Tasks
The sampling process of EBFlow can be accomplished through the inverse function or an MCMC process. The former is typical generation method adopted by flow-based models, while the latter is a more flexible sampling process that allows conditional generation without re-training the model. In the following paragraphs, we provide detailed explanations and visualized results of these tasks.
Inverse Generation. One benefit of flow-based models is that g can be directly adopted as gener- ator. While inverting the weight matrices in linear transformations typically demands time complexity of O(D* these inverse matrices are only required to be computed once 9 has converged, and can then be reused for subsequent inferences In this experiment, we adopt the Glow 4 model architecture and train using our method with LSSM on the MNIST dataset. We compare our visualized results with the current best flow-based model trained using the score matching objective (17). The results of 177 are generated using their officially released code with their best setup (i.e., FDSSM). As presented in Fig. 4 the results generated using model demonstrate significantly better visual quality than those of [77.
MCMC Generation. In comparison to the inverse generation method, the MCMC sampling process is more suitable for conditional generation tasks such as data imputation due to its flexibility [4I] For the imputation task, data vector is separated IM. The goal of imputation is to generate the masked part æ based on the observable part æo- To achieve this goal, one can perform Langevin MCMC process update æM according to the gradient of the energy function 2 E(x; 0). Given noise vector sampled from (0, I) and Ow
observable part zo and masked part

---

small step size the process iteratively updates aM based on the following equation:
(++1) M
(1) M
Fl a ():0) + V2az,
= x
(11)
(t) dæ M
where æ (2) represents aM at iteration E {1, (1....,.7. T, and i' the total number of iterations. MCMC generation requires overall cost of O(TD²L), potentially more economical than the O L) computation of the inverse generation method. Fig. deppccs the imputation results of the FC-based model trained using LDSM the CelebA [42] dataset (D = 3 x 64 64). In this example,
implement the masking part æM using the data from the KMNNST 43 and MNIST 19 datasets.
Conclusion
In this paper, we presented EBFlow, new flow-based modeling approach that associates the parameterization of flow-based and energy-basee bassed models. We showed that by optimizing EBFlow with score-matching objectives, the computation of Jacobian determinants for linear transformations can be bypassed, resulting in an improved training time complexity. In addition, we demonstrated that the training stability and performance can be effectively enhanced through the MaP and EMA techniques. Based on the improvements in both theoretical time complexity and empirical performance, our method exhibits superior training efficiency compared to maximum likelihood training.
Acknowledgement
The authors gratefully acknowledge the support from the National Science and Technology Council well as the financial
(NSTC) in Taiwan under grant number MOST 111 11-22223-E0007-004-MM support from MediaTTek Inc., Taiwan. The authors would also like to express their appreciation for the donation of the GPUs from NVIDIA Corporation and NVIDIA AI Technology Center (NVAITC) used in this work. Furthermore, the authors extend their gratitude to the National Center for High- Perforreffr Computing (NCHC) for providing the necessary computational and storage resources.
References
[1] L. Dinh, D. Krueger, and Y. Bengio. NICE: Non-linear Independent Components Estimation, 2015.
[2] L. Dinh, J. N. Sohl-Dickstein, and S. Bengio. Density estimation using Real NVP. In Proc. Int. Conf. on Learning Representations (ICLR), 2016.
[3] G. Papamakarios, I. Murray, and T. Pavlakou. Masked Autoregressive Flow for Density Estimation. In Proc. Conf. on Neural Information Processing Systems (NeurIPS), 2017.
[4] D. Kingma and P. Dhariwal. Glow: Generative Flow with Invertible 1x1 Convolutions. In Proc. Conf. Neural Information Processing Systems (NeurIPS). 2018.
[5] C. Durkan, A. Bekasov, I. Murray, and G. Papamakarios. Neural Spline Flows. In Proc. Conf. on Neural Information Processing Systems (NeurIPS), 2019.
[6] A. Hyvärinen and E. Oja. Independent Component Analysis: Algorithms and Applications. Neural Networks: the Official Journal of the International Neural Network Sociery, 13 4-5:411- 30, 2000.
[7] L. Gresele, G. Fissore, A. Javaloy, B. Schölkopf, and A. Hyvärinen. Relative Gradient Op- timization of the Jacobian Term in Unsupervised Deep Learning. In Proc. Conf. on Neural Information Processing Systems (NeurIPS), 2020.
[8] Y. Song, C. Meng, and S. Ermon. MintNet: Building Invertible Neural Networks with Masked Convolutions. In Proc. Conf. on Neural Information Processing Systems (NeurIPS), 2019.
[9] E. Hoogeboom, R. v. d. Berg, and M. Welling. Emerging Convolutions for Generative Normal- izing Flows. In Proc. Int. Conf. on Machine Learning (ICML), 2019.
10

---

[10] X. Ma and E. H. Hovy. MaCow: Masked Convolutional Generative Flow. In Proc. Conf. on Neural Information Processing Systems (NeurIPS), 2019.
[11] Y. Lu and B. Huang. Woodbury Transformations for Deep Generative Flows. In Proc. Conf. on Neural Information Processing Systems (NeurIPS), 2020.
[12] C. Meng, I.. Zhou, K. Choi, T. Dao, and S. Ermon. ButterfyyFlow: Building Invertible Layers with Butterfly Matrices. In Proc. Int. Conf. on Machine Learning (ICML), 2022.
[13] Y. LeCun, S. Chopra, R. Hadsell, A. Ranzato, and F, J. Huang. A Tutorial Energy-Based Learning. 2006.
[14] A. Hyvärinen. Estimation of Non-Normalized Statistical Models by Score Matching. Journal of Machine Learning Research (JMLR), 6(24):695-709 2005.
[15] P. Vincent. A Connection between Score Matching and Denoising Autoencoders. Neural computation, 23(7):1661-1674, 2011.
[16] Y. Song, S. Garg, J. Shi, and S. Ermon. Sliced Score Matching: A Scalable Approach to Density and Score Estimation. In Proc. Conf. Uncertainty in Arlificial Intelligence (UAT), 2019.
[17] T. Pang. K. Xu, C. Li, Y. Song, S. Ermon, and J. Zhu. Efficient Learning of Generative Models via Diffe-eDifeeeece Score Matching. Proc. Conf. on Neural Information Processing Systems (NeurIPS), 2020.
[18] S. Lyu. Interpretation and Generalization of Score Matching. In Proc. Conf. on Uncertainty in Artificial Intelligence (UAI), 2009.
[19] L. Deng. The MNIST Database of Handwritten Digit Images for Machine Learning Research. IEEE Signal Processing Magazine, 2916):141-142, 2012.
[20] G. Papamakarios, E. T. Nalisnick, D. J. Rezende, S. Mohamed, and B. Lakshminarayanan. Normalizing Flows for Probabilistic Modeling and Inference. Journal of Machine Learning Research (JMLR), 22:57:1-57:64, 2019.
[21] G. Roberts and R. L. Tweedie. Exponential Convergence of Langevin Distributions and Their Discrete Approximations. Bernoulli, 2(4):341 363,, 1996.
[22] G. O Roberts and J. S Rosenthal. Optimal Scaling of Discrete Approximations to Langevin Dif- fusions. Journal of the Royal Statistical Sociery: Series (Statistical Methodology). 60(1):255- 268, 1998.
[23] G. E. Hinton. Training Products of Experts by Minimizing Contrastive Divergence. Neural Computation, 14:1771-1800, 2002.
[24] Tieeleman. Training restricted Boltzmann machines using approximations to the likelihood gradient. In Proc. Int. Conf. on Machine Learning (ICML), 2008.
[25] Y. Du, Li, B. J. Tenenbaum, and I. Mordatch. Improved Contrastive Divergence Training of Energy Based Models. In Proc. Int. Conf. on Machine Learning (ICML). 2021.
[26] M. F Hutchinson. A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines. Communications in Statistics-Simulatin and Computation, 18(3): 1059- 1076, 1989.
[27] A. Neumaier. Introduction to Numerical Analysis. Cambridge University Press, 2001.
[28] E. Parzen. On Estimation of Probability Density Function and Mode. Annals of Mathematical Statistics, 33:1065-1076, 1962.
[29] J. M. Tomczak and M. Welling. Improving Variational Auto-Encoders using Houscholder Flow. ArXiv, abs/161 1.09630, 2016.
[30] T. Dao, A. Gu, M. Eichhorn, A. Rudra, and C. Ré. Learning Fast Algorithms for Linear Transforms Using Butterfly Factorizations. In Proc. Int. Conf. on Machine Learning (ICML., 2019.
11

---

[31] J. Behrmann, D. K. Duvenaud, and J.-H. Jacobsen. Invertible Residual Networks. In Proc. Int. Conf. on Machine Learning (ICML), 2018.
[32] R. T. Q. Chen, Behrmann, D. K. Duvenaud, and J.-H. Jacobsen. Residual Flows for Invertible Generative Modeling. In Proc. Conf. on Neural Information Processing Systems (NeurIPS), 2019.
[33] R. M. Neal. Annealed Importance Sampling. Statistics and Computing, 11:125-139, 1998.
[34] Y. Burda, R. B. Grosse, and R. Salakhutdinov. Accurate and Conservative Estimates of MRF Log-Likelihood using Reverse Annealing. volume abss|441..856, 2015.
[35] W. Grathwohl, R. T. Q. Chen, J. Bettencourt, Sutskever, and D. K. Duvenaud. FFJORD: Free-form Continuous Dynamics for Scalable Reversible Generative Models. In Int. Conf. on Learning Representations (ICLR), 2018.
[36] Y. Song and S. Ermon. Improved Techniques for Training Score-Based Generative Models. In Proc. Conf. Neural Information Processing Systems (NeurIPS), 2020.
[37] A. Krizhevsky. Learning Multiple Layers of Features from Tiny Images. 2009.
[38] A. V. D. Oord, N. Kalchbrenner, L. Espeholt, K. Kavukcuoglu, O. Vinyals, and A. Graves. Conditional Image Generation with PixeICNN Decoders. In Proc. Conf. on Neural Information Processing Systems (NeurIPS), 2016.
[39] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, I.. Antiga, A. Desmaison, A. Köpf, E. Yang, Z. Vito,, M. Raison, A. Tejani, S. Chintala. PyTorch: An Imperative Style, High-Performant Deep Learning Library. In Proc. Conf. on Neural Information Processing Systems (NeurIPS), 2019.
S. Chilamkurthy, B. Steiner, L. Fang. J. Bai, and
[40] A. Griewank and A. Walther. Evaluating Derivatives- Principles and Techniques of Algorithmic Differentiation, Second Edition. In Frontiers in applied mathematics, 2000.
[41] A. M Nguyen, J. Clune, Y. Bengio, A. Dosovitskiy, and J. Yosinski. Plug & Play Generative Networks: Conditional Iterative Generation of Images in Latent Space. In Proc. Int. Conf. Computer Vision and Pattern Recognition (CVPR), 2016.
[42] Z. Liu, P. Luo, X. Wang, and X. Tang. Deep Learning Face Attributes the Wild. In Proc. Int. Conf. Computer Vision (ICCV), December 2015.
[43] T. Clanuwat, M. Bober-Irizar, A. Kitamoto, A. Lamb, K. Yamamoto, and D. Ha. Deep Learning for Classical Japanese Literature. In Proc. Conf. on Neural Information Processing Systems (NeurIPS), 2018
[44] T. Anderson Keller, Jom W. T. Peters, Priyank Jaini, Emiel Hoogeboom, Patrick Forr'e, and Max Welling. Self Normalizing Flows. In Proc. Int. Conf. on Machine Learning (ICML), 2020.
[45] C.-H. Chao, W.-F. Sun, B.-W. Cheng, Y.C. Lo, C.-C Chang, Y.-L. Liu, Y.-. Chang, C.-P. Chen, and C.-Y. Lee. Denoising Likelihood Score Matching for Conditional Score-based Data Generation. In Proc. Int. Conf. on Learning Representations (ICLR), 2022.
[46] Diederik P. Kingma and Jimmy Ba. Adam: A Method for Stochastic Optimization. CoRR, abs/1412.699 2014.
[47] I. Loshchilov and F. Hutter. Decoupled Weight Decay Regularization. In Proc. Int. Conf. Learning Representations (ICLR). 2017.
[48] M. Ning, F. Sangineto, A. Porrello, S. Calderara, and R. Cucchiara. Input Perturbation Reduces Exposure Bias in Diffusion Models. In Proc. Int. Conf. Machine Learning (ICML), 2023.
12

---

[49] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, Goodfellow, A. Harp. G. Irving, M. Isard, Y. Jia, R. Joze- fowicz, I. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, Tucker, V. Vanhoucke, V. Va- sudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng. TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems, 2015. Software available from tensorflow.org.
[50] J. Sohl-Dickstein Two Equalities Expressing the Determinant of a Matrix in terms of Expecta- tions over Matrix-Vector Products, 2020.
[51] K. P. Murphy. Probabilistic Machine Learning: Advanced Topics, page 811. MIT Press, 2023.
[52] M. Zhang, O. Key, P. Hayes, D. Barber, B. Paige, and F.-X. Briol. Towards Healing the Blindness of Score Matching. Workshop on Score-Based Methods at Conf. on Neural Information Processing Systems (NeurIPS), 2022.
13

---

A Appendix
A.1 Derivations
In the following subsections, we provide theoretical derivations. In Section A.1.1 we discuss the asymptotic convergence properties as well the assumptions of score-matching methods. In Section A.1.2. we elaborate on the formulation of EBFlow (i.e., Eqs. (8) and (9), and provide explanation of their interpretation. Finally, in Section A.1.3 present theoretical analysis of KL divergence and Fisher divergence, and discuss the underlying mechanism behind the proposed MaP technique.
A.1.1 Asymptotic Convergence Property of Score Matching
In this subsection, we provide a formal description of the consistency property of score match- ing. The description follows [16] and the notations are replaced with those used in this paper. The regularity conditions for p(:;6) are defined in Assumptions |.
A.7 In the following para- Eax; 0)
graph, the parameter space is defined as O. In addition, s(æ;; 0) Ow logp(z;; ()
represents the score function. (0) A N
ΣN [xxx; 0) denotes an unbiased estimator of F(x;0) 0) = 0)|1 + Tr ( <(x; 0))
(0), where f(;; 0) - 2 and (xI,... .Xv) represents a collection of i.i.d. samples drawn from px- For notational simplicity, Ba h(x; 0) and (x; 0) A @
21 F(x;; 0)| Tr
0²
we denote @h(x; 0)
oz. (,;;;;;) 0), where (a(3; 0) denotes the j-th
element of h.
Assumption A.1. (Positiveness) p(a; 0) 0 and px(x)) 0, 40 O, Vo ER.
Assumption A.2. (Regularity of the score functions) The parameterized score function s(x; 0) and the true score function Fe log Px((x) are both continuous and differentiable. In addition. their expectations (m)) [s(æ;;) 6)] and (a) de log (x)))( finite. (i.e., (x) [s(x; and (x)( 80 logPx Px (n:] < o0)
(x)))(( 0) 4O E Assumption A.3. (Boundary condition) -
Assumption A.4. (Compactness) The parameter space is compact.
Assumption A.5. (Identifiability) There exists set of parameters 6* such that px(x)) = p(a: 8*),
where 8" E θ, Vo E R...
Assumption A.6. (Uniqueness) 0 + 0* + p(æ; 0) + p(x; 0*), where 0,0 0" O
Assumption A.7. (Lipschitzness of f) The function f is Lipschitz continuous w.r.t. 0, i.c., [[ 01) - f(:; < L(x)) |0 02112- 10,, O2 E O, where L(x) represents a Lipschitz con- stant satisfying (p-(#) [L(x]] < cc.
Theorem A.8. (Consistency of a score-matching estimator ON A argmingge is consistent, i.e.
[6] The score-matching estimator
On 2 0", as N
Assumptions A.1
A.3
the conditions that ensure Dr [px (ar) (p))) 0)] LsM (0). Assump-
tions A.4 A.7 lead to the uniform convergence property [16] of a score-matching estimator, which gives rise the consistency property. The detailed derivation can be found in Corollary (16). In the following Lemma A.9 and Proposition A. A. satisfy Assumption A.7
we examine the sufficient condition for and Pa to
Lemma A.9. (Sufficient condition for the Lipschitzness of f) The function f(x; = [1<(;; 0)| Tr De s(x; 0)) is Lipschitz continuous if the score function s(x; 0) satisfies the following conditions: vo, 0, 02 E Vi (I, (.,
||s(;; 0) 2 < (x))
| 01) s(2; < (œ) ||0 - 0z
(0,,(<r; 01) 0,s((r; 02)111 < (x)))) ||0 0g 0211g:
where L1, Ly, and Lg Lipschitz constants satisfying Epa(>) [Z]( Ep)()) [2(x:)] < co. and (x)( [La(x']
14

---

Proof. The Lipschitzness of f can be guaranteed by ensuring the Lipschitzness of ||s( 0) and Tr ((s((;;)))
Step (Lipschitzness of [s((2))11))
[| 01)112 - [|s( 12)
ss(;; (1)0)1))0000 s(;; 01) s(x; 02) (z) s(x;
(s(x;; (1) 01) - s(x; 81) )d s(x: (2)) (1)=;;)))s((;; B2) - s(x; =2x)=(<;;())))))))))))))<)
ss(x;; 01)² (s(z; 01) s(x: 02)) + (z: (s(z:; 01) - s(x; 82))
01) (s( s(x; 02) + ( 02) (s(x; 01) s(x; 02))
[| 01 )l1 [|<(2; 01) s(=; 0z) || + || (2;; 02) 112 [1= 01) s( 02)112
1,1 (a:) [s<(;; 01) 02x;02)11 + T (x)) (a: 01) - <2; 02)1l2
66) 2L L2) 10 - 02112
where (i) is based on triangle inequality, (ii) is due to Sauchy-Schwar inequality, and (itt) follows from the listed assumptions.
Step 2. (Lipschitzness of Tr (0s(x;)))
Tr (0s(2;00)))) Tr (0s(2; 02)) = |Tr (0s(2;;0) 01) as(a;; 02))
D [10s(<;;) 01) - 0s(z; 02)112
(40)
D
10:s(1:; 01) (12(1))1
(68) S Dy 2L 110 - 02112
DV DVDL (a:) 110 02l12
where (i) holds by Von Neumann's 's trace inequality. (di) is due to the property ||A|||2 S v V.,, la. where a is the column vector of A. (iti) holds by the listed assumptions. Based on Steps 1 and 2, the Lipschitzness of f is guaranteed, since
[s 01) - f(x; 02) =
2 |s( 0,)11 + Th
a
(ax; 01)
2 [1((000))1+
- TT
(x;; 02)
dx
2 ||( 0 0011 - 2 |x 02)11
+ 'T
s(x; s(x;01) 01)
Tr
s (x; 02)
Ox
Ox
S ||<(;; 0, 2112 || 02) 02)
+ +r
Oa: s2;; 01)
Tr
sax; 02)
Fx.
s L1 11: - 02112 + DVDLa(x)) 110 02112
(2)()) (L( DVDL 110 82ll2.
Proposition A.10. (Sufficient condition for the Lipschitzness of f) The function f is Lipschitz continuous if (x:; 8) has bounded first, second, and ohird-order derivatives, Le,, Vi,j E {1,-.. D), v0 E O.
(x; a 0)1, S (1x), 11, 0)1, < 11000, 0)ll, S lz(x), and smooth enough on O i.e., 01, O:
[1 02) ((x; 02) rx) 110 0z l
15

---

[x 01) Jg(x; (2x=;;2)1 02)112 5 (=) |0] - 82112
g 01) (2J,(2;;))] S r2()) ||0 8211.
| o (m; 01) (, (2:00)))) 02)11 S (r(2)) 02
In addition, it satisfies the following conditions:
J J [2; 0) - ((())) | |
(o: 0)1 s (2(1).
- 02)11 112 S (w) 10: 02112
||aJ, 01) a.J, (a; 02)112 S (r((a)) 1|0 02l12
where J represents the inverse matrix of g.- Furthermore, the prior distribution Pu satisfies:
||su(u)]] - t,, [10;su(())) tz ((n(())) su (22)))2 I La 1| u2 10:s((u)) (12)( 12 S La |lue u2 12:
where 8,(()) A log (w)) is the score function of Pu- The Lipschitz constants listed above (i.e., Tye l 1a. 13, l l22 and ri r22) have finite expectations.
Proof. We show that the sufficient conditions stated in Lemma A.9 can be satisfied using the conditions listed above.
Step (Sufficient condition of [|<(;; 0) L1 (xx))
Since [|s(;;; 0)12
2 log PogPa(g(a (z;; 0))
log det Jy (a; demonstrate
0)]l
<
dee log PogPoog(a0 (p( (1)12)
aee
log |dct Jy (z: 0)
we first
that
a
and log (u()((;)))l
Be
2
dm
Ø de
log |det (x; 0) are both bounded.
(1.1)
d logpa(g(>; a;; 0)) 2 is bounded:
log Pa x: 0)) (s(()(z;0))*
S [su(g(=;;)))) [I, S tt)((x).
dæ
2
(1.2) bn
log det (x; |||| is bounded:
(x:
a
log |et (x; 0)
det (x; 0)|
det J (x;
Zx
(del J. (2x; 0))
det J (z; 0)
dx
(det. J, (2;; 0)) del J, (x; (); 0)
[ 0)ll
where (i) derived using Jacobi's formula, and (x;; 0) (J (x: 0)0,J, 0)).
|< 0)1 =
Tr 0)a,), (z)
Σ D2 J (x; 0),J (x;
D2 || (x: 0>| |1a. (x; 0)11
(H6)
(x1((1))
= (=)((()))(
16

---

where (i) holds by Von Neumann's 's trace inequality, (ii) is due to the property of matrix norm, and (iti) is follows from the listed assumptions.
Step 2. (Sufficient condition of the Lipschitzness of s(x;; 0))
Since s(x; 8) dæ log Pu (u( 8)) + dæ log |dct Jy(::) 8)), we demonstrate that On log Pu (a(g(a; (x; 8)) 8 log det J, (2; 0) are both Lipschitz continuous on O
and 82
(2.1) Lipschitzness of de log log
a
log (u 81))
dæ logpu(g(z:
Ox
(sulg(<; 81))) J (x; (su(g(æ; 02)))² (z 82)
|| I| (z: (1) (z; 02)112 + ||su 01)) (2))l1; (T 02)
4i (m)) ||0: 02l12 124 (z:) [ 01) ((=:02)]1
(8)
2 (x) ||0: - 0211. + 124 (2)) 110 - 02ll2
(t1ri (x) (z)rooo) |01 82ll2
where (i) obtained using similar derivation Step I in Lemma A.9 while (ii) follows from the listed assumptions.
(2.2) Lipschitzness of am log det (z;
Let a( _
0). We first demonstrate that M is Lipschitz continuous:
(M, M
-|1 (1:
(a; (m; (0m 02)112
$ |3: [10.J 01) J.J, 02)) 12 + [J - -J (xx; 02) |J,J (x:; 02)
(0)
((0))r=(=) 0211, -2 (2)()rr 110 02l12
(4 (((0))((0))++())))i) 110 02l12
where (i) obtained by an analogous derivation of the step in Lemma A.9 and (ii) holds by the listed assumption.
The Lipschitzness of M leads to the Lipschitzness
ön log |dct. Jy (z; 0) since:
log /det J (x; 01)1 Zi log det 02)
dr
|| 01) v(x; 02)112
(Tr(M((;a;))) Tr (M(;, 02)))
(Tr (MM ( a(;, 0z)))²
Σ Da |M((,,;;,) Mi a;;
Σ D2 ([(()))(()++(())))))))) 2 10, - 0211
vD (((x))2(() + (())r((()) 110 02l12
17

---

where (1) holds by Von Neumann's 's trace inequality, (ži) is due to the Lipschitzness of M.
Step 3. (Sufficient condition of the Lipschitzness of 0;s(2;; 0))
(æ; 0) can be decomposed as (0,su(g(<;;)) 0))) (z;00)), (su(g(a;())) 0,J, (x; and 0, [0 0] as follows:
a, 0) 0, (su(o(<;; 0))) (z 0)
+ 0 [v(2;;)) 0)]
(0:su(9(2;; 0))) 0)]
[su(o(<: 0))) (2; 0)] a [v(x;;
(3.1) The Lipschitzness of (a;; 0))) " Jn (2:; 0) and (su(g(æ; 0))) a.J, (x; 0) can be derived using proofs similar to that in Step 2.1:
(0,sulg(2; 0())) J,(2; 01) - (D.sulg(a; 02))) Jg (x: 02)
- (12r1 +40l)))) 10 - 02112
(su(o(2;0))) (1)) (a3 - (su(g(a;;))) 02)))) , (z; 02)
- ( ((rr0(a))+ 1 - 02l12
(3.2) Lipschitzness of (v( 0)]:
Let 0: (x; 0)] A Tr (M(j x; = Tr (8.M(, x;; We first show that a,M(j, x;; 0) can be decomposed as:
0,M(j. x: 0) à (J x [ = (a.J, (a:J (z:
+ (J (x: (
The Lipschitz constant of 0 equals to (/2 (2)) + (())))((()) + (4(x))a((x) la)z)r
based on a similar derivation as in Step 3.1. The Lipschitzness of 8 M, x;; 0) leads to the Lipschitz- ness of 0 (a;; 0)]:
||0: (a;; 01)] - 0: (a; 02)1112
Tr (0.M( - Tv
Σ Tr (8)M(,))) - a;; 02))
D2 [0.M(j x
D2
(a(())((())) 110 - 02)1l²
D (((z))2(z) + +)(0)((()) + +())))()))++(())((())))+++++)+)) 101 - 02)1l2
where (1) holds by Von Neumann's trace inequality, (ti) is due to the Lipschitzness of O.M.
A.1.2 Derivation of Eqs. and
Energy-based models are formulated based on the observation that any continuous pdf p(æ; 0) can be expressed a Boltzmann distribution (xp 0)) --) [I3, where the energy function E(:;0) can be modeled as any scalar-valued continuous function. In EBFlow, the en- ergy function E (z; 0) is selected as - log(pu (g(x; 0)) EI, (det(J. (a,-1;0)) according
Eq.
This suggests that the normalizing constant Z(Ø) = f exp -((x;; 0)) dx is equal to (det(Jg. 9. (0)))) according to Lemma A.11
gg
Lemma A.11.
II det (o,(J) g.E..
[ S Pu 0))
II
[det(J,, (x)-1;0)) dæ.
(AI)
g.ES,
18

---

Proof.
[det(Jg, (xx-1; dx. plæ; 0)d
WERD
Pu (g(x;
|ct(JJg gy (x.-1:0))
II S.ES;
Idet (0)) dx
ZERD
S.ES,
/
[det(Jo 9. (0))
Pu (g(x;
II 9:ES,
[det(Jg. g. (x.-1:8)) dæ
EERD
si ES
By multiplying (Iorss, (d)) -1
both sides equation, we arrive the conclusion:
II get (x)))) -1
Pu (g(z; 9i ES
II
ZERD
9: ES
point-wise donsity evaluation
(x,
Px(x) Pa (x) M
(data)
PxL
Px,
Px
9L
95+1
9)
91
XL
Xj
X
(madel)
P
Pj
po
piint-wise
density evalistion
Pu(a))
(y(n)) = Pulge ( det(J,)))
2-311
KL divergence
DKL. Px. lpz)
Ds (Px, [Ips]
DRL [Px [o]
divergence (Lemmna A31
De [Px]lpp]]
general)
Dy Dp[[Ps,llps)
(in general)
DF[px| Po]
Figure A1: An illustration of the relationship between the variables discussed in Proposition 4.1 Lemma 1.12 and Lemma A.13 x represents a random vector sampled from the data distribution {g: is series of transformations. X, o ...... o g1 (x), and Px, is its pdf. py(m)) = PulgL g)+1(0,)) II /det (J,,)l, where Pu is a prior distribution. The properties of KL divergence and Fisher divergence presented in the last two rows derived in Lemmas A.12
Px-
and A.13
A.1.3 Theoretical Analyses of KL Divergence and Fisher Divergence
In this section, provide formal derivations for Proposition 4.1 Lemmal 12 and Lemma A.13
To ensure a clear presentation, we provide a visualization of the relationship between the variables used in the subsequent derivations in Fig. AT
Lemma A.12. Let Px, be the pdf of the latent variable xf & gjo goo---ogi (x) indexed by j. In addition,
let (.(.) be pdf modeled as Pulgi o ... ()+1(-)) [[= Idet (J,.), where j E {0,-.. L - 1}.
It follows that:
DKL [px, [pps] DEL (pxl1po) Vj E (1,... -- 1}.
(A2)
19

---

Proof. The equivalence DKL [px lxpllol = DKL [Px,, [|pb]
holds for any j fl,... L 1 since:
DKL [px Ipol
Pxx
(p)
log
(o(x)
Pxx (9j 0...... o (2))
i=1 det (Jg:)]
(x))
log
Pu o............................... (x))
1= det (J
Px, Px, ..... ...01.11000 g1()))))
log
Pu(gL 0.....0......................... (1())))
I1 det (J (Jor))
Px, (2t)
Bpn, (1z)
log
Pullr o... 9j+1 a (z))) I1 dee (J
DKL [px, [lps]:
where (i) is due to the property that (w) o gj o.... gi o....g.
Bpn, (a,)) for given
function f. Therefore, DKL [px, [pps] DEL [xxllpo]. Vj E {1..... , 1-
Lemma A.13. Let be the pdf the the latent variable of gjo...oo .ogi (x) indexed by j. In addition,
let (.(.) be pdf modeled as Pulg1 o (j+t I |det (Jg.), where j E (0, (0,... I - 1}. follows that:
1
0
Px, (x,) (T{))
DF [pxllpo) = Epp, (ars)
log
Jn.
vj e (1, ... ,.,,,.,------- 11. (A3)
2
Proof. Based on the definition, the Fisher divergence between px and po is written
De [px]lPo]
θ
Px (x(
(pp
log
2
Ox
(x))
Px, (gg o... o 91 (1(1))))))))) (x)) IT
dee (Jg.]
(pp(
log
Oa:
Pulgi o ..................... (x)
dee (J,.).
0
Px, Px; ....o 00....i.6.... gx()))
= (xe)
log
2
Ox
Pulge
o... gi (x:) [ det (Jg:))
a
Px, (gj 9n (x))
Og;
91 (x)
Epxe (#)
log
0gj ... o (x)
Pu(gL ..... ( II [7 |det (J,.))
dx
(
Px, (g o.................... (xx)
E Epe (20)
log
J g1
2 Oyg o o... (x))
Pulgl o... )z IT ldet (J
(i8)
Px x
Epx, (1ss)
log Dx.
J 9.
2
Pu(gL o .... go+1 (xx
I1 det (J₂.))
1=1
a
Pxx (2xz) (T(
= Bpx, (m1)
log
②
g: 1=1
where (i) is due to the chain rule, and (ii) is because Ep If o o o.... o...))(((x)) = Ep,, (e))) (xx,) for given function f.
0
20

---

Remark A.14. Lemma A.13
implies that Dp [px, [|pp]
* De [px][Po]] in general, as the latter contains
an additional multiplier Io as shown below:
Px, (zy) (z({z))
Dp [px||po] = Byy, (se,)
log
J g:
2
dxx
1
8
Pxx (xx((0s) py(2T)
De [Px, [1pp] = Byx, (oes)
log
2
J
Proposition 4.1. Let Px; be the pdf of the latent variable of Xj - 9j (x(x) indexed by
j. In addition, let (.).) be a pdf modeled as Pulgi o ...... 9 95-1 (.)) II= det (Jg,)l where
j E (0, -0. , 1-. It follows that:
Dp [xx, [pps] = =0 DF INx [po] = 0,Vi E (1,
(A4)
Proof. Based on Remark A. 14 the following holds:
Px, (z,)
Dr
Pxx
|py]
Pay (ze,)
log
=
Jæ.
i
Px, (xx,)
log
0
Fzy
Pi (2tz)
i
Prxy (1zz)
log
Jg,
(2;;)
2 Dr [pxl
1
θ
log Px3 (ax)) Py(m,)
Bapsy (oes)
J
0,
2
Ox,
where (i) and (ii) both result from the positiveness condition presented Assumption A.1 Specif-
ically, for (i), Px; (xzz) Px(9
o 9 (x;)) det. (J O, since px and
Ih det J | = -h= /det (Jg:)) 1=)1 > 0. Meanwhile (it) holds since IT /eet (Jg.)l
and thus all of the singular values of IT Jg. are non-zero.
A.2 Experimental Setups
In this section, elaborate on the experimental setups and provide the detailed configurations for the experiments presented in Section of of the main manuscript. The code implementation for the experiments is provided in the following repository: https: wig.t.u.ccom. chen hao- chao/ ebflow Our code implementation is developed based on
A.2.1 Experimental Setups for the Two-Dimensional Synthetic Datasets
Datasets. In Section 5.1 we present the experimental results three two-dimensional synthetic datasets: Sine, Swirl, and Checkerboard. The Sine dataset is generated by sampling data points from the set {(4w 2, sin( (12w 6))
[0, 1]}. The Swirl dataset is generated by sampling data points from the set { - ग / " cos ( ww). . ग / sin /i)) w E [o. 1]}. The Checkerboard dataset generated by sampling data points from the set {(4w 2, t 2s 44w 2] mod 2) [o, 1]. t
[o, 1]. s E {0, 1}}, where [.] is a floor function, and mod represents the modulo operation.
To establish px for all three datasets, we smooth Dirac function using Gaussian kernel. Specifically, define the Dirac function f()) & 1 M E 6(][ =() [|). where (=(() MM are M uniformly-
M
sampled data points. The data distribution is defined as Px & I [() 8² o @ The closed-form expressions for (x: and Dur log Px (xx)) can be ob-
M M (() tained using the derivation in [4). In the experiments, M is set as 50, 000, and fixed at 0.375 for all three datasets.
21

---

Implementation Details. The model architecture of g(:; ; ;0) consists of ten Glow blocks 4]. Each block comprises an actnorm 4 layer, a fully-com fully-connected layer, and an affine coupling layer. Table A2 provides the formal definitions of these operations. Pu(-) is implemented as an isotropic Gaussian with zero mean and unit variance. To determine the best hyperparameters, we perform grid search over the following optimizers, learning rates, and gradient clipping values based on the evaluation results in terms of the KL divergence. The optimizers include Adam [46], AdamW [47], and RMSProp.
The learning rate and gradient clipping values are selected from (5e-3, le-3, 5e-4, le-4) and (None, 2.5, 10.0), respectively. Table AT] summarizes the selected hyperparameters. The optimization processes of Sine and Swirl datasets require 50,000 training iterations for convergence, while that of the Checkerboard dataset requires 100,000 iterations. The batch size is fixed at 5,000 for all setups.
A.2.2 Experimental Setups for the Real-world Datasets
The experiments presented in Section 5.2 are performed on the MNIST 19 and CIFAR- 10 (37] datasets. The training and test sets of MNIST and CIFAR-10 contain 50,000 and 10,000 images, respectively. The data are smoothed using the uniform dequantization method presented in M The observable parts (i.e., 20) of the images in Fig. [ are produced using the pre-trained
Datasets.
model in 48
Implementation Details. In Sections 5.2 and 5.4, we adopt three types of model architectures:
FC-based [7]. CNN-based, and Glow models. The FC-based model contains two fully-connected layers and a smoothed leaky ReLU non-linearity 7 in between, which is identical to 7 The CNN- based model consists of three convolutional blocks and two squeezing operations 2 between every convolutional block. Each convolutional block contains two convolutional layers and a smoothed leaky ReLU in between. The Glow model adopted in Section [3.4]is is composed of 16 Glow blocks. Each of the Glow block consists of an actnorm (4) layer, a convolutional layer, and an affine coupling layer. The squeezing operation is inserted between every eight blocks. The operations used in these models are summarized in Table A2 The smoothness factor a of Smooth Leaky ReLU is set to 0.3 and 0.6 for models trained on MNIST and CIFAR-10, respectively. The scaling and transition functions s(-: 6) and t(:: :)) of the affine coupling layers are convolutional blocks with ReLU activation functions. The prior distribution p.(()) is implemented as an isotropic Gaussian with zero mean and unit variance. The FC-based and CNN-based models are trained with RMSProp using learning rate initialized at le-4 and a batch size of 100. The Glow model is trained with an Adam optimizer using a learning rate initialized at le-4 and a batch size of 100. The gradient clipping value is set to 500 during the training for the Glow model. The learning rate scheduler MultiStepLR in PyTorch ch is used for gradually decreasing the learning rates. The { used in DSM and FDSSM are selected based on a grid search over {0.05, 0.1, 0.5, 1.0). The selected { are {1.0, 1.0) and f0.1, for the MNIST and CIFAR-10 datasets, respectively. implemented using PyTorch (39].
The parameter m in EMA is set to 0.999. The algorithms The gradients w.r.t x' and 8 are both calculated using automatic differential tools [40] provided by PyTorch (39. The runtime is evaluated on Tesla V100 NVIDIA GPUs. In the experiments performed on CIFAR-10 and CelebA using score-matching methods, the energy function (i.e., E. (x) (F(x; 0))) is added as regularization loss with balancing factor fixed at 0.001 during the optimization processes. The results in Fig. (b) are smoothed with the exponential moving average function used in Tensorboard [49 i.c., w d:-1 (1 w) x d,, where is set 0.45 and d, represents the evaluation result at the i-th iteration.
Table A1: The
used in the two-dimensional synthetic example in Section 5.1
Dataset
ML
SML
SSM
DSM
FDSSM
Optimizor
Adam
AdamW
Adam
Adam
Adam
Sine
Learning Rate
Se-4
Se-4
1c-4
lc-4
le4
Gradient Clip
1.0
None
1.0
1.0
1.0
Optimizor
Adam
Adam
Adam
Adam
Adam
Swirl
Learning Rate
5c-3
le-4
le-4
le-4
lc-4
Gradient Clip
None
10.0
10.0
10.0
2.5
Optimizor
AdamW
AdamW
AdamW le-4
AdamW le-4
Adam
Checkerboard
Learning Rate
le-4
le-4
le-4
Gradient Clip
10.0
10.0
10.0
10.0
10.0
22

---

Table A2: The components g(::0) used in this paper. In this table, and y are the output and the input of a layer, respectively. 3 and represent the mean and variance of an actnorm layer. w is convolutional kernel, and w - Wy, where is a convolutional operator, and W is D x D matrix. W and represent the weight and bias in a fully-connected layer. a is a hyper-parameter for adjusting the smoothness of smooth leaky ReLU. In the affine coupling layer, & and y are split into two parts { ] and {ya: 3o respectively. s(;:0) and t(:;0) are the scaling and transition networks parameterized with 8. sig (y) = 1/(1 + exp (-y)) represents the sigmoid function. dim (+)
represents the dimension of the input vector. y represents the i-th element of vector y.
Layer
Function
Log Jacohian Determinant
Set
actnorm
x (y 3)/7
Di= log 1 1//gal log det W log Idet (W)
S
convolutional
+ b
Si
fully-connected
2 Wy +
S
smooth leaky ReLU 7
z = ay + (1 log(1 + exp (#))
E log aa + (1 - as)sig (ma)
Sn
affine coupling 4
Za s s3ss: 0) + fyy:: 0). 2
Σ -==1
diant vs log (<vo: ra
Sn
Table A3: The simulation results of Eq. A6 The error rate measured by Idurue dest /d/uuel where drrue and dest represent the true and estimated Jacobian determinants, respectively.
=50
-100
D =200
Error Rate (M =50)
0.004211
0.099940
0.355314 14
Error Rate (M =100)
0.003503
0.034608
0.076239
Error Rate (M =200)
0.002332
0.015411
0.011175
Results of the Related Works. The results of the relative gradient 77. SSM [16], and FDSSM [17] methods are directly obtained from their original paper. On the other hand, the results of the DSM method is obtained from [17]. Please note that the reported results of [16] and [17] differ from each other given that they both adopt the NICE I model. Specifically, the SSM method achieves respectively. Moreover, the DSM method achieves NLL= 4, 363 and NLL= 3, 398 in and respectively. Table 4
NLL= 3, 355 and NLL= 6, 234 in [16] and 17
we report the results with
lower NLL.
A.3 Estimating the Jacobian Determinants using Importance Sampling
Importance sampling is a technique used to estimate integrals, which can be employed to approximate the normalizing constant Z(0) in energy-based model. In this method, a pdf q with simple closed form that can be easily sampled from is selected. The normalizing constant can then be expressed as the following formula:
z(0) = exp(-E(z;0)) dx /f
cxp (- -E(x: 0))
q(x)
dx
2EER
JEERD
((x))
Tepp
(-E(x; 0))
exp (-E(#0);
Em(())
g(x))
M
((x())
(A5)
where M represents M i.i.d. samples drawn from q. According to Lemma و
A.11
the Jacobian
determinants of the layers in St can be approximated using Eq. (A5) follows:
M Pu (g(* 01) IL, (Jt (a
II E.
[et (0))
23
(A6)
M
9 (
j=1
23

---

Table A4: An overall comparison between EBFlow, the baseline method, the Relative Gradient method 77. and the methods that utilize specially designed linear layers 18-12 29]. The notations /I X in row "Unbiased" represent whether the models optimized according to unbiased target. On the other hand, the notations v/ X in row "Unconstrained" represent whether the models can be constructed with arbitrary linear transformations. (7) The approximation errors o(E) FDSSM controlled by its Syper-parameter Ç. (1) The error o(W) of the Relative Gradient method is determined by the values of model's weights.
KL-Divergenee-Based
Fisher-Divergee BineerDDiveggnee--Be
Baseline (ML) EBFlow (SML) Relative Grad. Special Linear EBFlow (SSM) EBFlow (DSM) EBFlow (FDSSM)
Complexity
O(D'L
O(D*L)
OLD*L)
OID*L
O(DLL)
O(DLL)
O(D*L
Unhiased
X
X!I)
Unconstrained
To validate this idea, provide simple simulation with Pa NO0 1= N(,, f), g(;; W) Wæ, æ. .M = (50, 100, 2003, and D = {50, 100, 200} Table A3 The results show that larger values of .M lead to more accurate estimation of the Jacobian determinants. Typically, the choice of is crucial to the accuracy of importance sampling. To obtain an accurate approximation, one can adopt the technique of annealed importance sampling (AIS) [33 or Reverse AIS Estimator (RAISE) 34]. which are commonly-adoppe
algorithms for effectively estimating (0)).
Eq. A66 can be interpreted as a generalization of the stochastic estimator presented in [50 where the distributions Pu and 9 are modeled isotropic Gaussian distributions, and g is restricted as linear transformation. For the further analysis of this concept, particularly in the context of determinant estimation for matrices, we refer readers to Section of [50], where a more sophisticated approximation approach and the corresponding experimental findings are provided.
A.4 Comparison among the Methods Discussed in this Paper
In Sections
2 3 and we discuss various methods for efficiently training flow-based models.
To provide comprehensive comparison of these methods, summarize their complexity and characteristics in Table A4
A.5 The Impacts of the Constraint of Linear Transformations on the Performance of Flow-based Model
In this section, examine the impact of the constraints of linear transformations the performance of a flow-based model. key distinction between constrained and unconstrained linear layers lies
Leamable Weight
Masked Weight (Not Leamable)
X
Lower Triangular Matrix (L)
Upper Triangular Matrix (U)
Lower & Upper Trinngular Matrices (LU)
Full
Figure A2: An illustration of the weight matrices in the F, L, U, and LU layers described in Section A.5
X(4
Figure A3: Visualized marginal distributions Pxcj) for i = 1,2, 3, and 4.
24

---

in how they model the correlation between each element in a data vector. Constrained linear trans- formations, such as those used in the previous works [2 29 impose predetermined correlations that are not learnable during the optimization process. For instance, masked linear layers [8I00] are constructed by masking either the upper or lower triangular weight matrix in linear layer. In contrast, unconstrained linear layers have weight matrices that are fully learnable, making them more flexible than their constrained counterparts.
To demonstrate the influences of the constraint on the expressiveness of a model, we provide performance comparison between flow-based models constructed using different types of linear layers. Specifically, we compare the performance of the models constructed using linear layers with full matrices, lower triangular matrices, upper triangular matrices, and matrices that are the multiplication of both lower and upper triangular matrices. These four types of linear layers are hereafter denoted as F, I., U, and L.U, respectively, and the differences between them are depicted in Fig. A2 Furthermore, to highlight the performance discrepancy between these models, we construct the target distribution ps based on an autoregressive relationship of data vector x. Let Xij denote the i-th element of x, and Pxyy represent its associated pdf. X is constructed based on the following equation:
[
if 1, Lanh(ug x 8) x (XK-1) + d x 2'),, i (2...., D}.
[
(A7)
where is sampled from an isotropic Gaussian, and and are coefficients controlling the shape and distance between each mode, respectively. In Eq. (A7). p the function tanh(-) can be intuitively where H(-) represents the Heaviside step
viewed as a smoothed variant of the function (-) function. In this context, the values of (X + d x 2¹) are multiplied by a value close to either -1 or 1, effectively transforming a positive number a negative one. Fig. A3 depicts number of examples of Рxµa constructed using this method. By employing this approach to design Px, where capturing Px is presumed to more challenging than modeling Px,i for any j < í, we can inspect how the applied constraints impact performance. Inappropriately masking the linear layers, like the U-type layer, anticipated to result in degraded performance, similar to the anti-casual effect explained in [51).
In this experiment, we constructed flow-based models using the smoothed leakyReLU activa-
LU
70
tion and different types of linear layers (i.e., F. L, U, and LU) with dimensionality of D = 10. The models are optimized according to Eq. The performance of these models is evaluated in terms of NI.L, and its trends are depicted in Fig. A4 It is observed that the flow-based model built with the F-type layers achieved the lowest NLL, indicating the advantage of using unconstrained weight matrices in linear layers, In addition, there is noticeable performance discrepancy between models with the L-type and U-type layers, indicating that imposing inappro- priate constraints on linear layers may negatively affect the modeling abilities of flow-based mod- els. Furthermore, even when both L-type and U-type layers were adopted, as shown in the red curve in Fig. A4 the performance remains inferior to those using the F-type layers. This experimental evidence suggests that linear layers constructed based on matrix decomposition (e.g., [4] 9 may not possess the same expressiveness as unconstrained linear layers.
60
45
2,000 4,000
6,000
8,000
10,000
Step
Figure A4: The evaluation curves in terms of NLL of the flow-based models constructed with the F- type, L-type. U-type, and LU-type layers. The curves and shaded area depict the mean and 95% confidence interval of three independent runs.
A.6 Limitations and Discussions
We noticed that score-mater matching methods sometimes exhibit difficulty in differentiating the weights be- tween individual modes within multi-modd distribution. This deficiency is illustrated in Fig. A5 (a).
where EBFlow fails to accurately capture the density of the Checkerboard dataset. This phenomenon bears resemblance to the blindness problem discussed in [52] While the solution proposed in (52 has the potential to address this issue, their approach is not directly applicable to the flow-based architectures employed in this paper.
25

---

2 5
2
1 4
FDSSM
epoch:
opoch:
Figure A5: (a) Visualized examples of EBFlow trained with SSM, DSM, and FDSSM the Checkerboard dataset. (b) The samples generated by the Glow model at the 40-th training epoch. (c) The samples generated by the Glow model the 80-th training epoch.
In addition, we observed that the sampling quality of EBFlow occasionally experiences significant reduction during the training iterations. This phenomenon illustrated in Fig. (b) and (c), where the Glow model trained using approach demonstrates decline in performance with extended training periods. The underlying cause of this phenomenon remains unclear, and we consider it potential avenue for future investigation.
26


Reviews:
[{"review_id": "vwPNb08GQg", "reviewer": "Reviewer_htQa", "summary": "This paper has proposed an energy-based normalizing flow model, where the computaiton of Jacobin determinants for linear transformations can be skiped with score-matching objectives. This could enable deeper layers of linear transformation and make the flow model more efficient.", "questions": "1. Could the author explain in detail why the complexity of calcuating Jacobin determinants is $O(D^3 \\cdot L)$ for non-linear transformation? Is it general in all kinds of non-linear functions?\n2. What does the sentence in line 62-63 mean?\n3. What is the FID score in generation tasks?\n4. What is the distribution of $p_{u}$?", "limitations": "This paper has propose an interesting view of taking flow model as an eneryg-based model with score matching objectives. However, the expriments lack comparison of many baseline models as listed in 'weakness' part and other eval metrics in generation. Additional, the author should rephrase some sentences and formulars to avoid misleading.", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 3, "strengths": "This paper has proposed an interesting idea about the relationship between flow model and energy-based model. The flow model could be interpreted as an energy-based model with a tractable normalizaiton term and the proposed score-matching objective could significantly improve the efficiency.", "weaknesses": "1. The experiments lack comparison with other energy-based flow model like [1][2][3] and recent flow models like [4][5][6].\n2. Is there any theoretic analysis for the complexity between different objective function?\n3. The author misses citations for the very first deep-learning-based energy-based models like [7][8].\n\nWriting weakness:\n1. The background in section 2.2 could break up into two sections, i.e., energy-based model and score-matching model.\n2. It is hard to see the difference between results in the first two rows in Figure 1.\n3. In line 51, even the tranformation function in flow model is invertible, the author should specify that $g()$ is the inverse transformation function.\n4. In line 51, using $g_i(\\cdot;\\theta)$ is misleading. Does different transformation functions share the same parameters?\n5. I think the author intends to model $g(\\theta)$. Therefore, it is better to avoid using $E(x;\\theta)$ to avoid confusion as it is usually used for representing a parametric energy model.\n6. For Table 2, the author could change the unit to batch/second for easier comparison.\n\nReference:\n[1] Xie, Jianwen, et al. \"A Tale of Two Latent Flows: Learning Latent Space Normalizing Flow with Short-run Langevin Flow for Approximate Inference.\" arXiv preprint arXiv:2301.09300 (2023).\n\n[2] Xie, Jianwen, et al. \"A tale of two flows: Cooperative learning of langevin flow and normalizing flow toward energy-based model.\" arXiv preprint arXiv:2205.06924 (2022).\n\n[3] Nijkamp, Erik, et al. \"Learning energy-based model with flow-based backbone by neural transport mcmc.\" arXiv preprint arXiv:2006.06897 2 (2020).\n\n[4] Chen, Ricky TQ, et al. \"Residual flows for invertible generative modeling.\" Advances in Neural Information Processing Systems 32 (2019).\n\n[5] Maaløe, Lars, et al. \"Biva: A very deep hierarchy of latent variables for generative modeling.\" Advances in neural information processing systems 32 (2019).\n\n[6] Vahdat, Arash, Karsten Kreis, and Jan Kautz. \"Score-based generative modeling in latent space.\" Advances in Neural Information Processing Systems 34 (2021): 11287-11302.\n\n[7] Xie, Jianwen, et al. \"A theory of generative convnet.\" International Conference on Machine Learning. PMLR, 2016.\n\n[8] Xie, Jianwen, et al. \"Cooperative learning of energy-based model and latent variable model via mcmc teaching.\" Proceedings of the AAAI Conference on Artificial Intelligence. Vol. 32. No. 1. 2018.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer’s time and effort spent on the review, and would like to respond to the reviewer’s questions as follows. Please note that parts of the responses are provided in the [global comment](https://openreview.net/forum?id=AALLvnv95q&noteId=9f9rk2L3tW) due to the word limit.\n\n---\n### **Comments**\n---\n**C1.** The experiments lack comparison with other energy-based flow model like [1][2][3] and recent flow models like [4][5][6]. \\\n**C2.** The author misses citations for the very first deep-learning-based energy-based models like [7][8].\n\n**Response:** The references offered by the reviewer are either not applicable to the context of this work or beyond the scope of this work. Therefore, attempting a comparison of their performance may be impractical.\n\n[1] discusses the utilization of an additional decoder structure to define and train flow-based models in latent space. On the other hand, [2] and [3] explore the training techniques of optimizing an energy-based model according to a flow-based model. It is important to note that the primary emphasis of [1-3] is not on accelerating the training process of flow-based models, which constitutes the central theme of this paper.\n\nResidual Flow [4] presents an efficient method for training i-ResNet. However, the concept of [4] differs from this work, as it is derived based on the assumption of the Lipschitz constraint. In particular, each residual-flow block in [4] is constrained to have a Lipschitz constant less than 1, which is not a requirement in the architecture adopted in this paper. Due to the fundamentally different premises, a comparison between our method and [4] would provide limited insights.\n\nRegarding [5,6], they primarily concentrate on VAE and score-based models, respectively, which are not directly related to the research field of this work. Moreover, it is crucial to highlight that the primary focus of this work lies in normalizing flows rather than energy-based models. As a result, references [7,8] are also independent to the scope of this study.\n\n---\n**C3.** Is there any theoretic analysis for the complexity between different objective function? \n\n**Response:** The energy function defined in Eq. (9) has a forward propagation time complexity of $O(D^2L)$. This is based on the observation that forward passing $g$ requires $O(D^2L)$ time, and calculating the determinant of the non-linear transformations requires $O(D^2L)$ time, as specified in Lines 55-63. This implies that the complexity of calculating FDSSM is $O(D^2L)$ time.\n\nMoreover, differentiation operations exhibit the same computational complexity as the forward propagation process. As a result, the SSM and DSM objectives, along with their gradients with respect to $\\theta$, can be computed in $O(D^2L)$ time. \n\nOn the other hand, the ML and SML objectives require a complexity of $O(D^3L)$ due to specific computational requirements. The ML objective involves computing the Jacobian determinants of the linear layers, while the SML objective necessitates the computation of inverse matrices for the linear transformations during training. These operations contribute to the increased time complexity to $O(D^3L)$ when compared to the other objectives.\n\n---\n**C4.** The background in section 2.2 could break up into two sections, i.e., energy-based model and score-matching model. \n\n**Response:** We appreciate the reviewer's suggestion. However, splitting Section 2.2 into two distinct subsections, “Energy-based Models” and “Score-matching Models”, might not be ideal, as it undermines the main goal of Section 2.2 of this paper: explicitly distinguishing between parameterization and training objectives. More specifically, models trained using score-matching objectives can also fall under the category of energy-based models. In the current presentation, we first discuss the parameterization (i.e., Eq. (3)) and then describe its training methods (i.e., Eq. (4)-(7)). This arrangement provides clarity in terms of distinguishing the concepts of parameterization and training methods. Therefore, we maintain that the current presentation in Section 2.2 is more suitable.\n\n---\n**C5.** It is hard to see the difference between results in the first two rows in Figure 1. \n\n**Response:** As specified in Lines 247-249 of the manuscript, the objective of Fig. 1 is to display comparable qualitative results between the baseline and EBFlow. It is important to note that the central aim of Section 5.1 is to motivate the adoption of EBFlow by demonstrating the ability to match its performance with the baseline method. The experimental results presented in Fig. 1 substantiate this claim by showcasing similar qualitative outcomes between the baseline and EBFlow.\n\n---\n**C6.** In line 51, even the tranformation function in flow model is invertible, the author should specify that $g()$ is the inverse transformation function. \n\n**Response:** The definition of flow-based models in Eq. (1) suggests that $g$ is an inverse transformation function, despite not being literally stated. In Eq. (1), the function $g$ transforms the data vector $x$ to the latent variable $u$. This indicates that $g$ is the inverse of the generator $g^{-1}$ which maps $u$ to $x$.\n\n---\n**C7.** In line 51, using $g_i(\\cdot; \\theta)$ is misleading. Does different transformation functions share the same parameters?\n\n**Response:** Employing separate notations, such as denoting $\\theta_i$ for each $g_i$ in Line 51, could provide greater precision. However, in the context of this paper, each $g_i$ is parameterized using a subset of $\\theta$. This formulation offers notational simplicity, which contributes to the conciseness of this work.\n\n---\n**C8.** I think the author intends to model $g(\\theta)$. (...).\\\n**C9.** For Table 2, (...).\n\n**Response:** (See the global comment)\n\n--- \n### **Questions**\n---\n\n(See the global comment)\n"}, {"review_id": "eCZMxmEapZ", "reviewer": "Reviewer_DEYZ", "summary": "This paper proposes a new flow based approach to generative modeling called EBFlow. EBFlow uses training objectives from EBMs to reduce training cost for flow based methods unlike previous approaches that mostly achieved this via different architectures or biased estimation methods. The main insight for EBFlow is that the change of variables formula for a series of transformations for a NF is broken in to two parts: one containing only linear flows term and the other containing terms from non-linear transformations. The linear flows term is then interpreted as the normalizing constant and this interpretation aids in faster training strategies for flow based methods.", "questions": "See above.", "limitations": "Yes.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The motivation for reducing training costs for NFs is discussed well and is a very relevant and difficult problem. Several studies have been conducted to address this and the authors here propose a nice and novel methodology to address this.\n\nThe exploration and discussion of the method as well as the motivation is thorough. \n\nThe empirical analysis considered is detailed.", "weaknesses": "I am not sure how to interpret the empirical results. It is obvious from the plots in Fig 2 that the proposed method is much faster than ML based objective. However, the performance in all the tasks for the ML based objective is usually better or very close to EBFlow. I was curious why this is? \n\nAlso, for the results reported in Tables 1 and 2, were the methods trained for the same wall clock time?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer’s time and effort spent on the review, and would like to respond to the reviewer’s questions as follows.\n\n---\n### **Comments**\n---\n\n**C1.** I am not sure how to interpret the empirical results. It is obvious from the plots in Fig 2 that the proposed method is much faster than ML based objective. However, the performance in all the tasks for the ML based objective is usually better or very close to EBFlow. I was curious why this is? \n\n**Response:** The use of KL-divergence-based objectives (i.e., ML and SML objectives) and Fisher-divergence-based objectives (i.e., SSM, DSM, and FDSSM objectives adopted by EBFlow) exhibits a tradeoff between NLL performance and speed. KL-divergence-based objectives are expected to achieve superior NLL performance, since they explicitly minimize the NLL metric during the optimization process. On the other hand, employing Fisher-divergence-based objectives offer the advantage of efficient training by circumventing the computation of Jacobian determinants of linear transformations. This increased speed may come at the cost of NLL performance compared to the KL-divergence-based objectives, since NLL is not explicitly minimized in Fisher-divergence-based objectives.\n\n---\n\n**C2.** Also, for the results reported in Tables 1 and 2, were the methods trained for the same wall clock time?\n\n**Response:** The training wall clock time for the methods in Table 1 and 2 are different, since the models are trained until convergence as mentioned in the caption of Figure 2. We appreciate the question raised by the reviewer and will add the details in the forthcoming revision.\n"}, {"review_id": "oQ3iMrdlue", "reviewer": "Reviewer_Bq9G", "summary": "The authors use a connection between flows and EBMs to devise a new training method for flows via score-matching.\nTheir training procedure improves efficiency via avoiding computation of the Jacobian determinants that are independant to values of the samples (i.e. only contribute to the normalizing constant of the PDF).\nAdditionally they identify two methods for improving training of the flows via score-matching, namely (1) match-after-preprossing whereby the score-matching occurs to the pre-processed variables avoiding numerical instability from the pre-processing layers, and (2) using an exponential moving average for the weight updates.\nIn their experiments, they show their approach significantly improves training efficiency.", "questions": "This is the first time I've seen imputation done with flows - is this a novel contribution from the paper?\nHow would sampling work for large $D$?", "limitations": "There is no discussion around inverting the flow when $D$ is large (see Weaknesses/Questions). ", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "- The paper identifies that computing the unnormalized PDF of a flow is significantly cheaper than the normalized PDF for linear flow layers, O(D^2) instead of O(D^3), noting that computing the Jacobian determinant is not necessary for the unnormalized PDF. \nThis contribution is novel and valuable. \n- The experiments of the paper are strong, and demonstrate the utility of their approach.", "weaknesses": "For sample generation by inverting the flow, the inverse matrix for the linear flow layers will cost $D^3$ - this is not mentioned in the paper. \nThis inversion may be performed once, and then re-used for sample generation but for large $D$ this could be prohibitive. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely thank the reviewer for the time and effort spent on the review, and would like to respond to the reviewer’s questions as follows.\n\n---\n### **Comments**\n---\n\n**C1.**\n**(a)** For sample generation by inverting the flow, the inverse matrix for the linear flow layers will cost $D^3$ - this is not mentioned in the paper. \n**(b)** This inversion may be performed once, and then re-used for sample generation but for large $D$ this could be prohibitive.\n\n**Response:**\n**(a)** We agree with the reviewer that indicating the $O(D^3L)$ sampling cost in the paper enhances the clarity. This will be included in the forthcoming revision.\n\n**(b)** If the dimensionality $D$ is excessively large such that the inverse matrix of linear layers cannot be explicitly calculated, a viable alternative is to adopt a stochastic sampler such as Langevin MCMC generator discussed in Lines 327-335. The associated computational cost would depend on the number of steps (denoted as $N$) required for convergence and the time complexity of deriving the gradients of $E(x;\\theta)$ (i.e., $O(D^2L)$). This results in an overall cost $O(D^2NL)$, potentially more economical than the $O(D^3L)$ computation.\n\n---\n\n### **Questions**\n\n---\n\n**Q1. (a)** This is the first time I've seen imputation done with flows - is this a novel contribution from the paper? **(b)** How would sampling work for large $D$?\n\n**Response:**\n**(a)** Imputation for flow-based models is a well-established concept, as demonstrated in the previous work [1], which explored such an application. The process of image inpainting becomes feasible when an energy function can be explicitly defined. In the case of flow-based models, a straightforward choice is selecting $E(x;\\theta)=-\\log p(x;\\theta)$ for this purpose. Such an application of energy function would be more prevalently observed in score-based and energy-based generative modeling literature (e.g., [2-4]).\n\n**(b)** In response to the reviewer's inquiry, we have conducted supplementary experiments, and presented the imputation results of the FC-based models trained with DSM on the CelebA [5] and STL-10 [6] datasets. Both of the datasets have data dimensionality $D=$64x64x3. These results demonstrate EBFlow's potential in generating images with good quality. Please refer to Figure 1 of the PDF file in the [global response](https://openreview.net/forum?id=AALLvnv95q&noteId=9f9rk2L3tW) above.\n\n---\n\n**References:**\\\n[1] Dinh *et al.* NICE: Non-linear Independent Components Estimation, ICLR Workshop, 2015.\\\n[2] Nguyen *et al.* Plug & Play Generative Networks: Conditional Iterative Generation of Images in Latent Space, CVPR, 2017.\\\n[3] Du *et al.* Implicit Generation and Modeling with Energy Based Models, NeurIPS, 2019.\\\n[4] Song *et al.* Score-Based Generative Modeling through Stochastic Differential Equations, ICLR, 2021.\\\n[5] Liu *et al.* Deep Learning Face Attributes in the Wild, ICCV, 2015.\\\n[6] Coates *et al.* An Analysis of Single Layer Networks in Unsupervised Feature Learning, AISTAT, 2011."}, {"review_id": "7Jeujekw8s", "reviewer": "Reviewer_tDWv", "summary": "Normalizing flows are generative models trained by maximizing the log-likelihood of the dataset given the model, which is expressed as the sum of the base probability of a transformed variable plus a series of log-Jacobian terms that keep track of the volume distortion factors. Evaluating the log-Jacobian is usually the main bottleneck of normalizing flow architectures, so much so that the most popular architectures are designed so as to keep the log-Jacobian computations tractable. \n\nThis paper introduces the use of techniques from the energy-based models literature in order to more efficiently train normalizing flow architectures. In a normalizing flow architecture, it is common to have both elementwise non-linear layers, which give data dependent but tractable log-Jacobian terms, and linear layers whose log-Jacobian is data-independent but expensive to compute. The main idea introduced in this paper is to split the log-likelihood of a normalizing flow into the product of two terms: 1) an \"energy function\" that contains all the data-dependent terms and 2) a normalization constant that collects all the data-independent log-Jacobians of the linear layers.\n\nThis mirrors the split in energy-based models, where the energy function is tractable but the normalization function is intractable and it is not directly evaluated. Consequently, the authors apply a series of training objectives developed in order to train energy-based models to this new setting, obtaining an efficient way of training normalizing flows with expensive linear layers. \n\n", "questions": "- Could you clarify in which situations this model should be preferred over continuous flows and score-based diffusion models? ", "limitations": "The limitations of the approach have not be extensively discussed. I do not see any direct negative social impact.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The main idea is clever and original. It follows from the realization that the intractable terms of the log-Jacobian are independent from the data. This offers an interesting and well-motivated reason for borrowing techniques from the energy-based literature.\n\nThe experiment section is adequate, with experiments both on toy data and on real high-dimensional datasets. While the results are not spectacular, they certainly support the claim that the energy-based training can be used to train normalizing flows with a substantial speed-up and without a major loss of performance.\n\nThe paper is well-written and the relevant literature is properly referenced. \n", "weaknesses": "Major points:\nIn general, I am not convinced that this approach solves a major problem in the current literature. In fact, most successful flow architectures are designed to have fast log-Jacobian evaluations. Moreover, both continuous flows and the more recent flow-matching models can be used to train architectures with intractable Jacobian terms. \n\nThe experiments only compare with the normalizing flow baseline and therefore they do not provide evidence that this approach has a competitive advantage over continuous flows and diffusion-like models such as flow-matching models. \n\nMinor points:\nI find the presentation somehow difficult to read, with several important equations and concepts scattered all over the paper. In particular, I think it would be useful to the reader to explicitly discuss the score-matching loss associated with Eq.9, for example, by writing the concrete form of Eq.5 as applied to Eq.9 in the main text. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for thoroughly summarizing our paper, and would like to respond to the reviewer’s questions as follows.\n\n---\n### **Comments**\n---\n**C1.**\nMajor points: **(a)** In general, I am not convinced that this approach solves a major problem in the current literature. In fact, most successful flow architectures are designed to have fast log-Jacobian evaluations. **(b)** Moreover, both continuous flows and the more recent flow-matching models can be used to train architectures with intractable Jacobian terms. The experiments only compare with the normalizing flow baseline and therefore they do not provide evidence that this approach has a competitive advantage over continuous flows and diffusion-like models such as flow-matching models. \n\n**Response: (a)** We respectfully disagree with the reviewer’s comment. Architectures that prioritize fast log-Jacobian evaluations often come at the expense of limiting learnable transformations, potentially restricting their capacity to capture complex feature representations, as highlighted in [1,2]. In contrast, our work introduces an alternative approach to address training inefficiencies. This approach does not impose restrictions on linear transformations, and offers more flexibility in architectural design.\n\n**(b)** Continuous normalizing flows (CNF) (e.g., FFJORD [3] and flow-matching models [4]) reside in a research domain distinct from the central focus (i.e., (non-continuous) normalizing flows) of our study.\n\nCNF differs from (non-continuous) normalizing flows due to its reliance on an ODE solver for simulating the sampling process. This characteristic renders CNF less suited for certain downstream applications that could be efficiently achieved through (non-continuous) normalizing flows. For instance, models leveraging the latent information of a flow-based model in their training objective (e.g., [5]) may require performing backward propagations through the inverse transformation $g$. In these scenarios, the adjoint method is required for calculating the gradients by backward propagating through the ODE solutions, which can be computationally intensive in terms of both time and memory. In addition, density estimation in CNF also necessitates the use of an ODE solver, which introduces additional errors for the approximated solution (e.g., the relative and absolute tolerances in the RK45 solver). This characteristic may render CNF less common in statistical analyses demanding accurate density estimation, such as [6] and [7]. While we acknowledge the ongoing efforts of enhancing CNF, we posit that CNF and (non-continuous) normalizing flows represent distinct research avenues, with each having its unique applications.\n\n---\n**C2.** Minor points: (see above)\n\n**Response:** We appreciate the reviewer's suggestion. However, substituting $E$ in Eq. (5) with the specific form of Eq. (9) would lead to overly lengthy equations. Separating Eq. (9) and Eq. (5) offers the advantage of preserving clarity throughout the presentation.\n\n---\n### **Questions**\n---\n**Q1.** Could you clarify in which situations this model should be preferred over continuous flows and score-based diffusion models?\n\n**Response:** As explained in the response to C1 (b), (non-continuous) normalizing flows are capable of (1) exact density evaluation and (2) offering computational benefits in certain downstream applications compared to CNF, which lacks these capabilities. In a similar vein, due to these unique characteristics, score-based diffusion models also fall short in these two areas, making them unsuitable replacements for (non-continuous) normalizing flows.\n\nMore specifically, to explain why (non-continuous) normalizing flows are preferred, we first inspect the score-based diffusion models and categorize them into three different types based on their continuity with respect to $t$ and the presence of added noises. This categorization serves to clarify the difference between (non-continuous) normalizing flows and them. \n\nFirst, when score-based diffusion models are interpreted as probability flow ODEs [8,9], they fall within the category of CNFs, thus sharing the characteristics and limitations discussed in the response to C1 (b). \n\nSecond, when diffusion models are interpreted as SDEs [10], they stand apart from normalizing flows due to their intrinsic stochasticity. The stochasticity introduces difficulty in exact density estimation, leading score-based diffusion models to rely on upper-bound approximations for their likelihood [9]. \n\nThird, discrete variants of diffusion models may also require approximations for density estimation. For example, DDPM [11] leverages upper-bound to approximate negative log likelihood, setting it apart from normalizing flows in this aspect.\n\nAs a result, in certain scenarios, (non-continuous) normalizing flows are the favored choice, and diffusion models cannot adequately replace them.\n\n---\n**References:**\\\n[1] Chen et al. Residual Flows for Invertible Generative Modeling, NeurIPS, 2019.\\\n[2] Gresele et al. Relative gradient optimization of the Jacobian term in unsupervised deep learning, NeurIPS, 2020.\\\n[3] Grathwohl et al. FFJORD: Free-form Continuous Dynamics for Scalable Reversible Generative Models, ICLR, 2019.\\\n[4] Lipman et al. Flow Matching for Generative Modeling, ICLR, 2023.\\\n[5] Nalisnick et al. Hybrid Models with Deep and Invertible Features, ICML, 2019.\\\n[6] Nalisnick et al. Do Deep Generative Models Know What They Don't Know?, ICLR, 2019.\\\n[7] Jiang et al. Revisiting Flow Generative Models for Out-of-distribution Detection, ICLR, 2022.\\\n[8] Karras et al. Elucidating the Design Space of Diffusion-Based Generative Models, NeurIPS, 2022.\\\n[9] Song et al. Maximum Likelihood Training of Score-Based Diffusion Models, NeurIPS, 2021.\\\n[10] Song et al. Score-Based Generative Modeling through Stochastic Differential Equations, ICLR, 2021.\\\n[11] Ho et al. Denoising Diffusion Probabilistic Models, NeurIPS, 2020."}, {"review_id": "2RRgAdEDuP", "reviewer": "Reviewer_fofM", "summary": "The authors consider normalizing flows with linear transformations. Computing the determinant of their Jacobian is expensive, limiting their applicability when training via maximum likelihood. Therefore, the authors introduced an energy-based method to train normalizing flows with linear transformation, which circumvents the need of computing the determinant of the Jacobian. It only has to be computed once for inference, which leads to a significant speed-up in training.\n\nCompared to the maximum likelihood-based training, the authors achieve similar results while getting a significant speed-up.", "questions": "* How does your method compare against Glow in terms of runtime?\n* How does using an estimator for the Jacobian's determinant compare against your method?", "limitations": "Limitations were mentioned in the weaknesses section of the review.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "The idea of combining energy-based training with normalizing flows by phrasing the flow as an energy-based model is novel and interesting. It circumvents the need of having to compute the determinant of the Jacobian of some layers, which can save a lot of compute. It still needs to be computed for inference; however, since it is constant in the cases the authors considered, it has to be computed only once and can be stored for further use.", "weaknesses": "My biggest concern is the relevance of this work. Normalizing flows typically use transformations that are designed in such a way that the Jacobian determinant is easy to compute. For that reason, linear transformations are not very popular, also because they are not very expressive. The claim of the authors, that they are used in Glow and subsequent flow models for images, is true, but the linear transformation was only applied across channels, and since the number of channels is much smaller than the total number of input dimensions, i.e. the number of all pixels across all channels, they are not very costly in terms of evaluating their Jacobian's determinant.\n\nThe authors apply their method to MNIST and CIFAR10, roughly matching the performance of their baseline with their method. However, both perform significantly worse than Glow and related methods. Hence, it is not clear whether their method gives them any advantage over these methods, either in terms of performance or speed.\n\nThere are also procedures to estimate the Jacobian's determinant to cut down cost, such as the Skilling-Hutchinson estimator used by residual flows. The authors should have compared their method with using such an estimator. Unfortunately, their method cannot be used for residual flows, as the determinant of the Jacobian cannot depend on the input to be treated as a normalization constant in an energy-based method.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer’s time and effort spent on the review, and would like to respond to the reviewer’s questions as follows. Please note that parts of the responses are provided in the [global comment](https://openreview.net/forum?id=AALLvnv95q&noteId=9f9rk2L3tW) due to the word limit.\n\n---\n### **Comments**\n---\n\n**C1. (a)** My biggest concern is the relevance of this work. Normalizing flows typically use transformations that are designed in such a way that the Jacobian determinant is easy to compute. For that reason, linear transformations are not very popular, also because they are not very expressive. **(b)** The claim of the authors, that they are used in Glow and subsequent flow models for images, is true, but the linear transformation was only applied across channels, and since the number of channels is much smaller than the total number of input dimensions, i.e. the number of all pixels across all channels, they are not very costly in terms of evaluating their Jacobian's determinant. \n\n**Response: (a)** We value the reviewer's perspective, but respectfully hold a different viewpoint. Architectures that prioritize efficient log-Jacobian evaluations often come at the expense of limiting learnable transformations, which might compromise their capacity to capture complex feature representations, as highlighted in previous studies [r1-r3]. In contrast, our work introduces an alternative approach to address training inefficiencies. This approach does not impose restrictions on the linear transformations, and offers more flexibility in architectural design. \n\nIn addition, linear transformations are well-recognized and extensively used in the current literature. This trend of utilizing linear transformations, specifically 1x1 convolution layers, was pioneered by Glow [r4], and subsequent studies [r5-r8] have gradually relaxed the constraints on the kernel sizes. The aforementioned literature has a common goal of incorporating linear layers with a more general form in model designs, and their experimental results have supported the notion that linear transformations play a vital role in enhancing the expressiveness of flow-based models. \n\n**(b)** In response to the reviewer's comment on the computational cost, it is important to highlight that the described scenario is a niche, pertaining exclusively to cases with few channels and limited receptive fields within the model. While earlier methods [r4-r8] have gradually loosened this constraint, their methods necessitate careful management of the number of channels and the receptive fields to ensure computational feasibility. In contrast, our paper enables the use of arbitrary linear transformation in flow-based architectures. This key distinction allows the extension of EBFlow to more complex architectures.\n\n---\n\n**C2.** The authors apply their method to MNIST and CIFAR10, roughly matching the performance of their baseline with their method. However, both perform significantly worse than Glow and related methods. Hence, it is not clear whether their method gives them any advantage over these methods, either in terms of performance or speed. \n\n**Response:** We would like to emphasize that our primary objective is not to achieve state-of-the-art performance. Instead, our goal is to deliver performance comparable to contemporary works while significantly enhancing training efficiency. A unique feature of our approach is its ability to incorporate arbitrary linear transformations in flow-based models without compromising this efficiency, setting our work apart from previous research. The results demonstrated in Section 5 validate that our proposed method can train flow-based models with linear layers efficiently. This leads to considerable speed improvements while maintaining competitive performance levels. We believe that these experimental results are sufficient to substantiate the main claim of our paper.\n\n---\n\n**C3. (a)** There are also procedures to estimate the Jacobian's determinant to cut down cost, such as the Skilling-Hutchinson estimator used by residual flows. The authors should have compared their method with using such an estimator. **(b)** Unfortunately, their method cannot be used for residual flows, as the determinant of the Jacobian cannot depend on the input to be treated as a normalization constant in an energy-based method.\n\n**Response: (a)** The Skilling-Hutchinson estimator used by residual flows is not generally applicable to our model. In residual flows, each transformation block adheres to the constraint that its Lipschitz constant is less than 1 (i.e., $Lip<1$). Under such a constraint, the Jacobian determinant of each block can be estimated with convergence guarantee through an infinity power series. In contrast, the flow-based architectures adopted in this paper do not enforce constraints on Lipschitzness. This suggests that the estimator presented in [r2] is not applicable to our models, since our models can exhibit $Lip\\geq 1$.\n\n**(b)** We value the insights provided by the reviewer. However, we would like to highlight that our method is compatible with residual flows, while its efficiency enhancement depends on the architecture employed. The training efficiency of residual flows can potentially be improved under the case where a number of residual-flow blocks consist solely of linear transformations. In such a case, the Jacobian determinant of these residual-flow blocks can be treated as $Z(\\theta)$ in EBFlow, allowing its computation to be bypassed. For instance, consider a residual-flow block $g_i(x;\\theta)=x+\\tilde{W}x$, where $\\tilde{W}$ is the weight matrix with its spectral norm less than 1 (i.e., $||\\tilde{W}||_2<1$). In this case, $J{g_i}=I+\\tilde{W}$ is formulated as $Z(\\theta)$, enabling its determinant calculation to be circumvented.\n\n---\n### **Questions**\n---\n\n(see the global comment)\n\n---\n**References:**\\\n(see the global comment)\n"}]

Author rebuttal:
This global comment includes extended discussions of the questions raised by reviewers htQa and fofM. The attached PDF file contains additional experimental results.

---
### **Reviewer htQa (Cont'd)**
---
**C8.** I think the author intends to model $g(\theta)$. Therefore, it is better to avoid using $E(x;\theta)$ to avoid confusion as it is usually used for representing a parametric energy model.

**Response:** Using the notation $E$ is necessary since this paper aims to explore the connections between flow-based and energy-based models. Introducing the symbol $E$ emphasizes the perspective of viewing a flow-based model as a parametric energy-based model.

---
**C9.** For Table 2, the author could change the unit to batch/second for easier comparison.

**Response:** We thank the reviewer for the suggestion. The following tables compare our results presented in batch/second and second/batch on the MNIST dataset.

- MNIST (FC-based)
||ML|SML|SSM|DSM|FDSSM|
|-|-|-|-|-|-|
|sec. / batch|1.25 e-1|8.15 e-2|3.02 e-2|1.50 e-2|7.68 e-3|
|batch / sec.|8.00|12.27|33.11|66.67|130.21|

---
**Q1.** Could the author explain in detail why the complexity of calcuating Jacobin determinants is $O(D^3L)$ for non-linear transformation? Is it general in all kinds of non-linear functions?

**Response:** It seems that there might be some misunderstanding regarding the premise of the flow-based architecture discussed in this paper. As explained in Lines 59-61, the non-linear transformation is defined to have a computational cost of $O(D^2L)$. This complexity is determined by the sparsity of the Jacobian of $g_i \in S_n$. For instance, in Neural Spline Flows (NSF), the Jacobian determinant calculation of the non-linear transformation can be simplified as multiplying the diagonal elements of the Jacobian, resulting in a complexity that satisfies the premise stated in Lines 59-61.

---
**Q2.** What does the sentence in line 62-63 mean? 

**Response:** In Lines 62-63, we present a number of examples of implementing the flow-based architecture satisfying the criteria stated in Lines 59-61.

---
**Q3.** What is the FID score in generation tasks? 

**Response:** The FID score is seldom employed in the literature on normalizing flows. Contemporary studies, including [r4~r8], prefer NLL (Negative Log-Likelihood) and Bits/Dim (Bits per Dimension) as the evaluation metrics.

Moreover, in the case of our related works [r3,r9], which are trained on the MNIST dataset, evaluating FID would necessitate additional adjustments, such as modifying the pre-trained Inception backbone to accommodate 28x28x1 inputs. Such modifications can significantly impact the numerical results and are rarely adopted in recent literature.

---
**Q4.** What is the distribution of $p_u$?

**Response:** $p_u$ can be any continuous distribution as long as its sampling process and density estimation can be easily implemented. In this paper, we follow many popular flow-based modeling methods (e.g., [r4,r10]) to choose $p_u$ as an isotropic Gaussian with zero mean and unit variance throughout all experiments.

---
### **Reviewer fofM (Cont'd)**
---
**Q1.** How does your method compare against Glow in terms of runtime? 

**Response:** The runtimes of different losses calculated based on the Glow architecture are demonstrated in the following tables. The Glow architecture is fixed to 4 blocks to ensure the computation is manageable, and the kernel size ($ks$) of the convolutional layers within Glow varies from 1x1 to 7x7. All of the experiments are performed on an NVIDIA TITAN V GPU, and the results are reported as the average runtime for 10 independent runs. 

As shown in the rightmost three columns of both tables, SSM, DSM, and FDSSM exhibit great scalability with respect to the kernel size, as increasing the value of $ks$ does not lead to rapid growth in the average runtime. In contrast, the computational cost of the ML objective (used in the original Glow paper) increases significantly with respect to the kernel size. This phenomenon arises from the fact that transformations with larger $ks$ tend to exhibit Jacobian matrices with several non-zero elements, and thus increase the computational burden associated with the ML objective.

- $D=$28x28x1
|$ks$|Baseline (ML)|SML|SSM|DSM|FDSSM|
|-|-|-|-|-|-|
|1x1|0.269|0.201|0.340|0.189|0.114|
|3x3|0.380|0.366|0.344|0.213|0.122|
|5x5|0.554|0.460|0.346|0.213|0.122|
|7x7|0.859|0.526|0.351|0.222|0.123|

- $D=$32x32x3
|$ks$|Baseline (ML)|SML|SSM|DSM|FDSSM|
|-|-|-|-|-|-|
|1x1|2.237|1.573|0.354|0.219|0.121|
|3x3|3.547|2.104|0.395|0.227|0.124|
|5x5|6.511|3.129|0.441|0.232|0.124|
|7x7|9.249|4.156|0.443|0.238|0.128|

Please note that the results correspond to the training time measured in sec. / batch.

---
**Q2.** How does using an estimator for the Jacobian's determinant compare against your method?

**Response:** We presume the reviewer is referring to the estimator discussed in C3 (a). However, this estimator is inapplicable without the enforcement of an additional Lipschitz constraint, as explained in the response to C3 (a).

---
**References:**\
[r1] Behrmann et al. Invertible Residual Networks, ICML, 2018.\
[r2] Chen et al. Residual Flows for Invertible Generative Modeling, NeurIPS, 2019.\
[r3] Gresele et al. Relative gradient optimization of the Jacobian term in unsupervised deep learning, NeurIPS, 2020.\
[r4] Kingma et al. Glow: Generative Flow with Invertible 1×1 Convolutions, NeurIPS, 2018.\
[r5] Hoogeboom et al. Emerging Convolutions for Generative Normalizing Flows, ICML, 2019.\
[r6] Ma et al. MaCow: Masked Convolutional Generative Flow, NeurIPS, 2019.\
[r7] Lu et al. Woodbury Transformations for Deep Generative Flows, NeurIPS, 2020.\
[r8] Meng et al. ButterflyFlow: Building Invertible Layers with Butterfly Matrices, ICML, 2022.\
[r9] Pang et al. Efficient Learning of Generative Models via Finite-Difference Score Matching, NeurIPS, 2020.\
[r10] Dinh et al. Density estimation using Real NVP, ICLR, 2017.

Meta-review: The current paper uses a connection between flow-based model and energy-based model to devise a new training framework for flow-based models via score-matching. Experiments show the proposed approach significantly improves training efficiency. After the rebuttal and discussion between the authors and reviewers, the paper received all positive ratings. All five reviewers reached a consensus to accept the paper due to its substantial contributions to the realm of flow-based generative modeling, as well as convincing experimental results. The AC agrees with the reviewers and recommends accepting the paper. To further improve the paper quality, AC suggests the authors to revise their paper by taking into account all the suggestions provided by the reviewers, including completing necessary prior works about EBMs and flow models, and integrating supplementary experimental results conducted during the rebuttal phase.
Decision: Accept (poster)