Towards Empirical Sandwich Bounds on the Rate-Distortion Function
Authors: Yibo Yang
Stephan Mandt
Venue: ICLR 2022 Poster (2022)
Keywords: information theory
deep generative modeling
lossy data compression

Abstract: Rate-distortion (R-D) function, a key quantity in information theory, characterizes the fundamental limit of how much a data source can be compressed subject to a fidelity criterion, by any compression algorithm. As researchers push for ever-improving compression performance, establishing the R-D function of a given data source is not only of scientific interest, but also reveals the possible room for improvement in existing compression algorithms. Previous work on this problem relied on distributional assumptions on the data source (Gibson, 2017) or only applied to discrete data (Blahut, 1972; Arimoto, 1972). By contrast, this paper makes the first attempt at an algorithm for sandwiching the R-D function of a general (not necessarily discrete) source requiring only i.i.d. data samples. We estimate R-D sandwich bounds for a variety of artificial and real-world data sources, in settings far beyond the feasibility of any known method, and shed light on the optimality of neural data compression (Ballé et al., 2021; Yang et al., 2022). Our R-D upper bound on natural images indicates theoretical room for improving state-of-the-art image compression methods by at least one dB in PSNR at various bitrates. Our data and code can be found at https://github.com/mandt-lab/RD-sandwich.

Published as a conference paper at ICLR 2022
TOWARDS EMPIRICAL SANDWICH BOUNDS ON THE RATE-DISTORTION FUNCTION
Yibo Yang, Stephan Mandt
Department of Computer Science, UC Irvine (yibo. yang, mandt @unddt@@uuc.eem edu
ABSTRACT
Rate-distortion (R-D) function, key quantity information theory, characterizes the fundamental limit of how much data source can be compressed subject to a fidelity criterion, by any compression algorithm. As researchers push for ever- improving compression performance, establishing the R-D function of given data source not only of scientific interest, but also reveals the possible room for improvement in compression algorithms. Previous work this problem relied on distributional assumptions on the data source (Gibson, 2017) or only applied discrete data. By contrast, this paper makes the first attempt at algorithm for sandwiching the R-D function of a general (not necessarily discrete) source requiring only i.i.d. data samples. We estimate R-D sandwich bounds for variety of artificial and real-worid data sources, in settings far beyond the feasibility of any known method, and shed light on the optimality of neural data compression (Ballé et al., 2021; Yang et al., 2022). Our R-D upper bound on natural images indicates theoretical room for improving state-of mfte-off-the-ar image compression methods by at least one dB in PSNR at various bitrates. Our data and code can be found here.
INTRODUCTION
From storing astronomical images captured by the Hubble telescope. to delivering familiar faces and voices over video calls, data compression, i.c., communication of the "same" information but with less bits, is commonplace and indispensable to our digital life, and even arguably lies at the heart of intelligence (Mahoney, 2009). While for lossless compression, there exist practical algorithms that can compress any discrete data arbitrarily close to the information theory limit (Ziv & Lempel, 1977; Witten et al., 1987), no such universal algorithm has been found for lossy data compression (Berger & Gibson, 1998), and significant research efforts have dedicated to lossy compression algorithms for various data. Recently, deep learning has shown promise for learning lossy compressors from raw data examples, with continually improving compression performance often matching exceeding traditionally engineered methods (Minnen et al., 2018; Agustsson et al., 2020; Yang et al., 2020a).
However, there are fundamental limits to the performance of any lossy compression algorithm, due to the inevitable trade-off between rate, the average number of bits needed to represent the data, and the distortion incurred by lossy representations. This trade-off is formally described by the rate-distortion (R-D) function, for a given source (i.e., the data distribution of interest; referred such in information theory) and distortion metric. The R-D function characterizes the best theoretically achievable rate- dise-iittoiion performance by any compression algorithm, which can be seen as compressooo counterpart and generalization of Shannon entropy lossless compression.
Despite its fundamental importance, the R-D function is generally unknown analytically, and estab- lishing it for general data sources, especially real world data, a difficult problem (Gibson, 2017). The default method for computing R-D functions, the Blahut-Arimote algorithm (Blahut, 1972; Arimoto, 1972), only works for discrete data with a known probability mass function and has complexity exponential in the data dimensionality. Applying it to unknown data source requires discretization (if it is continuous) and estimating the source probabilities by a histogram, both of which introduce errors and computationally infeasible beyond couple of dimensions. Previous work characterizing the R-D function of images and videos (Hayes et al., 1970; Gibson, 2017) all assumed statistical model of the source, making the results dependent the modeling assumptions.

---

Published as a conference paper at ICLR 2022
In this work, we make progress this long-standing problem in information theory using tools from machine learning, and introduce new algorithms for upper and lower bounding the R-D function of a general (i.e., discrete, continuous, neither), unknown memoryless source. More specifically,
1. Similarly to how a VAE with a discrete likelihood model minimizes an upper bound on the data entropy, we establish that any 3- VAE with a likelihood model induced by a distortion metric minimizes an upper bound on the data rate-distortion function. We thus open the deep generative modeling toolbox to the estimation of an upper bound on the R-D function.
2. We derive a lower bound estimator of the R-D function that can be made asymptotically exact and optimized by stochastic gradient ascent. Facing the difficulty of the problem involving global optimization, we restrict to a squared error distortion for a practical implementation.
3. We perform extensive experiments and obtain non-trivial sandwich bounds on various data sources, including H'N-ggenerate artificial sources and real-world data from speech and physics. Our results shed light on the effectiveness of neural compression approaches (Ballé et al., 2021; Minnen et al., 2018; Minnen & Singh, 2020), and identify the intrinsic (rather than nominal) dimension of data as a key factor affecting the tightness of our lower bound.
4. Our estimated R-D upper bounds on high-resolution natural images (evaluated on the standard Kodak and Tecnick datasets) indicate theoretical room for improvement of state-of- the-art image compression methods by at least one dB in PSNR, at various bitrates.
We begin by reviewing the prerequisite rate-distortion theory in Section 2, then describe our upper and lower bound algorithms in Section and Section 4, respectively. We diseuss related work in Section 5, report experimental results in Section 6, and conclude in Section 7.
2 BACKGROUND
(R-D) theory deals with the fundamental trade-ofi off between the average number of bits per sample (rate) used to represent data source X and the distortion incurred by the lossy representation Y. It asks the following question about the limit of lossy compression: for a given data source and a distortion metric (a.k.a., a fidelity criterion), what is the minimum number of bits (per sample) needed to represent the source at a tolerable level of distortion, regardless of the computation complexity of the compression procedure? The answer is given by the rate-distortion function R(D). To introduce it, let the source and its reproduction take values in the sets X' and y, conventionally called the source and reproduction alphabets, respectively. We define the data source formally by random variable X X following a (usually unknown) distribution Px, and assume distortion metric : X' x y + [0, oo) has been given, such as the squared error p(æ,, y) = || - yll². The rate-distortion function is then defined by the following constrained optimization problem.
inf R(D) Qvix: Elp() ((;Y),
(1)
where we consider all random transforms Qrix whose expected distortion is within the given threshold D > O, and minimize the mutual information between the source X and its reproduced representation Y Shannon's 's lossy source coding theorems (Shannon, 1948; 1959) gave operational significance to the above mathematical definition of R(D), as the minimum achievable rate with which any lossy compression algorithm can code i.i.d. data samples at a distortion level within D.
The R-D function thus gives the tightest lower bound the rate-distortion performance of any lossy compression algorithm, and can inform the design and analysis of such algorithms. If the operational distortion-rate performance of an algorithm lies high above the source R(D)-eurve (D, R(D)), then further performance improvement may be expected; otherwise, its performance is already close to theoretically optimal, and we may focus our attention on other aspects of the algorithm. As the R-D function does not have an analytical form in general, we propose to estimate it from data samples, making the standard assumption that various expectations w.r.t. the true data distribution Px exist and can be approximated by sample averages. When the source alphabet is finite, this assumption automatically holds, and R(D) also provides lower bound on the Shannon entropy of discrete data.
Both the expected distortion and mutual information terms are defined w.r.t. the joint distribution Px Qx We formally describe the general setting of the paper, including the technical definitions, in Appendix A.I.
2

---

Published as a conference paper at ICLR 2022
3 UPPER BOUND ALGORITHM
R-D theory (Cover & Thomas, 2006) tells us that every (distortion, rate) pair lying above the R(D)- curve is in theory realizable by (possibly expensive) compression algorithm. An upper bound on R(D) thus reveals what kind of compression performance is theoretically achievable. Towards this goal, we borrow the variational principle of the Blahut-Arimoto (BA) algorithm, but extend it to general (e.g.. non-discrete) source requiring only its samples. Our resulting algorithm optimizes a 3-VAE whose likelihood model is specified by the distortion metric, a common case being Gaussian likelihood with fixed variance. For the first time, we establish this class of models as computing model-agnostic upper bound on the source R-D function, as defined by a data compression task.
Variational Formulation. Following the BA algorithm (Blahut, 1972; Arimoto, 1972), we consider a Lagrangian relaxation of the constrained problem defining R(D), which has the variational objective
L(Qrix, Qy A) == Px I x [Qy]] + XEpx [p(X,
(2)
where Qy is new, arbitrary probability measure on y. The first (rate) term is a variational upper bound on the mutual information I(X;)) and the second (distortion) term enforces the distortion tolerance constraint in Eq. 1. For each fixed À > O, the BA algorithm globally minimizes L w.r.t. the variational distributions Qryx and Qy by coordinate descent; at convergence, the (distortion, rate) pair yields point on the R(D) curve (Csiszár, 1974a). Unfortunately, the BA algorithm only applies when X' and y are finite (hence discrete), and the source distribution known. Otherwise, a preprocessing step is required to discretize a continuous source and/or estimate source probabilities by a histogram, which introduces a non-negligible bias. This bias, along with its exponential complexity in the data dimension, also makes BA infeasible beyond a couple of (usually 3) dimensions.
Proposed Method. To avoid these difficulties, we propose to apply (stochastic) gradient descent on L. w.r.t. flexibly parameterized variational distributions Qy x and Qy- In this work we parameterize the distributions by neural networks, and predict the parameters of each Qrxx== by an encoder network (xz) as in amortized inference (Kingma & Welling, 2014). Given data samples, the estimates of rate and distortion terms of L yield point that in expectation lies on an R-D upper bound Ry (D). and we tighten this bound by optimizing C; repeating this procedure for various À traces out Ry (D).
The objective closely resembles the negative ELBO (NELBO) objective of a 3-V (Higgins et al., 2017) if we view the reproduction space ) as the "latent space". The connection is immediate when X is continuous and a squared error P specifies the density of Gaussian likelihood p(zly) x exp(-llz->) [)])) /²) However, unlike in data compression, where ) is determined by the application (and often equal to X for a full-reference distortion), the latent space in a (3-)VAE typically has a lower dimension than X, and decoder network is used parameterize a likelihood model in the data space. To capture this setup, we introduce a new, arbitrary latent space Z on which we define variational distributions Q2ix: Qz, and (possibly stochastic) decoder function w 2 + y. This results in an extended objective, resembling a 3 with likelihood density p(z|2)) x
J [Q2ixx Q,,w, X) == Ex [K L(Q21x=z11Q22) + XE lp(X, (Z))]-
(3)
This objective is closely related to the original data compression task and provides upper bound on the source R(D), as follows. Treating Z as the reproduction alphabet, we can define a new distortion P (I,, == p( w()), and w-induced R-D function, R.,. (D) = in EXE, for each choice of w. Our Theorem A.3 then guarantees that R. (D) Z R(D), for any w, and consequently the (distortion, rate) of J always lies above R ()) Moreover, (D) R(D) for bijective w, which offers some theoretical support for the use of invertible pixel-shuffe operations instead of upsampled convolutions in the decoder of image compression autoencoders (Theis et al., 2017: Cheng et al., 2020). We can now minimize the NEL BO-like objective J w.r.t. the parameters of (Qzjx, Qz,w) w) similar to training 3AVAE, knowing that we are optimizing an upper bound on the rate-distortion function of the source. This can be seen as the lossy counterpart to the lossless setting, where it is well-established that minimizing the NELBO minimizes an upper bound on the Shannon entropy of the source (Frey& & Hinton, 1997; MacKay, 2003), the limit of lossless data compression.
T(X; Z),
The extended objective offers the freedom to define variational distributions on any suitable latent space Z, rather than V, which we found to simplify the modeling task and yield tighter bounds. E.g.. even if y is high-dimensional and discrete, we can still work with densities continuous and
3

---

