Sum-of-Squares Polynomial Flow

Introduction

Neural density estimation methods are gaining popularity for the task of multivariate density estimation in machine learning [Kingma et al., 2016; Dinh et al., 2015, 2017; Papamakarios et al., 2017; Uria et al., 2016; Huang et al., 2018]. These generative models provide a tractable way to evaluate the exact density, unlike generative adversarial nets [Goodfellow et al., 2014] or variational autoencoders [Kingma & Welling, 2014; Rezende et al., 2014]. Popular methods for neural density estimation are autoregressive models [Neal, 1992; Bengio & Bengio, 1999; Larochelle & Murray, 2011; Uria et al., 2016] and normalizing flows [Rezende & Mohamed, 2015; Tabak & Vanden-Eijnden, 2010; Tabak & Turner, 2013]. These models aim to learn an invertible, bijective and increasing transformation $\mathbf{T}$ that pushes forward a (simple) source probability density (or measure, in general) to a target density such that computing the inverse $\mathbf{T}^{-1}$ and the Jacobian $|\mathbf{T}^{\prime}|$ is easy.

In probability theory, it has been rigorously proven that increasing triangular maps [Bogachev et al., 2005] are universal, i.e. any source density can be transformed into a target density using an increasing triangular map. Indeed, the Knothe-Rosenblatt transformation [Ch.1, Villani, 2008] gives a (heuristic) construction of such a map, which is unique up to null sets [Bogachev et al., 2005]. Furthermore, by definition the inverse and the Jacobian of a triangular map can be very efficiently computed through univariate operations. However, for multivariate densities computing the exact Knothe-Rosenblatt transform itself is not possible in practice. Thus, a natural question is: Given a pair of densities, how can we efficiently estimate this unique increasing triangular map?

This work is devoted to studying these increasing, bijective, and monotonic triangular maps, in particular how to estimate them in practice. In §2, we precisely formulate the density estimation problem and propose a general maximum likelihood framework for estimating densities using triangular maps. We also explore the properties of the triangular map required to push a source density to a target density.

Subsequently, in §3, we trace back the origins of the triangular map and connect it to many recent works on generative modelling. We relate our study of increasing, bijective, triangular maps to works on iterative Gaussianization [Chen & Gopinath, 2001; Laparra et al., 2011] and normalizing flows Tabak & Vanden-Eijnden ; Tabak & Turner ; Rezende & Mohamed . We show that a triangular map can be decomposed into compositions of one dimensional transformations or equivalently univariate conditional densities, allowing us to demonstrate that all autoregressive models and normalizing flows are subsumed in our general density estimation framework. As a by-product, this framework also reveals that autoregressive models and normalizing flows are in fact equivalent. Using this unified framework, we study the commonalities and differences of the various aforementioned models. Most importantly, this framework allows us to study the universality in a much cleaner and more streamlined way. We present a unified understanding of the limitations and representation power of these approaches, summarized concisely in Table 1 below.

In §4, by understanding the pivotal properties of triangular maps and using our proposed framework, we uncover a new neural density estimation procedure called the Sum-of-Squares polynomial flows (SOS flows). We show that SOS flows are akin to higher order approximation of $\mathbf{T}$ depending on the degree of the polynomials used. Subsequently, we show that SOS flows are universal, i.e. given enough model complexity, they can approximate any target density. We further show that (a) SOS flows are a strict generalization of the inverse autoregressive flow (IAF) of Kingma et al. , (b) they are interpretable; its coefficients directly control the higher order moments of the target density and, (c) SOS flows are easy to train; unlike NAFs [Huang et al., 2018] which require non-negative weights, there are no constraints on the parameters of SOS.

In §5, we report our empirical analysis. We performed holistic synthetic experiments to gain intuitive understanding of triangular maps and SOS flows in particular. Additionally, we compare SOS flows to previous neural density estimation methods on real-world datasets where it achieved competitive performance.

We summarize our main contributions as follows:

We study and propose a rigorous framework for using triangular maps for density estimation

Using this framework, we study the similarities and differences of existing flow based and autoregressive models

We provide a unified understanding of the limitations and representational power of these methods

We propose SOS flows that are universal, interpretable, and easy to train.

We perform several synthetic and real-world experiments to demonstrate the efficacy of SOS flows.

