1 Introduction
Every day, 2500 petabytes of data are generated. Clearly, there is a need for compression to enable efficient transmission and storage of this data. Compression algorithms aim to decrease the size of representations by exploiting patterns and structure in data. In particular, lossless compression methods preserve information perfectly–which is essential in domains such as medical imaging, astronomy, photography, text and archiving. Lossless compression and likelihood maximization are inherently connected through Shannon’s source coding theorem shannon1948mathematical , i.e., the expected message length of an optimal entropy encoder is equal to the negative loglikelihood of the statistical model. In other words, maximizing the loglikelihood (of data) is equivalent to minimizing the expected number of bits required per message.
In practice, data is usually highdimensional which introduces challenges when building statistical models for compression. In other words, designing the likelihood and optimizing it for high dimensional data is often difficult. Deep generative models permit learning these complicated statistical models from data and have demonstrated their effectiveness in image, video, and audio modeling kingma2018glow ; kumar2019videoflow ; prenger2019waveglow . Flowbased generative models dinh2014nice ; dinh2016density ; papamakarios2017masked ; kingma2018glow ; grathwohl2018ffjord ; hoogeboom2019emerging are advantageous over other generative models: i)
they admit exact loglikelihood optimization in contrast with Variational AutoEncoders (VAEs)
kingma2014auto and ii) drawing samples (and decoding) is comparable to inference in terms of computational cost, as opposed to PixelCNNs van2016pixel. However, flowbased models are generally defined for continuous probability distributions, disregarding the fact that digital media is stored discretely–for example, pixels from 8bit images have 256 distinct values. In order to utilize continuous flow models for compression, the latent space must be quantized. This produces reconstruction errors in image space, and is therefore not suited for lossless compression.
To circumvent the (de)quantization issues, we propose Integer Discrete Flows (IDFs), which are invertible transformations for ordinal discrete data–such as images, video and audio. We demonstrate the effectiveness of IDFs by attaining stateoftheart lossless compression performance on CIFAR10, ImageNet32 and ImageNet64. For a graphical illustration of the coding steps, see Figure 1. In addition, we show that IDFs achieve generative modelling results competitive with other flowbased methods. The main contributions of this paper are summarized as follows: 1) We introduce a generative flow for ordinal discrete data (Integer Discrete Flow), circumventing the problem of (de)quantization; 2) As building blocks for IDFs, we introduce a flexible transformation layer called integer discrete coupling; 3) We propose a neural network based compression method that leverages IDFs; and 4) We empirically show that our image compression method allows for progressive decoding that maintains the global structure of the encoded image. Code to reproduce the experiments will be released at a later stage.
2 Background
The continuous change of variables formula lies at the foundation of flowbased generative models. It admits exact optimization of a (data) distribution using a simple distribution and a learnable bijective map. Let be a bijective map, and a prior distribution on . The model distribution can then be expressed as:
(1) 
That is, for a given observation , the likelihood is given by evaluated at
, normalized by the Jacobian determinant. A composition of invertible functions, which can be viewed as a repeated application of the change of variables formula, is generally referred to as a normalizing flow in the deep learning literature
deco1995decorr ; tabak2010density ; tabak2013family ; rezende2015norm .2.1 Flow Layers
The design of invertible transformations is integral to the construction of normalizing flows. In this section two important layers for flowbased generative modelling are discussed.
Coupling layers are tractable bijective mappings that are extremely flexible, when combined into a flow dinh2016density ; dinh2014nice
. Specifically, they have an analytical inverse, which is similar to a forward pass in terms of computational cost and the Jacobian determinant is easily computed, which makes coupling layers attractive for flow models. Given an input tensor
, the input to a coupling layer is partitioned into two sets such that . The transformation, denoted , is then defined by:(2) 
where denotes elementwise multiplication and and may be modelled using neural networks. Given this, the inverse is easily computed, i.e., , and , where denotes elementwise division. For to be invertible, must not be zero, and is often constrained to have strictly positive values.
Factorout layers allow for more efficient inference and hierarchical modelling. A general flow, following the change of variables formula, is described as a single map . This implies that a
dimensional vector is propagated throughout the whole flow model. Alternatively, a part of the dimensions can already be
factoredout at regular intervals dinh2016density , such that the remainder of the flow network operates on lower dimensional data. We give an example for two levels () although this principle can be applied to an arbitrary number of levels:(3) 
where and . The likelihood of is then given by:
(4) 
This approach has two clear advantages. First, it admits a factored model for , , which allows for conditional dependence between parts of . This holds because the flow defines a bijective map between and . Second, the lower dimensional flows are computationally more efficient.
2.2 Entropy Encoding
Lossless compression algorithms map every input to a unique output and are designed to make probable inputs shorter and improbable inputs longer. Shannon’s source coding theorem shannon1948mathematical states that the optimal code length for a symbol is , and the minimum expected code length is lowerbounded by the entropy:
(5) 
where denotes the encoded message, is length, denotes entropy, is the data distribution, and is the statistical model that is used by the encoder. Therefore, maximizing the model loglikelihood is equivalent to minimizing the expected number of bits required per message, when the encoder is optimal.
Stream coders encode sequences of random variables with different probability distributions. They have nearoptimal performance, and they can meet the entropybased lower bound of Shannon
rissanen1979arithmetic ; moffat1998arithmetic . In our experiments, the recently discovered and increasingly popular stream coder rANS duda2013asymmetric is used. It has gained popularity due to its computational and coding efficiency. See Appendix A.1 for an introduction to rANS.3 Integer Discrete Flows
We introduce Integer Discrete Flows (IDFs): a bijective integer map that can represent rich transformations. IDFs can be used to learn the probability mass function on (highdimensional) ordinal discrete data. Consider an integervalued observation , a prior distribution with support on , and a bijective map defined by an IDF. The model distribution can then be expressed as:
(6) 
Note that in contrast to Equation 1, there is no need for renormalization using the Jacobian determinant. Deep IDFs are obtained by stacking multiple IDF layers , which are guaranteed to be bijective if the individual maps are all bijective. For an individual map to be bijective, it must be onetoone and onto. Consider the bijective map . Although, this map is a bijection, it requires us to keep track of the codomain of , which is impracticable in the case of many dimensions and multiple layers. Instead, we design layers to be bijective maps from to , which ensures that the composition of layers and its inverse is closed on .
3.1 Integer Discrete Coupling
As a building block for IDFs, we introduce integer discrete coupling layers. These are invertible and the set is closed under their transformations. Let be an input of the layer. The output is defined as a copy , and a transformation , where denotes a nearest rounding operation and is a neural network (Figure 2).
Notice the multiplication operation in standard coupling is not used in integer discrete coupling, because it does not meet our requirement that the image of the transformations is equal to . It may seem disadvantageous that our model only uses translation, also known as additive coupling, however, largescale continuous flow models in the literature tend to use additive coupling instead of affine coupling kingma2018glow .
In contrast to existing coupling layers, the input is split in 75%–25% parts for and , respectively. As a consequence, rounding is applied to fewer dimensions, which results in less gradient bias. In addition, the transformation is richer, because it is conditioned on more dimensions. Empirically this results in better performance.
Backpropagation through Rounding Operation As shown in Figure 2, a coupling layer in IDF requires a rounding operation (
) on the predicted translation. Since the rounding operation is effectively a step function, its gradient is zero almost everywhere. As a consequence, the rounding operation is inherently incompatible with gradient based learning methods. In order to backpropagate through the rounding operations, we make use of the Straight Through Estimator (STE)
bengio2013estimating . In short, the STE ignores the rounding operation during backpropagation, which is equivalent to redefining the gradient of the rounding operation as follows:(7) 
Lower Triangular Coupling
There exists a tradeoff between the number of integer discrete coupling layers and the complexity of the layers in IDF architectures, due to the gradient bias that is introduced by the rounding operation (see section 4.1). We introduce a multivariate coupling transformation called Lower Triangular Coupling, which is specifically designed such that the number of rounding operations remains unchanged.
This leads to marginally improved performance. For more details, see Appendix B.
3.2 Tractable Discrete distribution
As discussed in Section 2, a simple distribution is posed on in flowbased models. In IDFs, the prior is a factored discretized logistic distribution (DLogistic) kingma2016improved ; salimans2017pixelcnn++ . The discretized logistic captures the inductive bias that values close together are related, which is wellsuited for ordinal data.
The probability mass for an integer , mean , and scale is defined as the density assigned to the interval
by the probability density function of
(see Figure 3). This can be efficiently computed by evaluating the cumulative distribution function twice:
(8) 
where denotes the sigmoid, the cumulative distribution function of a standard Logistic. In the context of a factorout layer, the mean and scale are conditioned on the subset of
data that is not factored out. That is, the input to the th factorout layer is split into and . The conditional distribution on is then given as , where and are parametrized as neural networks.
Discrete Mixture distributions The discretized logistic distribution is unimodal and therefore limited in complexity. With a marginal increase in computational cost, we increase the flexibility of the latent prior on by extending it to a mixture of logistic distributions salimans2017pixelcnn++ :
(9) 
Note that as , the mixture distribution can model arbitrary univariate discrete distributions. In practice, we find that a limited number of mixtures () is usually sufficient for image density modelling tasks.
3.3 Lossless Source Compression
Lossless compression is an essential technique to limit the size of representations without destroying information. Methods for lossless compression require i) a statistical model of the source, and ii) a mapping from source symbols to bit streams.
IDFs are a natural statistical model for lossless compression of ordinal discrete data, such as images, video and audio. They are capable of modelling complicated highdimensional distributions, and they provide errorfree reconstructions when inverting latent representations. The mapping between symbols and bit streams may be provided by any entropy encoder. Specifically, stream coders can get arbitrarily close to the entropy regardless of the symbol distributions, because they encode entire sequences instead of a single symbol at a time.
In the case of compression using an IDF, the mapping is defined by the IDF. Subsequently, is encoded under the distribution to a bitstream using an entropy encoder. Note that, when using factorout layers, is also defined using the IDF. Finally, in order to decode a bitstream , an entropy encoder uses to obtain . and the original image is obtained by using the map , i.e., the inverse IDF. See Figure 1 for a graphical depiction of this process.
In rare cases, the compressed file may be larger than the original. Therefore, following established practice in compression algorithms, we utilize an escape bit. That is, the encoder will decide whether to encode the message or save it in raw format and encode that decision into the first bit.
4 Architecture
The IDF architecture is split up into one or more levels. Each level consists of a squeeze operation dinh2016density , integer flow layers, and a factorout layer. Hence, each level defines a mapping from to , except for the final level , which defines a mapping . Each of the integer flow layers per level consist of a permutation layer followed by an integer discrete coupling layer. Following dinh2016density , the permutation layers are initialized once and kept fixed throughout training and evaluation. Figure 5 shows a graphical illustration of a two level IDF. The specific architecture details for each experiment are presented in Appendix D.1. In the remainder of this section, we discuss the tradeoff between network depth and performance when rounding operations are used.
4.1 Flow Depth and Network Depth
The performance of IDFs depends on a tradeoff between complexity and gradient bias, influenced by the number of rounding functions. Increasing the performance of standard normalizing flows is often achieved by increasing the depth, i.e., the number of flowmodules. However, for IDFs each flowmodule results in additional rounding operations that introduce gradient bias. As a consequence, adding more flow layers hurts performance, after some point, as is depicted in Figure 5. We found that the limitation of using fewer coupling layers in an IDF can be negated by increasing the complexity of the neural networks part of the coupling and factorout layers. That is, we use DenseNets huang2017densely in order to predict the translation in the integer discrete coupling layers and and in the factorout layers.
5 Related Work
There exist several deep generative modelling frameworks. This work builds mainly upon flowbased generative models, described in rippel2013high ; dinh2014nice ; dinh2016density
. In these works, invertible functions for continuous random variables are developed. However, quantizing a latent representation, and subsequently inverting back to image space may lead to reconstruction errors
dewitte1997lossless ; calderbank1997lossless ; calderbank1998wavelet .Other likelihoodbased models such as PixelCNNs van2016pixel
utilize a decomposition of conditional probability distributions. However, this decomposition assumes an order on pixels which may not reflect the actual generative process. Furthermore, drawing samples (and decoding) is generally computationally expensive. VAEs
kingma2014auto optimize a lower bound on the log likelihood instead of the exact likelihood and can be used, in combination with bitsback coding, for lossless compression townsend2019practical ; kingma2019bit . However, the performance of this approach is bounded by the lower bound. Moreover, encoding the initial data example is inefficient and the extra bits should be random, which is not the case in practice and may lead to coding inefficiencies townsend2019practical .Nonlikelihood based generative models tend to utilize Generative Adversarial Networks goodfellow2014generative , and can generate highquality images. However, since GANs do not optimize for likelihood, which is directly connected to the expected number of bits in a message, they are not suited for lossless compression.
In the lossless compression literature, numerous reversible integer to integer transforms have been proposed ahmed1974discrete ; dewitte1997lossless ; calderbank1997lossless ; calderbank1998wavelet . Specifically, lossless JPEG2000 uses a reversible integer wavelet transform jp2_standard . However, because these transformations are largely handdesigned, they are difficult to tune for realworld data, which may require complicated nonlinear transformations.
Around time of submission, unpublished concurrent work appeared discrete2019tran that explores discrete flows. The main differences between our method and this work are: i) we propose discrete flows for ordinal discrete data (e.g. audio, video, images), whereas they are are focused on categorical data. ii) we provide a connection with the source coding theorem, and present a compression algorithm. iii) We present results on more largescale image datasets.
6 Experiments
To test the compression performance of IDFs, we compare with a number of established lossless compression methods: PNG png_standard ; JPEG2000 jp2_standard ; FLIF sneyers2016flif
, a recent format that uses machine learning to build decision trees for efficient coding; and BitSwap
kingma2019bit , a VAE based lossless compression method. We show that IDFs outperform all these formats on CIFAR10, ImageNet32 and ImageNet64. In addition, we demonstrate that IDFs can be very easily tuned for specific domains, by compressing the ER + BCa histology dataset. For the exact treatment of datasets and optimization procedures, see Section D.4.6.1 Image Compression
The compression performance of IDFs is compared with competing methods on standard datasets, in bits per dimension and compression rate. The IDFs and BitSwap are trained on the train data, and compression performance of all methods is reported on the test data in Table 1. IDFs achieve stateoftheart lossless compression performance on all datasets.
Even though one can argue that a compressor should be tuned for the source domain, the performance of IDFs is also examined on outofdataset examples, in order to evaluate compression generalization. We utilize the IDF trained on Imagenet32, and compress the CIFAR10 and ImageNet64 data. For the latter, a single image is split into four patches. Surprisingly, the IDF trained on ImageNet32 (IDF) still outperforms the competing methods showing only a slight decrease in compression performance on CIFAR10 and ImageNet64, compared to its sourcetrained counterpart.
As an alternative method for lossless compression, one could quantize the distribution and the latent space of a continuous flow. This results in reconstruction errors that need to be stored in addition to the latent representation , such that the original data can be recovered perfectly. We show that this scheme is ineffective for lossless compression. Results are presented in Appendix C.
Dataset  IDF  IDF  BitSwap  FLIF sneyers2016flif  PNG  JPEG2000 

