Invertible Generative Modeling using Linear Rational Splines

INTRODUCTION

Flow-based modeling, widely known as normalizing flows (Tabak and Turner,, 2013; Rezende and Mohamed,, 2015), is a novel approach used for density estimation problems. The main idea behind this method is to model the distribution of any arbitrary set of data as the mapping of a simple base random variable using a set of invertible transformations. In doing so, they make use of the well-known change of variables formula from probability theory (Papoulis and Pillai,, 2002). In particular, let $\mathbf{Z}$ denote a random vector with a simple distribution such as a standard normal. Furthermore, let $\mathbf{X}$ be the result of applying an invertible transformation $\mathbf{f}(\cdot)$ on $\mathbf{Z}$. Then, it can be shown that:

Eq. (1) can be used for maximum likelihood problems. This can be done by parameterizing the distribution of the data through an invertible mapping $\mathbf{f}(\cdot)$ and then finding the parameters of this transformation using an appropriate optimization algorithm.

However, a major bottleneck to using normalizing flows in high-dimensional scenarios is the complexity of computing the Jacobian determinant. In general, determinant calculation has an \mathcal{O}\big{(}D^{3}\big{)} complexity for an arbitrary $D\times D$ matrix. Hence, to adapt flow-based models to high-dimensional cases, this issue must be addressed. As a result, the main objective in designing normalizing flow algorithms is to come up with an invertible transformation whose Jacobian determinant is tractable.

This work presents a new transformation for use in normalizing flow algorithms. We show that the transformation used has an analytical inverse. Moreover, experiments done using this transformation indicate its competitive performance with existing state-of-the-art algorithms.

RELATED WORK

In this section, we review the existing methods for flow-based modeling.

The main idea behind this category of transformations is as follows. They first split the data into two parts. The first part is output without any change. Then, a transformation is constructed based on the first partition of the data. This transformation is then applied to the second part of the data to generate the output.

where $\mathbf{g}_{\boldsymbol{\theta}({\mathbf{x}_{1}})}(\cdot)$ is an invertible, element-wise transformation with parameters $\boldsymbol{\theta}({\mathbf{x}_{1}})$ which have been computed based on $\mathbf{x}_{1}$. The final output of the transformation is then $\mathbf{x}=[\mathbf{x}_{1},\mathbf{x}_{2}]$. It can be shown that such transformations have a lower triangular Jacobian whose determinant is the multiplication of the diagonal elements. Note that coupling layers do not change the first split of their inputs. Thus, to prevent that part of the data from remaining unchanged, it is necessary to use an alternating mask and switch the splits at two consecutive transformations to make sure that every dimension of the data has a chance to be changed.

The NICE algorithm (Dinh et al.,, 2015) uses a simple translation function as its transformation $\mathbf{g}_{\boldsymbol{\theta}({\mathbf{x}_{1}})}(\cdot)$. It sets $\mathbf{g}_{\boldsymbol{\theta}({\mathbf{x}_{1}})}(\mathbf{z}_{2})=\mathbf{z}_{2}+\boldsymbol{\theta}({\mathbf{x}_{1}})$ where $\boldsymbol{\theta}(\cdot)$ is a multi-layer perceptron (MLP) with Rectified Linear Units (ReLU) as the activation function.

Real NVP (Dinh et al.,, 2017) generalizes NICE by using an affine transformation as $\mathbf{g}_{\boldsymbol{\theta}({\mathbf{x}_{1}})}(\cdot)$. Here, the authors use two Residual Networks (ResNet) (He et al.,, 2016) to come up with the translation and scaling operations required for an affine mapping. As in NICE, Real NVP uses a simple alternating mask to involve every dimension of the data in the transformation. In Kingma and Dhariwal, (2018), it is suggested to linearly combine dimensions of the data before feeding it to each layer of transformation using an invertible $1\times 1$ convolution operation, resulting in a new algorithm called Glow. This permutation is simply a matrix multiplication. Thus, to make the Glow algorithm fast, the authors propose using an LU-decomposition to calculate the matrix mentioned above. It is shown that Glow (Kingma and Dhariwal,, 2018) improves Real NVP’s performance (Dinh et al.,, 2017) for generative image modeling.

