{'fraction_non_alphanumeric': 0.05481227775993911, 'fraction_numerical': 0.014960695470167802, 'mean_word_length': 3.8432208270936528, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 23, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 28, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Metropolized Knockoff Sampling

Model-X knockoffs is a wrapper that transforms essentially any feature importance measure into a variable selection algorithm, which discovers true effects while rigorously controlling the expected fraction of false positives. A frequently discussed challenge to apply this method is to construct knockoff variables, which are synthetic variables obeying a crucial exchangeability property with the explanatory variables under study. This paper introduces techniques for knockoff generation in great generality: we provide a sequential characterization of all possible knockoff distributions, which leads to a Metropolis-Hastingsformulation of an exact knockoff sampler. We further show how to use conditional independence structure to speed up computations. Combining these two threads, we introduce an explicit set of sequential algorithms and empirically demonstrate their effectiveness. Our theoretical analysis proves that our algorithms achieve near-optimal computational complexity in certain cases. The techniques we develop are sufficiently rich to enable knockoff sampling in challenging models including cases where the covariates are continuous and heavy-tailed, and follow a graphical model such as the Ising model.


Introduction
In modern science, researchers often have access to large data sets featuring comprehensive measurements about some phenomenon of interest. The question is then to discover meaningful relationships between an outcome and all the measured covariates. While it is often expected that only a small fraction of the covariates may be associated with the outcome, the relevance of any particular variable is unknown a priori. For instance, a researcher may be interested in understanding which of the thousands of gene-expression profiles may help determine the severity of a tumor. In such circumstances, the researcher often relies on statistical algorithms to sift through large data sets and find those promising candidates, making variable selection a topic of central importance in contemporary statistical research.
The knockoff filter (Barber and Candès 2015;Candès et al. 2018) has recently emerged as a useful framework for performing controlled variable selection, allowing the user to convert any black-box feature importance measure into a variable selection procedure while rigorously controlling the expected fraction of false positives. This means that the statistician can use essentially any blackbox importance measure to return a list of variables with the guarantee that, on the average, the ratio between the number of false positives-loosely speaking, a false positive is a variable that does not influence the response, see Candès et al. (2018)-and the total number of reported variables is below a user-specified threshold. The strength of this method is that the guarantees hold in finite samples and in situations where nothing can be assumed about the dependence between the response and the explanatory variables. Instead, the statistician must have knowledge of the distribution of the explanatory variables. When this happens to be the case, a remaining challenge is the ability to generate the knockoffs, a set of synthetic variables, which can essentially be used as negative controls; these fake variables must mimic the original variables in a particular way without having any additional predictive power. In sum, constructing valid knockoff distributions and sampling mechanisms across a wide range of covariate models is critical to deploying the knockoff filter in a number of applications.

Our contribution
This paper describes a theory for sampling knockoff variables and introduces a general and efficient sampler inspired by ideas from Markov chain Monte Carlo (MCMC). Before moving on, we pause to explicitly mention the two main considerations one should keep in mind when constructing knockoffs: Computation. How can we efficiently sample nontrivial knockoffs?
Statistical power. How can we generate knockoffs that will ultimately lead to powerful variable selection procedures? On this note, it has been observed that knockoffs that are less correlated with the original variables lead to higher power (Barber and Candès 2015;Candès et al. 2018) and, therefore, low correlation must be a design objective.
Having said that, our work makes several specific contributions.
1. Characterization of all knockoff distributions. We provide a sequential characterization of every valid knockoff distribution. Furthermore, we introduce a connection linking pairwise exchangeability between original and knockoff variables to reversible Markov chains, enabling the use of powerful sampling tools from computational statistics.

Related literature
This paper draws most heavily on Candès et al. (2018), which builds on Barber and Candès (2015) to introduce the model-X knockoff framework. In particular, the former reference proposes the Sequential Conditional Independent Pairs (SCIP) procedure for knockoff generation; this is the only known generic knockoff sampler to date, which shall serve as our starting point. The SCIP procedure, however, is only abstractly specified and prior to this paper, implementations were only available for Gaussian distributions and discrete Markov chains. Briefly, Sesia et al. (2018) developed a concrete SCIP algorithm for discrete Markov chains, and then leveraged this construction to sample knockoffs for covariates following hidden Markov models widely used in genome-wide association studies. Similarly relevant is the work of Gimenez et al. (2018), which developed a sampling strategy for a restricted class of Bayesian networks, most notably Gaussian mixture models. In contrast, we address here knockoff sampling for a much larger class of distributions, namely, arbitrary graphical models. We also describe the form of all valid knockoff sampling strategies, thereby providing a framework possibly enabling the construction of future knockoff sampling algorithms. Hence, our work may be of value to the increasing number of researchers deploying the knockoff framework for feature selection in a variety of applications including neural networks (Lu et al. 2018), time-series modeling ), Gaussian graphical model structure learning (Zheng et al. 2018), and biology (Xiao et al. 2017;Gao et al. 2018). Lastly, we close by emphasizing that our contribution is very different from a new strand of research introducing approximate knockoffs generated with techniques from deep learning (Romano et al. 2018;Jordon et al. 2019;Liu and Zheng 2018). While these approaches are tantalizing and demonstrate promising empirical performance in lowdimensional situations, they currently lack formal guarantees about their validity.
2 Characterizing knockoff distributions 2.1 Knockoff variables mixed second moments of (X,X) must match. Assume the second moments of X exist and write Σ = Cov(X). Then the covariance of the vector (X,X) must take the form Cov(X,X) = Γ(s) : where s ∈ R p is any vector such that the right-hand side is positive semi-definite. In other words, for each pair (i, j) with i = j, we have Cov(X i ,X j ) = Cov(X i , X j ).
We are interested in constructing knockoff variables and below we call a knockoff sampler a procedure that takes as inputs a distribution P and a sample X ∼ P and returnsX such that (1) holds. Nontrivial samplers have been demonstrated in a few cases, for instance, when X ∼ N (0, Σ) is multivariate Gaussian. In this case, Candès et al. (2018) show that if (X,X) is jointly Gaussian with mean zero and covariance Γ(s), then the entries ofX are knockoffs for X. One can say that appropriately matching the first two moments is sufficient to generate knockoffs in the special case of the multivariate normal distribution. However, this does not extend and matching the first two moments is in general not sufficient; to be sure, (1) requires that all moments match appropriately.
Gibbs measures. As a motivating example, consider the Ising model, a frequently discussed family of Gibbs measures first introduced in the statistical physics literature (Ising 1925). In this model, the random vector X ∈ {−1, 1} d 1 ×d 2 defined over a d 1 × d 2 X ∈ {−1, 1} d 1 ×d 2 grid has a probability mass function (PMF) of the form here, I = {(i 1 , i 2 ) : 1 ≤ i 1 ≤ d 1 , 1 ≤ i 2 ≤ d 2 } is the grid and α and β are parameters. As we have seen, knockoffsX for X must marginally follow the Ising distribution (3). Furthermore,X must be dependent on X in such a way that any vector of the form {(Z 1 , . . . , Z p ) : Z j = X j or Z j =X j , 1 ≤ j ≤ p} has PMF given by (3). It is tempting to naïvely define a joint PMF for (X,X) as Although the joint distribution is symmetric in X s andX s , the marginal distribution of X is not an Ising model! Hence, this is not a valid joint distribution. Other than the trivial constructioñ X = X, it is a priori unclear how one would construct knockoffs. Any distribution continuous or discrete factoring over a grid poses a similar challenge.

SCIP and its limitations
The only generic knockoff sampler one can find in the literature is SCIP from Candès et al. (2018), given in Procedure 1. While this procedure provably generates valid knockoffs for any input distribution, there are two substantial limitations. The first is that SCIP is only given abstractly; it is challenging to specify L(X j | X -j ,X 1:(j−1) ), b let alone to sample from it. As a result, it is b We use L(W1 | W2) to denote the conditional distribution of W1 given W2. We use the subscript 1 : 0 to mean an empty vector.
only known how to implement SCIP for very special models such as discrete Markov chains and Gaussian distributions. The second limitation is that SCIP is not able to generate all valid knockoff distributions. Recall that we want knockoffs to have low correlations with the original variables so that a feature importance statistic will correctly detect true effects. To achieve this goal, we might need a wider range of sampling mechanisms.
Procedure 1: Sequential Conditional Independent Pairs (SCIP) for j = 1 to p do SampleX j from L(X j | X -j ,X 1:(j−1) ), conditionally independently from X j end

Sequential formulation of knockoff distributions
We begin by introducing a sequential characterization of all valid knockoff distributions, which will later lead to a new class of knockoff samplers.
Theorem 1 implies that a sequential knockoff sampling algorithm faithful to these two conditions is as general as it gets. The challenge now becomes creating exchangeable random variables at each step (with a little caution on the dependence on the previous pairs of variables). In turn, this task happens to be equivalent to designing a time-reversible Markov chain, as formalized below.
Proposition 1. A pair of random variables (Z,Z) is exchangeable, i.e., (Z,Z) d = (Z, Z), with marginal distribution π for Z-and, therefore, forZ as well-if and only if there exists a timereversible Markov chain {Z n } ∞ n=1 such that Z 1 ∼ π is a stationary distribution of the chain, and (Z 1 , Z 2 ) d = (Z,Z).
Combining these two results gives SCEP (Procedure 2 below), which is a completely general strategy for generating knockoffs: at each step j, we design a time-reversible Markov chain with stationary distribution L(X j | X -j ,X 1:(j−1) ), and draw a sample by taking one step of this chain starting from X j . Proposition 1 implies that the conditional exchangeability (4) holds. Furthermore, the symmetry requirement on the transition kernel implies that SCEP does not break the exchangeability from previous steps; that is, the knockoff symmetry (5) also holds. Theorem 1 then implies that such a procedure produces valid knockoffs.

