Title: COMPRESSION WITHOUT QUANTIZATION

Abstract: Standard compression algorithms work by mapping an image to discrete code using an encoder from which the original image can be reconstructed through a decoder. This process, due to the quantization step, is inherently non-differentiable so these algorithms must rely on approximate methods to train the encoder and decoder end-to-end. In this paper, we present an innovative framework for lossy image compression which is able to circumvent the quantization step by relying on a non-deterministic compression codec. The decoder maps the input image to a distribution in continuous space from which a sample can be encoded with expected code length being the relative entropy to the encoding distribution, i.e. it is bitsback efficient. The result is a principled, end-to-end differentiable compression framework that can be straight-forwardly trained using standard gradient-based optimizers. To showcase the efficiency of our method, we apply it to lossy image compression by training Probabilistic Ladder Networks (PLNs) on the CLIC 2018 dataset and show that their rate-distortion curves on the Kodak dataset are competitive with the state-of-the-art on low bitrates.

1. INTRODUCTION:

The recent development of powerful generative models, such as Variational Auto-Encoders (VAEs) and their hierarchical extensions, such as Probabilistic Ladder Networks (PLNs) (Kingma & Welling, 2014; Sønderby et al., 2016; Higgins et al., 2017) has caused a great deal of interest in their application to lossy compression, notably Ballé et al. (2016) ; Theis et al. (2017) ; Rippel & Bourdev (2017) ; Ballé et al. (2018) ; Mentzer et al. (2018) ; Johnston et al. (2018) . The benefit of using these models as opposed to hand-crafted methods is that they can adapt to the statistics of their inputs much better, and hence allow significant gains in compression rate at a given quality setting, and in quality at a given compression rate. A second advantage is their easier adaptability to new media formats, such as light-field cameras, 360 • images, Virtual Reality (VR), video streaming, etc. for which classical methods are not currently applicable, or would perform pathologically badly.

Lossy compression codecs usually perform a lossy transformation on the continuous representation of the data, after which it is mapped to a discrete domain through an operation known as quantization, such that entropy coding may be used to encode this transformed representation. This continuous-to-discrete mapping will mean that quantization will necessarily have zero derivative almost everywhere. Since the above mentioned generative models are usually trained using some gradient-based technique, this means that they cannot be directly applied as compression models, as quantization will destroy the learning signal, and hence a workaround is needed. A nice remedy to this problem is to replace the quantities affected by quantization (or at least their derivatives) by some smooth approximation during training, as is done in recent approaches, e.g. Ballé et al. (2016) ; Theis et al. (2017) ; Ballé et al. (2018) ; Mentzer et al. (2018) ; Johnston et al. (2018) .

Similarly to Ballé et al. (2016) ; Theis et al. (2017) ; Ballé et al. (2018) , we aim to minimize the rate-distortion of our model directly, but unlike them we replace entropy coding with a more general coding technique that allows us to use probability densities for coding instead of probability masses. This then allows us to forgo the quantization operation altogether and we get an end-to-end differentiable model that can be optimized with standard gradient-based optimizers.

We select an appropriate generative model p(x, z) = p(x | z)p(z) with an approximate posterior q(z | x), where the latents z are going to serve as the transformed representation of x that we wish to compress. We show that -log p(x | z) corresponds to the distortion and KL [ q(z | x) || p(z) ] corresponds to the rate of the model. Hence, optimizing our model for the rate-distortion is as simple as minimizing -log p(

x | z) + βKL [ q(z | x) || p(z) ],

where β is a scalar hyperparameter controlling the compression rate.

To compress some data x with our trained model, instead of selecting z deterministically and then coding it, we map x to the posterior q(z | x) and use it to code a sample from it, using an adapted version of the importance sampling algorithm proposed by Havasi et al. (2018) . First, set K = KL [ q(z | x) || p(z) ] and stochastically draw independent samples (z 1 , . . . , z exp{K} ) from p(z). Next, select z c * for some 1 ≤ c * ≤ exp{K}. We can use the index c * as the compressed representation of z c * . Clearly, c * can be coded in K nats, and it can also be shown that z c * will be a low-bias sample from q(z | x). As the efficiency of this algorithm is KL [ q(z | x) || p(z) ], we will refer to it as relative entropy coding (REC), and note that REC algorithms achieve bits-back efficiency (Hinton & Van Camp, 1993) .

The contributions of this paper are as follows:

