#  Perfect Sampling from Pairwise Comparisons

Dimitris Fotakis111National Technical University of Athens,
<EMAIL_ADDRESS>Supported by the HFRI, project BALSAM, HFRI-FM17-1424.
NTUA  Alkis Kalavasis222National Technical University of Athens,
<EMAIL_ADDRESS>Supported by the HFRI, project BALSAM, HFRI-
FM17-1424.
NTUA  Christos Tzamos333University of Wisconsin-Madison & University of
Athens<EMAIL_ADDRESS>Supported by NSF grants CCF-2144298 and CCF-2008006.
UW Madison & UOA

###### Abstract

In this work, we study how to efficiently obtain perfect samples from a
discrete distribution $\mathcal{D}$ given access only to pairwise comparisons
of elements of its support. Specifically, we assume access to samples $(x,S)$,
where $S$ is drawn from a distribution over sets $\mathcal{Q}$ (indicating the
elements being compared), and $x$ is drawn from the conditional distribution
$\mathcal{D}_{S}$ (indicating the winner of the comparison) and aim to output
a clean sample $y$ distributed according to $\mathcal{D}$. We mainly focus on
the case of pairwise comparisons where all sets $S$ have size 2. We design a
Markov chain whose stationary distribution coincides with $\mathcal{D}$ and
give an algorithm to obtain exact samples using the technique of Coupling from
the Past. However, the sample complexity of this algorithm depends on the
structure of the distribution $\mathcal{D}$ and can be even exponential in the
support of $\mathcal{D}$ in many natural scenarios. Our main contribution is
to provide an efficient exact sampling algorithm whose complexity does not
depend on the structure of $\mathcal{D}$. To this end, we give a parametric
Markov chain that mixes significantly faster given a good approximation to the
stationary distribution. We can obtain such an approximation using an
efficient learning from pairwise comparisons algorithm (Shah _et al._ , JMLR
17, 2016). Our technique for speeding up sampling from a Markov chain whose
stationary distribution is approximately known is simple, general and possibly
of independent interest.

##  1 Introduction

Machine Learning questions dealing with pairwise comparisons (PCs) constitute
a fundamental and broad research area with multiple real-world applications
[SBB+16, Kaz11, LOAF12, SBGW16, FH10, WJJ13]. For instance, the preference of
a consumer to choose one product over another constitutes a pairwise
comparison between the two products. The two bedrocks of modern ML
applications are inference/learning and sampling. The area of learning from
pairwise comparisons has a vast amount of work [RA14, SBGW16, DKNS01, NOS12,
Hun04, JLYY11, APA18, LSR21, HOX14, SBB+16, VY16, NOS17, VYZ20]. In this
problem, there exist $n$ alternatives and an unknown weight vector
$\mathcal{D}\in\mathbb{R}^{n}$ where $\mathcal{D}(i)$ is the weight of the
$i$-th element with $\mathcal{D}(i)\geq 0$ and
$\sum_{i\in[n]}\mathcal{D}(i)=1$, i.e., $\mathcal{D}$ is a probability
distribution over $[n]$. Roughly speaking, the learner observes pairwise
comparisons (or more generally $k$-wise comparisons) of the form
$(w,\\{u,v\\})\in[n]\times[n]^{2}$ that compare the alternatives $u\neq v$ and
$w$ equals either $u$ or $v$ indicating the winner of the comparison. The most
fundamental distribution over pairwise comparisons is the BTL model [BT52,
Luc59] where the probability that the learner observes $w=u$ when the items
$u,v$ are compared is $\mathcal{D}(u)/(\mathcal{D}(u)+\mathcal{D}(v))$. The
learner’s goal is mainly to estimate the underlying distribution $\mathcal{D}$
using a small number of pairwise comparisons.

In this work, we focus on questions concerning sampling, the second bedrock of
modern ML. In fact, we will deal with a fundamental area of the sampling
literature, namely _perfect sampling_. In contrast to approximate sampling,
perfect (or exact) sampling from a probabilistic model requires an exact
sample from this model. Problems in this vein of research lie in the
intersection of computer science and probability theory and have been
extensively studied [PW96, PW98, Hub04, Hub16, FGY19, JSS21, AJ21, GJL19,
BC20, HSW21]. Our work poses questions concerning perfect sampling in the area
of pairwise comparisons. We ask the following simple yet non-trivial question:
Is it possible to obtain a perfect sample $y$ distributed as in $\mathcal{D}$
by observing only pairwise comparisons from $\mathcal{D}$?

Apart from the algorithmic and mathematical challenge of generating exact
samples, perfect simulation has also practical motivation, noticed in various
previous works [Hub16, JSS21, AJ21]. Crucially, its value does not stem from
the fact that the output of the sampling distribution is perfect; usually, one
can come up with a Markov chain whose evolution rapidly converges to the
stationary measure. Hence, the distance between the (approximate) sampling
distribution from the desired one can become small efficiently. This
approximate sampling approach has a clear drawback: there is no termination
criterion, i.e., the Markov chain has to be run for sufficiently long time and
this requires an a priori known bound on the chain’s mixing time, which may be
much larger than necessary. On the contrary, perfect sampling algorithms come
with a stopping rule and this termination condition can be attained well ahead
of the worst-case analysis time bounds in practice.

Second, the question that we pose appears to have links with the important
literature of truncated statistics [Gal97, DGTZ18]. Truncation is a common
biasing mechanism under which samples from an underlying distribution
$\mathcal{D}$ over $\mathbb{X}$ are only observed if they fall in a given set
$S\subseteq\mathbb{X}$, i.e., one only sees samples from the conditional
distribution $\mathcal{D}_{S}$. There exists a vast amount of applications
that involve data censoring (that go back to at least [Gal97]) where one does
not have direct sample access to $\mathcal{D}$. Recent work in this literature
[FKT20] asked when it is possible to obtain perfect samples from some specific
truncated discrete distributions (namely, truncated Boolean product
distributions). The area of pairwise comparisons can be seen as an extension
of the truncated setting where the truncation set can change (i.e., each
subset of alternatives corresponds to some known truncation set); while
perfect sampling has been a subject of interest in truncated statistics, the
neighboring area of pairwise comparisons lacks of efficient algorithms for
this important question.

##### Problem Formulation.

Let us define our setting formally (which we call Local Sampling Scheme). This
setting is standard (without this terminology) and can be found e.g., at
[RA14].

######  Definition 1 (see [RA14]).

Let $\mathcal{Z}$ be a finite discrete domain. Consider a target distribution
$\mathcal{D}$ supported on $\mathcal{Z}$ and a distribution $\mathcal{Q}$
supported on subsets of $\mathcal{Z}$. The sample oracle
$\mathrm{Samp}(\mathcal{Q};\mathcal{D})$ is called _Local Sampling Scheme
(LSS)_ and each time it is invoked, it returns an example $(x,S)$ such that:
(i) the set $S\subseteq\mathcal{Z}$ is drawn from $\mathcal{Q}$ and (ii) $x$
is drawn from the conditional distribution $\mathcal{D}_{S}$, where
$\mathcal{D}_{S}(x)=\mathcal{D}(x)\bm{1}\\{x\in S\\}/\mathcal{D}(S)$ and
$\mathcal{D}(S)=\sum_{x\in S}\mathcal{D}(x)$.

For simplicity, we mostly focus on the case where $\mathcal{Q}$ is supported
on a subset of $\mathcal{Z}\times\mathcal{Z}$, and often refer to
$\mathcal{Q}$ as the _pair distribution_. Then, $\mathcal{Q}$ naturally
induces an (edge weighted) undirected graph $G_{\mathcal{Q}}$ on
$\mathcal{Z}$, where $\\{u,v\\}\in E(G_{\mathcal{Q}})$, if $\\{u,v\\}$ is
supported on $\mathcal{Q}$, and has weight equal to the probability
$\mathcal{Q}(u,v)$ that $\\{u,v\\}$ is drawn from $\mathcal{Q}$. When we deal
with a pair distribution $\mathcal{Q}$, the above generative model lies in the
heart of the pairwise comparisons [FV86, Mar96, WJJ13]. The comparisons are
generated according to the Bradley-Terry-Luce (BTL) model (where the
comparisons have size $2$, [BT52, Luc59]): there is an unknown weight vector
$\mathcal{D}\in\mathbb{R}_{+}^{n}$ ($\mathcal{D}(u)$ indicates the quality of
item $u$) and, for two items $u,v$, the algorithm observes that $u$ beats $v$
with probability $\frac{\mathcal{D}(u)}{\mathcal{D}(u)+\mathcal{D}(v)}$.
Motivated by the importance of perfect simulation, we ask the following
question:

_(Q): Is there an efficient algorithm that draws i.i.d. samples $(x,S)$ from
$\mathrm{Samp}(\mathcal{Q};\mathcal{D})$ and generates a single sample
$y\sim\mathcal{D}$? _

To put our contribution into context, let us ask another question whose answer
is clear from previous work. What is the sample complexity of _learning_
$\mathcal{D}$ from a Local Sampling Scheme
$\mathrm{Samp}(\mathcal{Q};\mathcal{D})?$ Here, the goal is to estimate
$\mathcal{D}$ in some $L_{p}$ norm given such comparisons. Some works focus on
learning the re-parameterization $\bm{z}\in\mathbb{R}^{n}$ with
$z_{x}=\log(\mathcal{D}(x))$. Depending on the context, either the
normalization condition $\langle\bm{1},\mathcal{D}\rangle=1$ or
$\langle\bm{1},\bm{z}\rangle=0$ is used. [SBB+16] deal with the problem of
learning the re-parameterization $\bm{z}$ in the $L_{2}$ norm when $|S|=2$
(pairwise comparisons) or $|S|=k$ ($k$-wise comparisons with $k=O(1)$). The
sample complexity for learning $\bm{z}$ in $L_{2}$ norm from pairwise
comparisons is $n/(\lambda(\bm{Q})\epsilon^{2})$, where $\lambda(\bm{Q})$ is
the second smallest eigenvalue of a Laplacian matrix $\bm{Q}$, induced by the
pair distribution $\mathcal{Q}$ (the learning result is tight in a minimax
sense for some appropriate semi-norm [SBB+16]). In particular, the matrix
$\bm{Q}$ is defined as $\bm{Q}_{xy}=-\mathcal{Q}(x,y)$ and
$\bm{Q}_{xx}=\sum_{y\neq x}\mathcal{Q}(x,y)$. Hence, in the pairwise
comparisons setting, the sample complexity of learning incurs an overhead
associated with $\lambda(\bm{Q})$ (Table 1). The algorithm of [SBB+16] is a
pivotal component for our main result, as we will see later. [APA18] provide a
similar result for learning $\mathcal{D}$ in TV distance using a random-walk
approach. We remark that learning from pairwise comparisons requires some mild
conditions concerning the sampling process (LSS in our case). We review them
shortly below.

Sample Access |  Learning $(\epsilon>0)$ of $\bm{z}$ in $L_{2}$ |  Exact Sampling
---|---|---
Definition 1 with PCs |  $O\left(n/\epsilon^{2}\right)\cdot\frac{1}{\lambda(\bm{Q})}$ [SBB+16] |  $\widetilde{O}(n^{2})\cdot\frac{1}{\lambda(\bm{Q})}$ (Theorem 3)
Table 1: Learning and Exact Sampling from PCs of size 2.
$\widetilde{O}(\cdot)$ subsumes logarithmic factors.

##### Conditions for Local Sampling Schemes.

We provide a standard pair of conditions for Local Sampling Schemes (these or
similar conditions are also needed in the learning problem). We are interested
in LSSs that satisfy both of these conditions. However, some of our results
still hold when only the first condition is true; this will be more clear
later. 1 is a necessary information-theoretic condition for the pair
distribution $\mathcal{Q}$.

######  Assumption 1 (Identifiability).

The support of the pair distribution $\mathcal{Q}$ of the Local Sampling
Scheme contains a spanning tree, i.e., the induced graph $G_{\mathcal{Q}}$ is
connected. We define $\mathcal{E}:=E(G_{\mathcal{Q}}).$

Intuitively, the pair distribution $\mathcal{Q}$ needs to be supported on a
spanning tree. Equivalently, the associated Laplacian matrix $\bm{Q}$ should
have positive Fiedler eigenvalue $\lambda(\bm{Q})>0$. 2 is a condition about
the support of $\mathcal{Q}$ _and_ the target distribution $\mathcal{D}$; it
allows efficient learning of $\mathcal{D}$ from pairwise comparisons and is
related to the difficulty of estimating small $\mathcal{D}$-probabilities.

######  Assumption 2 (Efficiently Learnable).

There exists a constant $\phi>1$ such that target distribution $\mathcal{D}$
satisfies
$\frac{1}{\phi}\leq\max_{(x,y)\in\mathcal{E}}\frac{\mathcal{D}(x)}{\mathcal{D}(y)}\leq\phi\,,$
where $\mathcal{E}$ is the support of the distribution $\mathcal{Q}$.

From a dual perspective, 2 says that $\mathcal{Q}$ is supported only on edges
where the two corresponding $\mathcal{D}$-ratios are well controlled and
implies that any conditional distribution $\mathcal{D}_{\\{x,y\\}}$ has
bounded variance, i.e.,
$(1/\phi)/(1+\phi)^{2}\leq\frac{\mathcal{D}(x)\mathcal{D}(y)}{(\mathcal{D}(x)+\mathcal{D}(y))^{2}}\leq\phi/(1+1/\phi)^{2}$.
A quite similar property has been proposed in the problem of learning from
pairwise comparisons [SBB+16].

Contributions and Techniques. We provide an algorithmic answer to question
_(Q)_ regarding perfect sampling (under the above mild conditions) by showing
the following:

_There exists an efficient algorithm that draws
$\widetilde{O}(n^{2}/\lambda(\bm{Q}))$ samples from
$\mathrm{Samp}(\mathcal{Q};\mathcal{D})$ and outputs a perfect sample
$y\sim\mathcal{D}$._

The tool that makes this possible (and algorithmic) is a novel variant of the
elegant Coupling from the Past (CFTP) algorithm of [PW96]. CFTP is a technique
developed by Propp and Wilson [PW96, PW98], that provides an exact random
sample from a Markov chain, that has law the (unique) stationary distribution
(we refer to Section A.0.3 for an introduction to CFTP). To make our
contribution clear, we provide two theorems; 1 that directly and “naively”
applies CFTP (indicating the challenges behind efficient sample correction)
and a more sophisticated 2 that _combines CFTP and learning from pairwise
comparisons_ in order to get the result of Table 1. We show how to efficiently
generate “global” samples from the true distribution $\mathcal{D}$ from few
samples from a Local Sampling Scheme $\mathrm{Samp}(\mathcal{Q};\mathcal{D})$.
To this end, our novel algorithm combines the CFTP method [PW96, PW98], a
remarkable technique that generates exact samples from the stationary
distribution of a given Markov chain, with (i) an efficient algorithm that
estimates the parameters of Bradley-Terry-Luce ranking models from pairwise
(and, in general, $k$-wise) comparisons [SBB+16] and (ii) procedures based on
Bernoulli factories [DHKN17, MP05, NP05] and rejection sampling.

At a technical level, we observe (Section 2) that a Local Sampling Scheme
$\mathrm{Samp}(\mathcal{Q};\mathcal{D})$ naturally induces a canonical ergodic
Markov chain on $[n]$ with transition matrix $\bm{M}$. The probability of a
transition from a state $x$ to a state $y$ is equal to the probability
$\mathcal{Q}(x,y)$, that the pair $\\{x,y\\}$ is drawn from $\mathcal{Q}$,
times the probability $\mathcal{D}(y)/(\mathcal{D}(x)+\mathcal{D}(y))$, that
$y$ is drawn from the conditional distribution $\mathcal{D}_{\\{x,y\\}}$.
Then, the unique stationary distribution of the resulting Markov chain
coincides with the true distribution $\mathcal{D}$. Applying (a non-adaptive
version of) the Coupling from the Past algorithm, we obtain the following
result under 1 and 2 for a target distribution $\mathcal{D}$ over $[n]$. For
the matrix $\bm{M}$, we let $\Gamma(\bm{M})$ be its absolute spectral gap (see
Section 1.2 for extensive notation).

######  Informal Theorem 1 (Direct Exact Sampling from LSS).

With an expected number of
$\widetilde{O}\left(\frac{n^{2}}{\Gamma(\bm{M})}\right)$ samples from a Local
Sampling Scheme $\mathrm{Samp}(\mathcal{Q};\mathcal{D})$, there exists an
algorithm (Algorithm 1) that generates a sample distributed as in
$\mathcal{D}$.

The sample complexity of 1 is, for instance, attained by the instance of
Figure 1, which we will discuss shortly444For the instance of Figure 1 one can
reduce the sample complexity of our results by a factor of $n$, since the
state space has a natural ordering and one could apply monotone CFTP (it
suffices to coalesce two random processes and not $n$).. We remark that 1
holds more generally for Local Sampling Schemes that only satisfy 1 and not 2;
for a formal version of the above result, we refer to Theorem 2 and for a
discussion on the resulting sample complexity, we refer also to Remark 1. Let
us now see why 1 is quite unsatisfying.

In 1, the matrix $\bm{M}$ contains the transitions induced by
$\mathrm{Samp}(\mathcal{Q};\mathcal{D})$ and the transition $u\to v$ is
performed with probability
$\mathcal{Q}(u,v)\mathcal{D}(v)/(\mathcal{D}(u)+\mathcal{D}(v))$ and
$\Gamma(\bm{M})$ is the absolute spectral gap of the Markov chain’s transition
matrix $\bm{M}$ on which CFTP is performed. Note that both $\bm{M}$ and
$\Gamma(\bm{M})$ depend on $\mathcal{Q}$ and the target distribution
$\mathcal{D}$. In many natural cases, e.g., if $\mathcal{D}$ is multimodal or
includes many low probability points, the CFTP-based algorithm of 1 can be
quite inefficient! For instance, Figure 1 is a bad instance for the standard
CFTP algorithm. It corresponds to an instance where the support of
$\mathcal{Q}$ is the path graph over $[n]$ and the target distribution
$\mathcal{D}$ is bimodal, satisfying 1 and 2 with $\phi=2$. In this case,
$\bm{M}$ captures the transitions of a negatively biased nearest-neighbor
random walk on the path, i.e., the probability of moving to the boundary is
larger than moving to the interior. The coalescence probability of the two
extreme points is of order $\alpha^{n}$ for some $0<\alpha<1$. This
probability is connected to the maximum hitting time and, consequently, we get
that the mixing time is $\Omega(\alpha^{-n})$ [LP17] and so the coalescence
time is exponentially large. 1 depends on $\Gamma(\bm{M})$ and hence the
_exponential dependence on_ $n$ is reflected in the conductance of the chain,
which is exponentially small. So, can we obtain perfect $\mathcal{D}$-samples
using the oracle of Definition 1 efficiently?