CIFAR10  3.34 (2.40)  3.60 (2.22)  3.82 (2.09)  4.37 (1.83)  5.89 (1.36)  5.20 (1.54) 
ImageNet32  4.18 (1.91)  4.18 (1.91)  4.50 (1.78)  5.09 (1.57)  6.42 (1.25)  6.48 (1.23) 
ImageNet64  3.90 (2.05)  3.94 (2.03 )  –  4.55 (1.76)  5.74 (1.39)  5.10 (1.56) 
6.2 Tuneable Compression
Thus far, IDFs have been tested on standard machine learning datasets. In this section, IDFs are tested on a specific domain, medical images. In particular, the ER + BCa histology dataset janowczyk2018resolution is used, which contains 141 regions of interest scanned at , where each image is pixels (see Figure 6, left). Since current hardware does not support training on such large images directly, the model is trained on random px patches. See Figure 6, right for samples from the model. Likewise, the compression is performed in a patchbased manner, i.e., each patch is compressed independently of all other patches. IDFs are again compared with FLIF and JPEG2000, and also with a modified version of JPEG2000 that has been optimized for virtual microscopy specifically, named JP2WSI helin2018optimized . Although the IDF is at a disadvantage because it has to compress in patches, it considerably outperforms the established formats, as presented in Table 2.
Dataset  IDF  JP2WSI  FLIF sneyers2016flif  JPEG2000 