Despite its better performance compared to Real NVP, Glow may struggle to learn synthetic probability distributions with multiple separated modes (Grathwohl et al.,, 2019). Since all the previously mentioned algorithms use an affine transformation, this may be the cause of the degradation in their general expressiveness.

2 Autoregressive Transformations

The chain rule in probability theory states that if {\mathbf{Z}=\big{(}Z_{1},Z_{2},\dots,Z_{D}\big{)}} is a $D$-dimensional random variable, then the distribution of ${\mathbf{Z}}$ can be written as (Papoulis and Pillai,, 2002)

where $\mathbf{z}_{<i}$ is shorthand notation denoting all dimensions of vector $\mathbf{z}$ with an index less than $i$.

Autoregressive flows exploit this rule to build their transformations. In particular, they transform each one-dimensional variable $z_{i}$ conditioned on its previous dimensions $\mathbf{z}_{<i}$ using an invertible transformation. It can be shown that the Jacobian determinant of such transformations is again lower triangular and can be computed efficiently.

Inverse Autoregressive Flows (IAF) (Kingma et al.,, 2016) and Masked Autoregressive Flows (MAF) (Papamakarios et al.,, 2017) are among the first designs in this category. They use affine functions to build their autoregressive transformation. Later, it was argued that the use of affine functions might limit the expressiveness of such models (Huang et al.,, 2018). Hence, Neural Autoregressive Flows (NAF) (Huang et al.,, 2018) were introduced. NAFs build their transformation using a neural network with positive weights and monotonic activations to ensure invertibility. A hyper neural network, called the conditioner network, is trained to capture the dependency of the neural network parameters to the previously seen data. Later, Block Neural Autoregressive Flows (B-NAF) (Cao et al.,, 2019) suggest a more straightforward structure omitting the conditioner network. It was shown that B-NAFs reach the same performance as NAFs using orders of magnitude fewer parameters. Moreover, universal density approximation has been proved for both NAFs and B-NAFs. This theorem states that a structure of these models always exists that could accurately represent the data.

Despite their success, there are two important observations regarding NAFs and B-NAFs. First, universal density approximation only proves the existence of such models without any convergence guarantee. Moreover, these two methods are not analytically invertible, questioning their usage for generative probabilistic modeling tasks where one might want to compute the inverse.

In general, unlike coupling layer transformations, inverting an autoregressive flow cannot be done in a single pass, making their inverse computationally expensive. Thus, it would be better if one can come up with a coupling layer transformation whose performance can compete with autoregressive flows. There are also other autoregressive models to which we refer the interested reader (Germain et al.,, 2015; Chen et al.,, 2017; Oliva et al.,, 2018).

3 Other Methods

In addition to the coupling layer and autoregressive transformations, there are other methods that do not fall into any of the previous categories. We summarize the most well-known ones here.

Continuous Normalizing Flows (CNF) are flow-based models that are constructed upon Neural Ordinary Differential Equations (ODE) (Chen et al.,, 2018). FFJORD (Grathwohl et al.,, 2019) improves CNF’s computational complexity using the Hutchinson’s trace estimator (Hutchinson,, 1990). Both of these methods involve a system of first-order ODEs that replaces the change of variables formula in Eq. (1). Then, this ODE system is solved by a proper integrator resulting in a flow-based model. Here, the sampling process requires solving a system of ODEs, which may slow down sample generation.

A closely related method to FFJORD is invertible residual networks, or i-ResNets for short (Behrmann et al.,, 2019). In this paper, the authors first state a sufficient condition for making residual networks invertible. Then, a flow-based generative model is built using this invertible transformation. Realization of such models depends on calculating the Jacobian determinant of an entire residual network. In Behrmann et al., (2019), it is shown that an infinite power series can replace this Jacobian determinant. To compute this infinite series feasibly, the authors suggest truncating it after a finite number of terms $n$, which causes the Jacobian determinant estimator to become biased (Behrmann et al.,, 2019).

Residual Flows (Chen et al.,, 2019) address the bias issue of i-ResNets by using a so-called Russian roulette estimator (Kahn,, 1955). In short, instead of a deterministic $n$ used for truncation of the infinite series as in i-ResNets, a Russian roulette estimator models $n$ as a discrete random variable on natural numbers. Then, the infinite power series is replaced with the first $n$ terms divided by appropriate weights. Unlike i-ResNets, here $n$ is a realization of an arbitrary distribution on natural numbers. Since Residual Flows use an unbiased estimator for the Jacobian determinant, it is shown that they can achieve a better performance compared to i-ResNets (Behrmann et al.,, 2019).