Procedure 2: Sequential Conditional Exchangeable Pairs (SCEP)
for j = 1 to p do SampleX j by taking one step of a time-reversible Markov chain starting from X j . The transition kernel must be such that it depends only on X (j+1):p and the unordered pairs {X 1 ,X 1 }, . . . , {X j−1 ,X j−1 }, and admits L(X j | X -j ,X 1:(j−1) ) as a stationary distribution. end To rehearse the universality of SCEP, consider an arbitrary knockoff sampler producingX 1 , . . . ,X p . Then from Theorem 1 we know that X 1 andX 1 must be exchangeable conditional on X -1 . Therefore, X 1 may be sampled by taking one step of a reversible Markov chain starting at X 1 . Moving on to X 2 , Theorem 2 informs us that X 2 andX 2 are exchangeable conditional on {X 1 ,X 1 }, X 3 , . . . , X p , sõ X 2 can again be viewed as taking one step of a reversible Markov chain starting at X 2 . Continuing in this fashion for j = 3, . . . , p establishes our claim.
SCEP as stated remains too abstract to be considered an implementable algorithm, so we will next develop a concrete version of this procedure. Although this may not yet be clear, we would like to stress that formulating a knockoff sampler in terms of reversible Markov chains is an important step forward because it will ultimately enable the use of flexible MCMC tools.

The Metropolized knockoff sampler
We now demonstrate how one can generate knockoffs in a sequential manner by making proposals which are either accepted or rejected in a Metropolis-Hastings-like fashion as to ensure pairwise exchangeability.

Algorithm description
The celebrated Metropolis-Hastings (MH) algorithm (Metropolis et al. 1953;Hastings 1970) provides a general time-reversible Markov transition kernel whose stationary distribution is an arbitrary density function π. To construct a transition from x to y, MH operates as follows: generate a proposal x * from a distribution q(· | x) (any distribution depending on x) and set c y = x * with prob. α, This can be implemented even when the density π is unnormalized, as the normalizing constants cancel. In our setting, we shall make sure that the choice of the proposal distribution depends on the previously sampled pairs in a symmetric fashion, thereby remaining faithful to the knockoff symmetry condition (5) in Theorem 1. As such, we call such proposals faithful. Consider now running SCEP (Procedure 2) with the MH kernel, where at the jth step, the target distribution π is taken to be L(X j | X -j ,X 1:j−1 ). The issue with such a naïve implementation is that the target π cannot be readily evaluated. To understand why this is the case, set j = 2 and consider c More generally, we take as acceptance probability γ α with γ ∈ (0, 1]. In this work, γ is set to 1 as default, except in Section 3.3 and Appendix F.2, which are cases where tuning γ is recommended. L(X 2 | X -2 ,X 1 ). This distribution has density proportional to P(X = x)P(X 1 =x 1 | X = x), which is equal to The first term in the summation within the brackets corresponds to the acceptance case while the second corresponds to the rejection case. This latter term cannot be evaluated because of the integral over x * . Hence, the target density cannot be evaluated either.
We propose an effective solution to this problem: condition on the proposals and at step j, let the target distribution be L(X j | X -j ,X 1:j−1 , X * 1:j−1 ) rather than L(X j | X -j ,X 1:j−1 ). This has the effect of removing the integral and makes computing the rejection probability tractable. This is best seen by returning to our example where j = 2. Here, L(X j | X -j ,X 1:j−1 , X * 1:j−1 ) has density now proportional to We will show in Section 4 how such terms can be efficiently computed. Leaving aside implementation details for the moment, this strategy leads to Algorithm 1. Here and elsewhere, P denotes the density of the variables under study, or formally, the Radon-Nikodym derivative with respect to a common dominating measure.
for j = 1 to p do Sample X * j = x * j from a faithful proposal distribution q j . Accept the proposal with probability min 1, At this point, it should be clear that Metropolized knockoff sampling generates exact knockoffs, a fact we formally record below.
Proof. For the sake of the proof, let U j be the indicator of acceptance at step j, and Z j = (1−U j )X * j . We will prove pairwise exchangeability jointly with the U j 's and Z j 's; marginalizing out these variables will establish the claim. For 1 ≤ j ≤ p, let f j (x j , x -j ,x 1:(j−1) , u 1:(j−1) , z 1:(j−1) ) be the joint density function of (X j , X -j ,X 1:(j−1) , U 1:(j−1) , Z 1:(j−1) ), in this order. We will use induction to show that the density of (X,X 1:j , U 1:j , Z 1:j ) is symmetric in X k andX k for 1 ≤ k ≤ j. For 1 ≤ j ≤ p, the inductive hypothesis is that f j is symmetric in x k andx k for 1 ≤ k ≤ j − 1 (since f j is just the density of (X,X 1:(j−1) , U 1:(j−1) , Z 1:(j−1) ) after reordering the variables). For 1 ≤ j ≤ p, the density of (X,X 1:j , U 1:j , Z 1:j ) at (x,x 1:j , u 1:j , z 1:j ) = f j (x j , x -j ,x 1:(j−1) , u 1:(j−1) , z 1:(j−1) )× which is symmetric in the first j − 1 pairs by the inductive hypothesis. For the symmetry in the jth pair, when u j = 1, the density simplifies to which is invariant to swapping x j andx j ; when u j = 0, the delta function δ(x j − x j ) ensures x j =x j , and thus swapping them has no effect. Hence, when the algorithm terminates, all pairs are exchangeable and therefore remain exchangeable after marginalizing out the U j 's and Z j 's.
Anticipating possible future applications, we wish to remark that Metro can be easily adapted to sampling group knockoffs (Dai and Barber 2016); see Appendix E.

Covariance-guided proposals
Now that we have available a broad class of knockoff samplers, we turn to the question of finding faithful proposal distributions that will generate statistically powerful knockoffs. The overall challenge is to propose samples that are far away from X to make good knockoffs, but not as far that they are systematically rejected. A rejection at the jth step givesX j = X j , leading to a knockoff with poor contrast. Below, we shall borrow ideas from existing knockoff samplers for Gaussian models to make sensible proposals.
Suppose that X has mean µ and covariance Σ, and consider s ∈ R p with non-negative entries such that Γ(s) from (2) is positive semidefinite. Such a vector s can be found with techniques from Barber and Candès (2015) and from Candès et al. (2018). We have seen that if X were Gaussian, this covariance matrix would induce a multivariate Gaussian joint distribution over X andX with the correct symmetry. In non-Gaussian settings, our observation is that we can still make proposals as if the variables were Gaussian, but use the MH correction to guarantee exact conditional exchangeability. This can be viewed as a Metropolis-adjustment to the second-order knockoff construction of Candès et al. (2018). Concretely, the distribution q j for a covarianceguided proposal-used to generate a proposal X * j -is normal with mean and variance Γ here, X * 1:(j−1) is the sequence of already generated proposals, Γ (j) 11 = Γ 1:(p+j−1),1:(p+j−1) , Γ 22 = Γ p+j,p+j , Γ (j) 12 = Γ 1:(p+j−1),p+j , µ is the mean of X, and † stands for the pseudoinverse. The parameters of q j can be efficiently computed using the special structure of Γ; see Appendix D. The faithfulness of the proposal is shown in Appendix B.
The covariance-guided proposals are valid even when Σ is replaced by any other positive semidefinite matrix-any faithful proposal distribution will give valid knockoffs. This allows us to use an empirical estimate of Cov(X) based on simulated samples from L(X), or even to apply the covariance-guided proposals to discrete distributions by rounding each proposal X * j to the nearest point in the support of X j . These proposals will be most successful when X is well-approximated by a Gaussian density, indeed when X is exactly Gaussian and the true covariance is used, the covariance-guided proposals will never be rejected. Numerical simulations in a variety of settings can be found in Section 5.

Multiple-try Metropolis
A possibility for samplingX j "far away" from X j is to run multiple MH steps instead of a single one. The issue with this is that this would make the conditional distributions from Metro prohibitively complex at later steps. Longer chains also require conditioning later proposals on a longer sequence of proposals and acceptances or rejections, which will constrain those proposals to be closer to their corresponding true variables and thus reduce power. Instead, we use the multiple-try Metropolis (MTM) technique introduced in Liu et al. (2000).
The key idea of MTM is to propose a set of several candidate moves in order to increase the probability of acceptance. As in Qin and Liu (2001), we take the candidate set to be C m,t x = {x±kt : 1 ≤ k ≤ m}, where m is a positive integer and t is a positive number; see Figure 1 for an illustration. MTM proceeds by choosing one element x * from the set C m,t x , with probability proportional to the target density, i.e., This proposal is then accepted with probability γ min 1, , γ ∈ (0, 1), where γ is an additional tuning parameter explained in Appendix F.2. This parameter should be taken to be near 1 in most settings. If no element of C m,t x has positive probability, then one automatically rejects. MTM is a special case of MH with the proposal q(x * | x) distribution defined implicitly by the above rules, and furthermore, the proposals are faithful. Thus, MTM can be used in Metro. While there is no universally optimal combination of m and t, we provide guidance about default values based on our experimental results from Section 5. To understand the choice of parameters, first observe that with a fixed t, the knockoff distribution produced by the algorithm should eventually stabilize as m grows to infinity, since (non-pathological) π will vanish at positive and negative infinities and equations (8) and (9) will converge. Large values of m require more density evaluations, so we would like to choose the smallest value of m such that the distribution defined by equations (8) and (9) is nearly converged to its limit as m → ∞. Turning our attention to t, smaller values cause higher acceptance rates, and at the same time, encourageX j to be close to X j . Clearly, there is a trade-off. Based on our experiments, a sensible default setting is m = 4 and t j = 1.5 1/(Σ −1 ) jj where Σ = Cov(X). In the Gaussian case, Var(X j |X -j ) = 1/(Σ −1 ) jj for any observed value of X -j (Anderson 2009), hence this choice of scaling is intuitive. In the non-Gaussian case 1/(Σ −1 ) jj should be viewed as an approximation to the conditional variance. We have found that this parameter setting achieves nearly the best performance in most of our experiments.