Histology  2.42 (3.19)  3.04 (2.63)  4.00 (2.00)  4.26 (1.88) 
6.3 Progressive Image Rendering
In general, transferring data may take time because of slow internet connections or disk I/O. For this reason, it is desired to progressively visualize data, i.e., to render the image with more detail as more data arrives. Several graphics formats support progressive loading. However, the encoded file size may increase by enabling this option, depending on the format png_standard , whereas IDFs support progressive rendering naturally. To partially render an image using IDFs, first the received variables are decoded. Next, using the hierarchical structure of the prior and ancestral sampling, the remaining dimensions are obtained. The progressive display of IDFs for ImageNet64 is presented in Figure 8, where the rows use approximately 15%, 30%, 60%, and 100% of the bitstream. The global structure is already captured by smaller fragments of the bitstream, even for fragments that contain only 15% of the stream.
6.4 Probability Mass Estimation
In addition to a statistical model for compression, IDFs can also be used for image generation and probability mass estimation. Samples are drawn from an ImageNet 3232 IDF and presented in Figure 7. IDFs are compared with recent flowbased generative models, RealNVP dinh2016density , Glow kingma2018glow , and Flow++ in analytical bits per dimension (negative loglikelihood). To compare architectural changes, we modify the IDFs to Continuous models by dequantizing, disabling rounding, and using a continuous prior. The continuous versions of IDFs tend to perform slightly better, which may be caused by the gradient bias on the rounding operation. IDFs show competitive performance on CIFAR10, ImageNet32, and ImageNet64, as presented in Table 3. Note that in contrast with IDFs, RealNVP uses scale transformations, Glow has convolutions and actnorm layers for stability, and Flow++ uses the aforementioned, and an additional flow for dequantization. Interestingly, IDFs have comparable performance even though the architecture is relatively simple.
Dataset  IDF  Continuous  RealNVP  Glow  Flow++ 