A major drawback of models such as i-ResNets and Residual Flows, which make use of invertible residual networks as their building blocks, is their inversion. Although i-ResNets are guaranteed to be invertible, they do not have an analytical inverse. I-ResNets, as well as Residual Flows, need a fixed-point iterative algorithm to compute the inverse of each layer. Hence, the inversion process takes much more time (around 5-20x) than inference (Behrmann et al.,, 2019).

PROPOSED APPROACH

As we saw in Eq. (2.1), the only requirements that $\mathbf{g}_{\boldsymbol{\theta}({\mathbf{x}_{1}})}(\cdot)$ must satisfy is being invertible and differentiable. Also, we saw that Real NVP and Glow use a simple affine transformation, to build their algorithm and maintain analytical invertibility. In contrast, autoregressive methods either use a simple affine transformation like IAFs and MAFs, or use non-analytically invertible neural networks with positive weights as in NAFs and B-NAFs. However, there are a plethora of differentiable functions that are invertible and lie in-between: they can be more expressive than a simple affine mapping while having an easy-to-compute inverse.

A family of functions with such properties is splines. Splines are piecewise functions where each piece is expressed as a closed-form standard function. The most popular form of splines is those defined by polynomials.

In the context of flow-based modeling, the idea of replacing the affine transformation used in methods such as NICE with piecewise polynomial functions was first introduced by Müller et al., (2019). They use a piecewise linear or quadratic function as a replacement for the affine transformation in coupling layers. It is shown in that this change improves the expressiveness of methods such as Real NVP.

Later, this work was extended to the cases of cubic and rational quadratic splines (Durkan et al., 2019a, ; Durkan et al., 2019b, ). It is shown that the previous coupling layer methods such as Real NVP and Glow could benefit from this change. Also, simulation results demonstrate their competitive performance against the most expressive methods, such as NAF and B-NAF, without sacrificing analytical invertibility.

In this work, we explore the usage of another family of piecewise functions, namely linear rational splines, in the context of normalizing flows. We aim to seek more straightforward piecewise functions whose inversion can be done efficiently. Previous methods in this area use quadratic, cubic, and rational quadratic functions whose inversion is done after solving degree 2 or 3 polynomial equations. However, piecewise linear rational splines can perform competitively with these methods without requiring a polynomial equation to be solved in the inversion.

Next, we review linear rational splines and their application in constructing monotonically increasing functions. Then, based on this algorithm, we propose our flow-based model.

2 Monotonic Data Interpolation using Linear Rational Splines

Let \big{\{}\big{(}x^{(k)},~{}y^{(k)}\big{)}\big{\}}_{k=0}^{K} be a set of monotonically increasing points called knots. Furthermore, let {\big{\{}d^{(k)}>0\big{\}}_{k=0}^{K}} be a set of positive numbers representing the derivative of each point. Consider that we wish to find linear rational functionsIn the context of splines, they are also known as linear/linear rational functions. Also, they are sometimes referred to as homographic functions. of the form $y=\tfrac{ax+b}{cx+d}$ that fit the given points and their respective derivatives in each interval (also called a bin) \big{[}x^{(k)},x^{(k+1)}\big{]}. Also, consider that we require the function to be monotone.

In each bin, after satisfying function value constraints at the start and end points, we can write down the desired function as:

in which $w^{(k)}$ and $w^{(k+1)}$ are two arbitrary weights and \phi=\big{(}x-x^{(k)}\big{)}/\big{(}x^{(k+1)}-x^{(k)}\big{)}, which belongs to the interval $$. As can be seen in Eq. (4), we only have one degree of freedom left while still needing to satisfy two derivative constraints at the extreme points of the bin. Thus, we cannot use a single linear rational function and satisfy all the constraints of each bin.

Fuhr and Kallay, (1992) suggest to solve this issue by considering an intermediate point in each interval \big{(}x^{(k)},x^{(k+1)}\big{)}. This way, we can treat each bin as it was two. Since there are no constraints on this particular intermediate point in terms of the value and derivative, we can use its associated parameters to add more degrees of freedom to the existing ones. Then, we can fit one linear rational function to each one of these two intervals and satisfy the end point values and derivative constraints.