Graphical models and conditional independence
One outstanding issue is whether the Metropolized knockoff sampler can be run in reasonable time for cases of interest. We begin by showing why sequential knockoff sampling is prohibitively expensive without additional structure, and then turn our attention to a common type of structure that enables efficient sampling: graphical models. The central contribution of this section will be a complexity bound on Metro showing how the graphical structure affects the difficulty of sampling. To complete this line of investigation, we give a complexity lower bound for all knockoff samplers which shows that Metro is optimal in some cases.

Why do we need structure?
Consider running Metro for some input distribution P and sample X = x. In view of (7), at step j we need to evaluate P(X -j , X j = z j ,X 1:(j−1) , X * 1:(j−1) ) for z j ∈ {x j , x * j } up to a constant. d Metro defines a joint distribution on (X,X 1:(j−1) , X * 1:(j−1) ) implicitly, so the only way to evaluate this density is to compute it step by step, from 1 to j − 1, i.e., through the sequential decomposition P(X -j , X j = z j ,X 1:(j−1) , X * 1:(j−1) ) = P(X -j , X j = z j )× j−1 k=1 P(X k | X -j , X j = z j ,X 1:(k−1) , X * 1:k )P(X * k | X -j , X j = z j ,X 1:(k−1) , X * 1:(k−1) ) . (10) Consider the term P(X k | X -j , X j = z j ,X 1:(k−1) , X * 1:k ). By the definition of Metro, computing this term will require evaluating an acceptance probability of the form Now, to compute the terms in the acceptance probability, we must use the same sequential decomposition (10) for the terms P(X -(j,k) , X k = z k , X j = z j ,X 1:k−1 , X * 1:k−1 ) for z k ∈ {x k , x * k }. Considering k = j − 1, we see that step j is making two calls to the probability at step j − 1, each of which is in turn making two calls to the probability function at step j − 2 and so on. Thus, each evaluation of (10) will require Ω(2 j ) function calls. This behavior is not due to a shortcoming of Metro; any genuine knockoff sampler with access only to an unnormalized density will require time exponential in p. We will present the formal statement of this lower bound later in Theorem 3. d In this section, when not explicitly specified, a variable is set to its observed value, e.g., Although knockoff sampling with no restriction on the distribution is prohibitively slow, we will show how to avoid the exponential complexity when there is additional known structure. Consider a Markov chain, i.e., a density that factors as P(x) = p−1 j=1 φ j (x j , x j+1 ). In this case, the joint density (10) can be evaluated efficiently provided we proceed along the chain in the natural order. Assume for simplicity that the proposal distribution is fixed in advance so that the second term within the square brackets in (10) does not depend on any variables and can be ignored. Due to the Markovian structure, only the k = j − 1 term in the product depends on z j , so it suffices to compute the acceptance probability (11) for k = j − 1. Again using the Markovian structure, this simplifies to where a j−1 is the ratio which does not depend on z j by the Markov structure. The key here is that a j−1 was previously computed with z j = x j when samplingX j−1 . Thus, the acceptance probability can be computed in constant time. Putting this all together, for a Markov chain, each of the necessary joint probabilities (10) can be computed in constant time, and the time to sample the entire vectorX is linear in the dimension p. Markov chains are not the only case where computing the acceptance probability can be done quickly; for other distributions with conditional independence structure, we next develop a systematic way of computing (10), using the graphical structure to control the depth of the recursion and hence control the running time.

Time complexity of Metro for graphical models
We have seen that we must restrict our attention to a subset of distributions in order to efficiently sample knockoffs, so in this section we show how to implement Metropolized knockoff sampling for a very broad class of distributions: graphical models. Let X ∈ R p be a random vector whose density factors over a graph G: here, C is the set of maximal cliques of the graph G and Φ is an unnormalized version of P. The variables in X can be either discrete or continuous. All graphical models with positive density or mass take this form (Hammersley and Clifford 1971), and such distributions are known to have particular conditional independence properties. We refer the reader to Koller and Friedman (2009) for a general treatment.
In order to take advantage of the conditional independence structure of X, we use a graphtheoretic object known as a junction tree (Bertele and Brioschi 1972) which encodes properties of the graph G.
2. For every edge (j, k) in G, j ∈ V and k ∈ V for some vertex V .

(Running intersection property)
If the vertices V and V both contain a node of G, then every vertex in the unique path from V to V also contains this node. Figure 2 gives an example of a junction tree over a 2 × 3 grid. The size of the largest vertex of T minus one is known as the width of the junction tree T , and the smallest width of a junction tree over G is called the treewidth of G, a measure of graph complexity. Finding the junction tree of lowest width for a graph G is known to be NP-hard (Arnborg et al. 1987), but there exist efficient heuristic algorithms for finding a junction tree with small width (Kjaerulff 1990;Koller and Friedman 2009).
Given a junction tree T for the graph G, we will soon prove that Metro can be run with O(p2 w ) queries of the unnormalized density Φ, where w is the width of T . In view of (7), at step j of Metro we need to evaluate P(X -j , X j = z j ,X 1:(j−1) , X * 1:(j−1) ) for z j ∈ {x j , x * j } up to a constant as well as sample from and evaluate the proposal distribution q j (· | x j ). We can use the graphical model structure to make these operations tractable by both (1) sampling the variables in a specific order, and (2) choosing proposal distributions that are not unnecessarily complex. We formalize these two requirements below.
We first consider the order in which we sample the variables. Recalling (10), the complexity of the computations of P(X -j , X j = z j ,X 1:(j−1) , X * 1:(j−1) ) depends on the number of function calls implied by the recursion (10). For simplicity, assume that the proposal terms in the product, P(X * k | X -j , X j = z j ,X 1:(k−1) , X * 1:(k−1) ), never depends on z j ; this will be relaxed soon. In that case, we need only consider the terms in (10) of the form P(X k | X -j , X j = z j ,X 1:(k−1) , X * 1:k ) for k < j. When there is graphical structure, not all such terms will depend on z j , and the number of terms that do depend on z j determines the recursion depth. In particular, if at step j only r j terms depend on z j , then there will be O(2 r j ) function calls in the recursion. A desirable ordering of the variables is then one that minimizes the largest r j , and such an ordering can be extracted from a junction tree T using Algorithm 2.
Algorithm 2: Junction tree variable ordering for Metro Select a leaf node V of T active . V is connected to at most one other node V of T active because it is a tree. In any order, append each j ∈ V \ V to the end of the list J. If no V exists, append all j ∈ V to J in any order. Remove V from the active tree T active . end Return J Algorithm 2 is valid in that when a node is removed, no j ∈ J remains in any node in T active . e From now on we assume that the variables are numbered according to this ordering. Our second consideration is to create proposals that do not add unnecessary complexity. No matter which proposal distribution we choose, P(X j = z j | X -j ,X 1:(j−1) , X * 1:(j−1) ) will still depend on some X for > j due to dependencies among coordinates of X; we however restrict ourselves to proposal distributions that do not add any additional dependencies.
Definition 2 (Compatible proposal distributions). Let V j be the node of the junction tree when j is appended to J from Algorithm 2. SetV j = {1, . . . , j − 1} ∪ V j . We say that proposal distributions q j are compatible with a junction tree T if they depend only on XV j ,X 1:(j−1) , and X * 1:(j−1) . This definition is motivated by the property Thus, a proposal distribution at step j that violates the compatibility property and relies on X for some / ∈V j will result in additional non-one terms in the product in (10) at step , soV j is the largest set that the proposal can be allowed to depend on without increasing the number of function calls/runtime. Although not all proposals are compatible, it is a rich enough class to handle a broad range of knockoff distributions, including the distribution induced by SCIP.
With these two conditions in place, we now state our main result about the efficiency of knockoff sampling, giving an upper bound on the number of evaluations of the unnormalized density function Φ that is required by Metro when the graphical structure is known. Assuming the variable ordering from Algorithm 2 and faithful proposal distributions compatible for T such that sampling from and evaluating the proposal distributions does not require evaluating Φ, we reach the following result: Theorem 2 (Computational efficiency of Metro). Let X be a random vector with a density which factors over a graph G as in (12). Let T be a junction tree of width w for the graph G. Under the conditions above, Metro uses O(p2 w ) queries of Φ.
This result means that we can efficiently implement Metropolized knockoff sampling for many interesting distributions, and it shows precisely how the complexity of the conditional independence structure of X affects the complexity of the sampling algorithm. Furthermore, in the next section we will prove that this is the optimal complexity in some cases.

