Greedy Poisson Rejection Sampling
Authors: Gergely Flamich
Venue: NeurIPS 2023 poster (2023)
Keywords: channel simulation
relative entropy coding
reverse channel coding
rejection sampling
Poisson process
TL;DR: We propose rejection sampling procedures using Poisson processes, dual to A* sampling in a certain sense. Then, we use it to develop algorithms to optimally compress samples from any distribution, which is useful for neural data compression.

Abstract: One-shot channel simulation is a fundamental data compression problem concerned with encoding a single sample from a target distribution $Q$ using a coding distribution $P$ using as few bits as possible on average. Algorithms that solve this problem find applications in neural data compression and differential privacy and can serve as a more efficient and natural alternative to quantization-based methods. Unfortunately, existing solutions are too slow or have limited applicability, preventing their widespread adaptation. In this paper, we conclusively solve one-shot channel simulation for one-dimensional problems where the target-proposal density ratio is unimodal by describing an algorithm with optimal runtime. We achieve this by constructing a rejection sampling procedure equivalent to greedily searching over the points of a Poisson process. Hence, we call our algorithm greedy Poisson rejection sampling (GPRS) and analyze the correctness and time complexity of several of its variants. Finally, we empirically verify our theorems, demonstrating that GPRS significantly outperforms the current state-of-the-art method, A* coding.