In particular, let $x^{(m)}=(1-\lambda)x^{(k)}+\lambda x^{(k+1)}$. When $0<\lambda<1$, this equation denotes a point in the interval \big{(}x^{(k)},x^{(k+1)}\big{)}. As in Fuhr and Kallay, (1992), we aim to fit two linear rational functions like Eq. (4) to the intervals \big{[}x^{(k)},x^{(m)}\big{]} and \big{[}x^{(m)},x^{(k+1)}\big{]}. Here, the value, weight and derivative of the intermediate point are treated as parameters to satisfy the value and derivative constraints that we have at each bin.

Our desired function is required to be continuous and monotonic. Since the derivative of a linear rational function does not change its sign, when we interpolate it using positive derivatives, this constraint is always satisfied. Hence, we only need to take care of the continuity of the function and its derivative at $x^{(m)}$.

After satisfying all those constraints, the final algorithm for monotonic data interpolation using linear rational functions for each interval \big{[}x^{(k)},x^{(k+1)}\big{]} can be shown as in Algorithm 1 (Fuhr and Kallay,, 1992). We can then use this algorithm for all the bins, and end up having a piecewise monotonically increasing function, also known as a linear rational spline. Specifically, for each bin \big{[}x^{(k)},x^{(k+1)}\big{]} we have:

By having a monotonically increasing function whose derivatives exist at all points, we can use it as an alternative for the function used in Eq. (2.1) and construct a flow-based model. Unlike previously used piecewise functions, linear rational splines have the advantage of having a straightforward inverse that does not require solving a degree 2 or 3 polynomial equation. Moreover, the inverse of linear rational functions has the same format as its forward form, but with different parameters. Thus, the regular and inverse function evaluations cost the same. More interestingly, having one extra degree of freedom (namely $\lambda^{(k)}$Note that although $w^{(k)}$ is chosen freely in Algorithm 1, since $w^{(m)}$ and $w^{(k+1)}$ are a multiplication of $w^{(k)}$, this value does not provide any degree of freedom.) per bin provides the opportunity to manipulate the curve that fits through a set of possible knots and derivatives. This flexibility is shown in Figure 1. Here, a set of knot points with fixed derivatives are interpolated with different $\lambda^{(k)}$s. In this figure, the assumption is that a single function shares the same value for all $\lambda^{(k)}$s. However, this does not need to be the case, and one could manipulate $\lambda^{(k)}$ of each bin separately, resulting in greater flexibility.

Note that spline transformations are defined within a finite interval. To deal with unbounded data, there are two possible solutions. First, a transformation can be used to map the unbounded data to the desired range. For instance, Müller et al., (2019) and Durkan et al., 2019a suggest mapping the data into a $$interval to deal with this issue. A downside to this method is the numerical errors caused by this extra transformation (Durkan et al., 2019a, ). As an alternative approach, Durkan et al., 2019b propose to use linear tails outside of the interval defined by piecewise functions. This way, the transformation range covers all the real line, circumventing the need to force the data itself to be in the interval$$. We also prefer this way to deal with this problem.

3 Linear Rational Spline (LRS) Flows

Having an algorithm for monotone data interpolation, we can exploit and adapt it to coupling layers, and thus, come up with a normalizing flow. We follow the steps of neural spline flows (Durkan et al., 2019b, ).

First, $K$ and $B$, the number of bins and the boundary of the spline calculation are set. Then, in order to compute the transformation as in Eq. (2.1), a neural network such as a ResNet should be trained to determine the parameters of the transformation $\mathbf{g}_{\boldsymbol{\theta}({\mathbf{x}_{1}})}(\cdot)$. As we have seen in Algorithm 1, $4K-1$ parameters are required for each dimension of this transformation: a width, height, and $\lambda$ for each bin; and $K-1$ derivatives at all points except the start and end points. At these two points the derivative is set to be $1$ for consistency with linear tails. After determining the transformation parameters for each dimension of the data, we apply the linear rational spline algorithm to come up with a set of functions that later construct different dimensions of $\mathbf{g}_{\boldsymbol{\theta}({\mathbf{x}_{1}})}(\cdot)$. As in Glow, a $1\times 1$ convolution constructed by LU-decomposition is also used to linearly combine the data that is going to be fed into the coupling layer. Furthermore, a multi-scale architecture is used for image generation scenarios as in Real NVP and Glow.