Time complexity of general knockoff sampling
In the previous section we analyzed the runtime of Metro and showed that it will be tractable for graphs of sufficiently low treewidth. Now, we investigate the computational complexity of knockoff sampling in general. To formalize our investigation, we discuss a model of computation in which we have no information about the distribution of X beyond its graphical structure and the ability to query its (possibly unnormalized) density at any given point.
Oracle model. In this model, we are given as inputs (a) a p-dimensional vector X drawn from a density λΦ, where λ is a (possibly unknown) positive scalar so that we can think of Φ as an unnormalized density, (b) the support of Φ, and (c) a black box capable of evaluating Φ at arbitrary query points, and (d) a graph G for which the density is known to have the form (12). No other information about Φ is available.
We show that in the oracle model with the complete graph, i.e., when there is no graphical structure, knockoff sampling requires exponential time in the number of covariates, p. Please note that any complexity bound must take into account the quality of the generated knockoffs sincẽ X = X is a trivial knockoff that can be sampled in no time.
Theorem 3 (Complexity lower bound for knockoff sampling). Consider a procedure operating in the oracle model which makes a finite number of calls to the black box Φ and returnsX, thereby inducing a joint distribution (X,X) obeying the pairwise exchangeability (1) for all Φ. If G is the complete graph so that the procedure generates valid knockoffs for any input density, then the total number N of queries of Φ must obey N ≥ 2 #{j:X j =X j } − 1 a.s..
This result means that for any knockoff sampler, we cannot have both full generality and time efficiency. Put differently, in order to efficiently generate nontrivial knockoffs, we will need to restrict our attention to a subset of distributions for which we have structure. This fact justifies our decision to focus on distributions with graphical structure. We also derive a lower bound for the complexity of knockoff sampling for graphical models, stated next.
Corollary 2 (Complexity lower bound for graphical models). Consider the setting of Theorem 3. Fix a graph G with maximal cliques C. Suppose that for all Φ of the form Φ(x) = c∈C φ c (x c ), the procedure induces a joint distribution (X,X) obeying pairwise exchangeability (1). Then N ≥ max c∈C 2 #{j∈c:X j =X j } − 1 a.s..
This proposition shows that even after making some useful structural assumptions, there is still a trade-off between knockoff quality and computation. We next derive a byproduct, which proves that Metro is achieving a good runtime.
Proposition 2 (Optimality of Metro for chordal Gaussian graphical models). Consider continuous distributions of the form Φ(x) = c∈C φ c (x c ) over a chordal graph G. f On the one hand, for any input, Metro can be run with O(p 2 + p2 w ) queries of Φ. Furthermore, in the case where the distribution is Gaussian with zero mean and positive definite covariance (i.e., Φ(x) ∝ exp −xΣ −1 x /2 ), Metro can produce knockoffs with X j =X j for all j with probability 1. On the other hand, any general procedure that samples knockoffs such that X j =X j for all j with probability > 0 will require at least 2 w − 1 queries of Φ with probability at least .
Proposition 2 means that for chordal graphs, any general knockoff sampling algorithm such that P(X j =X j for all j) is bounded away from zero needs, in expectation, the same exponential order of queries as Metro (with the proviso that p is negligible compared to 2 w ).

Divide-and-conquer to reduce treewidth
Theorem 2 shows that Metro enables efficient computations for random vectors whose densities factor over a graph G of low treewidth. Not all graphs corresponding to random vectors of interest have low treewidth, however. A d 1 × d 2 grid, for example, has treewidth min(d 1 , d 2 ) (Diestel 2018). This section develops a mechanism for simplifying the graphical structure of a random vector X, allowing for faster computation of exact knockoffs at the cost of reduced knockoff quality.
To simplify graphical structure, we fix a set of variables C that separates the graph G into two subgraphs A and B. After fixing the variables in C, knockoffs can be constructed for the variables in A and B independently.
Proposition 3 (Validity of divide-and-conquer knockoffs). Suppose the sets A, B, C form a partition of {1, . . . , p} such that C separates A and B in the graph G, i.e., there is no path from some j ∈ A f A chordal graph is a graph such that any cycle of length 4 or larger has a chord. Figure 3: Two examples of conditioning to reduce the treewidth of a 6 × 6 grid from 6 to 3.
to some k ∈ B in G that does not contain some ∈ C. SupposeX is a random vector such that X C =X C a.s. and for all j A ∈ A and j B ∈ B, Furthermore, assume we construct the knockoffs for A and B separately, The divide-and-conquer technique can be applied recursively to split the graph into components of low treewidth until the junction-tree algorithm for constructing knockoffs can be used on each component. For example, for an arbitrary planar graph with p nodes, the planar separator theorem gives the existence of a subset of nodes C of size O( √ p) that separates the graph into components A and B with max(|A|, |B|) ≤ 2p/3 (Lipton and Tarjan 1979), suggesting that this technique will apply to many cases of interest. Figure 3 illustrates this technique for a d 1 × d 1 grid. We split the grid into rectangular ribbons of size d 1 × d 2 for small d 2 ; each resulting ribbon has treewidth d 2 .
The drawback of this approach is that for j ∈ C, we shall have X j =X j . When we think of deploying the knockoff framework in statistical applications, one should remember that we will work with multiple copies of X corresponding to distinct observations. We can then choose different separator sets for each observation so that in the end, X j =X j for most of the observations. For example, in the setting of Figure 3, one would randomly choose between the two choices of C for each observation. This technique is explored numerically in Section 5.2.3.

Discrete distributions
For discrete distributions with a small number of states for each coordinate X j , the junction tree techniques from Section 4.2 can be directly applied without using Metropolized knockoff sampling. When each variable X j can take on at most K values, the probability mass function P(X j | X -j ,X 1:(j−1) ) can be represented as a vector in R K , so at step j of the algorithm we simply need to evaluate P(X j = z j , X -j ,X 1:(j−1) ) for z j in the support of X j . This is the same quantity we computed in Section 4.2; see, e.g., (10). Once these probabilities have been computed, sampling from the resulting multinomial probability gives the SCIP procedure. In principle, this can be viewed as a special case of Metro, but for a practical implementation it is simpler to work directly with the probability vectors. A similar analysis to the proof of Theorem 2 then shows that the procedure requires O(pK w ) queries of the density Φ; see Appendix C.5 for details. For discrete distributions with infinite or large K, this is not tractable. However, Metro still applies and is much faster. Figure 4: An illustration of sampling knockoffs for an Ising model on a grid from Section 4.6. The blue dashed nodes represent the active variables of the junction tree when variable X 3,3 (shown in green) is being sampled. Gray nodes indicate variables that have already been sampled, and white nodes indicate variables that have not been sampled yet and are not in the active node of the junction tree.

Knockoffs for the Ising model
The tools from this section have the power to generate knockoffs for the Ising model on a grid (3). To construct an efficient knockoff sampler for this distribution, we need to find a junction tree of minimal width for the d 1 ×d 2 grid so that we can apply the technique from Section 4.5. A junction tree for the 2×3 grid of width 2 is shown in Figure 2, and the construction immediately generalizes to a junction tree of width min(d 1 , d 2 ) for the d 1 × d 2 grid, which is the optimal width. When d 1 ≥ d 2 , this leads to a knockoff sampler that proceeds from left to right, top to bottom; when variable X i,j is sampled, the other variables in the active node of the junction tree are X i,j+1:p and X i+1,1:j ; see Figure 4. Per our upper bound, this knockoff sampler will have runtime O(d 1 d 2 2 min(d 1 ,d 2 ) ). If min(d 1 , d 2 ) is large, this runtime may still be prohibitively long, but the divide-and-conquer technique from Section 4.4 greatly increases speed at the cost of slightly worse knockoffs than the impractical full procedure. We conduct a simulation experiment of both the small-grid and large-grid setting in Section 5.2.3.

Numerical experiments
We now empirically examine the Metropolized knockoff sampler, beginning with the few models where previously known samplers are available as a baseline, and then continuing on to cases with no previously known samplers. Condensed plots are presented in the main text, while more comprehensive versions can be found in Appendix F. We provide approximate runtimes with a singlecore g implementation in either R or Python. All source code is available from https://github. com/wenshuow/metro with interactive computing notebooks at http://web.stanford.edu/group/ candes/metro demonstrating the usage of the code and presenting further experimental results.

Measuring knockoff quality
The mean absolute correlation (MAC) is a useful measure of knockoff quality for a joint distribution of (X,X): We will use this as our measure of knockoff quality in our simulation experiments. Lower values of MAC are preferred. Let Γ be the correlation matrix of (X,X); pairwise exchangeability implies Γ g The hardware varies across simulations, but each CPU is between 2.5Ghz and 3.3Ghz.
is of the form (2). The MAC is then 1 p p j=1 |1 − s j |. Since Γ = Γ(s) has to be positive semidefinite, a lower bound on the MAC achievable by any knockoff-generation algorithm for a given distribution is the optimal value of the program This minimization problem can be solved efficiently with semidefinite programming (Barber and Candès 2015); we call the solution the SDP lower bound for the MAC. This lower bound can be achieved for Gaussian distributions (Candès et al. 2018). Valid knockoffs, however, must match all moments, not just the second moments, so this lower bound is not expected to be achievable in general; still it provides a useful goalpost in our simulations.

Gaussian Markov chains
We first apply our algorithm to Gaussian Markov chains and compare with the SDP Gaussian knockoffs, whose MAC achieves the SDP lower bound exactly, and SCIP knockoffs, both from Candès et al. (2018). We take p = 500 features such that X 1 ∼ N (0, 1) and X j+1 | X 1:j ∼ N (ρ j X j , 1 − ρ 2 j ). First, since the model is multivariate Gaussian, the covariance-guided proposal with s computed by the SDP method (14) will be identical to the SDP Gaussian knockoffs, so already a clever implementation of Metro is as good as a method specifically designed for Gaussian distributions, and since both achieve the SDP lower bound, one cannot do better in terms of MAC. Thus, we only investigate the MTM-proposals for implementing Metro. Note that the Gaussian knockoffs from Candès et al. (2018) do not use the Markovian structure of this problem, but instead rely on operations on 2p × 2p matrices, whereas the MTM knockoffs from this work utilize the Markovian structure to achieve time complexity linear in p.
The results are presented in Figure 5. Following Section 3.3, we vary the number of proposals and the step size. We find that choosing the step size for X j to be proportional to 1/(Σ −1 ) jj gives consistent results across different sets of ρ j 's. The MTM consistently outperforms the SCIP procedure, and is reasonably close to the SDP procedure. It is observed that the defaults from Section 3.3 of eight proposals (m = 4) and t j = 1.5 1/(Σ −1 ) jj performs nearly the best in all settings. Confirming our reasoning in Section 3.3, we find that the performance stabilizes as m grows and the step size should not be too large or too small, although for sufficiently large m the MAC is fairly stable to the choice of t. In this setting, it takes around 1 second for MTM to sample one knockoff vector with m = 4 and t j = 1.5 1/(Σ −1 ) jj .