$1/2$$1$$1/4$$2$$1/8$$3$$1/16$$4$$1/8$$5$$1/4$$6$$1/2$$7$ Figure 1: A bad
instance for 1.

As a technical contribution, we provide a novel exact sampling algorithm whose
convergence does not depend on $\mathcal{D}$. We show how to improve the
efficiency of CFTP by removing the dependence of its sample (and
computational) complexity on the structure of the true distribution
$\mathcal{D}$. In particular, in Section 3, we show the following under 1 and
2:

######  Informal Theorem 2 (Exact Sampling from LSS using Learning).

There exists an algorithm (Algorithm 2 & Algorithm 3) that draws an expected
number of $\widetilde{O}\left(\frac{n^{2}}{\lambda(\bm{Q})}\right)$ samples
from a Local Sampling Scheme $\mathrm{Samp}(\mathcal{Q};\mathcal{D})$, and
generates a sample distributed as in $\mathcal{D}$.

In the above, $\lambda(\bm{Q})$ is the second smallest eigenvalue (a.k.a., the
Fiedler eigenvalue) of $\bm{Q}$, the weighted Laplacian matrix of the graph
$G_{\mathcal{Q}}$ induced by $\mathcal{Q}$; this is the same matrix that
appears in the learning result of Table 1. _Crucially, the sample complexity
in the above result is polynomial in $n$_, provided that the pair distribution
$\mathcal{Q}$ is reasonable. For the path graph of Figure 1 with $\mathcal{Q}$
the uniform distribution, 1 yields an exponential sample complexity but the
spectrum of the matrix of 2 is
$\mathrm{spec}(\bm{Q})=\\{\frac{1}{n}\Theta(1-\cos(\pi i/n))\\}_{i=0}^{n-1}.$
Hence, since $1-\cos(x)=\sin^{2}(x/2)\approx x^{2}$, we get that the sample
complexity is $\mathrm{poly}(n).$ We also note that the runtime of the
algorithm is sample polynomial.

##### Technical Overview.

A first key idea towards establishing 2 is to exploit _learning_ in order to
accelerate the convergence of CFTP. We first apply (a modification of) the
gradient-based algorithm of [SBB+16] to the empirical log-likelihood
objective, which estimates the parameters of BTL ranking models from PCs, with
samples from $\mathrm{Samp}(\mathcal{Q};\mathcal{D})$. As a result of this
learning algorithm (whose efficiency requires 2, as in [SBB+16]), we get a
probability distribution $\widetilde{\mathcal{D}}$, which approximates
$\mathcal{D}$ with small relative error, in the sense that
$\mathcal{D}(x)\in(1\pm\epsilon)\cdot\widetilde{\mathcal{D}}(x)$, for any $x$
in $\mathcal{D}$’s support. As we saw in Table 1, the sample complexity will
be roughly of order $O(n/(\lambda(\bm{Q})\epsilon^{2})))$. To be precise, this
learning phase gives us a sample complexity term of order
$O(n^{2}/\lambda(\bm{Q}))$ and corresponds to a single execution of the above
learning algorithm for the target $\mathcal{D}$ in relative error with
accuracy $\epsilon=\Theta(1/\sqrt{n})$ (the accuracy’s choice will be more
clear later).

The new idea is that we can use the output of the learning algorithm
$\widetilde{\mathcal{D}}$ and a Bernoulli factory mechanism to transform the
Markov chain induced by $\mathrm{Samp}(\mathcal{Q};\mathcal{D})$ into an
ergodic chain with almost uniform stationary probability distribution, whose
transition probabilities (and mixing time) essentially do not depend on
$\mathcal{D}$! To do so, for each pair of states $x,y$, with
$(x,y)\in\mathcal{E}$ and
$\widetilde{\mathcal{D}}(x)>\widetilde{\mathcal{D}}(y)$, we can use Bernoulli
downscaling (which constitutes an essential building block for general
Bernoulli factories, see e.g., [DHKN17, MP05, NP05]) to transform the
(unknown) probability of a transition from $y$ to $x$ to
$\mathcal{D}(x)/(\mathcal{D}(x)+\mathcal{D}(y))\cdot\widetilde{\mathcal{D}}(y)/\widetilde{\mathcal{D}}(x)$
(conditional that the edge $(x,y)$ is drawn). Since
$\widetilde{\mathcal{D}}(x)\approx\mathcal{D}(x)$, for all states $x$, this
makes the probability of a transition from $x$ to $y$ almost equal to the
probability of a transition from $y$ to $x$, which, in turn, implies that _the
stationary distribution of the modified Markov chain is almost uniform_. Thus,
we speed up the CFTP algorithm by removing the modes due to the landscape of
the target $\mathcal{D}$ and can get an exact sample from the stationary
distribution of the modified Markov chain with
$O\left(n\log(n)/\lambda(\bm{Q})\right)$ samples. Namely, the sample
complexity of this parameterized variant of the CFTP algorithm does not depend
on the structure of $\mathcal{D}$ anymore (of course, it does depend on the
size $n$ of $\mathcal{D}$’s support). The above constitute our second key idea
and the intuition behind this re-weighting comes from the method of simulated
tempering [MP92, GLR18b, GLR18a].

Our third idea deals with the output of the algorithm. Since we want an
exact/clean sample from $\mathcal{D}$, we cannot simply print the output of
the above process. As a last step, we use rejection sampling, guided by the
distribution $\widetilde{\mathcal{D}}$, and apply the “reverse” modification,
bringing the distribution of the resulting sample back to $\mathcal{D}$ (see
Algorithm 2). With high probability, this last rejection sampling will succeed
after $\Theta(n)$ executions of the parameterized CFTP algorithm. However,
note that we do not have to execute the learning algorithm again. Hence, the
total sample complexity is equal to the cost of the learning phase (that is
performed only once and its output is given as input to the parameterized CFTP
algorithm; see Algorithm 2) and the cost of the modified CFTP multiplied by
$\Theta(n)$ (due to rejection sampling), which yields a total of
$\widetilde{O}(n^{2}/\lambda(\bm{Q}))$ samples. The above modification of the
Markov chain induced by $\mathrm{Samp}(\mathcal{Q};\mathcal{D})$, which
improves possibly significantly its mixing time, is simple and general and the
parameterized variation of the standard CFTP algorithm (Algorithm 2) might be
of independent interest. In Appendix G, we discuss how to extend the analysis
of the main Algorithm 2 to sets of size larger than $2$.

###  1.1 Related Work

Our results are related to and draw ideas from several areas of theoretical
computer science:

##### Exact Sampling.

Exact sampling comprises a well-studied field [Ken05, Hub16] with numerous
applications in theoretical computer science (e.g., to approximate counting
[Hub98]). One of our main tools, the Coupling From the Past algorithm555 We
would like to shortly mention an interesting example that shows how perfect
sampling from pairwise comparisons is possible from an algebraic perspective
which is different from the standard random-walk viewpoint: consider the case
where $n=3$ and the discrete distribution $\mathcal{D}$ is the vector
$(\mathcal{D}(a),\mathcal{D}(b),\mathcal{D}(c))$. Assume that $\mathcal{Q}$ is
supported on $\\{a,b\\}$ and $\\{b,c\\}$ with equal probability. The following
process generates an exact sample from $\mathcal{D}$: given a sample $x$ from
$\mathcal{D}_{\\{a,b\\}}$ and $y$ from $\mathcal{D}_{\\{b,c\\}}$, we repeat
until $(x,y)\neq(a,c)$. In this case, $b$ will be observed in at least one
coordinate and the algorithm outputs the other coordinate. We claim that this
process outputs a perfect sample from $\mathcal{D}$. An algebraic intuition
behind this random process is that $(a+b)(b+c)=b(a+b+c)+ac$, where $+$
corresponds to the logical or; the bad event is the term $ac$, while the good
event contains $b$ and the perfect sample ($a$ or $b$ or $c$). This is the
generating polynomial of the above sampling process. Formally, the probability
that this process outputs $a$ is
$\Pr[a]\propto\frac{\mathcal{D}(a)}{\mathcal{D}(a)+\mathcal{D}(b)}\frac{\mathcal{D}(b)}{\mathcal{D}(b)+\mathcal{D}(c)}\sum_{i=0}^{\infty}\left(\frac{\mathcal{D}(a)}{\mathcal{D}(a)+\mathcal{D}(b)}\frac{\mathcal{D}(c)}{\mathcal{D}(b)+\mathcal{D}(c)}\right)^{i}$
and note that $\Pr[x]\propto\mathcal{D}(x)$ for any $x\in\\{a,b,c\\}$ and all
have the same normalization constant. This verifies our claim. , is the
breakthrough result of [PW96, PW98], which established the applicability of
perfect simulation, and made efficient simulation of challenging problems
(e.g., Ising model) feasible. After this fundamental result, the literature of
perfect simulation flourished and various exact sampling procedures have been
developed, see e.g., [PW96, GM99, Wil00, HS00, FH00, Hub04], and the
references therein.

##### Learning and Ranking from Pairwise Comparisons.

There is a vast amount of literature on ranking and parameter estimation from
pairwise comparisons (see e.g., [FV86, Mar96, Cat12] and the references
therein). A significant challenge in that field is to come up with models that
accurately capture uncertainty in pairwise comparisons. The standard ones are
the Bradley-Terry-Luce (BTL) model [BT52, Luc59] and the Thurstone model
[Thu27]. Our work is closely related to previous work on the BTL model [RA14,
SBGW16, DKNS01, NOS12, Hun04, JLYY11, APA18, LSR21] and our sample complexity
rates depend on the Fiedler eigenvalue of a pairwise comparisons’ matrix, as
in [HOX14, SBB+16, VY16, NOS17, VYZ20]. For more details on ranking
distributions, see e.g., [LM18, BFFSZ19, FKS21, FKP22].

##### Truncated Statistics.

Our work falls in the research agenda of efficient learning and exact sampling
from truncated samples. Recent work has obtained efficient algorithms for
efficient learning from truncated [DGTZ18, KTZ19, DGTZ19, IZD20, DRZ20,
BDNP21, NP19, FKT20, DKTZ21] and censored [MMS21, FKKT21, Ple21] samples.
Closer to our work, [FKT20] also deal with the question of exact sampling from
truncated samples.

##### Correction Sampling.

Sample correction is a field closely related to our work. This area of
research deals with the case where input data are drawn from a noisy or
imperfect source and the goal is to design _sampling correctors_ [CGR18] that
operate as filters between the noisy data and the end user. These algorithms
exploit the structure of the distribution and make “on-the-fly” corrections to
the drawn samples. Moreover, the problem of local correction of data has
received much attention in the contexts of self-correcting programs, locally
correctable codes, and local filters for graphs and functions (see e.g.,
[BLR93, SS10, JR13, ACCL04, BGJ+10, GLT20]).

##### Conditional Sampling.

Conditional sampling deals with the _adaptive_ learning analog of LSS, where
the goal is again to learn an underlying discrete distribution from
conditional/truncated samples, but the learner can change the truncation set
on demand; see e.g., [Can15, FJO+15, ACK15b, BC18, ACK15a, GTZ17, KT19,
CCK+19].

###  1.2 Notation

In this section, we provide the basic notation we are going to use in the
technical part.

##### General Notation.

We define $[n]:=\\{1,\dots,n\\}$. We denote vectors with small bold letters
$\bm{x}$ and matrices with bold letters $\bm{P}=[\bm{P}_{xy}]$. For a vector
$\bm{x}$, we let $x_{i}=x(i)$ be its $i$-th entry. For $p\in\\{1,2,\infty\\}$,
we denote the $L_{p}$ norm by $\|.\|_{p}$. For a matrix $\bm{A}$, let
$\|\bm{A}\|_{2}$ denote its spectral norm (largest singular value). We
consider graphs $G=(V,E)$, with $|V|=n$ vertices, whose associated Laplacian
matrix is denoted $\bm{\Lambda}(G)$. The edge set of the graph $G$ is the set
$E(G)$. In Section A.0.2, we have included some basic definitions and facts
about random walks and Markov chains. We define the weighted Laplacian matrix
of $G_{\mathcal{Q}}$ for a distribution $\mathcal{Q}$ as

$\bm{Q}_{xy}=-\mathcal{Q}(x,y)~{}\text{for any $x\neq y$ and
}~{}\bm{Q}\bm{1}=\bm{0}\,,$ (1.1)

##### Distributions and Distances.

The probability simplex is denoted by $\Delta^{n}$. We let
$\mathrm{supp}(\mathcal{D})$ denote the support of a distribution
$\mathcal{D}$. The total variation distance between two discrete distributions
$\mathcal{D}_{1},\mathcal{D}_{2}$ over $[n]$ is equal to
$\mathrm{TV}(\mathcal{D}_{1},\mathcal{D}_{2})=\max_{S\subseteq[n]}(\mathcal{D}_{1}(S)-D_{2}(S))=\frac{1}{2}\|\mathcal{D}_{1}-\mathcal{D}_{2}\|_{1}$.
We let $\mathcal{D}(S)=\sum_{x\in S}\mathcal{D}(x)$ and $\mathcal{D}_{S}$ is
the conditional distribution on the set $S$, i.e.,
$\mathcal{D}_{S}(x)=\frac{\mathcal{D}(x)}{\mathcal{D}(S)}\bm{1}\\{x\in S\\}.$
The vector $\bm{z}\in\mathbb{R}^{n}$ with $z_{i}=\log(\mathcal{D}(i))$ is
called the natural parameter vector of $\mathcal{D}$.

##### Random Walks.

For a reversible transition matrix $\bm{P}$, let its (real) spectrum be
$1=\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{n}\geq-1$. We define the
_absolute spectral gap_ of $\bm{P}$ to be the difference
$\Gamma(\bm{P})=1-\max\\{\lambda_{2}(\bm{P}),|\lambda_{n}(\bm{P})|\\}$. If
$\bm{P}$ is aperiodic and irreducible, then $\Gamma(\bm{P})>0$. One could also
define the _spectral gap_ $1-\lambda_{2}(\bm{P})$. The two gaps are equal when
the chain is lazy, i.e., when $\forall i\in[n]$, $\bm{P}_{ii}\geq 1/2$. We let
$\lambda(\bm{L})$ be the second _smallest_ eigenvalue of a Laplacian matrix
$\bm{L}$, a.k.a., its Fiedler eigenvalue. The Laplacian matrix
$\bm{L}\in\mathbb{R}^{n\times n}$ is positive semi-definite and induces a
semi-norm on $\mathbb{R}^{n}$ with
$\|\bm{v}\|_{L}:=\sqrt{\bm{v}^{T}\bm{L}\bm{v}}$ for $\bm{v}\in\mathbb{R}^{n}$.
Recall that a semi-norm differs from a norm in that the semi-norm of a non-
zero element is allowed to be zero. In Algorithm 1, Algorithm 2 and Algorithm
7, we use the notation $F_{t}(x)\leftarrow F_{t+1}(w)$ for some $t<0$. This
means that we append the path of $x$ at time $t$ by adding the transition
$x\to w$ and then the path of $w$ from time $t+1$ up to $0$.

##### Matrix Operators.

We denote by $\bm{A}\odot\bm{B}=(\bm{A}_{ij}\bm{B}_{ij})_{1\leq i\leq j\leq
n}$ the standard Hadamard matrix product and by $\bm{A}\circ\bm{B}$ a
variation of the Hadamard matrix product, where the off-diagonal entries are
equal to those of the standard Hadamard product $\bm{A}\odot\bm{B}$, but the
diagonal entries correspond to the diagonal matrix with entries $(-\sum_{j\neq
i}\bm{A}_{ij}\bm{B}_{ij})_{i\in[n]}$. Finally, we let
$\bm{A}\otimes\bm{B}=\bm{I}-\bm{A}\circ\bm{B}.$

##  2 An Exact Sampling Algorithm using CFTP

Consider a discrete target distribution $\mathcal{D}$. For pair distributions
$\mathcal{Q}$, the Local Sampling Scheme
$\mathrm{Samp}(\mathcal{Q};\mathcal{D})$ can be regarded as a graph $G=(V,E)$
with $V=\mathrm{supp}(\mathcal{D})$ and
$\mathcal{E}=\mathrm{supp}(\mathcal{Q})$. We let $V=[n]$. Let $\bm{M}$ denote
the canonical transition matrix of the Markov chain, associated with the
oracle $\mathrm{Samp}(\mathcal{Q};\mathcal{D})$. The entries of
$\bm{M}=[\bm{M}_{xy}]_{x,y\in[n]}$ are defined as:

$\bm{M}_{xy}=\mathcal{Q}(x,y)\frac{\mathcal{D}(y)}{\mathcal{D}(x)+\mathcal{D}(y)}~{}\textnormal{for
any $x\neq y$ and }~{}\bm{M}_{xx}=1-\sum_{y\neq x}\bm{M}_{xy}~{}\textnormal{
for }x\in[n]\,.$ (2.1)

Observe that the transition from $x$ to $y$ is performed when both the edge
$\\{x,y\\}$ is drawn from $\mathcal{Q}$, with probability proportional to
$\mathcal{Q}(x,y)$ and $y$ is drawn from the conditional distribution
$\mathcal{D}_{\\{x,y\\}}$ . The Markov chain, whose transitions are described
by $\bm{M}$, has some notable properties: it has $\mathcal{D}$ as stationary
distribution and is ergodic. We can invoke the fundamental CFTP algorithm over
the Markov chain of [PW96, PW98] to obtain a perfect sample from
$\mathcal{D}$. We extensively discuss CFTP in Section A.0.3 for the interested
reader. CFTP yields the following (unsatisfying) result, whose proof can be
found at the Appendix B.

Algorithm 1 Exact Sampling from Local Sampling Schemes of 1 & Theorem 2

1:procedure ExactSampler() $\triangleright$ _Sample access to the LSS oracle
$\mathrm{Samp}(\mathcal{Q};\mathcal{D})$._

2: $t\leftarrow 0$

3: $F_{0}(x)\leftarrow x$, for any $x\in[n]$

4: $\textbf{while}~{}F_{t}~{}\text{has not coalesced}~{}\textbf{do}$
$\triangleright$ _While no coalescence has occured._

5: $~{}~{}~{}~{}~{}~{}~{}~{}t\leftarrow t-1$

6: Draw sample $(i,j,w)$ with $(i,j)\in\mathcal{E}$ and $w\in\\{i,j\\}$

7: for $x=1\ldots n$ do$\triangleright$ _In order to update state $x$._

8: if $x\notin\\{i,j\\}~{}\mathrm{or}~{}x=w$ then $F_{t}(x)\leftarrow
F_{t+1}(x)$ $\triangleright$ _Append $F_{t}$ in the past._

9: else $F_{t}(x)\leftarrow F_{t+1}(w)$ $\triangleright$ _New transition with
probability
$\mathcal{Q}(x,w)\frac{\mathcal{D}(w)}{\mathcal{D}(x)+\mathcal{D}(w)}$._

10: end

11: end

12: Output $F_{t}(1)$ $\triangleright$ _Output the perfect sample._

######  Theorem 2 (Direct Exact Sampling from LSS).