Note that all spline transformations can also be used in an autoregressive fashion. However, coupling layers are preferred as they can perform the inversion in a single pass.

Moreover, as Durkan et al., 2019b suggest, we perform a transformation on the first partition of the data using trainable parameters. In particular, let $\boldsymbol{\phi}$ be a set of trainable parameters that does not depend on the data. Then, Eq. (2.1) can be re-written as

where $\mathbf{g}_{\boldsymbol{\psi}}(\cdot)$ is an invertible, element-wise linear rational spline with parameters $\boldsymbol{\psi}$. This way all the variables are transformed while the Jacobian determinant still remains lower triangular.

Before considering the simulation results, it is worthwhile to highlight the differences between this work and previous approaches. First, here the inverse has a straightforward relationship and does not require solving degree 2 or 3 polynomial equations. Second, as in Real NVP and Glow, the inverse of this transformation is given with a similar format to its forward form: both of them are linear rational splines. This property can be useful in theoretical analysis of this flow-based model. For instance, it is sufficient to investigate mathematical properties (such as the bi-Lipschitz property that can be used for stability guarantees as in i-ResNets (Behrmann et al.,, 2019)) for a rational linear spline. Then, this property can hold for both the forward and the inverse transformations as both of them are linear rational splines. In contrast, rational quadratic splines need slightly fewer parameters. This issue can be avoided by fixing $\lambda$ in linear rational splines.

SIMULATION RESULTS

In this section, we review our simulation results. We see that the proposed method can perform competitively despite using a lower order polynomial. The code is available online at: https://github.com/hmdolatabadi/LRS_NF.

As a first experiment, we studied a density estimation scenario on synthetic 2-d distributions. This task involves the reconstruction of a continuous probability distribution given a set of its samples.

Figure 2 compares the performance of Glow (Kingma and Dhariwal,, 2018), i-ResNets (Behrmann et al.,, 2019) and our proposed method under this scenario. Qualitatively, linear rational spline (LRS) flows can reconstruct the underlying distribution precisely, and outperform the other two models.

In fact, linear rational splines can perform density estimation tasks on more complicated distributions. Figure 3 shows the result of a density estimation task, which involves a highly sophisticated distribution. In both of these experiments, our model consists of two coupling layers constructed by linear rational splines. For detailed information on the configuration used in these simulations, refer to Appendix B.1.

2 Density Estimation of Real-world Data

For the next experiment, we apply density estimation using maximum likelihood on standard benchmark datasets. Four of these datasets (Power, Gas, HEPMASS, and MiniBooNE) are tabular data from the UCI machine learning repositoryhttp://archive.ics.uci.edu/ml. Furthermore, BSDS300 is a dataset containing patches of natural images (Martin et al.,, 2001). We used the preprocessed data of Masked Autoregressive Flows (Papamakarios et al.,, 2017), which is available onlinehttps://doi.org/10.5281/zenodo.1161203.

Table 1 shows the test set log-likelihood comparison of our proposed method and FFJORD (Grathwohl et al.,, 2019), Glow (Kingma and Dhariwal,, 2018), Masked Autoregressive Flow (MAF) (Papamakarios et al.,, 2017), Neural Autoregressive Flow (NAF) (Huang et al.,, 2018), Block Neural Autoregressive Flow (B-NAF) (Cao et al.,, 2019), and Neural Spline Flows (NSF). Note that NSF models either use quadratic (Q) piecewise polynomials as in Müller et al., (2019) or rational quadratic (RQ) splines as in Durkan et al., 2019b . Moreover, recall that piecewise polynomials can be used in either coupling or autoregressive modes. These two modes are indicated in this table and the following ones with (C) and (AR), respectively. As suggested in Durkan et al., 2019b , ResMADE architecture (Durkan and Nash,, 2019) was used for autoregressive transformations.

It can be seen in Table 1 that linear rational spline flows perform competitively with RQ-NSF despite using lower degree polynomials. Furthermore, although here we set the $\lambda^{(k)}$ of each bin to be a single parameter for itself, one could consider a unified parameter to be set for the $\lambda^{(k)}$ of all intervals. While this action can slightly degrade the performance, the results are still comparable to the other methods.

3 Generative Modeling of Image Datasets