CIFAR10  3.32  3.31  3.49  3.35  3.08 
ImageNet32  4.16  4.13  4.28  4.09  3.86 
ImageNet64  3.90  3.85  3.98  3.81  3.69 
7 Conclusion
We have introduced Integer Discrete Flows, flows for ordinal discrete data that can be used for deep generative modelling and neural lossless compression. We show that IDFs are competitive with current flowbased models, and that we achieve stateoftheart lossless compression performance on CIFAR10, ImageNet32 and ImageNet64. To the best of our knowledge, this is the first lossless compression method that uses invertible neural networks.
References
 [1] Nasir Ahmed, T Natarajan, and Kamisetty R Rao. Discrete cosine transform. IEEE transactions on Computers, 100(1):90–93, 1974.
 [2] Yoshua Bengio, Nicholas Léonard, and Aaron Courville. Estimating or propagating gradients through stochastic neurons for conditional computation. arXiv preprint arXiv:1308.3432, 2013.
 [3] A Robert Calderbank, Ingrid Daubechies, Wim Sweldens, and BoonLock Yeo. Lossless image compression using integer to integer wavelet transforms. In Proceedings of International Conference on Image Processing, volume 1, pages 596–599. IEEE, 1997.
 [4] AR Calderbank, Ingrid Daubechies, Wim Sweldens, and BoonLock Yeo. Wavelet transforms that map integers to integers. Applied and computational harmonic analysis, 5(3):332–369, 1998.
 [5] Gustavo Deco and Wilfried Brauer. Higher Order Statistical Decorrelation without Information Loss. In G. Tesauro, D. S. Touretzky, and T. K. Leen, editors, Advances in Neural Information Processing Systems 7, pages 247–254. MIT Press, 1995.
 [6] Steven Dewitte and Jan Cornelis. Lossless integer wavelet transform. IEEE signal processing letters, 4(6):158–160, 1997.
 [7] Laurent Dinh, David Krueger, and Yoshua Bengio. NICE: Nonlinear independent components estimation. 3rd International Conference on Learning Representations, ICLR, Workshop Track Proceedings, 2015.
 [8] Laurent Dinh, Jascha SohlDickstein, and Samy Bengio. Density estimation using Real NVP. 5th International Conference on Learning Representations, ICLR, 2017.
 [9] Jarek Duda. Asymmetric numeral systems. arXiv preprint arXiv:0902.0271, 2009.
 [10] Jarek Duda. Asymmetric numeral systems: entropy coding combining speed of huffman coding with compression rate of arithmetic coding. arXiv preprint arXiv:1311.2540, 2013.
 [11] International Organization for Standardization. JPEG 2000 image coding system. ISO Standard No. 154441:2016, 2003.
 [12] International Organization for Standardization. Portable Network Graphics (PNG): Functional specification. ISO Standard No. 15948:2003, 2003.
 [13] Ian Goodfellow, Jean PougetAbadie, Mehdi Mirza, Bing Xu, David WardeFarley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. Generative adversarial nets. In Advances in neural information processing systems, pages 2672–2680, 2014.
 [14] Will Grathwohl, Ricky TQ Chen, Jesse Betterncourt, Ilya Sutskever, and David Duvenaud. Ffjord: Freeform continuous dynamics for scalable reversible generative models. 7th International Conference on Learning Representations, ICLR, 2019.
 [15] Henrik Helin, Teemu Tolonen, Onni Ylinen, Petteri Tolonen, Juha Näpänkangas, and Jorma Isola. Optimized jpeg 2000 compression for efficient storage of histopathological wholeslide images. Journal of pathology informatics, 9, 2018.
 [16] Emiel Hoogeboom, Rianne van den Berg, and Max Welling. Emerging convolutions for generative normalizing flows. Proceedings of the 36th International Conference on Machine Learning, 2019.