Published as a conference paper at ICLR 2022
lower-dimensional 2 and draw upon tools such as normalizing flows (Kobyzev et al., 2021). We can also treat as the concatenation of sub-vectors [Z1, .2.,.. ZL}, and parameterize Qz in terms of simpler component distributions Qz = [ Q2ize (similarly for Qz x) as in a hierarchical VAE.
4 LOWER BOUND ALGORITHM
Without knowing the tightness of an R-D upper bound, we could be wasting time and resources trying to improve the R-D performance of compression algorithm, when it in fact already close to the theoretical limit. This would be avoided if we could find a matching lower bound on R(D). Unfortunately, the problem turns out to be much more difficult computationally. Indeed, every compression algorithm, or every pair of variational distributions (Qy. Qy1x) yields point above R(D). Conversely, establishing lower bound requires disproving the existence of any compression algorithm that can conceivably operate below the R(D) curve. this section, we derive an algorithm that can in theory produce arbitrarily tight R-D lower bounds. However, as an indication of its difficulty, the problem requires globally maximizing family of partition functions. By restricting to a continuous reproduction alphabet and squared error distortion, we make some progress on this problem and obtain useful lower bounds especially on data with low intrinsic dimension (see Sec. 6).
Dual characterization of R(D). While upper bounds on R(D) arise naturally out of its definition as a min- imization problem, variational lower bound requires expressing R(D) through a maximization problem. For this, we introduce "dual" function the optimum of the Lagrangian Eq. 2 (Qy is eliminated by replacing the rate upper bound with the exact mutual information T(X; Y)):
rale
F())
R(D)
F(A) AD
F(A) = inf (() XBIp(X,
(4)
distortion
Qy
-o kg
As illustrated by Fig. F(A) is the maximum R-axis intercept of a straight line with slope -X, among all such lines that lie below langent
The geometry of the R-D lower bound problem. For a given slope -λ, we seek to maximize the rate-axis in- tercept, E[- log g()), over all g 2 functions admissible according to Eq. 6.
Figure 1:
R(D); the R- D curve can then be found by taking the upper enve- lope of lines with slope -X and R-intercept F FA)) i.e., R(D) max F(A) AD. This key result is captured mathematically by Lemma A. and the following theorem:
Theorem 4.1. (Csiszár, 1974b) (As follows, all the expectations with respect to the data source Lv. X Px.) Under basic conditions (e.g., satisfied by bounded p; see Appendix A.2), it holds that
F(X)
max [ maa{E] g(z)
log (())),,
(5)
where the maximization is over all non-negative functions x + [0, oo) satisfying the constraint
exp(- (pp(-AppXX,1)) v))
exp( (1-xppx,)))) v)) dPx I 1, Vy y
E
9(X)
g(x)
In other words, every admissible g yields a lower bound of R(D), via an underestimator of the intercept E|- log ((X) F()). We give the origin of this result in related work in Section 5,
Proposed Unconstrained Formulation. The constraint in Eq. 6 is concerning-in is a family of possibly infinitely many constraints, one for each y. To make the problem easier work with, we propose to eliminate the constraints by the following transformation. Let g be defined in terms of another function u(x) 2 0 and scalar e depending on u, such that
eepp - Ap(XX
g(x) = (z)) where c:= sup (v),), and V (.. = E yEy
x(X)
This reparameterization of g is without loss of generality, and can be shown to always satisfy the constraint in Eq. 6. While this form of y bears a superficial resemblance energy-based model
4

---

Published as a conference paper at ICLR 2022
(LeCun al., 2006), with 1 resembling a normalizing constant, there is important difference: c= sup, (3) is in fact the supremum of a family of "partition functions" (.(y)) indexed by y, we thus refer to as the nup-partition function. Although all these quantities have X-dependence, omit this from our notation since À is a fixed input parameter (as in the upper bound algorithm).
Consequently, F is now the result of unconstrained maximization over all functions, and we obtain a lower bound on it by restricting u to a subset of functions with parameters @ (e.g., neural networks),
F(A)
max El log uf (X) log sup V (1)) 2 mg max log u (X) log sup ()(y)} 4
yEy
yey
Define the 0-parameterized objective {()) == E[- log (X)] - c(0), with c(0) = sup V Voly). Given samples of X, we can in principle maximize {(O) by (stochastic) gradient ascent. However, computing the sup-partition function c()) poses serious computation challenges: even evaluating Ve(y) for a single y value involves a potentially high-dimensional integral w.r.t. Px; this is only exacerbated by the need globally optimize w.r.t. y, an NP-hard problem even in one-dimension.
Proposed Method. To tackle this problem, we propose an over-estimator of the sup-partition function inspired by IWAE (Burda et al., 2015). Fix 0 for now; we denote the inteeara in Eq. 7 by (0,,y) y) = exp -Aplaaaa. (so c = suppey Blyb (X, v)), and omit the dependence on to simplify notation. Given k 2 1 ii.d. random variables X define the estimator C2 sup 1 Σ, (X,, 3). We prove in Theorem A.4 that E[C, 2 ... c, i.e. Cg is in expectation over-estimator of the sup-partition function Similarly the Importance- Sreprhance-Weightedd ELBO (Burda et al.. 2015), the bias of this estimator decreases monotonically as + 00, and asymptotically vanishes under regularity assumptions. In light of this, replace by E E(Ca] and obtain a sampll under-estimator of the objective E(O) (which in turn underestimates F(A)):
a(x)
(0) == E - log (X)) log E[Cx]; moreover, & (0) S (2( S... S E()).
In order to apply stochastic gradient ascent, overcome two more technical hurdles. First, each draw of C3 requires solving global maximization problem. We note that by restricting to a squared-error p and y/ = X, Cs be computed by finding the mode of Gaussian mixture density; for this we use the method of Carreira-Perpinan (2000), essentially by hill-climbing from each of the & centroids. Second, turn- log [C] into an expectation, we follow Poole al. (2019) and underestimate it by linearizing - log around a scalar parameter a > 0, resulting in the following lower bound objective:
(0) = EE[ log u X)] EE(Ck]/ a
log 1.
(8)
(0) can finally be estimated by sample averages, and yields a lower bound on the optimal intercept F(A), via (f) I (0) S E(O) < F()). A trained model ug- then yields an R-D lower bound, (D) = -D + (0*). We give more detailed derivation and pseudocode in Appendix A.4.
5 RELATED WORK
The past few years have seen significant progress in applying machine learning lossy data compression. Theis et al. (2017); Ballé et al. (2017) first showed that a particular type of 3-VAE can be trained to perform data compression using the same objective as Eq. 3. The variational distributions in such a model have shape restrictions to simulate quantization and entropy coding (Ballé et al., 2017). Our upper bound is directly inspired by this line of work, and suggests that such a model can in principle compute the source R-D function when equipped with sufficiently expressive variational distributions and "rich enough" decoder (see Sec. 3). We note however not all compressive autoencoders admit a probabilistic formulation (Theis et al., 2017); recent work has found training with hard quantization to improve compression performance (Minnen & Singh, 2020), and methods have been developed (Agustsson & Theis, 2020; Yang et al., 2020b) to reduce the gap between approximate quantization at training time and hard quantization at test time. Departing from compressive autoencoders, Yang et al. (2020c) and Flamich et al. (2020) use Gaussian B VAES for data compression and exploit the flexibility of variable-width Gaussian posteriors. Flamich al. (2020)'s 's method, and more generally, reverse channel coding (Theis& & Yosri, 2021), can transmit a sample of Qzix with a rate close to that optimized by our upper bound model in Eq. 3, more precisely, T(X; Z) log(T(X: Z) + 1) + O(1). i.e., in this ohe-shot setting (which is standard for
Machine Learning:
5

---

Published as a conference paper at ICLR 2022
neural image compression), R(D) is no longer achievable; rather, the achievable R-D performance is characterized by R(D) + loog(o(())))) + +)+ (11). Therefore, our R-D bounds can be shifted upwards by this logarithmic factor to give an estimate of the achievable R-D performance in this setting.
Information theory has also broadly influenced unsupervised learning and representation learning. The Information Bottleneck method (Tishby et al., 2000) directly motivated by, and extends R-D theory and the BA algorithm. Alemi et al. (2018) analyzed the relation between generative modeling and representation learning with a similar R-D Lagrangian to Eq. 2, but used abstract, model-dependent distortion p(y, x) := log (3]l) with an arbitrary y and without considering a data compression task. Recently, Huang et al. (2020) proposed to evaluate decoder-based generative models by computing restricted version of Rw (with Qy fixed); our result in Sec. 3 (R,, 2 R(D)) gives a principled way to interpret and compare these model-dependent R-D curves.
Information Theory: While the BA algorithm (Blahur, 1972; Arimoto, 1972) computes the R(D)
of discrete source with known distribution, no tool currently exists for the general and unknown case. Riegler et al. (2018) share our goal of computing R RDD) of a general source, but still require the source to be known analytically and supported on known reference measure. Harrison & Kontoyiannis (2008) consider the same setup as ours of estimating R(D) of an unknown source from samples, but focus on purely theoretical aspects, assuming prefect optimization. They prove statistical consistency of such estimators for a general class of alphabets and distortion metrics, assuring that our stochastic bounds on R(D) optimized from data samples, when given unlimited computation and samples, can converge to the true R DDD.. Perhaps closest in spirit to our work is by Gibson (2017), who estimates lower bounds on R RDD) of speech and video using Gaussian autoregressive models of the source. However, the correctness of the resulting bounds depends on the modeling assumptions.
variational lower bound on R(D) was already proposed by Shannon (1959), and later extended (Berger, 1971) to the present version similar to Theorems 4.1 and A.2. In the finite-alphabet case, the maximization characterization of R(D) follows from taking the Lagrange dual of its standard definition in Eq. I; the dual problem can then be solved by convex optimization (Chiang & Boyd, 2004), but faces the same computational difficulties as the BA algorithm. Csiszár (1974b) proved the general result in Theorem 4. applicable to abstract alphabets (in particular, with source X taking values in an arbitrary probability space), by analyzing the fixed-point conditions of the BA algorithm.
6 EXPERIMENTS
We estimate the R-D functions of variety of artificial and real-world data sources, in settings where the BA algorithm is infeasible and no prior known method has been applied. On Gaussian sources, our upper bound algorithm is shown to converge to the exact R-D function, while our lower bounds become increasingly loose in higher dimensions, an issue we investigate subsequently. We obtain tighter sandwich bounds on particle physics and speech data than on similar dimensional Gaussians, and compare with the performance of neural compression. We further investigate looseness in the lower bound, experimentally establishing the intrinsic dimension of the data as a much more critical contributing factor than the nominal/ambient dimension. Indeed, we obtain tight sandwich bounds on high- nigh-dimension G-N-geneaaeed images with a low intrinsic dimension, and compare with popular neural image compression methods. Finally, we estimate bounds on the R-D function of natural images. The intrinsic dimension is likely too high for our lower bound to be useful, while our upper bounds on the Kodak and Tecnick datasets imply at least one dB (in PSNR) of theoretical room for improving state-of-the-aa image compression methods, at various bitrates. We also validate our bounds against the BA algorithm over the various data sources, using 2D marginal of the source to make BA feasible. We provide experimental details and additional results in Appendix A.5 and A.6.
6.1 GAUSSIAN SOURCES
We start by applying our algorithms to the factorized Gaussian distribution, one of the few sources with an analytical R-D function. We randomly generate the Gaussian sources in increasing dimensions.
For the upper bound algorithm, we let Qy and Qy x be factorized Gaussians with learned parameters, predicting the parameters of Qyyx by a 1-layer MLP encoder. As shown in Fig. 2a-top, on n = 1000 dimensional Gaussian (the results are similar across all the we tested), our upper bound (yellow

---

Published as a conference paper at ICLR 2022
RUG
Poin (rats per sample dimession)
ReTR) ree Diee
Rate (rats per pampan pameee
Rate (eats pe Smmp
tre HCM
Bellé 2001
propoed Ae1
propoved (D)
Elauu-AAiiote R1A)
RaDi w
propoee (a
propore (i.!))
RatP)
0.25
0.58
LaMiE
13 Distortion (mean squared arros)
55
(d) visualizations for
Datortion (meen aquared weroo)
Distortion (mean squared arror)
physies and speech data
(a) Gaussian
(b) particle physics
(c) speech
Figure 2: 2a top: R-D upper bound estimates a randomly generated = 1110000imeesiiii Gaus- sian source; bottom: R-D lower bound estimates on standard Gaussians with increasing dimensions (the result of the BA algorithm for 2 is also shown for reference). 2b top: estimated R-D bounds and the R-D performance of Ballé et al. (2021) ) on the particle physics dataset; bottom: the same experiment but on a 2D marginal distribution of the data, to compare with the BA algorithm. 2c: the same set of experiments as in 2b, repeated on the speech dataset. 2d top: histogram of the 2D marginal distribution of the physics data; bottom: example spectrogram computed on speech clip.
curve) accurately recovers the analytical R(D). We also optimized the variational distributions in latent space Z with varying dimensions, using an MLP decoder to map from 2 to y (see Sec. 3). The resulting bounds are similarly tight when the latent dimension matches or exceeds the data dimension n (green, brown), but become loose otherwise (red and purple curves), demonstrating the importance of a rich enough latent space for a tight R-D bound, as suggested by our Theorem A.3.
For the lower bound algorithm, we parameterize log u by a 2-layer MIP,, and study the effect of source dimension n and the number of samples & used in our estimator Cg (and objective (.). To simplify comparison of results across different source dimensions, here we consider standard Gaussian sources, whose R-D curve does not vary with v we scale the rate by 1 the results randomly generated Gaussians shows that the resulting bounds quickly become loose with increasing source dimension. This is due to the bias of our estimator Cj for the sup-partition function, which causes under-estimation in the objective l.- While C is defined similarly to M-estimator (Van der Vaart, 2000), analyzing its convergence behavior is not straightforward, as it depends on the function u being learned. In this experiment, we observe the bias of C& is amplified by an increase large l is required for effective training. In another experiment, we estimate that the & needed to close the gap in the lower bound increases exponentially in n; see results Fig. 4, and detailed discussion on this, in Appendix A.5.2. Fortunately, as we see in Sec. 6.3, the bias in our lower bound appears to depend on the intrinsic rather than (often much higher) nominal dimension data, giving us a more favorable trade-off between computation and a tighter lower bound as controlled by ht
n (i.e., rate per sample per dimension); similar. First, fix k 1024; Fig. 2a-bottom
X, such that an increasingly
6.2 DATA FROM PARTICLE PHYSICS AND SPEECH
The quickly deteriorating lower bound on higher (even 16) dimensional Gaussians may seem dis- heartening. However, the Gaussian is also the hardest continuous source to compress under squared error (Gibson, 2017), and real-world data often exhibits considerably more structure than Gaussian noise. In this subsection, we experiment on data from partiele physics (Howard et al., 2021) and speech (Jackson et al., 2018), and indeed obtain improved sandwich bounds compared to Gaussians.
First, we consider the Z-boson decay dataset from Howard et al. (2021), containing n==1--dieesiooal vectors of four-momenta information from independent particle decay events. We ran the neural compression method from (Ballé et al., 2021), 1). as well as our bounding algorithms with similar configurations before, except we fix k 2048 for the lower bound, and use normalizing flow for Qz in the upper bound model for better expressiveness. Fig. 2b top shows the resulting estimated R-D bounds (sandwiched region colored in red) and the operational R-D curve for Ballé al. (2021)
7