Discrete Markov chains
For discrete Markov Chains there is one previously-known knockoff sampler, which is an implementation of the SCIP procedure (Sesia et al. 2018). We consider here Metro with MTM proposals. (The covariance-guided proposals would require ad-hoc rounding so we do not consider this here.) We take a simple Markov Chain with K ∈ {5, 10} states with uniform initial distribution and transition probabilities Q(j, j ) defined as Step size (t) We examine α from 0 (independent coordinates) to 0.5 (strong dependence between adjacent coordinates), with p = 500 features. We examine the MTM methods across a range of values of the tuning parameters, and the results are presented in Figure 6. Full simulation results are given in Appendix F. Note that the cases with K = 5 and α ≤ 0.15 are tuned with the additional parameter γ from (9), as detailed in Appendix F.2. We find that the best-tuned MTM method outperforms the SCIP method and achieves MAC near the lower bound for all dependence levels α. It takes around 0.5 seconds and 0.7 seconds respectively, to run MTM (m = 4 and t = 1) for K = 5 and K = 10.

Heavy-tailed Markov chains
As an example of a heavy-tailed distribution, we consider a Markov chain with t-distributed tails.
for j = 1, . . . , p = 500 where t ν represents the Student's t-distribution with ν > 2 degrees of freedom (note this is not a multivariate t-distribution). We try both the covariance-guided proposal with s provided by the SDP method (14) and the MTM proposals. We set ν = 5 and use the same ρ j 's as in the Gaussian setting. As in Section 5.1.1, a step size of 1.5 1/(Σ −1 ) jj again performs well. The covariance-guided proposals also perform well, although unlike the Gaussian case, there is now a gap between the lower bound and the performance of the covariance-guided proposals.
In this setting, it takes around 1.6 seconds for MTM to sample one knockoff vector with m = 4 (eight proposals) and t j = 1.5 1/(Σ −1 ) jj . For the covariance-guided proposals, it takes around 12.5 seconds for the one-time computation of the parameters (excluding time used for computing s, which varies depending on the method) and then 0.3 seconds to sample each knockoff vector. Step size (t)

Asymmetric Markov chains
As an example of asymmetric, continuous distributions, we take a standardized equal mixture of Gaussian and exponential random variables and then form a Markov chain. Explicitly, where Y G ∼ N (0, 1), Y E ∼ Expo(1) and I ∼ Bern(1/2) are independent. The parameters µ and σ are chosen so that Z j has mean 0 and variance 1. We then take X 1 = Z 1 , X j+1 = ρ j X j + 1 − ρ 2 j Z j+1 for j = 2, . . . , p.
We examine both the covariance-guided proposal with s provided by the SDP (14) and the multipletry proposals. We use the same ρ j 's as in the Gaussian setting. As in the previous case, m = 4 Step size (t) (eight proposals) and t j = 1.5 1/(Σ −1 ) jj performs essentially as well as any other MTM parameter choices, and significantly outperforms the covariance-guided proposals. The timing results are the same as in the heavy-tailed Markov chains.

Ising model
In this section, we consider an Ising model over a square grid (3). We generate knockoffs with the method for discrete random variables from Section 4.5 combined with the divide-and-conquer technique, the combination of which was described for Ising models in Section 4.6; no other exact knockoff samplers are known for the Ising model. Although our sampling procedures for the Ising model do not explicitly use the Metropolis-Hastings step, as explained in Section 4.5, we will refer to the sampler as "Metro" in this section for simplicity. First, we take a 10 × 10 grid and set all β i,j,i ,j = β 0 and all α i,j = 0. The results are presented in Figure 9. The left panel shows how the MAC increases-or, the quality decreases-as the dependence between adjacent variables-β 0 -increases. We see that the procedure is close to the lower bound for large β 0 . In the middle panel, we plot cor(X j,k ,X j,k ) across different coordinates (j, k). We see that on the edges of the grid, especially on the corners, knockoffs have lower correlation with their original counterparts. These variables are less determined by the values of the rest of the grid, so this is expected. In this setting, it takes about 12 seconds to sample a knockoff.
Next, we demonstrate the divide-and-conquer technique from Section 4.4. Here we consider the Ising model from above on a 100 × 100 grid, for a total dimension of 10, 000. The 100 × 100 grid has treewidth 100, so Metro would not be tractable without the the divide-and-conquer technique. We divide the graph into subgraphs of width w, by fixing entire columns as in Figure 3. To measure the effect of the slicing, we compute the MAC on the interior points and compare this to the MAC of the interior points of a smaller grid for a procedure without slicing, see Appendix F.3 for details. We find that the quality of the knockoffs increases as we take larger slices, as expected. Furthermore, even modest values of w such as w = 5 result in a procedure that achieves a MAC close to that of the baseline. Recall that the complexity of Metro scales as 2 w , so fixing w = 5 dramatically reduces the computation time compared to w = 100. With w = 5, it takes about 2.5 minutes to generate one knockoff for the 100 × 100 grid.

Gibbs measure on a grid
Lastly, we demonstrate the MTM proposals simultaneously with the junction tree techniques for complex dependence structure. Consider a Gibbs measure on {1, . . . , K} d×d , with a probability mass function and note that like the Ising model, this density factors over the grid. For our experiment, we take a 10 × 10 grid and examine different dependence levels β 0 with K = 20 possible states for each variable. We apply Metro with the MTM proposals and the divide-and-conquer technique on the grid, tuning the procedure across a range of parameters as detailed in Appendix F. The condensed results are given in Figure 10. We do not know of another knockoff sampler in this setting. Having said this, we observe that our procedure has MAC close to the lower bound. We also observe that in the case where w = 3, with as few as two proposals, our procedure performs well and takes about half a second to generate a knockoff copy; when we increase the number of proposals to ten, the compute time is around 2 minutes. When w is set to 5, the slowest setting is m = t = 1, which takes less than 4 minutes.

Discussion
This paper introduced a sequential characterization of all valid knockoff-generating procedures and used it along with ideas from MCMC and graphical models to create Metropolized knockoff sampling, an algorithm which generates valid knockoffs in complete generality with access only to X's unnormalized density. Although we proved in Theorem 3 that no algorithm (including Metro) can sample exact knockoffs efficiently for arbitrary X distributions, we characterized one way out of this impossibility result: conditional independence structure in X. An interesting future direction Step size (t) would be to establish other sufficient conditions on a model family that would allow one to sample knockoffs efficiently. Another way out of the lower bound in Theorem 3 is to forgo exact knockoffs and settle for approximations. Although this arguably is a tall order, it would be interesting to establish theoretical guarantees on the approximation quality of these or other approximate knockoff constructions, and better understand the tradeoff between knockoff approximation quality and time complexity.

A Junction tree lemmas
This section includes several important lemmas which will be used in some proofs in the appendix.
Lemma 1. If the variables are ordered by Algorithm 2, then for each 1 ≤ j ≤ p, any node in the junction tree that contains j is an element of the set {V 1 , V 2 , . . . , V j }. In addition, if j ∈ V k for some k > j, then V k = V j .
Proof of Lemma 1. According to Algorithm 2, when a node V is selected, all variables in V \ Vhere, V is the unique neighbor of V in the remaining junction tree-are sampled before the next node is selected. Recall that V j is the selected node when j is sampled; let k ≥ j be the last sampled variable when V j is selected (in this case, we have V j = V j+1 = · · · = V k by definition). After k is sampled, V j is removed, and by Algorithm 2, this means j, j + 1, . . . , k do not appear in V j 's only remaining neighbor. Now we claim no remaining node contains j. Otherwise, if there is a node V * j which contains j and still remains, by the running intersection property, all nodes on the unique path between V j and V * j contain j. This would imply that V j 's remaining neighbor in the junction tree also contains j since the unique path must pass through V j 's only remaining neighbor. This is a contradiction. Now we know that any node that contains j is some V with 1 ≤ ≤ k. Since V j = V j+1 = · · · = V k , the lemma follows.
Lemma 2. If the variables are ordered by Algorithm 2, then for any j > such that j ∈V , we haveV ⊆V j .
Proof of Lemma 2. Consider any k ∈V . Now we show such a k must be inV j . The case k ≤ j is trivial, since k ∈V j by definition. If k > j, then j, k ∈ V . Assume V = V j ; otherwise there is nothing to prove. Before j-and, therefore, k-are sampled, each time a node V containing j and k (e.g., V ) is selected, j and k appear in V 's neighbor. By a recursive argument, before V j is selected, each time the node containing j and k is selected, j does not get sampled and neither does k (k is sampled after j). Hence, before V j is selected, there is always at least one remaining node that contains both j and k. By Algorithm 2 and Lemma 1, V j is the last selected node j appears in, so it has to contain both j and k (otherwise no node contains both j and k at this point). This means that k ∈ V j ⊆V j .
Lemma 3. If the variables are ordered by Algorithm 2, then for any j = k, if j and k are connected in G, we have j ∈V k and k ∈V j .
Proof of Lemma 3. Without loss of generality, we assume k > j, so j ∈V k by definition. Now we show k ∈V j also holds. By the second property of the junction tree, k co-appears with j at least once in some node V . But j does not appear in any node after V j by Lemma 1, so there is some ≤ j such that {j, k} ⊆ V . If j = , then we already have k ∈ V = V j ⊆V j ; otherwise, j > , so by Lemma 2, k ∈ V ⊆V ⊆V j .