• We show how any generative model p(x, z) = p(x | z)p(z) with approximate latent posterior q(z | x) can be used for bits-back efficient lossy transform coding, using an adapted version of Minimal Random Code Learning (MIRACLE, Havasi et al. (2018) ).

• We present the concept of relative entropy coding (REC) algorithms and attainable bounds on their worst-case performance. We note their current intractability in high-dimensional settings and propose an approximate solution that is tractable, though the rate-distortion quality is impacted.

• We demonstrate the effectiveness by applying our method to lossy image compression.

Concretely, we train a Probabilistic Ladder Network (Sønderby et al., 2016) on the CLIC (2018) dataset and show that the rate-distortion curve of our method is competitive with the state-of-the-art on the Kodak dataset (Eastman Kodak Company, 1999) on low bitrates.

2. LITERATURE REVIEW:



2.1. TRANSFORM CODING:

Neural network-based approaches lend themselves easiest to transform coding as a realisation of lossy compression. In transform coding (Goyal, 2001) , given some data x, it is first transformed using the analysis transform z = f a (x) and then quantized ẑ = [z], where [•] could denote rounding, or any other quantizer in general. The distribution of this quantized representation is modelled by some probability mass function P (ẑ), which is used to encode ẑ using some entropy code. Note, that this is different from the true distribution of ẑ, P (ẑ), which is entirely determined by the distribution of [f a (x)], where x ∼ p(x). During decompression, once ẑ is decoded, the synthesis transform x = f s (z) is applied to it to obtain a reconstruction of x.

Given some metric d(•, •), the distortion of the codec is

EQUATION

and its rate is the expected codelength of the entropy code, which is

EQUATION

The equality holds, since ẑ is a deterministic transformation of x.

2.2. BITS-BACK EFFICIENCY:

If instead of mapping some data we wish to compress x to a single deterministic point z and entropy coding it by using some prior p(z), we map it to an approximate posterior q(z | x), then we may select a random representative z * from q(z | x) in the sense that z * has a short Minimum Description Length (MDL, Grünwald et al. (2007) ). Bits-back efficient methods are such that this MDL is approximately KL [ q(z | x) || p(z ]. An important distinction should be made in whether we wish to communicate the distribution q(z | x) or a sample z * drawn from it, as the communication problem from which the MDL can be obtained, its precise value and the method required to obtain it are different.

The concept originates in the work of Hinton & Van Camp (1993) . Given data D and model weights w, they derive the MDL of the posterior weight distribution q(w | D) of a Bayesian Neural Network (BNN) with prior p(w) to be K = KL [ q(w | D) || p(w) ]. However, they did not give a method for achieving this, in fact, in their argument they merely show that in a longer message only K nats are used to code the posterior.

In our setting, we are only interested in communicating a single sample z * from q(z | x). We assume that we have access to a generative model p(x | z)p(z) and a posterior q(z | x). Then, the problem we must answer is what is the minimum expected number of nats T [x : z] that is required to be communicated, such that we could draw a sample z ∼ q(z | x). This question was answered by Harsha et al. (2007) , who proved the following theorem:

Theorem 1 (Harsha et al., 2007) Let the setting be as described above, and assume that the coder and the decoder share a common source of randomness (e.g. by using the same pseudo-random number generator with the same seed), then the minimum expected number of nats T [x : z] that need to be communicated in order for the decoder to be able to sample a z ∼ q(z | x) can be bounded in terms of the mutual information I[x : z] as

EQUATION

They prove the above claim constructively, by giving a rejection sampler that achieves this. This, to our knowledge, is the first example of what we call a single-shot relative entropy coding (REC) algorithm, and we will examine it in further detail in Section 4. Bits-back efficiency comes from the fact that

I[x : z] = E x∼p(x) [KL [ q(z | x) || p(z) ]].

Relative entropy coding subsumes entropy coding in the sense that any entropy code is a relative entropy code as well, since for a discrete ẑ ∼ P (ẑ) and the point-mass δ ẑ on ẑ,

EQUATION

See Havasi et al. (2018) for details on the derivation. In this sense, methods that use quantization and then entropy coding are also bits-back efficient, but are restricted to point-masses for their latent posteriors. Thus, one of the advantages of our method lies in that it allows a much wider class of posteriors to be used.

2.3. RELATED WORKS:

Minimal Random Code Learning (Havasi et al., 2018) is a powerful coding technique that allows us to code a random sample from a continuous distribution, rather than coding a deterministic sample, which requires quantization, as any fixed sample would have 0 density. They use their technique to compress the weights of BNNs. The structure of the weight space is much different from the structure of the latent space of a generative model trained for the compression of some medium. Hence, we argue that their practical compression procedure is infeasible for our setting, the random fixed-size block partitioning as well as the retraining between the coding of blocks is incompatible with our setting as well as computationally infeasible. Townsend et al. (2019) were the first to develop a bits-back efficient image compression algorithm using VAEs, called Bits-back with ANS (BB-ANS). However, they develop a lossless compression algorithm, and hence the performance their approach is not comparable to ours. The sense in which their method is "bits-back" also differs from ours. Townsend et al. (2019) stick to the original bitsback setting. Concretely, they show a clever way of utilizing the overhead of the method of Hinton & Van Camp (1993) for efficient coding, our work on the other hand, avoids sending it in the first place. Their coding efficiency collapses onto entropy efficiency when coding a single image and bits-back efficiency is only achieved when coding larger batches.

Lately, there has been an explosion of interest in using powerful deep generative models for lossy image compression, and thus all previous methods had to deal with the issue of the non-differentiability of quantization. These approaches can generally be put into three different categories, depending on how they avoid this problem. There are approaches that waive the paradigm of end-to-end training, e.g. by not incorporating the feedback from the gradient of the rate term into the optimization of the entropy coder, either by using inflexible distributions without any parameters (Toderici et al., 2017) , or by optimizing them post-training (Johnston et al., 2018; Rippel & Bourdev, 2017) . The issue with this, as pointed out in Ballé et al. (2018) is that not training end-to-end might severly limit the efficiency of these methods.

Another approach of getting aroud this issue has been to learn to generate samples from the image distribution p(x), using a Generative Adversarial Network (GAN, Goodfellow et al. (2014) ), by learning the synthesis transform first and using an auxiliary discriminator. Then, fixing the discriminator, the analysis transform is trained in order to minimize the reconstruction error d(x, f s (f a (x))).

The noise provided f s during the initial GAN training should be a continous relaxation of quantization noise, such as a uniform distribution (Ballé et al., 2016) .

Finally 2017)), though they only derive it for a smaller class of distributions than we do, as well as through a different route, as we arrive at it through the MDL principle (Grünwald et al., 2007) . Crucially, these methods are still only equivalent to VAEs during training, and the continuous relaxations are switched back to the original discrete operations thereafter.