Let $\mathcal{D}_{\mathrm{min}}$ be the minimum value of the target
distribution $\mathcal{D}$. Under 1 and for any $\delta>0$, Algorithm 1 draws,
with probability at least $1-\delta$,
$N={O}\left(\frac{n\log(1/\mathcal{D}_{\mathrm{min}})}{\Gamma(\bm{M})}\log(\frac{1}{\delta})\right)$
samples from a Local Sampling Scheme $\mathrm{Samp}(\mathcal{Q};\mathcal{D})$
over $[n]$, runs in time polynomial in $N$ and outputs a sample distributed as
in $\mathcal{D}$.

Note that in the above result (Theorem 2), the confidence parameter $\delta$
corresponds to the number of samples required and not the quality of the
output sample (the sample is perfect with probability $1$). We comment that
the proof of the above Theorem is nearly straightforward and follows from the
analysis of the CFTP algorithm; the purpose of Theorem 2 is to indicate that a
direct application of CFTP would potentially lead to an exponential number of
samples due to the shape of the target distribution. Our main technical
contribution begins in the next section: We give an exact sampler that
surpasses the dependence on the target’s structure using learning and
rejection sampling.

##  3 Improving the Exact Sampling Algorithm using Learning

The main drawback of Theorem’s 2 algorithm is that its sample complexity
depends on both $\bm{Q}$ and $\bm{D}$ and consequently in many natural
scenaria the sample complexity may be exponentially large in $n$. Next, we
present Algorithm 2, which will allow us to remove the dependence on $\bm{D}$
for the Markov chain and reduce the sample complexity to polynomial in the
domain size. The crucial differences compared to Algorithm 1 are the  orange
lines. Algorithm 2 is a parameterized extension of CFTP and its input is a
vector $\bm{p}\in[0,1]^{n}$. We remark that when $\bm{p}$ is the all-ones
vector $\bm{1}$, we obtain the previous Algorithm 1.

Algorithm 2 Parameterized Extension of Algorithm 1

1: Input: $\bm{p}\in[0,1]^{n}$ $\triangleright$ _For instance, this may be the
output of Algorithm 5._

2:procedure ExactSampler($\bm{p}$) $\triangleright$ _Sample access to the LSS
oracle $\mathrm{Samp}(\mathcal{Q};\mathcal{D})$._

3: $t\leftarrow 0$, $F_{0}(x)\leftarrow x$, for any $x\in[n]$

4: $\textbf{while}~{}F_{t}~{}\textnormal{has not coalesced}~{}\textbf{do}$
$\triangleright$ _While no coalescence has occured._

5: $~{}~{}~{}~{}~{}~{}~{}~{}t\leftarrow t-1$

6: Draw sample $(i,j,w)$ with $(i,j)\in\mathcal{E}$ and $w\in\\{i,j\\}$

7: for $x=1\ldots n$ do$\triangleright$ _In order to update state $x$._

8: if $x\notin\\{i,j\\}~{}\mathrm{or}~{}x=w$ then $F_{t}(x)\leftarrow
F_{t+1}(x)$

9: else
$F_{t}(x)\leftarrow\left\\{\begin{array}[]{lr}\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}F_{t+1}(w),&\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\textnormal{with
probability
}\min\\{p(x)/p(w),1\\}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\\\
\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}F_{t+1}(x),&\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\textnormal{otherwise
}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color<EMAIL_ADDRESS>
10: end

11: end

12: Draw $C\sim\mathrm{Be}(p(F_{t}(1)))$ $\triangleright$ _$\mathrm{Be}(p)$ is
a $p$-biased Bernoulli coin_.

13: if $C=1$ then Output $F_{t}(1)$ else Output $\perp$ $\triangleright$
_Output the perfect sample or Fail._

We perform Algorithm 2 using as input parameter $\bm{p}$ an estimate
$\widetilde{\mathcal{D}}$ of the target $\mathcal{D}$ with (relative) error of
order $\epsilon=1/\sqrt{n}$ (this estimate is obtained via Theorem 4 &
Algorithm 5 using $O(n^{2}/\lambda(\bm{Q}))$ samples from
$\mathrm{Samp}(\mathcal{Q};\mathcal{D})$, as we will see later). Under the
perspective of LSSs as random walks on an irreducible Markov chain, Algorithm
2 proceeds by executing CFTP. Recall that in each draw from
$\mathrm{Samp}(\mathcal{Q};\mathcal{D})$, a transition can be realized by the
matrix $\bm{M}$ . The fact that this transition depends on the target
$\mathcal{D}$ is pathological, as we observed previously. The algorithm,
instead of performing this transition (as the non-adaptive CFTP algorithm of
Theorem 2 does), performs a Bernoulli factory mechanism (downscaling), using
the estimates $\widetilde{\mathcal{D}}$ in order to make each transition
almost uniform (see Line 9 of Algorithm 2). This change introduces some
(known) bias to the random walk and changes the stationary distribution from
$\mathcal{D}$ to an almost uniform one. By performing the CFTP method, the
algorithm iteratively reconstructs the past of the infinite simulation, until
all the simulations have coalesced at time $t=0$. When coalescence occurs, the
algorithm has an exact sample from the biased (almost uniform) stationary
distribution. Since the introduced bias is known (it corresponds to the
inverse of the estimates $\widetilde{\mathcal{D}}$), the algorithm accepts the
sample using rejection sampling, guided by $\widetilde{\mathcal{D}}$, so that
the final sample is distributed according to the target distribution
$\mathcal{D}$ (see Line 13 of Algorithm 2). Specifically, we show (for the
proof, see Appendix C) that:

######  Theorem 3 (Exact Sampling from LSS using Learning).

Under 1 and 2, for any $\delta>0$, there exists an algorithm (Algorithm 3)
that draws, with probability at least $1-\delta$,
$N={O}\left(n^{2}\log^{2}(n)/\lambda(\bm{Q})\cdot\log(1/\delta)\right)$
samples from a Local Sampling Scheme $\mathrm{Samp}(\mathcal{Q};\mathcal{D})$,
runs in time polynomial in $N$, and outputs a sample distributed as in
$\mathcal{D}$.

Algorithm 3 Exact Sampling from Local Sampling Schemes of 2 & Theorem 3

1:procedure ExactSamplerWithLearning($\delta$) $\triangleright$ _The algorithm
of Theorem 3_.

2:
$\widetilde{\mathcal{D}}\leftarrow\textsc{Learn}(\epsilon:=1/\sqrt{n},\delta)$
$\triangleright$ _Learn $\mathcal{D}$ in relative error with Algorithm 5
(Theorem 4)._

3: x $\leftarrow\perp$

4: while x $=\perp$ repeat

5: x $\leftarrow$ ExactSampler$(\widetilde{\mathcal{D}})$ $\triangleright$
_Call Algorithm 2._

6: Output x

We proceed with the analysis of Algorithm 3, which can be decomposed into two
parts: the Learning Step (Line 2) and several iterations of Algorithm 2 (Line
5).

##### Learning Phase.

For two distributions $\mathcal{D},\widetilde{\mathcal{D}}$ with ground set
$[n]$, we introduce the sequence/list (of length $n$)
$1-\mathcal{D}/\widetilde{\mathcal{D}}:=(1-\mathcal{D}(x)/\widetilde{\mathcal{D}}(x))_{x\in[n]}.$
Observe that the pair of sequences
$(1-\mathcal{D}/\widetilde{\mathcal{D}},1-\widetilde{\mathcal{D}}/\mathcal{D})$
captures the relative error between the two distributions. The sample
complexity of the task of learning $\mathcal{D}$ in $\epsilon$-relative error
is summarized by the following:

######  Theorem 4 (Learning Phase).

For any $\epsilon,\delta>0$, there exists an algorithm (Algorithm 5) that
draws
$N=O\left(\frac{n}{\lambda(\bm{Q})\epsilon^{2}}\log(\frac{1}{\delta})\right)$
samples from $\mathrm{Samp}(\mathcal{Q};\mathcal{D})$ satisfying Assumptions 1
and 2, runs in time polynomial in $N$, and, with probability at least
$1-\delta$, computes an estimate $\widetilde{\mathcal{D}}$ of the target
distribution $\mathcal{D}$, that satisfies the relative error guarantee
$\max\big{\\{}\|1-\mathcal{D}/\widetilde{\mathcal{D}}\|_{\infty},\|1-\widetilde{\mathcal{D}}/\mathcal{D}\|_{\infty}\big{\\}}\leq\epsilon$.

We execute _once_ the learning algorithm (Algorithm 5) that estimates
$\mathcal{D}$ in relative error and we choose accuracy $\epsilon=1/\sqrt{n}$.
We would like to shortly comment on the choice of the accuracy value: An
accuracy of constant order would still potentially yield a biased random walk.
After the re-weighting step, our goal is to show that the target distribution
becomes almost uniform and that the absolute spectral gap of the reformed
matrix is of order $\lambda(\bm{Q})$ (see the proof of Lemma 15). In order to
ensure this, any accuracy that vanishes with $n$ is sufficient. We chose to
set $\epsilon=1/\sqrt{n}$ so that the complexity of learning
$(\frac{n}{\lambda(\bm{Q})\epsilon^{2}})$ matches that of the exact sampling
procedure $(\frac{n^{2}}{\lambda(\bm{Q})})$. This learning step costs
$n^{2}/\lambda(\bm{Q})$ samples and the estimate $\widetilde{\mathcal{D}}$ is
used as input to Algorithm 2. For the proof of Theorem 4, we refer to the
Appendix E. We can now move to Line 5.

##### Executing Algorithm 2.

We next focus on a single execution of the parameterized CFTP algorithm. Our
goal is to transform the Markov chain of the Local Sampling Scheme defined in
Equation (2.1) to a modified Markov chain with an almost uniform stationary
distribution. We are going to provide some intuition. Applying Theorem 4, we
can consider that, for any $x\in[n]$, there exists a coefficient
$\widetilde{\mathcal{D}}(x)\approx\mathcal{D}(x)$. The idea is to use
$\widetilde{\mathcal{D}}(x)$ and make the stationary distribution of the
modified chain close to uniform. Intuitively, this transformation should
speedup the convergence of the CFTP algorithm. The downscaling method of
Bernoulli factories gives us a tool to do this modification.

We can implement the modified Markov chain via downscaling as follows:
Consider an edge $\\{x,y\\}$ with transition probability pair
$(p_{xy},p_{yx})$. Without loss of generality, we assume that
$\widetilde{\mathcal{D}}(y)>\widetilde{\mathcal{D}}(x)$ (which intuitively
means that we should expect that $p_{xy}>p_{yx}$). Then, the downscaler leaves
$p_{yx}$ unchanged and reduces the mass of $p_{xy}$ to make the two
transitions almost balanced. Consider the LSS transition matrix $\bm{M}$ with
$\bm{M}_{xy}=\mathcal{Q}(x,y)\frac{\mathcal{D}(y)}{\mathcal{D}(x)+\mathcal{D}(y)}$
and $\bm{M}_{xx}=1-\sum_{y\neq x}\bm{M}_{xy}$. Also, let
$\widetilde{\mathcal{D}}$ be an estimate for the distribution $\mathcal{D}$
(Theorem 4). For the pair $(x,y)$, we modify the transition probability
$p_{xy}:=\bm{M}_{xy}$, only if $p_{xy}>p_{yx}$, to be equal to the following:
$\widetilde{p}_{xy}=\mathcal{Q}(x,y)\frac{\mathcal{D}(y)}{\mathcal{D}(x)+\mathcal{D}(y)}\frac{\widetilde{\mathcal{D}}(x)}{\widetilde{\mathcal{D}}(y)}\approx\mathcal{Q}(x,y)\frac{\mathcal{D}(x)}{\mathcal{D}(x)+\mathcal{D}(y)}=p_{yx}\,,$
where we use that $\widetilde{\mathcal{D}}(x)\approx\mathcal{D}(x)$ and
$\widetilde{\mathcal{D}}(y)\approx\mathcal{D}(y)$. The transition probability
from $x$ to $y$ corresponds to a Bernoulli variable $\mathrm{Be}(p_{xy})$,
which is downscaled by
$\widetilde{\mathcal{D}}(x)/\widetilde{\mathcal{D}}(y)<1$. The modified
transition probability $\widetilde{p}_{xy}$ can be implemented by drawing a
$\Lambda\sim\mathrm{Be}(\widetilde{\mathcal{D}}(x)/\widetilde{\mathcal{D}}(y))$
and then drawing a $P\sim\mathrm{Be}(p_{xy})$ (from
$\mathrm{Samp}(\mathcal{Q};\mathcal{D})$) and, finally, realizing the
transition from $x$ to $y$ only if $\Lambda P=1$. This implementation is valid
since the two sources of randomness are independent. So, we perform a
downscaled random walk based on the transition matrix $\widetilde{\bm{M}}$
with $\widetilde{\bm{M}}_{xy}=\widetilde{p}_{xy}$ (the “small” probability
$p_{yx}$ remained intact). We call $\widetilde{\bm{M}}$ the
$\widetilde{\mathcal{D}}$-scaling of $\bm{M}$. Lemma 5 summarizes the key
properties of $\widetilde{\bm{M}}$. The definition of the mixing time, used
below, can be found at the notation Section 1.2 and the full proof can be
found at the Appendix D.

######  Lemma 5 (Properties of Rescaled Random Walk).

Let $\mathcal{D}$ be a distribution on $[n]$ and consider an
$\epsilon$-relative approximation $\widetilde{\mathcal{D}}$ of $\mathcal{D}$,
as in Theorem 4 with $\epsilon=1/\sqrt{n}$. Consider the transition matrix
$\bm{M}$ of the Local Sampling Scheme $\mathrm{Samp}(\mathcal{Q};\mathcal{D})$
(see Equation (2.1)), and let $\widetilde{\bm{M}}$ be the
$\widetilde{\mathcal{D}}$-scaling of $\bm{M}$. Then, the transition matrix
$\widetilde{\bm{M}}$ (i) has stationary distribution $\widetilde{\pi}_{0}$
with $\widetilde{\pi}_{0}(x)=\Theta(1/n)$ for all $x\in[n]$ and (ii) the
mixing time of the transition matrix $\widetilde{\bm{M}}$ is
$T_{\mathrm{mix}}(\widetilde{\bm{M}})=O\left(\frac{\log(n)}{\lambda(\bm{Q})}\right)$.

###### Proof Sketch.

Set $J_{xy}=\mathcal{D}(x)+\mathcal{D}(y).$ For a pair $(x,y)\in\mathcal{E}$,
before the downscaling, the pairwise $(x,y)$-comparison corresponds to the
random coin $\left(\mathcal{D}(y)/J_{xy},\mathcal{D}(x)/J_{xy}\right)$ and let
$\mathcal{D}(y)>\mathcal{D}(x)$. After the downscaling, this coin becomes
(locally) almost fair, i.e.,
$\left(\mathcal{D}(x)/J_{xy}+\xi_{xy},\mathcal{D}(x)/J_{xy}\right).$ Hence,
the modified transition matrix $\widetilde{\bm{M}}$ can be written as
$\widetilde{\bm{M}}=\widetilde{\bm{D}}\otimes\bm{Q}+\bm{Q}\circ[\xi_{xy}]\,,$
where $\widetilde{\bm{D}}$ is a symmetric matrix with
$\widetilde{\bm{D}}_{xy}=\mathrm{min}\\{\mathcal{D}(x),\mathcal{D}(y)\\}/(\mathcal{D}(x)+\mathcal{D}(y))=\widetilde{\bm{D}}_{yx}$,
$\bm{Q}$ is the Laplacian matrix of the graph $G_{\mathcal{Q}}$ and
$[\xi_{xy}]$ is a matrix with entries $\xi_{xy}$ where $\xi_{xy}$ is only non-
zero in the modified transitions $x\to y$, i.e., if $p_{xy}>p_{yx}$ then
$\widetilde{\bm{M}}_{xy}=\mathcal{Q}(x,y)\cdot(\mathcal{D}(x)/J_{xy}+\xi_{xy})$;
for this pair, we set $\xi_{yx}=0$. Finally, $\bm{Q}\circ[\xi_{xy}]$ denotes
the modified Hadamard product666We have that
$\Big{(}\bm{Q}\circ\bm{A}\Big{)}_{xy}=\mathcal{Q}_{xy}A_{xy},\Big{(}\bm{Q}\circ\bm{A}\Big{)}_{xx}=-\sum_{y\neq
x}\mathcal{Q}_{xy}A_{xy}.$ between $\bm{Q}$ and $[\xi_{xy}]$. Due to the
accuracy $\epsilon$ of the learning algorithm, we have that
$|\xi_{xy}|=O(1/\sqrt{n})$.

We analyze (i) the stationary distribution and (ii) the mixing time. For Part
(i), the chain $\widetilde{\bm{M}}$ remains irreducible and the detailed
balance equations of $\widetilde{\bm{M}}$ satisfy
$\widetilde{p}_{xy}/\widetilde{p}_{yx}=\mathcal{D}(y)\cdot\widetilde{\mathcal{D}}(x)/(\widetilde{\mathcal{D}}(y)\cdot\mathcal{D}(x)).$
Since $\widetilde{\mathcal{D}}$ is an $\epsilon$-relative approximation of
$\mathcal{D}$ with $\epsilon=O(1/\sqrt{n})$, it holds that for any $x\in[n]$:
$\mathcal{D}(x)/\widetilde{\mathcal{D}}(x)\in[1-\epsilon,1+\epsilon]$. Hence,
the unique stationary distribution $\widetilde{\pi}_{0}$ satisfies:
$\widetilde{\pi}_{0}(x)/\widetilde{\pi}_{0}(y)\in[1-\Theta(1/\sqrt{n}),1+\Theta(1/\sqrt{n})]$
for any $x,y\in[n]$. So, it must be the case that
$\widetilde{\pi}(x)=\Theta(1/n)$.

For Part (ii), in order to show that the mixing time is of order
$O(\log(n)/\lambda(\bm{Q}))$, it suffices to show (by Lemma 9) that (i) the
minimum value of the stationary distribution is $O(1/n)$ and (ii) the absolute
spectral gap $\Gamma(\widetilde{\bm{M}})$ of $\widetilde{\bm{M}}$ satisfies
$\Gamma(\widetilde{\bm{M}})=\Omega(\lambda(\bm{Q}))$, where $\lambda(\bm{Q})$
is the minimum non-zero eigenvalue of the Laplacian matrix $\bm{Q}$. The first
property follows directly from the closeness of the stationary distribution of
$\widetilde{\bm{M}}$ to the uniform one. We now focus on the latter. Our goal
is to control the absolute spectral gap $\Gamma(\widetilde{\bm{M}})$. In fact,
this matrix can be seen as a perturbed version of the matrix
$\widetilde{\bm{D}}\otimes\bm{Q}$ as we discussed in the beginning of the
proof sketch. The perturbation noise matrix $\bm{E}=\bm{Q}\circ[\xi_{xy}]$ is
a matrix whose entries contain (among others) the approximation errors
$[\xi_{xy}]$. Hence, we can use Weyl’s inequality [Tao12] in order to control
the absolute spectral gap of $\widetilde{\bm{M}}$. This inequality gives a
crucial perturbation bound for the singular values for general matrices. Since
the matrix $\widetilde{\bm{M}}$ can be seen as a perturbation (as we discussed
previously), we can (informally speaking) use Weyl’s inequality to control the
singular values. After some algebraic manipulation, we derive that
$\Gamma(\widetilde{\bm{M}})=\Omega(\lambda(\bm{Q}))$. ∎