Density estimation through triangular map

In this section we set up our main problem, introduce key definitions and notations, and formulate the general approach to estimate density functions using triangular maps.

Let $p,q$ be two probability densityAll of our results can be extended to two probability measures satisfying mild regularity conditions. For simplicity and concreteness we restrict to probability densities here. functions (w.r.t. the Lebesgue measure) over the source domain $\mathsf{Z}\subseteq\mathds{R}^{d}$ and the target domain $\mathsf{X}\subseteq\mathds{R}^{d}$, respectively. Our main goal is to find a deterministic transformation $\mathbf{T}:\mathsf{Z}\to\mathsf{X}$ such that for all (measurable) set $B\subseteq\mathsf{X}$,

In particular, when $\mathbf{T}$ is bijective and differentiable [e.g. Rudin, 1987], we have the change-of-variable formula $\mathbf{x}=\mathbf{T}(\mathbf{z})$ such that

where $|\mathbf{T}^{\prime}(\mathbf{z})|$ is the (absolute value) of the Jacobian (determinant of the derivative) of $\mathbf{T}$. In other words, by pushing the source random variable $\mathbf{z}\sim p$ through the map $\mathbf{T}$ we can obtain a new random variable $\mathbf{x}\sim q$. This “push-forward” idea has played an important role in optimal transport theory [Villani, 2008] and in recent Monte carlo simulations [Marzouk et al., 2016; Parno & Marzouk, 2018; Peherstorfer & Marzouk, 2018].

Here, our interest is to learn the target density $q$ through the map $\mathbf{T}$. Let $\mathcal{F}$ be a class of mappings and use the KL divergenceOther statistical divergences can be used as well. to measure closeness between densities. We can formulate the density estimation problem as:

When we only have access to an i.i.d. sample $\lbag\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\rbag\sim q$, we can replace the integral above with empirical averages, which amounts to maximum likelihood estimation:

Conveniently, we can choose any source density $p$ to facilitate estimation. Typical choices include the standard normal density on $\mathsf{Z}=\mathds{R}^{d}$ (with zero mean and identity covariance) and uniform density over the cube $\mathsf{Z}=^{d}$.

Computationally, being able to solve (5) efficiently relies on choosing a map $\mathbf{T}$ whose

inverse $\mathbf{T}^{-1}$ is “cheap” to compute;

Jacobian $|\mathbf{T}^{\prime}|$ is “cheap” to compute.

Fortunately, this is always possible. Following Bogachev et al. we call a (vector-valued) mapping $\mathbf{T}:\mathds{R}^{d}\to\mathds{R}^{d}$ triangular if for all $j$, its $j$-th component $T_{j}$ only depends on the first $j$ variables $x_{1},\ldots,x_{j}$. The name “triangular” is derived from the fact that the derivative of $\mathbf{T}$ is a triangular matrix functionThe converse is clearly also true if our domain is connected.. We call $\mathbf{T}$ (strictly) increasing if for all $j\in[d]$, $T_{j}$ is (strictly) increasing w.r.t. the $j$-th variable $x_{j}$ when other variables are fixed.