3. COMPRESSION WITHOUT QUANTIZATION:

Our method relies on two observations. First, we note that the sole purpose of quantization is that once we have obtained a deterministic representation ẑ, it has non-zero probability mass, and hence entropy coding may be used to compress it. Hence, if we could replace entropy coding with a method that could use a probability density p(z) instead of the mass P (ẑ), then we could forgo the quantization step and thus our training objective becomes end-to-end differentiable.

Second, as we have seen in Section 2.2, if in addition have access to an approximate latent posterior q(z | x), then it is enough to communicate approximately

EQUATION

nats such that the decoder may obtain a sample z ∼ q(z | x). Applying f s to z could then be used to reconstruct x.

These two could be combined by noting that our setting fulfils all premises required for both observations to hold. This enables the following development of a lossy transform codec:

Training After fixing an appropriate architecture for f a (x) = q(z | x) and f s such that the likelihood is p(x = f s (z) | z), we train our model using the β-Evidence Lower Bound (Higgins et al., 2017) :

EQUATION

Note, that if the distortion metric d(•, •) is such, that exp{-d(x, x )} can be normalized, then the negative log-likelihood can be identified with the distortion D.

E x∼p(x) [KL [ q(z | x) || p(z) ]] = I[x : z]

can also be identified with the rate R. Hence, in this case the negative β-ELBO is equivalent to training for the rate-distortion loss

EQUATION

which was our aim to train for in the first place. This gives us a new way to interpret β in Equation 6: it allows the model to converge on different points on its rate-distortion curve (see Ballé et al. (2016) for the full argument).

Note, that even if normalizing exp{-d(x, x )} is intractable β can be adjusted in Equation 6to compensate for this unknown constant during training.

Compression and Decompression:

For given x that we wish to compress, we calculate q(z | x) using the analysis transform. Then, we use relative entropy coding (REC) to code a sample z ∼ q(z | x). The decoder can then re-obtain z and apply the synthesis transform to it to reconstruct x.

4. RELATIVE ENTROPY CODING:

As mentioned in the previous section, our method relies on the existence of relative entropy coding (REC) in order to be able to code the latent representation of a given input x with bits-back efficiency. To our knowledge, the first such algorithm was developed by Harsha et al. (2007) , as a variant of rejection sampling. This algorithm is unfortunately intractable in more than a few dimensions, as it requires to keep track of probabilities across the entire range of ẑ.