So, a _single iteration_ of the parameterized CFTP algorithm with transition
matrix $\widetilde{\bm{M}}$ (see Line 5 of Algorithm 3) guarantees that with
an expected number of
$N=\widetilde{O}(nT_{\mathrm{mix}}(\widetilde{\bm{M}}))=\widetilde{O}\left(n\log(n)/\lambda(\bm{Q})\right)$
samples, it reaches Line 12 with a sample $y=F_{t}(1)$ that is generated by
the parameterized CFTP algorithm (the iterations in Line 4). We then have to
study the _quality of the output sample $y$_. Crucially, by the utility of the
standard CFTP algorithm, the law of $y$ is the normalized measure that is
proportional to $\mathcal{D}/\widetilde{\mathcal{D}}$, i.e.,
$\widetilde{\pi}_{0}(y)=(\mathcal{D}(y)/\widetilde{\mathcal{D}}(y))/\sum_{x\in[n]}\mathcal{D}(x)/\widetilde{\mathcal{D}}(x)$
(this is the stationary distribution of the modified random walk). Hence, we
cannot simply output $y.$ For this reason, we perform rejection sampling.
Specifically, as we will see, a single execution of Algorithm 2 either
produces a single sample distributed as in $\mathcal{D}$ with some known
probability $p$ or fails to output any sample with probability $1-p$ (this is
because we do not sample from $\mathcal{D}$ but from an almost uniform
distribution). If the rejection sampling process is successful, then the while
loop (Line 4) of Algorithm 3 will break and the clean sample will be given as
output. One can show that $\Theta(n)$ iterations of the parameterized CFTP
algorithm suffice to generate a perfect sample. Hence, the total sample
complexity is of order $O(n^{2}\log(n)\log(n/\delta)/\lambda(\bm{Q}))$ (at
each step $\delta^{\prime}=\Theta(\delta/n))$. We continue with a sketch
behind the rejection sampling process of Line 13 of Algorithm 2.

######  Claim 1 (Sample Quality and Rejection Sampling).

Algorithm 2, with input the output of the Learning Phase, outputs a state
$x\in[n]$ with probability proportional to $\mathcal{D}(x)$ or outputs
$\perp$. Moreover, Algorithm 3 outputs a perfect sample from $\mathcal{D}$
with probability 1.

To see this, set $A=\sum_{y\in[n]}\mathcal{D}(y)/\widetilde{D}(y)$. At the end
of Line 11, Algorithm 2 outputs $x\propto\mathcal{D}(x)/\widetilde{D}(x)$,
since the unique stationary distribution is
$\widetilde{\pi}(x)=(\mathcal{D}(x)/\widetilde{\mathcal{D}}(x))/A$ and the
utility of the CFTP guarantees that the outputs state follows the stationary
distribution. For the sample $x$, we perform a rejection sampling process,
with acceptance probability $\widetilde{\mathcal{D}}(x)$. Algorithm 2 has
$n+1$ potential outputs: It either prints $x\in[n]$ or $\perp$ indicating
failure. The arbitrary sample $x\in[n]$ is observed with probability
$\widetilde{\pi}(x)\cdot\widetilde{\mathcal{D}}(x)=\mathcal{D}(x)/A$. The
remaining probability mass is assigned to $\perp.$ We claim that the output of
Line 6 of Algorithm 3 has law $\mathcal{D}$. Observe that the whole stochastic
process of Algorithm 3 outputs $x\in[n]$ with probability
$\sum_{i=0}^{\infty}\Pr[\mathrm{Reject}]^{i}\Pr[x\textnormal{ is
Accepted}]=\mathcal{D}(x)/A\cdot\sum_{i=0}^{\infty}\Pr[\mathrm{Reject}]^{i}$,
where
$\Pr[\mathrm{Reject}]=\frac{1}{A}\sum_{y\in[n]}\frac{\mathcal{D}(y)}{\widetilde{\mathcal{D}}(y)}(1-\widetilde{\mathcal{D}}(y))=1-\frac{1}{A}\in(0,1)\,.$
As a result, we get that the probability that $x$ is the output of the above
stochastic process is $\mathcal{D}(x).$

##  4 Conclusion

We close with an open question: Is the bound
$\widetilde{O}(n^{2}/\lambda(\bm{Q}))$ tight? The main difficulty is that it
is not clear how to obtain a lower bound against _any possible_ exact sampling
algorithm. We would like to mention that the dependence on $1/\lambda(\bm{Q})$
is expected and intuitively unavoidable. This term is connected to the mixing
time. Let us become more specific. The (sample) complexity of our proposed
method is $\widetilde{O}(n^{2}/\lambda(\bm{Q}))$. The term $1/\lambda(\bm{Q})$
appears due to the mixing time of the (transformed) Markov chain. The
complexity of CFTP is $\widetilde{\Theta}(nT_{\mathrm{mix}})$. Hence, any
CFTP-based exact sampler (or more broadly random walk based algorithm) will
have that kind of spectral dependence.

Similarly, a dependence on the support’s size $n$ is necessary in general.
Note that the structure of the target distribution $\mathcal{D}$ is crucial:
if the target is uniform (or generally “flat”), then a bound of order
$\widetilde{O}(n/\lambda(\bm{Q}))$ is obtained by the naive CFTP method; for
distributions with modes, our algorithm achieves quadratic dependence on $n$.
The sample complexity of a potential (not necessarily random-walk based) exact
sampling algorithm will depend on $\mathcal{Q}$, but might achieve a sample
complexity not directly dependent on $\lambda(\bm{Q})$. The unclear part is
our algorithm’s dependence on $n$. While our algorithm attains a bound of
$\widetilde{O}(n^{2})$, it is not obvious whether an algorithm with sample
complexity $o(n^{2})$ exists.

Similarly, it is interesting to ask whether there are (efficient) ways to
transform the chain in order to get a faster sample corrector. Our learning
method finds a vector $\bm{p}$ (input of Algorithm 2) that makes the re-
weighted distribution close to uniform. Is it possible to find another
$\bm{p}$ (without our learning algorithm) such that (i) when transforming the
chain with $\bm{p}$, we obtain a fast mixing chain; and (ii) the success
probability in the rejection sampling step will be sufficiently high?

This work provided a theoretical understanding to perfect sampling from
pairwise comparisons. We believe it is an interesting direction to
experimentally evaluate our proposed methodology.

##### Acknowledgments.

We would like to thank the anonymous reviewers of NeurIPS 2022 for their
feedback and comments.

## References

  * [ACCL04] Nir Ailon, Bernard Chazelle, Seshadhri Comandur, and Ding Liu.  Property-preserving data reconstruction.  In International Symposium on Algorithms and Computation, pages 16–27. Springer, 2004.
  * [ACK15a] Jayadev Acharya, Clément L. Canonne, and Gautam Kamath.  Adaptive Estimation in Weighted Group Testing.  In Proceedings of the 2015 IEEE International Symposium on Information Theory, ISIT ’15, pages 2116–2120. IEEE Computer Society, 2015.
  * [ACK15b] Jayadev Acharya, Clément L. Canonne, and Gautam Kamath.  A Chasm Between Identity and Equivalence Testing with Conditional Queries.  In Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques., RANDOM ’15, pages 449–466, 2015.
  * [AJ21] Konrad Anand and Mark Jerrum.  Perfect sampling in infinite spin systems via strong spatial mixing.  arXiv preprint arXiv:2106.15992, 2021.
  * [APA18] Arpit Agarwal, Prathamesh Patil, and Shivani Agarwal.  Accelerated spectral ranking.  In International Conference on Machine Learning, pages 70–79. PMLR, 2018.  URL: https://proceedings.mlr.press/v80/agarwal18b/agarwal18b.pdf.
  * [BC18] Rishiraj Bhattacharyya and Sourav Chakraborty.  Property Testing of Joint Distributions using Conditional Samples.  Transactions on Computation Theory, 10(4):16:1–16:20, 2018.
  * [BC20] Siddharth Bhandari and Sayantan Chakraborty.  Improved bounds for perfect sampling of k-colorings in graphs.  In Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing, pages 631–642, 2020.
  * [BDNP21] Arnab Bhattacharyya, Rathin Desai, Sai Ganesh Nagarajan, and Ioannis Panageas.  Efficient statistics for sparse graphical models from truncated samples.  In International Conference on Artificial Intelligence and Statistics, pages 1450–1458. PMLR, 2021.
  * [BFFSZ19] Róbert Busa-Fekete, Dimitris Fotakis, Balázs Szörényi, and Manolis Zampetakis.  Optimal Learning of Mallows Block Model.  CoRR, abs/1906.01009, 2019.
  * [BGJ+10] Arnab Bhattacharyya, Elena Grigorescu, Madhav Jha, Kyomin Jung, Sofya Raskhodnikova, and David P Woodruff.  Lower bounds for local monotonicity reconstruction from transitive-closure spanners.  In Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques, pages 448–461. Springer, 2010.
  * [BLR93] Manuel Blum, Michael Luby, and Ronitt Rubinfeld.  Self-testing/correcting with applications to numerical problems.  Journal of computer and system sciences, 47(3):549–595, 1993.
  * [BT52] R.A. Bradley and M.E. Terry.  Rank Analysis of Incomplete Block Designs: I. The Method of Paired Comparisons.  Biometrika, 39:324, 1952.
  * [Can15] Clément L. Canonne.  Big Data on the Rise? - Testing Monotonicity of Distributions.  In Proceedings of the 42nd International Colloquium on Automata, Languages, and Programming, ICALP ’15, pages 294–305, 2015.
  * [Cat12] Manuela Cattelan.  Models for paired comparison data: A review with emphasis on dependent data.  Statistical Science, pages 412–433, 2012.
  * [CCK+19] Clément L Canonne, Xi Chen, Gautam Kamath, Amit Levi, and Erik Waingarten.  Random Restrictions of High-Dimensional Distributions and Uniformity Testing with Subcube Conditioning.  CoRR, abs/1911.07357, 2019.
  * [CGR18] Clément L Canonne, Themis Gouleakis, and Ronitt Rubinfeld.  Sampling correctors.  SIAM Journal on Computing, 47(4):1373–1423, 2018.
  * [DGTZ18] Constantinos Daskalakis, Themis Gouleakis, Christos Tzamos, and Manolis Zampetakis.  Efficient Statistics, in High Dimensions, from Truncated Samples.  In 59th Annual IEEE Symposium on Foundations of Computer Science (FOCS), pages 639–649. IEEE, 2018.
  * [DGTZ19] Constantinos Daskalakis, Themis Gouleakis, Christos Tzamos, and Manolis Zampetakis.  Computationally and Statistically Efficient Truncated Regression.  In Conference on Learning Theory (COLT), pages 955–960, 2019.
  * [DHKN17] Shaddin Dughmi, Jason D Hartline, Robert Kleinberg, and Rad Niazadeh.  Bernoulli factories and black-box reductions in mechanism design.  In Proceedings of the 49th Annual ACM SIGACT Symposium on Theory of Computing, pages 158–169, 2017.  URL: https://arxiv.org/pdf/1703.04143.pdf.
  * [DKNS01] Cynthia Dwork, Ravi Kumar, Moni Naor, and Dandapani Sivakumar.  Rank aggregation methods for the web.  In Proceedings of the 10th international conference on World Wide Web, pages 613–622, 2001.
  * [DKTZ21] Constantinos Daskalakis, Vasilis Kontonis, Christos Tzamos, and Emmanouil Zampetakis.  A statistical taylor theorem and extrapolation of truncated densities.  In Conference on Learning Theory, pages 1395–1398. PMLR, 2021.
  * [DRZ20] Constantinos Daskalakis, Dhruv Rohatgi, and Manolis Zampetakis.  Truncated linear regression in high dimensions.  arXiv preprint arXiv:2007.14539, 2020.
  * [FGY19] Weiming Feng, Heng Guo, and Yitong Yin.  Perfect sampling from spatial mixing.  arXiv preprint arXiv:1907.06033, 2019.
  * [FH00] James Allen Fill and Mark Huber.  The randomness recycler: a new technique for perfect sampling.  In Proceedings 41st annual symposium on foundations of computer science, pages 503–511. IEEE, 2000.
  * [FH10] Johannes Fürnkranz and Eyke Hüllermeier.  Preference learning and ranking by pairwise comparison.  In Preference learning, pages 65–82. Springer, 2010.
  * [FJO+15] Moein Falahatgar, Ashkan Jafarpour, Alon Orlitsky, Venkatadheeraj Pichapati, and Ananda Theertha Suresh.  Faster Algorithms for Testing under Conditional Sampling.  In Proceedings of the 28th Annual Conference on Learning Theory, COLT ’15, pages 607–636, 2015.
  * [FKKT21] Dimitris Fotakis, Alkis Kalavasis, Vasilis Kontonis, and Christos Tzamos.  Efficient algorithms for learning from coarse labels.  In Conference on Learning Theory, pages 2060–2079. PMLR, 2021.
  * [FKP22] Dimitris Fotakis, Alkis Kalavasis, and Eleni Psaroudaki.  Label ranking through nonparametric regression.  In Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 6622–6659. PMLR, 17–23 Jul 2022.
  * [FKS21] Dimitris Fotakis, Alkis Kalavasis, and Konstantinos Stavropoulos.  Aggregating incomplete and noisy rankings.  In International Conference on Artificial Intelligence and Statistics, pages 2278–2286. PMLR, 2021.
  * [FKT20] Dimitris Fotakis, Alkis Kalavasis, and Christos Tzamos.  Efficient parameter estimation of truncated boolean product distributions.  In Conference on Learning Theory, pages 1586–1600. PMLR, 2020.
  * [FV86] Michael A Fligner and Joseph S Verducci.  Distance Based Ranking Models.  Journal of the Royal Statistical Society. Series B (Methodological), pages 359–369, 1986.
  * [Gal97] Francis Galton.  An examination into the registered speeds of American trotting horses, with remarks on their value as hereditary data.  Proceedings of the Royal Society of London, 62(379-387):310–315, 1897.
  * [GJL19] Heng Guo, Mark Jerrum, and Jingcheng Liu.  Uniform sampling through the lovász local lemma.  66(3), 2019.
  * [GLR18a] Rong Ge, Holden Lee, and Andrej Risteski.  Beyond log-concavity: provable guarantees for sampling multi-modal distributions using simulated tempering langevin monte carlo.  In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pages 7858–7867, 2018.
  * [GLR18b] Rong Ge, Holden Lee, and Andrej Risteski.  Simulated tempering langevin monte carlo ii: An improved proof using soft markov chain decomposition.  arXiv preprint arXiv:1812.00793, 2018.
  * [GLT20] Evangelia Gergatsouli, Brendan Lucier, and Christos Tzamos.  Black-box methods for restoring monotonicity.  In International Conference on Machine Learning, pages 3463–3473. PMLR, 2020.
  * [GM99] Peter J Green and Duncan J Murdoch.  Exact sampling for bayesian inference: towards general purpose algorithms.  Bayesian statistics, 6:301–321, 1999.
  * [GTZ17] Themistoklis Gouleakis, Christos Tzamos, and Manolis Zampetakis.  Faster Sublinear Algorithms Using Conditional Sampling.  In Proceedings of the 28th Annual ACM-SIAM Symposium on Discrete Algorithms, SODA ’17, pages 1743–1757. SIAM, 2017.
  * [HJ12] Roger A. Horn and Charles R. Johnson.  Matrix analysis.  Cambridge university press, 2012.
  * [HOX14] Bruce Hajek, Sewoong Oh, and Jiaming Xu.  Minimax-optimal inference from partial rankings.  Advances in Neural Information Processing Systems, 27:1475–1483, 2014.
  * [HS00] Olle Häggström and Jeffrey E Steif.  Propp–wilson algorithms and finitary codings for high noise markov random fields.  Combinatorics, Probability and Computing, 9(5):425–439, 2000.
  * [HSW21] Kun He, Xiaoming Sun, and Kewen Wu.  Perfect sampling for (atomic) lov$\backslash$’asz local lemma.  arXiv preprint arXiv:2107.03932, 2021.
  * [Hub98] Mark Huber.  Exact sampling and approximate counting techniques.  In Proceedings of the thirtieth annual ACM symposium on Theory of computing, pages 31–40, 1998.
  * [Hub04] Mark Huber.  Perfect sampling using bounding chains.  The Annals of Applied Probability, 14(2):734–753, 2004.
  * [Hub16] Mark L. Huber.  Perfect simulation, volume 148.  CRC Press, 2016.
  * [Hun04] David R. Hunter.  Mm algorithms for generalized bradley-terry models.  The annals of statistics, 32(1):384–406, 2004.
  * [HW71] David Lee Hanson and Farroll Tim Wright.  A bound on tail probabilities for quadratic forms in independent random variables.  The Annals of Mathematical Statistics, 42(3):1079–1083, 1971.
  * [IZD20] Andrew Ilyas, Emmanouil Zampetakis, and Constantinos Daskalakis.  A theoretical and practical framework for regression and classification from truncated samples.  In International Conference on Artificial Intelligence and Statistics, pages 4463–4473. PMLR, 2020.
  * [Jer03] Mark Jerrum.  Counting, sampling and integrating: algorithms and complexity.  Springer Science & Business Media, 2003.
  * [JLYY11] Xiaoye Jiang, Lek-Heng Lim, Yuan Yao, and Yinyu Ye.  Statistical ranking and combinatorial hodge theory.  Mathematical Programming, 127(1):203–244, 2011.
  * [JR13] Madhav Jha and Sofya Raskhodnikova.  Testing and reconstruction of lipschitz functions with applications to data privacy.  SIAM Journal on Computing, 42(2):700–731, 2013.
  * [JSS21] Vishesh Jain, Ashwin Sah, and Mehtaab Sawhney.  Perfectly sampling $k\geq(8/3+o(1))\delta$-colorings in graphs.  In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing, pages 1589–1600, 2021.
  * [JSV04] Mark Jerrum, Alistair Sinclair, and Eric Vigoda.  A polynomial-time approximation algorithm for the permanent of a matrix with nonnegative entries.  Journal of the ACM (JACM), 51(4):671–697, 2004.
  * [Kaz11] Gabriella Kazai.  In search of quality in crowdsourcing for search engine evaluation.  In European conference on information retrieval, pages 165–176. Springer, 2011.
  * [Ken05] Wilfrid S. Kendall.  Notes on perfect simulation.  Markov chain Monte Carlo: innovations and applications, 7, 2005\.
  * [KT19] Gautam Kamath and Christos Tzamos.  Anaconda: A Non-Adaptive Conditional Sampling Algorithm for Distribution Testing.  In Proceedings of the Thirtieth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 679–693. SIAM, 2019.
  * [KTZ19] Vasilis Kontonis, Christos Tzamos, and Manolis Zampetakis.  Efficient Truncated Statistics with Unknown Truncation.  In 260th Annual IEEE Symposium on Foundations of Computer Science (FOCS), pages 1578–1595. IEEE, 2019.
  * [LM18] Allen Liu and Ankur Moitra.  Efficiently learning mixtures of mallows models.  In 2018 IEEE 59th Annual Symposium on Foundations of Computer Science (FOCS), pages 627–638. IEEE, 2018.
  * [LOAF12] Miguel Angel Luengo-Oroz, Asier Arranz, and John Frean.  Crowdsourcing malaria parasite quantification: an online game for analyzing images of infected thick blood smears.  Journal of medical Internet research, 14(6):e2338, 2012.
  * [LP17] David A. Levin and Yuval Peres.  Markov chains and mixing times, volume 107.  American Mathematical Soc., 2017.  URL: https://www.stat.berkeley.edu/users/aldous/260-FMIE/Levin-Peres-Wilmer.pdf.
  * [LSR21] Wanshan Li, Shamindra Shrotriya, and Alessandro Rinaldo.  The performance of the mle in the bradley-terry-luce model in $\ell_{\infty}$-loss and under general graph topologies.  arXiv preprint arXiv:2110.10825, 2021.  URL: https://arxiv.org/abs/2110.10825.
  * [Luc59] R.D. Luce.  Individual Choice Behavior.  Wiley, 1959.
  * [Mar96] John I. Marden.  Analyzing and modeling rank data.  CRC Press, 1996.
  * [MMS21] Ankur Moitra, Elchanan Mossel, and Colin P Sandon.  Learning to sample from censored markov random fields.  In Mikhail Belkin and Samory Kpotufe, editors, Proceedings of Thirty Fourth Conference on Learning Theory, volume 134 of Proceedings of Machine Learning Research, pages 3419–3451. PMLR, 15–19 Aug 2021.
  * [MP92] Enzo Marinari and Giorgio Parisi.  Simulated tempering: a new monte carlo scheme.  EPL (Europhysics Letters), 19(6):451, 1992.
  * [MP05] Elchanan Mossel and Yuval Peres.  New coins from old: computing with unknown bias.  Combinatorica, 25(6):707–724, 2005.
  * [NOS12] Sahand Negahban, Sewoong Oh, and Devavrat Shah.  Iterative ranking from pair-wise comparisons.  In Advances in neural information processing systems, pages 2474–2482, 2012.
  * [NOS17] Sahand Negahban, Sewoong Oh, and Devavrat Shah.  Rank centrality: Ranking from pairwise comparisons.  Operations Research, 65(1):266–287, 2017.
  * [NP05] Şerban Nacu and Yuval Peres.  Fast simulation of new coins from old.  Annals of Applied Probability, 15(1A):93–115, 2005.
  * [NP19] Sai Ganesh Nagarajan and Ioannis Panageas.  On the Analysis of EM for truncated mixtures of two Gaussians.  In 31st International Conference on Algorithmic Learning Theory (ALT), pages 955–960, 2019.
  * [Ple21] Orestis Plevrakis.  Learning from censored and dependent data: The case of linear dynamics.  In Mikhail Belkin and Samory Kpotufe, editors, Proceedings of Thirty Fourth Conference on Learning Theory, volume 134 of Proceedings of Machine Learning Research, pages 3771–3787. PMLR, 15–19 Aug 2021.
  * [PW96] James Gary Propp and David Bruce Wilson.  Exact sampling with coupled markov chains and applications to statistical mechanics.  Random Structures & Algorithms, 9(1-2):223–252, 1996.  URL: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.46.3607&rep=rep1&type=pdf.
  * [PW98] James Gary Propp and David Bruce Wilson.  How to get a perfectly random sample from a generic markov chain and generate a random spanning tree of a directed graph.  Journal of Algorithms, 27(2):170–217, 1998.
  * [RA14] Arun Rajkumar and Shivani Agarwal.  A statistical convergence perspective of algorithms for rank aggregation from pairwise data.  In International Conference on Machine Learning, pages 118–126, 2014.
  * [RV13] Mark Rudelson and Roman Vershynin.  Hanson-wright inequality and sub-gaussian concentration.  Electronic Communications in Probability, 18, 2013.
  * [SBB+16] Nihar B Shah, Sivaraman Balakrishnan, Joseph Bradley, Abhay Parekh, Kannan Ramchandran, and Martin J Wainwright.  Estimation from pairwise comparisons: Sharp minimax bounds with topology dependence.  The Journal of Machine Learning Research, 17(1):2049–2095, 2016\.  URL: https://jmlr.org/papers/volume17/15-189/15-189.pdf.
  * [SBGW16] Nihar Shah, Sivaraman Balakrishnan, Aditya Guntuboyina, and Martin Wainwright.  Stochastically transitive models for pairwise comparisons: Statistical and computational issues.  In International Conference on Machine Learning, pages 11–20. PMLR, 2016.
  * [SS10] Michael Saks and Comandur Seshadhri.  Local monotonicity reconstruction.  SIAM Journal on Computing, 39(7):2897–2926, 2010.
  * [Tao12] Terence Tao.  Topics in random matrix theory, volume 132.  American Mathematical Soc., 2012.  URL: http://www.stat.ucla.edu/~ywu/research/documents/BOOKS/TaoRandomMatrix.pdf.
  * [Thu27] Louis L Thurstone.  A law of comparative judgment.  Psychological review, 34(4):273, 1927.
  * [Tro15] Joel A Tropp.  An introduction to matrix concentration inequalities.  Foundations and Trends® in Machine Learning, 8(1-2):1–230, 2015.  URL: https://arxiv.org/pdf/1501.01571.pdf.
  * [VY16] Milan Vojnovic and Seyoung Yun.  Parameter estimation for generalized thurstone choice models.  In International Conference on Machine Learning, pages 498–506, 2016.
  * [VYZ20] Milan Vojnovic, Se-Young Yun, and Kaifang Zhou.  Convergence rates of gradient descent and mm algorithms for bradley-terry models.  In International Conference on Artificial Intelligence and Statistics, pages 1254–1264. PMLR, 2020.  URL: https://arxiv.org/abs/1901.00150.
  * [Wil00] David Bruce Wilson.  How to couple from the past using a read-once source of randomness.  Random Structures & Algorithms, 16(1):85–113, 2000.
  * [WJJ13] Fabian Wauthier, Michael Jordan, and Nebojsa Jojic.  Efficient ranking from pairwise comparisons.  In International Conference on Machine Learning, pages 109–117. PMLR, 2013.
  * [ZCKM20] Tom Zahavy, Alon Cohen, Haim Kaplan, and Yishay Mansour.  Unknown mixing times in apprenticeship and reinforcement learning.  In Jonas Peters and David Sontag, editors, Proceedings of the 36th Conference on Uncertainty in Artificial Intelligence (UAI), volume 124 of Proceedings of Machine Learning Research, pages 430–439. PMLR, 03–06 Aug 2020.  URL: https://proceedings.mlr.press/v124/zahavy20a.html.