[17]
Gao Huang, Zhuang Liu, Laurens Van Der Maaten, and Kilian Q Weinberger.
Densely connected convolutional networks.
In
Proceedings of the IEEE conference on computer vision and pattern recognition
, pages 4700–4708, 2017.  [18] Andrew Janowczyk, Scott Doyle, Hannah Gilmore, and Anant Madabhushi. A resolution adaptive deep hierarchical (radhical) learning scheme applied to nuclear segmentation of digital pathology images. Computer Methods in Biomechanics and Biomedical Engineering: Imaging & Visualization, 6(3):270–276, 2018.
 [19] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. 3rd International Conference on Learning Representations, ICLR, 2015.
 [20] Diederik P Kingma, Tim Salimans, Rafal Jozefowicz, Xi Chen, Ilya Sutskever, and Max Welling. Improved variational inference with inverse autoregressive flow. In Advances in Neural Information Processing Systems, pages 4743–4751, 2016.
 [21] Diederik P Kingma and Max Welling. AutoEncoding Variational Bayes. In Proceedings of the 2nd International Conference on Learning Representations, 2014.
 [22] Durk P Kingma and Prafulla Dhariwal. Glow: Generative flow with invertible 1x1 convolutions. In Advances in Neural Information Processing Systems, pages 10236–10245, 2018.
 [23] Friso H Kingma, Pieter Abbeel, and Jonathan Ho. Bitswap: Recursive bitsback coding for lossless compression with hierarchical latent variables. 36th International Conference on Machine Learning, 2019.
 [24] Manoj Kumar, Mohammad Babaeizadeh, Dumitru Erhan, Chelsea Finn, Sergey Levine, Laurent Dinh, and Durk Kingma. Videoflow: A flowbased generative model for video. arXiv preprint arXiv:1903.01434, 2019.
 [25] Alistair Moffat, Radford M Neal, and Ian H Witten. Arithmetic coding revisited. ACM Transactions on Information Systems (TOIS), 16(3):256–294, 1998.
 [26] George Papamakarios, Iain Murray, and Theo Pavlakou. Masked autoregressive flow for density estimation. In Advances in Neural Information Processing Systems, pages 2338–2347, 2017.