The simplest way to extend the univariate algorithm to the multivariate setting is to note the independence assumption for the dimensions of ẑ and thus concatenating univariate samples for each dimension will given an exact multivariate sample. Hence, we could use REC to encode the dimensions of ẑ individually, i.e. for each ẑi obtain some code c i , and build an empirical distribution P (c) over these. Then, we can encode each c i using entropy coding. The issue with this approach is that while the upper bound in Theorem 1 is quite tight if we apply it to the whole multivariate distribution over ẑ, when applied to each of its dimensions ẑi , the O(1) cost is incurred for each dimension, which makes this approach very inefficient. Hence, we focus on multivariate REC.

4.1. ADAPTIVELY GROUPED IMPORTANCE SAMPLING:

Algorithm 1 Importance sampling algorithm proposed by Havasi et al. (2018) Inputs:

P -Proposal distribution Q -Target Distribution x i ∼ Q | i ∈ N -Shared sequence of random draws from Q procedure IMPORTANCE-SAMPLER(P, Q, x i ∼ Q | i ∈ N ) K ← exp{KL [ Q || P ]} wi ← log Q(xi) P (xi)

∀i = 1, . . . K j ← arg max i wi return j, x j end procedure

The basis of our approach is the approximate importance sampler proposed by Havasi et al. (2018) . For given x, it proceeds by letting K = KL [ q(z | x) || p(z) ] and drawing S = e K samples (z 1 , z 2 , . . . , z S ) from p(z). Then the index of the maximal importance weight

c * = arg max c∈{1,2,...S} q(z c | x) p(z c ) (8)

is used as the encoding of x. The decoder can recover z c * by drawing c * samples from p(z) and then can use it to reconstruct x, assuming it can generate the same sequence of samples. Clearly, by virtue of the choice for S, c * can be coded in KL [ q(z | x) || p(z) ] nats, and Havasi et al. (2018) also show that z c * has negligible bias. The algorithm is depicted in Algorithm 1, note that Havasi et al. (2018) originally sample the softmax distribution of the wi . We found it better to select samples biased towards the mode by always selecting the most likely sample. Though now we do not need to keep track of the probabilities over the entire range of z, the run-time is still going to be exponential in K, and hence becomes prohibitive for K larger than 15 -20. In our experiments we had to deal with KLs well above 10 5 , and hence solving this issue is critical for the viability of our method.

To solve this issue, Havasi et al. (2018) randomly partition the vector of interest w into fixed size blocks, in the hope that on average each individual block will have approximately some C nats of total KL. They then code each the blocks sequentially. The concern is that the next block to be coded is a suboptimal choice, i.e. its total KL is significantly higher or lower than the desired C nats. To circumvent this problem, they introduce a blockwise KL penalty to the loss that forces each block to have equal KL, and they intersperse coding a block with several rounds of retraining. They show that this solves the suboptimal choice problem, and their coding procedure also does not impact the quality of the coded w.

Unfortunately, this fix is not applicable in our case. Havasi et al. (2018) can use the retraining trick to compensate for suboptimal block choices. However, if we hope to obtain a useful compression codec, we should avoid extremely expensive operations to code a single item, such as performing gradient descent for each individual input x to our algorithm.

To remedy these issues, instead of selecting fixed-size blocks randomly we impose a group size constraint of G bits and total group KL constraint of B bits. Then, we partition the dimensions of z such that each set in the partition (what we call a group) has size at most 2 G and each group has total KL at most K bits. The details are depicted in Algorithm 2.

Algorithm 2 Adaptively grouped importance sampler Inputs:

K -Maximum individual KL allowed G -Maximum group size B -Bit budget per group P -Proposal distribution Q -Target Distribution x i ∼ Q | i ∈ N -Shared sequence of random draws from Q procedure ADAPTIVE-IMPORTANCE-SAMPLER(K, G, B, P, Q, x i ∼ Q | i ∈ N ) Γ ← () Initialize empty list of group sizes k i ← KL [ Q i || P i ] ∀i = 1, . . . N Get KLs for each dimension Q , P , O ← REMOVE-OUTLIERS({k i } N i=1 , K) γ ← 0, k ← 0

Current group size and group KL

for i ← 1, . . . dim(Q ) do if k + kl i > B or γ + 1 > G then Append γ to Γ k ← kl i , γ ← 1 else k ← k + kl i , γ ← γ + 1 end if end for Append γ to Γ