Next, we perform invertible generative modeling on benchmark image datasets including MNIST (LeCun,, 1998), CIFAR-10 (Krizhevsky and Hinton,, 2009), ImageNet $32\times 32$ and $64\times 64$ (Deng et al.,, 2009; Chrabaszcz et al.,, 2017). We measure the performance of our proposed method in bits per dimension, and then compare it with other existing methods. These include Real NVP (Dinh et al.,, 2017), Glow (Kingma and Dhariwal,, 2018), FFJORD (Grathwohl et al.,, 2019), i-ResNets (Behrmann et al.,, 2019), residual flow (Chen et al.,, 2019) and rational quadratic spline flows (Durkan et al., 2019b, ). The results are given in Table 2.

The results in Table 2 demonstrate the competitive performance of linear rational splines with respect to other methods despite their simplicity. Note that Glow uses twice as many parameters as used in neural spline flows including linear rational functions. Furthermore, although residual flows perform better than our proposed method, they require much more computation time in the sampling process where they have to compute the inverse of each invertible ResNet using a fixed point iterative method. In contrast, linear rational splines have an analytic inverse which is also a linear rational spline. This means both the forward and inverse have the same cost. In particular, our experiments indicate that the sampling process in linear rational splines is faster than residual flows by an order of magnitude.

For a detailed information on the experiments and randomly generated sample images of the model, see Appendices B.3 and C.1, respectively. Also, more simulation results can be found in Appendix D.

4 Variational Auto-encoders

Finally, we test the performance of our proposed model in a variational auto-encoder (VAE) (Kingma and Welling,, 2014) setting. In short, in a VAE data points are modeled as realizations of a random variable whose distribution is assumed to be the result of marginalization over a lower-dimensional latent variable. To make this process tractable, the type of prior and approximate posterior random variables need to be determined carefully. Like other flow-based models, linear rational spline flows can be used as effective models for priors and approximate posteriors in a VAE setting. For a detailed explanation on normalizing flows in the context of VAEs, we refer the interested reader to (Rezende and Mohamed,, 2015).

Table 3 shows the quantitative results of VAE simulation using linear rational splines. As can be seen, our proposed model performs almost as well as other models such as Glow (Kingma and Dhariwal,, 2018) and rational quadratic splines (Durkan et al., 2019b, ). For more details on the experimental configuration and model samples see Appendices B.4 and C.2, respectively.

CONCLUSION

In this paper, we investigated the use of monotonic linear rational splines in the context of invertible generative modeling. We saw that this family of piecewise functions have the advantage that their inverse is straightforward, and does not require solving degree 2 or 3 polynomial equations. Furthermore, since the same family of functions defines both the forward and inverse, investigation of the mathematical properties of these models is more straightforward. Also, we showed that despite their simplicity, they could perform competitively with more complicated methods in a suite of experiments on synthetic and real datasets.

We would like to thank the reviewers for their valuable feedback on our work, helping us to improve the final manuscript.

This research was undertaken using the LIEF HPC-GPGPU Facility hosted at the University of Melbourne. This Facility was established with the assistance of LIEF Grant LE170100200.

References

Appendix A MONOTONIC LINEAR RATIONAL SPLINES

Using the quotient rule for derivatives, the derivative of a linear rational spline function (as $g(\phi)$ in Eq. (5)) can be computed as:

To calculate the derivative with respect to $x$, we only need to divide Eq. (7) by $\delta^{(k)}=x^{(k+1)}-x^{(k)}$. As can be seen, the derivative of the function $g(x)$ never changes sign, even outside the interval \big{[}x^{(k)},x^{(k+1)}\big{]}.

A.2 Inverse Computation

Unlike rational quadratic splines which require computing the root of a degree two polynomial, linear rational splines have a straightforward closed-form inverse. This function is again a linear rational spline, but with different parameters. The inverse of Eq. (5) can be computed as:

Again, this function gives us the value of $\phi$ in each interval. We should calculate $x=\delta^{(k)}\phi+x^{(k)}$ to translate this into the interval \big{[}x^{(k)},x^{(k+1)}\big{]}.

A.3 Inverse Derivative Computation

The derivative of the inverse can be computed using the following relationship:

This function captures the change of inverse with respect to $\phi$ in each interval. To translate this into $x$, we should multiply this derivative by $\delta^{(k)}$. As in the forward pass, we can see that the derivative of the inverse does not change its sign even outside the interval $0\leq\phi\leq 1$.