---

Published as a conference paper at ICLR 2022
2
Kodak
Mixcc 2018 (hpeaprke)
Mirees 2000
00
16
proyood Rx(D)
(2)
10
Rate (eats pe Semain
Quality (PSNR)
35
-4
prppoeed Rp( (Bedie-VAA)
proposed for(D) (Minnen 2020 (AVAFF
82
VTW
AV11
Mirns 2018 -mppeprr-
28
Minnen 2018 (hppeppior)
BPG
10
JFFC2000
26
1,000 0L00ES
2,0000
EUMG
0L2000
Rate (BPP, bits per sample per pixel)
Distortion
(meon squared error)
Figure 3: Left: 128x 128 GAN-generated images, with intrinsic dimension d = 4. Middle: Bounds on R(D) of GAN images, for d = 2 (top) and d = 4 (bottom). Right: Quality-rate curves of ours and state-of-the-an image compression methods on Kodak (1993), corresponding to R-D upper bounds.
(blue). The resulting sandwich bounds appear significantly tighter here than on the Gaussian source with equal dimension (n = 16, bottom curve in Fig. 2a bottom), and the neural compression method operates with a relatively small average overhead of 0.5 nat/sample relative to our upper bound. To also compare to the ground-truth R(D) as estimated by the BA algorithm, we created a 2-dimensional marginal data distribution (plotted in Fig. 2d top), so that the BA algorithm can be feasibly run with a fine discretization grid. As shown in Fig. 2b, the BA estimate of R (D) (green) almost overlaps with our upper bound on the 2D marginal, and is tightly sandwiched from below by our lower bound.
We then repeat the same experiments on speech data from the Free Spoken Digit Dataset (Jackson et al., 201 8). We constructed our dataset by pee-prccessing the audio recordings into spectrograms (see, e-g., Fig. 2d bottom), then treating the resulting n = 33--iimesionna frequency feature vectors as independent across time. As shown in Fig. 2c, the gap in our R-D bounds appears wider than on the physics dataset (top), but the results on the corresponding 2D marginal appear similar (bottom).
6.3 THE EFFECT OF INTRINSIC V.S. NOMINAL DIMENSION OF DATA
It is known that learning manifold has a complexity that depends on the intrinsic dimension of the manifold, but not on its ambient dimension (Narayanan & Mitter, 2010). Our experiments suggest a similar phenomenon for our lower bound, and show that we can still obtain tight sandwich bounds on data with a sufficiently low intrinsic dimension despite the high ambient dimension.
First, explore the effect of increasing the ambient dimension, while keeping the intrinsic dimension of the data fixed. We borrow the 2D banana- sanana-source from Ballé et al. (2021), and randomly embed sandwich bounds on
in R As shown in Fig. 7 and Fig. (in Appendix due space constraint), the 2D source appear tight, and closely agree with BA (similar to the results in Fig. 2); moreover, the tightness appears unaffected by the increase in ambient dimension this for up to 1000). Unlike in the Gaussian experiment, where increasing required seemingly exponentially larger k for a good lower bound, here constant k: 1024 worked well for all n.
4, 16, and 100 (we verified
Next, we experiment on high-dimension GAN-generated images with varying intrinsic dimension, and obtain R-D sandwich bounds that help assess neural image compression methods. Following Pope et al. (2021), we generate 128 x 128 images of basenji from pre-trained GAN, and control the intrinsic dimension by zeroing out all except d dimensions of the noise input to the GAN. As shown in Fig. 3-Left, the images appear highly realistic, showing dogs with subtly different body features and in various poses. We implemented 6-layer ResNet-VAL VAE (Kingma et al., 2016) for our upper bound model, and simple convolutional network for our lower bound model. Fig. 3-Middle plots our resulting R-D bounds and the sandwiched region (red), along with the operational R-D curves (blue, green) of neural image compression methods (Minnen et al., 2018; Minnen & Singh, 2020) trained on the GAN images, for d = 2 and = 4. We see that despite the high dimensionality estimate R(D) to
(n 128 1 128 3), the images require few nats to compress; e.g., for d
8

---

Published as a conference paper at ICLR 2022
be between 4 and 8 nats per sample at D 0.001 (30 dB in PSNR). Notably, the R-D curve of Minnen & Singh (2020) stays roughly within a factor of 2 above our estimated true R(D) region. The neural compression methods show improved performance as d decreases from 4 to 2, following the same trend as our R-D bounds. This demonstrates the effectiveness of learned compression methods at adapting to low-dimensional structure in data, in contrast traditional methods such as JPEG, whose R-D curve on this data does not appear to vary with d and lies orders of magnitude higher.
6.4 NATURAL IMAGES
To establish upper bounds on the R-D function of natural images, we define variational distributions (Qz, ^z)x) x on a Euclidean latent space for simplicity, and parameterize them as well as a learned decoder w via hierarchical VAEs. We consider two VAE architectures: 1. we borrow the convolutional autoencoder architecture of state-of-the-art image compression model (Minnen Singh, 2020), but use factorized Gaussians for the variational distributions (we still keep the deep factorized hyperprior, but no longer convolve it with a uniform prior); 2. we also reuse our ResNet- VAE architecture with stochastic layers from the GAN experiments (Sec. 6.3). We trained the models with mean-squared error (MSE) distortion and various on the COCO 2017 (Lin et al., 2014) images, and evaluated them on the Kodak (1 (1993) and Tecnick (Asuni & Giachetti, 2014) datasets. Following image compression conventions, we report the rate in bits-per-pixel, and the quality (i.e., negative distortion) in PSNR averaged over the images for each (A, model) pair 2, The resulting quality-rate (Q-R) curves can be interpreted as giving upper bounds on the R-D functions of the image-generating distributions. We plot them in Fig. 3, along with the Q-R performance (in actual bitrate) of various traditional and learned image compression methods (Ballé et al., 2017; Minnen et al., 2018; Minnen & Singh, 2020), for the Kodak dataset (see similar results on Tecnick in Appendix Fig. 11). Our 3 3-VAE version of (Minnen & Singh, 2020) (orange) lies on average 0.7 dB higher than the Q-R curves of the original compression model (red) and VTM (green). With a deeper latent hierarchy, our ResNet- (essNee-(--VAA gives higher Q-R curve (blue) that is on average 1.1 dB above the ofate-of-the-art Q-R curves (gap shaded in cyan). We leave it to future work to investigate which choice of autoencoder architecture and variational distributions are most effective, as well as how the theoretical R-D performance of such a 3- VAE can be realized by a practical compression algorithm (see discussions in Sec. 5).
7 DISCUSSIONS
In this work, we proposed machine learning techniques to computationally bound the rate-distortion function of a data source, a key quantity that characterizes the fundamental performance limit of all lossy compression algorithms, but is largely unknown. Departing from prior work in the information theory community (Gibson, 2017; Riegler et al., 2018), our approach applies broadly to general data sources and requires only i.i.d. samples, making it more suitable for real-world application.
Our upper bound method is a gradient descent version of the classic Arimmoo algorithm, and closely relates (and extends) variational autoencoders from neural lossy compression research. Our lower bound method optimizes a dual characterization of the R-D function, which has been known for some time but seen little application outside of theoretical work. Due to difficulties involving global optimization, our lower bound currently requires a squared eттor distortion for tractability in the continuous case, and only tight on data sources with a low intrinsic dimension. We hope that a better understanding of the lower bound problem will lead to improved algorithms in the future.
To properly interpret bounds on the R-D function, we emphasize that the significance of the R-D function is two-fold: 1. for a given distortion tolerance D. no coding procedure can operate with a rate less than R(D), and that 2. this rate is asymptotically achievable by potentially expensive block code. Thus, while a lower bound makes universal statement about what performance is "too good to be true", , the story is more subtle for the upper bound. The achievability proof relies on a random coding procedure that jointly compresses multiple data samples in arbitrarily long blocks (Shannon, 1959). When compressing at a finite block length b, R(D) is generally no longer achievable due to O ) rate overhead (Kontoyiannis, 2000; Kostina & Verdú, 2012). Extending our work to the settings of finite block lengths or non-memoryless sources may be additional future directions.
Technically, to estimate an R-D upper bound with MSF as p. one should compute the distortion by averaging MSF images; however, the results of many image compression baselines are only available in average PSNR.

---

Published as a conference paper at ICLR 2022
ACKNOWLEDGEMENTS
Yibo Yang acknowledges support from the Hasso Plattner Institute UCL. This material is in part based upon work supported by DARPA under Contract No. HR001 12000211 Stephan Mandt acknowledges support by the National Science Foundation (NSF) under the NSF CAREER Award 2047418: NSF Grants 1928718, 2003237 and 2007719; the Department of Energy under grant DESCO022331, as well as Intel, Disney, and Qualcomm. Any opinions, findings and conclusions or recommendations expressed in this material those of the author(s) and do not necessarily reflect the views of DARPA or NSF.
REFERENCES
E. Agustsson and L. Theis. Universally Quantized Neural Compression. In Advances in Neural Information Processing Systems 33, 2020.
Eirikur Agustsson, David Minnen, Nick Johnston, Johannes Balle, Sung Jin Hwang, and George Toderici. Scale-space flow for end-to-end optimized video compression. Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 3503-8512 2020.
Alexander Alemi, Ben Poole, Ian Fischer, Joshua Dillon, Rif A Saurous, and Kevin Murphy. Fixing broken ELBO. In International Conference on Machine Learning, pp. 159-168. PMLR, 2018.
Suguru Arimoto. An algorithm for computing the capacity of arbitrary discrete memoryless channels. IEEE Transactions on Information Theory, 18(1):14-20, 1972.
N. Asuni and A. Giachetti. 'ESTIMAGES: A large-scale archive for testing visual devices and basic image processing algorithms (SAMPLING 1200 RGB set). In STAG: Smart Tools and Apps for Graphics, 2014. URL httpp::/souureeeeeeeeeee OLD/OLD_ SAMP WNteLINNN//eeeteeeeeeeeeeeeeeeee http://s/...ooeeeee....eeee//eee
zip.
J. Ballé, V. Laparra, and E. P. Simoncelli. End-to-end Optimized Image Compression. In International Conference on Learning Representations, 2017.
Johannes Ballé, David Minnen, Saurabh Singh, Sung Jin Hwang, and Nick Johnston. Variational Image Compression with Scale Hyperprior. ICLR, 2018
J. Ballé, P. A. Chou, D. Minnen, S. Singh, N. Johnston, E. Agustsson, S. J. Hwang, and G. Toderici. Nonlinear transform coding. JEEE Trans. on Special Topics in Signal Processing, 15, 2021.
T Berger. Rate distortion theory, mathematical basis for data compression (prentice-hall, Inc. Englewood Cliffs, New Jersey, 1971.
Toby Berger and Jerry D Gibson. Lossy source coding. IEEE Transactions on Information Theory, 2466):269,,,2 1998.
R. Blahut. Computation of channel capacity and rate-distortion functions. IEEE Transactions on Information Theory, 18(4):4604477 1972. doi: 10.1109//T..1972115553
Andrew Brock, Jeff Donahue, and Karen Simonyan. Large scale gan training for high fidelity natural image synthesis. arXiv preprint arXiv: 1096, 2019.
Yuri Burda, Roger Grosse, and Ruslan Salakhutdinov, Importance weighted autoencoders. arXiv preprint arXiv:1509.005/2 2015.
Miguel A. Parreira-Perpinan. Mode-finding for mixtures of gaussian distributions. IEEE Transactions on Pattern Analysis and Machine Intelligence, (2(11))1318-122 2000.
Miguel A. Carreira-Perpinan, Gaussian mean-shift is an em algorithm. IEEE Transactions on Pattern
Analysis and Machine Intelligence, 2965)7777-776, 2007. doi: 10.1 109/TPAM1.2007.10557
Miguel A. Carreira-Perpina How many modes can a Gaussian mixture have, 2020. URL https: / fuculty.... ucmerced edu aodd/mmarreii racterrrerrrpppeeeee ra-perrinnn//rrreerrrhhh hew..che////......eee html.
10

---

Published as a conference paper at ICLR 2022
Zhengxue Cheng. Heming Sun, Masaru Takeuchi, and Jiro Katto. Learned image compression with discretized gaussian mixture likelihoods and attention modules. arXiv preprint arXiv:2001.01568, 2020.
Mung Chiang and Stephen Boyd. Geometric programming duals of channel capacity and rate distortion. IEEE Transactions on Information Theory, 50(2):245-22) 258,, 2004.
T. M. Cover and J. A. Thomas. Elements of Information Theory, volume 2. John Wiley & Sons, 2006.
I. Csiszár. On the computation of rate-distortion functions (corresp.). IEEE Transactions on Information Theory, 20(1):122-224, 1974a. doi: 10.1109/TIT.197.10065 1055146.
Imre Csiszár. On extremum problem of information theory. Studia Scientiarum Mathematicarum Hungarica, 9, 01 1974b.
Thomas Ferguson. A course in large sample theory. Routledge, 2017.
G. Flamich, M. Havasi, and J. M. Hernández-Lobato. Compressing Images by Encoding Their Latent Representations with Relative Entropy Coding, 2020. Advances in Neural Information Processing Systems 34.
Brendan J. Frey and Geoffrey E. Hinton. Efficient stochastic source coding and an application to bayesian network source model. The Computer Journal, 3)) 1997.
Jerry Gibson. Rate distortion functions and rate distortion function lower bounds for real-world sources. Entropy, 19611)6604, 2017.
Robert M Gray. Entropy and information theory. Springer Science & Business Media, 2011.
Matthew T. Harrison and Ioannis Kontoyiannis. Estimation of the rate distortion function. IEEE
Transactions on Information Theory, 54(3):375733762, Aug 2008. ISSN 0018-9448 doi: 10.1109/ tit.2008.926388 URL http:: / dx ..d.xddoo. doi org/orggg///oooo 110. 0. 1109/1111TT TIT 2008. 926387.
J. Hayes, A. Habibi, and P. Wintz. Rate-distortion function for gaussian source model of images (corresp.). IEEE Transactions on Information Theory, 16(4):507-509 1970. doi: 10. 109/TIT.
Irina Higgins, Loic Matthey, Arka Pal, Christopher Burgess, Xavier Glorot, Matthew Botvinick, Shakir Mohamed, and Alexander Lerchner. beta-vae: Learning basic visual concepts with constrained variational framework. Iclr, 2(5):6, 2017.
Emiel Hoogeboom, Jorn Peters, Rianne van den Berg, and Max Welling. Integer discrete flows and lossless compression. In Advances Neural Information Processing Systems, pp. 12134-12144, 2019.
Jessica N. Howard, Stephan Mandt, Daniel Whiteson, and Yibo Yang. Foundations of fast, data- driven, machine-learned simulator, 2021.
Sicong Huang, Alireza Makhzani, Yanshuai Cao, and Roger Grosse. Evaluating lossy compression rates of deep generative models. International Conference Machine Learning, 2020.
Zohar Jackson, César Souza, Jason Flaks, Yuxin Pan, Hereman Nicolas, and Adhish Thite. http:////eeeesppeeeeeiii--itiiie 5281/ 1342401.
v1.0.8, August 2018. URL https:://bii.rrg/////
Eric Jang, Shixiang Gu, and Ben Poole. Categorical reparameterization with gumbel-softmax. arXiv preprint /61 1.0/ 144, 2016.
D. Kingma and M. Welling. Auto-encoding variational Bayes. In International Conference on Learning Representations, 2014.
Durk P Kingma, Tim Salimans, Rafal Jozefowicz, Xi Chen, Ilya Sutskever, and Max Welling. Improved variational inference with inverse autoregressive flow, Advances in neural information processing systems, pp. 4743 2016.
11