Append the last group size S = (), I = (), g ← 0 Samples, sample codes and current group index for γ in Γ do Now importance sample each group i, s ← IMPORTANCE-SAMPLER(P g:g+γ , Q

g:g+γ , x i ∼ Q | i ∈ N ) Append i to I, append s to S g ← g + γ end for return I, S, O end procedure

While the concern of intractability due to high KL is resolved this way but we are faced with two inefficiencies. First, for groups whose size is much smaller than 2 G , dedicating G bits to coding the group size is wasteful. Second, many groups will have size 2 G but total KL much less then K, and hence dedicating K bits to code the sample is wasteful. To solve these issues, we built an empirical distribution over group sizes using our training dataset, which is tractable since we chose G to be not too large (2 ≤ G ≤ 6), and then used arithmetic coding (Rissanen & Langdon, 1981) to code the group sizes. This worked very effectively, to the extent that for appropriately selected G the cost of communicating the group sizes is effectively negligible. where the W × H × C/S implies a W × H convolution kernel, with C target channels and S gives the downsampling rate (given a preceding letter "d") or the upsampling rate (given a preceding letter "u"). If the slash is missing, it means that there is no up/downsampling. All convolutions operate in same mode with zero-padding. GDN / IGDN: these are the non-linearities described in Ballé et al. (2016) . Leaky ReLU: elementwise non-linearity defined as max{x, αx}, where we set α = 0.2. Sigmoid: Elementwise non-linearity defined as The relative entropy codes of the groups were also further compressed using Elias δ-coding (Elias, 1975) . This is the main point of inefficiency of this method, as the empirical distribution of sample indices has a much heavier tail than the distribution implied by δ-coding. Concretely, we found that using δ-coding the resulting code lengths are approximately 1.6 -1.8 times the theoretically possible REC code lengths. However, building an empirical distribution over the sample indices and using it for entropy coding is more challenging, as the number of possible indices for a "useful" K is between 2 15 ≤ 2 K ≤ 2 20 . This means that to obtain a good estimate of the tail mass, i.e. of indices above 2 14 -2 15 , we need a lot of data, or some reasonable smoothing should be used, and we leave this for future work. is displayed for each dimension, where the µ i and σ i are the means and standard deviations of the dimensions of p(z (1) | z (2) ), respectively. We see that relative to the conditional prior, the latent structure is effectively captured.

5. EXPERIMENTS:

We applied our framework to lossy image compression. We followed the work of Ballé et al. (2016) and Ballé et al. (2018) . Not only do their results represent the current state-of-the-art, but their work is also the closest to ours. In particular, as mentioned before, to our knowledge they were the first to derive the connection of the weighted rate-distortion loss with the β-ELBO as the training objective for a VAE. However, we stress that they derived this for a VAE whose latent distribution was a continuous relaxation of the quantization error. Therefore, their results were restricted to a narrow set of distributions, and neither did they make the connection to bits-back efficientcy. Most notably though, they only used the continuous relaxation during the training of the model, thereafter switching back to quantization and entropy coding, which, they show does not impact the predicted performance, and hence confirming that their relaxation during training is reasonable. This is in striking contrast with our method, in which not only do we allow any latent distributions for which the KL can be calculated, but our architecture also remains the same after training.

Nonetheless, we based our architecture entirely on the one they present in Ballé et al. (2018) , appropriately modified for our setting. Most importantly, where they use i.i.d. uniform distributions for the posterior and complicated non-parametric priors, we opted for Gaussians, a much more standard choice for VAEs. We also opted for a similar hierarchical VAE structure, which in our case is a Probabilistic Ladder Network (PLN, Sønderby et al. ( 2016)). The reason for this is as shown by Ballé et al. (2018) and also confirmed by us empirically, the full independence assumption in the latent distribution of a VAE is very limiting, as there is a lot of topological structure even on the latent level. One solution to introduce a rich dependence structure is to introduce further latent variables, conditioned upon which the rest are assumed to be independent (Bishop, 1998) . Concretely, we break up our generative model p(x, z) = p(x | z (1) )p(z (1) | z (2) )p(z (2) ), and our approximate latent posteriors now become q(z (2) | x) and q(z (1) | x, z (2) ) ∝ q(z (1) | x)p(z (1) | z (2) ). Further, our architecture also only uses convolutions and deconvolutions as non-linearities. This allows us to compress arbitrary sized images, and the latent space can naturally grow with the dimensions of the images. Full details of the architecture can be seen in Figure 1 . A confirmation that the PLN's conditional independence structure efficiently removes the topological structure in the latent space on the first stochastic level is shown in Figure 2 . Note, that in order for us to code a latent sample z ∼ q(z | x), we first code z (2) ∼ q(z (2) | x) using p(z (2) ) and then code z (1) | z (2) ∼ q(z (1) | x, z (2) ).