##  Appendix A Preliminaries

In this section, we provide some preliminaries and some useful tools about (i)
concentration of random matrices, (ii) random walks and (iii) CFTP. The reader
may skip this section.

####  A.0.1 Random Matrices

We continue with some definitions for random matrices, needed for the proof of
Lemma 17. The following can be found at [Tro15].

######  Definition 6.

Let $(\Omega,\mathcal{F},\mu)$ be a probability space. A random matrix
$\bm{X}$ is a measurable map $\bm{X}:\Omega\rightarrow\mathcal{M}^{n_{1}\times
n_{2}}$.

The entries of $\bm{X}$ may be considered complex random variables that may or
may not be correlated with each other. A finite sequence $\\{\bm{X}_{i}\\}$ of
random matrices is independent when

$\mu(\bm{X}_{k}\in F_{k}\text{ for each }k)=\prod_{k}\mu(\bm{X}_{k}\in
F_{k})\,,$

for any collection $\\{F_{k}\\}$ of Borel subsets of $\mathcal{M}^{n_{1}\times
n_{2}}$.

######  Proposition 7 (Hermitian Matrix Chernoff Bounds (see [Tro15])).

Consider a finite sequence $\\{\bm{X}_{i}\\}$ of $m$ independent, random,
Hermitian square matrices with common dimension $n$. Assume that
$0\leq\lambda_{\min}(\bm{X}_{i})\leq\lambda_{\max}(\bm{X}_{i})\leq M$ for any
$i\in[m]$. Let $\bm{Y}=\sum_{i}X_{i}$. Then, for $\epsilon\in[0,1)$:

$\Pr\Big{[}\lambda_{\min}(\bm{Y})\leq(1-\epsilon)\lambda_{\min}(\mathbb{E}\bm{Y})\Big{]}\leq
n(e^{-\epsilon}/(1-\epsilon)^{1-\epsilon})^{\lambda_{\min}(\mathbb{E}\bm{Y})/M}\leq
n\exp\Big{(}-\epsilon^{2}\lambda_{\min}(\mathbb{E}\bm{Y})/(2M)\Big{)}\,.$

####  A.0.2 Random Walks and Markov Chains

This section is mostly based on [LP17] and we refer the reader to [LP17] for a
thorough exposition.

##### Markov Chains.

Let $\Omega$ be a finite state space. A Markov chain is a sequence of random
variables $X_{0},X1,\dots$ that satisfy the Markov property, i.e.,

$\Pr[X_{t+1}=x|X_{0}=x_{0},\dots,X_{t}=x_{t}]=\Pr[X_{t+1}=x|X_{t}=x_{t}]\,.$

A Markov chain is called time-homogeneous, if the RHS of the above equation
does not depend of $t$. Such a chain is associated with a transition matrix
$\bm{P}=\\{\bm{P}(x,y)\\}$, where $(x,y)\in\Omega\times\Omega$. It holds that

$\bm{P}(x,y)=\Pr[X_{t+1}=y|X_{t}=x]\text{ for all
}x,y\in\Omega,t\in\mathbb{N}\,.$

A Markov chain is _ergodic_ if there exists a time $t\in\mathbb{N}$ such that
$\bm{P}^{(t)}(x,y)>0$ for any $x,y\in\Omega$, i.e., there exists a finite time
$t$ so that the probability of going from any vertex to any other in $t$ steps
is positive. For finite state Markov chains, ergodicity is equivalent to
irreducibility and aperiodicity. A Markov chain is _irreducible_ if for any
two states $x,y\in\Omega$, there exists a time step $t$ such that
$\bm{P}^{(t)}(x,y)>0$ and is _aperiodic_ if, for any state $x$, it holds that
$\mathrm{gcd}\\{t|\bm{P}^{(t)}(x,x)>0\\}=1$. Observe that a Markov chain with
self-loops is aperiodic.

##### Stationary Distribution.

A stationary distribution $\bm{\pi}\in\Delta^{n}$, i.e., the
$(n-1)$-dimensional simplex whose vertices are the $n$ standard unit vectors,
is defined as the fixed point of the transition matrix $\bm{P}$, that is
$\bm{\pi}^{T}\bm{P}=\bm{\pi}^{T}$. An ergodic Markov chain has a _unique
stationary distribution_ $\bm{\pi}$ and, as $t$ increases, it converges to it.
An interesting property of Markov chains is the reversibility. A Markov chain
is _reversible_ if there exists a distribution $\bm{\pi}$ that satisfies the
detailed balanced equations:

$\bm{\pi}(x)\bm{P}(x,y)=\bm{\pi}(y)\bm{P}(y,x)\text{ for all }x,y\in\Omega\,.$
(A.1)

In this case, we can verify that $\bm{\pi}$ is a stationary distribution. We
can simply write $\pi$ for the stationary distribution with associated
probability vector $\bm{\pi}$.

##### Mixing Time.

It is important to understand the convergence time of a Markov chain to its
stationary distribution $\pi$. A crucial random variable for convergence time
is the _mixing time_.

######  Definition 8 (Mixing Time).

Let $0<\epsilon<1/2$. Let $M$ be an ergodic Markov chain on a finite state
space $\Omega$ with stationary distribution $\pi$. Then, the mixing time with
accuracy $\epsilon$ of $M$ equals:

$T_{\mathrm{mix}}(\bm{P};\epsilon)=\inf\\{t>0:\max_{x\in\Omega}\mathrm{TV}(\bm{P}_{x}(X_{t}),\pi)\leq\epsilon\\}\,,$

where $\bm{P}_{x}(X_{t})$ is the distribution of the state $X_{t}$ at time $t$
for starting state $x\in\Omega$.

For a reversible transition matrix $\bm{P}$, let its spectrum be:

$1=\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{n}\geq-1\,.$

Note that $\lambda_{2}<1$ if and only if the chain is irreducible (exactly one
connected component) and $\lambda_{n}>-1$ if and only if the chain is
aperiodic (e.g., not a bipartite graph). We define the _absolute spectral gap_
of $\bm{P}$ to be the difference:
$\Gamma(\bm{P})=1-\max\\{|\lambda_{2}(\bm{P})|,|\lambda_{n}(\bm{P})|\\}$. It
holds that if $\bm{P}$ is aperiodic and irreducible, then $\Gamma(\bm{P})$ is
strictly positive. One could also define the spectral gap
$\lambda(\bm{P})=1-\lambda_{2}(\bm{P})$. The two gaps are equal when the chain
is lazy, i.e., when for any state $x\in\Omega,$ it holds that $\bm{P}(x,x)\geq
1/2$.

######  Lemma 9 (Bounding $T_{\mathrm{mix}}$, see [LP17]).

Let $0<\epsilon<1$. Assume that the transition matrix $\bm{P}$ is aperiodic,
irreducible and reversible with respect to $\pi$. Then, it holds that:

$(t_{\mathrm{rel}}-1)\log\Big{(}\frac{1}{2\epsilon}\Big{)}\leq
T_{\mathrm{mix}}(\bm{P};\epsilon)\leq\log\Big{(}\frac{1}{\epsilon\pi_{\mathrm{min}}}\Big{)}t_{\mathrm{rel}}\,,$

where $t_{\mathrm{rel}}=1/\Gamma(\bm{P})$ is the relaxation time of the Markov
chain, i.e., the inverse absolute spectral gap.

##### Coalescencing Random Walks.

Let $\Omega=[n]$. In a coalescing random walk, a set of $n$ particles perform
independent discrete time random walks on an undirected connected graph
$\mathcal{G}=(V,E)$ with $|V|=n$, with each particle initially placed at a
single (distinct) vertex $x\in V$. In each time step, all particles move
simultaneously. Whenever two or more particles meet at a vertex, they unite to
form a single particle which then continues to make a random walk through the
graph. The _coalescence time_ $T_{\mathrm{coal}}$ is a random variable and is
the first time when all particles coalesce. More formally, we can define
_coalescence_ as follows.

######  Definition 10 (Coalescence of Stochastic processes (see [Hub16])).

Let $\mathcal{A}$ be a collection of stochastic processes, defined over a
common index set $\mathcal{I}$ and common state space $\Omega$. If there is an
index $i\in\mathcal{I}$ and state $x\in\Omega$ such that, for all stochastic
processes $X\in\mathcal{A}$, it holds that $X_{i}=x$, then we say that the
stochastic processes have coalesced.

##### Coalescence Time of Random Walks.

Consider the random variable $X^{(i)}_{t}$, that indicates the position of the
$i$-th particle at time $t$. The coalescence time is equal to:

$T_{\mathrm{coal}}=\inf_{t>0}\Big{\\{}X_{t}^{(i)}=X_{t}^{(j)}\text{ for any
}i\neq j,i,j\in[n]\Big{\\}}\,.$

We conclude this section with a folklore result, that deals with the possible
locations in the complex plane of the eigenvalues of a square matrix
$\bm{A}\in\mathcal{M}_{n}$.

######  Lemma 11 (Geršgorin’s Theorem (see [HJ12])).

Let $\bm{A}=[\bm{A}_{ij}]\in\mathcal{M}_{n}$. For any $i\in[n]$, let

$R_{i}(\bm{A})=\sum_{j\neq i}|\bm{A}_{ij}|\,,$

and consider the $n$ Geršgorin disks

$B_{i}=\\{z\in\mathbb{C}:|z-\bm{A}_{ii}|\leq R_{i}(\bm{A})\\}\,.$

Then, the eigenvalues of $\bm{A}$ are in the union of the Geršgorin disks,
i.e.,

$\mathrm{spec}(\bm{A})\subseteq\bigcup_{i=1}^{n}B_{i}\,.$

####  A.0.3 Exact Sampling and CFTP Preliminaries