[27]
Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary
DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer.
Automatic differentiation in PyTorch.
In NIPS Autodiff Workshop, 2017.  [28] Ryan Prenger, Rafael Valle, and Bryan Catanzaro. Waveglow: A flowbased generative network for speech synthesis. In ICASSP 20192019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3617–3621. IEEE, 2019.
 [29] Danilo Rezende and Shakir Mohamed. Variational Inference with Normalizing Flows. In Proceedings of the 32nd International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pages 1530–1538. PMLR, 2015.
 [30] Oren Rippel and Ryan Prescott Adams. Highdimensional probability estimation with deep density models. arXiv preprint arXiv:1302.5125, 2013.
 [31] Jorma Rissanen and Glen G Langdon. Arithmetic coding. IBM Journal of research and development, 23(2):149–162, 1979.
 [32] Tim Salimans, Andrej Karpathy, Xi Chen, and Diederik P Kingma. PixelCNN++: Improving the pixelcnn with discretized logistic mixture likelihood and other modifications. 5th International Conference on Learning Representations, ICLR, 2017.
 [33] Claude Elwood Shannon. A mathematical theory of communication. Bell system technical journal, 27(3):379–423, 1948.
 [34] Jon Sneyers and Pieter Wuille. Flif: Free lossless image format based on maniac compression. In 2016 IEEE International Conference on Image Processing (ICIP), pages 66–70. IEEE, 2016.
 [35] EG Tabak and Cristina V Turner. A family of nonparametric density estimation algorithms. Communications on Pure and Applied Mathematics, 66(2):145–164, 2013.
 [36] Esteban G Tabak, Eric VandenEijnden, et al. Density estimation by dual ascent of the loglikelihood. Communications in Mathematical Sciences, 8(1):217–233, 2010.
 [37] James Townsend, Tom Bird, and David Barber. Practical lossless compression with latent variables using bits back coding. 7th International Conference on Learning Representations, ICLR, 2019.
 [38] Dustin Tran, Keyon Vafa, Kumar Agrawal, Laurent Dinh, and Ben Poole. Discrete flows: Invertible generative models of discrete data. ICLR 2019 Workshop DeepGenStruct, 2019.