In our experiments, we compare the performance of our method with the method of Ballé et al. (2018) and JPEG as a baseline on their rate-distortion curves, using the Peak Signal-to-Noise Ratio (PSNR) and Multiscale Structural Similarity Index (MS-SSIM) perceptual metrics (Huynh-Thu & Ghanbari, 2008; Wang et al., 2003) as distortion distances. Our models were trained using β ∈ {300, 600, 1200, 2400, 4800} and the Mean Absolute Error (MAE) loss. The reason for the MAE loss is due to Zhao et al. (2015) , who show that using it greatly improves image reconstruction quality at no significant extra cost. The model of Ballé et al. (2018) 

:



Figures and Tables:

Figure 1: PLN network architecture. The blocks signal data transformations, the arrows signal the flow of information. Block descriptions: Conv2D: 2D convolutions along the spatial dimensions,where the W × H × C/S implies a W × H convolution kernel, with C target channels and S gives the downsampling rate (given a preceding letter "d") or the upsampling rate (given a preceding letter "u"). If the slash is missing, it means that there is no up/downsampling. All convolutions operate in same mode with zero-padding. GDN / IGDN: these are the non-linearities described inBallé et al. (2016). Leaky ReLU: elementwise non-linearity defined as max{x, αx}, where we set α = 0.2. Sigmoid: Elementwise non-linearity defined as

. We ran all experiments presented here with N = 196, M = 128, F = 128, G = 24.

Figure 2: Demonstration of the effectiveness of the PLN's conditional independence assumption p(z) = p(z (1) | z (2) )p(z (2) ). a) Original image from the Kodak dataset. b) Randomly selected channel of a latent sample z (1) ∼ q(z (1) | z (2) , x), clearly there is a lot of structure present in the sample. c) The same channel of the latent sample standardized according to p(z (1) | z (2) ), i.e. we z (1) i -µi σi

Figure 7: Rate-distortion curves for kodim05

Figure 8: Rate-distortion curves for kodim08

, there are methods that stick to Variational Auto-Encoders (VAEs,Kingma & Welling (2014);Higgins et al. (2017)), optimize the weighted rate-distortion directly, and replace the nondifferentiable quantities with smooth approximations during training.Theis et al. (2017) introduce compressive auto-encoders, and replace the derivative of the quantization operation with the identity function, and the rate of the discrete ẑ by a differentiable upper bound which is used in the loss. Closest to our work are the works ofBallé et al. (2016; 2018);Minnen et al. (2018), who all model the quantization noise as uniform noise centered on the quantized representation.Ballé et al. (2016) are the first to derive the connection between the weighted rate-distortion loss and the β-Evidence Lower Bound (β-ELBO,Higgins et al. (

Table content: <table/>

Peer review (iclr2020):
The authors propose a new image compression method that does not require quantizing the encoded bits in an auto-encoding style image compression model. The method builds on a VAE with image x, code z, and posterior q(z | x). Instead of directly encoding z, the proposed method samples z_{c^star} from posterior and store c^star as the compressed representation. A decoder then reconstruct an approximation of x from z_{c^star}. The authors show that this framework is bits-back efficient and draw connections to prior theoretical results. Experiments were conducted on the Kodak dataset based on the model of Balle et al. The proposed method works only slightly worse than Balle et al. at low bit-rate region, but the gap becomes larger in higher bit rate regions. The method is technically sound and the paper is clearly written. My main concerns fall in practical aspects. In Figure 3, for 3 out of 4 images, the theoretical upper bound of the proposed method still do not outperform Balle et al. This suggests limitations of the proposed method. Discussion on the limitations of the method is limited. The results also beg the question: How much does quantization in existing methods impact performance, and how much will fixing this benefit the overall system. Finally, in my opinions, in some sense the sampling of z_{c^star} is also a form of quantization. Does drawing different samples from posterior leads to different reconstructed images? If it does, doesn*t it also suffer from similar limitations as existing *quantization* methods? Overall, I think the direction of the proposed method has good potential, but it also leaves important questions unanswered. I think this paper will benefit from additional revisions.