For any two densities $p$ and $q$ over $\mathsf{Z}=\mathsf{X}=\mathds{R}^{d}$, there exists a unique (up to null sets of $p$) increasing triangular map $\mathbf{T}:\mathsf{Z}\to\mathsf{X}$ so that $q=\mathbf{T}_{\#}p$. The sameMore generally on any open or closed subset of $\mathds{R}^{d}$ if we interpret the monotonicity of $\mathbf{T}$ appropriately [Alexandrova, 2006]. holds over $\mathsf{Z}=\mathsf{X}=^{d}$.

Conveniently, to compute the Jacobian of an increasing triangular map we need only multiply $d$ univariate partial derivatives $|\mathbf{T}^{\prime}(\mathbf{x})|=\prod_{j=1}^{d}\frac{\partial T_{j}}{\partial x_{j}}.$ Similarly, inverting an increasing triangular map requires inverting $d$ univariate functions sequentially, each of which can be efficiently done through say bisection. Bogachev et al. further proved that the change-of-variable formula (3) holds for any increasing triangular map $\mathbf{T}$ (without any additional assumption but using the right-side derivative).

Thus, triangular mappings form a very appealing function class for us to learn a target density as formulated in (4) and (5). Indeed, Moselhy & Marzouk already promoted a similar idea for Bayesian posterior inference and Spantini et al. related the sparsity of a triangular map with (conditional) independencies of the target density. Moreover, many recent generative models in machine learning are precisely special cases of this approach. Before we discuss these connections, let us give some examples to help understand Theorem 1.

Consider two probability densities $p$ and $q$ on the real line $\mathds{R}$, with distribution function $F$ and $G$, respectively. Then, we can define the increasing map $T=G^{-1}\circ F$ such that $q=\mathbf{T}_{\#}p$, where $G^{-1}:\to\mathds{R}$ is the quantile function of $q$:

Let $p$ be uniform over $$and$q\sim\mathcal{N}(\mu,\sigma^{2})$ be normal distributed. The unique increasing transformation

Similar as above but we now find a map $S$ that pushes $q$ to $p$:

where $\Phi$ is the cdf of standard normal. As shown by Medvedev , $S$ must be the inverse of the map $T$ in Example 2. We observe that the derivative of $S$ is no longer a sum of squares of polynomials, but we prove later that it is approximately so. If we truncate at $k=0$, we obtain

where the leading term is also the inverse of the leading term of $T$ in (9).

We end this section with two important remarks.

If the target density has disjoint support e.g. mixture of Gaussians (MoGs) with well-separated components, then the resulting transformation will need to admit sharp jumps for areas of near zero mass. This follows by analyzing the transformaiton $T(z)=G^{-1}\circ F$. The slope $T^{\prime}(z)$ of $T(z)$ is the ratio of the quantile pdfs of the source density and the target density. Therefore, in regions of near zero mass for target density, the transformation will have near infinite slope. In Appendix B, we demonstrate this phenomena specifically for well-separated MoGs and show that a piece-wise linear function transforms a standard Gaussian to MoGs. This also opens the possibility to use the number of jumps of an estimated transformation as the indication of the number of components in the data density.

So far we have employed the (increasing) triangular map $\mathbf{T}$ explicitly to represent our estimate of the target density. This is advantageous since it allows us to easily draw samples from the estimated density, and, if needed, it results in the estimated density formula (3) immediately. An alternative would be to parameterize the estimated density directly and explicitly, such as in mixture models, probabilistic graphic models and sigmoid belief networks. The two approaches are conceptually equivalent: Thanks to Theorem 1, we know choosing a family of triangular maps fixes a family of target densities that we can represent, and conversely, choosing a family of target densities fixes a family of triangular maps that we can implicitly learn. The advantage of the former approach is that given a sample from the target density, we can infer the “pre-image” in the source domain while this information is lost in the second approach.

Connection to existing works

The results in Section 2 suggest using (5) with $\mathcal{F}$ being a class of triangular maps for estimating a probability density $q$. In this section we put this general approach into historical perspective, and connect it to the many recent works on generative modelling. Due to space constraint, we limit our discussion to work that are directly relevant to ours.

Origins of triangular map: Rosenblatt , among his contemporary peers, used the triangular map to transform a continuous multivariate distribution into the uniform distribution over the cube. Independently, Knothe devised the triangular map to transform uniform distributions over convex bodies and to prove generalizations of the Brunn-Minkowski inequality. Talagrand , unaware of the previous two results and in the process of proving some sharp Gaussian concentration inequality, effectively discovered the triangular map that transforms the Gaussian distribution into any continuous distribution. The work of Bogachev et al. rigorously established the existence and uniqueness of the triangular map and systematically studied some of its key properties. Carlier et al. showed surprisingly that the triangular map is the limit of solutions to a class of Monge-Kantorovich mass transportation problems under quadratic costs with diminishing weights. None of these pioneering works considered using triangular maps for density estimation.

Iterative Gaussianization and Normalizing Flow: In his seminal work, Huber developed the important notion of non-Gaussianality to explain the projection pursuit algorithm of Friedman et al. . Later, Chen & Gopinath , based on a heuristic argument, discovered the triangular map approach for density estimation but deemed it impractical because of the seemingly impossible task of estimating too many conditional densities. Instead, Chen & Gopinath proposed the iterative Gaussianization technique, which essentially decomposesThis can be made precise, much in the same way as decomposing a triangular matrix into the product of two rotation matrices and a diagonal matrix, i.e. the so-called Schur decomposition. the triangular map into the composition of a sequence of alternating diagonal maps $\mathbf{D}_{t}$ and linear maps $\mathbf{L}_{t}$. The diagonal maps are estimated using the univariate transform in Example 1 where $G$ is standard normal and $F$ is a mixture of standard normals. Later, Laparra et al. simplified the linear map into random rotations. Both approaches, however, suffer cubic complexity w.r.t. dimension due to generating or evaluating the linear map. The recent work of [Tabak & Vanden-Eijnden, 2010; Tabak & Turner, 2013] coined the name normalizing flow and further exploited the straightforward but crucial observation that we can approximate the triangular map through a sequence of “simple” maps such as radial basis functions or rotations composed with diagonal maps. Similar simple maps have also been explored in Ballé et al. . Rezende & Mohamed designed a “rank-1” (or radial) normalizing flow and applied it to variational inference, largely popularizing the idea in generative modelling. These approaches are not estimating a triangular map per se, but the main ideas are nevertheless similar.

(Bona fide) Triangular Approach: Deco & Brauer (see also Redlich ), to our best knowledge, is among the first to mention the name “triangular” explicitly in tasks (nonlinear independent component analysis) related to density estimation. More recently, Dinh et al. recognized the promise of even simple triangular maps in density estimation. The (increasing) triangular map in [Dinh et al., 2015] consists of two simple (block) components: $T_{1}(\mathbf{x}_{1})=\mathbf{x}_{1}$ and $T_{2}(\mathbf{x}_{1},\mathbf{x}_{2})=\mathbf{x}_{2}+m(\mathbf{x}_{1})$, where $\mathbf{x}=(\mathbf{x}_{1},\mathbf{x}_{2})$ is a two-block partition and $m$ is a map parameterized by a neural net. The advantage of this triangular map is its computational convenience: its Jacobian is trivially 1 and its inversion only requires evaluating $m$. Dinh et al. applied different partitions of variables, iteratively composed several such simple triangular maps and combined with a diagonal linear mapThey also considered a more general coupling that may no longer be triangular.. However, these triangular maps appear to be too simple and it is not clear if through composition they can approximate any increasing triangular map. In subsequent work, Dinh et al. proposed the extension where $T_{1}(\mathbf{x}_{1})=\mathbf{x}_{1}$ but $T_{2}(\mathbf{x}_{1},\mathbf{x}_{2})=\mathbf{x}_{2}\odot\exp(s(\mathbf{x}_{1}))+m(\mathbf{x}_{1})$, where $\odot$ denotes the element-wise product. This map is again increasing triangular. Moselhy & Marzouk employed triangular maps for Bayesian posterior inference, which was further extended in [Marzouk et al., 2016] for sampling from an (unknown) target density. One of their formulations is essentially the same as our eq. 4.

Autoregressive Neural Models: A joint probability density function can be factorized into the product of marginal and conditionals:

In his seminal work, Neal proposed to model each (discrete) conditional density by a simple linear logistic function (with the conditioned variables as inputs). This was later extended by Bengio & Bengio using a two-layer nonlinear neural net. The recent work of Uria et al. proposed to decouple the hidden layers in Bengio & Bengio and to introduce heavy weight sharing to reduce overfitting and computational complexity. Already in [Bengio & Bengio, 1999], univariate mixture models were mentioned as a possibility to model each conditional density, which was further substantiated in [Uria et al., 2016]. More precisely, they model the $j$-th conditional density as:

where $C_{j}$ is the so-called conditioner network that outputs the parameters for the (univariate) mixture distribution in (14). According to Example 1 there exists a unique increasing map $S_{j}(\cdot~{};\boldsymbol{\theta}_{j})$ that maps a univariate standard normal random variable $z_{j}$ into $x_{j}$ that follows (14). In other words,

where the last equality follows from induction, using the fact that $\boldsymbol{\theta}_{j}=C_{j}(x_{j-1},\ldots,x_{1})$. Thus, as already pointed out in Remark 2, specifying a family of conditional densities as in (14) is equivalent as (implicitly) specifying a family of triangular maps. In particular, if we use a nonparametric family such as mixture of normals, then the induced triangular maps can approximate any increasing triangular map. The special case, when $k=1$ in (14), was essentially dealt with by Kingma et al. : for $k=1$ the map $S_{j}(\boldsymbol{\theta}_{j})=\mu_{j}+\sigma_{j}z_{j}$ hence the triangular map

Obviously, not every triangular map can be written in the form (17), which is affine in $z_{j}$ when $z_{<j}$ are fixed. To address this issue, Kingma et al. composed several triangular maps in the form of (17), hoping this suffices to approximate a generic triangular map. In contrast, Huang et al. proposed to replace the affine form in (17) with a univariate neural net (with $z_{j}$ as input and $\mu_{j}$ and $\sigma_{j}$ serve as weights). Lastly, based on binary masks, Germain et al. and Papamakarios et al. proposed efficient implementations of the above that compute all parameters in a single pass of the conditioner network. It should be clear now that (a) autoregressive models implement exactly a triangular map; (b) specifying the conditional densities directly is equivalent as specifying a triangular map explicitly.

Other Variants. Recurrent nets have also been used in autoregressive models (effectively triangular maps). For instance, Oord et al. used LSTMs to directly specify the conditional densities while MacKay et al. chose to explicitly specify the triangular maps. The two approaches, as alluded above, are equivalent, although one may be more efficient in certain applications than the other. Oliva et al. tried to combine both while Kingma & Dhariwal used an invertible $1\times 1$ convolution. We note that the work of Ostrovski et al. models the conditional quantile function, which is equivalent to but can sometimes be more convenient than the conditional density.

Non-Triangular Flows. Sylvester Normalizing Flows (SNF) [Berg et al., 2018] and FFJORD [Grathwohl et al., 2019] are examples of normalizing flows that employ non-triangular maps. They both propose efficient methods to compute the Jacobian for change of variables. SNF utilizes Sylvester’s determinant theorem for that purpose. FFJORD, on the other hand, defines a generative model based on continuous-time normalizing flows proposed by Chen et al. and evaluates the log-density efficiently using Hutchinson’s trace estimator.

Sum-of-Squares Polynomial Flow

In Section 2 we developed a general framework for density estimation using triangular maps, and in Section 3 we showed the many recent generative models are all trying to estimate a triangular map in one way or another. In this section we give a surprisingly simple way to parameterize triangular maps, which, when plugged into (5), leads to a new density estimation algorithm that we call sum-of-squares (SOS) polynomial flow.

Our approach is motivated by some classical result on simulating univariate non-normal distributions. Let $z$ be univariate standard normal. Fleishman proposed to simulate a non-normal distribution by fitting a degree-3 polynomial:

where the coefficients $\{a_{l}\}$ are estimated by matching the first 4 moments of $x$ with those of empirical data. This approach was quite popular in practice because it allows researchers to precisely control the moments (such as skewness and kurtosis). However, three difficulties remain: (1) with degree-3 polynomial one can only (approximately) simulate a (very) strict subset of non-normal distributions. This can be addressed by using polynomials of higher degrees and better quantile matching techniques [Headrick, 2009]. (2) The estimated coefficients $\{a_{l}\}$ may not guarantee the monotonicity of the polynomial, making inversion and density evaluation difficult, if not impossible. (3) Extension to the multivariate case was done through composing a linear map [Vale & Maurelli, 1983], which can be quite inefficient.

We show that all three difficulties can be overcome using SOS flows. First, let us recall a classic result in algebra:

A univariate real polynomial is increasing iff it can be written as:

where $c\in\mathds{R}$, $r\in\mathds{N}$, and $k$ can be chosen as small as 2.

Note that a univariate increasing polynomial is strictly increasing iff it is not a constant. Theorem 2 is obtained by integrating a nonnegative polynomial, which is necessarily a sum-of-squares, see e.g. [Marshall, 2008]. Now, by applying (19) to model each conditional density in (16) we effectively addressed the last two issues above. Pleasantly, this approach strictly generalizes the affine triangular map (17) of [Kingma et al., 2016], which amounts to truncating $r=0$ in (19). However, by using a larger $r$, we can learn certain densities more faithfully (especially for capturing higher order statistics), without significantly increasing the computational complexity. Additionally, implementing (19) in practice is simple: It can be computed exactly since it is an integral of univariate polynomials.

Lastly, we prove that as the degree $r$ increases, we can approximate any triangular map. We prove our result for the domain $\mathsf{Z}=\mathsf{X}=\mathds{R}^{d}$, but the same result holds for other domains if we slightly modify the proof in the appendix.

Let $\mathsf{C}$ be the space of real univariate continuous functions, equipped with the topology of compact convergence. Then, the set of increasing polynomials is dense in the cone of increasing continuous functions.

Since the topology of pointwise convergence is weaker than that of compact convergence (i.e. uniform convergence on every compact set), we immediately know that there exists a sequence of increasing polynomials of the form (19) that converges pointwise to any given continuous function. This universal property of increasing polynomials allows us to prove the universality of SOS flows, i.e. the capability of approximating any (continuous) triangular map.

SOS flow consists of two parts: an increasing (univariate) polynomial $\mathfrak{P}_{2r+1}(z_{j};\mathbf{a}_{j})$ of the form (19) for modelling conditional densities and a conditioner network $C_{j}(z_{1},\ldots,z_{j-1})$ for generating the coefficients $\mathbf{a}_{j}$ of the polynomial $\mathfrak{P}_{2r+1}(z_{j};\mathbf{a}_{j})$. In other words, the triangular map learned using SOS flows has the following form:

If we choose a universal conditioner (that can approximate any continuous function), such as a neural net, then combining with Theorem 2 and Theorem 3 we verify that the triangular maps in the form of (20) can approximate any increasing continuous triangular map in the pointwise manner. It then follows that the transformed densities will converge weakly to any desired target density (i.e. in distribution). This solves the first issue mentioned before Theorem 2. We remark that our universality proof for SOS flows is significantly shorter and more streamlined than the previous attempt of Huang et al. , and it can be seemingly extended to analyze other models summarized in Table 1.

As pointed out by Papamakarios et al. we can also construct conditioner networks $C_{j}$ that take inputs $x_{1},\ldots,x_{j-1}$, instead of $z_{1},\ldots,z_{j-1}$. They are equivalent in theory but one can be more convenient than the other, depending on the downstream application. Figure 1 illustrates the main components of a single-block SOS flow, where we implement the conditioner network in the same way as in [Papamakarios et al., 2017]. To get a higher degree approximation, we can either increase $r$ or stack a few single-block SOS flows, as shown in Figure 2. The former approach appears to be more general but also more difficult to train due to the larger number of parameters. Indeed, the effective number of parameters for SOS flows obtained by stacking $L$ blocks with $k$ polynomials of degree $2r+1$ is $L\cdot k\cdot(r+1)$ whereas achieving the same representation with a single block wide SOS flow would require $1/2\cdot k\cdot((2r+1)^{L}-1)$ parameters. In Section 5.1 we perform simulated experiments to compare deep vs. wide SOS flows.

SOS flow is similar to the neural autoregressive flow (NAF) of Huang et al. in the sense that both are capable of approximating any (continuous) triangular map hence learning any desired target density. However, SOS flow has the following advantages:

As mentioned before, SOS flow is a strict generalization of the inverse autoregressive flow (IAF) of Kingma et al. , which corresponds to setting $r=0$.

SOS flow is more interpretable, in the sense that its parameters (i.e. coefficients of the polynomials) directly control the first few moments of the target density.

SOS flow may be easier to train, as there is no constraint on its parameters $\mathbf{a}$. In contrast, NAF needs to make sure the parameters are nonnegativeA typical remedy is to re-parameterize through an exponential transform, which, however, may lead to overflows or underflows..

Experiments

We evaluated the performance of SOS flows on both synthetic and real-world datasets for density estimation, and compare it to several alternative autoregressive models and flow based methods.

We performed a host of experiments on simulated data to gain in-depth understanding of SOS flows.

In Figure 3 we demonstrate the ability of SOS flows to represent transformations that lead to multi-modal densities by generating data from a mixture of Gaussians for two cases - well-connected and disjoint support. The true transformation can be computed exactly following Example 1. We show three transformations learned by SOS flows for each case corresponding to a deep SOS flow, wide SOS flow and wide-deep SOS flow. As is evident, SOS flows were fairly successful in learning the transformations. We further estimated the parameters of these simulated densities using Gaussian mixtures trained using maximum likelihood under three cases - exact (same number of components as target density), under-specified (lesser number of components) and over-specified. Subsequently, we plot the resulting transformation in each case following Example 1. While, GMMs with exact components work well as expected, the transformations learned by under-specified and over-specified models are not as good. This experiment also goes on to show that using a parameterized density to model conditionals is equivalent to implicitly learning a transformation. We also performed experiments to study the effect of relative ordering of variables for the conditioner network and the representational power of deep and wide SOS flows. Finally, we tested SOS flows on a suite of 2D simulated datasets – Funnel, Banana, Square, Mixture of Gaussians and Mixture of Rings. Due to space constraints we defer the figures and explanations to Appendix A.

2 Real-World Datasets

We also performed density estimation experiments on 5 real world datasets that include four datasets from the UCI repository and BSDS300. These datasets have been previously considered for comparison of flows based methods [Huang et al., 2018].

The SOS transformation was trained using maximum likelihood method with source density as standard normal distribution. We used stochastic gradient descent to train our models with a batch size of 1000, learning rate = 0.001, number of stacked blocks = 8, number of polynomials ($k$) = 5 and, degree of polynomials ($r$) = 4 with number of epochs for training = 40. We compare our method to previous works on normalizing flows and autoregressive models which include MADE-MoG [Germain et al., 2015], MAF [Papamakarios et al., 2017], MAF-MoG [Papamakarios et al., 2017], TAN [Oliva et al., 2018] and NAFs [Huang et al., 2018]. In Table 2, we report the average log-likelihood obtained using 10 fold cross-validation on held-out test sets for SOS flows. The performance reported for other methods are those reported in [Huang et al., 2018]. The results show that SOS flows are able to achieve competitive performance as compared to other methods.

Conclusion

We presented a unified framework for estimating complex densities using monotone and bijective triangular maps. The main idea is to specify one-dimensional transformations and then iteratively extend to higher-dimensions using conditioner networks. Under this framework, we analyzed popular autoregressive and flow based methods, revealed their similarities and differences, and provided a unified and streamlined approach for understanding the representation power of these methods. Along the way we uncovered a new sum-of-squares polynomial flow that we show is universal, interpretable and easy to train. We discussed the various advantages of SOS flows for stochastic simulation and density estimation, and we performed various experiments on simulated data to explore the properties of SOS flows. Lastly, SOS flows achieved competitive results on real-world datasets. In the future we plan to carry out the analysis indicated in Table 1, and to formally establish the respective advantages between deep and wide SOS flows.

Acknowledgement

We thank the reviewers for their insightful comments, and Csaba Szepesvári for bringing [Mulansky & Neamtu, 1998] to our attention, which allowed us to reduce a lengthy proof of Theorem 3 to the current slim one. We would also like to thank Ilya Kostrikov for the code which we adapted for our SOS Flow implementation. We thank Junier Oliva for pointing out an oversight about TAN in a previous draft and Youssef Marzouk for bringing additional references to our attention. Finally, we gratefully acknowledge support from NSERC. PJ was also supported by the Cheriton Scholarship, Borealis AI Fellowship and Huawei Graduate Scholarship.

References

Appendix A Simulated Experiments

Here, we explore the effect of relative ordering for the conditioner network for SOS flows as well as mixture of Gaussians. We again generated two sets of 2D densities given by $p(x_{1},x_{2})=\mathcal{N}(x_{2}~{};0,4)\mathcal{N}(x_{1}~{};0.25x_{2}^{2},1)$ and $p(x_{1},x_{2})=\mathcal{N}(x_{2}~{};2,2)\mathcal{N}(x_{1}~{};1/3x_{2}^{3},1.5)$. However, we trained both SOS flows and GMMs with the reverse order i.e. $x_{1},x_{2}$. For SOS flows we again tested using both deep and wide flows whereas for MoGs we tested with varying number of components for each conditional. We present the plots in Figure 4. The best performance here is by a deep SOS flow. Furthermore, while a flat SOS flow is able to achieve almost the same geometrical shape as the target density, its learned density still differs from the true density. For mixture of Gaussians, a large number of components for each conditional improved the performance of the resulting model.

We also test the representational power of deep and wide SOS flows and the results are given in Figure 5. In the first row, the true transformation was simulated by stacking multiple blocks of SOS transformation. Subsequently, we generated the target density using this transformation and estimated it using a deep flow, wide flow, wide-deep flow and mixture of Gaussians. In the second row, we simulated the true transformation using a single block SOS transformation and performed the same experiment as before. In both simulations, we tried to break our model by adding random noise to the coefficients of simulated transformation. As the figure shows, however, both deep and wide variants performed equally well in terms of representation. As expected however, the training time for wider flows was significantly longer than that for deeper flows.

Finally, we tested SOS flows on a suite of 2D simulated datasets – Funnel, Banana, Square, Mixture of Gaussians and Mixture of Rings. These datasets cover a broad range of geometries and have been considered before by Wenliang et al. . For these experiments, we constructed our model by stacking 3 blocks with each block being a sum of two polynomials each of degree four. We plot the log density function learned by SOS flow and the true model in Figure 6. The model is able to capture the true log density of datasets like Funnel and Banana. The true densities of Funnel and Banana are a simple linear transformation of Gaussians. Hence, flow based models that learn a continuous and smooth transformation are expected to perform well on these datasets. However, SOS demonstrates certain artifacts at the sharp corners of the Square although it is able to capture the overall density nicely. These three datasets – Funnel, Banana, and Square – were part of the unimodal simulated datasets.

The multimodal datasets included Mixture of Gaussians (MoG) and Mixture of Rings (MoR). As discussed earlier in Remark 1, when the target distribution has regions of near zero mass, the learned transformation admits sharp jumps to capture such regions. Flow based models by virtue of being invertible and smooth are often unable to learn such sharp jumps. SOS flows performs reasonably well for mixture of Gaussians although there are certain artifacts in the model that try to connect the two components. Similarly, there are some artifacts connecting the rings for the Mixture of Rings datasets. However, this issue of separated components can be dealt with relative ease in practice using clustering.

Appendix B Transformation for Mixture of Gaussians

The slope $T^{\prime}(z)$ of $T$ at any point $z$ is given by

i.e. the slope $T^{\prime}(z)$ is the ratio of probability density quantiles (pdQs) of the source random variable and the target random variable.

We now analyze the transformation required to transform a standard normal distribution to mixture of normal distributions. Figure 7 shows three columns of plots: In the leftmost column, the top plot is the source distribution ($\mathsf{Z}$) which is standard normal. The bottom plot is the target distribution for the random variable $\mathsf{X}$ which is a Gaussian mixture model with two components. The means are $-10$ and $10$, the variance is 1 and weights are 0.5 for each component respectively. The middle plot shows the transformation $T$ required to push forward a standard normal distribution to the target. In the second column of plots, we now transform a standard normal distribution to a mixture distribution but with means as -20 and 20, i.e. the components are more separated. Finally, in the plots given in the rightmost column, we transform a standard normal distribution to a mixture of three Gaussians with means -20, -5, and 15. The variances are 1 and weights are $\frac{1}{3}$ respectively.

We make the following observations here: In all three plots for the transformation, we notice that the transformation admits jumps (close to being vertical) i.e. the slope at these points is large and close to infinity. This is expected since the regions where the target has almost zero mass but the source has finite mass would lead to a slope with such behavior. In the plots, this is the region in between the components where the mass of the target density approaches zero. Furthermore, the larger this area, the longer is the height of this jump (see plots on column one and column two). With densities that have two such areas, the transformation as expected has two jumps (plots on column three). The slope of $T$ on the extremes is a constant and is equal to the standard deviation of the component on that extreme. This is because:

Appendix C Proofs

Let us define $\mathsf{P}$ to be the space of polynomials, and $\mathsf{I}$ the space of increasing functions. We need only prove on any compact set $K$, the set of polynomials of the form (19), i.e. $\mathsf{I}\cap\mathsf{P}$ thanks to Theorem 2, is dense in $\mathsf{C}(K)\cap\mathsf{I}$. By Weierstrass’ theorem we know $\mathsf{P}$ is dense in $\mathsf{C}(K)$. Moreover, the convex subset $\mathsf{I}\cap\mathsf{C}(K)$ has nonempty interior (take say a linear function with positive slope). Applying Lemma 1 above completes the proof. ∎