---

Published as a conference paper at ICLR 2022
Ivan Kobyzev, Simon J.D. Prince, and Marcus A. Brubaker. Normalizing flows: An introduction and review of current methods. IEEE Transactions on Pattern Analysis and Machine Intelligence, doi: 10.1 10.1109/pami.220022292222
43(11):3964-3979, Nov 2021. ISSN 1939-3539.
URL
http::/dxx. doi. wtt.ooootoooo..oooooo 1109/TPAMM 2020. 2992934.
Kodak. PhotoCD PCD0992, 1993. URL http: / rōkk us htt///graphhhkokoeee
loannis Kontoyiannis. Pointwise redundancy lossy data compression and universal lossy data compression. IEEE Transactions on Information Theory, 46(1):136-152, 2000.
Victoria Kostina. When is shannon's lower bound tight at finite blocklength? In 2016 54th Annual Allerton Conference on Communication, Control, and Computing (Allerton), pp. 982-989. IEEE, 2016.
Victoria Kostina and Sergio Verdú. Fixed-length lossy compression in the finite blocklength regime. IEEE Transactions on Information Theory, (116:3330--3338 2012.
Yann LeCun, Sumit Chopra, Raia Hadsell, M Ranzato, and F Huang. A tutorial on energy- based learning. Predicting structured data, 100), 2006.
Jasper H. Lee, Jerry Li, Christopher Musco, Jeff M. Phillips, and Wai Ming Tai. Finding the mode of a kernel density estimate, 2019.
Tsung-Y! Yi Lin, Michael Maire, Serge Belongie, James Hays, Pietro Perona, Deva Ramanan, Piotr Dollár, and C Lawrence Zitnick. Microsoft coco: Common objects in context. In European conference on computer vision, pp. 740-755. Springer, 2014.
David JC MacKay. Information theory, inference and learning algorithms. Cambridge university
press, 2003.
Chris Maddison, Andriy Mnih, and Yee Whye Teh. The concrete distribution: A continuous relaxation of discrete random variables. arXiv preprint arXiv:1611.00712, 2016.
Matt Mahoney. Rationale for large text compression benchmark. Retrieved (Oet. 1st, 2021) from:
http://mathhaaooeeeeedddccoooooo htmml 2009.
Paul Milgrom and Iya Segal. Envelope theorems for arbitrary choice sets. Econometrica, 70(2): 583-601, 2002.
D. Minnen and S. Singh. Channel-wise autoregressive entropy models for learned image compression. In IEEE International Conference on Image Processing (ICIP), 2020.
D. Minnen, J. Ballé, and G. D. Toderici. Joint Autoregressive and Hierarchical Priors for Learned Image Compression. In Advances in Neural Information Processing Systems 31. 2018.
Hariharan Narayanan and Sanjoy Mitter. Sample complexity of testing the manifold hypothesis. Advances in neural information processing systems, 23, 2010.
George Papamakarios, Theo Pavlakou, and lain Murray. Masked autoregressive flow for density estimation. In Advances in Neural Information Processing Systems, pp. 2338 2237, 2017.
Polyanskiy and Y Wu. Lecture notes on information theory. 2014.
Ben Poole, Sherjil Ozair, Aaron Van Den Oord, Alex Alemi, and George Tucker. vari- ational bounds of mutual information. In Kamalika Chaudhuri and Ruslan Salakhutdinov (eds.), Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 5171-5180. PMLR, 09-15 Jun 2019. URL html.
Phillip Pope, Chen Zhu, Ahmed Abdelkader, Micah Goldblum, and Tom Goldstein. The intrin- sic dimension of images and its impact on learning. In International Conference Learning Representations, 2021. URL https: / opwnreview net/ forum? wtt.://ooooummmddeeeee -d=ijkixikk 19x219222222.
Erwin Riegler, Günther Koliander, and Helmut Böleskei. Rate-distortion theory for general sets and measures. arXiv preprint arXiv: / 1804.0884000 2018
12

---

Published as a conference paper at ICLR 2022
T. Salimans, A. Karpathy, X. Chen, and D. P. Kingma. PixeICNN++: A pixelenn implementation with discretized logistic mixture likelihood and other modifications. International Conference on Learning Representations, 2017.
C. E. Shannon. Mathematical Theory of Communication. Bell System Technical Journal, 27: 379-423, 1948.
CE Shannon. Coding theorems for discrete source with a fidelity criterion. IRE Nat. Conv. Rec., March 1959, 4:142-163, 1959.
L. Theis and N. Yosri. Algorithms for the communication of samples. preprint, 2021. URL https://arxte httpa////rxivo.../aa.......
org//ogg/ss....ooooo /tt.b////110 2210.11011111 12805.
L. Theis, W. Shi, A. Cunningham, and F. Huszár. Lossy Image Compression with Compressive Autoencoders. In International Conference on Learning Representations, 2017.
Naftali Tishby, Fernando C Pereira, and William Bialek. The information bottleneck method. arXiv preprint physicx/0004077, 2000.
Aad W Van der Vaart. Asymptotic statistics, volume Cambridge university press, 2000.
lan H Witten, Radford M Neal, and John G Cleary. Arithmetic coding for data compression. Communications of the ACM, 30(6):20--540, 1987.
Ruihan Yang. Yibo Yang, Joseph Marino, and Stephan Mandt. Hierarchical autoregressive modeling for neural video compression. In International Conference on Learning Representations, 2020a.
Yibo Yang, Robert Bamler, and Stephan Mandt. Improving inference for neural image compression. In Neural Information Processing Systems (NeurIPS), 2020, 2020b.
Yibo Yang, Robert Bamler, and Stephan Mandt. Variational Bayesian Quantization. International
Conference Machine Learning, 2020c.
Yibo Yang, Stephan Mandt, and Lucas Theis. An introduction neural data compression. arXiv preprint arXiv:2202.06533. 2022.
Jacob Ziv and Abraham Lempel. A universal algorithm for sequential data compression. IEEE Transactions on information theory, 233))337-333, 1977.
13

---