Markov chain Monte Carlo (MCMC) methods constitute a class of algorithms for
sampling from probability measures and arise naturally in various fields of
science such as theoretical computer science (e.g., approximation algorithms
for #$\mathrm{P}$-complete problems (see [Jer03, JSV04])), statistical physics
(e.g., in order to understand phase transition phenomena for Ising models) and
statistics. However, the theory of MCMC, in terms of exact sampling, can be
seen as an asymptotic analysis (i.e., in the limit, the total variation
distance vanishes). Perfect simulation is analogous to MCMC and deals with
techniques for designing algorithms that return exact draws from the target
distribution, instead of long-time approximations. Exact sampling comprises a
well-studied field (see [Ken05, Hub16]) and has numerous applications in
computer science (e.g., in approximate counting [Hub98]). The importance of
exact simulation gave rise to various procedures in order to generate perfect
samples. See [PW96, GM99, Wil00, HS00, FH00, Hub04] for a small sample of this
line of research.

Figure 2: Execution of the CFTP Algorithm for a random walk on the path graph.
Observe that it suffices to execute the random procedure only for the two
extreme points, since their trajectories dominate any intermediate point of
the path. After a series of executions which did not coalesce (see the right
part of Figure 2), we observe that the two stochastic evolutions, the blue and
the red one in the left part of Figure 2, coalesce and continue as a single
trajectory until $t=0$ (see the green trajectory). A sample execution code can
be found at https://warwick.ac.uk/fac/sci/statistics/staff/academic-
research/kendall/personal/perfect_programs/

Coupling from the past (CFTP) is a technique developed by Propp and Wilson
[PW96, PW98], that provides an exact random sample from a Markov chain, that
has law the (unique) stationary distribution. The algorithm assumes that a
particle has been running on the Markov chain since time $-\infty$
(arbitrarily long in the past) and we are concerned with the location of the
particle at (the fixed) time $t=0.$ The fact that the stopping time is
deterministic and not a random variable (as in MCMC) is crucial. Since the
particle performs a random walk in the Markov chain infinitely long, one would
intuitively believe that the particle at time $t=0$ is distributed according
to the stationary distribution.

Algorithm 4 Coupling From The Past Algorithm

1:procedure CouplingFromThePast()$\triangleright$ _Assuming (adaptive) access
to the Markov oracle_.

2: $t\leftarrow 0$

3: $F_{(t,0]}(i)\leftarrow i,\text{ for any }i\in[n]$

4: while $F_{(t,0]}$ has not coalesced do

5: $t\leftarrow t-1$

6: $f_{t}\leftarrow\mathrm{Markov}()$

7: $F_{(t,0]}\leftarrow F_{(t+1,0]}\circ f_{t}$

8: end

9: return $F_{(t,0]}(1)$

We define $F_{(t,0]}$ as the $t$-step evolution of the Markov chain from time
$t\in\mathbb{Z}_{\leq 0}$ to the fixed time $T=0$ (evolution _from the past_).
Note that we can decompose the evolution of the chain into $t$ independent
applications of the random function $f$, which encodes the information of the
random walk, i.e., $\Pr[f(x)=y]=\bm{P}_{xy}$ for any pair of states
$x,y\in\Omega$. Hence, we have that

$F_{(t,0]}=f_{-1}\circ f_{-2}\circ\ldots\circ f_{t}\,.$

We can think of the data structure $F_{(t,0]}$ as a list of $n$ values, each
one capturing the evolution of the Markov chain, initiated in a different
state at time $t$ and stopped at time $0$. Similarly, the oracle call
$\mathrm{Markov}()$ corresponds to $n$ adaptive calls, each one giving a
transition from a fixed state $x\in\Omega$, where the next state is
distributed according to $\bm{P}(x,\cdot)$. Hence, $f_{t}$ is an
$n$-dimensional vector with values in $\Omega$. We refer the reader to Figure
2 for an execution of the CFTP algorithm.

The following fact is the main result behind CFTP.

######  Fact 1 ([PW96]).

With probability $1$, the Coupling from the Past algorithm returns a value,
and this value is distributed according to the stationary distribution of the
Markov chain.

The following proposition provides an upper bound on the convergence of the
CFTP algorithm.

######  Proposition 12 ([ZCKM20]).

Let $\pi$ be the stationary distribution of an ergodic Markov chain $\bm{P}$
with $n$ states. We run $n$ simulations of the chain each starting at a
different state. When two or more simulations coalesce, we merge them into a
single simulation. With probability at least $1-\delta$, all $n$ chains have
merged after at most $O(nT_{\mathrm{mix}}(\bm{P};1/4)\log(1/\delta))$
iterations.

##  Appendix B The Proof of Theorem 2 (Direct Exact Sampling from LSS)

In this section, we provide the proof of Theorem 2. Let us begin with the
algorithm. We remark that in the proof we will use the more technical notation
$F_{(a,b]}$ to denote the evolution of the chain from $a$ to $b$. In Algorithm
1, this is abbreviated as $F_{t}$ to denote $F_{(t,0]}$.

###### Proof.

At each iteration, the algorithm draws a single sample from
$\mathrm{Samp}(\mathcal{Q};\mathcal{D})$ Then, the algorithm performs the
following procedure for each state $v\in[n]$: let the drawn sample be
$(i,j,w),$ where $(i,j)\in\mathcal{E}$ and $w\in\\{0,1\\}$ indicating the
winning node between $i$ and $j$, i.e., the node $i$ stays in $i$ (if $w=1$)
or node $i$ moves to $j$ (if $w=0$). For any $v\in[n]$, the algorithm checks
whether $v\in\\{i,j\\}$. If not, the state $v$ remains at $v$; otherwise,
assume that $i=v$. Then, the transition $v\to j$ is implemented only if $w=0$,
otherwise it remains at $v$. We claim that this process simulates the
transitions of the matrix $\bm{M}$ and has unique stationary distribution
$\mathcal{D}$. The probability that the transition $v\to j$ occurs is

$\displaystyle\Pr[v\to j]$ $\displaystyle=\Pr[(v,j,w)\text{ is drawn and
}w=0]=\Pr[(v,j)\text{ is drawn}]\cdot\Pr[w=0|(v,j)\text{ is drawn}]$
$\displaystyle=\mathcal{Q}(v,j)\cdot\frac{\mathcal{D}(j)}{\mathcal{D}(v)+\mathcal{D}(j)}=\bm{M}_{vj}\,.$

Hence, for $v\in[n]$ and by the above process, the transition probabilities
are given by the matrix $\bm{M}$, that is a transition matrix with unique
stationary distribution $\mathcal{D}$. We will see that this random process
gives a perfect sample.

Consider the state space $\Omega=[n]$ and let $T$ be a positive integer. Let
$(X_{t})_{t\in\mathbb{Z}}$ be a random walk on $[n]$ with transition
probabilities as in $\bm{M}$ (see Equation (2.1)). Let
$F_{(a,b]}:[n]\rightarrow[n]$ be such that $F_{(-t_{0},T]}(x_{0})$ denotes the
state of the walk at time $t=T$, i.e., $X_{T}\in[n]$, when the random walk
begun at time $t=-t_{0}$ (i.e., _in the past_) with initial state
$X_{-t_{0}}=x_{0}$. Specifically, for each $t_{0},T$, we have that
$F_{(-t_{0},T]}:[n]\rightarrow[n]$ is a random map and $F_{(-t_{0},T]}(x_{0})$
is a realization of the Markov chain, with initial state $x_{0}$, executed
from time $-t_{0}$ to time $T$. We call $F_{(-t_{0},T]}$ a _stochastic flow_.

The CFTP method is initialized at time $t=-t_{0}$ into the past and places $n$
particles (i.e., simulates $n$ concurrent random walks), one at each state
$x\in[n]$. CFTP terminates when the stochastic flow $F_{(-t_{0},0]}$ is
constant (it has a single point random image, i.e.,
$F_{(-t_{0},0]}(x)=F_{(-t_{0},0]}(y)$, for all $x,y$). Namely, CFTP terminates
when all the random walks that start at time $t=-t_{0}$ have coalesced at time
$t=0$ at a single $z\in[n]$. Then, the algorithm outputs the (common) state
$z=F_{(-t_{0},0]}(x)$, for any $x\in[n]$.

Since the transitions of each random walk are performed with respect to
$\bm{M}$, the resulting sample $z$ is distributed according to the stationary
distribution $\mathcal{D}$. This can be shown as follows (see also [PW96,
Hub16]): Let $\mathcal{L}(X)$ denote the distribution of the random variable
$X$. Since the chain is aperiodic and irreducible, coalescence occurs almost
surely. Specifically, assume that coalescence occurs when the simulation
starts at a finite time $t=-t_{0}$. First, observe that if coalescence occurs
at time $-t_{0}$, then if CFTP starts at any time $-t<-t_{0}$, we end up in
the same state (because the randomness in $[-t_{0},0]$ remains the same),
i.e., $F_{(-t,0]}(x)=F_{(-t_{0},0]}(x)$ for any $-t<-t_{0}$ and any state
$x\in[n]$. Since, coalescence has occurred, we can work with a fixed starting
state $x_{0}$ and let $F_{(-t,0]}=F_{(-t,0]}(x_{0})$. Hence, it holds that
$\mathcal{L}(F_{(-t,0]})=\mathcal{L}(F_{(-t_{0},0]})$ and so

$\mathrm{TV}(\mathcal{L}(F_{(-t_{0},0]}),\mathcal{D})=\lim_{t\rightarrow\infty}\mathrm{TV}(\mathcal{L}(F_{(-t,0]}),\mathcal{D})\,.$

We note that the distribution of the time-backward coalescing random walk is
equal to the distribution of the time-forward random walk, i.e.,
$\mathcal{L}(F_{(-t,0]})=\mathcal{L}(F_{(0,t]})$. Thus,
$\mathrm{TV}(\mathcal{L}(F_{(-t_{0},0]}),\mathcal{D})=\lim_{t\rightarrow\infty}\mathrm{TV}(\mathcal{L}(F_{(0,t]}),\mathcal{D})=0$.
So, the output sample of the CFTP algorithm is distributed according to
$\mathcal{D}$. Using Proposition 12, we get that, with high probability,
coalescence occurs after $O(nT_{\mathrm{mix}}(\bm{M};1/4))$ iterations. From
Lemma 9 with $\epsilon=1/4$, we get that the mixing time of the Markov chain
is $O(\log(1/\mathcal{D}_{\mathrm{min}})/\Gamma(\bm{M}))$. Finally, note that
in order to perform a single iteration of the CFTP method, it requires to draw
a single sample from the Local Sampling Scheme oracle. This concludes the
analysis of the sample complexity. The time complexity is also polynomial in
the number of samples [PW96]. ∎

######  Remark 1.

As we mentioned in the discussion after 1, in order to obtain Theorem 2 only 1
is required. Under only this assumption, a quite large amount of Local
Sampling Schemes is allowed and the resulting sample complexity is
$\widetilde{O}\left(\frac{n\log(1/\mathcal{D}_{\mathrm{min}})}{\Gamma(\bm{M})}\right)$.
However, if we additionally assume that 2 holds, as we did in the statement of
1, we accept only LSSs whose target distribution satisfies
$\mathcal{D}_{\mathrm{min}}=\Omega(2^{-n}).$ This yields the (non-tight) bound
$\widetilde{O}(n^{2}/\Gamma(\bm{M}))$ of 1. We preferred to omit this detail
in 1 and emphasize only on the exponential dependence due to $\Gamma(\bm{M})$
for many natural instances.

##  Appendix C The Proof of Theorem 3 (Exact Sampling using Learning)

The main result of this section is the convergence analysis of the Coupling
From the Past algorithm when the algorithm performs the rescaling
transformation internally in each transition.

######  Theorem 13 (Convergence of CFTP).

Assume that 1 and 2 both hold. Algorithm 2, assuming access to a Local
Sampling Scheme $\mathrm{Samp}(\mathcal{Q};\mathcal{D})$ with weighted
Laplacian matrix $\bm{Q}$ (see Equation (1.1)) associated with the pair
distribution $\mathcal{Q}$ over $[n]\times[n]$, satisfies the following:

  * $(i)$

Algorithm 2 reaches Line 12 with probability $1$.

  * $(ii)$

For any $\delta\in(0,1)$, a single execution of of Algorithm 2 uses
$\widetilde{O}\left(\frac{n\log(n)\log(1/\delta)}{\lambda(\bm{Q})}\right)$
samples with probability $1-\delta$ and the running time is polynomial in the
number of samples.

  * $(iii)$

Algorithm 2, with input the output of the Learning Phase, outputs a state
$x\in[n]$ with probability proportional to $\mathcal{D}(x)$ or outputs
$\perp$. Moreover, Algorithm 3 outputs a perfect sample from $\mathcal{D}$.

Using the above properties, for any $\delta>0$, Algorithm 3 draws, with
probability at least $1-\delta$,
$N={O}\left(n^{2}\log^{2}(n)/\lambda(\bm{Q})\cdot\log(1/\delta)\right)$
samples from a Local Sampling Scheme $\mathrm{Samp}(\mathcal{Q};\mathcal{D})$,
runs in time polynomial in $N$, and outputs a sample distributed as in
$\mathcal{D}$.

###### Proof _of_ Theorem 13.

Consider a target distribution $\mathcal{D}$, supported on $[n]$ and sample
access to a Local Sampling Scheme $\mathrm{Samp}(\mathcal{Q};\mathcal{D})$ for
a pair distribution $\mathcal{Q}$. Let $\bm{M}=\bm{Q}\otimes\bm{D}$ be the
transition matrix of the Markov chain associated with the Local Sampling
Scheme graph. From Theorem 4 with accuracy $\epsilon>0$, there exists an
algorithm that, with high probability, computes the rescaling of the Markov
chain that corresponds to $\bm{M}$, which we denote by $\widetilde{\bm{M}}$,
the transition matrix of the rescaled (by $1/\widetilde{\mathcal{D}}$) Markov
chain. The learning algorithm uses $\widetilde{O}(n^{2}/(\lambda(\bm{Q})))$
with probability $1-\delta$ and is executed only once. We now analyze a single
execution of Algorithm 2.

######  Claim 2 (Termination of CFTP).

Algorithm 2 reaches Line 12 with probability $1$.

###### Proof.

Since the chain $\widetilde{\bm{M}}$ is ergodic, there is a $T$ such that for
any pair of states $(x,y)$, the probability $\widetilde{\bm{M}}^{T}(x,y)>0$.
Similarly to the proof sketch of the naive exact sampling algorithm (see also
Theorem 2), let $(X_{t})_{t\in\mathbb{Z}}$ be a random walk on $[n]$ and
define the mapping $F_{(a,b]}:[n]\rightarrow[n]$, so that
$F_{(-t_{0},T]}(x_{0})$ is the state of the random walk at time $t=T$, i.e.,
$X_{T}\in[n]$, where the random walk begun at time $t=-t_{0}$ (in the past)
with starting state $X_{-t_{0}}=x_{0}$. Then, the $n$ random walks have
coalesced after $T$ steps if and only if $|\mathrm{Image}(F_{(-T,0]})|=1$. By
the ergodicity of the chain, each one of the events
$E_{T}=\\{|\mathrm{Image}(F_{(-T,0]})|=1\\}$ for $T\in\mathbb{Z}_{\geq 0}$ has
strictly positive probability to occur. Since the events are independent, with
probability $1$, there will be an event $E_{T^{\star}}$ that occurs. Then, for
any $T>T^{\star},$ the desired property holds. ∎

######  Claim 3 (Sample Complexity of a Single Iteration).

It holds that each execution of Algorithm 2 uses
$O\left(\frac{n\log(n)\log(1/\delta)}{\lambda(\bm{Q})}\right)$ samples from
the Local Sampling Scheme with probability $1-\delta$.

###### Proof.

The CFTP draws a single sample for each iteration and performs a step
according to the matrix $\widetilde{\bm{M}}$ due to the rescaling. First, the
unique stationary distribution is almost uniform (Lemma 15) and so
$\log(1/\mathcal{D}_{\mathrm{min}})=\Theta(\log(n))$. This implies, by Lemma 9
and Lemma 15, that the mixing time of $\widetilde{\bm{M}}$ for some constant
accuracy $\epsilon_{0}=1/4$ is
$T_{\mathrm{mix}}(\widetilde{\bm{M}};1/4)=O(\log(n)/\lambda(\bm{Q}))$. Hence,
by Proposition 12, we get the desired result. ∎

Moreover, the update time of each execution of the CFTP algorithm is
polynomial in $n$. Note that, since during the Learning Phase, the algorithm
uses $O(n^{2}/\lambda(\bm{Q}))$ and is executed once, the sample complexity of
Algorithm 3 follows, since Algorithm 2 will be executed $\Theta(n)$ times.

######  Claim 4 (Rejection Sampling).

Algorithm 2, with input the output of the Learning Phase, outputs a state
$x\in[n]$ with probability proportional to $\mathcal{D}(x)$ or outputs
$\perp$. Moreover, Algorithm 3 outputs a perfect sample from $\mathcal{D}$.

###### Proof.

Let us set $A=\sum_{y\in[n]}\mathcal{D}(y)/\widetilde{D}(y)$. At the end of
Line 11, the parameterized CFTP algorithm (Algorithm 2) outputs a sample
$x\propto\mathcal{D}(x)/\widetilde{D}(x)$. This is because the unique
stationary distribution is
$\widetilde{\pi}(x)=(\mathcal{D}(x)/\widetilde{\mathcal{D}}(x))/A$ and the
utility of the CFTP mechanism guarantees that when the Algorithm initiates a
CFTP simulation from time $T=-T^{\star}$, the generated sample
$x=F_{(-T^{\star},0]}=F_{(-\infty,0]}$ will be distributed according to
$\widetilde{\pi}$. This follows from the observation that the distribution of
the state $F_{(-\infty,0]}$ is equal to the distribution of the state
$F_{[0,+\infty)}$ obtained by running the simulation up to the limit
$T\rightarrow+\infty$ forward in time. Hence, we have that $x$ has law the
unique stationary distribution $\widetilde{\pi}$. We now have to remove the
bias induced by the learning step. For the sample $x$, we perform a rejection
sampling process, with acceptance probability $\widetilde{\mathcal{D}}(x)$.
The Algorithm 2 has $n+1$ potential outputs: It either prints $x\in[n]$ or
$\perp$ indicating failure. The sample $x$ is the output of the algorithm with
probability
$\frac{\mathcal{D}(x)/\widetilde{\mathcal{D}}(x)}{A}\cdot\widetilde{\mathcal{D}}(x)=\mathcal{D}(x)/A$.
This holds for any $x\in[n]$. The remaining probability mass is assigned to
$\perp.$ Hence, a single execution of Algorithm 2 outputs a point $x\in[n]$
with probability $\mathcal{D}(x)/A$ and outputs “reject” with probability
$1-1/A$. As we will see, the accepted sample is distributed according to
$\mathcal{D}$. This is because, conditional on acceptance, its mass is exactly
that assigned by $\mathcal{D}$. Specifically, we claim that the output of Line
6 of Algorithm 3 has law $\mathcal{D}$. Observe that the whole stochastic
process of Algorithm 3 outputs $x\in[n]$ with probability

$\displaystyle\sum_{i=0}^{\infty}\Pr[\mathrm{Reject}]^{i}\Pr[x\text{ is
Accepted}]$
$\displaystyle=\sum_{i=0}^{\infty}\Pr[\mathrm{Reject}]^{i}\cdot\frac{(\mathcal{D}(x)/\widetilde{\mathcal{D}}(x))\cdot\widetilde{\mathcal{D}}(x)}{\sum_{y\in[n]}\mathcal{D}(y)/\widetilde{\mathcal{D}}(y)}$
$\displaystyle=\frac{\mathcal{D}(x)}{\sum_{y\in[n]}\mathcal{D}(y)/\widetilde{\mathcal{D}}(y)}\sum_{i=0}^{\infty}\Pr[\mathrm{Reject}]^{i}\,.$

We have that
$\Pr[\mathrm{Reject}]=\frac{1}{A}\sum_{y\in[n]}\frac{\mathcal{D}(y)}{\widetilde{\mathcal{D}}(y)}(1-\widetilde{\mathcal{D}}(y))=1-\frac{1}{A}\in(0,1)\,.$
Hence, we have that the whole stochastic process of Algorithm 3 outputs
$x\in[n]$ with probability

$\displaystyle\sum_{i=0}^{\infty}\Pr[\mathrm{Reject}]^{i}\Pr[x\text{ is
Accepted}]$
$\displaystyle=\frac{\mathcal{D}(x)}{\sum_{y\in[n]}\mathcal{D}(y)/\widetilde{\mathcal{D}}(y)}\cdot
A=\mathcal{D}(x)\,,$

since $\sum_{i\geq 0}\lambda^{i}=\frac{1}{1-\lambda}$ for $\lambda\in(-1,1).$
As a result, we get that the probability that $x$ is the output of the above
stochastic process is $\mathcal{D}(x).$ ∎

These claims complete the proof. Specifically, the total sample complexity is
derived as follows: Since we have to execute the CFTP iterations $\Theta(n)$
times, we should call each CFTP process with $\delta^{\prime}=C\delta/n$ for
some constant $C$. Proposition 12 guarantees that the CFTP process will
terminate using $O(n\log(n)/\lambda(\bm{Q})\cdot\log(n/\delta))$ samples with
probability $1-\Theta(\delta/n)$. Hence, applying the the union bound for the
$\Theta(n)$ CFTP calls and the single call of the learning algorithm, gives
that, with probability $1-\delta$, a number of
$O(n^{2}\log^{2}(n)/\lambda(\bm{Q})\cdot\log(1/\delta))$ samples suffices in
order to get a perfect sample from $\mathcal{D}$. ∎

##  Appendix D The Proof of Lemma 5 (Properties of Rescaled Random Walk)

In the next section, we discuss some useful steps as a warm-up for the proof
of Lemma 5. The proof can be found at the Section D.2.

###  D.1 Sketch of the Idea

By downscaling the transition probabilities (as we will see below), we can
decouple the Markov chain from $\mathcal{D}$. Then, the transition matrix of
the modified Markov chain only depends on the pair distribution $\mathcal{Q}$,
and the convergence of the CFTP algorithm is determined by $\lambda(\bm{Q})$.
Hence, we can sample efficiently even from target distributions $\mathcal{D}$
that may be multimodal or has many low probability points. E.g., in case of a
multimodal stationary distribution $\mathcal{D}$, the spectral properties of
the walk would remind these of a disconnected graph and the sample complexity
of Algorithm 1 would be quite high. We use the estimate
$\widetilde{\mathcal{D}}$ (which we obtained from our learning algorithm) of
the target distribution $\mathcal{D}$ to transform the Markov chain of the
Local Sampling Scheme defined in Equation (2.1) to a modified Markov chain
with an almost uniform stationary distribution. This transformation can be
viewed as a downscaling mechanism:

######  Definition 14 (Downscaling).

Let $p\in(0,1)$ and let $X\sim\mathrm{Be}(p)$ be a Bernoulli random variable.
Let $\lambda\in(0,1).$ The random variable $Y\sim\mathrm{Be}(\lambda p)$ is
called a $\lambda$-downscaler of $X$.

Applying Theorem 4, we can consider that, for any $x\in[n]$, there exists a
coefficient $\widetilde{\mathcal{D}}(x)\approx\mathcal{D}(x)$. The idea is to
use $\widetilde{\mathcal{D}}(x)$ and make the stationary distribution of the
modified chain close to uniform. Intuitively, this transformation should
speedup the convergence of the CFTP algorithm.

######  Step 1.

Implementation of the downscaling and the rescaled matrix
$\widetilde{\bm{M}}$.

We can implement the modified Markov chain via downscaling as follows:
Consider an edge $\\{x,y\\}$ with transition probability pair
$(p_{xy},p_{yx})$. Without loss of generality, we assume that
$\widetilde{\mathcal{D}}(y)>\widetilde{\mathcal{D}}(x)$ (which intuitively
means that we should expect that $p_{xy}>p_{yx}$). Then, the downscaler leaves
$p_{yx}$ unchanged and reduces the mass of $p_{xy}$ to make the two
transitions almost balanced. Our exact sampling algorithm will perform this
downscaling phase to the matrix $\bm{M}$ (see Equation (2.1)).

Consider the transition matrix $\bm{M}$ with
$\bm{M}_{xy}=\mathcal{Q}(x,y)\frac{\mathcal{D}(y)}{\mathcal{D}(x)+\mathcal{D}(y)}$
and $\bm{M}_{xx}=1-\sum_{y\neq x}\bm{M}_{xy}$. Also, let
$\widetilde{\mathcal{D}}$ be an estimate for the distribution $\mathcal{D}$
with some sufficiently small accuracy $\epsilon$, to be chosen (see also
Theorem 4). For the pair $(x,y)$, we modify the transition probability
$p_{xy}:=\bm{M}_{xy}$, only if $p_{xy}>p_{yx}$, to be equal to the following

$\widetilde{p}_{xy}=\mathcal{Q}(x,y)\frac{\mathcal{D}(y)}{\mathcal{D}(x)+\mathcal{D}(y)}\frac{\widetilde{\mathcal{D}}(x)}{\widetilde{\mathcal{D}}(y)}\approx\mathcal{Q}(x,y)\frac{\mathcal{D}(x)}{\mathcal{D}(x)+\mathcal{D}(y)}=p_{yx}\,,$

where we use that $\widetilde{\mathcal{D}}(x)\approx\mathcal{D}(x)$ and
$\widetilde{\mathcal{D}}(y)\approx\mathcal{D}(y)$. The transition probability
from $x$ to $y$ corresponds to a Bernoulli variable $\mathrm{Be}(p_{xy})$,
which is downscaled by
$\widetilde{\mathcal{D}}(x)/\widetilde{\mathcal{D}}(y)<1$. The modified
transition probability $\widetilde{p}_{xy}$ can be implemented by drawing a
$\Lambda\sim\mathrm{Be}(\widetilde{\mathcal{D}}(x)/\widetilde{\mathcal{D}}(y))$
and then drawing a $P\sim\mathrm{Be}(p_{xy})$ (from
$\mathrm{Samp}(\mathcal{Q};\mathcal{D})$) and, finally, realizing the
transition from $x$ to $y$ only if $\Lambda P=1$. This implementation is valid
since the two sources of randomness are independent.

The modified transition matrix $\widetilde{\bm{M}}$ can be written as:

$\widetilde{\bm{M}}=\widetilde{\bm{D}}\otimes\bm{Q}+\bm{Q}\circ[\epsilon_{xy}]\,,$
(D.1)

where $\widetilde{\bm{D}}$ is a symmetric matrix with
$\widetilde{\bm{D}}_{xy}=\mathrm{min}\\{\mathcal{D}(x),\mathcal{D}(y)\\}/(\mathcal{D}(x)+\mathcal{D}(y))=\widetilde{\bm{D}}_{yx}$
and $\bm{Q}$ is the Laplacian matrix of the graph associated with the
$\mathrm{Samp}(\mathcal{Q};\mathcal{D})$ and $\bm{Q}\circ[\epsilon_{xy}]$
denotes the modified Hadamard product777We have that
$\Big{(}\bm{Q}\circ\bm{A}\Big{)}_{xy}=\mathcal{Q}_{xy}A_{xy},\Big{(}\bm{Q}\circ\bm{A}\Big{)}_{xx}=-\sum_{y\neq
x}\mathcal{Q}_{xy}A_{xy}.$ between the Laplacian $\bm{Q}$ and the matrix with
the estimation error $\epsilon_{xy}$ is only non-zero in the modified
transitions $x\to y$, i.e., if $p_{xy}>p_{yx}$ then
$\widetilde{\bm{M}}_{xy}=\mathcal{Q}(x,y)\cdot(\frac{\mathcal{D}(x)}{\mathcal{D}(x)+\mathcal{D}(y)}+\epsilon_{xy})$.
For this pair, we set $\epsilon_{yx}=0$. We denote this error matrix with
$[\epsilon_{xy}]$.

Specifically, for a pair $(x,y)$, before the downscaling, the pairwise
$(x,y)$-comparison corresponds to the random coin
$\left(\frac{\mathcal{D}(y)}{\mathcal{D}(x)+\mathcal{D}(y)},\frac{\mathcal{D}(x)}{\mathcal{D}(x)+\mathcal{D}(y)}\right)$
and let $\mathcal{D}(y)>\mathcal{D}(x)$. After the downscaling, this coin
becomes (locally) almost fair, i.e.,
$(\frac{\mathcal{D}(x)}{\mathcal{D}(x)+\mathcal{D}(y)}+\epsilon_{xy},\frac{\mathcal{D}(x)}{\mathcal{D}(x)+\mathcal{D}(y)}).$

######  Step 2.

Obtaining a simpler matrix $\widetilde{\bm{M}}$ with absolute spectral gap of
the same order.

Since the coins are locally fair, we can work with the following matrix that
has also almost uniform stationary distribution and conductance of the same
order (since making each transition $(x,y)$ more lazy (i.e., increasing the
probability of $x\to x$ and $y\to y$ by some _constant_) and still keeping the
transitions $x\to y$ and $y\to x$ (almost) equal cannot significantly affect
the conductance, i.e., the desired symmetry is preserved). Let
$\widetilde{\bm{M}}=\bm{I}-c\cdot\bm{Q}+\bm{Q}\circ[\epsilon_{xy}]$ for some
_constant_ $0\leq c\leq 1/2$. This matrix with
$c=\min_{(x,y)\in\mathcal{E}}\mathcal{D}(x)/(\mathcal{D}(x)+\mathcal{D}(y))$
can be obtained by further performing downscaling at each transition (and $c$
is constant due to 2); then each transition will be equal to the minimum
global transition probability (potentially with some $O(1/\sqrt{n})$ noise
term). For simplicity, we let $c=1/2$. We get the following matrix:

$\widetilde{\bm{M}}_{xy}=\mathcal{Q}(x,y)\Big{(}\frac{1}{2}+\epsilon_{xy}\Big{)},~{}\widetilde{\bm{M}}_{yx}=\mathcal{Q}(x,y)\Big{(}\frac{1}{2}+\epsilon_{yx}\Big{)}\,,$
(D.2)

and

$\widetilde{\bm{M}}_{xx}=1-\sum_{y\neq
x}\widetilde{\bm{M}}_{xy}=1-\frac{1}{2}\sum_{y\neq
x}\mathcal{Q}(x,y)-\sum_{y\neq x}\epsilon_{xy}\mathcal{Q}(x,y)\,.$

Also, we observe that the generator of the Markov chain
$\bm{I}-\widetilde{\bm{M}}$ is close up to scaling to the Laplacian $\bm{Q}$.
This explains why the convergence time of the Markov chain with transition
matrix $\widetilde{\bm{M}}$ spectrally depends only on the distribution
$\mathcal{Q}$. Specifically, the Laplacian $\bm{Q}$ corresponds to the
dominant component for the convergence of the algorithm and the Hadamard
product is the low-order noise induced by the chain transformation. As we will
see, the larger the smallest non-zero eigenvalue of the Laplacian matrix and,
hence, the larger the spectral gap of the transition matrix of the transformed
chain, the faster the Markov chain converges to its stationary distribution.

###  D.2 The Proof of Lemma 5

The following lemma summarizes the key properties of the downscaled random
walk.

######  Lemma 15 (Properties of Rescaled Random Walk).

Let $\mathcal{D}$ be a distribution on $[n]$ and consider an
$\epsilon$-relative approximation $\widetilde{\mathcal{D}}$ of $\mathcal{D}$,
as in Theorem 4 with $\epsilon=O(1/\sqrt{n})$. Consider the transition matrix
$\bm{M}$ of the Local Sampling Scheme $\mathrm{Samp}(\mathcal{Q};\mathcal{D})$
(see Equation (2.1)), and let $\widetilde{\bm{M}}$ be the
$\widetilde{\mathcal{D}}$-scaling of $\bm{M}$ (see Equation (D.1)). Then, the
following hold.

  * $(i)$

The transition matrix $\widetilde{\bm{M}}$ has stationary distribution
$\widetilde{\pi}_{0}(x)=\Theta(1/n)$ for all $x\in[n]$.

  * $(ii)$

The absolute spectral gap $\Gamma(\widetilde{\bm{M}})$ of $\widetilde{\bm{M}}$
and the minimum non-zero eigenvalue $\lambda(\bm{Q})$ of the Laplacian matrix
$\bm{Q}$ satisfy $\Gamma(\widetilde{\bm{M}})=\Omega(\lambda(\bm{Q}))$.

  * $(iii)$

For any $\epsilon_{0}\in(0,1)$, the mixing time of the transition matrix
$\widetilde{\bm{M}}$ is
$T_{\mathrm{mix}}(\widetilde{\bm{M}};\epsilon_{0})=O\left(\frac{\log(n/\epsilon_{0})}{\lambda(\bm{Q})}\right)$.

In the above, we can choose $\epsilon_{0}=1/4$. The proof of the Lemma 15
follows.

###### Proof _of_ Lemma 15.

We remark that it suffices to work with the matrix of Equation (D.2), since
the results will only change by at most some constant. Observe that the
spectral gap of the matrix of Equation (D.1) is at least as high as one of the
matrix $\bm{I}-c\cdot\bm{Q}+\bm{Q}\circ[\epsilon_{xy}]$ with constant
$c=\min_{(x,y)\in\mathcal{E}}\mathcal{D}(x)/(\mathcal{D}(x)+\mathcal{D}(y))=\Theta(\phi)$
and this matrix has spectral gap of the same order as the matrix of Equation
(D.2). We break the proof into three claims.

######  Claim 5.

The transition matrix $\widetilde{\bm{M}}$ has stationary distribution
$\widetilde{\pi}_{0}$, that satisfies $\widetilde{\pi}_{0}(x)=\Theta(1/n)$ for
any $x\in[n]$.

###### Proof.

The chain $\widetilde{\bm{M}}$ remains irreducible and the detailed balance
equations of the matrix $\widetilde{\bm{M}}$ satisfy:

$\frac{\widetilde{p}_{xy}}{\widetilde{p}_{yx}}=\frac{\mathcal{Q}(x,y)\mathcal{D}(y)/\widetilde{\mathcal{D}}(y)}{\mathcal{Q}(y,x)\mathcal{D}(x)/\widetilde{\mathcal{D}}(x)}\,.$

Since $\widetilde{\mathcal{D}}$ is an $\epsilon$-relative approximation of
$\mathcal{D}$ with $\epsilon=O(1/\sqrt{n})$, it holds that for any $x\in[n]$:

$\mathcal{D}(x)/\widetilde{\mathcal{D}}(x)\in[1-\epsilon,1+\epsilon]\,,$

and hence the unique stationary distribution $\widetilde{\pi}_{0}$ satisfies:
$\widetilde{\pi}_{0}(x)/\widetilde{\pi}_{0}(y)\in[1-\epsilon,1+\epsilon]$ for
any $x,y\in[n]$. So, this holds for
$x=\operatorname*{argmax}\widetilde{\pi}_{0}$ and
$y=\operatorname*{argmin}\widetilde{\pi}_{0}$ and, since it should hold that
$\sum_{x}\widetilde{\pi}_{0}(x)=1$, it must be the case that
$\widetilde{\pi}(x)=\Theta(1/n)$. ∎

######  Claim 6.

For the absolute spectral gap $\Gamma(\widetilde{\bm{M}})$ of
$\widetilde{\bm{M}}$ and the minimum non-zero eigenvalue $\lambda(\bm{Q})$ of
the Laplacian matrix $\bm{Q}$, it holds that
$\Gamma(\widetilde{\bm{M}})=\Omega(\lambda(\bm{Q}))$.

###### Proof.

For $x\neq y$, there exist error estimates $\epsilon_{xy}$ and
$\epsilon_{yx}$; one of them is zero and the other’s absolute value is of
order $O(1/\sqrt{n}).$ Our goal is to control the absolute spectral gap
$\Gamma(\widetilde{\bm{M}})$. Recall that the transition matrix can be written
as follows

$\widetilde{\bm{M}}=\bm{I}-\frac{1}{2}\bm{Q}+\bm{Q}\circ[\epsilon_{xy}]\,,$

where $\epsilon_{xy}\in[-\epsilon,\epsilon]$ for some
$\epsilon=O(1/\sqrt{n})$. Note that $[\epsilon_{xy}]$ denotes the $n\times n$
error matrix and the operator $\circ$ denotes the standard Hadamard product.
Note that the matrix $\bm{I}-\frac{1}{2}\bm{Q}$ is symmetric, while
$\bm{Q}\circ[\epsilon_{xy}]$ needs not to be.

Let $\lambda_{1}>\lambda_{2}\geq\dots\geq\lambda_{n}$ be the spectrum of the
transition matrix $\widetilde{\bm{M}}$. Also, the matrix is right stochastic
with $\widetilde{\bm{M}}\bm{1}=\bm{1}$ and $\lambda_{1}=1$. Let
$\Gamma(\widetilde{\bm{M}})$ be the absolute spectral gap, that is strictly
positive by aperiodicity and irreducibility. The Weyl’s Inequality for general
matrices ([HJ12, Tao12]) describes a multiplicative majorization between the
ordered absolute eigenvalues and singular values of $\bm{A}$ and gives a
useful perturbation bound. Recall that
$\sigma_{i}^{2}(\bm{A})=\lambda_{i}(\bm{A}\bm{A}^{*})=\lambda_{i}(\bm{A}^{*}\bm{A})$
for an arbitrary matrix $\bm{A}\in\mathbb{C}^{m\times n}$.

######  Fact 2 (Weyl’s Inequality (see [HJ12])).

Consider the matrices $\bm{A},\bm{E}\in\mathbb{C}^{n\times n}$. Define the
singular value of $\bm{A}$ in decreasing order (counting multiplicity)
$\sigma_{1}(\bm{A})\geq\ldots\geq\sigma_{n}(\bm{A})\geq 0$. Then, for
$k=1,\dots,n$, the following hold

  * •

$\prod_{i\in[k]}|\lambda_{i}(\bm{A})|\leq\prod_{i\in[k]}\sigma_{i}(\bm{A})$
with
$|\lambda_{1}(\bm{A})|\geq|\lambda_{2}(\bm{A})|\geq\dots|\lambda_{n}(\bm{A})|$.

  * •

$|\sigma_{k}(\bm{A}+\bm{E})-\sigma_{k}(\bm{A})|\leq\|\bm{E}\|_{2}$ .

For the matrix $\widetilde{\bm{M}}$, we have that $|\lambda_{1}|=1$ and the
second largest absolute eigenvalue is $\max\\{|\lambda_{2}|,|\lambda_{n}|\\}$.
We can apply the first property for the singular values from the Fact 2 and
get that

$\Gamma(\widetilde{\bm{M}})=\lambda_{1}-\max\\{|\lambda_{2}|,|\lambda_{n}|\\}\geq
1-\sigma_{1}(\widetilde{\bm{M}})\sigma_{2}(\widetilde{\bm{M}})\,,$

where $\sigma_{1},\sigma_{2}$ correspond to the two largest singular values.
Observe that the largest singular value is
$\sigma_{1}(\widetilde{\bm{M}})=\max_{\|\bm{v}\|_{2}=1}\|\widetilde{\bm{M}}\bm{v}\|=1$
and, using the second property of Fact 2, we can control the perturbation of
the second largest singular value

$\sigma_{2}\left(\bm{I}-\frac{1}{2}\bm{Q}\right)-\|\bm{Q}\circ[\epsilon_{xy}]\|_{\mathrm{2}}\leq\sigma_{2}(\widetilde{\bm{M}})\leq\sigma_{2}\left(\bm{I}-\frac{1}{2}\bm{Q}\right)+\|\bm{Q}\circ[\epsilon_{xy}]\|_{\mathrm{2}}\,.$

In order to lower bound the absolute spectral gap, it suffices to upper bound
the second largest singular value. But, for the second largest singular
eigenvalue of the real symmetric matrix $\bm{I}-\bm{Q}/2$, it holds that

$\sigma_{2}\left(\bm{I}-\frac{1}{2}\bm{Q}\right)=\max\left\\{\left|\lambda_{2}\left(\bm{I}-\frac{1}{2}\bm{Q}\right)\right|,\left|\lambda_{n}\left(\bm{I}-\frac{1}{2}\bm{Q}\right)\right|\right\\}\,.$

Since the matrix $\bm{I}-\bm{Q}/2$ is symmetric, we have that

$\sigma_{2}\left(\bm{I}-\frac{1}{2}\bm{Q}\right)=\max\left\\{\Big{|}\max_{\bm{v}\perp\bm{1},\|\bm{v}\|_{2}=1}\bm{v}^{\top}\left(\bm{I}-\frac{1}{2}\bm{Q}\right)\bm{v}\Big{|},\Big{|}\min_{\|\bm{v}\|_{2}=1}\bm{v}^{\top}\left(\bm{I}-\frac{1}{2}\bm{Q}\right)\bm{v}\Big{|}\right\\}\,,$

and, hence, if we let $0=\mu_{1}\leq\mu_{2}\leq\ldots\leq\mu_{n}$ be the
spectrum of the Laplacian matrix $\bm{Q}$, we have that

$\sigma_{2}\left(\bm{I}-\frac{1}{2}\bm{Q}\right)=\max\Big{\\{}\Big{|}1-\frac{1}{2}\mu_{2}(\bm{Q})\Big{|},\Big{|}1-\frac{1}{2}\mu_{n}(\bm{Q})\Big{|}\Big{\\}}\,.$

For the Laplacian matrix $\bm{Q}$, observe that $|\bm{Q}_{xx}|=|\sum_{y\sim
x}\mathcal{Q}(x,y)|\leq 1$. Hence, from Gershgorin’s Theorem (see Lemma 11),
we get that all the eigenvalues of the Laplacian lie on the real axis (since
the matrix is symmetric) and all the disks $B(x,r)$ are centered at points
$|x|\leq 1$ and the radii $r$ are upper bounded by $1$. Hence, it holds that
$0\leq\mu_{2}(\bm{Q})\leq\mu_{n}(\bm{Q})\leq 2$. So, we can take that
$\sigma_{2}(\bm{I}-\frac{1}{2}\bm{Q})=1-\frac{1}{2}\mu_{2}(\bm{Q})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}1-\frac{1}{2}\lambda(\bm{Q})\geq
0$, where $\lambda(\bm{Q})$ denotes the smallest non-zero eigenvalue of the
Laplacian matrix $\bm{Q}$.