Appendix B DETAILS OF SIMULATION RESULTS

For the density estimation task on the Rings dataset in Figure 2, we generated a set of 350,000 data points. Then, we used batches of size 512 to train our model, which is a linear rational spline (LRS) flow in the coupling layer mode. The number of coupling layers is 2. For the LRS function of each layer, we used 64 bins and a tail bound of 5. For optimization, we used the Adam (Kingma and Ba,, 2015) optimizer, with a learning rate of $0.0005$ and cosine annealing (Loshchilov and Hutter,, 2017). Finally, a 2-d standard normal was used as the starting probability distribution.

Note that sometimes, it is common to use an infinite data generator, which generates a different set of data at each iteration. We performed our simulation under this condition, too. The results of our method after only 50,000 iterations are depicted in Figure 4.

For the results depicted in Figure 3, we used almost the same setup as for the Rings dataset. Here, however, we used a set of $1$M data samples and $1.5$M iterations. Also, we used a uniform random variable as the starting probability distribution.

B.2 Density Estimation of Real-world Data

Also, in Table 6 we have included the results of our model under the hyper-parameters set for rational quadratic spline flows (Durkan et al., 2019b, ).

B.3 Generative Modeling of Image Datasets

For the generative modeling tasks, we used the Adam optimizer with cosine annealing of the learning rate. The initial learning rate was set to $0.0005$. For all datasets, we used batches of size $256$, and trained the model for $200$k iterations. We followed the multi-scale architecture of Dinh et al., (2017) as used in rational quadratic splines (Durkan et al., 2019b, ) and Glow (Kingma and Dhariwal,, 2018). As in Glow, each layer consists of multiple stacked steps of basic transformations, which are built by using an actnorm, a 1x1 convolution, and a coupling layer. Here, we used rational linear spline functions to build the coupling layer transformation of each layer. Moreover, a ResNet with batch normalization was used to determine the parameters of each layer’s linear rational spline functions. The detailed configuration used for the simulation of each dataset is given in Table 7.

B.4 Variational Auto-Encoders

For variational auto-encoders, we follow the same procedure as neural spline flows (Durkan et al., 2019b, ). First, a linear warm-up multiplier is used for the KL-divergence term in the cost function. This multiplier starts at the value $0.5$, and then linearly increases to $1$ as $10$% of the training set passes. A ResNet with $2$ blocks determines the parameters of the linear rational splines used in either coupling (C) or autoregressive (AR) transformations. The dimension of the latent space is set to $32$, and $64$ context features are computed by the encoder.

As before, the Adam optimizer with cosine annealing of an initial $0.0005$ learning rate is used for optimization. We use batches of size $256$, and train the model for $150$k iterations. Model selection is made using a validation set of $10$k and $20$k samples for MNIST and EMNIST, respectively.

Appendix C IMAGE SAMPLES

C.2 Randomly Generated VAE Samples

Appendix D Further Simulation Results

To highlight the improvements in the current work, we perform a new set of image generation experiments on the MNIST (LeCun,, 1998) dataset. Other than using a different family of splines, all of the hyperparameters of the models (summarized in Table 8) are fixed to be the same. For a given depth, the experiment is performed for 8 different seeds. We then train the model and repeat the same procedure for 5 various depths. In each of the experiments we pick the best flow model using a validation set. Finally, we measure the log-likelihood on the test set in BPD.

Figure 7 shows the simulation results. The top-left figure shows the average of log-likelihood on the 8 seeds. As shown, linear rational splines consistently perform better than rational quadratic splines despite using lower degree polynomials. The top-right figure shows the standard deviation of the results across different seeds. As the figure shows, the standard deviation of linear rational splines consistently decreases as the depth increases. However, the results of rational quadratic splines show fluctuations, and for the depth of 32 their standard deviation gets worse. This might be an indication of the fact that since they are using higher degree polynomials, they require more numerical accuracy as the depth increases. In contrast, as a composition of linear rational splines is still a linear rational spline, the standard deviation of our method’s results consistently decreases. Finally, you can see the number of parameters, and its relative change in percentages for these simulations in the bottom figures. As the figures show, the increase in number of parameters is only 0.23% for this set of simulations which is negligible.