[39]
Rianne van den Berg, Leonard Hasenclever, Jakub M Tomczak, and Max Welling.
Sylvester normalizing flows for variational inference.
ThirtyFourth Conference on Uncertainty in Artificial Intelligence, UAI
, 2018. 
[40]
Aaron Van Oord, Nal Kalchbrenner, and Koray Kavukcuoglu.
Pixel recurrent neural networks.
In International Conference on Machine Learning, pages 1747–1756, 2016.
Appendix A Additional background
a.1 Asymmetric Numeral Systems
Asymmetric Numeral Systems (ANS) [9] is a recent approach to entropy coding. The rangebased variant: rANS, is generally used as a faster replacement for arithmetic coding, because a state is only represented by a single number and fewer mathematical operations are required [10].
The encoding function of rANS encodes a symbol into a code given the so far existing code :
(10) 
where is a large integer that functions as the quantization denominator. Integers are chosen for such that , where denotes the probability of symbol . Each symbol is associated with a unique interval , where , as depicted in Figure 9.
The decoding function needs to retrieve the encoded symbol , and the previous state from the new code . First consider the term , which is equal to the last two terms of the encoding function: . This term is guaranteed to lie in the interval , . Therefore, the symbol can be retrieved by finding:
(11) 
Consequently with the knowledge of , the previous state can be obtained by computing:
(12) 
In practice, is chosen as a power of two (for example ). As such, multiplication and division with reduces to bit shifts and modulo reduces to a binary masking operation.
Appendix B Lower Triangular Coupling
There exists a tradeoff between the number of integer discrete coupling layers and the complexity of the layers in IDF architectures, due to the gradient bias that is introduced by the rounding operation. For this reason, it is desired to increase the flexibility of layers without increasing the number of rounding operations. We introduce a multivariate coupling transformation called Lower Triangular Coupling, which is specifically designed such that the number of rounding operations remains unchanged. The transformation of is formed by multiplication with a strictly lower triangular matrix which is conditioned on :
(13) 
The main trick is to round the sum of all transformations, such that no additional gradient bias is introduced. This transformation is guaranteed to be invertible, and the inverse can be found with a modified version of forward substitution:
(14) 
where denotes the th element of , and and are still conditioned on , however, this notation is dropped for clarity. The continuous case can even be solved analytically by using the inverse .
In practice we restrict the computational cost on feature maps by parametrizing a local triangular matrix. That is, the transformation can be computed in parallel spatially, and is defined as: , where denote spatial coordinates, and are conditioned on , and denotes the number of channels in . Since the dimensions of are small, relative to the neural networks parametrizing them, the inverse can be found in iterations using spatially parallelized matrix operations.
Appendix C Quantizing a Continuous Flow
To test the lossless compression performance of continuous flows, the latent space is quantized to a linear spaced bins. Because the latent space is quantized, the reconstructions may contain errors. To enable lossless compression, FLIF is used to encode the errors in reconstruction. Hence, given the quantized latent variables and the reconstruction errors, the original input can be obtained.
The performance of the quantized flow is shown in Figure 10. When the bin size is large (), encoding the latent representation requires relatively few bits, because the probability area is larger. However, the residuals are higher, and require more bits to be modelled. Analogously, when the bin size is small (), encoding the latent representation requires more bits, but the residual can be modelled using fewer bits. Although the bits required for the residual or the quantized latents may be small individually, their sum is always large. In total the quantized flow performs poorly on lossless compression.
Appendix D Experimental details
d.1 Networks
The coupling and factor out layers are parametrized using neural networks. These networks are DenseNets [17]. Specifically we use intermediate channels and a depth . In contrast with standard DenseNets, we do not use normalization layers. A single layer in the densenet consists of:
d.2 IDF architecture
The exact architecture for experiments is specified in Table 4. All models are trained using Adamax [19] with standard parameters. Furthermore, the learning rate is computed as:
{epoch}
. We follow the preprocessing procedure for CIFAR10 as described in [22]. For ImageNet32 and ImageNet64, we do use additional preprocessing. For the ER + BCa dataset, we employ random horizontal and vertical flips during training.Dataset  densenet depth  densenet channels  batchsize  patchsize  train examples  lr decay  epochs  