Now, we have that
$\|\bm{Q}\circ[\epsilon_{xy}]\|_{\mathrm{2}}=\max_{\|\bm{v}\|_{2}=1}\|(\bm{Q}\circ[\epsilon_{xy}])\bm{v}\|_{2}=O(1/\sqrt{n})\cdot\mu_{n}(\bm{Q})=O(1/\sqrt{n})$
and so we get
$\sigma_{2}(\widetilde{\bm{M}})\leq\left(1-\Theta(\lambda(\bm{Q}))\right)+O(1/\sqrt{n})$.
This implies that

$\sigma_{2}(\widetilde{\bm{M}})\lesssim\max\\{1-\Theta(\lambda(\bm{Q})),O(1/\sqrt{n})\\}\,,$

and so

$\Gamma(\widetilde{\bm{M}})\gtrsim
1-\max\\{1-\Theta(\lambda(\bm{Q})),O(1/\sqrt{n})\\}$

This implies that $\Gamma(\widetilde{\bm{M}})=\Omega(\lambda(\bm{Q}))$ for
sufficiently large $n$. The above proof will be similar for any matrix
$\widetilde{\bm{M}}=I-c\cdot\bm{Q}+\bm{Q}\cdot[\epsilon_{xy}]$ for $0\leq
c\leq 1/2$ which is obtained by the discussion of Step 2. Such a matrix has an
absolute spectral gap that is, on the one side, lower bounded by
$\lambda(\bm{Q})$ (as we showed above) and, on the other side, upper bounded
by the absolute spectral gap of the original matrix of Equation (D.1). This
completes the proof. ∎