Published as a conference paper at ICLR 2022
A APPENDIX
A.1 TECHNICAL DEFINITIONS AND PREREQUISITES
In this work consider the source be represented by a random variable measurable function on an underlying probability space F,, F. P), and Px is the image measure of under X. We suppose the source and reproduction spaces are standard Borel spaces, Axx and (V' Ay), equipped with sigma-algebras Ax and Ay, respectively. Below we use the definitions of standard quantities from Polyanskiy & Wu (2014).
X. i.c.
The notation Qrjx denotes arbitrary conditional distribution (also known as a Markov kernel), i.., satisfies
Conditional distribution
1. For any X, Qylx=a()) ()) is probability measure
2. For any measurable set B E Ay, x Qrxx==(B) a measurable function on X.
Induced joint and marginal measures Given source distribution Px, each test channel distri- bution Qy defines a joint distribution PxQrıx on the product space y (equipped with the usual product sigma algebra, Ax x Ay) as follows:
PxQrix(E) [ Px((x) / 1((x,) E E)Q
for all measurable sets F. E Ax Ay. The induced y-marginal distribution Py is then defined by
Pr(B) f Qy)x=
for all measurable sets VB E Ay.
KL Divergence We use the general definition Kullback-Leiler (KL) divergence between two probability measures P, Q defined on a common measurable space:
aq log dP,
de
i
KL(PIIQ)
o0,
otherwise.
P P&Q Q denotes that P is absolutely continuous w.r.t. (i.e. for all measurable F,, Q(F) 0 P(E) 0). dP denotes the Radon-Nikodym derivative of P w.r.t. Q: for discrete aQ
distributions, we can simply take be the ratio of probability mass functions; and for continuous distributions, we can simply take to be the ratio of probability density functions.
Mutual Information Given Px and Qyix. the mutual information 1(X;)) defined as
KL(P PI Pv) PxmPy [I=(Qy)x==IIIx)]
where Py is the --magiinal of the joint Px Qxxx Px 2 Py denotes the usual product measure, and (-(--ll) the KL divergence.
For the mutual information upper bound, it's easy show that
To (Qxix Qy) = Exy [KT(QVI=== Qr)] I(X;))++++++))))) =<4 so the bound is tight when Py = Qr.
Obtaining R(D) through the Lagrangian. For each À 2 O, we define the Lagrangian by incor- porating the distortion constraint in the definition of R(D) through a linear penalty:
[(Qy = ((X;Y) xE (p(X,Y)),
(10)
and define its infimum w.r.t. Qyi x by the function
F()) = inf T(X;Y) + XE(p(X,
(11)
Qyix
14

---

Published as a conference paper at ICLR 2022
Geometrically, F(A) is the maximum of the R-axis intercepts of straight lines of slope X, such that they have no point above the R(D) curve (Csiszár, 1974b).
Define Dmn D' RR(D)) so). Since R(D) is convex, for each D > D min there exists 2 such that the line of slope through (D, R(D)) is tangent to the R(D) curve, ie..
R(D') AD' 2 R(D) X(D) F()), AD'.
When this occurs, we say that is associated to D.
Consequently, the RD) curve is then the envelope of lines with slope -A and R-axis intercept F(A). Formally, this can be stated as:
Lemma A.1. (Lemma 1.2. Csiszár (1974b): Lemma 9.7, Gray (2011)). For every distortion tolerance D >> min where D Dinn = inf ( R(D') < oo), it holds that
R(D) max F(A) XD X
(12)
We can draw the following conclusions:
1. For each >> Dmins the maximum above is attained iff associated to D.
2. For a fixed >, if f/Qxxx achieves the minimum of L(,, A), then associated to the point
(p(QVix) xx ie., there exists line with slope -À that is tangent to the R curve (p(O' )) I()))))) where we defined the shorthand p(Qrix) EEp((X Y)]
and InQv(x) = (X; as induced by PxQx1x.
A.2 FULL VERSION OF THEOREM 4.1
Theorem A.2. (Theorem 2.3, Csiszár (1974b); Theorem Kostina (2016).) Suppose that the following basic assumptions are satisfied.
1. R(D) finite for some D. i.e. Dmin :m R(D) oo) < oo;
2. The distortion metric p such that there exists a finite set C y such that
E[min ((X,)) y) < 00 WEE
Then, for each D Dains holds that
R(D) max g(z),)
[
- AD}
(13)
where the maximization is over g(x) 2 and 2 satisfying the constraint
(xxp(-Ap(X,pp) 8))
exp(-Ap(e: v))
E
=
dxx(x) <LYy
(14)
g(X)
g(x)
Note: the basic assumption is trivially satisfied when the distortion is bounded from above; the maximization over g(x) 2 0 can be restricted only g(x) 2 0. Unless stated otherwise, we use log base in this work, so the R(D) above is in terms of nats (per sample).
Theorem 2 can be seen as consequence of Theorem 4. in conjunction with Lemma A.I. We implemented an early version of our lower bound algorithm based on Theorem A.2, generating each R-D lower bound by fixing a target value and optimizing over both À and g as in Equation 13. However, the algorithm often diverged due to drastically changing A. We therefore based our current algorithm on Theorem 1. producing R-D lower bounds by fixing and only optimizing over g (or in our formulation).
A.3 THEORETICAL RESULTS
Theorem A.3, (A decoder-induced R-D function bounds the source R-D function from above). Let X Px be memoryless source subject to lossy compression under distortion p. Let Z be any
15

---

Published as a conference paper at ICLR 2022
measurable space ("latent space" in VAE), and w Z ++ y any measurable function ( ("decoder" in a VAE). This induces a new lossy compression problem with being the reproduction alphabet, under distortion function Pw X x Z + 00), (z) 2) p(xxw(z)). Define the corresponding w-dependent rate-distortion function
(D) =
inf 12x::1p
[)x =
inf x Blx(X
I(X; Z).
Then for any D 2 0, (D) 2 R(D). Moreover, the inequality is tight if is bijective (so that there "no information loss")
Proof. Fix D. Take any admissible conditional distribution Q2ix that satisfies Elp(X, D in the definition of R (D). Define a new kernel Qr between X' and y by Qyxx=r := Qz|x= w- Ve E X, i.., Qyjx= the image measure Qz)x== induced byw. w. Applying data processing inequality the Markov chain X Qzix, Z R Y, we have [XX: 2 (()
Moreover, since Qrjx is admissible in the definition of R(D), i.c.
(oeevixl((XXX)) ExxQ2 x 1 (X, w [2)] < D
therefore have
T(X: 2 1(X;Y) R(D) =
inf T(X
Finally, since T(X; 2) 2 R(D) holds for any admissible Qzix, taking infimum over such Qzix gives R.,. (D) 2 R(D).
To prove (D) = R(D) if w is bijective, it suffices to show that R(D) > (D) We use the same argument as before. Take any admissible Qyix in the definition of R(D). We can
then construct a Qzix by the process X Qvix Y Z. Then by DPI have I(X 2 I(X; Z). Morever, Qz x is admissible: ExxQzIx [p(X, (())) Exx (n(, (Y)) E [o(X, Y)] D. So T(X;Y) 2 T(X;%) 2 (D). Taking infimum such Qrix concludes the proof.
Remark. Although (w)() = Tdentity(o) for any injective w, our construction of Qzix in the last argument requires the inverse of w to exist, so that additionally w(w- Y. Several learned image compression methods have advocated for the use of sub-pixel convolution, i.e. convolution followed by (invertible) reshaping of the results, over upsampled convolution in the decoder, in order to produce better reconstructions (Theis et al., 2017: Cheng et al., 2020). This can be seen making the decoder more bijective, therefore reducing the gap of (D) over R(D), in light of our above theorem.
Corollary A.3.1. (A suitable 3-VAE defines an upper bound on the source R-D function). Let X Px be a data source, p x y + [0, oo) a distortion function, and be any latent space. Consider any 3-VAE consisting of prior distribution Q2 Z, and encoder which specifies conditional distribution Q for each z. Suppose further that the likelihood (observation) model chosen to have density p(a|2) x exp(-p(x,, w(2) for some decoder function w 2 y (a common example being isotropic Gaussian p(z|2) with mean w(z) and fixed variance, specified by a squared error distortion).
Then the two terms of the negative ELBO specifically, the "KL term"
1 [1T(22)===1122)
R = Ex
and the "log- likeliiol term" up to a constant,
D = ==EEx (p(X, w(2)] = E Exx l- + const,
define point (D, R) that lies above the source R(D)-curve; i.e. R R(D).
16

---

Published as a conference paper at ICLR 2022
Proof. R2 I2X (X; 2
inf 2221x Pxxxx (x. Px (X; =: R (D) 2 R(D)
The first inequality is the variational upper bound on mutual information (Eq. 9), the second inequality follows from the definition of the w-induced R-D function (and the fact that EpxQ2 (X,, 2)] =; D), and the last inequality is Theorem A.3.
Remark. As also remarked by Theis et al. (2017), not all distortion functions lead to the 3-VAE interpretation above. For this to work, the function expi-p(z, (()) must be normalizable in æ (w.r.t. suitable base measure on X) for every 2, and the normalizing constant may not depend on
Theorem A.4. (Basic properties of the proposed estimator C₁ of the sup-partition function.) Let X1, X2,... ~ Px be sequence of i..dd. random variables. Let ab : X x y - >R+ be a measurable function. For each k. define the random variable Ck sup, k E, (Xx,, v). Then
1. Ck is an overestimator of the sup-partition function c, ie., (X,,, y)] 2 sup, [xyb(X, v)] C; [[[Ca] = Exx....,, x, ( Te E.
2. The bias Cr decreases with increasing k. i.e., [[C1] 2 E[C2] 2..]] 2 E[Ck+1] 2 sup, (l = c;
3. If l/bb,, y) is bounded and continuous in and if is compact, then C2 is strongly consistent, i.e., Ck converges to e almost surely (as well as in probability, i.e, limg- PP(CC. - cl >
c) 0, Ve > 0), and limg- E[Ck] c.
Proof. We prove each in turn:
1. [[Ca] = Elsup, 1 E 26 (X,, w)1 sup, EI E, (X,, w)1 sup, E|b 8)]
2. First, note that [[C]] 2 [[Cb. since
[[C1] = E[sup 2 [1, 3) = E k
Σ
sup w)] 2 Elsup k
(x, 5) = E[C]] Σ
We therefore have
x+1
1
E[Ck+1] = Elsup
4 (X,,
k
im1
1
E|supf
v/ (X v) + ++1
2 (+(x+1.))) v)
k
+
k
1
S B|suut +
vxxi, y)} + supt k +
2 ((x+1.)))
1=1
E[Ch]
+
B[G1]
1+1 $ [[C]]
k
3. Define the shorthand Mx(y) 1 Σ d v). (vy) = E(((X,))) Denote the sup norm on the set of continuous bounded functions on y by [l.( == sup lf()). Then it holds that
sup My) sup (v)) sup [M - sup (1))
[1] (-)lol
[| (V(-)l
sup (M - (2)),
17

---

Published as a conference paper at ICLR 2022
where made use of the fact that vb 2 O in the first equality, and the reverse triangle inequality in the second-to-last step. By the uniform strong law large numbers (Ferguson, 2017),
lim sup Mv(v) (vy)l = k-000 y
almost surely. Therefore, by the inequality we just showed, we also have
lim | sup Mx(x) - sup (r)) - lim sup [M((y) ((vi] = 10000
k-
almost surely, i.e.,
lim Ck lim sup Mx(y) sup (X, == k=20
k++00
almost surely. Then trivially follows that Cg also converges to in probability, and that c by the bounded convergence theorem.
limg-+o E[Ca]
A.4 DETAILED LOWER BOUND METHOD
Approximations for stochastic gradient ascent. Recall the lower bound objective obtained by replacing the sup-partition function by the mean of the k-sample underestimator C₃:
Eg (0) == E[ log us (X)] log E[Ca]
Unfortunately, we still cannot apply stochastic gradient ascent to &g, two more difficulties remain, which we solve by further approximations. First, C₂ is still hard to compute, as it is defined through a global maximization problem. We note that by restricting to a symmetrie p and y = X' the maximization objective of Cg has the form of reweighted kernel density estimate,
(2,,))
(ap(-~p((x,,0000 v)), u(x;)
K
where each data point is reweighted by u(a) For a squared error distortion, this becomes a Gaussian mixture density, v E, () v) Σ, m. exp(-xXpX)) | y)(²), with centroids defined by the samples Ti, ... Ik. and mixture weights ग. (ku((x))) )1 The global mode Gaussian mixture can generally be found by hill-climbing from each of the & centroids, except in rare artificial examples (Carreira-Pena Perrpiaa,e,r,,nn, 2000; 2020); we therefore use this procedure compute Ck. but note that other methods exist (Lee et al., 2019; Carreira-Perpinan Perpinan, 2007).
The second difficulty is that even we could estimate [[C,] (with samples of C computed by global optimization), the objective f,. requires estimate of its logarithm; naive application of Jensen's 's inequality log [C[ck] E[- log Cx) results in an over-estimator (as does the IWAE estimator), whereas we require a lower bound. Following Poole et al. (2019), we use the following basic lower bound of log(ar) by its linearization around some point a, i.c.,
log(x)) 2 -//a- - log(a) 1, Vz, >
with equality achieved by This results in the final lower bound objective we optimize:
(0) E log un (X)] E[Cx]/a log
(15)
Since the linear underestimator of log B[Ca] is tightest adaptively to an estimate of B[Ck]] (separately from the [[[C]] term we estimate in the objective function).
E[Cx), ]. in algorithm set
The algorithm. We give a pseudo-code implementation of the algorithm in Algorithm 1. The code largely follows Python semantics, and assumes an auto-iiffernntiation package is available, such as Gradient as provided in Tensorflow We use 7h denote the function in the supremem definition of Ck i.e., 7)())= E exp( - Ap(are, v)} / (z.)-
18

---

Published as a conference paper at ICLR 2022
Example implementation of the proposed algorithm for estimating rate- distortion lower bound (D).
Algorithm I:
Requires:
O, model ug (e.g., a neural network) parameterized by 0, batch sizes k m, and gradient ascent step size
2 while l not converged do
// Draw 2m samples of Cei use the last m samples to set a and the first m samples to estimate E[Ca] in 4.
for j := 2m do
Draw k data samples, (at,..., 3¹, O.CC == coiputommpettipc Ca (0, A, fII,...,
end
m + Oj
2m
// Compute ob jective (0) and update f by gradient ascent.
with GradientTape() as tape:
E = E 7k 0, t.t, + // Estimate of [[C[]
9
l = - Σ log (x E loga +1
10
gradient = tape-apeeeeeee aeee-grdddent(e... 0)
:= gradient
end
Subroutine (y. 0, X, (r1,...
// Evaluate the global maximization objective at y.
Return 2/2 E. (xp(-Ap(z,,)))/u
15
16
Subroutine compute (0, A, (5i, #x]
// Compate the global opt imum of mul),, assuming squared distortion p(x,8) u1..
17
opt.loss -00
for i := to do // Run gradient ascent from the ith mixture component. y = 2;
20
while y not converged do tape:
with GradientTape()
22
loss := 7 (y, 0, x, {II,..., xx3))
gradient tape. aop..gradient(loss, y)
24
: + gradient
25
end
if loss ops_ then
27
opt == loss
28
y = y
end
end
31
= opt.los
return (k)
For given À > O and model ug - RR, the lower bound algorithm works by (stochastic) gradient ascent on the objective (0) (Eq. 8) w.r.t. the model parameters 0. To compute the gradient, we need an estimate of [E[ca] (ultimately log E(Cx)) as part of the objective, which requires us to draw m samples of Cj. In the pseudo-code, separate set of samples of C₂ are also drawn to set the linear expansion point a. In our actual implementation, we set to an exponential moving average [[[]]] estimated from previous iterations (e.g., replacing line 7 of the pseudo-code by 0.2o 0.8E). so only m samples of C need to be drawn for each iteration of the algorithm. We did not find this approximation significantly affect the training results.
Recall C₁ (0) is defined by the result of maximizing w.r.t y. When computing the gradient with respect to 0, differentiate through the maximization operation by evaluating Cs. mv()")
19

---

Published as a conference paper at ICLR 2022
the forward pass, using the argmax y" found by the subroutine compute.Ck. This is justified by appealing to a standard envelope theorem (Milgrom & Segal, 2002).
By Lemma A. each ug trained with given value of À yields linear under-estimator of R(D).
(D) = AD + (8*).
We obtain the final (D) lower bound by taking the upper envelope of all such lines, i.e.,
Rr. = max (D),
AEA
where A is the set of À values we trained with.
Tips and tricks:
To avoid numerical issues, we always parameterize log and all the operations involving C₁ and are performed in the log domain; e.g., we optimize log " instead of YK (which does not affect the result since log is monotonically increasing), and use logsumexp when taking the average of multiple Cs samples in the log domain.
During training, we only run an approximate version of the global optimization subroutine com- pute_Ck for efficiency, as follows. Instead of hill-climbing from each of the l mixture component centroids, we only do from the highest-valued centroids under the objective Yks for small number t (say 10). This approximation is not used when reporting the final results.
A.5 ADDITIONAL. EXPERIMENTAL DETAILS AND RESULTS
A.5.1 COMPUTATIONAL ASPECTS
Our methods are implemented using the Tensorflow library. Our experiments with learned image compression methods additionally used the wwwssoofooo--cpreessse experiments on images were run on Titan RTX GPUs, while the rest of the experiments were run on Intel(R) Xeon(R) CPUs. We used the Adam optimizer for gradient based optimization in all our experiments, typically setting the learning rate between le - and 5e - 4. Training the 8- BAVS for the upper bounds required from a few thousand gradient steps on the lower-dimension problems (under an hour), to a few million gradient steps on the image compression problem (a couple of days; similar to reported Minnen Singh (2020)). With our approximate mode finding procedure (starting hill climbing only from a small number of centroids, see Sec. A.4), the lower bound models required less than 10000 gradient steps to converge in all our experiments.
library. Our
A.5.2 GAUSSIANS
Data Here the data source is an n-dimensional Gaussian distribution with a diagonal covariance matrix, whose parameters are generated randomly. We sample each dimension of the mean uniformly randomly from 0-0.5, 0.5] and variance from [0, 2]. The ground-trui -ouund-truth R-D function of the source computed analytically by the reverse water-filling algorithm (Cover & Thomas, 2006).
Upper Bound Experiments For the Z = y (no decoder) experiment, we chose Qy and Qyix x to be factorized Gaussians; we let the mean and variance vectors of Qy be trainable parameters, and predict the mean and variance of Qy by one fully connected layer with 2n output units, using softplus activation for the n variance components.
For our experiments involving a decoder, we parameterize the variational distributions Qz and Qzix similarly to before, and an MLP decoder with one hidden layer (with the number of units equal the data dimension n. and leaky ReI.U activation) to map from the 2 to y space. We observe the best performance with a linear (or identity) decoder and simple linear encoder; using deeper networks with nonlinear activation required more training iterations for SGD to converge, and often to poorer upper bound. In fact, for factorized Gaussian source, it can be shown analytically that the optional Qyix= is Gaussian whose mean depends linearly on and an identity (no) decoder is optimal.
20

---

Published as a conference paper at ICLR 2022
14
true F.())
FIN
12
las
R(D) intercept at slope
10
1128
8-
6
11024
2
10
16
20
32
Dimension (11)
Figure 4: R-axis intercept estimates at (negative) slope = 2 from our lower bound algorithm, trained with increasing n and k. The R-axis intercept of optimally tight linear lower bound, (5), increases linearly as a function of n; however, we see that the best achievable R-axis intercept achieved with particular setting of k la plateaus out at log, increases.
Lower Bound Experiments In our lower bound experiments, we parameterize log u by an MLP with 2 hidden layers with 20n hidden units each and SeLU activation, where is the dimension of the Gaussian source. We fixed = 1024 for the results with varying n in Figure 2a-bottom. We vary both & and n in the additional experiment below.
To investigate the necessary k needed to achieve a tight lower bound, and its relation to À or the Gaussian source dimension n, on increasingly high dimensional standard Gaussian sources, each time setting À = - For an n-dimensional standard Gaussian, the true R-intercept, Fz (A). has an analytical form; in particular, (²) 2' In Figure 4, we plot the final objective estimate, Lx, using the converged MLP model, one for each k and n. As we can see, the maximum achievable & plateaus to the value log k as we increase n, and for sufficiently high dimension (e.g., n here), doubling h only brings out constant (log 2) improvement in the objective. This phenomenon relates to the eetimsstintionn bias of Cs when k is too low compared or the "spread" the data, and can be understood as follows.
series of experiments with A ranging from 64 1024
For given sample size k; there is a ""catestrophie" regime produced by increasingly large À (cor- responding to an exceedingly narrow Gaussian kernel), or increasingly high (intrinsic) dimension of the data, that the & data samples appear very far apart. Numerically, this is exhibited by very quick termination of the hill-clibing runs when computing C' (subroutine compute_Ci, well separated, and the mixture centroids are nearly stationary points of the mixture density. The value of the mixture density at each centroid r, be approximated as (x;) = [(. + Ejzı cxp(->p(z;; xr,}// ug (22)) = ku(z, since exp(-Ap(<ı, x;)} x 0 for all 전 i. The maximization problem defining Cj then essentially reduces to checking which mixture component is the highest, and returning the corresponding centroid (or This implies
in Algorithm. 1), since the mixture components
point very close to it), so that C₁ Rs sup, I....h ku(z.)
[[Cc] R E sup iml ku (X,)
sup S Elsup ku(x) x
ku(x))
therefore
44 = E[- log u( (X)] log E(Ck]
(16)
< E| log u( (X) log 1e inf log u(x)
(17)
log k + El- u(X) inf log u(x < logk.
(18)
Thus, when these approximations hold, the maximum achievable lower bound objective & cannot exceed log k. On sources such as high-ding high-dimension (e.g.. n = 10000) Gaussians, or 256 256 patches of natural images, À needs to be the order of 10 higher target even the low-distortion region
21

---

Published as a conference paper at ICLR 2022
of the R-D curve, and the above analysis well describes the behavior of the lower bound algorithm for any value of k we feasibly experimented with.
A.5.3 PARTICLE PHYSICS
Data We borrowed the Z - dataset from Howard et al. (2021), containing 331699 data observations. Each n = 16-dimensional data vector corresponds to independent Z-boson decay event, and is the concatenation of the four-momenta of an electron (e") and positron (e+) both in the theoretical prior space and observed space. The 2D marginal distribution (where we additionally compare with the BA algorithm) is formed by taking the two data coordinates with the least absolute value of correlation. A histogram is given in Figure 5a.
Unlike in the Gaussian experiments, where we generate new random samples from the source on the fly for training and evaluation, here the true underlying source only approximated by a finite number of samples in the given dataset. Since our R-D upper/lower bounds have the form of expectations of learned functions w.r.t. the true underlying source (explained in more detail in Section. A.6), we create a separate held-out test set of 10000 samples for our final estimates of R-D bounds using the trained models. We also use a small subset of the training data as a validation set for model selection, typically using the most expressive model we can afford without overfitting. We use the same train/test split and the same model architectures on the 2D marginal data as on the n = 16 data.
Upper Bound Experiments For our R-D upper bound model, we use a VAE with MLP decoder, a normalizing flow "prior" Qz. and an MLP encoder computing a factorized Gaussian Q21x. We set the latent space to have the same dimensionality n as the data. The encoder and decoder MLPs have three layers with [500, 500, !] units in each, where / = n for the deocder MLP and / = 2n for the encoder MIP (as computes both a mean and variance vector for Q2(x): we use softplus activation in each hidden layer. Qz is parameterized as a base Gaussian distribution transformed by three stacks of MAF (Papamakarios et al., 2017).
For the Nonlinear Transform Coding (NTC) (Ballé et al., 2021) model, we use VAE with similar 500-500-n MLP architecture for the encoder and decoder (thus the latent space has the same dimensionality as the data), but modify the variational distributions to simulate quantization and end-to-end compression. As in (Ballé et al., 2018; Ballé et al., 2021), we let Qz x be a factorized uniform distribution with unit width, and Qz be factorized neural density model convolved with uniform noise.
For the BA algorithm on the 2D marginal, we obtain discretized alphabets and source probabilities as follows. First, we determine the smallest box set S in IR² containing all the data samples (simply the Cartesian product of sample ranges along the and y axes); next, we finely tile up S using small rectangle bins, and compute a histogram of data samples over them. We then let the source and reproduction alphabets be the bin centers, and set the source distribution to be the frequency of samples in each bin. Finally, for computation efficiency. we remove the bin centers from the source alphabet associated with samples; this does not affect the results.
We ran maximum of 2000 steps of the BA algorithm on the test set, terminating early if the objective has not improved by more than 0.0001 in consecutive steps.
Lower Bound Experiments We parameterize our lower bound model log by MLP with 3 hidden layers with 200 hidden units each and SeL.U activation. We fixed k- = 2048.
A.5.4 SPEECH
Data We use the Free Spoken Digit Dataset (Jackson et al., 2018) consisting of recordings of spoken digits from 'zero' to 'nine' at a sampling rate of 8kHz. We use subset of the recordings with a randomly chosen speaker id ('theo"), and keep the original train-test split" consisting of 450 and 50 audio files for the train and test sets, respectively. From each audio file, extracted Short- time Fourier Transform (STFT) coefficients using Tensorflow's 's signal. 63-sample FFT window (frame length-63) in I-sample steps (frame step-1), took the magnitude and converted to log domain, common pre-processing step for computing spectrograms
routine with
hhttps://gi com wwwo///aakoookkk///oooko http://ttt....o.oooooooooooo///e free- www.esssiiiieoeeeeteeeeeee-eee-e /uuaae
22

---

Published as a conference paper at ICLR 2022
0.30
0.000
0.005
0.004
-0.004
0.20
0.003
0.003
0.15
0.002
8
0.002
0.10
0..01
-10
0.001
0.05
0.000
0.000
-0.05
0.00
0.05
12.5
10.0
7.5
5.0
2,5
(a)
(b)
Figure 5: Normalized histograms of the 2D marginals of the particle physics (5a) and speech (5b) datasets. The quantization bins and the associated empirical PMFs are also used define the discretized 2D source distributions which the BA algorithm is
and Mel-frequency cepstral coefficients. The resulting 33-dimensionall Fourier feature vectors are then shuffled across time steps and treated as i.id. samples in our experiments, with over million samples in the training set and 120000 samples in the test set. We follow the same procedure as in the physics experiment, creating 2D marginal from the two coordinates with the least absolute correlation; histogram is given in Figure 5b.
Experiments We the same model architectures and experimental procedureshyppmaraaneess as in the physics experiment.
A.5.5 MANANA-SHAPEL SOURCE
Data The 2D banana aanana-shaped source (Ballé et al., 2021) is RealNVP transform of 2D Gaussian See Figure 6 for plot of its density function.
We embed the banana source in a higher dimension n by multiplying with randomly sampled n x matrix. This is done by simply passing samples of the 2D banana source through linear MLP layer with n output units and no activation, with weight matrix randomly drawn from a Gaussian by the Glorot initialization procedure.
Upper Bound Experiments We reproduced the NTC experiment using the compressive autoencoder as in Ballé et al. (2021), with a 2-dimensional latent space and two-layer MLPs for both the encoder and decoder, with 100 hidden units and softplus activation each hidden layer. For our upper bound model, we use the same MLP architecture as in Ballé et al. (2021), but parameterize Qz by a MAF (Papamakarios al., 2017) and Qzix by factorized Gaussian (doubling the number of output of the encoder), similar on previous experiments.
For the higher-dimension embeddings the banana-shaped source, we use the same 3- 3-VAE as on the 2D source, but set the number of hidden units in each hidden layer to 100n, where n is the dimension of the source. We cap the number of hidden units per layer at 2000 for n > 20, and do not find this to adversely affect the results.
Lower Bound Experiments We MLP with three hidden layers and SeLU activations for the log u model; as in the upper bound models, here we set the number hidden units in each layer to be 100n, and cap it 1000. In all the experiments set k 1024.
SSource code taken from the authors' GitHub repo ht tptsps:s :/ll://////e /ttpttbttttt////////////ttttthtt www/tensorrlow/ fba 9995
htt.e//ssaaooblaa//e//////afafaa sources//ty_scurcces ipynb.
23

---

Published as a conference paper at ICLR 2022
Ball et al. (2021) Rs(P)
Rate (oas pe: sampely
propoued (2
proppoed (2 (n-4)
Blahut Arinato R.D)
Rate: (nots simple)
proppsed Ra(0) (x-4)
propposed Ra(D) (n-30)
propored R2
proppoeed Rg(D) (s-15)
proposed (g(7)) (e0)
proposed
Rc(P) (wm1000)
Distortion (menn squared error)
0.9
11.5
1.0
1.5
Divtortion (mean squarsd smro))
Figure 7: R-D sandwich bounds and the R-D performance of non- linear transform coding (Ballé et al., 2021) the banana- shaped source. For reference, the BA algorithm is also run fine discretization of the source.
Figure 8: R-D sandwich bounds on higher-dimension embed- dings of the banana-shaped source. The tightness of our bounds appear unaffected by the increasing n.
Figure 6: Density plot of the 2D banana-shaped distribution from Ballé et al. (2021).
Figure 9: Random samples of basenji i from BigGAN, d = 2.
Figure 10: Random samples of basenji from a BigGAN, d d= 4.
As illustrated by Figure 8, the tightness of our sandwich bounds do not appear to be affected by the dimension of the data (we verified this up to 1000), unlike on the Gaussian experiment where for a fixed the lower bound became increasingly loose with higher n.
A.5.6 GAN-GENERATED IMAGES
Data We adopt the same setup as in the GAN experiment by Pope et al. (2021). We use BigGAN (Brock al., 2019) pretrained on ImageNet at 128 128 resolution (so that 128 128 3 49152), downloaded from https: generate an image an ImageNet category, we provide the corresponding one-hot class vector as well as noise vector to the GAN. Following Poole et al. (2019), we control the intrinsic dimension d of the generated images by setting all except the first dimensions of the 128-dimension noise vector to zero, and use truncation level of 0.8 for the sample diversity factor. Since the GAN generates values in [-1, 1]. we rescale them linearly to [0, 1] to correspond to images, so that the alphabets X' y [0. 1 As in (Pope al., 2021). experimented with images from the ImageNet class basenji. i. In Figure and 10, we plot additional random samples, for 2 and d
fhub. dev/
/1. To
Upper Bound Experiments We implemented a version of ResNet- VAE following the appendix of Kingma et al. (2016), using the 3 x 3 convolutions of Cheng et al. (2020) in the encoder and decoder our model. Our model consists of 6 layers of latent variables, and implements bidirectional inference using factorized Gaussian variational distributions at each stochastic layer. During inference pass, an input image goes through stages of convolution followed by downsampling (cach time reducing the height and width by 2), and results in stochastic tensor each stage. In encoding
24

---

Published as a conference paper at ICLR 2022
Tecnick
Quality (PSNR)
88
36
34
proposed (D( (ResNett-VA)
32
proposed Rs (D) (Minnen 2020 (-VAE)
Minnen 2018 (hyperprior)
30
0.2
0.3
01S5
0.5
0.6
0.7
0.8
0.9
1.0
Rate (BPP, bits per sample per pixel)
Figure 11: Quality-rate curves on the Tecnick dataset.
order, the stochastic tensors have decreasing spatial extent and an increasing number of channels equal to 4, 16, 32, 64, and 128. The latent tensor Zo at the topmost generative hierarchy flattened and modeled by MAF prior Qz. (aided by the fact that all the images have fixed shape).
Lower Bound Experiments We parameterize the log u model by simple feedforward convo- lutional neural network. It contains three convolutional layers with 4, 8, and 16 filters, each time downsampling by 2, followed by fully connected layer with 1024 hidden units. Again we use Sel.U activation in the hidden units, and set h 2048 in all the experiments.
A.5.7 NATURAL IMAGES
Data For signals such images or video, which can have arbitrarily high spatial dimension, it is not immediately clear how to define i.i.d. samples. But since our focus is on image compression, follow the common practice of training convolutional autoencoders random 256x 256--ixeei crops high resolution images as in learned image compression research (Ballé et al., 2017), noting that current methods cannot effectively exploit spatial redundancies larger than this scale (Minnen & Singh, 2020). We can then use the trained models to estimate R-D upper bounds on the image- generating distributions underlying standard benchmark datasets such Kodak Kodak (1993) and Tecnick Asuni & Giachetti (2014).
As representative dataset of natural images, we used images from the COCO 2017 (Lin al., 2014) training set that are larger than 512x512, and downsampled each by a random factor in [0.6, 1] to remove potential compression artifacts from JPEG. We trained image compression models from (Minnen et al., 2018; Minnen & Singh, 2020) on our dataset and were able to closely reproduce their performance.
Upper Bound Experiments Instead of working with variational distributions over the set of pixel values y X = {0, 1, 1,.. 255}". which would require specialized tools for representing and optimizing high-dimensional discrete distributions (Salimans al., 2017; Hoogeboom al., 2019; Maddison et al., 2016; Jang et al., 2016), here we parameterize the variational distributions in Fuclidean latent spaces for simplicity (see Sec. 3).
We applied our upper bound algorithm based on two VAE architectures:
1.
The learned image compression model from Minnen & Singh (2020) 6 This corresponds to a hierarchical VAE with two layers of latent variables, with the lower level latents modeled
'https://gt wwtpbb..mtteess///////o f1 ow www.o//oroossooo/////o/oc http////looo//ccc/////// f9edo 199a55358 ffa fa (farfa5aaacd377escceeeeee models/ms20200
25

---

Published as a conference paper at ICLR 2022
by a channel-wis autoregressive prior conditioned the top-level latents. We adapt it for our purpose, by using factorized Gaussians for the variational posterior distributions Qzix (instead of factorized uniform distributions used to simulate rounding in the original model), and no longer convolving the variational prior Qz with uniform distribution as in the original model.
2. The ResNet- VAE from our experiment on GAN images. Following learned image compres- sion models, we maintain the convolutional structure of the topmost latent tensor (rather than flattening it in the GAN experiments) in order to make the model applicable to images variable sizes. Instead a MAF, we model the topmost latent tensor by deep factorized (hyper-)prion Qz, proposed by Ballé et al. (2018).
For the 3-VAE based on Minnen & Singh (2020), dim( = 0.28 di'()"); and for our 3-ResNet- VAE, dim(2) 0.66 di'(X)). We did not observe improved results by simply increasing the number of latent channels/dimenssonn in the former model. We follow the same model training practice as reported by Minnen & Singh (2020), training both models around 5 million gradient steps.
In our problem setup, we consider the reproduction alphabet to be the discrete set of pixel values A' = 1, ..., 255)". the decoder network w needs to discretize the output of the final convolutional layer (in Tensor f flow, this can be implemented by satu saturate cast to uint B). However, since the discretization operation has no useful gradient, we skip during training, and only apply for evaluation.
The operational quality-rate curves of various compression methods are taken from tensorf] flsorlloo fww.wwooommm cww-compees inn? and compressA libraries. The results of Minnen & Singh (2020) represent the current state-of-the-art of neural image compression, to the best of our knowl- edge.
Lower Bound Experiments When running our lower bound algorithm on 256 x 256 crops of natural images, we observed similar behavior as on high-dimension Gaussians, with the loss & quickly plateauing value around log k. This issue persisted when tried different model architectures. See a discussion on this phenomenon in Section A.5.2.
A.6 MEASURES OF VARIABILITY
As alluded to Sec. 3, the proposed R-D bounds are estimated as empirical averages from samples, which only equal the true quantities in expectation. Before we characterize the variability of the bound estimates, we first define the exact form of the estimators used.
Upper Bound Estimator For the proposed upper bound, consider the variational problem solved by the Lagrangian Eq. 3, and define the random variables from the rate and distortion terms of the Lagrangian:
R(Z, x) == (
log g()),
and
D(2,X)= (xx,)(Xxw)))
Above, X and Z follow the joint distribution Px Px the decoder function, and q(>(>)) and g(=) are the Lebesgue density functions of absolutely continuous Qz x=r and Q2 variational distributions (we only worked with absolutely continuous variational distributions for simplicity). Then is clear that the expectations of the two random variables equal the rate and distortion terms, i.c.,
E [R, X)) Eyg- Px [L 1Qz)l.
and
E[D(Z,X)) = Epx Q (((X,(()))
Recall for any given decoder and variational distributions, the above rate and distortion terms define point (E[)), that lies upper bound R, (D) of the source R(D) curve. Given i.i.d.
besutts/image_comppeee httpm//ooooo//ccmmmresssoo/o///e F322 9386188 (essb118113331e3333
7 Phttps://ght wttp://sooo//eoocceee/oeeooo/tot
httpg/ittaattnntttaeeeeeeeee hcc/CompressAl//reee bac8ef00000F12111000a22222222222
26

---

Published as a conference paper at ICLR 2022
[RR.
8R.
[[R] CI
[AD
8D
[B[D] CI
300
0.278
0.8
(0.15, 0.41)
0.006
0.00287
(0.00553, 0.00647)
500
1.25
1.27
(1.04, 1.46)
0.00351
0.00207
(0.00317, 0.00385)
le+03
2.69
1.21
(2,49, 2.89)
0.00157
0.00103
(0.0014, 0.00174)
4c+03
4.77
1.25
(4.56, 4.97)
0.000352
0.000285
(0.000305, 0.000399)
8c+03
5.65
1,31
(5,44, 5.87)
0.000196
0.000152
(0.000171, 0,000221)
1.6e+04
6.43
1.32
(6.21, 6.64)
0.000126
9.68e-05
(0.00011, 0.000142)
Table 1: Statistics for the estimated R-D upper bound on GAN-generated images with d 2,
computed with m 100 samples of (R. D). "CT" denotes 95% large-sample confidence interval.
A
PR
SR
[[ CI
JtD
$D
B[D) CI
300
0.265
0.689
(0.15, 0.38)
0.00978
0.00376
(0.00916, 0.0104)
500
1.44
1.24
(1.23, 1.64)
0.00693
0,00284
(0.00646, 0.00739)
le+03
4.08
1.73
(3.80, 4.36)
0.00332
0.00152
(0.00307, 0.00357)
2e+03
6.59
1.88
(6.28, 6.90)
0.00151
0.000752
(0.00139, 0.00164)
4e+03
8.66
1.89
(8.35, 8.97)
0.000774
0.000389
(0.00071, 0.000838)
8e+03
10.4
1.89
(10,06, 10.68)
0.000476
0.000229
(0.000438, 0.000514)
1.6e+04
12
1.83
(11.69, 12.29) 0.000323
0.000147
(0.000299, 0.000347)
Table 2: Statistics for the estimated R-D upper bound GAN-generated images with d 4.
computed with m 100 samples of (R, D). "CI" denotes 95% large-sample confidence interval.
pairs of (21, X1), X2), (2m Xm) PxQxxx. estimate such point by (PD, (RD, using the usual sample mean estimators PR == mm R(Z, X,) and up 1 Enn D(Z X,).
Lower Bound Estimator For the proposed lower bound, given i.d. X1, ..., X ~ Px. define the random variable
1 Σ log u(X.) { :
Ck logg +
a
where the various quantities are defined in Sec. 4. Then it is clear that its expectation equals the lower bound objective Eq. 8, i.e., [[{] la == Ei-log u (X)) E/Ca]/a log a + Given m i..dd. samples of &, we form the usual sample mean estimator #ę == 1 ES Ej, and form an R-D lower bound estimate by
R2(D) == -DD + Me.
In computing μęç, we fix the constant a to a sample mean estimate of [[Ca] computed separately (as in line of Algorithm 1). As explained in Sec. A.4, repeat this procedure with different X, and our overall R-D lower bound obtained by taking the point-wise maximum of such linear under-estimators over
Results Below we give measures of variability of for the sandwich bound estimates obtained on GAN generated images (Section 6.3). We report detailed statistics the following tables, and plot them in Figure 12. Our results on other experiments have comparable variability.
For the upper bound, given each trained (AVEE, we compute m samples of (R, essentially passing samples of X through the autoencoder and collect the rate and distortion values; then report the sample mean (μ), sample standard deviation (s), and 95% large-sample confidence interval for the population mean (the true expected value) for R and D, respectively, in Tables and 2.
For the lower bound, given each trained log u model, we compute m samples of &, then report the sample mean (ug), sample standard deviation (sg). and 90% large-sample confidence lower bound for the population mean (the true expected value E[{] fx) in Tables and
27

---

Published as a conference paper at ICLR 2022
X
Agg
8g
[ LCB
99.63
0.74
0.01
0.74
199.90
1.46
0.01
1.45
299.15
2.04
0.02
2.03
499.69
2.78
0.02
2.78
999.01
3.84
0.03
3.83
1999.70
4.58
0.05
4.57
2024.00
4.60
0.06
4.58
2999.30
4.97
0.05
4.95
3036.80
4.96
0.08
4.93
4000.00
5.18
0.05
5.17
Table 3: Statistics for the estimated R-D lower bound on d GAN-generate images, com- puted with m = 30 samples of &. "LCB" denotes 90% large-sample lower confidence bound.
12
Minnen 2018 (hyperprior)
Minnen 2020
10
% proposed (D)
Rate (nats per sample)
proposed R.(D)
90% conf. lower bound of proposed (D)
0.000
0.002
0.004
3,006
0.008
0.010
Distortion (mean squared error)
Figure Zoomed-in version of Figure Middle-Bottom (GAN-generated images with in- trinsic dimension ple mean estimates of R-D upper bound points (LDD, (R) bracketed by 95% confidence intervals (red), well lower bound estimates based on the sample mean jlg (maroon, dashed) and its 90% confidence lower bound (pink line, beneath maroon).
4), showing the sam-
X
Ag
8g
E[ LCB
99.63
1.09
0.01
1.09
199.90
2.18
0.02
2.17
299.15
3.08
0.02
3.07
499.69
4.44
0.04
4.43
999.01
5.97
0.06
5.95
1999.70 70
6.69
0.08
6.67
2024.00
6.69
0.07
6.66
2999.30
6.87
0.08
6.84
3036.80
6.85
0.09
6.83
4000.00
6.93
0.08
6.91
Table 4: Statistics for the estimated R-D lower bound on d GAN-generated images, com- puted with m = 30 samples of &. "LCB" denotes 90% large-sample lower confidence bound.
28


Reviews:
[{"review_id": "tnZW0Zq9qSc", "reviewer": "Reviewer_9Pwc", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I cannot find any strength. I will outline major weakness, not covering all the waeknesses (as all pages are weak).\n\n(1) The paper claims to propose upper and lower bounds on the R-D function of unknown memoryless information source. There are always nonzero probability of the \"upper bound\" being lower than the true R-D function, because the data samples are probabilistic. The standard tool for handling such a situation is one-sided confidence interval (see any undergraduate textbook on frequentiest statistics).\nThere is no reason for not using the standard statistical framework.\nIf one wants to use the Baysian counterpart of the interval estimation, one must have a prior probability of the unknown memoryless information source, but such a prior probability is not provided in the manuscript.\n\n(2) The paper gives the big claim in the first page that \"bounding the R-D function of a general (i.e. discrete, absolutely continuous or neither), unknown memoryless source\". There are several subtleties, for example: The expectations of random variables do not\nalways exist, for example the Cauchy distribution. What is the definition of the R-D function for such probability distributions?\nThe expectation is used in almost every displayed equations of the manuscript. What if the expectation is undefined, like the case of the\nCauchy distribution? If the authors assume the existence of expectation in every mathematical expression of the manuscript,\nwhat kind of the memoryless information source can be handled by the proposal? It is clear that the discrete and Gaussian sources\nare OK, but they can be handled by the Arimoto-Bluhat algorithm in 1970s, after estimating the probability distribution by standard\nstatistical methods.\n\n(3-1) The manuscript uses mathematical theorems that have not been proved for \"general\" memoryless information source.\nFor example, Lemma A.1 cites Csiszar and Gray, which proved Lemma A.1 for discrete and Gaussian sources.\nThe manuscript must prove Lemma A.1 for the class of information source that are considered.\n\n(3-2) Similar for Theorem A.2. Please provide a proof of Theorem A.2 for the \"general\" information source.\n\n(4) With any conficence interval estimation algorithm, there always exists a positive probability of the event that the true parameter does not belong\nto the confidence interval generated by the algorithm with data samples obtained. The confidence level is one minus such a \"failure\"\nprobability. Thus, every practical interval estimation algorithm accompanies with (desirably theoretical) derivation of the confidence level.\nThe manuscript does provide confidence levels that can handle \"general\" memoryless information source.\nIt is acceptable that confidence levels of the proposed \"upper\" and \"lower\" bounds are also algorithmically produced given\ndata samples.\n\n(5) Confidence levels should be evaluated with the finite number n of data samples, otherwise one cannot know what is the probabilities\nof the proposed \"lower\" and \"upper\" bounds fail in a practical situation.\n\n(6) If the authors want to retain their claim that their proposed algorithm can handle any probability distribution, please provide\nthe results of running their proposed lower and upper bounds algorithms on the Cauchy distribution as defined at \nhttps://en.wikipedia.org/wiki/Cauchy_distribution\n\n(7) Related to (6), it is unclear what kind of random variables/probability distributions can be handled by the proposed algorithms.\nIt must be clearly stated. Experiments are only conducted for the Gaussian distribution and so-called the \"banana\" source.\nAre the proposed algorithms can be applied to other classes of probability distributions? If the authors claim wider applicability\nof their proposed algorithm, it must be supported by experimental or theoretical evidences. If the proposals are only applicable\nto a limited class of distributions, it is unsuitable for ICLR.", "summary_of_the_paper": "The paper claims to propose \"upper bound\" and \"lower bound\" on the R-D function in lossy data compression.", "summary_of_the_review": "Mathematical reasoning in the paper is doubtful, and what kind of situations can be handled by the proposed\nalgorithm is unclear, see \"main review\" for the detail.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "1: The contributions are neither significant nor novel."}, "extra_text": {}}, {"review_id": "jXqraypjBA7", "reviewer": "Reviewer_zsSF", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is well-motivated and has valuable results for the lossy compression community. In more detail:\n\n$\\textbf{Strenghts:}$ \n- The paper is well-written and summarizes the prior work well for readers from both information theory and generative models backgrounds. \n\n- The authors provide a tight upper bound on the RD function of arbitrary sources without limiting the sources to discrete data or restricting the dimensions. I believe this is a useful contribution to the lossy compression field. For instance, as the authors pointed out in the paper, we can assess the success of lossy compression algorithms by comparing their performance with tight RD bounds.\n\n-  The authors use the dual formulation of the RD function to find a lower bound, which, to the best of my knowledge, is a novel approach.\n\n- The experimental results on random Gaussian samples show the exactness of their upper bound (for this particular source). For the banana-shaped source and its high-dimensional projection, the authors show that upper and lower bounds are tighter than what was proposed in prior work. These two sets of experiments provide empirical support for the theoretical results in the paper. \n\n- I find the experimental results on natural images interesting as the authors essentially set a (somewhat loose) limit for future work in lossy image compression without actually providing a compression algorithm. \n\n$\\textbf{Weaknesses:}$ \n- Despite the flexibility of the proposed method for finding upper bounds on RD function, the authors restrict the source to have low effective dimension with continuous reproduction alphabet. This makes the current lower bound not applicable to high-dimensional data such as images.\n\n- The lower bound is evaluated only on a banana-shaped source. It would be nice to see the tightness of the upper and lower bounds on a few other sources. \n\n- I could not access the code in the given link. The google drive folder seems empty.  \n\n- A few minor points:\n    - Section 4, Dual Characterization of R(D) paragraph, \"..., a variational lower would require...\" misses the word \"bound\"\n    - Section 6.1, last sentence of the second paragraph: optional - - > optimal\n    - There are a few problems with the Appendix. The text for the last section A.5.3 Natural Images seems missing. The figure number in the last paragraph is also missing.   ", "summary_of_the_paper": "This paper aims to provide stronger upper and lower bounds for the RD function of arbitrary sources. The authors handle unknown sources by requiring only i.i.d. samples. Specifically, the authors derive an upper bound to the RD function using a $\\beta$-VAE-like generative model, which has some similarities to the Blahut-Arimoto (BA) algorithm with two distinctions: 1) they do not restrict the source to be a low-dimensional discrete data, and 2) they use gradient descent to predict the parameters of the variational distributions $Q_{\\hat{X}|X}$ and $Q_{\\hat{X}}$. They also derive a lower bound to the RD function using the dual form of the RD function. Finally, the authors provide experimental results to show that the proposed upper bound is, in fact, exact for random Gaussian samples. For more complex data such as banana-shaped sources and their high-dimensional projections, they show that the proposed upper and lower bounds are tight. Furthermore, the proposed upper bounds on the RD function of natural images indicate that the state-of-the-art image compression methods can still be improved by at least one dB PSNR.", "summary_of_the_review": "Although the proposed method has limitations, I find some of the contributions valuable. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "icmyguqc7iw", "reviewer": "Reviewer_T3TM", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Novelty: The classical Blahut-Arimoto algorithm requires the restrictive assumption of a discrete source as well as a known source distribution. This paper aims to develop a technique for unknown sources with possible continuous support. While there has been some work on extending this in the past (Riegler et al, 2018; Harrison & Kontoyiannis, 2008; Gibson, 2017), none have tackled the problem to the breadth of this paper. The paper is well written and easy to follow.\n\nSignificance: On the technical side the main contribution of the paper is on computing lower bound for rate-distortion function, as the upper bound is a rather natural application of the ELBO analysis in VAEs. While the experiments on the Gaussian and Banana sources give credence to the tightness of the bounds, these are rather simple examples. In particular, the variational distribution used for the experiments is Gaussian and contains the source distribution in its family. Although the banana source is a nonlinear transform of a Gaussian source, by visual inspection the underlying inverse transform seems to be easily learnable by an encoder in principle.  Do the authors have examples of synthetic experiment in which the source distribution is significantly different from the variational/latent distributions (Q_{Z|X}, Q_Z) as well? \n\nThe results on image compression suggesting a possible 1 PSNR improvement seem to be significant. \nAs the authors have mentioned, it is unclear whether we can actually close this gap due to the nature of it being an overestimate. Furthermore, the lower bound computation, which seem to be the main technical innovation does not scale to such datasets. As a result the reviewer feels that the technical contribution of the paper in realistic image compression applications may be limited due to these reasons. \n\n", "summary_of_the_paper": "This paper considers the problem of estimating the rate-distortion function R(D) for arbitrary sources with unknown distribution. An upper bound is established using a variational distribution and is estimated using an iterative coordinate descent algorithm inspired by the Blahut-Arimoto algorithm. A lower bound is established using a parameterization of the dual form of the rate-distortion function by Cisszar. These are then evaluated using neural network encoder/decoders on Gaussian and Banana sources, and tested on natural images suggesting that there is room for improved compression rates achieved by current state-of-the-art image compression algorithms.\n\n", "summary_of_the_review": "Due to the weaknesses in the \"significance\" part of the main review I am marking the score as a 5.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "IZJpaJGSJqV", "reviewer": "Reviewer_tH6q", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper proposes to use ML to establish lower and upper bounds on rate distortion for general sources, thus going beyond the Blahut-Arimoto algorithm (which assumes discrete sources with known PMF, and has complexity exponential in the dimension of the data). Some theoretical results are provided, extending prior results. Some numerical results are provided.\n\nTo construct the upper bound, the paper starts with the unconstrained variational objective of Blahut-Arimoto.  The distributions Q(\\hat{x} | x) and Q(\\hat{x}) are parametrized by the parameters of a NN; minimizing the objective function provides a point on the upper bound of the R(D) curve. A key insight is the observation that this objective function is similar to NELBO in \\beta-VAE.  Based on this a projection to a lower-dimensional space is considered. The paper establishes R_\\omega(D) \\ge R(D), where \\omega is the decoder function. While this seems straightforward, there is no discussion of the complexity in finding each R(D) point, and of course the required sweep over \\lambda to find the entire UB curve. There is little discussion on the choice of \\omega; in particular what choices are likely to lead to tighter bounds. \nThe establishment of the lower bound is less clear. The paper notes that this much harder than the UB case. Again there is little detail on the architecture for computing u_\\theta(x). At the end of section 4, we are told that Appendix A has the detailed algorithm, but Appendix A says that an outline will be provided\n\nNumerical section: insufficient details are provided so that it would not be possible for an informed reader to reproduce the results in the paper. . While some detail is provided for the Gaussian sources (results shown in Fig 2), very little detail is provided for the banana sources (Fig 3) or the natural images (Fig 4). \n\nSection A.5.3 “Natural Images” is a blank section \n\nOther: \n\nWhat does “1 PSNR” (used in the abstract and several other places mean)? Does this mean 1 dB improvement in PSNR? \n\nI assume a diagonal Gaussian process is a vector Gaussian process in which the components are independent so that the covariance matrix is diagonal? \n\n", "summary_of_the_paper": "The paper proposes to use ML to establish lower and upper bounds on rate distortion for general sources, thus going beyond the Blahut-Arimoto algorithm (which assumes discrete sources with known PMF, and has complexity exponential in the dimension of the data). Some theoretical results are provided, extending prior results. Some numerical results are provided ", "summary_of_the_review": "The use of ML to address complex problems in information theory is appealing. The authors consider the problem of obtaining upper and lower bounds for the rate distortion curve for source distributions that may be continuous, discrete, or mixed, and may not be known (only data are available). Some theoretical results are provided. Two major issues: insufficient details about the learning architectures are provided, so that an informed reader will not be able to reproduce the results in the paper. Second, there is very little discussion of the complexity of the algorithm. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "Not applicable"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper proposes an algorithmic approach to estimating upper and lower bounds of the rate-distortion (R-D) function of a data source on the basis of samples drawn from it. The proposed upper bound is based on the variational objective employed in the Blahut-Arimoto algorithm, whereas the proposed lower bound is based on the dual characterization of the R-D function. In both bounds neural networks trained with samples are utilized. Experimental results on four sources (Gaussian, banana-shaped, GAN-generated images, and natural images) are provided.

The four review scores were initially two positives and two negatives. Some reviewers evaluated positively the argument on the lower bound of the R-D function. On the other hand, one reviewer showed his/her concern about lack of the argument on statistical confidence of the obtained bounds. In response, the authors have addressed it in Section A.6 in the Supplementary Materials (SM) of the revised manuscript with the experiment using GAN-generated images (high-dimensional data with low intrinsic dimension), and the results are summarized in Tables 1-4 and Figure 10 in SM. The authors have in their revision also provided a specification for the range of the sources to which the proposal would be applicable. Description of the experiments on images, which was missing in the initial version as pointed out by some reviewers, has been added in the revised manuscript. Still, as some reviewers mentioned, the main weakness is that the proposal failed to demonstrate its usefulness to estimate the lower bound in settings where the data dimension is truly high, as in the experiment with natural images, where statistical confidence analysis was not conducted either. Three reviewers have revised their respective scores upward after the author response.

Despite some weaknesses I think that this paper provides a novel and interesting algorithmic approach to estimating the rate-distortion function. I would therefore like to recommend acceptance of this paper, and would like to encourage the authors to perform confidence analysis also for the experiments on natural images.
Decision: Accept (Poster)