CIFAR10  3  8  12  512  256  32  40000  0.999  2000 
ImageNet32  3  8  12  512  256  32  1230000  0.99  100 
ImageNet64  4  8  12  512  64  64  1230000  0.99  20 
ER + BCa  4  8  12  512  50  80  114  0.99999  50000 
In our implementation, instead of using integers in , we use the equivalent representation , which we found to work better with standard weight initialization and optimization methods. Despite the fact that this implementation does not use integers, it is functionally equivalent to the method presented in the main text.
d.3 Dataset preparation
The dataset for CIFAR10 originally consists of 50000 train images and 10000 test images. We use the last 10000 images for validation which results in 40000 train, 10000 validation and 10000 test images. ImageNet32 and ImageNet64 originally contain approximately 1250000 train and 50000 validation images. The validation images are used solely for testing, and 20000 images are randomly selected as a new validation set. This results in roughly 1230000 train, 20000 validation and 50000 test images.
The ER + BCa dataset [18] ^{1}^{1}1http://andrewjanowczyk.com/wpstatic/nuclei.tgz is split into train images and test images such that specific patients IDs only occur in one of the two sets. The test patient identifiers are:
8915 8959 9023 9081 9256 9382 10264 10301 12749 16532 12818 12871 12884 12908 12931 12949 13106 13459 13459 13617 13694 14154 14305 16661 17117 17643 25289 25617