Greedy Poisson Rejection Sampling
Gergely Flamich Department of Engineering University of Cambridge gf3320cam.c ac.
Abstract
One-shot channel simulation is fundamental data compression problem concerned with encoding a single sample from a target distribution Q using a coding distribu- tion P using as few bits possible average. Algorithms that solve this problem find applications in neural data compression and differential privacy and can serve methods. Sadly, existing solu-
more efficient alternative to quantization-bast tions are too slow or have limited applicability, preventing widespread adoption. In this paper, we conclusively solve one-shot channel simulation for one-dimensional problems where the target-proposal density ratio is unimodal by describing an algo- rithm with optimal runtime. We achieve this by constructing rejection sampling procedure equivalent to greedily searching over the points of a Poisson process. Hence, we call our algorithm greedy Poisson rejection sampling (GPRS) and ana- lyze the correctness and time complexity of several of its variants. Finally, we em- pirically verify our theorems, demonstrating that GPRS significantly outperforms the current state-of state-of-the of-the art method, A coding. Our code is available at httpps:/ github com/ wwwmcomm////eeee herggerrggglggylee wlly-fllmicii ch/ cageeedy//eeyyyyyy/e poisson rejection pereeciion-samplnng..
Introduction
It is a common misconception that quantization is essential to lossy data compression; in fact, it is merely a way to discard information deterministically. In this paper, we consider the alternative, that is, to discard information stochastically using one-shot channel simulation. To illustrate the main idea, take lossy image compression as an example. Assume we have generative model given by a joint distribution Px.y images y and latent variables x, e.g. might have trained variational autoencoder (VAE; Kingma & Welling 2014) on a dataset of images. To compress a new image Pxly as its stochastic lossy representation. The decoder can obtain a lossy reconstruction of y by decoding x and drawing sample ý ~ Pylx (though in practice, for VAE we normally just take the mean predicted by the generative network).
y, encode single sample from its posterior
Abstracting away from our example, in this paper we will be entirely focused on channel simulation for pair of correlated random variables x, Px.y: given source symbol y ~ Py we wish to encode a single sample x ~ Pxly- A simple way to achieve this is to encode x with entropy coding using the marginal Px. whose average coding cost approximately the entropy H[x]. Surprisingly, however, we can do much better by using a channel simulation protocol, whose average coding cost approximately the mutual information / /[x; y) (Li El Gamal 2018). This remarkable, since not only [x: y] < H[[x], but in many cases / [x; y] might be finite even though [H[[x] is infinite, such as when x is continuous. Sadly, most existing protocols place heavy restrictions on Px.y or their runtime scales much worse than O(/ y)), limiting their practical applicability (Agustsson & Theis 2020
In this paper, we propose family of channel simulation protocols based on new rejection sampling algorithm, which we can apply to simulate samples from a target distribution Q using proposal distribution P over an arbitrary probability space. The inspiration for our construction comes from an exciting recent line of work which recasts random variate simulation as a search problem over set
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

15
x
x
x
x
10
X
x
X
5
x
*
x
x
x
X
x
x
x
0
-1
0
0
X
X
X
Figure I: Illustration of three GPRS procedures for Gaussian target Q - (1, 0.25*) and Gaussian proposal distribution = (0, 1), with the time axis truncated to the first 17 units. All three variants indicated
find the first arrival of the same (1, P)-Poisson process II under the graph of 4p by the thick dashed black line in each plot. Here, r = dQ is the target-propose density ratio, and a is given by Equation ( Left: Algorithm sequentially scarching through the points of II. The green circle shows the first point of II that falls under yp and is accepted. All other points rejected, indicated by red crosses (X). In practice, Algorithm does not simulate points of IT that arrive after the accepted arrival. Middle: Parallelized GPRS (Algorithm 4 searching through two independent (1/2, P)-Poisson processes II; and ΠI2 in parallel. Blue points are arrivals in II and orange points are arrivals in ll2. Crosses (X) indicate rejected, and circles indicate accepted points by each thread. In the end, the algorithm accepts the earliest arrival across all processes. which in this case is marked by the blue circle (0). Right: maanchandd-oundd GPRS (Algorithm) 55) when o is unimodal. The shaded red areas are never searched or simulated by the algorithm since, given the first two rejections, we know points in those regions cannot fall under yp
of randomly placed points, specifically a Poisson process (Maddison 2016). The most well-known examples of sampling-as-searde are the Gumbel-max trick and A sampling (Maddison et al. 2014). Our algorithm, which we call greedy Poisson rejection sampling (GPRS), differs significantly from all previous approaches in terms of what it is searching for, which
succinctly summarise as: "GPRS searches for the first arrival of a Poisson process 11 under the graph of an appropriately defined function yp" The first and simplest variant of GPRS is equivalent to an exhaustive search over all points of II in time order. Next, we show that the exhaustive search is embarrassingly parallelizable, leading to parallelized variant of GPRS. Finally, when the underlying probability space has more structure, we develop branch-ar anddbouudd variants of GPRS that perform a binary search over the points of II. See Figure [ for an illustration of these three variants.
While GPRS is an interesting sampling algorithm on its own, we also show that each of its variants induces a new one-shot channel simulation protocol. That is, after we receive y ~ Py. we can set Q + Pxly and P + Px and use GPRS to encode a sample x ~ Pxly at an average bitrate of a little more than the mutual information / [[x; y]. In particular, on one-dimensional problems where the density ratio / /dPx x is unimodal for all (which is often the case in practice), branch-and- bound GPRS leads to a protocol with an average runtime of (()x which is optimal. This is considerable improvement over A coding (Flamich et al. 2022). the current state-of-the- th-the-art method.
In summary, our contributions are as follows:
We construct new rejection sampling algorithm called greedy Poisson rejection sampling. which we can construe as greedy search over the points of Poisson process (Algorithm We propose a parallelized (A (Algorithm + and a branch-and-bour 6) of GPRS. We analyze the correctness and runtime of these algorithms.
variant (Algorithms
We show that each variant of GPRS induces one-shot channel simulation protocol for correlated random variables x, y ~ Px.y. achieving the optimal average codelength of 7[x; + log, [x: 1) + (1) bits.
We prove that when x is a REvalued random variable and the density ratio dPxty/ /dPx is always unimodal, the channel simulation protocol based on the binary search variant of GPRS achieves [(([x;y]) runtime, which is optimal.
We conduct toy experiments on one-dimensional problems and show that GPRS compares favourably against A* coding, the current state-of-the-an art channel simulation protocol.

---

2 Background
The sampling algorithms we construct in this paper are search procedures on randomly placed points in space whose distribution is given by a Poisson process 11. Thus, in this section, we first review the necessary theoretical background for Poisson processes and how we can simulate them on a computer. Then, we formulate standard rejection sampling as a search procedure over the points of a Poisson process serve as prototype for our algorithm in the next section. Up to this point, our exposition loosely follows Sections 2, 3 and 5.1 of the excellent work of Maddison (2016), peppered with few additional results that will useful for analyzing our algorithm later. Finally, we describe the channel simulation problem, using rejection sampling as rudimentary solution and describe its shortcomings. This motivates the development of greedy Poisson rejection sampling in Section 3
2.1 Poisson Processes
A Poisson process II is a countable collection of random points in some mathematical space Ω. In the main text, we will always assume that II is defined over Ω R+ x Rd and that all objects involved are measure-theoree theoretically well-behaved for simplicity. For this choice of Q the positive reals represent time, and TRd represents space. However, most results generalize to settings when the spatial domain RRd is replaced with some more general space, and we give general measure-theoretic construction in Appendix [ where the spatial domain is an arbitrary Polish space.
Basic properties of TI: For a A C Ω2 let ((A) def [IT A| denote the number of points of IT falling in the set A, where |-| denotes the cardinality of a set. Then, II is characterized by the following two fundamental properties (Kingman 1992). First, for two disjoint sets A, B S Ω2 A n B = 0, the number of points of II that fall in either set are independent random variables: ((A) 1 N(B)
Second, (A) is Poisson distributed with mean measure ((A) def [[)
Time-ordering the points of II: Since we assume that l + x Rd has a product space structure, we may write the points of II as a pair of time-space coordinates: II = we can order the points in II with respect to their time coordinates and index them accordingly, i.c. for i < j we have T; < T Hence, we refer to (Ta, Xx)) as the nihh arrival of the process.
Furthermore,
As a slight abuse of notation, we define N(t) def N(10,) x Rd) and s(t) def E[N()), i.e. these quantities measure the number and average number of points of II that arrive before time i, respectively. In this paper, we assume jut) has derivative n'(t), and as- sume for each t 2 0 there is a conditional probability distribution PxiT=r with density p(x t), such that we can write the mean measure as m(A) = SAP(E | t)) da dt.
Simulating II: A simple method to simulate II on com- puter is to realize it in time order, i.e. at step n, sim- ulate Tn and then use it to simulate X,, PXIT=T. We can find the distribution of II's first arrival by not- ing that no point of II can come before it and hence where the
Generating a (), )xir)-Poisson process.
Algorithm 1:
Input :Time rate X, Spatial distribution Px|r
P(T, 2 t] = P[N(t) = 0] = exp(-p(t) second equality follows since N(t) Poisson distributed with mean (ut). A particularly important case is when II is tmm-hoongeneouu, i.e. (u(t) At for some O, in which case Ti ~ Exp(A) is an exponential random variable with rate In fact, all of IT's inter-arrival times A, Tn T Tn-1 for > share this simple distribution, end where we set To = 0. To see this, note that
To < 0 for n 1,2, do
An Exp(A)
Ta + Tn- + An
Xn ~ PXIT-T, yield Xn)
PIA, > 1 T... 1== = [N (TTn
R2) = Tw-1] = exp( exp(-Xt), Xt),
Le. all An Tu-1 Exp(A). Therefore, can the above procedure to simulate time- homogeneous Poisson processes, described in Algorithm We will refer to a homogeoomges Poisson process with time rate À and spatial distribution Px|r (A, ) process.
Rejection sampling using IT: Rejection sampling is a technique to simulate samples from a target distribution Q using a proposal distribution P. assuming can find an upper bound M > 0 for
3

---

Algorithm 2: Standard rejection sampler.
Algorithm 3: Greedy Poisson rejection sampler.
Input Proposal distribution P. Density ratio r = =QQ/PP, Upper bound M for r.
Input :Proposal distribution P.
Density ratio r =
Stretch function c.
Output :Samle X ~ Q2 and its index N.
Output :Samle X ~ Q and its index N.
// Generator for a (1, (Poisson process using Algorithm 1
// Generator for a (Posssno process using Algorithm 1 II SimulatePP(1,P)) for = 1.2 do
II + SimulatePP(1,
for n = 1,2...... do
XA, t (Ix
(T,, f- [I)
Un ~ ([nn((((
T,, < a (r(X,)) then return X,,
U, < <(Xx)/M then return X,, end end
end
end
their density ratio == dQ/dPpQQpppaoo (technically, the Radon-Nikodym derivative). We can formulate this procedure using a Poisson process: we simulate the arrivals Xn) of (1, P)-Poisoon process II. but we only keep them with probability r (M, otherwise, we delete them. This algorithm is described in Algorithm
[ its correctness guaranteed by the thinning theorem (Maddison) 2016))
Rejection sampling is suboptimal: Using Poisson processes to formulate rejection sampling high- lights subtle but crucial inefficiency: it does not make use of IT's temporal structure and only uses the spatial coordinates. GPRS fixes this by using a rejection criterion that does depend on the time variable. As we show, this significantly speeds up sampling for certain classes of distributions.
2.2 Channel Simulation
The main motivation for our work to develop a one-shot channel simulation protocol using the sampling algorithm we derive in Section 3 Channel simulation is of significant theoretical and practical interest. Recent works used to compress neural network weights, achieving otate-of-the-art
performance (Havasi et al. 2018); to perform image compression using variational autoencoders (Flamich et al. (020)) and diffusion models with perfect realism (Theis al.) 2022); and to perform
differentially private federated learning by compressing noisy gradients (Shah et al. (022).
One-shot channel simulation is a communication problem between two parties, Alice and Bob, sharing a joint distribution Px,y over two correlated random variables x and y, where we assume that Alice and Bob can simulate samples from the marginal Px. In single round of communication. Alice receives a sample ~ Py from the marginal distribution over y. Then, she needs to send the minimum number of bits to Bob such that he can simulate a single sample from the conditional distribution x ~ Pxly- Note that Bob does not want to learn Pxty: he just wants to simulate a single sample from Surprisingly, when Alice and Bob have access shared randomness, e.g. by sharing the seed of their random number generator before communication, they can solve channel simulation very efficiently. Mathematically, in this paper we will always model this shared randomness by some time-homogeneou Poisson process (or processes) II, since given a shared random seed, Alice and Bob can always simulate the same process using Algorithm [ Then, the average coding cost of x given IT is its conditional entropy H|[x TI] and, surprisingly. y. it is always upper bounded by H[[x II] < 1 [[x:y yl + log, [x; y) + 1) + (1), where [[x; y] is the mutual information between x and y (Li & El Gamal 2018). This is especially curious result, given that in many cases H[[x] is infinite while 1 [x;y] yJ is finite, e.g. when x is a continuous variable. In essence, this result means that given the additional structure Pxxy. channel simulation protocols can "offload" an infinite amount of information into the shared randomness, and only communicate the finitely many "necessary" bits.
An example channel simulation protocol with rejection sampling: Given II and y Py. Alice sets Q <- P as the target and P <- Px as the proposal distribution with density ratio r = dQQddP and run the rejection sampler in Algorithm 2 to find the first point of IT that was not deleted. She counts the number of samples N she had to simulate before acceptance and sends this number to Bob. He can then decode sample Pxly by selecting the spatial coordinate of the Nith arrival of II.

---

Unfortunately, this simple protocol is suboptimal. To see this, let Dxo [Q][ def supxen (log2 r(x) denote Rényi oo-divergence of Q from P. and recall two standard facts: (I) the best possible upper bound Alice can use for rejection sampling is Moye = exp2 (D [Q[[)))) where expy(#) = 2 and (2), the number of samples N drawn until acceptance is geometric random variable with mean Maps (Maddison
We now state the two issues with rejection sampling that GPRS solves.
Problem 1: Codelength. By using the formula for the entropy of geometric random variable and assuming Alice uses the best possible bound Mopt in the protocol, we find that
(a) H[[ IT] Ry [HIN yy] 2 Ry-r, [Dxo Pxly [P]] 2 T[x;y,,
(2) see Appendix [orr the derivation. Unfortunately, inequality (a) can be arbitrarily loose, hence the average codelength scales with the expected oc-divergence instead of [xx; y). as would be optimal.
Problem 2: Slow runtime. We are interested in classsfying the time complexity of our protocol. As saw, for a target Q and proposal P. Algorithm 2 draws Mopt exp2 (D (QI[P)) samples average. Unfortunately, under the computational hardness assumption RP & NP, Agustsson & Theis (2020) showed that without any further assumptions, there is sampler that scales polynomially in DKL Q PI. However, with further assumptions, we can do much better, as we show in Section 3.1
3 Greedy Poisson Rejection Sampling
We now describe GPRS; its pseudo-code is shown in Algorithm 3 This section assumes that Q and P are the target and proposal distributions, respectively, and r = ddPdP is their density ratio. Let II be a (1, PJ-Poisson process. Our proposed rejection criterion is now embarrassingly simple: for an appropriate invertible function of the composite function & = a o as illustrated in the left plot in Figure | We refer to o as the stretch function for r, as its purpose is to stretch the density ratio along the time-axis towards so.
R + R+ accept the first arrival of II that falls under the graph
Deriving the stretch function (sketch, see Appendix for details): Let yp = o r, where for now is an arbitrary invertible function R+, let U [(1,x) E Ω 1 S [p()}} be the set of points under the graph of yp and let Il 11 n U. By the restriction theorem (Kingman 1992). II is also a Poisson process with mean measure ((A)) = ((AN)) Let (ī,X) be the first arrival of II, i.e. the first arrival of IT under oo and let Q, be the marginal distribution of X. Then, as we show in Appendix Al the density ratio dQ=/dP is given by
dQ.
dP (z)
-
P[T 2 P[T21]d
(3)
Now we pick a such that Q. Q, for which we need to ensure that dQ. =. Substituting
t= ((xr) into Equation (3). and differentiating,
get (o-1) (1) P 2 t]. Since (T, X) falls
under the graph of yo by definition, we have PP 2 f| = P[ 2 (X)) 2 o (())]. By expanding the definition of the right-hand side, obtain time-invariant ODE for a -1.
(o-1)== = we (o- - o -1 wp (a-1), with o -() 0,
where (h) daf Py PZ)p(r(Z) 2 h] and we defined analogously. Finally, solving the ODE, we get a (Th)
1 wo(1) - 17 wp(n) dn.
Remember that picking a according to Equation (5) ensures that GPRS is correct by construction.
To complete the picture, in Appendix A we prove that
always exists and Algorithm terminates with probability 1. We now turn our attention to analyzing the runtime of Algorithm 3 and surprisingly, we find that the expected runtime of GPRS matches that of standard rejection sampling; see Appendix B.2 for the proof. Note that in our analyses, we identify GPRS's time complexity with the number of arrivals of IT it simulates. Moreover, we assume that we can evaluate a in O(1) time.
Theorem 3.1 (Expected Runtime). Let Q and P be the target and proposal distributions for Al- gorithm 3 respectively, and r = dQ/dP their density ratio. Let N denote the number of samples simulated by the algorithm before it terminates. Then,
E[N]= = exp2 (D)) and V[N] 2 exp2 (Doo PQ)(
(6)
where V[:] denotes the variance of random variable.
5

---

Furthermore, using ideas from the work Li & Verdu (0118), we can show the following upper bound on the fractional moments of the index; see Appendix [B.3 for the proof.
Theorem 3.2 (Fractional Moments of the Index). Let Q and P be the target and proposal distributions for Algorithm 3 respectively, and dQ/ dP their density ratio. Let N denote the number of samples simulated by the algorithm before it terminates. Then, for a E (0,1) we have
Do (211))
B[N - 1 1+
exp2 (aD (aD,=IQQPP))
log2 e
where D=[Q T log2 Bz-q [r(Z)a-1] is the Rényi a-divergence of Q from P.
GPRS induces a channel simulation protocol: We use an analogous solution to the standard rejection sampling- based scheme in Section 2.2 to obtain channel simulation protocol for pair of correlated random variables x. y Px.y. First, we assume that both the encoder and the decoder have access to the same realization of a (1, Px)-Poisson process II, which they will use as their shared randomness. In a single round of communication, the encoder receives y Py. sets Q + Pxly as the target and P <- Px as the proposal distribution, and runs Algorithm 3 to find the index of the accepted arrival in II. Finally, following Li & EI Gamal (018)), they use a zeta distribution <nn s) x n with s [xx: y] + logg to encode N using entropy coding. Then, the decoder simulates a I Pxly- -isttribtt sample by looking up the spatial coordinate XN of ll's Nth arrival. The following theorem shows that this protocol is optimally efficient; see Appendix B.4 for the proof.
Theorem 3.3 (Expected Codelength). Let Px,y be a joint distribution over correlated random variables x and y. and let II be the (1, Px)-Poisson process used by Algorithm 3 Then, the above protocol induced by the the algorithm is optimally efficient up to a constant, in the sense that
H x
II| < y[ + loga y] + 1) + 6.
(8)
Note, that this result shows we can use GPRS as an alternative to the Poisson functional representation to prove of the strong functional representation lemma (Theorem 1: Li & EI Gamal 2018)
GPRS in practice: The world applicability of Algorithm 3 is limited by the following three computational obstacles: (1) GPRS requires exact knowledge of the density ratio /QPPPP to evaluate s2, which immediately excludes a many practical problems where dQ/dP is known only up to a normalizing constant. (2) To compute a, we need to evaluate Equation (4) or (5). which requires us to compute wpp and wiq, which is usually intractable for general dlQ/aPP.Fortumaeely
Fortunately, in some special provide the expressions for discrete, uniform, triangular, Gaussian, and Laplace distributions. (3) However, even when we can compute wp and wQ, analytically evaluating the integral in Equation (5) is usually not possible. Moreover, solving it numerically is unstable as is unbounded. Instead, we numerically solve for o- using Equation in practice, which fortunately turns out to be stable.
cases they can be computed analytically: in Appendix G
However, we emphasize that GPRS's ruison d' ê'êre is to perform channel simulation using the protocol we outline above or the ones we present in Section 3.1 not general-purposs random variate simulation. Practically relevant channel simulation problems, such as encoding a sample from the latent posterior of VAE (Flamich et al., 2020) or encoding variational implicit neural representations (Guo et al.) 2023 He et al., (223) usually involve simple distributions such as Gaussians. In these cases, we have readily available formulae for wp and uq and can cheaply evaluate 0- via numerical integration.
GPRS as greedy search: The defining feature of GPRS is that its acceptance criterion at each step is local, since if at step the arrival (Ta,X,) in II falls under y,, the algorithm accepts it. Thus is greedy search procedure. This is in contrast with A sampling (Maddison et al. 2014). whose acceptance criterion is global, the acceptance of the with arrival (7.,, X,, depends on all other points of II in the general case. In other words, using the language of of & Verdu (2018). GPRS is a causal sampler, as its acceptance criterion each step only depends on the samples the algorithm has examined in previous steps. On the other hand, A* sampling is non-causal as its acceptance criterion in each step depends on future samples as well. Surprisingly, this difference between the search criteria does not make difference in the average runtimes and codelengths in the general case. However, as we show in the next section, GPRS can be much faster in special cases.
3.1 Speeding up the greedy search
This section discusses two ways to improve the runtime of GPRS. First, we show how we can utilize available parallel computing power to speed up Algorithm 3 Second, we propose advanced search

---

Algorithm 4: Parallel GPRS with J available threads.
Algorithm 5: Branch-and-bound GPRS I with unimodal r
Input Proposal distribution P. Density ratio r = dQ/dP. Stretch function Ø, Number of parallel threads J.
Input : :Proposal distribution P. Density ratio r = dQ/dP, Stretch function Location xx' of the mode r.
Output:Sampll X ~ Q and its code Ny,, T", ,",
Output:Sample :aample X ~ Q and its heap index H.
N-- ,N= < oc, mil, nil,nil nil in parallel for j 1,...,J. J do
To: H, B <- (0, 1,)
11; 4 SimulatePP(1/4,9) for ni 2,,..... do
for d = 1.2....de do
Xa ~ PlB
(T() x(1) + next (II,)
Aa ~ Exp (P(B))
By
Ta + Ta-1+ Ad
i < T then terminate thread j.
if < o(r(Xx)) then
return X2, end
end
if T < a (r (xx² (i)
if Xa 2 x² then
then
B - B (- H t 2H
T ,", TMp x11 j.ny terminate thread j.
else
end end
B - Bn (Xd. oc)
H 2H +1
end
end
return x", (N"
end
strategy when the spatial domain 52 has more structure and obtain a asuper-exponential improvement in the runtime from exp; (Doo PI) ((II)) in certain cases.
Parallel GPRS: The basis for parallelizing GPRS is the superposition theorem Kingman 1992)
j=1 NT is (1, P)-Poisson process. As shown in Algorithm 4 this result makes parallelizing GPRS very simple given parallel threads: First, we independently look for the first arrivals of II,..., II, under the
For a positive integer J, let I . II, all be [1 P)- P)-Poisson processes; then, II =U=-
N. x respectively, where the Nj corresponds to the index of ITy's 's first arrival under sp. Then, we select the candidate with the earliest arrival time, i.e.
graph of ip yielding (T()) x!))
N,
j = arg min min N2- Now, by the superposition theorem, (TM) x())) the first arrival of II under şp, and hence x(1")
TT))
N.- ~ Q. Finally, Alice encodes the sample via the tuple (j,,,,,)) i.e. which of the J processes the first arrival occurred in, and the index the arrival in ny.. See the middle graphic in Figure for example with J parallel threads.
Our next results show that parallelizing GPRS results a linear reduction in both the expectation and variance of its runtime and more favourable codelength guarantee for channel simulation.
Theorem 3.4 (Expected runtime of parallelized GPRS). Let Q, P and r be defined as above, and let V denote the random variable corresponding to the number of samples simulated by thread j in Algorithm 4 Jusinn J threads. Then, for all j.
(9) Theorem 3.5 (Expected codelength of parallelized GPRS). Let Px.y be joint distribution over correlated random variables x and y, and let II,....,III II, be Px)J,Px)-Ponl -Poisssnn processes. Then, assuming log, J < y| parallelized GPRS induces a channel simulation protocol such that
E[py]] = (exp, (Dxc
and Viy1 2 (exp, 1P=(
H[ .I. [I,] < [[x:y] + loog (f[x,y] - logg J +
(10)
See Appendix C for the proofs. Note that we can use the same parallelisation argument with the appropriate modifications speed up A sampling/ coding too.
Branch-and-boume GPRS on We briefly restrict our attention to problems on when the density ratio> is unimodal, can exploit this additional structure to more efficiently search for II's 'ss first
7

---

arrival under s. Consider the example in the right graphic in Figure[] we simulate the first arrival (Ti,X)) .x). and reject it, since Ti (. Since r is unimodal by assumption and increasing, op is also unimodal. Hence, for x < X1 we must have s(()) < oo ()))), ), while the arrival time of any of the later arrivals will be larger than Ti. Therefore, none of the arrivals to the left of X1 will fall under yp either! Hence, it is sufficient to simulate IT to the right of X1, which we can easily do using generalized inverse transform sampling ! We repeat this argument for all further arrivals to obtain an
efficient search procedure for finding II's first arrival under y, described in Algorithm 5 We simulate the next arrival (TB, XB) of II within some bounds B, starting with = S2. If S (((B) we accept: otherwise, where Xg falls relative to yp's mode. We repeat these two steps until we find the first arrival.
truncate the bound to B t- B n (-oc, XB), B B n (Xg, 00) based on
Since Algorithm 5 does not simulate every point of II, we cannot use the index N of the accepted arrival to obtain a channel simulation protocol as before. Instead, we encode the search path, i.e. whether we chose the left or the right side of our current sample at each step. Similarly to A coding. we encode the path using its heap index (Flamich et al. 2022): the root has index Hroot = 1, and for node with index H, its left child is assigned index 2// and its right child 2// + 1. As the following theorems show, this version of GPRS is, in fact, optimally efficient; see Appendix E for proofs.
Theorem 3.6 (Expected Runtime of GPRS with binary search), Let Q.P P and be defined as above, and let D denote the number of samples simulated by Algorithm 5 and X = 1/ 1og2(4/3) Then, B S X P + O1))
(11)
Theorem 3.7 (Expected Codelength of GPRS with binary search). Let Px, be a joint distribution over correlated random variables x and y, and let II be (1, P=)- process. Then, GPRS with binary search induces a channel simulation protocol such that
H[x II] S I y[ + 2 log2 (fxxy) + +1.
(12)
Branch-and-boumt GPRS with splitting functions: With some additional machinery, Algorithm 5 can be extended to more general settings, such as RD and cases where is not unimodal, by
introducing the notion of a splitting function. For a region of space, B C S2 a splitting function split simply returns binary partition of the set, i.e. {L, Rt} split(B), such that L R = 0 and L J R- B. In this case, we can perform a similar tree search to Algorithm [ captured in Algorithm 6 Recall that (T, X) is the first arrival of I1 under s. Starting with the whole space B = f2 we simulate the next arrival (T, X) of 11 in B at each step. we reject it, we partition / into two parts, L and R. using split. Then, with probability P|X E R X E B,T T), we continue searching through the arrivals of II only in R, and only in L otherwise. We show the correctness of this procedure in Appendix F and describe how to compute PIX E R X E B, 2 T] and in practice. This splitting procedure is analogous to the general version of A* sampling/coding (Maddison et al. 2014: Flamich et al., (022)). which is parameterized by the same splitting functions. Note that Algorithm is special case of Algorithm where Ω R. is unimodal, and at each step for an interval bound B- = (a, b) and sample X E (a, b) we split B into (a, X),(X,)) (X,b)).
Experiments
We compare the average and one-shot case efficiency of our proposed variants of GPRS and couple of other channel simulation protocols in Figure 2 See the figure's s caption and Appendix for details of our experimental setup. The top two plots in Figure 2 demonstrate that our methods expected runtimes and codelengths align well with our theorems predictions and compare favourably to other methods. Furthermore, we find that the mean performance is a robust statistic for the binary search-basee variants of GPRS in that lines up closely with the median, and the interquartile range is quite narrow. However, we also see that Theorem 3.7 is slightly loose empirically. We conjecture that the coefficient of the (Tg,([[x;; + 1) term in Equation 12 could be reduced to 1, and that the constant term could be improved as well.
On the other hand, the bottom plot in Figure 2 demonstrates the most salient property of the binary search variant of GPRS: unlike all previous methods, its runtime scales with PP and not Doo [Q][P]. Thus, we can apply it to larger family of channel simulation problems, where DRI. 1Q [P] is finite but P! is large even infinite, and other methods would not terminate.
'Suppose we have a real- real-valued random variable X with CDF F' and inverse CDF F-1, Then, we can simulate X truncated to jo, b) by simulating U ~ Un) F(b) and computing X (U).
8

---

Algorithm 6: Branch-and- GPRS with splitting function.
40
GPRS
# of steps
30
PER
Input : Proposal distribution P. Density ratio r /QPdP. Stretch function 0, Splitting function split. Output: Output:Sample X ~ Q and its
AS*
sGPRS
10
dGPRS
T[x; [4]
heap index H To. H. B < (0, R)
UB
Codelength
15
for d = 1, ..
X, PlB
Aa ~ Exp (P())
5
Ta < Ta- Ad
0
iT < (r((a)) then
2
4
8
10
return Xa, H
[[x: 14]
else
Bp, t split(B)
30
sGPRS
gt
P[X E B1 X E B,A Tal
# of steps
20
AS*
DKL
3 + Bernoull(p))
10
Doo
H 2H +
B + Bg
end end
10
15
20
25
Dz 1Q11PP
Figure 2: Left: Binary search GPRS with arbitrary splitting function. Right: Performance compari- son of different channel simulation protocols. In each plot, dashed lines indicate the mean, solid lines the median and the shaded areas the 25 75 percentile region of the relevant performance metric. We computed the statistics over 1000 runs for each setting. Top right: Runtime comparison on a ID Gaussian channel simulation problem Px,,. plotted against increasing mutual information [[x; μ4]. Alice receives fa ev (0, o²) and encodes sample μu ~ ((,,, 1) to Bob. The abbreviations in the legend are: GPRS A Algorithm 33 PFR Poisson functional representation Global-bound A* oplit-on-sample A coding (Flamich iet al. 2022); sGPRS split-on-sample GPRS (Algorithm 5 dGPRS - GPRS with dyadic split (Algorithm 6). Middle right: Average codelength comparison of our proposed algorithms on the same channel simulation problem as above. UB in the legend corresponds to an upper bound of 1) bits. We estimate the algorithms' expected codelengths by encoding the indices returned by GPRS using Zeta distribution ((n X) x n with À = 1 1/1 y;y in each case, which is the optimal maximum entropy distribution for this problem setting (Li & El Gamal 2018). Bottom right: One-shot runtime comparison of sGPRS with AS* coding. Alice encodes samples from target Q = N(m, s²) using P = (N 1) as the proposal. We computed m bits for each problem, but Doc [Q[[P] increases. GPRS' runtime stays fixed as it scales with DRL P), while the runtime of A keeps increasing.
coding (Li El Gamal 2018 Flamich al. 2022). AS* * A5**
7[x: |4] [
and s2 such that DKL [?
5 Related Work
Poisson processes for channel simulation: Poisson processes were introduced the channel simulation literature by Li & El Gamal (2018) via their construction of the Poisson functional representation (PFR) their proof of the strong functional representation lemma. Flamich al. (2022) observed that the PFR construction is equivalent to certain variant of A sampling (Maddison) all.222114))) Maddison 2016). Thus, they proposed an optimized version of the PFR called A coding. which achieves O(D [Q)P)) runtime for one- dideemmmeiioll unimodal distributions. GPRS was mainly inspired by A* coding, and they are dual constructions to each other in the following sense: Depth- A * coding can be thought of as imporlance sampler, i.e. a Monte Carlo algorithm that returns an approximate sample in fixed time. On the other hand, GPRS a rejection sampler, i.e. Las Vegas algorithm that returns exact sample in random time.

---

Channel simulation with dithered quantization: Dithered quantization (DQ: Ziv 1985 is an alternative to rejection and importance sampling-based approaches to channel simulation. DQ exploits that for any E IR and U, U ~ Unif(-1/2, 1/2), the quantities [ 01 + U and + U" are equal in distribution. While DQ has been around for decades as a tool to model and analyze quantization error, Aussssson & Theis (2020) reinterpreted it as a channel simulation protocol and used it to develop a VAE- neural image compression algorithm. Unfortunaeely, basic DQ only allows uniform target distributions, limiiing its applicability. As partial remedy, Theis & Yosri (2022) showed DQ could be combined with other channel simulation protocols to speed them up and thus called their approach hybrid coding (HQ). Originally, HQ required that the target distribution be compactly supported, which was lifted by Flamich & Theis (203)), who developed an adaptive rejection sampler using HQ. In a different vein, Henazyy & Li (2022) generalize and analyze method proposed in the appendix of Aguutsssn & Theis] (2020) and show that DQ can be used to realize channel simulation protocols for one-dimensional symmetric, unimodal distributions.
Concurrently to our work, Flamich et al. (2023) introduce greedy rejection coding (GRC) which generalizes Harsha et al. ss rejection sampling algorithm to arbitrary probability spaces and arbitrary splitting functions, extending the work Flamich & Theis (2023). Furthermore, they also utilize a space-partitioning procedure to speed up the convergence of their sampler and prove that a variant of their sampler also achieves optimal runtime for -imeennnsi problems where dQ/ /d dPP is unimodal. However. the construction of their method differs significantly from ours. GRC is a direct generalization of Harsha et al. (2007)'s 's algorithm and, thus, a more "conventional" rejection sampler, while we base construction on Poisson processes. Thus. our proof techniques are also significantly different. It is an interesting research question whether GRC could be formulated using Poisson processes, akin to standard rejection sampling in Algorithm 2 as this could be used to connect the two algorithms and improve both.
Greedy Rejection Coding:
6 Discussion and Future Work
Using the theory of Poisson processes, we constructed greedy Poisson rejection sampling. We proved the correctness of the algorithm and analyzed its runtime, and showed that it could be used to obtain channel simulation protocol. We then developed several variations on it, analyzed their runtimes, and showed that they could all be used to obtain channel simulation protocols. As the most significant result of the paper, we showed that using the binary search variant of GPRS we can achieve (( [Q[[P) Pl) runtime for arbitrary mindimensional, unimodal density ratios, significantly improving upon the previous best ((Dx<QII)) bound by A coding.
There are several interesting directions for future work. From a practical perspective, the most pressing question is whether efficient channel simulation algorithms exist for multivariate problems; finding an efficient channel simulation protocol for multivariate Gaussians would already have far- reaching practical consequences. We also highlight interesting theoretical question. We found that the most general version of GPRS needs to simulate exactly eDpp(Dxo (Q)[ (QI samples expectation, which matches A* sampling. From Agustsson & Theis (020),, we know that for general problems, the average number of simulated samples has to be at least on the order of exp, (DxL. Q IP)) Hence we ask whether this lower bound can be tightened to match the expected runtime of GPRS and A* sampling or whether there exists an algorithm that achieves a lower runtime.
7 Acknowledgements
We are grateful to Lennie Wells for the many helpful discussions throughout the project, particularly for suggesting the nice argument for Lemma D.2, Peter Wildemann for his advice on solving some of the differential equations involved in the derivation of GPRS and Stratis Markou for his help that aided our understanding of some high-level concepts so that we could derive Algorithm 6 We would also like to thank Lucas Theis for the many enjoyable early discussions about the project and for providing helpful feedback on the manuscript; and Cheuk Ting Li for bringing Liu & Verduls work to our attention. Finally, we would like to thank Daniel Goe for proofreading the manuscript and and providing feedback on the early version of
Lorenzo Bonito for helping to design Figures [ and [ the manuscript. GF acknowledges funding from DeepMind.
10

---

References
Firikur Agustsson and Lucas Theis. Universally quantized neural compression. In Advances in Neural Information Processing Systems, volume 33, 2020.
Harry Bateman and Arthur Erdélyi. Higher transcendental functions, volume 1. McGraw-Hill book company, 1953.
Robert M. Corless, Gaston H. Gonnet, David E. G. Hare, David Jeffrey, and Donald E. Knuth. On the Lambert W function. Advances in Computational mathematics, 5:329-359 1996.
Gergely Flamich and Lucas Theis. Adaptive greedy rejection sampling. In IEEE International Symposium on Information Theory (ISIT), pp. 454-459, 2023.
Gergely Flamich, Marton Havasi, and José Miguel Hernández Lobato. Compressing images by encoding their latent representations with relative entropy coding, Advances in Neural Information Processing Systems, 33, 2020.
Gergely Flamich, Stratis Markou, and José Miguel Hernández Lobato. Fast relative entropy coding with A coding. In International Conference on Machine Learning, pp. 6548-6577, 2022.
Gergely Flamich, Stratis Markou, and José Miguel Hernández. Lobato. Faster relative entropy coding with greedy rejection coding. In Advances in Neural Information Processing Systems, 2023.
Zongyu Guo, Gergely Flamich, Jiajun He, Zhibo Chen, and José Miguel Heráández Lobato. Compres- sion with Bayesian implicit neural representations. Advances in Neural Information Processing Systems, 2023.
Prahladh Harsha, Rahul Jain, David McAllester, and Jaikumar Radhakrishnan. The communica- tion complexity of correlation. In Twenty-Secowe Anmual IEEE Conference on Computational Complexity (CCC'07), pp. 10-23. IEEE, 2007.
Marton Havasi, Robert Peharz, and José Miguel Hernández Lobato. Minimal random code learning: Getting bits back from compressed model parameters. International Conference on Learning Representations, 2018.
Jiajun He, Gergely Flamich, Zongyu Guo, and José Miguel Hernández Lobato. RECOMBINER: Robust and enhanced compression with Bayesian implicit neural representations. arXiv preprint MXXi:23009.17182 2023.
Mahmoud Hegazy and Cheuk Ting Li. Randomized quantization with exact error distribution. In 2022 IEEE Information Theory Workshop (ITW), pp. 350-355. IEEE, 2022.
Diederik P. Kingma and Max Welling. Auto-encoding variational Bayes. International Conference on Learning Representations, 2014.
John F. C. Kingman. Poisson Processes. Oxford Studies in Probability. Clarendon Press, 1992. ISBN 9780191591242.
Cheuk Ting Li and Abbas El Gamal. Strong functional representation lemma and applications to coding theorems. IEEE Transactions on Information Theory, 64( (4(11))6967-6978 2018.
Jingbo Liu and Sergio Verdu. Rejection sampling and noncausal sampling under moment constraints. In 2018 IEEE International Symposium on Information Theory (ISIT), pp. 1565-1569. IEEE, 2018
Chris Maddison. Poisson process model for Monte Carlo. Perturbation, Optimization, and Statistics, pp. 193-232, 2016.
Chris J. Maddison, Daniel Tarlow, and Tom Minka. A sampling. Advances in Neural Information Processing Systems, volume 27. pp. 2086-3094, 2014.
Pat Muldowney, Krzysztof Ostaszewski, and Wojciech Wojdowski. The Darth Vader rule. Tatra Mountains Mathematical Publications, 52(1):33-63, 2012.
Jorma Rissanen and Glen Langdon. Arithmetic coding. IBM Journal of research and development, 23(2):149-162, 1979.

---

Abhin Shah, Wei-Ning Chen, Johannes Balle, Peter Kairouz, and Lucas Theis. Optimal compression of locally differentially private mechanisms. In International Conference on Artificial Intelligence and Statistics, pp. 7680-7723. PMLR, 2022.
Lucas Theis and Noureldin Yosri. Algorithms for the communication of samples. In International Conference on Machine Learning, 2022.
Lucas Theis, Tim Salimans, Matthew D. Hoffman, and Fabian Mentzer. Lossy compression with Gaussian diffusion. arXiv:2206.08889, 2022.
Jacob Ziv. On universal quantization. IEEE Transactions on Information Theory, 31(3):344-347, 1985.
12

---

A Measure-theoretic Construction of Greedy Poisson Rejection Sampling
In this section, we provide construction of greedy Poisson rejection sampling (GPRS) in its greatest generality. However, we first establish the notation for the rest of the appendix.
Notation: In what follows, will always denote GPRS's target distribution Q and its proposal distribution as P. We assume that Q and P are Borel probability measures over some Polish space Q
with Niioon derivative def dQ/dP. We denote the standard Lehesgue measure on Ik" by X. All logarithms are assumed to be to the base 2 denoted by logg- Similarly, we use the less common notation of exp2 to denote exponentiation with base 2, i.e. exp2(x) = 2". The relative entropy Q from defined as P def Ja log2 r(z) (1) and the Rényi oo-divergence as Doo [Q[PP] def ess suszeR (zogr(2)), where the essential supremum is taken with respect to P.
Restricting Poisson process "under the graph of a function": We first consider the class of functions "under which" we will be restricting our Poisson processes. Let
+ IR+ be a measurable function, such that (o(())] ([()], i.e. the image of yp has the same cardinality as the image of = dQ/dP. This implies that there exist bijections between r(2) and (12). Since both images are subsets of RR+ both images have linear order. Hence, a monotonically increasing bijection ā exists such that o = ā o Since ā is monotonically increasing, we can extend it to a monotonically increasing continuous function ø E + R This fact significant, to reduce questions about functions of interest from Q to IR to invertible functions on the positive reals. Since a is continuous and monotonically increasing on the positive reals, it is invertible on the extended domain. We now consider restricting a Poisson process under the graph of and work out the distribution of the spatial coordinate of the first arrival.
it allows us
Let II = (TT,, X,, )m=
be a Poisson process over R+ x f2 with mean measure μ = x P. Let def ((t.z) E Ω S (()) be the set of points "under the graph" of y. By the restriction theorem
KKingman 1992), π def INn U is a Poisson process with mean measure ū(AA) = ((An U) for an arbitrary Borel set A. As a slight abuse of notation, we define the shorthand ja() (() x ?) for the average number of points in II up to time Let Nç denote the counting measure of Il and let (T, X) be the first arrival of Î1, i.e. the first point of the original process that falls under the graph of sp, assuming that exists. To develop GPRS, we first work out the distribution of X for an arbitrary ip.
Remark: Note, that the construction below only holds the first arrival guaranteed exist, otherwise the quantities below do not make sense. From formal point of view, we would have
to always condition on the event (xc) def lim (t) > O, e.g. we would need to write P[T t (co) 0j. However, make "leap of faith" instead and assume that our construction valid to obtain our "guesses" " for the functions a and y. Then, we show that for our specific guesses the derivation below is well-defined, i.e. namely that P[Ni > 0] 1, which
then implies P t (oo) > = P( t). To do this, note that
P|N > 0] = 1 P|N (00) = 0]
(13)
lim -a()) t-
(14)
Hence, for the well-definition of our construction it is enough to show that for our choice of (i(t) + xc as t + sc.
Constructing a: We wish to derive
dP[ = x
000 d x = x]
dt.
(15)
dP
d(A P)
J0
13

---

To do this, recall that we defined (t) def N (10, t) x 12) and define (s, t) def [n((s, x (2).
Similarly, let p((t) def (t)) and jis,t) t) def
n b) - ju(t) - fi(s). Then,
P E A. T E [s t]]
P[X E A, A, E dt)
lim 8-22
(16)
s-s
P[X E A. n (s, 1) (,)) - 0]
lim 8-2t
(17)
8-
P[X E N (s, 1JPIN IJP[N 0
lim
(18)
1-8
where the last equality holds, because (n(+)) 1 Niss, (), since (0, x S2 and [s. x Q are disjoint. By Lemma from Maddison (2016),
P[X A (s, 1] ([(TTXX] [s.t] x A (s,ttt
(19)
[ t| x A)
(20)
m(s,t)
Substituting this back into Equation (18 and using the fact that N is Poisson, we find
p| A, E di)
lim s-2t
M(ss t| A) i ()
Mis D) ((at)) E E - 8)
/(1-s)
(21)
[ū((s,) 1] x A) ( e-B(o)) /(1 /(--s)
lim
(22)
c[[ t| X A)
=e-((t) alt)) lim g-3
(23)
t-s
/
-(t) 1 dP(r),
(24)
where the last equation holds by the definition of the derivative and the fundamental theorem of calculus. From this, by inspection we find
dPP t = X x] d(A x P)
1lt s (()))( A(0) = 1/t S (()) 2 t].
(25)
Finally, by marginalizing out the first arrival time, we find that the spatial distribution is
dP [X
dP T t, X x
dt
(26)
dP
d(> x P)
(p(())
P[T dt.
(27)
Deriving by finding a: Note that Equation algorithm, we need to set it such that dP[&=x]
27 holds for any However, to get a correct
dP aP r. Since we can write o = a r by our earlier argument, this problem is reduced to finding appropriate invertible continuous function
a + + R Thus, we wish find a such that
(r(x))
Ve E S2
r(x)
P( 2 1] di.
(28)
Now, introduce (f(())), so that (r) r(æ). Note that this substitution only makes sense for T E r(2). However, since we are free to extend o to R in any we like so long as it is monotone, we may require that this equation hold for all 7 E R.+ Then, we find
/"" F(T 2 Adt
o )
(29)
= (o-1))
PTT 2 r]
(30)
14

---

with 0 000) 0. Before we solve for J, we define
wp(h) def Py~plr(Y)) 2 h]
(31)
w def (y-olr(Y) 2 h]
(32)
Note that wp and wQ are supported on [0, r² where r" d = Now, we can rewrite Equation (30)
0 10 (t) = P(T 2 t]
(33)
= P(T 2 (X) 2 $ + 21, <
(34)
=0 due mutual exclusivity
PIT 2 0, (x) 2 t
(35)
=P[p(X)22) since? T 2 always
-(r( 2 1 (0) P -,
(36)
= Lo - o (t) (t)
(37)
where the second term in the last equation follows by noting that
P(r(X) 2 o (t) 2 L.L. 1r(() 2 o 2
2 7] dr dP(x)
1[((a)2 (O), since TIC
(38)
= 1|r(x) 2 (+)] ( P[ 2 r] dr dP(x)
(39)
29
= (t)up( (t)).
(40)
Now, by the inverse function theorem, we get
o' d'(h) (h(h) h wp(h)
(41)
thus we finally find
bb Jo (Q( 77 wp(n) dn.
1
(hh)
(42)
Thus, to recapitulate, setting o = a o r, where a is given by Equation 42 will ensure that the spatial
distribution of the first arrival of II under yp is the target distribution
An alternate form for P[T 2 t|: We make use of the following integral representation of the indicator function to present an alternative form of the survival function of 7:
for a > 00: 1[a - x]
5(y a) dy,
(43)
where d is the Dirac delta function. Then, we have
we(h) (h) = / 1/r(x) 2 dxx))
(44)
1 %'m r(z)0(y- r(x)) r(x)) dyd)(x)
= -1-
(45)
L/" [ 6(y r(x)) dl(z) dy
(46)
- Iv(1 wo(w)) + S (1 (a-wr(u)) dy
(47)
1 + wp(h))) - // wP(() dy.
(48)
15

---

Then, from this identity, we find that
P(T 2 t = w@ ( (t)) - - (o
(49)
1 (t)
wp(y) dy.
(50)
Well-definition of GPRS and Algorithm 3 terminates with probability 1: To ensure that our construction is useful, we need to show that the first arrival under the graph T exists and is almost surely finite, so that Algorithm terminates with probability First, note that for any L > O, we have j((t) S p(t) P (t) = 0) > O. Second, note that for h E [0,r")
oc. Since P[ i (t) 0] = e-((t) this implies that for any finite time
wh(h) - h wp(h)
1[r(x) 2 (1(r(a)-1)oal) - k)dP(x),
(51)
from which
ww(r") - r'wp(r") =
1|r(x) 2 r'((((x) - r²) dP(x)
(52)
1/((x) -+)
(53)
0.
(54)
Thus, in particular, find that
lim e e((((())) = lim (n))) 0] lim (wq (h) - h wp(h)) 0.
(55)
h-r
h-rr
By continuity, this only hold if ((((()) monotonically increasing, this also implies that ohh + oc as - r". Thus, we have shown a couple of facts:
h - r". Since ja()) is finite for all t and
0o h + r-, hence & is always unbounded points in Q that achieve the supremum r". Furthermore, this implies that
(h)
o (t)
cc.
(56)
Since - is increasing, this implies that is bounded from above by r".
((t) < oc for all t > 0, but (it) + 0o as t + co. Hence, by Equation 14 we have P[N i (c)) > 0] 1. Thus, the first arrival of II exists almost surely, and our construction is well- In particular, it is meaningful write P[T 2 t).
P > 1 P(N (1) 0] + 0 / > co, which shows that the first arrival time is finite with probability 1. In turn, this implies that Algorithm 3 will terminate with probability one, desired.
Note, that wp and uq can be computed in many practically relevant cases, see Appendix G
B Analysis of Greedy Poisson Rejection Sampling
Now that we constructed a correct sampling algorithm in Appendix A we turn our attention to deriving the expected first arrival time T in the restricted process II, the expected number of samples N before Algorithm 3 terminates and bound on Nss variance, and an upper bound on its coding cost. The proofs below showcase the true advantage of formulating the sampling algorithm using the language of Poisson processes: the proofs are all quite short and elegant.
B.1 The Expected First Arrival Time
At the end of the previous section, we showed that T is finite with probability one, which implies that E[T] < co. Now, we derive the value of this expectation exactly. Since T is a positive random variable, by the Darth Vader rule (Muldowney et al. 2012). may write its expectation
E 7 =
P 2 t dt
lim o -1(())
56
(57)
100
16

---

B.2 The Expectation and Variance of the Runtime
Expectation of N: First, let Ii def II \ II be the set of points in II above the graph of yp Since II and II are defined on complementary sets, they independent. Since by definition T' is the first arrival of II and the Nith arrival of I1, it must mean that the first N - I arrivals of II occurred in i1; Thus, conditioned on T, we have N - H X) mean
C T T is Poisson distributed with
E N -
1
1
/
1 2 ((x)) dP(x) di
(58)
I-" J
T
1 - 1 < x(x)] dP(x) dt
(59)
= 1 T
(60)
By the law of iterated expectations,
BIN - 1| = Ej E
N - 1 [
(61)
= E -
E j
(62)
Focusing on the second term, we find
/ ju(t) Pf E dt) dt
E
j0
T
(63)
p(t) p (t)e mit)) dt
(64)
ue- du
(65)
where the third equality follows by substituting u = ja(t). Finally, plugging the above and Equa- tion (57p into Equation 62 we find
E N = 1| = 1 - -1='"
(66)
Variance of N: We now show that the distribution of N is super-Poissonian, i.e. [N[] S V[N]. Similarly to the above, we begin with the law of iterated variances to find
V[N] V[N 1] Ey
[
+ +
V 7
1
T
(67)
T
(68)
where the second equality follows from the fact that the variance of N - 1 matches its mean conditioned on T, since it is a Poisson random variable. Focussing on the last term, we find
T 2 E {T)
T
E
(69)
j(tt) PP[ E dt) dt
(70)
w' dt - 1=
(71)
where the third equality follows from a similar u-substitution as in Equation 65 Thus, putting everything together, we find
V[N] = ---r--++ -1+ +1=++1+-0 E[N].
(72)
17

---

B.3 The Fractional Moments of the Index
In this section, we follow the work Liu Verdu (2018) and analyse the fractional moments the index N and prove Theorem 3.2 As before, set Q and P as GPRS's ' target and proposal distribution, respectively, and let r ap- To begin, we define the Rényi a-divergence between two distributions for a E (0, 1) U o0) as
dQ
Da[QIIP] = log2 Ez~Q [+(2)a-]]
(73)
By taking limits, we can extend the definiton to a = 1, 00 as
P| def
lim D. Po[ PKIPQ
(74)
P def
lim P) = css sup (log2 r()) =E
(75)
a0000
Now fix a E (0, and x E Q2 and note that by Jensen's inequality,
E[ X = x] s EIN X x]"
(76)
Now, continue, we will show the following bound:
E[N X = x] wp(r(x))
(77)
Let II be the (1, PoiPossonn process realized by Algorithm
B and let
M = min(m E N (Tm, II, r(Xm) 2 r(z)},
(78)
Le. the index of the first arrival in 1I such that the density ratio of the spatial coordinate exceeds r(z).
Since X; ~ P are i..dd.
M Geom (Pzwe(r(Z)) (r())),
(79)
hence
E[M] Pz~ppr 2 r(z)] wr(r(x))
(80)
Now, since geometrically distributed random variables are memoryless, we find
E[M] = E[M - N M N]
(81)
S E M > N).
(82)
which then implies
E[M M S N < E[M]
(83)
Before we proceed, we first require the following result
Lemma B.1. Let all quantities be as defined above and let 2; y E Q such that r(x) I r(y). Then,
N x = y stochastically dominates N = x, i.e. for all positive integers n we have
P(N 2 n X = x] < P(N 2 n x= y).
(84)
Proof. The proof follows via direct calculation. Concretely,
PIN < X = z] = PIN 1 15 1 X x7]
(85)
E Ej PrIX=
P[N - 1 n- x, T
(86)
(n 1)1 B Bi
[(n, [(m,f-ū((T))
(87)
where [(s, .x) is the upper incomplete r-fucction defined as
del [ 1 dt.
I(s
(88)
jz
18

---

Equation (87) follows from the fact that conditioned on T we have that N - 1 is a Poisson random variable with mean T j(()). Now, focussing on the expectation term in Equation (87 we let
(o(()))
=(r(x)) def E Bi r(n, (oor--(()))
(o- (t) dt. (89)
r(x)
Substituting h r(z), we show that is decreasing function, which proves the lemma. To this end, regardless of the image of the density ratio r(12), extend the definition of e to the entirety of (0, ess sup r(:)). Now, observe, that
=1 by inverse function theorem
de 1 1/ h-TTm,,(h)- (M)())) dh h?
(o(h) -S000 110000/000
1o (o(h)) d'(h) -
(o-1) (t) dt
(90)
Focussing on the integral, via integration by parts we find
(b)
T(n,t-a(t)) (o (1)
(91)
h- h- - ((((h)))) + 1 (t- (e- -1---000 wp(o (t)))o w()o= (t) dt.
poth)
(92)
Plugging this back into the above calculation, we get de J
(h))
(t- (00
(t -
(93)
dh
Note, that since the integrand is positive and h arbitrary, de is negative for every h in its domain, which proves the lemma.
di
Returning Equation 183 we now examine the cases M < N and M N separately.
Case M N: In this case, we have the following inequalities:
r(x) S r(XM) by the definition of M. Applying a to both sides we find (px)) S (p(XM) (Xar) < TM, since M < N and by definition N is the index the first arrival of under the graph of yp.
Now, note that since M < N conditioned on Tas, the first M points of IT are all rejected and hence are all in IT II \ ПI. Thus, M (M N. TM) is Poisson distributed with mean Th (TTM). Hence,
E M NJ B[M - M NJ
(94)
+ ETU~PaA MeN [E[M 1 1 M < N,TMl
(95)
+ ET MeN [MM - [T(TM)]
(96)
2 + BTM~Pry MeN (x(x)) ()(((x)))
(97)
+ ((()) ()(p(z)).
(98)
On the other hand, since the N-th point is the index of first arrival of TI under the graph of yp hence T s xx()). Furthermore, conditioning X
get T ppx)). Similarly to the above N is Poisson distributed with mean T P(T)). Hence, we get
E[N X = x] = 1 + E x -
EIN 118 = xx, 21
(99)
[ T --(]
= + Ejji=
(100)
2 1 + Eyjx== (4((x) - -(p(()))
(101)
= + (xx) - ()(((x))
(102)
19

---

Putting the above two inequalities, we find
E[N
X = x| S ((z) ((v(()) < B[M M N].
(103)
Case M N: In this case, we have
B[M N
M] E[N r(Xx)= 2 r(x)]
(104)
> E XN = X m]
(105)
where Equation (105) follows by applying Lemma B.1
Finally, putting these two cases together we find
1 W (r())))
E[M]
(106)
E N > M]
(107)
z E X z||,
(108)
which establishes Equation 77
Now, taking expectation over X, we get
dQ(x))
No
(109)
2
wr(r(x))
r(x)
dP(x)
(110)
wp(((x))a
y
8(y (((()))) dy
(111)
wp(y)"
y
8(y r())) dP (z) dy
(112)
y . (y((y)
1- wp(y) 1-a dy
(113)
wP(y)* 1- dy.
(114)
In the above, Equation 113 follows by integration by parts and from the fact that
/ 1(r(x) 2 (P:)(x))
wr(h)
(115)
5 h / 8(m r(x)) dP(x) dn.
1 L/a
(116)
We now note that we can bound up using Markov's inequality:
wr(h) = P Pz~rlp (r( 2 h]
(117)
= PPap
TT 2 -I
(118)
S shTi- i2. Ez~p 1r()) 12
(119)
= -h -h Ez~Q (((2)===
(120)
= h - exp2
a
D 1 [Q1PP]]
(121)
1-a
20

---

Putting the above in Equation
114
we get
1
B[Nd]
-
wp (y) 1 dy
(122)
a
1
exxp (Q11)))
+
wp(y) ( dp
(123)
a
1
(xxp((==1111))
dy expa (AD (2D==(Q1P))
14
(124)
1-a
1
1
Doo [Q11
exp2
(aD (
(125)
1-a
log2
B.4 The Codelength of the Index
As discussed in Section 3 Alice and Bob can realize one-shot channel simulation protocol using GPRS for a pair of correlated random variables x.y Px,y if they have access to shared randomness and can simulate samples from Px. In particular, after receiving y ~ Py. Alice runs GPRS with proposal distribution Px and target Pxly. and the shared randomness to simulate the Poisson process II. She then encodes the index N of the first arrival of TI under the graph of 4p. Bob can decode Alice's sample X ~ Pxly by simulating the same N samples from IT using the shared randomness. The question is. how efficiently can Alice encode N? We answer this question by following the approach Li & El Gamal (2018). Namely, we first bound the conditional expectation E|log2 N y = vl. after which we average over y to bound Bllog, N]. Then, we use the maximum entropy distribution subject to the constraint of fixed Ellog, N] to bound H[[N]. Finally, noting that x is function of II and N, get
H[x, N | II] = H[x N, II
HI [N
11]
(126)
H
x, II
+H[[ II]
(127)
=0
from which
H|[x M] = H(N I s H[N].
(128)
which will finish the proof.
Bound on the conditional expectation: Fix y and set Q Pxly and P = Px as GPRS's target and proposal distribution, respectively. Let (t,z) E R+ x l be a point under the graph of p. i.c. (t) S <zz). Then,
r(x) 2 (t)
(129)
[ 2 r] dr
(130)
[P/ 2 T|
2 t- inf
(131)
0.4
P(T 2 t].
(132)
From this, we get
r(x) +1 r(x) +1 IP[T 2
1 1+15
(133)
P[T 2 f|
Next, conditioning on the first arrival (T, we get
Ellog, N T = t, x = œ²| < log E[N - 1 ī =, -
(134)
3 60
log2 (t - (t) + +1)
(135)
< log. (t + 1)
(136)
eq. 183
log, (r(x:) + 1) - log, P(T 2 1
(137)
= loga (r(x) + 1) + f (t) log2 e,
(138)
21

---

where the first inequality follows by Jensen's 's inequality. Now, by the law of iterated expectations, find
E [log2 NN-
EE x
E 1X
Ellog2
N T, [X.]]]
(139)
0q
x
logg
rr x
+ +1]
+ Eup
j T
loga
(140)
1 by eq. (65)
Dku[Q|PP|+E,
log2
1 r
+ log2 e
(141)
{
}]
- DkI[Q||P]
loge
e (x(
logg e + logg
(142)
D P| + logg e,
(143)
where the second inequality follows from switching the natural base and using first-order Taylor approximation for e..
Bound on the marginal expectation and entropy: Equation (143) is a one-shot bound, which yields
Elogg N yl S PKL[Pxt] [Px] + 2- loga c.
(144)
Taking expectation over y Py, we get
Elogg N] < y] + 2- logg
(145)
Finally, following Proposition 4 in Li & El Gamal (018), the maximum entropy distribution for N
subject constraint on Ellog, N obeys
H[N]
Ellog2 N] log (Bllog, N] + 1) 1.
(146)
Plugging in Equation
145 into the above, we get
H[[N] < [x; y] + logg [[[ + log, e + 1) + log, e +
(147)
< 7 [xx y] + loga [X] + 1) + 2 loga + 1 + loga loga + 1)
(148)
[X: y] + log, y[ + 1) 5.84
(149)
y] + loga [X: + 1) 6.
(150)
Similarly to Li & El Gamal 2018). can encode N using Zeta distribution <(n 8) n with
def
(151)
y| + 2 2ogg
With this choice of the coding distribution, the expected codelength of N is upper bounded by [| (/ + 1) + 7 bits.
C Analysis of Parallel GPRS
Parallel GPRS (PGPRS) is general technique to accelerate GPRS when parallel computing power is available and is presented in Section 3.1 Assuming that J parallel threads are available, for target distribution Q and proposal P. PGPRS simulates // Poisson processes II,,... .i....,II,, II, in parallel, all of them with mean measure μ4 def } P. Assuming (TT). x
(7()) x are the first arrivals of II...., II, under i, respectively, PGPRS selects the arrival with the overall smallest arrival time J = arg min
( N, ) By the superposition theorem (Kingman
(J)) N₂ 1992). UI = II is Poisson process with mean measure μ À x P. hence (T X will be the first arrival of under the graph of yp and the theoretic construction in Appendix A therefore guarantees the correctness of PGPRS.
22

---

Expected runtime: We will proceed similarly to the analysis in Appendix B.2 Concretely, let T TV be the first arrival time of II under the graph of yp. Let
Vj def { (T. 1T < 7, (x))
(152)
be the number of points in I1, that are rejected before the global first arrival occurs. By an indepen- dence argument analogous to the one given in Appendix [.2 the V2,..., My are independent Poisson random variables, each distributed with mean
- 1 (f-a (T))
R
(153)
Hence, by the law of iterated expectations, we find that we get
62
(154)
rejections in each of the J threads on average before the global first arrival occurs. Now, thread terminates either when the global first arrival occurs in it or when the thread's s current arrival time T!/) provably exceeds the global first arrival time. Thus, on average, each thread will have one more rejection compared to Equation [154). Hence the average runtime of a thread
Bly 1| -1
J +1,
(155)
and the average number of samples simulated by PGPRS across all its threads is
J ( +1 --1
(156)
Variance of Runtime: Once again, similarly to Appendix B.2 we find by the law of iterated variances that the variance of the runtime in each of the j threads is
Vi, + 1]
Viy
J J +
(v
(157)
meaning that make (()/J) reduction in the variance of the runtime compared to regular GPRS.
Codelength: PGPRS can also realize a one-shot channel simulation protocol for pair of correlated random variables x, y Px.y. For a fixed y ~ Py. Alice applies PGPRS to the target Q = Fxly and proposal P Px. and encodes the two-part code (J+, N. Bob can then simulate Nr from y and recover Alice's ' sample.
J- samples
Encoding J": By the symmetry of the setup, the global first arrival will occur with equal probability in each subprocess īT Hence J follows a uniform distribution on [1,....J]. Therefore, Alice can encode J optimally using [log2 /| bits.
Encoding N We can develop a bound using an almost identical argument to the one in Ap- pendix B. 4 In particular, by adapting the conditional bound in Equation (188) appropriately using
Equation 155). we get
E logg No i = =x= ] (rg((((2) + 1) + p(t) logg logg J.
(158)
Then, using this conditional bound and adapting Equation 143) we find E [logg N, yy, J'] < Dkr. | Px]
2 logg e - log2 J
(159)
to obtain one-shot bound. Taking expectation over y Py. we get [logz Nr JJ] S l[x;y] + 2 log, e logg J, hence, by adapting the maximum entropy bound in Equation (150). we find
(160)
H- J*] < y[x log, J + logg( (1822[[[x [[xx y Thus, we finally find that the entropy of the two-part code (J* N, is upper bounded by
log, J 1) + 6.
(161)
H < I [x; + log, yx logz J + 1) +7.
(162)
Using Zeta distribution ((n s) x n - to encode N J with
def
(163)
1x: yj 2 log, e log, J we find that the expected codelength of the two-part code is upper bounded by y] + log, [x] log2 J + 1) + 8 bits.
(164)
23

---

D
Simulating Poisson Processes Using Tree-Structured Partitions
Algorithm 7: Simulating tree construction with another
Input
Prrpossl distribution
Target splitting function split larget>
Simulating splitting function split Target heap index Harget
Output: : (T, X)) with heap index Hurget in Tearget and heap index Hsim of the arrival in
Tim
P PriorityQueue
B + Ω
K + [log2 Husgel
X P
T Exp(1)
P. PPpuhhTT X, $2, 1)
for k ..... K do
repeat
T, T,X, C, + Popp()
L + splitsim sim (C)
BnL + then
A.. ~ EExp(((L)) Ti T X Ple P.push(TT, X1, L, 2H) end
if BBBR 0 then AR Ex)()(R) Tn + AnAn Xn PlR P.push (Tn. Nn. RH + 1)
end
until X B
/* Exit loop when we find first arrival in B. */
if h: < K then
[Ba, B1} < spli targer (B)
d +-
Huga 2K-E
mod? 2
B t Ba
end
end
return 7,X, H
In this section, we examine an advanced simulation technique for Poisson processes, which required to formulate the binary search-based variants of GPRS. We first recapitulate the tree-based simulation technique from (Flamich et al., (022)) and some important results from their work. Then, finally formulate the optimally efficient variant of Note: For simplicity, we present the ideas for Poisson processes whose
present Algorithm 7 using which
GPRS in Appendix spatial measure P is non-atomic. These ideas can also be extended to atomic spatial measures with appropriate modifications.
Let II be spatiotemporal Poisson process on some space S with non-atomic spatial measure P. In this paragraph, we will deal with itself yet, but the space on which is defined and the measure P. Now, assume that there is function split, which for any given Borel set St produces (possibly random) P-essential purtition of B consisting of two Borel sets, i.e.
Splitting functions:
split(B) = (L, R}, such that L n R and R) = P(B).
(165)
The condition simply allows split to discard some points from that will not included in either the left right split. For example, if we wish to design splitting function for subset of

---

B _ Rd to split B along some hyperplane, we do not need to include the points on the hyperplane in either set. The splitting function that always exists is the trivial splitting function, which just returns the original set and the empty set:
split = 03.
(166)
A more interesting example when f2 = R is the on-sample splitting function, where for a R-valued random variable X ~ Pla it splits B as
split sump (B)
B).
(167)
This function is used by Algorithm 5 and AS* coding (Flamich et al. 2022). dyadic splitting function operating on bounded interval (a, b), splitting as
Another example is the
{( a. 2 a+b) a ++ (1t² ()
split syad ((a, b)) =
(168)
We call this dyadic because when we apply this splitting function to (0,1) and subsets pro- duced by applying it, we get the set of all intervals with dyadic endpoints on (0, 1), i.c. {(0. 1), (0, 1/2), (1/2, 1), (0, 1/4)....... This splitting function is used by Algorithm 6 and AD* coding. As one might imagine, there might be many more possible splitting functions than the two examples we give above, all of which might be more less useful in practice.
Every splitting function on a space Q induces an infinite set of subsets by repeatedly applying the splitting function to the splits it produces. These sets can be organised into an infinite binary space partitioning free (BSP-tree), where each node in the tree is represented by set produced by split and an unique index. Concretely, let the root of the tree be represented by the whole space 1 and the index Hoo = 1. Now we recursively construct the rest of the tree as follows: Let (B, H) be node in the tree, with B Borel set and H its index, and let {L, R} = split(B) be the left and right splits of B. Then, we set (L, 211) as the node's 's left child and (R, 2H 1) its right child. We refer to each node's 's associated index their heap index.
Split-induced binary space partitioning tree:
Heap-indexing the points in II and strict heap invariant: As we saw in Section 2.1 each point in IT can be uniquely identified by their time index N, i.e. if we time-order the points of II, N represents the Nth arrival in the ordered list. However, we can also uniquely index each point in II using a splicting imnction-induce BSP-tree as follows.
We extend each node in the BSP-tree with a point from 11, such that the extended tree satisfies strict heap invariant: First, we extend the root node (S2, 1) by adjoining II's first arrival (Ti, X1) and the first arrival index 1 to get (Ti, X2, Z, 1, 1). Then, for every other extended node (T, X, B, H, N) with parent node (T', X', B', H/2), M) we require that (T,X) is the first arrival in the restricted process II ((T', oo) x B). This restriction enforces that we always have that T T and, therefore, M < N. Furthermore, it is strict because there are no other points of II in B between those two arrival times.
Notation for the tree structure: Let us denote the extended BSP tree T on II induced by split. Each node
T is 5-tuple coordinate X, its bounds B, split-induced heap index H and time index N. As slight overload of notation, for a node v. we use as subscript to refer its elements, e.g. T, is the arrival time associated with node v. We now prove the following important result, which ensures that we do not lose any information by using T to represent II. An essentially equivalent statement was first proven for Gumbel processes (Appendix, "Equivalence under partition" subsection; Maddison et al. 2014).
(T, X, B, H, N) consisting of the arrival time T, spatial
Lemma D.1. Let I1, P, split and T be as defined above. Then, P-almost surely, T contains every point of II.
Proof. We give an inductive argument.
Base case: the first arrival of is by definition contained in the root node of T. Inductive hypothesis: assume that the first N arrivals of II are all contained in T.
Case for N ++1, 1: We begin by observing that if the first N nodes are contained in T. they must form connected subtree Ty. To see this, assume the contrary, i.e. that the first N arrivals form a subgraph I' C T with multiple disconnected components. Let e I' be node contained in component
25

---

of N' that is disconnected from the component of N'v containing the root node. Since T is tree, there is unique path x from the root node in T, and since v and the root are disconnected in I' must contain node E N. However, since is an ancestor of v, by the heap invariant of ' we must have that the time index of e No < N, I N hence T N. a contradiction.
Thus, let TN represent the subtree of T containing the first N arrivals of I1. Now, let FN represent the frontier of TN, i.e. the leaf nodes' children:
FN def
(vE v & TN, parent(v) E Tvh.
(169)
where parent retrieves the parent of node in TT. Let
ΩN def U B
(170)
VEFN
be the union of all the bounds of the nodes in the frontier. A simple inductive argument shows that for all N the nodes in FN provide P-essential partition of f2 from which P(fx) 1. Let Ty be the Nth arrival time of II. Now, by definition, the arrival time associated with every node in the frontier FN must be later than Ty- Finally, consider the first arrival time across the nodes in the frontier:
def
TV
min Tw.
(171)
VEIN
Then, conditioned on TN. TN is the first arrival of 11 restricted to (I', oo) x Q thus it is P-almost surely the N + 1st arrival 11, as desired.
connection between the time index an heap index of a node: Now that we have two ways of uniquely identifying each point in II it is natural to ask whether there is any relation between them. In the next paragraph we adapt an argument from Flamich et al. (2022) to show that under certain cireumstances, the answer is yes.
First, we need to define two concepts, the depth a node in an extended BSP-tree and contractive splitting function. Thus. let T be extended BSP-tree over II induced by split. Let E T. Then, the depth D of / in T is defined as the distance from V to the root. A simple argument shows that the heap index H, of every node obtain the depth of a node y from the heap index via the formula Dj = llgg, H, Next, we say that split is contractive if all the bounds it produces shrink on average. More formally, let 1, and for a node v E T, let Ar denote the set of ancestor nodes of v. Then, split is contractive if for every non-root node E T we have
depth D in T is between 2D < H < 2D+ Thus, we can easily
TT D. [P(b,)]
(172)
where TA, and denotes the subtree of T containing the arrivals of the ancestors of v.
Note that 2 1/2 for any split function. This is because if {L,, R) = split(B) for some set B, then PR P(B)= - P(L). Thus, if P(R)- = (B), then P(L) = (1 P and by definition max a). which is minimized when 1/2, from which = 1/2.
For example, spli ty is contractive with 1/2, while split not contractive. By Lemma 1 from Flamich et al. (2022). split also contractive with e == 3/4. We now strengthen this result using a simple argument, and show that split on-samp is contractive with e = 1/2.
Lemma D.2. Let v. D be defined as above, let P be non-atomic probability measure over R with CDF Fp. Let II (1, P.. Poisson process and T be the BSP- over induced by split sampppp Then,
ETA, (PBB)) -
(173)
Proof. Fix D. = d, and let Na = (n E Dn =d3 be the set of nodes in T whose depth is d.
Note, that [Nal = 2d. We will show that
Vn, m E Na:
P(B,) 1 P(Bm).
(174)
26

---

i.c. that the distributions of the bound sizes are all the same. From this, we will immediately have B[)( [[,( for every n E Na. Then, using the fact, that the nodes N2 for P-almost partition of f2, we get:
B.) = E
-E
1 -
P
(UBB) B.
U
P(B.)
Σ E[P(B.)] Wal [[P(B.]] (175) nEN₂
Lne
Dividing the very left and very right by Wal = 2d yields the desired result.
To complete the proof, we now show that by symmetry, Equation 174 holds. We begin by exposing the fundamental symmetry of on for a node v with left child L and right child Rt the left
and right bound sizes are equal in distribution:
P(B - P(B.) P(B.).
(176)
First, note that by definition, all involved bounds will be intervals. Namely, assume that B, (a. b) for some b and X. is the sample associated with Then, B,, (a, X²) and BR (X,, b) and
hence P(Br) = Fp(X,) (X FF((a). Since X, Pl by the probability integral transform, we have
F(X,) ~ Unif(Fp(a)) Fp())), from which P(Br) Uni (b) (b)=)((a)) Unif P(B))). Since P(BR) = P(B)) - P(B)), we similarly have P(BR) Unif P(B,)), which establishes our claim.
Now, to show Equation (74). fix and fix n E Na. Let An denote the ancestor nodes of n. As we saw in the paragraph above,
P(B,) P(Bparent(n)) d P(Bparent(()) - U, U ~ Un, 1),
(177)
regardless of whether n is a left or right child of its parent. We can apply this d times the ancestors of n. find the marginal:
P(B,) 4 II U, U.~ Unif(0, 1).
(178)
Since the choice of n was arbitrary, all nodes in Na have this distribution, which is what wanted to show.
Now, we have the following result.
Lemma D.3. Let split be a contractive splitting function for some c E [1/2, 1). Then, for every node in T with time index N, and depth Dy, have
ED,IN, [Do] - log, Ny.
(179)
Proof. Let examine the case N, first. In this case, is the root of T and has depth D, by definition. Thus, < - log, holds trivially.
Now, fix E T with time index N, = N > 1. Let TN-1 be the subtree of T containing the first 1 arrivals, FN-1 be the frontier of TN and TN-1 the (N - (1)st arrival time. Then, as we saw in Lemma D.I the Nith arrival in 11 occurs in one of the nodes in the frontier FN-1, after T' -1- In particular, conditioned on TN- the arrival times associated with each node f E FN-1 will be shifted exponentials T Tx-1 + Exp(P(By)) and the Nth arrival time in II is the minimum of these: T, Tn = minsef=-1 Ty. It a standard fact (see e.g. Lemma Maddison (016)) that the index of the minimum
N N-1
Fx arg min JEFN
(180)
is independent of T. = Tey and P[FN that the number of nodes on the frontier FN-11 = N. Thus, we have a simple upper bound on the entropy of F:
f Tv = P(By). A simple inductive argument shows
B l- log, P(B,) = H|[Fy IN = =N] S log, N.
(181)
27

---

Thus, taking expectation over TN-1. we find
[8]
log2 N
E Fx.Tx TN 1|N, -IN,-NI- log2 P(Bv)]
(182)
En.IN, BFN. 11D,,N, - x- log2 P(B,)]
(183)
ED.IN [- logg BF F (P(L.. PPB,)]
(184)
172 0q.
2
E.,.NN. log2 P.]
(185)
logg E) IN.
(186)
The second inequality holds by Jensen's inequality. In the third inequality, apply Equation and one might worry about conditioning on N, here. However, this is not an issue because
172
E
[P(
(187)
1/d s 1
Σ fEFN-1
P(FN fID,
(188)
1[d - N 1]
E
é P.PNFN = fl =d
(189)
SEFN-1
= 1|d3 - -c
(190)
Thus, we finally get
(-logz e) ED.IN,= [D,] S log, N
ED.IN [Do] S llgg, No
(191)
by dividing both sides by - -log2* and we obtain the desired result.
Converting between different heap indices: Assume now that we have two splitting functions, spl ielit,, jarget and splitsim which induce their own BSP-ordering II, Tuaret and Tsim Now, given spl it taggg---gdde heap index H, Algorithm 7 presents method for simulating the appropriate node E Tagge by simulating nodes from Tim In other words, given node with some heap index induced by a splitting function, Algorithm 7 lets find the heap index of the same arrival induced by different splitting function. The significance of Algorithm 7 is that it lets us develop convenient search methods using given splitting function, but might be more efficient to encode the heap index induced by another splitting function.
Theorem D.4. Let II, splittagger split,im Turger and Tsim be as above, Let E Tranzet with and let (T Xsim: Himm) be the arrival and its heap index output by Algorithm given the above as input as well as Hy as the target index. Then, Algorithm| i correct, in the sense that
T, Tsim and X d X
(192)
and Haim is the heap index of (Taims Xmm) in Tsim-
Proof. First, observe that when split ungei = split sim Algorithm AlgoaithnnTcollapsss onto just the extended BSP tree construction for II and simply returns the arrival with the given heap index Hurget in T Tuaget- In particular, the inner loop will always exit after one iteration, and every time one and only one of the conditional blocks will be executed. In other words, in this case, the algorithm becomes equivalent to Algorithm 2 in Flamich et al. (2022).
Let us now consider the case when split * split,im Denote the depth of the required node by D2 [log2 J. Now, we give an inductive argument for correctness.
Base case: Consider D, = 0. In this case, the target bounds = Ω2 and the first sample we draw X ~ P is guaranteed fall in Ω. Hence, for D, 0 the outer loop only runs for one iteration. Furthermore, during that iteration, the inner loop will also exit after one iteration, and Algorithm ? returns the sample (T, X) sampled before the outer loop with heap index 1. Since T ~ Exp(1) and X P, this will be a correctly distributed output with the appropriate heap index.
Inductive hypothesis: Assume Algorithm 7 is correct heap indices with depths up to D, d.
28

---

Case D, = 1: Let E Turger be the parent node of v with arrival XppXXp) Then, D, = d. hence by the inductive hypothesis, Algorithm 7 will correctly simulate a branch Turgen up to node p. At the end of the dth iteration of the outer loop Algorithm 7 sets the target bounds B <- B.. Then, in the final, d + 1st iteration, the inner loop simply realizes Tsim and accepts the first sample after X that falls inside B, whose time index T Tp- Due to the priority queue, the loop simulates the nodes of Tam in time order; hence the accepted sample will also be the one with the earliest arrival time. Furthermore, Agooithmm7onl ever considers nodes of Tsim whose bounds intersect the target bounds By hence the inner loop is guaranteed to terminate, which finishes the proof.
E GPRS with Binary Search
We now utilise the machinery we developed in Appendix D to analyze Algorithm 5
Correcntess of Algorithm 5 Observe that Algorithm 3 constructs the extended BSP tree for the on-sample splitting function. Thus, we will now focus on performing a binary tree search using the extended BSP tree induced by the on-sample splitting function, which we denote by T. The first step of the algorithm matches GPRS's first step (Algorithm Hence it is correct for the first step. Now consider the algorithm in arbitrary step k before termination, where the candidate sample (T, X) is rejected, i.e. T > (xX). By assumption, the density ratio r is unimodal, and since is monotonically increasing, 1p Assume for now that X < x*, , the case X > x² follows
r is unimodal too. Thus, let x* E Ω be such that
r(x") = =r", where r" by symmetric argument. By the unimodality assumption, since X < x=, it must hold that for all y < X. we have p(y) < (xX). Consider now the arrival (T2, X2) of II in the current node's 's left child. Then, we will have T < Tr. and X X by construction. Thus, finally, we get
(Dx
x(X)) ((X) < <T,
(193)
meaning that the current node's left child is also guaranteed to be rejected. This argument can be easily extended to show that any left-descendant of the current node will be rejected, and is sufficient to search its fightddescendants only. By similar argument, when X > x²,, find that it is sufficient only to check the current node's left-descendants. Finally, since both Algorithms [nd simulate Il and search for its first arrival under şp by the construction in Appendix A the sample returned by both algorithms will follow the desired target distribution.
Analyzing the expected runtime and codelength: Since Algorithm simulates a single branch of Ton-sample. its runtime is equal to the depth D of the branch. Furthermore, let N denote the time index of the accepted arrival and H its heap index in Ton-sample. which means D [log, H] ]. It is tempting to use lemmas D.2 and D33 to in our analysis of N, H and D. However, they do not apply, because they assume no restriction
D and N, while in this section, we condition on the fact that N. H
and D identify the arrival in II that was accepted by Algorithm[5] Hence, we first prove a lemma analogous to lemma D.2 in this restricted case, based on Lemma 4 of Flamich al. 20022) Note the additional assumptions we need for our result.
Lemma E.1. Let Q and P be distributions over R with unimodal density ratio r = dQ/ dP,, given to Algorithm als the target and proposal distribution as input, respectively. Assume P has CDF Fp. Let H denote the accepted sample's 's heap index, let D = [log2 HJ, and let B1s denote the bounds associated with the accepted node. Then,
3 D
[[P(BH)
H1 <
(194)
Proof. We prove the claim by induction. For H 1, D = 0 the claim holds trivially, since B₁ R. hence P(B.)= = 1. Assume now that the claim holds for D = [log =]=dd =dj=d- 1 for d > 1, and we prove the statement when D = [loge H| = d. Let Xo..... X X. denote the samples simulated along the branch that leads to the accepted node, i.e. Xo is the sample associated with the root node and X₂ is the accepted sample. By the definition of the on-sample splitting function and law of iterated expectations, we have
E Ex [P(BH) H = [xx (xx.(P(Br)] | H]
when == 1
(195)
Exo.a Exa_1 (X [P(BH) H] H] when d 2 2.
29

---

Below, we only examine the inner expectation in the d 2 2 case; computing the expectation in the 1 case follows mutatis mutandis. First, denote the bounds of the accepted node's parent as B B = (a, for some real numbers < B and define = Fr(a)) b (f) and U = FP(Xa-1). Since Xa-1 X ~ L by the generalized probability integral transform we have U Unif(a, b), where Unif(a (a. b) denotes the uniform distribution on the interval (a, b). The two possible intervals from which A will choose (a, X2-1) and (X2-1,18). whose measures are P((a, Xa-1)) P((Xa-1;18) b U. Then, the worst case, the algorithm always picks the bigger of the two intervals, thus P(BH) S max (U- B - U, from which we obtain the bound
Pr(Xa-1)
Fr(()) U a and similarly
3 3 (b (b-o) a) P(Buraa)) (196)
Ixa- H] S Ey [max(U a, UH 4
Plugging this into Equation 195 we get
3 H S Ex [P ((B11)))) H
E Ex (P(BH)
(197)
< 3
(198)
where the second inequality follows from the inductive hypothesis, which finishes the proof.
Fortunately, lemma does not depend on this condition, so we can modify it mutatis mutandis so that it still applies this latter case as well, and we restate here for completeness:
Lemma F.2. Let N, D, H, Bu be defined above. Then
Bl-log P(Bn)) & log, N
(199)
and
DD] -
(200)
Proof. Equation (199) follows via the same argument as Equation [81) and Equation (200) follows
via the same argument as lemma D.3 combined with lemma E.1
Expected runtime: Combining Equation (200) with Equation 143 yields
DKL[Q| P] + 2 log,
B[D
(201)
(4/3)
which proves Theorem 3.6 In the channel simulation case, taking a pair of correlated random variables x, y ~ Px,y and setting Q t -xly and P t Px, and taking expectation over we get
B[D] < 1[x; yl + loga
(202)
log, 1og,(4/3) (4/3)
Entropy coding the heap index: Let E Ton-sample now denote the node accepted and returned by Algorithm We will encode v using its heap index H. As we saw in Appendix D can be decomposed into two parts: the depth D and the residual R::
H = + +R, 0 R PD
(203)
Then, have
l|ll < H|D, R H( D H[D]
(204)
(205)
Dealing with the second term first, we note that the maximum entropy distribution for random variable Z over the positive integers with the moment constraint [[[ 1/p is the geometric distribution with rate p. Thus, based on Equation 2022.. we set
logg (4/3)
p [[x;y]
(206)
loga
30

---

from which we get
H[D] I H [Geom (p)]
(207)
1 1 P log2(1 - -p)
= log2 p
(208)
P
< logg( (xx y] + log2 e) -log2 log2 (4/3) logz
(209)
< logg (1[x; + 1) logg (2 e) logg logg (4/3) + logg e;
(210)
where Equation 208
follows from the fact that the term 1-9 loga (1 p) is bounded from above by log2 which it attains [x: ->
For the [[R DJ term, we begin by noting that there is particularly natural encoding for given a particular run of the algorithm. Let B(D) B be the bounds associated with the depth-D nodes of Ton-sample D, the peiit-on-samplle BSP tree given D. Then, since B(D) ..... 22 P-essential partition of S2
B form a
encode R by encoding the bounds associated with the accepted
sample, (( Here, we can borrow a technique from arithmetic coding (Rissanen & Langdon 1979; ; encode the largest dyadic interval that is subset of B which can be achieved using logg ) ++1 bits. Note, that Algorithm Jonly ever realises single branch of san-se hence, in this section only, let 1 denote the infinite sequence of samples along this branch, starting with X1 P at the root. Furthermore, from here onwards, let Ba and Ta denote the bounds and arrival time associated with Xa in the sequence, respectively. Then, by the above reasoning, have
H[ D] S Ex Tho [D log, P(Bp)]] +1
(211)
Σ E Ex [P(D= = X1s-1,, Ti: log, P(BB) +1.
(212)
We now develop a bound on P(BB). To this end, note, that by the definition of GPRS,
P ++ XI.A-1 Tisk D 2 k] = P N (TTx Ta) x BB) 0 X TIk D k).
Hence,
P[D k X1:k-1, Tisk - X X1. D
(213)
= IP [N 7h. Bx)
2 1 X Ti:k 1
k]
(214)
E [ k-11 ( Th. B₂)
X Ti:-1, k]
(215)
S E [N (TT Tal Bx)
k]
(216)
= (TB Tk-1),
(217)
where eq. (215) follows from Markov's inequality, and eq. 216 follows from the fact that the mean measure of is just restricted version of II's mean measure. For convenience, we now write A₁ = Ti - Th-1. Then, rearranging the above inequality, we get
log, P(Bx) < - llog P P k X Tik-I, A D k] + log, Ac.
(218)
Substituting this into eq. 212)). get
H R D
(219)
E Ex - P
T1a-1,T7.--1,44/ - log, P = X T 11 Ax]
E 01 A.
P(D k
X1 Ti A log, AA] +
(220)
k=
31

---

We deal with the above two sums separately. In the first sum, we apply Jensen's inequality on the concave function I + -x log, x in each summand, from which we get
Ex T.n P = k X1:k-1, T1x--1, As] logg P = k X1:e-1, T1:-11, Ax]
< P[D = k) log P.D = k]
(221)
k=1 =]
(222)
For the second sum, consider the kth summand first:
P(D XI 1 -1,, Ak, D 2 k] logg Ak
(223)
dPID k. X X1 Ti:k-I, AR, D k]
loga A (a))))
dP
Be
1[r(x) 2 o
(T&)) log2 A dP(x).
P(Bg) BB
To continue, we now examine the integrand:
0q.
r(x) 2 o (Tk)
P Tx]
Tk 2 P N T
As
(224)
from which
log2 A S log2 r(x) log, PIT N Ty]
(225)
Substituting this back, we get
P[D k XI:k-1, 1, Ti:-- AR, D 2 k] logz As
(226)
1
/
$ P(B+) Ba
1[r(x) 2 20 log, r(x) dP(x)
TP[ k X Ti:--1, Ar. D 2 k] log P( 2 Tx]-
(227)
Therefore, for the second summand, we get
Ex -oak-ITha--
[F[D- = k 1
log2 Ax]
(228)
Ex-P
dPP
D
k: Xx X
logg r(x)
Σ P = k] logg P 2 Tkl.
dP
(229)
For the first term, by Fubini's theorem and by the correctness of GPRS, we get
00
dP
k, X2 X]
dP k X2 X
E
X
log, r(x)
= ExaP
log2 r(x)
dP
dP
(230)
Exap [r(X) log, r(X)
(231)
DPL.Q||P..
(232)
For the second term, note that
P = kj = P( = Tk].
(233)
Hence,
P.D kj loga
P( > Th) =
P(D k] log2 P > k]
(234)
k=1
1
(P(D 2 k||
P(D A: + 1))
logg
(235)
P > kl
1=1
< - J,
log, ddp
(236)
loga e.
(237)
32

---

Thus, we finally get that In the average case, this yields
H D) PKL[Q||| TH| [ log2 + 1.
(238)
HI[R DJ s [x + [[] log2 + 1.
(239)
Thus, we finally get that the joint entropy bounded above by
Hi R. DJ < y| + 2 logg [|[x 1)
(240)
2 (log2 logg e) log2 log2 (4/3) log2 e) log2 +
(241)
//x + 2 logg x +) 11,
(242)
which finally proves Theorem 3.7
F General GPRS with Binary Search
We finally present a generalized version spaces and remove the requirement that be unimodal.
aancch-an-boundd GPRS (Algorithm 6 for more general
Decomposing Poisson process into a mixture of processes: Let II be a process over R* x P as before, with spatial measure P. and let Q be target measure, = o and U = {(t,z) E
Rt x Ω t < (o())} as before. Let 11 be II 11 n U restricted under yp and X1) its first arrival. Let L, R] form an P-essential partiton of f2, i.. L n R = 0 and P(L U R) = 1, and let 11, = 1 I x L and 1R 11 g + x R be the restriction of 11 to L and R, respectively. Let Ii = IT nU and i = UI U be the restrictions of the two processes under yp as well. Let (I) and (t((t) be the mean measures of these processes. Thus, the first arrival times in these processes have survival functions
PIT 1 e (t)
(243)
PIT 2 - -ix(t)
(244)
Note, that by the superposition theorem, M U = II, hence the first arrival (Ti. X2) occurs in either II, or lIR. Assume now that we have already searched through the points of II up to time т without finding the first arrival. At this point, we can ask: will the first arrival occur in 11,, given that Ti 2 r? Using Bayes rule, we find
PIX E 2
DRIELITIEN
(245)
PIT 2
More generally, assume that the first arrival of II occurs in some set A C Ω and we know that the first arrival time larger than 7. Then, what is the probability that the first arrival occurs in some set B C A? Similarly to the above, we find
A = P E B, 2 27,8, A
P E B 2
X1 2
(246)
PI 2 2, A]
p r, 8 B]
(247)
P 2 2 A]
Let &1) be the first arrival of II, Then, the crucial observation is that
X T.
(248)
This enables us search for the first arrival of 11 under the graph of &o using an extended BSP tree construction. At each node, if we reject, we draw a Bernoulli random variable b with mean equal to the probability that the first arrival occurs within the bounds associated with the right child node. Then, if b = 1, we continue the search along the right branch. Otherwise, we search along the left branch.
Note, however, that in a restricted process AT the spatial measure no longer integrates to Further- more, target Radon-Nikodym derivative is r(z). 1/x AA)Q(A). This means we need to change
33

---

the graph to some new graph A to ensure that the spatial distribution of the returned sample is still correct. Therefore, for set we define the restricted versions of previous quantities:
(A(I) def [. 1/7 S (a(l)))P(()) di
(249)
wp(h A) def S 1 S <((x)] dp(x)
(250)
Then, via analogous arguments to the ones in Appendix A
we find
dP = x, x]
1[r Al1lt S pa(x))P(TA t]
(251)
d(A P)
dP [ A = x]
= 1 A] / (
PITA 2 de
(252)
dP
PA = ooor.
(253)
Similarly, setting
1sx Al ((A)) and setting x)(((x))) we get
dP
( =Q(A) J PŞTA 2 t]
(254)
(o7)) (T) = P r| = =P -T, E A).
(255)
(256)
From this, again using similar arguments to the ones in Appendix A we find
(T))) ()))) (0( A) 7 (r(o7 A).
(257)
G Necessary Quantities for Implementing GPRS in Practice
Ultimately, given target- prgep-proposal pair (Q, P) with density ratio r, we would want easy-to- evaluate expression for the appropriate stretch function a or -1 to plug directly into our algorithms. Computing requires computing the integral in Equation (5) and finding o requires solving the non-linear ODE in Equation (4, neither of which is usually possible in practice. Hence, we usually resort to computing
numerically by using an ODE solver for Equation (4)
In any case, we need to compute wp and uq. which are analytically tractable in the practically relevant cases. Hence, in this section, we give closed-form expressions for wp and wq for all the examples we consider and give closed-form expressions for a and - for some of them. If we do not give closed-form expression of or we use numerical integration to compute
instead.
Note on notation: In this section only, we denote the natural logarithm as In and exponentiation with respect to the natural base exp.
G.1 Uniform- Uniform Case
Let P be the uniform distribution over an arbitrary space Q and Q uniform distribution supported on some subset X' C Ω, with P(x)) C' for some C' S Then,
r(x) 1 C
1 x]
(258)
wp(h) = C
(259)
(h)h) = 1
(260)
1
(hh)
C log (1 - -C.
(261)
1
- = (1 exp(- ppp C hh))
(262)
C
Note that using GPRS in the uniform-uniform case is somewhat overkill, as this case, it is simply equivalent to standard rejection sampling.
34

---

G.2
Triangular-Unifrra Case
Let = (n 1) and for some numbers <<<<<< < 1, let be the triangular distribution, defined by the PDF
if < a if S <c
2(x-a)
(b (b-a](c-a)
((x))
2
if =
(263)
2(b-x)
(b-o][b-c)
if 1b< x.
In this parameterization, the distribution has support on (a, b). Moreover, it is linearly increasing on (a, e) and linearly decreasing on (c, b), with the mode located at c. For convenience, let & a.
Then.
r(x) ((x))
(264)
2
1
[ log2
loga
(265)
2
[o] log, l
(266)
log2 (1 + a)
Da[QIIP] loga
(267)
wp(h) = l 2
h
(268)
2
" (h) 4
1/
(269)
2h
o(h) 2 3-6 f h 2t
(270)
o (1) 2 t
(271)
Additionally, in this case we find
4
PP[T 2 1]
(272)
(2 1)2
8 l
P( E di)
(273)
(2 + t)
2G
22 12
0
0,
p9(8)
(274)
28s sign(s)
(1- |s|
(2C1
sin 2/s
+
π
2Si (2)²) (²x)) cos 2|
[2|
(275)
(2 20s 4s2e (s(s)= =
25/
for s S 0,
(276)
g
where 29 and M denote the characteristic and the moment-generatie functions of Î, respectively. Furthermore, Si, Ci and are the sine, cosine and exponential integrals, defined as
Si(x) f t sin(t) dt
(277)
Ci(x) 1- cos(1) di
(278)
Bi(r) = t-ell dt.
(279)
35

---

Finally, G is the Meijer-G function (Sections 5.1-3; Bateman & Erdélyi 1953) Interestingly, we can see that this case, while the expectation of 7 is finite, we find that V[T] co. Thus, by Equation 72 we find that V[N] = well, where N is the number of samples simulated by Algorithm
G.3 Finite Discrete-Discrete/ Piecewise Constant Case
Without loss of generality, let Q be discrete distribution over finite set with K elements defined by the probability vector (r1 < 12 <... тк) and p= Unif([ K]) the uniform distribution on K elements. Note that any discrete distribution pair can be rearranged into this setup. Now, embed this distribution into [0, 1] by mapping P P Unif(O, 1) the uniform distribution [0, and mapping Q to Q with density
K
((x) = Σ 1 K 1). K.rr.
(280)
Then, we can draw a sample a from Q by simulating a sample & Q and computing = K [&/K] 1. Hence, we can instead reduce the problem to sampling from piecewise constant distribution. Thus, let us now instead present the more general case arbitrary piecewise
constant distributions over [0, 1], with Q defined by probabilities (91 < < <... < =g) and corre-
sponding piece widths (wi < W2 < wx). We require, that Σk=1 9h = Σm= k= = 1. Then, the density
((x)
Σ j=1
w S x
x
Σ
Wi
(281)
Wk
1=1
Define rk = 9k/uuu Uk. Then,
K (p)h) Σ
1 S rx]] wk
(282)
k=1
K
(h)
1[h < TM]]
(283)
k=
K
Av BBTE-1
(hh) Σ
[ 1
In
(284)
Bx
B& min(rk, h}
k=1
where we defined
K
K
A
Wirr Σ dj
(285)
J=k
j=
K
Bg
Σ
Wi
(286)
j=k
G.4 Diagonal Gaussian-Gausiian Case
Without loss of generality, let Q N(u, o2) and P = N/O, be d-dimensional Gaussian distributions with diagonal covariance, As slight abuse of notation, let A(xx μ,, o², denote the probability density function of a Gaussian random variable with mean and covariance o² evaluated
36

---

at z. Then, when o² < 1, we have
r(z) = 2 N (x m, 11)
(287)
m
(288)
x²
8 o²
(289)
x²
(1 0²)-
Z = N(n 0, (1 02)1)
(290)
(h) P x2 (d, (|mml?)) S -22 + c]
(291)
wh(h) P x² d m 2
2s hn
(292)
C= 2 (2 Z - d (2
(293)
Unfortunately, in this case, it is unlikely that we can solve for the stretch function analytically, so in our experiments, we solved for it numerically using Equation
G.5 One-dimensional Laplace-L.aplace Case
Without loss of generality, let Q = L(u, b) and P= = (0, 1) be Laplace distributions. As a slight abuse of notation, let f(x jr, b) denote the probability density function of Laplace random variable with mean and scale b evaluated at x. Then, for < 1 we have
11|
r(x)
exp
+
|x]|
(294)
1 exp
bla(bk) I-b
cosh [
In I -1 6
wp(h)-
exp
6s In(bh)) |p|
1 1-62
sinh
b(ln(bh)-
1 1-b²
2/2
(bh) F exp
(bh) निm exp (-1%)) In
-Itt
la h b.
I+b
(295)
1 exp
1-b lu(hb) - m1 cosh 1-
In(hb) cosh (
if in S
5-1
Iu|
5
Im
(h)
1 exp
((In(hB)-lul) 1-6"
} (bli) I'M exp 1
(bh)
rie exp (-141
|u|
if In >
Ial
In 6
5
(296)
from which
1 1- o + (br b
-
cosh ( (o-1)
if In S -lui Jul - -lnn b
(a-1)²=
20 [0.1]
(297)
1 レー b
exp
|4| (0-1)
b r
b
exp
(a-1)
+
otherwise
Similarly to the Gaussian case, we resorted to numerical integration using Equation to solve for
H Experimental Details
Comparing Algorithm 3 versus Global-bound A coding: We use setup similar to the one used
by Theis Yosri (2022). Concretely, we assume the following model for correlated random variables
x, JC
P, = 1)
(298)
P Pxin= = (2)
(299)
37

---

From this, we find that the marginal on must be Px = =V(0, o² + 1). The mutual information is
7[x; [i] = log (1 + 2²) bits, which is what we plot as [xx; 14] in the top two panels in Figure 2
For the bottom panel in Figure 2 we fixed standard Gaussian prior
(O 1), fixed K DKL Q P| and linearly increased R = Do [Q] To find a target that achieves the desired values for these given divergences, we set mean and variances as
a² exp (W(A exp(B)) B
(300)
/ 2K o² + on +1
(301)
A = 2ln R 2K -
(302)
B = 2l R 1,
(303)
where W is the principal branch of the Lambert W W-function (Corless et al. 1996).
We can derive this formula by assuming wish find μ and o² such that for fixed numbers K and R, and ((x) = <x μ. 1,,²),,(z) = (xx 0, 1). Then, we have that
DKI.(ale)] K and sup ((x)) = R. p(x)
(304)
We know that
1
K = DKL[gllp]
[nt + - llog 0² --11 -1]
2
(305)
((x))
112
log R log sup zER
llog
p(x))
2(1 02)
From these, get that
1²2 = 2K- - 0² + log 0² +
(306)
" = 2(1 - -0o) 12 + log σ)
Setting these equal to each other
2K - 0² + log 2 + = 2 log R + log 0² 20 log R - o2log o²
2² logo - 0² + 20 log R = 2 log R 2K -
a² log + (2 R - 1) A o² ( 0² + B) A 0² log(o22
A
(307)
B 3 2 2) = A
Ing(a'(0) log 2e Ae
log 2 = W(Ae
0² W(de")
where we made the substitutions A = 2 log R 2K and B = 2og R 1.
I Rejection sampling index entropy lower bound
Assume that we have pair of correlated random variables x, ~ Px.y and Alice and Bob wish to realize channel simulation protocol using standard rejection sampling as given by, e.g. Algorithm 2 Thus, when Alice receives a source symbol y ~ Py. she sets Q = Pxty as the target and P = Px as the proposal for the rejection sampler. Let N denote the index of Alice's 's accepted sample. which is also the number of samples she needs to draw before her algorithm terminates. Since each acceptance decision is an independent Bernoulli trial in standard rejection sampling. N follows geometric distribution whose mean equals the upper bound M used for the density ratio (Maddison 2016). The lower bound on the optimal coding cost for N is given by its entropy
11
1
H N - (M - 1) loga
1
+ log, M 2 log, M,
(308)
M
38

---

where the inequality follows by omitting the first term, which is guaranteed to be positive since I ト 12- log, x is positive on (0, 1). Hence, by using the optimal upper bound on the density ratio M Doc Pl) and plugging it into the formula above, we find that
1x1Q11 Pl & H[N]
(309)
Now, taking expectation over y, we find
By~P, [D PPly [|] S H[N]
(310)
39


Reviews:
[{"review_id": "psdsNDW8iT", "reviewer": "Reviewer_sb37", "summary": "This paper proposes a new algorithm for lossy compression, using ideas from Poisson rejection sampling.\n\nGiven a sample $y \\sim P_y$, Alice wants to communicate the smallest number of bits possible such that Bob can simulate $x \\sim P_{x | y}$, when Alice and Bob have access to the distribution $P_{x, y}$ (and shared randomness to allow simulation).\n\nThe proposed algorithm achieves an almost optimal codelength by rejection sampling -- by introducing a temporal random variable t, Alice can communicate the hitting time $t$ for which $(t, X_t)$ lies under an appropriately defined function $\\phi$. Variations of this algorithm are suggested: one variation is to parallelize using $L$ servers to reduce the runtime by a factor of $L$, and another binary search variant which gives an exponential improvement when the Radon-Nikodym derivative of $dP_{x|y} / dP_{x}$ is bounded.", "questions": "Questions are listed in the weaknesses section.\n\nMy main concern is that this is probably a poor fit for NeurIPS. The theorems are nice, but the impact to a NeurIPS audience is probably extremely limited, given that the algorithms are only feasible for explicit likelihood models and scalar random variables.", "limitations": "There needs to be more discussion about the limitations of this work, especially if the authors are going to mention connections between the proposed algorithm and compression of neural networks. The connections to generative modelling are extremely weak.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- The paper has several nice ideas from information theory, and seems to be the first that can provably achieve runtimes predicted by Li and El Gamal.\n\n- The algorithm is built on A* coding, but is more robust, in the sense that while the runtime of A* coding depends on the Renyi entropy between $P_{x|y}$ and $P{x}$, the runtime of the proposed algorithm is proportional to the KL divergence between them, which may be much smaller and finite even when the Renyi entropy is infinite.\n", "weaknesses": "- The rejection sampler requires access to likelihoods of the conditional distribution and marginal, and hence it's not clear how much this can be generalized.\n\n- The inverse function in Eqn (2) is solved numerically. How would this be computed for distributions with harder CDFs? \n\n- All experiments are toy experiments on Gaussian, Binomial variables etc.\n\n- The experiments are for scalar random variables, and extending to higher dimensions would be extremely non-trivial. The authors state that extending this to multivariate Gaussians is an open problem.\n\n- The paper tries to establish connections with areas like generative modelling by arguing that source compression is useful in VAEs/ latent diffusion models, but this is probably not correct. The works of Theis et al for example use generative models for compression, not the other way around.\n\nMinor:\n- The name GPRS for a communication protocol is probably a poor choice given that it's an existing standard.\n- $\\mathbb{V}$ in Theorem 3.1 is undefined.\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their nice comments and valuable feedback and attempt to address the concerns they raise below.\n\n> The rejection sampler requires access to likelihoods of the conditional distribution and marginal, and hence it's not clear how much this can be generalized.\n\nWe believe the reviewer means that in our definition of GPRS, we require exact knowledge of $dQ/dP$ (i.e. an unnormalized version will not suffice). Furthermore, we need to be able to compute the $w_P$ and $w_Q$ quantities just to numerically evaluate $\\sigma^{-1}$, which makes applying GPRS difficult to more complex problems. These concerns are completely fair, and we agree that GPRS, in its current form, would be a terrible general-purpose sampling algorithm! However, we are interested in applying GPRS to solve channel simulation, which we would use as a core component in a neural data compression pipeline; where the most frequently used channel distributions are uniform, Gaussian and Laplace distributions.\n\nWe will improve our presentation in the camera-ready version in two ways based on the reviewer's comment: \n1. We will more clearly delineate between the relevance of using GPRS as a general-purpose rejection sampler (for which it is terrible) and for channel simulation (for which it yields theoretically optimal and SOTA results)\n2. It is an exciting question if we could GPRS extended to give a general-purpose sampling algorithm by lifting the requirements on the exact knowledge of $dQ/dP$ and being able to compute $w_P$ and $w_Q$ analytically.\n\n> The inverse function in Eqn (2) is solved numerically. How would this be computed for distributions with harder CDFs?\n\nWe believe here the author is asking about cases where even if we know the derivative exactly and can compute $w_P$ and $w_Q$ analytically, evaluating them might still be computationally expensive. This is a very good question; sadly, we lack a good general answer. Thankfully, this is not an issue in the practically relevant cases (i.e. uniform, Gaussian, Laplace); these quantities are all cheap to evaluate (see Appendix G for the precise analytic formulae). \n\nThus, based on the reviewer's comment, we will discuss the concerns mentioned above regarding the practicality of the numerical integration in the camera-ready version of the paper.\n\n> The experiments are for scalar random variables, and extending to higher dimensions would be extremely non-trivial. The authors state that extending this to multivariate Gaussians is an open problem.\n\nThis is correct, though we'd like to specify that standard GPRS can also be applied to multivariate Gaussians. The open problem is concerned with finding a fast algorithm with $\\mathcal{O}(D_{KL}[Q || P])$ runtime for multivariate Gaussians.\n\n> The paper tries to establish connections with areas like generative modelling by arguing that source compression is useful in VAEs/ latent diffusion models, but this is probably not correct. The works of Theis et al for example use generative models for compression, not the other way around.\n\nWe think this might be a misunderstanding; our work is fully in line with the work of Theis et al. (2022); in particular, we could apply \"standard\" GPRS to obtain a similar diffusion-model-based lossy compression algorithm. Could the reviewer please elaborate on which of our arguments they are referring to and believe to be incorrect?\n\n> The name GPRS for a communication protocol is probably a poor choice given that it's an existing standard.\n\nWe thank the reviewer for pointing this out! Fortunately, we believe that the two concepts are sufficiently far apart (one being a rejection sampler, the other being a cellular communication standard) that the chances of much confusion arising should be minimal.\n\n> $\\mathbb{V}$ in Theorem 3.1 is undefined.\n\nWe thank the reviewer for noting this, $\\mathbb{V}$ simply denotes the variance of a random variable in analogy to the expectation operator $\\mathbb{E}$. We will clarify this in the manuscript. \n\n> My main concern is that this is probably a poor fit for NeurIPS. The theorems are nice, but the impact to a NeurIPS audience is probably extremely limited, given that the algorithms are only feasible for explicit likelihood models and scalar random variables.\n\nWe appreciate the concern of the reviewer. However, we argue that our work is, in fact, relevant to the learned data compression community at NeurIPS, as evidenced by many related papers at NeurIPS, ICML, and ICLR, such as Havasi et al. (2018), Flamich et al. (2020), Agustsson & Theis (2020) and Theis et al. (2022). All of these works use channel simulation/relative entropy coding. Their primary limitation is the exponential scaling of the runtime of their algorithms, which could be improved by using a faster algorithm. Thus, our paper presents a significant step towards making all the abovementioned methods more practical.\n\nHowever, we must communicate this point more clearly in the paper. Hence, in the camera-ready version, we will include a paragraph or two explaining and discussing this to help the reader put our work into its appropriate context. \n\n## References\n\n - Agustsson, E., & Theis, L. (2020). Universally quantized neural compression. In NeurIPS 2020\n - Flamich, G., Havasi, M., & Hernández-Lobato, J. M. (2020). Compressing images by encoding their latent representations with relative entropy coding. In NeurIPS 2020\n - Havasi, M., Peharz, R., & Hernández-Lobato, J. M. (2018). Minimal random code learning: Getting bits back from compressed model parameters. In ICLR 2018.\n - Theis, L. & Yosri, N. (2022) Algorithms for the communication of samples. In ICML 2022.\n - Theis, L., Salimans, T., Hoffman, M. D., & Mentzer, F. (2022). Lossy compression with Gaussian diffusion."}, {"review_id": "PrVIM2rcdD", "reviewer": "Reviewer_HQuL", "summary": "This paper investigates the problem of one-shot channel simulation, which can be used as lossy compression without involving quantization.  \nA new rejection sampling algorithm called greedy Poisson rejection sampling (GPRS) is proposed. Then,  a parallelized and a branch-and-bound variant is proposed. Those algorithms are analyzed regrading both correctness and runtime.  \nToy experiments on one-dimensional problems show that GPRS compares\nfavourably against A* coding, the current state-of-the-art channel simulation protocol.   ", "questions": "L331-334, please explain in more detail about the application scenario of efficient channel simulation protocol for multivariate Gaussians. ", "limitations": "see above.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The proposed GPRS is properly positioned against previous and concurrent works, and contains sufficient novelty.    \nThe proposed method has much better runtime complexity compared with previous SOTA A* coding.  \nThe manuscript is well written and enough background information is provided for general readers.  ", "weaknesses": "Though this manuscript is a theoretical contribution, it would be better to discuss more about promising application scenarios and current gaps regarding both performance and efficiency, which can encourage more future works in this field and facilitate more applications of similar techniques.  ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their valuable comments and the concerns they raise, to which we respond below.\n\n> Though this manuscript is a theoretical contribution, it would be better to discuss more about promising application scenarios and current gaps regarding both performance and efficiency, which can encourage more future works in this field and facilitate more applications of similar techniques.\n\nWe agree with the reviewer that we need to discuss promising applications for computationally efficient channel simulation/relative entropy coding to put our work into a better context. Indeed, channel simulation already has a few exciting non-trivial applications, such as model compression (Havasi et al., 2018), lossy data compression with perfect realism (Blau & Micali, 2018; Theis et al., 2022), and differential privacy (Shah et al., 2020). Our work in this area is particularly relevant as the practicality of all these works is limited precisely by the undesirable exponential scaling of the runtime of the channel simulation algorithms they use.\n\nHence, in the camera-ready version of the paper, we will add a paragraph or two explaining and discussing this to help the reader put our work into its appropriate context. Furthermore, we will add more interesting possible future directions in the paper's final section to encourage future works, as the reviewer suggested.\n\n> L331-334, please explain in more detail about the application scenario of efficient channel simulation protocol for multivariate Gaussians.\n\nOf course! In fact, the example we give below is the primary motivation for our work; for a more complete discussion, please refer to (Flamich et al. 2020; 2022) and Theis et al. (2022). \n\nThe relevance of multivariate Gaussians is their widespread use in machine learning, in particular for variational inference and diffusion models. Thus, let us assume we wish to develop a neural image compression algorithm using a variational autoencoder (VAE) and channel simulation. We choose a standard Gaussian prior and a mean-field Gaussian variational posterior for our VAE. Note that with a powerful-enough network architecture, this is not a very restrictive choice for the distributions, and the VAE will be able to adapt to the distribution of images very well. Now, for simplicity, we can train this VAE with the standard beta-ELBO on an image dataset to fit the network parameters, or we can also incorporate an adversarial loss term if we are interested in compression with realism constraints.\n\nOnce we fit the model, we can compress a new image as follows, assuming that besides sharing common randomness and the VAE's prior $p$, the decoder also knows the VAE's generative network. The encoder receives a new image $x$ and uses the VAE's inference network to obtain the image's latent variational posterior $q(z | x)$. Then, the encoder uses a channel simulation protocol to encode a single multivariate Gaussian sample $z \\sim q(z | x)$ using the VAE's prior $p$ as the proposal distribution. The decoder obtains a lossy reconstruction of the image $x$ by decoding the encoder's sample and pushing it through the VAE's generative network.\n\nThe practicality of this scheme hinges upon the efficiency of the channel simulation protocol. Previous works have all used inefficient, exponentially scaling encoding algorithms and had to resort to small-scale experiments only or develop tricks to make the encoding procedure feasible. In contrast, we could break up encoding a single multivariate sample into a sequence of one-dimensional samples and apply branch-and-bound GPRS dimensionwise. This would be fast but incur at least 1 bit of overhead per dimension, so the codelength would be somewhat suboptimal. Thus, solving the multivariate Gaussian problem would yield a maximally efficient solution that does not have the dimensionwise codelength overhead.\n\n## References\n- Blau, Y., & Michaeli, T. (2018). The perception-distortion tradeoff. In Proceedings of the IEEE conference on computer vision and pattern recognition\n- Flamich, G., Havasi, M., & Hernández-Lobato, J. M. (2020). Compressing images by encoding their latent representations with relative entropy coding. In NeurIPS 2020\n- Flamich, G., Markou, S., & Hernández-Lobato, J. M. (2022). Fast relative entropy coding with A* coding. In ICML 2022\n- Havasi, M., Peharz, R., & Hernández-Lobato, J. M. (2018). Minimal random code learning: Getting bits back from compressed model parameters. In ICLR 2018\n- Shah, A., Chen, W. N., Balle, J., Kairouz, P., & Theis, L. (2022). Optimal compression of locally differentially private mechanisms. In AISTATS 2022\n- Theis, L., Salimans, T., Hoffman, M. D., & Mentzer, F. (2022). Lossy compression with Gaussian diffusion."}, {"review_id": "GNVLZLBIWk", "reviewer": "Reviewer_4NF6", "summary": "This paper focuses on the channel simulation problem, which finds applications in stochastic lossy compression. In contrast to the importance sampling approach A* coding (Flamich et.al. 2022) for channel simulation with Poisson processes, this work adopts a rejection sampling method. They demonstrate that the standard sampling approach may yield suboptimal code lengths. As an alternative, the authors propose a novel greedy rejection sampling algorithm called the Greedy Rejection Sampling (GPRS) algorithm. This algorithm composites the density ratio $r = \\frac{dQ}{dP}$ with an invertible function referred to as the \"stretch\" function $\\sigma$ that utilizes the temporal structure of the Poisson process without changing the target distribution.  Theoretical findings establish the correctness and optimality of the proposed GPRS algorithm in terms of code length. Furthermore, the expected runtime of GPRS matches the runtime of standard rejection sampling.\n\nIn addition to the vanilla GPRS, the authors introduce parallel and Branch-and-bound variants to further enhance the runtime of the algorithm. Specifically, the Branch-and-bound GPRS demonstrates a provable runtime improvement over A* coding when the density ratio exhibits unimodal characteristics. To validate the theoretical results, the authors conduct several toy experiments.\n\n\n\n\n", "questions": "- GPRS is thoroughly compared to A* coding both theoretically and experimentally. However, it is unclear why the baseline standard rejection sampling (Algorithm 2) is not empirically compared in Section 4 to validate the theoretical assertion that GPRS enhances the code length efficiency compared to standard RS.\n- Additionally, while the more recent works dithered quantization methods (DQ) and greedy rejection coding (GRC) are briefly mentioned in the related works section, there is a lack of theoretical or empirical comparisons with GPRS. Although these methods may have distinct formulations and proof techniques, conducting some comparisons would be valuable in justifying the significance of GPRS and its potential advantages over alternative approaches.\n- GPRS needs to numerically solve an ODE for $\\sigma^{-1}$ before the time steps. What about the extra cost in runtime? ", "limitations": "N/A", "rating": 6, "confidence": 1, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "- Mathematical derivations look rigorous. \n- This paper is well-written and easy to follow.\n- Theoretical results seem to be compelling. The authors successfully demonstrate, through rigorous analysis, that the GPRS algorithm outperforms both standard rejection sampling and A* coding in terms of either runtime or code length in some cases.", "weaknesses": "- Some baselines and closely related works are not compared theoretically or empirically (See the first two bullets in Questions).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their nice comments; we address the reviewer's questions below.\n\n> GPRS is thoroughly compared to A* coding both theoretically and experimentally. However, it is unclear why the baseline standard rejection sampling (Algorithm 2) is not empirically compared in Section 4 to validate the theoretical assertion that GPRS enhances the code length efficiency compared to standard RS.\n\nWe omitted standard RS from the codelength comparison because its codelength is provably suboptimal (as we outline in Problem 1 in Section 2.2 and show in Appendix I), and hence including it in Figure 2 would be uninformative. On the other hand, Figure 2 does more than verify our theorems regarding the expected runtime of the different algorithms; it also demonstrates that our upper bound on GPRS's codelength appears quite tight, and empirically it tends to heavily concentrate around its mean and that the mean appears to be quite robust in that it closely lines up with the median.\n\nIn the camera-ready version of the paper, we will add clarifying sentences based on our response above in the experiments section (Section 4).\n\n> Additionally, while the more recent works dithered quantization methods (DQ) and greedy rejection coding (GRC) are briefly mentioned in the related works section, there is a lack of theoretical or empirical comparisons with GPRS. Although these methods may have distinct formulations and proof techniques, conducting some comparisons would be valuable in justifying the significance of GPRS and its potential advantages over alternative approaches.\n\n**Regarding GRC:** We agree with the reviewer that in the current version of the paper, the comparison between GRC and GPRS in terms of high-level theoretical properties and empirical performance is too loose. We decided to omit GRC from Figure 2 because we felt it made the figure too cluttered; however, in the camera-ready version, we will add another section in the appendix with more empirical comparisons and specifically comparisons to different variants of GRC there. Furthermore, we will extend the discussion of GRC in the related works section and compare their theoretical properties more precisely by contrasting the exact bounds on their runtime, for example.\n\n**Regarding DQ:** Unfortunately, we are uncertain how to compare DQ with GPRS best as its formulation is significantly different. In fact, there isn't currently a complexity-theoretic framework that would cover all channel simulation protocols. We are actively working on this as part of our research agenda. In our paper, since we developed a rejection sampler, the number of samples simulated by the algorithm is a natural proxy for the computational time complexity that lends itself to theoretical analysis. We assume all other aspects (such as simulating a single sample from the proposal or evaluating sigma) take $\\mathcal{O}(1)$ time. However, this proxy is inappropriate for DQ, as it is not a rejection or importance sampler and only needs to simulate one sample from the proposal. Furthermore, DQ only applies to cases where both the target and proposal distributions are uniform, and in this case, it enjoys a clear advantage over GPRS.\n\nIn the camera-ready version, we will clarify that we are using the number of samples simulated by GPRS as a proxy for the algorithms' time complexity and assume that all other operations can be carried out in constant time. Based on our answer above, we will also discuss the difficulties of comparing DQ with GPRS.\n\n> GPRS needs to numerically solve an ODE for $\\sigma^{-1}$ before the time steps. What about the extra cost in runtime?\n\nThis is an important practical question, and currently, we lack a good answer to it. As discussed above, in our runtime analysis, we assumed that evaluating $\\sigma$ or its inverse could be done in constant time and disregarded practical implementation details as it was not a hindrance to carrying out our experiments (we used Scipy's `odeint` function to solve the ODE which was sufficiently fast in practice). That said, developing a cheap approximation to $\\sigma$ or its inverse is an important future direction, and we are actively working on it.\n\nThus, in the camera-ready version, we will discuss the practical considerations for the numerical integration and point out that analyzing the complexity of evaluating $\\sigma$ and developing cheap-to-evaluate approximations is an exciting future research direction."}, {"review_id": "M5ip8GSZ2t", "reviewer": "Reviewer_prqL", "summary": "The paper addresses the problem of representing a target distribution using the least possible number of bits. The authors refined the idea of encoding a sample from the target distribution as the first sample from the proposal distribution that passes a rejection sampling condition. The refined approach is shown to achieve runtime optimality.", "questions": "- What are the requirements the target and proposal distributions should fulfil? Why do the distributions need to be 1-dimensional? Is this only for theoretical purposes (i.e. to estimate the runtime bounds)? Or would it be also a practical constraint?\n- How expansive is it to solve the ODE for computing \\sigma? Should this be included in the total runtime?\n- How is the density ratio estimated? Is the method feasible if the target distribution is unavailable?\n- What are the key differences between the proposed algorithm, A*, and the approach described in Maddison 16? Why is rejection sampling supposed to be better than importance sampling? It would be good to include an intuitive explanation of why and when GPRS is faster.\n", "limitations": "The authors do not discuss the limitation of the proposed approach. It is unclear what assumptions are required for the proofs. The theoretical runtime of related methods is not reported.", "rating": 5, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "The idea of using a Poisson process to encode a target distribution is powerful. The strategy is well known but further investigations and attempts to make it more feasible could have an impact in various domains.", "weaknesses": "The authors' technical contribution. could have been outlined more precisely. It is hard to understand why the new strategy is expected to be better than the A* algorithm or naive rejection sampling. The authors should have explained why \"a criterion that does depend on the time variable\" is better.\n\nIn the introduction, the authors say that the distribution is assumed to be unimodal and 1 dimensional. This is a strong assumption and it is unclear why the proposed method requires it. Similar assumptions are made in Flemich 22. The authors may have specified whether these are also required in Maddison 16.\n\nThe authors should have found a way of testing the algorithm on real-world data.\n\nThe text could be more curated and structured.\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their comments and feedback on our work. We would like to begin by clarifying a crucial point that, in our experience, is the most common source of confusion regarding channel simulation. \n> The paper addresses the problem of representing a target distribution using the least possible number of bits.\n\nThis is not quite correct; in channel simulation, our goal is **not** to encode the target distribution so that the decoder can recover it. Please see lines 25-27 and 136-139 in our paper that outline this. As the reviewer notes in the second sentence of their summary, we only seek to encode a sample from the target distribution.\n\nWe now address the reviewer's concerns and questions.\n> The authors' technical contribution. could have been outlined more precisely.\n\nWe outline our technical and empirical contributions on lines 56-67, visually depict all GPRS variants in Figure 1, sketch GPRS's construction in Section 3, and construct and analyse it rigorously in the Appendices. Did the reviewer find these too vague or mathematically imprecise or simply that the writing needed to be improved? We are happy to incorporate any concrete feedback the reviewer has!\n> It is hard to understand why the new strategy is expected to be better than the A* algorithm or naive rejection sampling.\n\nThis is a tricky question that would take more space to answer than we have available for this rebuttal; if the reviewer is interested, we are happy to follow up with more detail!\n\nNote that the expected runtimes (in terms of the number of samples simulated) of GPRS and naive rejection sampling (NRS) are all equal to $2^{D_\\infty[Q || P]}$ and the only thing that differs is the codelength. Here, GPRS performs better than NRS because by only accepting samples below the graph of $\\phi$ and rejecting everything above, GPRS concentrates the acceptance probability of samples in the earliest step possible. We see this best in the discrete case, where we might draw the same sample $X_k = X_{k+1}$ from the proposal twice in a row. In this case, GPRS will always accept $X_k$, whereas NRS might reject $X_k$ but accept $X_{k+1}$.\n\nAt a high level, GPRS's branch-and-bound variant is better than A* coding's because it immediately terminates once it finds the sample it will accept. On the other hand, when A* sampling finds the sample it will eventually accept, it does some \"extra useless work\" by checking more samples due to its non-greedy acceptance criterion.\n\nIn the camera-ready version, we will add an explanation of why GPRS outperforms NRS and A* coding based on our answers above.\n\n> The authors should have found a way of testing the algorithm on real-world data.\n\nWe developed a rejection sampler; what kind of test using real-world data does the reviewer have in mind?\n> What are the requirements the target and proposal distributions should fulfil? Why do the distributions need to be 1-dimensional? Is this only for theoretical purposes (i.e. to estimate the runtime bounds)? Or would it be also a practical constraint?\n\nAs we mentioned in lines 82-84, we can define the standard version of GPRS over very general probability spaces. We only require that $dQ/dP$ be bounded, meaning that GPRS is theoretically as widely applicable as NRS. On the other hand, the branch-and-bound version requires 1-dimensional distributions over $\\mathbb{R}$ because it can be considered a stochastic bisection search, which is naturally defined over $\\mathbb{R}$ and extending it to more general spaces is highly non-trivial.\n> How expansive is it to solve the ODE for computing \\sigma? Should this be included in the total runtime?\n\nOur analysis in the paper assumes that we can evaluate $\\sigma$ in $\\mathcal{O}(1)$ time; we will clarify this in the camera-ready version. Of course, in practice, we need to ensure that we can cheaply do the numerical integration, but we leave this for future work.\n> How is the density ratio estimated? Is the method feasible if the target distribution is unavailable?\n\nGPRS wouldn't work very well as a general-purpose sampling algorithm in its current form. In this paper, we are only interested in applying it to channel simulation, where it is reasonable to assume that the encoder knows the density ratio exactly.\n> What are the key differences between the proposed algorithm, A*, and the approach described in Maddison 16?\n\nThe approach proposed in Maddison 16 _is_ A* sampling, and channel simulation is not discussed. For a comparison of GPRS and A* _coding_, please see lines 207-214 and 291-300 in our paper.  \n> Why is rejection sampling supposed to be better than importance sampling?\n\nThe answer to this question depends on what problem we are applying the sampler to and what the reviewer means by \"better\". Rejection sampling yields an exact sample from the desired target distribution, but its termination time is stochastic. In contrast, importance sampling has a fixed runtime but only produces an approximate sample.\n> It would be good to include an intuitive explanation of why and when GPRS is faster.\n\nSee our response above for a brief discussion regarding the speed comparison of the algorithms. \n> The authors do not discuss the limitation of the proposed approach. It is unclear what assumptions are required for the proofs.\n\nAs noted in the paper, we intentionally left a lot of discussion in the main text informal and only sketched the constructions of the different variants of GPRS in Section 3 to illustrate the important, high-level ideas. In each section, we refer the reader to the appropriate parts of the appendix for the fully rigorous development of each result. Could the reviewer elaborate on what they found lacking and needing improvement, please? \n\n_We will gladly answer any further questions from the reviewer. On the other hand, should we have answered the reviewer's concerns adequately, we kindly ask the reviewer to consider raising their score._"}, {"review_id": "9LaSM59OTY", "reviewer": "Reviewer_tLU5", "summary": "The paper proposes a new relative entropy coding algorithm for compression without quantization. Compression without quantization is an exciting line of research that tries to eliminate training-test mismatches in learned compression by avoiding discrete representations altogether. Consequently, one can losslessly compress data in a continuous latent variable model, such as a VAE, without quantizing the latent representation. The downside of existing algorithms for relative entropy coding is that their runtimes are intolerably slow. The paper makes a big step towards more runtime efficiency by designing new efficient rejection sampling approaches for transmitting samples from a variational posterior, especially in cases where the likelihood ratio between prior and posterior is uni-modal. ", "questions": "* Could the authors discuss in depth the commonalities and differences to a recent paper with a similar title, “adaptive greedy rejection sampling” (https://arxiv.org/pdf/2304.10407.pdf) that the authors already mentioned? \n* The presentation at the beginning of Section 3 is very dense. For example, for an arbitrary P and Q, how is w_p approximated—analytically or based on samples? How does that affect the numerical solvability of Eq. 2?\n* It appears that the stretching function sigma is learned anew for every datum. How does this aspect affect the runtime, and is it respected in the scaling analysis? A discussion would be adequate. \n* Parallelization and branch & bound: it appears that the parallelization affects the cost to *encode* data only, but not to *decode* data. Provided this is correct, could the authors explain why the encoding step is a bigger concern than decoding? In practice (e.g., video streaming), one mostly cares about reducing *decoding* complexity. \n* Could the authors also comment on how they propose to simulate from the Poisson process restricted to subregions? It is unclear why this saves in runtime performance since sampling from a spatially-restricted Poisson process again amounts to rejection sampling. \n\nWriting suggestion:\n* It took me a long time to understand why a temporal poisson process construction was used, where the original communication problem is time-independent. Maybe just say early on that including a temporal poisson process replaces rejection sampling with a deterministic criterion, a trick that was previously proposed for numerical integration (Maddison).", "limitations": "Limitations in terms of scaling and speed are sufficiently addressed. Societal impacts are a less important topic for general-purpose data transmission protocols. ", "rating": 9, "confidence": 5, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "+ Relative entropy coding is a somewhat under-appreciated but up-and-coming field of research at the intersection between information theory and machine learning. There are only a few papers in this domain. (For example, recent advances made compression in diffusion models possible [Theis et al., 2022].) \n+ The paper significantly improves runtime efficiency for these algorithms, paving the way toward its scalable deployment. Compared to runtime complexities exponential in the Renyi divergence, it achieves a runtime *linear* in the KL divergence between the target and proposal distributions. Notably, the involved algorithm is conceptually much simpler (and in practice faster) than relative entropy coding with A* coding, which achieves the same asymptotic scaling. \n+ Beyond its base version, the paper also proposes two computationally much more efficient extensions: parallel GPRS and branch&bound GPRS. Both approaches rely on non-trivial mathematical insights, e.g., additivity of a Poisson process. \n+ Moreover, in one-dimensional channel simulation problems where the density ratio is unimodal, the algorithm displays the theoretically-optimal runtime (unlike existing approaches such as A* sampling)\n+ The paper shows a nice combination of solid empirical results and theory. For example, most claims on the runtime complexities and theoretical guarantees (finite runtime) are proven rigorously. ", "weaknesses": "See questions below. There are no immediate weaknesses evident, but some clarifying questions need to be addressed. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their glowing review of our work; we are delighted that the reviewer shares our excitement for relative entropy coding/channel simulation!\n\nWe answer the reviewer's questions below and will gladly answer any further questions.\n> Could the authors discuss in depth the commonalities and differences to a recent paper with a similar title, \"adaptive greedy rejection sampling\"?\n\nOf course! The \"obvious\" commonality between GPRS and adaptive greedy rejection sampling (AGRS) is that they are both rejection samplers for channel simulation, i.e. they encode exact samples with optimal expected codelength, but their runtime is stochastic. This contrasts with the depth-limited version of A* sampling, an importance sampler whose runtime is fixed but only returns approximate samples. Moreover, GPRS and AGRS are greedy algorithms, though they differ in _how_ they are greedy. As we explain in lines 207-210, we call GPRS greedy because it performs a greedy search over the points of a Poisson process instead of the non-greedy search performed by A* sampling/coding. On the other hand, AGRS is greedy in the sense defined by Harsha et al. (2009): it greedily maximizes the acceptance probability of a proposed sample under a certain correctness constraint; see Section IV in Harsha et al. (2009) for the precise explanation. These two notions of greediness are likely related, and we are actively investigating this link. \n\nFurthermore, Flamich & Theis (2023) also propose a space partitioning method that empirically appears to achieve an exponential speedup similar to GPRS's branch-and-bound variant. However, the authors' space partition differs significantly from ours in several ways: The authors' method is based on dithered quantization (DQ; Agustsson & Theis, 2020) and hybrid coding (Theis & Yosri, 2022), while we based GPRS's branch-and-bound variant on a bisection search method, inspired by the AS* and AD* variants of A* coding. Moreover, our method comes with theoretical guarantees on the correctness and runtime of the algorithm, while Flamich & Theis do not even prove the correctness of their space partitioning variant, let alone provide guarantees on the runtime; all their results are empirical. That said, it is interesting whether their DQ-based space partitioning method could be adapted to be used by GPRS too.\n\n> ... for an arbitrary P and Q, how is w_p approximated—analytically or based on samples? How does that affect the numerical solvability of Eq. 2?\n\nThis is a good question; unfortunately, we lack a good answer for the general case. Appendix G provides analytic solutions for $w_P$ and $w_Q$ for the most relevant practical cases, i.e., the uniform, Gaussian and Laplace distributions. Hence, numerical integration is fine in these cases. Evaluating $w_P$ and $w_Q$ even approximately and studying its effects on the numerical integration in more general cases is an interesting future direction for research. As a first step, we are currently investigating if there is a natural family of distributions (e.g. natural exponential families) for which $w_P$ and $w_Q$ can always be analytically computed.\n\nIn the camera-ready version, we will add a clarifying discussion on this and note the challenges the reviewer highlighted as interesting future research directions.\n> It appears that the stretching function sigma is learned anew for every datum. How does this aspect affect the runtime, and is it respected in the scaling analysis?\n\nWe are not quite sure what the reviewer means by \"datum\". If they mean target distribution, then this is correct; we need to perform the numerical integration during each encoding procedure; unfortunately, there doesn't seem to be a good way to precompute $\\sigma$. In our runtime analysis, we assumed that we could evaluate $\\sigma$ in $\\mathcal{O}(1)$ time, and we used the number of samples simulated by the algorithm as the proxy for the algorithm's runtime, as it is the dominant factor. However, reducing the complexity of the numerical integration is an important practical concern, and we will note this as future work in the camera-ready version.\n> Parallelization and branch & bound: it appears that the parallelization affects the cost to encode data only, but not to decode data. Provided this is correct, could the authors explain why the encoding step is a bigger concern than decoding?\n\nIn practice, we use the step count as the random seed (mixed with a global seed) for the proposed samples. Hence, the decoder only needs to set the step count they received as their PRNG seed and simulate a single sample from the proposal distribution. In the camera-ready version, we will add a section to the appendix discussing the basic implementation details for GPRS.\n> Could the authors also comment on how they propose to simulate from the Poisson process restricted to subregions? It is unclear why this saves in runtime performance since sampling from a spatially-restricted Poisson process again amounts to rejection sampling.\n\nGreat question! The reviewer is correct in the general case; simulating from arbitrarily truncated distributions is computationally hard. However, we only truncate to intervals in the one-dimensional branch-and-bound variant of GPRS. We can use generalised inverse transform sampling in this case: Suppose we have a real-valued random variable $X$ with CDF $F$ and inverse CDF $F^{-1}$. Then, we can simulate $X$ truncated to $[a, b]$ by simulating $U \\sim \\mathrm{Uniform}(F(a), F(b))$ and computing $X = F^{-1}(U)$.\n\n## References\n\n- Agustsson, E., & Theis, L. (2020). Universally quantized neural compression. In NeurIPS 2020\n\n- Flamich, G., & Theis, L. (2023). Adaptive Greedy Rejection Sampling. arXiv preprint\n\n- Harsha, P., Jain, R., McAllester, D., & Radhakrishnan, J. (2009). The Communication Complexity of Correlation. IEEE Transactions on Information Theory\n\n- Theis, L., & Ahmed, N. Y. (2022). Algorithms for the communication of samples. In ICML"}]

Author rebuttal:
We thank all the reviewers for their valuable feedback on our paper, which will help us improve it significantly. We are delighted that all reviewers agree that our contributions are significant, that most reviewers (sb37, HQuL, 4NF6 and tLU5) found our exposition well-written and easy to follow, and that Reviewers prqL and tLU5 highlighted the potential impact of our work.

However, the reviewers also identified several aspects of the paper that need improvement. Thus, we highlight the most common concerns and the most significant changes we propose for the camera-ready version of our paper.

## Delineate GPRS as a sampling algorithm from its use in channel simulation

We presented GPRS in its full generality and did not comment much on its limitations as a general sampling algorithm, detached from its intended use for channel simulation. Hence, Reviewers sb37, prqL and tLU5 have rightly raised concerns regarding cases that often arise in practice when applying a sampling algorithm, e.g. that the Radon-Nikodym derivative is only known up to a constant or that the functions $w_P$ and $w_Q$ could be troublesome to evaluate.

This question is especially interesting because A* coding, based on the A* sampling algorithm, is more generally applicable in practice just as a sampling algorithm.

On the other hand, since we are ultimately motivated by applying our GPRS to neural compression via channel simulation, the set of distributions of practical interest is much smaller, consisting of the uniform, Gaussian and Laplace distributions. This is because, in neural compression, the data distribution (e.g. the distribution of natural images) is usually approximated by a latent representation following a simple distribution (e.g. Gaussian) that is then transformed by powerful neural networks to match the distribution statistics (e.g. in the case of VAEs and diffusion models).

Hence, in the camera-ready version of the paper, we will add clarifying discussion regarding the limitations of applying GPRS as a general-purpose sampling algorithm and note that extending it to more general cases is future work.

## Applications of GPRS
Reviewers sb37, HQuL and prqL expressed concern that the motivation or application of our method (and perhaps channel simulation in general) needs to be communicated more clearly.

Developing a fast channel simulation algorithm is of great practical interest, as there have been several papers at NeurIPS, ICML, and ICLR, such as Havasi et al. (2018), Flamich et al. (2020), Agustsson & Theis (2020) and Theis et al. (2022) utilizing it. All of these works use channel simulation/relative entropy coding. Their primary limitation is the exponential scaling of the runtime of their algorithms, which could be improved by using a faster algorithm. Thus, our paper presents a significant step towards making all the abovementioned methods more practical.

Hence, in the camera-ready version, we will include a paragraph or two explaining and discussing this to help the reader put our work into its appropriate context.

## Concerns regarding the implementation and numerics
Related to the first point, Reviewers sb37,  4NF6, prqL and tLU5 expressed concerns regarding the implementation details and numerics of GPRS, such as:
- how $w_P$ and $w_Q$ can be computed or approximated
- the numerical stability and time complexity of evaluating $\sigma$ or $\sigma^{-1}$
- the intricacies of the encoding and decoding procedures

While there are many subtle details, here we would like to emphasize that:
 - for the practically relevant channel simulation cases $w_P$ and $w_Q$ have analytic forms which are stated in Appendix G
 - for the runtime complexity analysis we assumed that $\sigma$ or its inverse can be evaluated in $\mathcal{O}(1)$ time
 - the time complexity of the numerical integration is essentially independent of the number of samples simulated, and we found it to be negligible in our experiments, as the number of samples GPRS simulates heavily dominates the time complexity anyways

To address all the above and more, in the camera-ready version of the paper, we will add another section to the Appendix to detail all relevant considerations for the practical implementation of each GPRS variant.

__We thank the reviewers again for their time reviewing our work and providing important feedback to improve our work. We will gladly address any further questions the reviewers might have during the discussion phase. Should we have answered a reviewer's concerns adequately, we kindly invite them to consider raising their score.__

## References
- Agustsson, E., & Theis, L. (2020). Universally quantized neural compression. In NeurIPS 2020
- Flamich, G., Havasi, M., & Hernández-Lobato, J. M. (2020). Compressing images by encoding their latent representations with relative entropy coding. In NeurIPS 2020
- Havasi, M., Peharz, R., & Hernández-Lobato, J. M. (2018). Minimal random code learning: Getting bits back from compressed model parameters. In ICLR 2018.
- Theis, L. & Yosri, N. (2022) Algorithms for the communication of samples. In ICML 2022.
- Theis, L., Salimans, T., Hoffman, M. D., & Mentzer, F. (2022). Lossy compression with Gaussian diffusion.


Meta-review: The authors solve the problem of 1-shot channel simulation (a fundamental data compression problem) in the 1-dimensional case where the target-proposal density ratio is unimodal, by giving a provable algorithm based on Poisson processes with optimal runtime. Empirical results show improvement over SOTA algorithms.
All reviewers recommend acceptance after discussion. I recommend acceptance.
Decision: Accept (poster)