B Covariance-guided proposals
This section includes details on the covariance-guided proposal introduced in Section 3.2, and proofs of its faithfulness and compatibility. See Appendix D for details on how to do the necessary linear algebra computations for the covariance-guided proposals efficiently.
We first recall the definition of the covariance-guided proposals. Let Z be a 2p-dimensional vector drawn from N ((µ, µ), Γ(s)), where Γ(s) is as in (2). Let q j be the probability density function of Z p+j conditional on Z 1:(p+j−1) . The proposal distribution at the jth step q j (x * j | x 1:p , x * 1:(j−1) ), is defined to be the conditional density of Z p+j given Z 1:(p+j−1) = (x 1:p , x * 1:(j−1) ). h Proposition 4. The covariance-guided proposals are faithful (see Section 3.1) in that the proposal distribution at the jth step depends on (X k ,X k ) in a symmetric way for 1 ≤ k < j.
Proof of Proposition 4. Consider the proposal at step j. To see how the proposal distribution depends on (X k ,X k ) for k < j, note that we are using the distribution of Z j+p conditional on . (17) We only need to check if the proposal density changes when swapping x k andx k for k ≤ j − 1. Note that if we rejected at the kth step, x k =x k , so there is no effect of swapping the two; if we accepted at the kth step, the dependence is symmetric in (x k ,x k ) because of the structure of the covariance matrix.
Proposition 5. If (Σ −1 ) ij = 0 whenever i = j and (i, j) is not an edge in the graph G, then the covariance-guided proposals are compatible with G (see Definition 2).
Proof of Proposition 5. We wish to show, L(X * j | X, X * 1:j−1 ) only depends on X k if k ∈V j . To do this, we use induction over j. For j = 1, since X * 1 | X is a Gaussian distribution whose conditional variance does not depend on X, it suffices to show that E[X * 1 | X] depends on X k only if k ∈V 1 . Note that is a linear function of X 1 and µ cond . It is easy to see that µ cond = E[X 1 | X -1 ] depends on X k only if k and 1 co-appear in some node of the junction tree. By Lemma 1, this node can only be V 1 and, therefore, k ∈ V 1 . Thus the base case j = 1 holds.
Suppose the claim is true up to j − 1. By the same argument on the conditional distribution (X j , X * j ) | X -j , X * 1:(j−1) , it suffices to show that E[X j | X -j , X * 1:(j−1) ] only depends on X k if k ∈V j . We have h It seems equally plausible to use x1:p andx 1:(j−1) (i.e., (x1, x2, . . . , xp,x1,x2, . . . ,xj−1) instead of equation (17)). However, we find that its empirical performance is not as good as the version presented in the main text.
By the induction hypothesis, for < j, P(x * | x j , x -j , x * 1:( −1) ) does not depend on x j unless j ∈V (which implies j ∈ V since j > ), so the th term in the product can be removed from the numerator and the denominator if j / ∈ V . Thus, we now have Now we will prove that P(x j | x -j ) : <j, j∈V P(x * | x j , x -j , x * 1:( −1) ) depends on x k only if k ∈V j , which will conclude case j. Consider first the terms P(x * | x j , x -j , x * 1:( −1) ) in the product: such a term depends on x k only if k ∈V by the induction hypothesis, and by Lemma 2,V ⊆V j . Next, the term P(x j | x -j ) only depends on x k if k = j or k is connected to j in G. If k = j, k ∈V j by definition; otherwise, k ∈V j follows from Lemma 3.
We have established that as long as Σ reflects the structure of G, i.e., for i = j, Σ −1 ij = 0 only if i and j are connected in G, the covariance-guided proposals are compatible. For these proposals, sampling and evaluating the proposal density can be done without without querying the density Φ, so Theorem 2 implies that Metro with covariance-guided proposals requires O(p2 w ) queries of Φ.
It is easy to see that if X is Gaussian and γ = 1, we always accept because the acceptance ratio is always 1.

C.1 Necessity of the knockoff symmetry condition
At first glance, it might not be directly clear why we need the symmetry condition (5) in Theorem 1. To illustrate why this is necessary, consider the following example: let p be the density function Each pair (X j ,X j ) is unexchangeable unless j = p. However, marginalizing out any coordinate would yield the uniform distribution. In other words, in any sequential construction, we would have that X j andX j are conditionally independent and, therefore, exchangeable up until the last step. Since X p andX p are exchangeable conditionally on everything else, conditional exchangeability would hold. This example shows that if we require conditional exchangeability only, we would not necessarily end up with valid knockoffs. To press this point further, imagine that in the SCIP algorithm, we change the last step: instead of conditional independence we simply require conditional exchangeability. Then we are not guaranteed to get valid knockoffs. Violation of the symmetry condition in just one step is, in general, not allowed.
Assume now that (4) and (5) hold. We prove by induction that (18) holds for j = 1, 2, . . . , p; when j = p, we achieve pairwise exchangeability (1). When j = 0, there is nothing to prove. Assume (18) holds up until j − 1. The distribution of (X 1 , . . . , X p ,X 1 , . . . ,X j ) can be decomposed into the marginal distribution of (X -j ,X 1:(j−1) ) and the conditional distribution (X j ,X j ) | X -j ,X 1:(j−1) . The former is symmetric in X k andX k , 1 ≤ k ≤ j − 1 as seen by taking the induction hypothesis and marginalizing out X j . The latter is symmetric in X j andX j because of (4), and symmetric in X k andX k for 1 ≤ k ≤ j − 1 because of (5).
Proof of Proposition 1. Let Z 1 ∼ π, and the Markov kernel be given by the law L(Z | Z). Then the chain has π as a stationary distribution. Also, (Z 1 , Z 2 ) d = (Z,Z). Time reversibility also holds since The converse is a direct consequence of time reversibility.