We continue with the mixing time result of the transformed Markov chain. Since
we have shown that the absolute spectral gap is $\Omega(\lambda(\bm{Q}))$, we
can directly get the desired result, whose proof relies on the analysis of
Lemma 9 and can be found, e.g., in [LP17].

######  Claim 7.

The $\epsilon_{0}$-mixing time of the transition matrix $\widetilde{\bm{M}}$
is equal to

$T_{\mathrm{mix}}(\widetilde{\bm{M}};\epsilon_{0})=O(\log(n/\epsilon_{0})/\lambda(\bm{Q}))\,.$

In the above claim, we choose $\epsilon_{0}=1/4$.

Combining these claims, the proof is completed. ∎

##  Appendix E  The Proof of Theorem 4 (Main Result of Learning Phase)

The first step of Algorithm 3 is to learn the target distribution
$\mathcal{D}$ in $\epsilon$-relative error for some $\epsilon>0$ and pass it
as input to Algorithm 2. In this section, we will provide the learning results
for abstract $\epsilon$; however, our algorithm applies these results with
$\epsilon=1/\sqrt{n}$. For two distributions
$\mathcal{D},\widetilde{\mathcal{D}}$ with ground set $[n]$, we introduce the
sequence/list (of length $n$)
$1-\mathcal{D}/\widetilde{\mathcal{D}}:=(1-\mathcal{D}(x)/\widetilde{\mathcal{D}}(x))_{x\in[n]}.$
Observe that the pair of sequences
$(1-\mathcal{D}/\widetilde{\mathcal{D}},1-\widetilde{\mathcal{D}}/\mathcal{D})$
captures the relative error between the two distributions. The sample
complexity of the task of learning $\mathcal{D}$ in $\epsilon$-relative error
is summarized by the following (restatement of Theorem 4):

######  Theorem.

For any $\epsilon,\delta>0$, there exists an algorithm (Algorithm 5) that
draws
$N=O\left(\frac{n}{\lambda(\bm{Q})\epsilon^{2}}\log(\frac{1}{\delta})\right)$
samples from a Local Sampling Scheme $\mathrm{Samp}(\mathcal{Q};\mathcal{D})$
satisfying Assumptions 1 and 2, runs in time polynomial in $N$, and, with
probability at least $1-\delta$, computes an estimate
$\widetilde{\mathcal{D}}$ of the target distribution $\mathcal{D}$, that
satisfies the following relative error guarantee
$\max\big{\\{}\|1-\mathcal{D}/\widetilde{\mathcal{D}}\|_{\infty},\|1-\widetilde{\mathcal{D}}/\mathcal{D}\|_{\infty}\big{\\}}\leq\epsilon$.

Algorithm 5 Learn using Shifting from (Pairwise) Local Sampling Schemes

1:procedure Learn-Shift($\epsilon,\delta$) $\triangleright$ _Sample access to
oracle $\mathrm{Samp}(\mathcal{Q};\mathcal{D})$_.

2: Set
$N=\Theta\left(\frac{n}{\lambda(\bm{Q})\epsilon^{2}}\log(\frac{1}{\delta})\right)$.

3: Draw $N$ samples of the form
$((x_{i},y_{i}),q_{i})\in\mathcal{E}\times\\{0,1\\}.$

4: Obtain $\widehat{\bm{z}}$ using Algorithm 6 (the estimation vector
satisfies $\langle\bm{1},\widehat{\bm{z}}\rangle=0$ and must be shifted).

5: Compute $C=\log\left(\sum_{x\in[n]}\exp(\widehat{z}_{x})\right)$

6: Set $z_{x}^{\prime}=\widehat{z}_{x}-C$ for any $x\in[n]$ $\triangleright$
_See Section F.2_.

7: Output $\bm{z}^{\prime}$ $\triangleright$ _The output satisfies
$\|\bm{z}-\bm{z}^{\prime}\|_{\infty}\leq\epsilon$_.

Let us sketch the proof of Theorem 4. For an arbitrary weight vector
$\bm{w}\in\mathbb{R}^{n}_{>0}$, we use the re-parameterization
$z_{i}=\log(w_{i})$. When in addition $\bm{w}\in\Delta^{n}$, i.e., $\bm{w}$ is
a probability distribution over $[n]$ and is usually denoted by $\mathcal{D}$,
we call $\bm{z}$ the natural parameter vector of $\mathcal{D}$. Recall that
$\bm{Q}$ is a Laplacian matrix where $\bm{Q}_{xy}=-\mathcal{Q}(x,y)$, i.e., it
is the Laplacian matrix of the graph $G_{\mathcal{Q}}$ weighted by the mass
assigned by $\mathcal{Q}$. The proof goes as follows: we draw i.i.d. samples
from the LSS. We first apply the following result which is a modification of
the results of [SBB+16].

######  Theorem 16 (Variant of [SBB+16]).

Let $G_{\mathcal{Q}}$ be the graph of the support $\mathcal{E}$ of
$\mathcal{Q}$ satisfying 1 with $|V|=n$ vertices and associated weighted
Laplacian matrix $\bm{Q}$ (see Equation (1.1)). Let $\bm{L}$ be the associated
empirical Laplacian matrix with $\mathbb{E}[\bm{L}]=\bm{Q}$. Consider a vector
$\bm{z}\in\mathbb{R}^{n}$ with $\langle\bm{1},\bm{z}\rangle=0$ satisfying the
constraint $\max_{x,y\in\mathcal{E}}|z_{x}-z_{y}|\leq\log(\phi)$ for some
constant $\phi>1$. Let $\mathrm{Ex}(\mathcal{Q};\bm{z})$ be the oracle that
generates the example $\\{(x,y),q\\}$ as follows: the edge $(x,y)$ is chosen
with probability $\mathcal{Q}(x,y)$ and the bit $q$ is set to $1$ with
probability $\exp(z_{x})/(\exp(z_{x})+\exp(z_{y}))$; otherwise it is set to
$0$. For any $\epsilon,\delta>0$, there is a maximum likelihood estimator of
$\bm{z}$ (Algorithm 6) which draws
$N=O\left(\frac{n}{\lambda(\bm{L})\epsilon^{2}}\log(\frac{1}{\delta})\right)$
samples from $\mathrm{Ex}(\mathcal{Q};\bm{z})$ and computes, in time
polynomial in the number of samples $N$, an estimate $\widehat{\bm{z}}$ such
that $\|\bm{z}-\widehat{\bm{z}}\|_{2}\leq\epsilon$, with probability at least
$1-\delta$.

The proof can be found at Section F.1. Some comments are in order:

  1. 1.

In the work of [SBB+16], the target is the (re-parameterized) weights vector
$\bm{z}^{\star}\in\mathbb{R}^{n}$ satisfying the conditions
$\langle\bm{1},\bm{z}^{\star}\rangle=0$ and $\|\bm{z}^{\star}\|_{\infty}\leq
B$ for some constant $B$. The provided algorithm minimizes the empirical log-
likelihood and, using $O(n/(\lambda(\bm{Q})\cdot\epsilon^{2})))$ samples,
computes an estimate $\widehat{\bm{z}}$ so that
$\|\bm{z}^{\star}-\widehat{\bm{z}}\|_{2}\leq\epsilon$.

  2. 2.

Our provided Algorithm 6 has the same guarantees as the algorithm of [SBB+16]
but the target weight vector satisfies the conditions
$\langle\bm{1},\bm{z}^{\star}\rangle=0$ and $|z^{\star}_{x}-z^{\star}_{y}|\leq
B$ for some constant $B$. The algorithm draws sufficiently many samples and
then minimizes the empirical negative log-likelihood objective over an
appropriately selected constrained set (see Section F.1). Observe that for a
Local Sampling Scheme with distribution $\mathcal{D}$, the natural parameter
vector $\bm{z}^{\star}$ does not satisfy
$\langle\bm{1},\bm{z}^{\star}\rangle=0$ and so Theorem 16 and Algorithm 6
cannot be directly applied.

  3. 3.

We will discuss (Section F.2) how to apply Algorithm 6 to distributions, i.e.,
target weight vectors with $\langle\bm{1},\bm{z}^{\star}\rangle\leq 0.$ This
results in learning the target vector $\bm{z}^{\star}$ in $L_{\infty}$ norm
(which is sufficient for Algorithm 3).

The algorithm of Theorem 16 follows.

Algorithm 6 Learning from (Pairwise) Local Sampling Schemes

1:procedure Learn($\epsilon,\delta$) $\triangleright$ _Sample access to oracle
$\mathrm{Ex}(\mathcal{Q};\bm{z})$_.

2: Set
$N=\Theta\left(\frac{n}{\lambda(\bm{Q})\epsilon^{2}}\log(\frac{1}{\delta})\right)$.

3: Draw $N$ samples of the form
$((x_{i},y_{i}),q_{i})\in\mathcal{E}\times\\{0,1\\}.$

4: Compute the empirical negative log-likelihood objective
$L_{N}(\bm{z};\\{((x_{i},y_{i}),q_{i})\\}_{i\in[N]})=-\frac{1}{N}\sum_{i=1}^{N}\bm{1}\\{q_{i}=1\\}z_{x_{i}}+\bm{1}\\{q_{i}=0\\}z_{y_{i}}-\log(\exp(z_{x_{i}})+\exp(z_{y_{i}}))\,.$

5: Minimize $L_{N}$ using gradient descent in the subspace $\Omega_{\phi}$.
$\triangleright$ _See Section F.1_.

6: Output the guess vector $\widehat{\bm{z}}$.

To complete the proof of Theorem 4, we combine the upcoming Lemma 17 and Lemma
18 with the main result of [SBB+16].

Lemma 17 states that, after drawing sufficiently many samples from
$\mathrm{Samp}(\mathcal{Q};\mathcal{D})$, we can construct an empirical
Laplacian matrix $\bm{L}$, whose Fiedler eigenvalue is of the same order as
the one of the unknown matrix $\bm{Q}$ of the Local Sampling Scheme. Moreover,
notice that Lemma 17 implies that, with high probability, the graph induced by
the matrix $\bm{L}$ is connected, i.e., $\lambda(\bm{L})>0$.

######  Lemma 17 (Concentration of Empirical Fiedler Eigenvalue).

Let $\epsilon>0$ and $V=\\{\bm{v}\in\\{\pm
1,0\\}^{n}:\bm{v}=\bm{e}_{i}-\bm{e}_{j}\text{ for any }i<j\\}$ and let
$\mathcal{Q}$ be a distribution over $V$. Let $\bm{Q}$ be the Laplacian matrix
of $\mathcal{Q}$ (see Equation (1.1)). There exists an algorithm that uses
$O(\log(n/\delta)/(\lambda(\bm{Q})\epsilon^{2}))$ samples from $\mathcal{Q}$
and computes a matrix $\bm{L}$ that satisfies
$|\lambda(\bm{L})-\lambda(\bm{Q})|\leq\epsilon\lambda(\bm{Q})$, with
probability at least $1-\delta$, where $\lambda(\cdot)$ denotes the second
smallest eigenvalue of a Laplacian matrix.

###### Proof.

Let $\\{\bm{X}_{i}=\bm{v}_{i}\bm{v}_{i}^{T}\\}$ be a finite sequence of $m$
independent symmetric square matrices (of common dimension $n$) with
$\bm{v}_{i}\sim\mathcal{Q}$, let $\bm{L}=\frac{1}{m}\sum_{i=1}^{m}\bm{X}_{i}$,
and let $M=\lambda_{\max}(\bm{X}_{i})=\Theta(1)$. To deal with the second
smallest eigenvalue of the Laplacian matrix, we can project the Laplacian
matrix $\bm{L}$ to the orthogonal complement of the vector $\bm{1}$. Under
this transformation, we can obtain a matrix $\bm{L}^{\prime}$ that has two
crucial properties. First, it still is a sum of independent positive
semidefinite terms and its minimum eigenvalue coincides with the Fiedler
eigenvalue of $\bm{L}$. To achieve this, as in [Tro15], we introduce the
transformation $\bm{R}\in\mathbb{R}^{(n-1)\times n},$ that satisfies
$\bm{R}\bm{R}^{T}=\bm{I}_{n-1}$ and $\bm{R}\bm{1}=\bm{0}$. Hence, it holds
that $\bm{L}^{\prime}=\bm{R}\bm{L}\bm{R}^{T}$ and
$\mathbb{E}[\bm{L}^{\prime}]=\bm{R}\bm{Q}\bm{R}^{T}$. Then,

$\Pr_{\bm{v}_{1..m}\sim\mathcal{Q}^{\otimes
m}}\Big{[}\lambda(\bm{L})\leq(1-\epsilon)\lambda(\bm{Q})\Big{]}=\Pr_{\bm{v}_{1..m}\sim\mathcal{Q}^{\otimes
m}}\Big{[}\lambda_{\min}(\bm{L}^{\prime})\leq(1-\epsilon)\lambda_{\min}(\bm{R}\bm{Q}\bm{R}^{T})\Big{]}\,.$

Using Proposition 7, we get that:

$\Pr_{\bm{v}_{1..m}\sim\mathcal{Q}^{\otimes
m}}\Big{[}\lambda(\bm{L})\leq(1-\epsilon)\lambda(\bm{Q})\Big{]}\leq
n\left(e^{-\epsilon}/(1-\epsilon)^{1-\epsilon}\right)^{\lambda_{\min}(\bm{R}\bm{Q}\bm{R}^{T})/M}\leq
n\exp\left(-\epsilon^{2}m\frac{\lambda(\bm{Q})}{2M}\right)\,.$

To upper bound this probability by $\delta$, it suffices to draw
$O(\log(n/\delta)/(\epsilon^{2}\lambda(\bm{Q})))$ samples. So, we get that,
with probability at least $1-\delta$, the Fiedler eigenvalue of the empirical
matrix $\bm{L}$ satisfies $\lambda(\bm{L})\in(1\pm\epsilon)\lambda(\bm{Q})$.
Moreover, the second smallest eigenvalue of the empirical Laplacian matrix
$\bm{L}$ is strictly positive and thus, the induced graph is connected (recall
that the number of connected components in the graph is the dimension of the
nullspace of the Laplacian matrix). ∎

The next lemma guarantees that learning the natural parameters of any
distribution in $L_{\infty}$ norm is sufficient for learning the distribution
with small relative error over its support.

######  Lemma 18 (Stability of Relative Error).

Let $\mathcal{D}_{1},\mathcal{D}_{2}$ be discrete distributions supported on
the ground set $[n]$ and let $\bm{z}_{1},\bm{z}_{2}$ be the corresponding
natural parameter vectors. For any sufficiently small accuracy parameter
$\epsilon>0$, if it holds that
$\|\bm{z}_{1}-\bm{z}_{2}\|_{\infty}\leq\epsilon~{}$, then the distributions