C.3 Proof of Theorem 3
Proof. Suppose we are given a procedure K that always generates valid knockoffs for X given an unnormalized density function Φ for X and (implicitly) the support of Φ (or the dominating measure). Below, the symbols P Φ , L Φ , etc., indicate that we are working in the probability space defined by Φ (together with its support) and, implicitly, the procedure K; (X,X, N ) are jointly defined on this space. Let π be the normalized density, which is defined by Φ through We abuse notation slightly and for a Borel set M write Φ(X ∈ M ) for P Φ (X ∈ M )/λ Φ . We will consider the conditional probability P Φ (N ≥ 2 #{j:X j =X j } − 1 | X,X) and prove it is almost surely one. Conditioning on (X,X) makes the probability easier to analyze because it fixes the exponent #{j : X j =X j }. We will basically identify which 2 #{j:X j =X j } − 1 points have to be queried: any point obtained by changing x j tox j for j in any non-empty subset of {j : x j =x j }. We will prove the theorem for both discrete and continuous distributions. Analysis of the conditional probability can be done directly in the discrete case, while in the continuous case we cover the possible values of (X,X) by a countable union of sets, and then prove the conditional probability of interest is one on each of the sets. Although more technical, the proof for the continuous case shares the same structure as that for the discrete case.
Discrete case. Let (x,x) be some pair of input and output, respectively, of K. For any S ⊆ {1, 2, . . . , p}, define x ch(S) as x except with x j changed tox j for all j ∈ S, and vice versa forx ch(S) , so that (x,x) swap(S) = (x ch(S) ,x ch(S) ).
We will prove that as long as x ch(S) = x, then K must have queried the oracle at x ch(S) ; now assume x ch(S) = x. We also assume P Φ (X = x,X =x) > 0 (otherwise (x,x) is not a possible pair of input and output), which implies π(x ch(S) ) > 0 for any S ⊆ {1, 2, . . . , p} by pairwise exchangeability. Let . Also by pairwise exchangeability, Let A x be the event that the input vector is x, so Let Bx be the event that the output of K isx, so Let C S be the event that K does not query the oracle at x ch(S) . Dividing (19) by λ Φ gives Equation (20) holds for any Φ such that Φ(x), Φ(x) > 0. Consider a new unnormalized density where η ∈ (0, 1]. This new density has the same support as Φ, so Φ S η (x ch(S ) ) > 0 for any S ⊆ {1, 2, . . . , p}, and thus (20) also holds for Φ S η . Now consider The first probability does not depend on η because the conditioning on C S means changing the oracle only at x ch(S) does not affect the procedure K in any way. The second probability does not depend on η either, for the points that K queries can only depend on η after K queries the oracle at x ch(S) . Therefore, for any η ∈ (0, 1]. Thus, we get from equation (20) that (recall we are assuming Since Φ(x) > 0 (i.e., x is a possible input), we conclude by letting , given x as an input,x is a possible output), we conclude P Φ (C S | A x , Bx) = 0. That is, if K generatesx from input x, it must have queried x ch(S) . Thus, combining the results for all the S's that make x = x ch(S) , we can claim that given x as input, K outputsx only if it has queried the oracle at least at the set of points Mathematically, when P Φ (A x , Bx) > 0, Since there are 2 #{j:x j =x j } − 1 points in H (x,x) , we have After marginalizing out x andx, this leads to the claimed a.s. inequality: Continuous case. The proof works similarly for the continuous case. Loosely speaking, we will construct non-overlapping hypercubes around the points in H (x,x) defined previously, and show they all contain points of query with probability one. Concretely, consider a hypercube around z, defined as The denominator 3 in the definition of g is not essential, as long as it is large enough so F (z,z) does not overlap with F (z,z) swap(j) if z j =z j . Let E (z,z) be a joint hypercube around (z,z), which is defined through F (z,z) and F (z,z) as We use the fact that the rational points are dense in R 2p to cover the entire space using a countable collection of sets. Hence, if we can prove the conditional probability P Φ (N ≥ 2 #{j:X j =X j } − 1 | X,X) = 1 on every set in this collection, we can claim the corresponding equality holds unconditionally. Formally, we first consider E (r,r) , where (r,r) ∈ Q 2p and r k =r k , 1 ≤ k ≤ q for some positive integer q ≤ p. We will show P Φ (N ≥ 2 q − 1 | (X,X) ∈ E (r,r) ) = 1 as long as P Φ ((X,X) ∈ E (r,r) ) > 0. Now suppose P Φ ((X,X) ∈ E (r,r) ) > 0, which implies P Φ (X ∈ F (r,r) swap(S) ) > 0 for any S ⊆ {1, 2, . . . , p}. Define Let S be any non-empty subset of {1, 2, . . . , q}, so (r,r) swap(S) = (r,r). Now we have and P Φ ((X,X) ∈ E (r,r) ) = P Φ ((X,X) swap(S) ∈ E (r,r) ) = P Φ ((X,X) ∈ E (r,r) swap(S) ) = P Φ (X ∈ F (r,r) swap(S) ,X ∈ F (r,r) swap(S) ) = P Φ (X ∈ F (r,r) swap(S) )P Φ (X ∈ F (r,r) swap(S) | X ∈ F (r,r) swap(S) ) ≤ P Φ (X ∈ F (r,r) swap(S) ).
Now, similar to the discrete case, we consider a new unnormalized density for η ∈ (0, 1], which has the same support/dominating measure as Φ. It is easy to check that equation (22) also holds for Φ S η . Let C S be the event that K does not query Φ at any point in F (r,r) swap(S) . To use the same trick as in the discrete case, we next prove P Φ S η (C S ∩ B (r,r) | A (r,r) ) does not depend on η. We have Φ S η (X ∈ F (r,r) ) = Φ(X ∈ F (r,r) ) because F (r,r) ∩ F (r,r) swap(S) is empty (so Φ S η = Φ on F (r,r) ). Note that by definition, and P Φ S η (C S ∩ B (r,r) | X), as a function of X, does not depend on η; this holds since all Φ S η 's for η ∈ (0, 1] have the same support, and we can follow the same argument as in the discrete case. Specifically, P Φ S η (C S ∩ B (r,r) | X) = P Φ (C S ∩ B (r,r) | X). In addition, since the unnormalized density is the same on the set F (r,r) , we have A (r,r) ). The last two equations together with equations (23) imply P Φ (C S ∩B (r,r) | A (r,r) ) = P Φ S η (C S ∩B (r,r) | A (r,r) ). By (22), The left hand side of (24) goes to 0 as η → 0. Thus, recall that we are assuming P Φ ((X,X) ∈ E (r,r) ) > 0 and Φ(X ∈ F (r,r) ) > 0, and so we get P Φ (C S ∩ B (r,r) | A (r,r) ) = 0. By Bayes' rule, since P Φ (A (r,r) ∩ B (r,r) ) > 0, we have That is, unless (X,X) ∈ E (r,r) happens with zero probability, with probability one at least one point in F (r,r) swap(S) is queried conditional on (X,X) ∈ E (r,r) . Combining the results for all the S's that make (r,r) swap(S) = (r,r), we get 2 q − 1 disjoint sets, each of which must contain at least one point of query. Hence, now we can claim that for any E (r,r) , where (r,r) ∈ Q 2p and r k =r k , 1 ≤ k ≤ q, either P Φ ((X,X) ∈ E (r,r) ) = 0 or P Φ (N ≥ 2 q − 1 | (X,X) ∈ E (r,r) ) = 1.
Note that this is equivalent to P Φ (N ≥ 2 q − 1 | X,X) = 1 almost surely on E (r,r) . i These two equations imply that, as a function of (X,X), the conditional probability satisfies because the union is over a countable index set. There is nothing special about choosing the r k =r k on first q coordinates, so we have i The conditional probability PΦ(N ≥ 2 q − 1 | X,X) = 1 is almost surely one on a set U means PΦ(N ≥ 2 q − 1 | (X,X) = (x,x)) = 1 for (x,x) ∈ U \ V , where V is some set satisfying PΦ((X,X) ∈ V ) = 0.
Note the case D = ∅ does not follow the exact same proof, but no proof is needed in this case, since N ≥ 2 |∅| − 1 = 0 holds trivially. We shall thus keep in mind that (25) holds for any Φ. Now we go back to the conditional probability we are interested in, mathematically defined as We want to show f Φ (X,X) = 1, a.s. Let therefore, D is the set of coordinates where x andx could differ, and 1/n measures the minimum difference between these original and knockoff coordinates. Since any point ( This is also a countable union, so in order to show f Φ (X,X) = 1 a.s., we only have to show that f Φ (X,X) = 1 a.s. for any T n,D that has positive probability of containing (X,X). Note that since there are exactly |D| coordinates that differ for x andx in the set T n,D , So now we only need to show P Φ (N ≥ 2 |D| − 1 | X,X) = 1, a.s. on T n,D , which would be implied by (25) if we can show that (r,r)∈Q 2p r k =r k ,k∈D E (r,r) ⊇ T n,D .

C.4 Divide-and-conquer knockoffs
Proof of Proposition 3. Consider the distribution of X conditional on X C = x C . Since C separates A and B in the graph G, we have The assumptions of the proposition then imply that and sinceX C = X C a.s., Lastly, we note that conditional exchangeability implies marginal exchangeability, so as claimed.

C.5 Complexity proofs for Metropolized knockoff sampling
Lemma 4. When the proposal distributions are compatible for the junction tree T , for 1 ≤ j ≤ p, P(X j , X * j | X,X 1:(j−1) , X * 1:(j−1) ), depends on X k only if k ∈V j .
Proof of Lemma 4. We use induction over j. For the base case j = 1, P(X * 1 | X) depends on X k only if k ∈V 1 by our assumption of compatible proposals. And P(X 1 | X, X * 1 ) is a function of the acceptance probability The first term depends only on X k ∈V 1 by assumption of compatible proposals. The second term depends on X k if k is connected to 1 in G (or k = 1). Since the variables are ordered by Algorithm 2, 1 only appears in V 1 by Lemma 1, so k has to appear in V 1 if 1 and k are connected. The base case is thus proved.
Suppose the claim is true for 1, . . . , j − 1. First we have P(X * j | X,X 1:(j−1) , X * 1:(j−1) ) depends on X k only if k ∈V j by our assumption of compatible proposals. Now P(X j =x j | X,X 1:(j−1) , X * 1:j ) is a function of the acceptance probability, which is computed from the ratio of the proposal densities (which depends only on X k for k ∈V j by assumption of compatible proposals) and the following ratio P(X j = x * j , X -j ,X 1:(j−1) , X * 1:(j−1) ) P(X j = x j , X -j ,X 1:(j−1) , X * 1:(j−1) ) = P(X j = x * j | X -j ) P(X j = x j | X -j ) P(X 1:(j−1) , X * 1:(j−1) | X j = x * j , X -j ) P(X 1:(j−1) , X * 1:(j−1) | X j = x j , X -j ) . (26) We first consider the second term in the right hand side of the above. Consider the decomposition The th term in this product depends on X j only if j ∈V , by the induction hypothesis. Lemma 2 then implies that for such ,V ⊆V j . Any term in the product that does not depend on X j will be identical in the numerator and denominator of (26) and so will cancel. Together, this shows that the second term on the right hand side of (26) depends only on k for k ∈V j . Next, the numerator and denominator of the first term of the right hand side of (26) only depends on X k if k = j or k is connected to j in G. If k = j, then k ∈V j by definition; otherwise k ∈V j by Lemma 3. Now we have showed (26) depends only on X k for k ∈V j and the proof of the lemma is complete.
Proof of Theorem 2. Take γ = 1 for simplicitly, and take a proposal distribution q j that can be sampled from and evaluated without an evaluation of Φ (e.g., an independent Gaussian proposal in the continuous setting). We will show how Algorithm 2 uses the conditional dependence structure encoded in the graph G to make computations of (10) simpler.
Define 1) ). By the definition of Metro, we can write F j as the product of the proposal density and the acceptance/rejection probability: are the proposal terms and is the acceptance probability. Sequentially decomposing the ratio of probabilities in the term α, we get where all we did in the second equality was cancel the normalizing constants in the first ratio. In light of Lemma 4, the ratio is one unless j ∈V k , since the value of X j is the only one that differs between the numerator and denominator. Recall thatV k = V k ∪ {1, 2, . . . , k} and note j > k, so j ∈V k is equivalent to j ∈ V k .
Thus, (27) gives . Now we will show that the numerators in the product satisfy . By the definition of F k , we need to show that ). Inspecting the equation, we note that the only difference between the two quantities is the value of the variables that are being conditioned on, and only the values of X V j \V k are are different. Thus, we only need to show this set of values do not affect the conditional density. By Lemma 4, we just have to show that V j \ V k does not overlap withV k . To see this, take any ∈V k = V k ∪ {1, 2, . . . , k}. If ∈ V k , then certainly / ∈ V j \ V k . Now we consider the case where ∈ {1, 2, . . . , k} and thus less than j. If ∈ V j \ V k , then it must be in V j . By Lemma 1, we have V j = V , which means variables , . . . , j are all sampled when V is selected; specifically, k is sampled when V is selected, so V k = V = V j . But in this case certainly V j \ V k = ∅, which is a contradiction.
Similarly, we also have that the corresponding denominators in the product satisfy Combining all these together, the acceptance probability becomes α = min 1, Note that the only difference between c (1) and c (2) is changing x * j to z j . From this expression, we see that F j (X V j = z V j ) : z ∈ {x , x * } for all ∈ V j can be computed in terms of we have that the total number of queries of Φ is O(p2 w ). Having access to the F k (X V k = z V k )'s is sufficient to run the algorithm, because at each step of the algorithm, it is clear that the acceptance ratio α can be computed from these F k (X V k = z V k )'s, so the proof is complete.
More generally, if an evaluation of Φ costs a units of computation and a floating point operation requires 1 unit, then the same proof shows that the algorithm takes O(p(p + a)2 w ) since computing the F j (X V j = z V j ) requires a total of O(p2 w ) evaluations of Φ and each F j (X V j = z V j ) requires O(p) floating point operations to compute c (1) and c (2) .

MTM
For the MTM method of Section 3.3, the proposal distribution requires evaluations of Φ, so the requirements of Theorem 2 do not hold. The proof of Theorem 2 can easily be adapted to apply to MTM, however. Instead, for the MTM method we see that at step j we need access to By an analysis similar to the proof of Theorem 2, it suffices to have access to In order to compute F j (X V j = z V j ) for all z V j such that z ∈ C m,t x ∪ C m,t x * for all ∈ V j for use in later steps, we additionally need to compute Thus, MTM requires O(p(3m + 1) w ) queries of Φ, where 3m + 1 is an upper bound on |C m,t x ∪ C m,t x * | for all .

Discrete distributions with small support
Consider the direct methods for discrete distributions with small support in Section 4.5, which can be viewed as a Metro algorithm that never rejects. The proof of Theorem 2 can easily be adapted to apply to the this case. Note that since the procedure never rejects, we have X * =X and we can omit writing terms of X * in all of the following discussion.
Let C j be the support of X j and suppose that |C j | ≤ K for all j. Then at step j, the method requires access to P(X j = z j , X -j = x -j ,X 1:(j−1) =x 1:(j−1) ) for all z j ∈ C j . By an analysis similar to the proof of Theorem 2, it suffices to have access to In order to compute F j (X V j = z V j ) for all z V j such that z ∈ C for all ∈ V j for use in later steps, we additionally need to compute Thus, the rejection-free procedure requires O(pK w ) queries of Φ.

C.6 Lower bounds for graphical models
Proof of Proposition 2. We only have to prove that we can design a Metro algorithm that meets the requirement in the proposition. The rest of the proposition is implied by Corollary 2, proved next. We focus on the continuous case (and when the support is all of R p ), because we are interested in the Gaussian case. Since Metro can only learn about the distribution by making queries to the oracle, we describe an algorithm which first makes about p 2 /2 queries to attempt to recover the covariance matrix. This can be done in such a way that if the model is Gaussian as described in Proposition 2, we recover the correct covariance matrix. Thus, covariance-guided proposals will accept at every step and get us a knockoff vector that differs with the input original vector at every coordinate (with probability one). Let e j be the jth vector of the canonical basis. Assume the model is N (0, Σ) let Φ be the unnormalized density. Then where W is an unknown positive constant. Hence, if we query the oracle at the above 1 + p(p + 1)/2 points, we can always solve for a potential precision matrix. The algorithm's next step depends on the solution to these equations.
• If the matrix formed by the solution to this system of equations is positive definite, and reflects the structure of the graph G, i.e., the (i, j)th entry is non-zero only if i = j or i and j are connected by an edge in G, then the algorithm inverts this solution matrix to get Σ, and then proceeds with covariance-guided proposals with any positive s which makes Γ(s) in (2) positive definite.
• Otherwise, the model must not be a multivariate Gaussian distribution with zero mean, positive definite covariance matrix and have the required conditional independence structure. In that case, the algorithm will just choose any proposal distribution (e.g., independent Gaussian proposals).
Since running Metro with the covariance-guided proposal requires O(p2 w ) queries of Φ, this algorithm in total requires O(p 2 + p2 w ) queries of Φ. If indeed Φ(x) ∝ exp −xΣ −1 x /2 , the algorithm will recover the right covariance matrix Σ, and therefore will never reject, and produce a knockoffX such that X j =X j for all j (with probability one).

D Efficient matrix inversion for covariance-guided proposals
Let Σ j be the matrix composed of the first (p + j) rows and columns of We want to find the inverses of Σ 0 , Σ 1 , . . . , Σ p−1 (assuming Σ p−1 is invertible). Note that where σ 2 j+1 is the (j + 1)th diagonal element of Σ, and γ j+1 is the truncated (p + j)th column of Γ. We have And by the Sherman-Morrison formula, With all the elements of recursion in place, we can invert Σ 0 and recursively calculate the inverse matrices of Σ 1 , . . . , Σ p−1 . We have made code available that implements this recursion efficiently.

E Group knockoffs
We can easily generalize our work to the group knockoff filter first presented in Dai and Barber (2016) to control the group false discovery rate. As in that work, let {I 1 , I 2 , . . . , I k } be a partition of {1, 2, . . . , p}, and suppose we want to constructX such that for each j = 1, . . . , k, (X I 1 , X I 2 , . . . , X I k ,X I 1 ,X I 2 , . . . ,X I k ) d = (X I 1 , X I 2 , . . . , X I k ,X I 1 ,X I 2 , . . . ,X I k ) swap(I j ) .

F.2 Effect of γ
The tuning parameter γ for the MTM procedure introduced in Section 3.3 may appear mysterious to the reader and warrants an explanation. In most of our MTM simulations, γ is set to be 0.999, but there are cases where it is necessary to tune γ for improved performance. For example, in the discrete Markov chain experiment in Figure 11, we sometimes observe that MAC increases with the number of proposals. This is surprising, and upon closer inspection, we find that the reason is that many pairs (X j ,X j ) have negative correlations, which leads to increased MAC. In this case, what is happening is that we are proposing and accepting points that are so far away from X j that they become negatively correlated with X j , which is undesirable. To shift the negative correlations toward zero, γ can be decreased so that X j =X j more frequently. We illustrate this in Figure 13. For example, in the setting where have independent coordinates taking on K = 5 possible states, the best performance is obtained with γ = 4/5, since with this value of γ there will be a probability 1 − γ = 1/5 of rejection, which makes X j andX j independent. Tuning γ may also enable fewer rejections at later stages of the algorithms, since the knockoffs at later coordinates will be less constrained. j Based on our simulation results, we only recommend tuning γ when the variables are discrete with small support and the dependence between variables is weak.

F.3 Ising model simulation details
Sampling Ising variables X is done with a Metroplis-Hastings sampler implemented in the bayess R package.

F.3.1 Divide-and-conquer simulation details
We setX i 1 ,i 2 := X i 1 ,i 2 for all (i 1 , i 2 ) such that i 1 belongs to the set C of columns defined as C = {1 ≤ i ≤ 100 : i = a 0 + b(w − 1), b ∈ N}; the spacing w is a fixed constant (see Figure 3 for an illustration) and the offset a 0 is chosen uniformly from {2, . . . , w + 1}. This implies that for all sites i 1 , i 2 ), P(X i 1 ,i 2 =X i 1 ,i 2 ) < 1.
The SDP lower bound is not available in this case, because it would require computing a 10000× 10000 covariance matrix and then solving the SDP, which is intractable. Instead, to evaluate the quality of the knockoffs, we compare to Ising model knockoffs on a smaller grid that does not require the divide-and-conquer technique. Our baseline is thus the MAC evaluated at interior nodes 1 < i 1 , i 2 < 10 achieved by the SCIP procedure. We consider interior nodes because we recall that correlations on the edges of the grid are smaller. We compare this figure of merit to the MAC of the interior variables of the 100 × 100 grid. Without the divide-and-conquer strategy, the MAC of the two procedures would be very similar-hence, this is a sensible baseline.

F.4 Gibbs model simulation details
In Figure 10, the best MTM specification is taken from {(γ, m, t, w) : γ = 0.999, 1 ≤ m ≤ 5, 1 ≤ t ≤ 7, w = 3} ∪ {(γ, m, t, w) : γ = 0.999, m = 1, 1 ≤ t ≤ 10, w = 5} for β = 0.07, 0.1, 0.3, and from {(γ, m, t, w) : γ = 0.999, 1 ≤ m ≤ 5, 1 ≤ t ≤ 7, w = 3} ∪ {(γ, m, t, w) : γ = 0.999, m = 1, 1 ≤ t ≤ j Choosing γ less than 1 means that no matter what the proposal X * j is, it will be rejected with positive probability. While we typically want to avoid rejections, rejecting at early stages in the algorithm my lead to better performance by enabling higher quality knockoffs at later steps in the algorithm. In particular, with γ = 1, at some step k > j, it might be the case that none of the points in the proposal set have positive probability, because any point in the proposal set is inconsistent with a rejection that occurred previously in step j. When γ < 1, however, any proposed value at step k is consistent with a rejection at step j, because there is always at least a 1 − γ chance of rejecting at step j, so we avoid the undesirable situation described above. 10, w = 3, 5} for β = 0.003, 0.01, 0.02, 0.05. The best m = 1, w = 3 MTM for each β is taken from the same set intersecting {(γ, m, t, w) : m = 1, w = 3}.