# Learning generative models for valid knockoffs using novel multivariate-rank
based statistics

Shoaib Bin Masud1, Shuchin Aeron1 $1$ \- Tufts University, Dept. of ECE,
Corresponding authors: Shoaib Bin Masud, email<EMAIL_ADDRESS>
###### Abstract

We consider the problem of generating valid knockoffs for knockoff filtering
which is a statistical method that provides provable false discovery rate
guarantees for any model selection procedure. To this end, we are motivated by
recent advances in multivariate distribution-free goodness-of-fit tests
namely, the rank energy (RE), that is derived using theoretical results
characterizing the optimal maps in the Monge’s Optimal Transport (OT) problem.
However, direct use of use RE for learning generative models is not feasible
because of its high computational and sample complexity, saturation under
large support discrepancy between distributions, and non-differentiability in
generative parameters. To alleviate these, we begin by proposing a variant of
the RE, dubbed as soft rank energy (sRE), and its kernel variant called as
soft rank maximum mean discrepancy (sRMMD) using entropic regularization of
Monge’s OT problem. We then use sRMMD to generate deep knockoffs and show via
extensive evaluation that it is a novel and effective method to produce valid
knockoffs, achieving comparable, or in some cases improved tradeoffs between
detection power Vs false discoveries.

##  1 Introduction

Selection of significant features with control over the false discovery rates
(FDR), is one of the most important problems in application of machine
learning methods (Jović et al., 2015; Rietschel et al., 2018) for problems
such as untargeted metabolomics (Heinemann, 2019) and genome-wide association
study (GWAS) (Frommlet et al., 2012). Due to the large number of features
involved, it is of utmost importance to provide provable guarantees in
selecting the true underlying features that can explain certain phenotypic
conditions.

Classical methods of FDR control depend on the assumptions on how the features
and the responses are related (Benjamini and Hochberg, 1995; Gavrilov et al.,
2009). Barber and Candes in their seminal paper (Candes et al., 2016),
proposed a novel FDR control approach, called the Model-X knockoff that can be
used as a statistical wrapper around any machine learning method that can
select features. Model-X knockoff does not rely on the nature of the
relationship between the features and responses and therefore is model-free.
In order to control the FDR, the Model-X framework generates a synthetic set
of features called knockoffs, which mimic the original features but are
conditionally independent of the responses given the original features.

#### Related works.

Existing methods of knockoff generation either (i) assume the distribution of
the features, or (ii) incorporate a generative model to learn the feature
distribution from data. Second-order knockoff (Candes et al., 2016) assumes
that the distribution of the features is jointly Gaussian. Another knockoff
generation approach based on the Hidden Markov Model (HMM) (Sesia et al.,
2019) characterizes the feature distribution with the help a Markov chain.
Conditional independence knockoff (CIK) (Liu and Rigollet, 2019) can produce
valid knockoffs only if the Gaussian graphical model associated with the
features is a tree. A more flexible approach that samples knockoffs from a
Bayesian network is discussed in (Gimenez et al., 2019), where the covariates
are modeled as the observed variables of the network.

Methods such as KnockoffGAN (Salimans et al., 2016), Deep Knockoffs (Romano et
al., 2020), and Auto-Encoding knockoff (Liu and Zheng, 2018) focus on learning
a deep generative model to produce valid knockoffs. To create valid knockoffs,
KnockoffGAN employs a generative adversarial network (GAN) that requires
simultaneous training of four different but interconnected neural networks and
is therefore computationally expensive. Auto-Encoding knockoff uses latent
variables to reconstruct the covariate distribution using a variational
autoencoder (Kingma and Welling, 2019). The performance of this method depends
on the dimension of the latent space considered. Higher dimensional latent
space can be used to make the model better but may have diminished power if
the covariates violate the low-dimensional approximation. Deep-Knockoffs uses
the celebrated two-sample goodness-of-fit statistic called maximum mean
discrepancy (MMD) (Gretton et al., 2012) as the loss function in a generative
model to produce knockoffs. Though training a generative model with MMD is
less expensive (Li et al., 2015), in high dimension MMD has less power (Ramdas
et al., 2015). Another method, called DDLK (Sudarshan et al., 2020) generates
knockoffs by first maximizing the likelihood of the features and then
minimizing the KL (Kullback-Liebler) divergence between the joint distribution
of the features and the knockoffs and any possible swaps between them. KL
divergence cannot provide a useful gradient when the supports of the two
distributions are disjoint (Bińkowski et al., 2018), which is a long-
recognized problem for KL divergence-based discrepancy measures.

#### Summary of main contributions:

We introduce a new statistic called soft Rank Energy $(\text{sRE})$, which is
heavily inspired by the recent development in multivariate distribution-free
goodness-of-fit tests based on Rank Energy $(\text{RE})$ (Deb and Sen, 2019),
that in turn is based on the fundamental results in optimal transportation
theory (McCann et al., 1995). sRE is the extension of RE that is obtained by
entropic regularization of the optimal transport problem. Similar to (Cuturi
et al., 2019), this makes sRE, when used as loss function, a differentiable
function of the generative model parameters. We highlight on the properties of
sRE which make it a desirable candidate for measuring two-sample goodness-of-
fit. We also introduce the kernel variant of sRE called soft Rank Maximum Mean
Discrepancy $(\text{sRMMD})$. We inspect the behaviour of sRE and sRMMD w.r.t.
the sample size, dimensions and the regularization parameters. We show that
with an appropriate entropy regularizer, and sample size, sRMMD based
generative models do not suffer from mode collapse. We use sRMMD as a loss
function in a generative model to produce valid knockoffs. We demonstrate that
knockoffs generated by our proposed method can keep the FDR under control,
with either comparable or increased detection power compared to existing
baselines in the cases considered.

#### Notations:

We use bold-math capital letters $\bm{X}$ for multivariate random variables,
bold-face capital letters $\mathbf{X}$ for matrices and maps, lower-case bold
face bold-math $\bm{x}$ for vectors. We denote by
$\mathcal{P}(\mathbb{R}^{d})$\- the set absolutely continuous measure in
$\mathbb{R}^{d}$. $\overset{d}{=}$ refers to the equality in distribution.
$\mathcal{N}(\mu_{d},\mathbb{I}_{d})$ is referred to as the Gaussian
distribution, where $\mathbb{I}_{d}$ is a $d$-dimensional identity matrix. The
rest of the notations is standard and should be clear from the context.

##  2 Overview of the knockoff filter

Given that $\bm{X}=(X_{1},\dots,X_{d})\in\mathbb{R}^{d}$ and $y\in\mathbb{R}$
are the multivariate random variable and the response variable, respectively.
Assume that the distribution associated with $\bm{X}$ is known and denoted as
$\mathsf{F}_{\bm{X}}$, whereas there exists no knowledge about the conditional
distribution, $\mathsf{F}_{y|\bm{X}}$. A variable $X_{j}$ is said to be “null
or unimportant” if and only if $y$ is conditionally independent on $X_{j}$,
given the other variables, $y\perp\\!\\!\\!\perp X_{j}|X_{-j}$, where
$X_{-j}=\\{X_{1},\dots,X_{d}\\}\setminus X_{j}$. Otherwise, it will be
considered as a relevant or important variable. Under these assumptions, the
primary goal of the Model-X Knockoff filter (Candes et al., 2016) is to
discover as many relevant variables as possible while keeping the FDR under
control. FDR is defined as following,

$\text{FDR}=\mathbb{E}\left[\frac{|\hat{\mathcal{S}}\cap\mathcal{H}_{0}|}{\text{max}(1,\;\;|\hat{\mathcal{S}}|)}\right],$

where $\mathcal{H}_{0}$ denotes the true set of null or unimportant variables
and $\hat{\mathcal{S}}\subseteq\\{1,2,\dots,d\\}$ is a subset of variables
that is selected by any variable selection method.

To control FDR, Model-X knockoff (Candes et al., 2016) generates knockoff,
$\bm{\tilde{X}}$\- a random vector that satisfies the following properties:

  * (a)

Exchangeability:
$(\bm{X},\tilde{\bm{X}})_{\text{swap}(B)}\overset{d}{=}(\bm{X},\tilde{\bm{X}})$,

  * (b)

Conditional independence: $y\perp\\!\\!\\!\perp\tilde{\bm{X}}\;|\;\bm{X}$.

Exchangeability property ensures that for any subset
$B\subset\\{1,2,\dots,p\\}$, the joint distribution remains unchanged when the
variables and their corresponding knockoffs exchange their positions. That is,
for any random vector, $\bm{X}=(X_{1},X_{2},X_{3})$, and a set, $B=\\{2,3\\}$,
in order to satisfy the exchangeability condition, we must have
$(X_{1},\tilde{X}_{2},\tilde{X}_{3},\tilde{X}_{1},X_{2},X_{3})\overset{d}{=}(X_{1},X_{2},X_{3},\tilde{X}_{1},\tilde{X}_{2},\tilde{X}_{3})$.
The second property ensures that knockoffs are generated without the knowledge
of the response variable.

To use knockoffs in a controlled variable selection procedure, it is important
to construct a knockoff statistic, $W_{j}$ that can satisfy the flip-sign
property. That is, $W_{j}=w_{j}([\bm{X},\tilde{\bm{X}}],y)$, for some function
$w_{j}$ and $B$, must satisfy the following,

$\displaystyle
w_{j}([\bm{X},\tilde{\bm{X}}]_{\text{swap}(B)},y)=\begin{cases}w_{j}([\bm{X},\tilde{\bm{X}}],y)\;\text{if}\;j\not\in
B,\\\ -w_{j}([\bm{X},\tilde{\bm{X}}],y)\;\text{if}\;j\in B.\end{cases}$ (1)

Equation (1) says that, swapping the position of $j$-th variable with its
knockoff will change the sign of $W_{j}$.

Given the knockoff statistics $W_{j}$ for $j=1,\dots,d$, a guaranteed control
of FDR is achieved by the knockoff filter at a level $q\in(0,1)$ by selecting
$j$-th variable, $j\in\\{1,\dots,d\\}$, such that $W_{j}\geq\tau$, where
$\tau$ is achieved via the following,

$\tau=\min_{t>0}\Big{\\{}t:\frac{1+|\\{j:W_{j}\leq-t\\}|}{|\\{j:W_{j}\geq
t\\}|}\leq q\Big{\\}}.$ (2)

The power and FDR control of a knockoff filter depends on how the
exchangeability condition is satisfied. A naive approach called Second-order
knockoff (Candes et al., 2016) only approximates first two moments to satisfy
the exchangeability condition assuming that the covariates follow a
multivariate Gaussian distribution. However, Second-order knockoffs provide
insufficient guarantees for FDR control if the assumption does not hold. On
the other hand, Deep knockoffs (Romano et al., 2020) use a generative model
based on the maximum mean discrepancy (MMD) (Gretton et al., 2012) statistic
to produce higher-order knockoffs that can achieve provable FDR guarantee
under general conditions.

We now briefly describe the knockoff generation procedure used in Deep
Knockoffs (Romano et al., 2020).

###  2.1 Deep knockoffs

Deep knockoff (Romano et al., 2020) employs a deep neural network,
$f_{\theta}(\bm{X},\bm{V})$ that generates a knockoff
$\tilde{\bm{X}}\in\mathbb{R}^{d}$ from original input variable
$\bm{X}\in\mathbb{R}^{d}$, and a noise vector
$\bm{V}\sim\mathcal{N}(0,\mathbb{I}_{d})\in\mathbb{R}^{d}$. Here, $\theta$
denotes the parameters of the network. In order to satisfy the exchangeability
condition without any assumption on the covariate distribution, Deep knockoff
(Romano et al., 2020) uses the unbiased MMD (Li et al., 2015) estimate as a
loss function- that can match the higher-order moments including the first two
moments. In particular, for a design matrix $\mathbf{X}\in\mathbb{R}^{n\times
d}$ with $n$ observations, and the corresponding knockoff matrix
$\mathbf{\tilde{X}}\in\mathbb{R}^{n\times d}$, Deep knockoff computes the loss
in the following manner:

$\displaystyle\ell_{\text{MMD}}(\mathbf{X},\tilde{\mathbf{X}})$
$\displaystyle={\text{MMD}}\Big{[}(\mathbf{X}^{\prime},\tilde{\mathbf{X}^{\prime}}),(\tilde{\mathbf{X}^{\prime\prime}},\mathbf{X}^{\prime\prime})\Big{]}+{\text{MMD}}\Big{[}(\mathbf{X}^{\prime},\tilde{\mathbf{X}^{\prime}}),(\mathbf{X}^{\prime\prime},\tilde{\mathbf{X}^{\prime\prime}})_{\text{swap}(B)}\Big{]},$
(3)

where
$\mathbf{X}^{\prime},\mathbf{X}^{\prime\prime},\mathbf{\tilde{X}}^{\prime},\mathbf{\tilde{X}}^{\prime\prime}\in\mathbb{R}^{n/2\times
d}$ are obtained by randomly splitting the design matrix in half and $B$ is a
uniformly chosen random subset of $\\{1,\dots,d\\}$, such that $j\in B$ with
probability $1/2$. To increase the power of the knockoff filter, a
regularization term is added to the loss function that penalizes the large
pairwise correlation,

$\displaystyle\text{D}_{\text{corr}}(\mathbf{X},\tilde{\mathbf{X}})=\|\text{diag}(\hat{\mathbf{G}}_{\bm{X}\tilde{\bm{X}}})-1+s^{*}_{\text{SDP}}(\hat{\mathbf{G}}_{\bm{X}\bm{X}})\|^{2},$
(4)

where
$\mathbf{\hat{G}}_{\bm{X}\bm{\tilde{X}}},\mathbf{\hat{G}}_{\bm{X}\bm{X}}\in\mathbb{R}^{d\times
d}$ are the empirical covariance matrices and $s^{*}_{\text{SDP}}$ is the
solution to the following semidefinite programming,

$\displaystyle s^{*}_{\text{SDP}}(\Sigma_{\bm{X},\bm{X}})$
$\displaystyle=\arg\min_{s\in[0,1]^{d}}\sum_{j=1}^{p}|1-s_{j}|,\,\,\,\text{s.t.}\;\;2\Sigma\succeq\text{diag}(s).$

Deep knockoffs also adds a second-order loss term (Equation 7 in (Romano et
al., 2020)) which is claimed to be effective in reducing the training time.
Total loss of the generative model is given as,

$\displaystyle\ell(\mathbf{X},\mathbf{\tilde{X}})\\!\\!=\\!\\!\ell_{\text{MMD}}(\mathbf{X},\tilde{\mathbf{X}})+\lambda\ell_{\text{s-o}}(\mathbf{X},\mathbf{\tilde{X}})+\delta\text{D}_{\text{corr}}(\mathbf{X},\tilde{\mathbf{X}}).$
(5)

For further insights on the generative model, we refer the reader to (Romano
et al., 2020) and reference therein.

##  3 Background: OT based multivariate rank energy

Rank-based goodness-of-fit tests have been studied comprehensively in 1-D
e.g., Kolmogorov-Smirnov test (Smirnov, 1939), Wilcoxon signed-rank test
(Wilcoxon, 1947), Wald-Wolfowitz runs test (Wald and Wolfowitz, 1940). Unlike
1-D, due to the lack of canonical ordering in $d$\- dimensional space, for
$d\geq 2$, the notion of rank cannot be defined in a straightforward way.
Recently several authors (Hallin et al., 2017, 2021; Chernozhukov et al.,
2017; Deb and Sen, 2019; Hallin et al., 2020; Shi et al., 2020b, a) studied
the notion of multivariate rank based on OT theory. In this paper, we will
consider the setting in (Deb and Sen, 2019) explicitly and build upon the
ideas therein.

#### Ranks and Quantiles for univariate distributions.

Let $X$ be a univariate random variable with c.d.f.
$\mathsf{F}:\mathbb{R}\rightarrow[0,1]$. It is a standard result that when
$\mathsf{F}$ is continuous, the random variable
$\mathsf{F}(X)\sim\mathsf{U}[0,1]$ \- the uniform distribution on $[0,1]$. For
any $x\in\mathbb{R}$, $\mathsf{F}(x)$ is referred to as the rank-function. For
any $0<p<1$, the _quantile function_ is defined by
$\mathsf{Q}(p)=\inf\\{x\in\mathbb{R}:p\leq\mathsf{F}(x)\\}$. When $\mathsf{F}$
is continuous, the quantile function $\mathsf{Q}=\mathsf{F}^{-1}$.

#### Ranks and Quantiles for multivariate distributions.

Since there exists no natural ordering in $\mathbb{R}^{d}$, defining ranks and
quantiles in high dimension is not straightforward. To extend the notion of
rank in $\mathbb{R}^{d}$, theory of Optimal Transport (OT) has been used to
propose meaningful and useful notions of multivariate rank and quantile
functions (Hallin et al., 2017; Chernozhukov et al., 2017; Deb and Sen, 2019;
Hallin et al., 2020). In it’s most standard-setting, given two distributions,
a source distribution $\mu\in\mathcal{P}(\mathbb{R}^{d})$ and a target
distribution $\nu\in\mathcal{P}(\mathbb{R}^{d})$, OT aims to find a map
$\mathbf{T}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ that pushes $\mu$ to
$\nu$ with a minimal cost. That is, given $\bm{X}\sim\mu$, and
$\bm{Y}\sim\nu$, OT finds a map $\mathbf{T}$,

$\displaystyle\inf_{\mathbf{T}}\int\|\bm{x}-\mathbf{T}(\bm{x})\|^{2}d\mu(\bm{x})\;\;\text{subject
to.}\;\;\bm{Y}=\mathbf{T}(\bm{X})\sim\nu.$ (6)

Note that if $\mathbf{T}(\bm{X})\sim\nu$ when $\bm{X}\sim\mu$, we write
$\nu=\mathbf{T}_{\\#}\mu$. The key insight in using the theory of OT to
multivariate ranks and quantiles comes from noticing that in case of $d=1$,
the optimal transport map is given by
$\mathbf{T}=\mathsf{F}_{\nu}^{-1}\circ\mathsf{F}_{\mu}$, where
$\mathsf{F}_{\mu}$ and $\mathsf{F}_{\nu}$ are the distribution functions for
$\mu$ and $\nu$, respectively. When $\nu=\mathsf{U}[0,1]$, this gives the rank
function $\mathsf{F}_{\mu}$. Following (Deb and Sen, 2019), McCann’s theorem
(McCann et al., 1995) stated below, is used to extend the notion of rank to
the multivariate setting.

######  Theorem 1 ((McCann et al., 1995)).

Assume $\mu,\nu\in\mathcal{P}(\mathbb{R}^{d})$ be absolutely continuous
measures, then there exists transport maps $\mathbf{R}(\cdot)$ and
$\mathbf{Q}(\cdot)$, that are gradients of real-valued $d$-variate convex
functions such that $\mathbf{R}_{\\#}\mu=\nu,\;\;\mathbf{Q}_{\\#}\nu=\mu$,
$\mathbf{R}$ and $\mathbf{Q}$ are unique and
$\mathbf{R}\circ\mathbf{Q}(\bm{X})=\bm{X}$,
$\mathbf{Q}\circ\mathbf{R}(\bm{Y})=\bm{Y}$.

Based on this result, the authors in (Deb and Sen, 2019) give the following
definitions for the rank and quantile functions in high dimensions.

######  Definition 1 ((Deb and Sen, 2019)).

Given an absolutely continuous measure $\mu\in\mathcal{P}(\mathbb{R}^{d})$ and
$\nu=\mathsf{U}[0,1]^{d}$ \- the uniform measure on the unit cube in
$\mathbb{R}^{d}$, the ranks and quantile maps for $\mu$ are defined as the
maps $\mathbf{R}(\cdot)$ and $\mathbf{Q}(\cdot)$, respectively as defined in
Theorem 1.

###  3.1 Rank Energy

To state the definition of Rank Energy (Deb and Sen, 2019), we begin with a
brief introduction on the Energy distance (Baringhaus and Franz, 2004; Székely
and Rizzo, 2013). Energy distance is a multivariate, two-sample goodness-of-
fit measure. Given two independent multivariate random variables
$\bm{X}\in\mathbb{R}^{d}\sim\mu_{\bm{X}},\bm{Y}\in\mathbb{R}^{d}\sim\mu_{\bm{Y}}$,
where $\mu_{\bm{X}},\mu_{\bm{Y}}\in\mathcal{P}(\mathbb{R}^{d})$, energy
distance is defined via:

$\displaystyle\text{E}(\bm{X},\\!\bm{Y})\\!=\\!\gamma_{d}\\!\int_{\mathbb{R}}\\!\int_{\mathcal{S}^{d-1}}\\!\Big{(}\\!\mathbb{P}(\bm{a}^{\top}\\!\bm{X}\\!\leq\\!t)-\mathbb{P}(\bm{a}^{\top}\\!\bm{Y}\\!\leq\\!t)\Big{)}^{2}\\!\\!\\!d\kappa(\bm{a})dt,$
(7)

where $\gamma=(2\Gamma(d/2))^{-1}\sqrt{\pi}(d-1)\Gamma\big{(}(d-1)/2\big{)}$
for $d>1$,
$\mathcal{S}^{d-1}\overset{.}{=}\\{\bm{x}\in\mathbb{R}^{d}:\|\bm{x}\|=1\\}$,
and $\kappa(\cdot)$ denotes the uniform measure on $\mathcal{S}^{d-1}$. It is
self-evident that for all $\bm{a}\in\mathcal{S}^{d-1}$ and $t\in\mathbb{R}$,
the two measures will be equal if and only if
$\mathbb{P}(\bm{a}^{\top}\bm{X}\leq t)=\mathbb{P}(\bm{a}^{\top}\bm{Y}\leq t)$.

Deb (Deb and Sen, 2019) proposed a rank-based version of the energy measure
and defined it as the rank energy measure.

######  Definition 2 ((Deb and Sen, 2019)).

Suppose that $\bm{X}\sim\mu_{\bm{X}}$ and $\bm{Y}\sim\mu_{\bm{Y}}$. Let
$\mathbf{R}_{\lambda}$ denote the population rank corresponding to a mixture
distribution, $\lambda\mu_{\bm{X}}+(1-\lambda\mu_{\bm{Y}})$, for some
$\lambda\in(0,1)$. Then the population rank energy
$\mathop{}\\!\mathrm{RE}_{\lambda}$ is defined as:

$\displaystyle\mathop{}\\!\mathrm{RE}_{\lambda}(\bm{X},\bm{Y})=\gamma_{d}\int_{\mathbb{R}}\int_{\mathcal{S}^{d-1}}\Big{(}\mathbb{P}(\bm{a}^{\top}\mathbf{R}_{\lambda}(\bm{X})\leq
t)-$ $\displaystyle\mathbb{P}(\bm{a}^{\top}\mathbf{R}_{\lambda}(\bm{Y})\leq
t)\Big{)}^{2}d\kappa(\bm{a})dt.$ (8)

In other words, rank energy between $\bm{X}$ and $\bm{Y}$ is the energy
measure between the corresponding ranks, $\mathbf{R}_{\lambda}(\bm{X})$ and
$\mathbf{R}_{\lambda}(\bm{Y})$.

$\mathop{}\\!\mathrm{RE}_{\lambda}(\cdot,\cdot)$ is a distribution-free
multivariate two-sample goodness-of-fit measure. A nice feature of
$\text{RE}_{\lambda}$ is that, in 1-D, it is equivalent to the widely used
Cramér-Von Mises (Anderson, 1962) statistic for two-sample distribution
testing, (Deb and Sen, 2019).

$\displaystyle\frac{1}{2}\text{RE}_{\lambda}(X,Y)=\int(\mathsf{F}_{X}(t)-\mathsf{F}_{Y}(t))^{2}\text{d}\mathsf{H}_{\lambda}(t),$
(9)

where $\mathsf{F}_{X},\mathsf{F}_{Y}$ and $\mathsf{H}_{\lambda}$ are
absolutely continuous distribution functions corresponding to the probability
measures $\mu_{X},\mu_{Y}$ and $\lambda\mu_{X}+(1-\lambda)\mu_{Y}$,
respectively. Based on this result, we note the following lemma.

######  Lemma 1.

Let $\mu_{X}$ and $\mu_{Y}$ are supported on the interval $[a,b]$ and
$[a+s,b+s]$, respectively, and their mixture measure denoted as
$\mu_{H}=\lambda\mu_{X}+(1-\lambda)\mu_{Y}$, for any $\lambda\in(0,1)$. Also,
let $\mathsf{F}_{X},\mathsf{F}_{Y}$ and $\mathsf{H}_{\lambda}$, denote the
cumulative distribution functions for $\mu_{X},\mu_{Y}$ and $\mu_{H}$,
respectively. Then $\mathop{}\\!\mathrm{RE}_{\lambda}$ as a function of $s$,
is a constant independent of $s$ for $s\geq b-a$ and $s\leq a-b$.

###### Proof.

Please see Supplementary A.1. ∎

We also verify this phenomenon empirically for higher dimensions in Section 5.

#### Sample rank energy.

To define sample rank energy (hereafter referred to as RE), we begin with the
definition of a sample rank map. Given a set of i.i.d. samples
$\\{\bm{X}_{1},\dots,\bm{X}_{m}\\}\in\mathbb{R}^{d}$ with empirical measure
$\mu_{m}^{\bm{X}}=m^{-1}\sum_{i=1}^{m}\delta_{\bm{X}_{i}}$ and a set of Halton
sequence
$\mathcal{H}_{m}^{d}:=\\{\bm{h}_{1},\dots,\bm{h}_{m}\\}\subset[0,1]^{d}$ (Deb
and Sen, 2019) with empirical measure
$\nu_{m}^{\mathbf{H}}=m^{-1}\sum_{i=1}^{m}\delta_{\bm{h}_{i}}$ that weakly
converges to $\mathsf{U}[0,1]^{d}$, to compute a sample rank map
$\mathbf{\widehat{R}}_{m}$, one solves for the following discrete optimal
transport problem,

$\displaystyle\mathbf{\hat{P}}=\arg\min_{\mathbf{P}\in\Pi}\sum_{i,j=1}^{m}\mathbf{C}_{i,j}\mathbf{P}_{i,j},$
(10)

where
$\mathbf{C}_{i,j}=\|\bm{X}_{i}-\bm{h}_{j}\|^{2},\Pi=\\{\mathbf{P}:\mathbf{P}\bm{1}=\frac{1}{m}\bm{1},\bm{1}^{\top}\mathbf{P}=\frac{1}{m}\bm{1}^{\top}\\}$.
It is well known that the solution to this problem, under the given set-up, is
one of the scaled permutation matrices. Consequently, one obtains a map
$\mathbf{\hat{R}}_{m}(\bm{X}_{i})=\bm{h}_{\sigma(i)}$, where $\sigma(i)$ is
the non-zero index in the $i$-th row of $\mathbf{\hat{P}}$.

Now to compute RE between two sets of i.i.d. samples
$\\{\bm{X}_{1},\dots,\bm{X}_{m}\\}\sim\mu_{m}^{\bm{X}}$ and
$\\{\bm{Y}_{1},\dots,\bm{Y}_{n}\\}\sim\mu_{n}^{\bm{Y}}$, a joint sample rank
map $\mathbf{\widehat{R}}_{m,n}$ is computed via solving (10) between the
joint empirical measure
$\mu_{m,n}^{\bm{X},\bm{Y}}=(m+n)^{-1}(m\mu_{m}^{\bm{X}}+n\mu_{n}^{\bm{Y}})$
and the target empirical measure
$\nu_{m,n}^{\mathbf{H}}=(m+n)^{-1}\sum_{i=1}^{m+n}\delta_{\bm{h}_{i}}$. Given
the sample ranks corresponding to $\bm{X}_{i}$’s and $\bm{Y}_{j}$’s, RE is
defined as,

$\displaystyle\text{{RE}}\doteq$
$\displaystyle\frac{2}{mn}\sum_{i=1}^{m}\sum_{j=1}^{n}\|\widehat{\mathbf{R}}{m,n}(\bm{X}_{i})-\widehat{\mathbf{R}}_{m,n}(\bm{Y}_{j})\|-\frac{1}{m^{2}}\sum_{i,j=1}^{m}\|\widehat{\mathbf{R}}{m,n}(\bm{X}_{i})-\widehat{\mathbf{R}}{m,n}(\bm{X}_{j})\|$
$\displaystyle-\frac{1}{n^{2}}\sum_{i,j=1}^{n}\|\widehat{\mathbf{R}}{m,n}(\bm{Y}_{i})-\widehat{\mathbf{R}}{m,n}(\bm{Y}_{j})\|.$
(11)

RE is distribution-free under the null for fixed sample size. However,
discrete OT (3.1) has large sample complexity $\mathcal{O}(n^{-1/d})$ (Genevay
et al., 2019) and also computationally expensive
$\mathcal{O}(n^{3}\text{log}(n))$, given the sample size $n$. Moreover, RE is
constant in the region where the supports of the two distributions are
completely disjoint (see Supplementary).

##  4 Proposed soft rank energy

Towards developing the notion of soft rank energy we first state the
Kantorovich relaxation of the Monge problem, where instead of a map, one seeks
an optimal coupling $\bm{\pi}$ between a source distribution
$\mu\in\mathcal{P}(\mathbb{R}^{d})$ and a target distribution
$\nu\in\mathcal{P}(\mathbb{R}^{d})$. That is, given $\bm{X}\sim\mu$ and
$\bm{Y}\sim\nu$, Kantorovich relaxation (Santambrogio, 2015) solves for,

$\displaystyle\operatorname*{arg\,min}_{\bm{\pi}\in\Pi(\mu,\nu)}\int
c(\bm{x},\bm{y})\text{d}\bm{\pi}(\bm{x},\bm{y}),$ (12)

where $\Pi(\mu,\nu)$ is the set of joint probability measure over the product
set $\mathbb{R}^{d}\times\mathbb{R}^{d}$ with marginals $\mu$ and $\nu$, and
typically $c(\bm{x},\bm{y})$ is a symmetric positive cost function satisfying
$c(\bm{x},\bm{x})=0$. It is well-known that the Kantorovich solution coincides
with the Monge solution, whenever the latter exists (Santambrogio, 2015).

Adding an entropic regularization term to (12) allows one to get the transport
plan in closed form (Feydy et al., 2019) and makes the OT problem
differentiable everywhere w.r.t. the weights of the input measures (Blondel et
al., 2018). For any regularizer $\varepsilon>0$, the primal formulation of the
entropy regularized OT is given by:

$\displaystyle\operatorname*{arg\,min}_{\bm{\pi}\in\Pi(\mu,\nu)}\int
c(\bm{x},\bm{y})\text{d}\bm{\pi}(\bm{x},\bm{y})+\varepsilon\text{KL}(\bm{\pi}(\bm{x},\bm{y})|\mu(\bm{x})\otimes\nu(\bm{y})),$
(13)
$\displaystyle\text{where}\;\text{KL}(\pi(\bm{x},\bm{y})|\mu(\bm{x})\otimes\nu(\bm{y}))\doteq\int\text{ln}\Big{(}\frac{\text{d}\bm{\pi}(\bm{x},\bm{y})}{\text{d}\mu(\bm{x})\text{d}\nu(\bm{y})}\Big{)}\text{d}\pi(\bm{x},\\!\bm{y}).$

The unconstrained dual of (13) is achieved by using the Fenchel-Rockafellar’s
duality theorem (Clason et al., 2021),

$\displaystyle\operatorname*{arg\,max}_{f,g}\int
f(\bm{x})\text{d}\mu(\bm{x})+\int
g(\bm{y})\text{d}\nu(\bm{y})-\\!\\!\varepsilon\\!\\!\int\\!\\!\int\\!\text{exp}\big{(}\frac{1}{\varepsilon}(f(\bm{x})\\!\\!+\\!\\!g(\bm{y})\\!\\!-\\!\\!c(\bm{x},\bm{y}))\big{)}\text{d}\mu(\bm{x})\text{d}\nu(\bm{y})\\!+\\!\varepsilon,$
(14)

where $(f,g)$ is a pair of continuous functions. Given the optimal dual
solutions $f,g$, the optimal plan $\bm{\pi}^{\varepsilon}$ is given as (Seguy
et al., 2017),

$\displaystyle\text{d}\bm{\pi}^{\varepsilon}(\bm{x},\bm{y})=\text{exp}\Big{(}\\!\frac{1}{\epsilon}(f(\bm{x})+g(\bm{y})-c(\bm{x},\bm{y}))\\!\Big{)}\text{d}\mu(\bm{x})\text{d}\nu(\bm{y}).$

Unlike $\mathbf{T}$, $\bm{\pi}^{\varepsilon}$ is not a map, but a diffused
coupling, where the degree of diffusion is directly proportional to
$\varepsilon$. Given $\bm{\pi}^{\varepsilon}$, we now define the soft rank.

######  Definition 3.

Given an absolutely continuous measure $\mu\in\mathcal{P}(\mathbb{R}^{d})$ and
$\nu=\mathsf{U}[0,1]^{d}$ \- the uniform measure on the unit cube in
$\mathbb{R}^{d}$, the soft rank is defined as:

$\displaystyle\mathbf{R}^{\varepsilon}(\bm{x})\overset{.}{=}\int_{\bm{y}}\frac{\bm{y}\,\text{d}\bm{\pi}^{\varepsilon}(\bm{x},\bm{y})}{\int_{\bm{y}}\text{d}\bm{\pi}^{\varepsilon}(\bm{x},\bm{y})}.$
(15)

In other words, soft rank $\mathbf{R}^{\varepsilon}(\cdot)$ is the conditional
expectation of $\bm{y}$ under $\bm{\pi}^{\varepsilon}$, given $\bm{X}=\bm{x}$.

We note that the soft rank as defined above corresponds to what is referred to
as the barycentric projection of optimal transport plan (Seguy et al., 2017;
Deb et al., 2021). Based on Definition 3, we now define the population version
of soft rank energy.

######  Definition 4.

Given two independent multivariate random variables $\bm{X}\sim\mu_{\bm{X}}$
and $\bm{Y}\sim\mu_{\bm{Y}}$, where
$\mu_{\bm{X}},\mu_{\bm{Y}}\in\mathcal{P}(\mathbb{R}^{d})$. Also, let
$\mathbf{R}^{\varepsilon}_{\lambda}$ denote the population soft rank
corresponding to a mixture distribution
$\lambda\mu_{\bm{X}}+(1-\lambda\mu_{\bm{Y}})$, for some $\lambda\in(0,1)$.
Then the population soft rank energy is defined as:

$\displaystyle\mathop{}\\!\mathrm{sRE}^{\varepsilon}_{\lambda}(\bm{X},\bm{Y})=\gamma_{d}$
$\displaystyle\int_{\mathbb{R}}\int_{\mathcal{S}^{d-1}}\Big{(}\mathbb{P}(\bm{a}^{\top}\mathbf{R}_{\lambda}^{\varepsilon}(\bm{X})\leq
t)-\mathbb{P}(\bm{a}^{\top}\mathbf{R}_{\lambda}^{\varepsilon}(\bm{Y})\leq
t)\Big{)}^{2}d\kappa(\mathbf{a})dt.$ (16)

In other words, soft rank energy is the energy measure between the soft ranks
corresponding to $\bm{X}$ and $\bm{Y}$.

Based on this definition, we state the following lemma.

######  Lemma 2.

Under the assumptions of Definition 4, we state the following properties of
$\mathop{}\\!\mathrm{sRE}^{\varepsilon}_{\lambda}$.

  * (a)

Given $\bm{X}_{1},\bm{X}_{2}$ i.i.d. $\sim\mu_{\bm{X}}$ and given
$\bm{Y}_{1},\bm{Y}_{2}$ i.i.d. $\sim\mu_{\bm{Y}}$, following (Székely and
Rizzo, 2013), the soft rank energy as defined via equation (16) is equivalent
to

$\displaystyle\mathop{}\\!\mathrm{sRE}_{\lambda}^{\varepsilon}(\bm{X}_{1},$
$\displaystyle\bm{Y}_{1})\\!\\!=\\!\\!2\mathbb{E}\|\mathbf{R}^{\varepsilon}_{\lambda}(\bm{X}_{1})-\mathbf{R}^{\varepsilon}_{\lambda}(\bm{Y}_{1})\|\\!\\!-\\!\\!\mathbb{E}\|\mathbf{R}^{\varepsilon}_{\lambda}(\bm{X}_{1})-\mathbf{R}^{\varepsilon}_{\lambda}(\bm{X}_{2})\|\\!\\!-\\!\\!\mathbb{E}\|\mathbf{R}^{\varepsilon}_{\lambda}(\bm{Y}_{1})\\!\\!-\\!\\!\mathbf{R}^{\varepsilon}_{\lambda}(\bm{Y}_{2})\|.$
(17)

  * (b)

(Symmetric)
$\mathop{}\\!\mathrm{sRE}^{\varepsilon}_{\lambda}(\bm{X}_{1},\bm{Y}_{1})=\mathop{}\\!\mathrm{sRE}^{\varepsilon}_{\lambda}(\bm{Y}_{1},\bm{X}_{1})$,

  * (c)

$\mathop{}\\!\mathrm{sRE}^{\varepsilon}_{\lambda}(\bm{X}_{1},\bm{Y}_{1})=0$ if
$\bm{X}_{1}\overset{d}{=}\bm{Y}_{1}$,

  * (d)

$\mathop{}\\!\mathrm{sRE}^{\varepsilon}_{\lambda}$ will converge to
$\mathop{}\\!\mathrm{RE}_{\lambda}$ as $\varepsilon\rightarrow 0$.

###### Proof.

Please see Supplementary A.2. ∎

Note that one has the flexibility to employ any multivariate two-sample test
to measure the goodness-of-fit between distributions. That is, in similar
fashion, one can define $\text{RMMD}_{\lambda}$ or
$\text{sRMMD}^{\varepsilon}_{\lambda}$ by measuring the maximum mean
discrepancy (MMD) (Gretton et al., 2012) between the rank-transformed random
variables $\mathbf{R}_{\lambda}(\bm{X})$ and $\mathbf{R}_{\lambda}(\bm{Y})$ or
$\mathbf{R}_{\lambda}^{\varepsilon}(\bm{X})$ and
$\mathbf{R}_{\lambda}^{\varepsilon}(\bm{Y})$.

We now state the sample soft rank energy $(\text{sRE}^{\varepsilon})$ and
sample soft rank maximum mean discrepancy $(\text{sRMMD}^{\varepsilon})$. For
the sake of brevity, hereafter, we refer to them as sRE and sRMMD,
respectively.

Figure 1: $\text{RMMD}(\varepsilon=0)$ and sRMMD (y-axis) between two uniform
distributions with bounded supports. $n$ and $d$ denote the number of samples
and dimension. RMMD achieves a constant maximum in the region, where the two
supports are completely disjoint, for all $n,d$. sRMMD do not saturate for
increased $d$ if $\varepsilon$ and $n$ is larger.

###  4.1 sRE and sRMMD

Using the similar setting as described in Section 3.1, we now define the
sample soft rank map $\mathbf{\widehat{R}}_{m}^{\varepsilon}$. To compute
$\mathbf{\widehat{R}}_{m}^{\varepsilon}$, first an entropy regularized OT is
solved via Sinkhorn (Peyré et al., 2019)-a iterative fixed point algorithm for
an entropy regularizer, $\varepsilon>0$,

$\displaystyle\mathbf{P}^{\varepsilon}=\arg\min_{\mathbf{P}\in\Pi}\sum_{i,j=1}^{m}\mathbf{C}_{i,j}\mathbf{P}_{i,j}-\varepsilon
H(\mathbf{P}),$ (18)

where
$\mathbf{C}_{i,j}=\|\bm{X}_{i}-\bm{h}_{j}\|^{2},\Pi=\\{\mathbf{P}:\mathbf{P}\bm{1}=\frac{1}{m}\bm{1},\bm{1}^{\top}\mathbf{P}=\frac{1}{m}\bm{1}^{\top}\\}$,
and $H(\mathbf{P})=-\sum_{i,j}\mathbf{P}_{i,j}\log\mathbf{P}_{i,j}$ is the
entropy functional. Given $\mathbf{P}^{\varepsilon}$, sample soft rank map is
defined via,

$\displaystyle\widehat{\mathbf{R}}_{m}^{\varepsilon}(\bm{X}_{i})=\sum_{j=1}^{m}\frac{\mathbf{P}_{i,j}^{\varepsilon}}{\sum_{j=1}^{m}\mathbf{P}_{i,j}^{\varepsilon}}\bm{h}_{j}.$
(19)

Now, to define sRE and sRMMD, we first compute a joint sample soft rank map
$\hat{\mathbf{R}}_{m,n}^{\varepsilon}$ via solving (18) between
$\mu_{m,n}^{\bm{X},\bm{Y}}$ and $\nu_{m,n}^{\mathbf{H}}$. Given the sample
soft ranks corresponding to $\bm{X}_{i}$’s and $\bm{Y}_{j}$’s, we define sRE
as,

$\displaystyle\text{{sRE}}\doteq$
$\displaystyle\frac{2}{mn}\sum_{i=1}^{m}\sum_{j=1}^{n}\|\widehat{\mathbf{R}}^{\varepsilon}_{m,n}(\bm{X}_{i})-\widehat{\mathbf{R}}_{m,n}^{\varepsilon}(\bm{Y}_{j})\|-\frac{1}{m^{2}}\sum_{i,j=1}^{m}\|\widehat{\mathbf{R}}^{\varepsilon}_{m,n}(\bm{X}_{i})-\widehat{\mathbf{R}}^{\varepsilon}_{m,n}(\bm{X}_{j})\|$
$\displaystyle-\frac{1}{n^{2}}\sum_{i,j=1}^{n}\|\widehat{\mathbf{R}}^{\varepsilon}_{m,n}(\bm{Y}_{i})-\widehat{\mathbf{R}}^{\varepsilon}_{m,n}(\bm{Y}_{j})\|$
(20)

and sRMMD, with a characteristic kernel function $k(\cdot,\cdot)$, as,

$\displaystyle\text{{sRMMD}}\doteq$
$\displaystyle\frac{1}{m(m-1)}\sum_{i,j\neq
i}^{m}k(\widehat{\mathbf{R}}_{m,n}^{\epsilon}(\bm{X}_{i}),\widehat{\mathbf{R}}_{m,n}^{\epsilon}(\bm{X}_{j}))-\frac{2}{mn}\sum_{i,j=1}^{m,n}k(\widehat{\mathbf{R}}_{m,n}^{\epsilon}(\bm{X}_{i}),\widehat{\mathbf{R}}_{m,n}^{\epsilon}(\bm{Y}_{j}))$
$\displaystyle+\frac{1}{n(n-1)}\\!\\!\sum_{i,j\neq
i}^{n}\\!\\!k(\widehat{\mathbf{R}}_{m,n}^{\epsilon}(\bm{Y}_{i}),\widehat{\mathbf{R}}_{m,n}^{\epsilon}(\bm{Y}_{j})),$
(21)

Employing entropic OT to define rank maps make sRE and sRMMD differentiable
everywhere (Seguy et al., 2017). Use of Sinkhorn iterations (Peyré et al.,
2019) allows one to exploit the power of GPU. The computational complexity of
entropic OT is
$\mathcal{O}(\varepsilon^{-2}n^{2}\log\,n\|\mathbf{C}\|_{\infty}^{2})$ (Peyré
et al., 2019) and with increasing $\varepsilon$, complexity decreases.
Moreover, under some mild assumptions, namely sub-Gaussianity of the measures,
the estimation of entropic OT does not suffer from the curse of dimensionality
for a sufficiently large $\varepsilon$ (Mena and Niles-Weed, 2019).

###  4.2 Proposed Knockoff generation

We now propose to use sRMMD as the loss function in a generative model to
produce knockoffs. Given the design matrix $\mathbf{X}\in\mathbb{R}^{n\times
d}$ and its knockoff matrix $\mathbf{\tilde{X}}\in\mathbb{R}^{n\times d}$, the
loss function is defined via,

$\displaystyle\ell(\mathbf{X},\mathbf{\tilde{X}})\\!\\!=\\!\\!\ell_{\text{sRMMD}}(\mathbf{X},\tilde{\mathbf{X}})\\!\\!+\\!\\!\lambda\ell_{\text{s-o}}(\mathbf{X},\mathbf{\tilde{X}})\\!\\!+\\!\\!\delta\text{D}_{\text{corr}}(\mathbf{X},\tilde{\mathbf{X}}).$
(22)

$\ell_{\text{sRMMD}}(\mathbf{X},\tilde{\mathbf{X}})$ is computed in a similar
way as (3),

$\displaystyle\ell_{\text{sRMMD}}(\mathbf{X},\tilde{\mathbf{X}})$
$\displaystyle=\text{{sRMMD}}\Big{[}(\mathbf{X}^{\prime},\tilde{\mathbf{X}^{\prime}}),(\tilde{\mathbf{X}^{\prime\prime}},\mathbf{X}^{\prime\prime})\Big{]}+\text{{sRMMD}}\Big{[}(\mathbf{X}^{\prime},\tilde{\mathbf{X}^{\prime}}),(\mathbf{X}^{\prime\prime},\tilde{\mathbf{X}^{\prime\prime}})_{\text{swap}(S)}\Big{]}.$

For optimal performance, the hyperparameters should be tuned to the specific
data distribution at hand.

##  5 Experiments

In this section, we empirically observe the properties of sRMMD as a loss
function and evaluate the performance of the proposed knockoff generation
method on synthetic and real datasets.

####  sRMMD as a loss function.

To investigate the quality of sRMMD as a loss function in a generative model,
we consider two Uniform distributions having bounded supports,
$\mathsf{U}[0,1]^{d}$ and $\mathsf{U}[s,s+1]^{d}$ and measure the sRMMD
between them for $s\in[-10,10]$. Figure 1 empirically shows the saturation
behavior of sRMMD.

This saturation limits the use of sRMMD for learning a generative model, as
the saturation in those regimes will imply vanishing gradients for learning.
As one can see from the plots that either higher values of $\varepsilon$ or
larger sample size is needed with increasing dimensions in order to prevent
the saturation behavior of sRMMD. Plots showing similar behavior for sRE can
be found in Supplementary C.2.

Figure 2: Mode reconstruction for a Gaussian mixture model (upper left) with
sRMMD for various $\varepsilon$

.

(a)

(b)

(c)

(d)

Figure 3: FDR vs amplitude and power vs. amplitude. FDR level is set to 0.1
(black dotted line).

####  sRMMD for a generative model.

To check this, we sample
$\bm{X}\sim\sum_{k=1}^{3}m_{k}\mathcal{N}(\mu_{k},\Sigma_{k})$, which is a
Gaussian mixture model with three modes. $\Sigma_{k}$ is a $d$-dimensional
covariance matrix whose $(i,j)$ entry is $\rho^{|i-k|}$. We use
$(\rho_{1},\rho_{2},\rho_{2})=(0.6,0.4,0.2)$, cluster mean
$(\mu_{1},\mu_{2},\mu_{3})=(0,20,40)$ and mixture proportion
$(m_{1},m_{2},m_{3})=(0.4,0.2,0.4)$ for $d=16$. We use (22) with
$(\lambda,\gamma)=(0,0)$ in the generative model described in (Romano et al.,
2020).

Figure 2 shows that the generative model suffers from mode collapse for
$\varepsilon=1$, which is expected (Figure 1). On the other hand,
$\varepsilon=5,10$, perfectly captures the mode, covariance, and proportions
of the Gaussian mixture model, hence establish sRMMD as a suitable loss
function for the knockoff generation.

#### Synthetic dataset.

We evaluate the performance of sRMMD based knockoffs on four different
distributional settings adapted from (Romano et al., 2020). Below, we briefly
describe each setting with the optimal hyperparameters used in (22) to train
the model.

  1. (a)

_Multivariate Gaussian_ : An AR(1) model, with
$\bm{X}=(X_{1},X_{2},\dots,X_{d})\sim\mathcal{N}(0,\Sigma)$, where $\Sigma$\-
a $d$ dimensional covariance matrix, whose $(i,j)$ entry is $\rho^{|i-j|}$.
$\rho$ is set to $0.5$. Hyperparameters $(\lambda,\delta)$ were set to $(1,1)$
in the loss function (22).

  2. (b)

_Gaussian Mixture Model_ : Mixture of three multivariate Gaussian
distributions, $\mathcal{N}(0,\Sigma_{1})$, $\mathcal{N}(0,\Sigma_{2})$ and
$\mathcal{N}(0,\Sigma_{1})$, each of equal probability.
$\Sigma_{1},\Sigma_{2},\Sigma_{3}$ have the same first-order autoregressive
structure. Covariance matrices are
$\Sigma_{1}=0.3^{|i-j|},\Sigma_{2}=0.5^{|i-j|}$, and $\Sigma_{3}=0.7^{|i-j|}$.
$(\lambda,\delta)=(1,1)$.

  3. (c)

_Multivariate Student’s t distribution_ :
$\bm{X}=\sqrt{\frac{(\nu-2)}{\nu}}\frac{\bm{Z}}{\sqrt{\Gamma}}$, where $\nu=3$
is the degrees of freedom, $\bm{Z}\sim\mathcal{N}(0,\Sigma)$ and $\Gamma$ is
independently drawn from a Gamma distribution with shape and rate parameters
both equal to $\nu/2$. $(\lambda,\delta)=(1,0.1)$.

  4. (d)

_Sparse Gaussian Variables_ : For a scalar random variable
$\eta\sim\mathcal{N}(0,1)$, and random subset $A\in\\{1,\dots,p\\}$ of size
$|A|=L$, sparse Gaussian variable is defined as,

$\displaystyle
X_{j}=\sqrt{\frac{\binom{L}{p}}{\binom{L-1}{p-1}}}.\begin{cases}\eta,\;\;\text{if}\;\;j\in
A,\\\ 0,\;\;\text{otherwise}\end{cases}\;\;$
$\displaystyle\text{and,}\;\;\Sigma_{i,j}=\begin{cases}1,\;\;\text{if}\;\;i=j,\\\
\frac{L-1}{p-1}\;\;\text{otherwise.}\end{cases}$

with $L=30$, and $(\lambda,\delta)=(1,5)$.

Experimental setup: We consider a fully connected neural network with 6 hidden
layers, each consists of 5$d$ nodes as the generative architecture. This
architecture is almost similar to the structure used in (Romano et al., 2020).
The network takes $\bm{X}\in\mathbb{R}^{d}$ and a noise vector
$V\in\mathbb{R}^{d}\sim\mathcal{N}(\bm{0},\mathbb{I}_{d})$ as inputs and
produce the knockoff vector $\tilde{\bm{X}}\in\mathbb{R}^{d}$. For further
insights, we refer the reader to (Romano et al., 2020).

We train the model with Stochastic gradient descent with $n=2500$ training
samples, each of $d=100$ dimension. The minibatch size, learning rate and the
number of epochs were set to 500, $0.01$, and 100, respectively. We use
$\varepsilon=10$, to compute the soft ranks. We employ a Gaussian mixture
kernel
$k(x,y)=\frac{1}{8}\sum_{i=1}^{8}\text{exp}[-\|x-y\|]_{2}^{2}/(2\sigma^{2})$,
with $\sigma=(1,2,4,8,16,32,64,128)$ to compute sRMMD.

We design a matrix $\mathbf{X}\in\mathbb{R}^{m\times d}$, with $m=100$ random
samples unknown to the training dataset and generate knockoffs
$\mathbf{\tilde{X}}\in\mathbb{R}^{m\times d}$. For each sample
$i\in\\{1,2,\dots,m\\}$, we simulate the response according to a Gaussian
linear model, $y_{i}=\bm{X}_{i}^{\top}\bm{\beta}+z$, where
$z\sim\mathcal{N}(0,1)$. $\bm{\beta}\in\mathbb{R}^{d}$ denotes the coefficient
vector which has 30 non-zero entries each having an amplitude equal to
$a/\sqrt{m}$. To find the knockoff statistics, we fit a LASSO-regression
(Friedman et al., 2010) model on the augmented matrix
$[\mathbf{X},\mathbf{\tilde{X}}]\in\mathbb{R}^{m\times 2d}$ with the response
vector $\bm{y}\in\mathbb{R}^{m}$ and compute the coefficients
$[\hat{\bm{\beta}},\hat{\bm{\beta}}^{K}]\in\mathbb{R}^{2d}$ via solving,

$\displaystyle(\hat{\bm{\beta}},\hat{\bm{\beta}}^{K})=\arg\min_{(\bm{\beta},\bm{\beta}^{K})}$
$\displaystyle\frac{1}{2m}\|y-\mathbf{X}\bm{\beta}-\mathbf{\tilde{X}}\bm{\beta}^{K}\|_{2}^{2}+\alpha(\|\bm{\beta}\|_{1}+\|\bm{\beta}^{K}\|_{1}),$
(23)

where $\hat{\bm{\beta}}$ and $\hat{\bm{\beta}}^{K}$ denotes the coefficient
vectors corresponding to the original variables and the knockoffs,
respectively, and $\alpha$ is a hyperparameter chosen carefully for each
distributional setting. Absolute LASSO coefficient difference is taken as the
knockoff statistics,
$W_{j}=|\hat{\bm{\beta}}_{j}|-|\hat{\bm{\beta}}_{j}^{k}|$, for each
$j\in\\{1,\dots,d\\}$. We apply (2) on the statistics $W_{j}$ for the FDR
level $q=0.1$. We repeat the experiment 500 times for different values of $a$
and take the average FDR and power.

We compare our method to other knockoff generation techniques, namely Second-
order knockoff (Candes et al., 2016), Deep Knockoffs (Romano et al., 2020),
KnockoffGAN (Jordon et al., 2018) and DDLK (Sudarshan et al., 2020). For all
comparisons, we use publicly available implementations of the code (if
available) and used their recommended configurations and hyperparameters. For
each method, we follow the exact procedure to compute knockoff statistics as
described above.

#### Performance on synthetic data.

Figure 3 shows the average FDR vs. amplitude and power vs. amplitude curves
for different knockoff generation methods. For the multivariate Gaussian
distribution, we observe that, all methods showcase very similar power and FDR
control over the entire amplitude region. For the Gaussian mixture model,
KnockoffGAN shows higher power, but cannot control the FDR, whereas sRMMD and
other methods demonstrate comparable power and keep the FDR below the FDR
level. We notice that for the multivariate Student’s t distribution, Second-
order and DDLK knockoffs achieves comparably higher power in the low amplitude
region but fail to control the FDR. We also observe that, though KnockoffGAN
and Deep knockoffs can control the FDR but they have low power in the small
amplitude region. On the other hand, sRMMD gains higher power than the
KnockoffGAN and Deep knockoff at this region and is also able to control the
FDR. For the sparse Gaussian setting, we notice that DDLK gains very high
power in the small amplitude region, but has poor control over FDR. Second-
order knockoff also fails to control the FDR in this case. On the other hand,
the plot indicates that knockoffGAN and Deep knockoff control the FDR but have
less power in the small amplitude region, whereas sRMMD gains higher power and
more control over the FDR in this region.

#### Evaluation on real data benchmark

We apply the proposed knockoff filter to a publicly available metabolomics
dataset in order to discover important biomarkers. In absence of the ground
truth, to qualitatively analyze the performance, we cross-reference the
selected metabolites with published literature. We use a study titled
_Longitudinal Metabolomics of the Human Microbiome in Inflammatory Bowel
Disease_ (Lloyd-Price et al., 2019) which is available at the NIH’s common
Fund’s National Metabolomics Data Repository (NMDR) website, the Metabolomics
Workbench, https://www.metabolomicsworkbench.org/ under the project DOI:
10.21228/M82T15. The study is related to Inflammatory Bowel Disease (IBD),
that describes the conditions like ulcerative colitis (UC) and Crohrn’s
disease (CD). We use _C18 Reverse-Phase negative mode_ dataset which was
collected under this study. There are $546$ samples in this dataset, each
having $91$ metabolites on average. Each sample belongs to one of the three
different classes e.g., CD, UC, non-IBD and we assign response $y$ to either
$\\{0,1,\text{or}\;2\\}$, respectively. Before applying the knockoff filter,
we preprocess the dataset in multiple steps, (i) keeping the metabolites that
have at least $80\%$ filled values (ii) missing value imputation using
K-nearest neighbour (KNN) (Gromski et al., 2014) (iii) standardization. After
generating the knockoffs, we apply Random Forest classifier (Trainor et al.,
2017) on the augmented matrix, and take the difference between the out-of-bag
(OOB) (Trainor et al., 2017) scores corresponding to the original and
knockoffs as the knockoff statistics. We repeat the whole procedure $100$
times, since the generated knockoffs are random. We select those metabolites
that appear at least 70 times out of $100$ instances, setting the FDR level at
0.05. 22 metabolites are found as significant. Among them 19 metabolites are
shown to have relation to IBD, which are cross referenced by published
literature.

Additional experiments and the list of the selected metabolites are to be
found in Supplementary C, and, B.

##  6 Conclusion and Future work

In this paper, we introduced a new statistic called soft-rank Maximum Mean
Discrepancy (sRMMD) and used it to generate valid knockoffs. We demonstrate
through a series of experiments that sRMMD is a valid loss function that can
be used in a generative model. We also showed that knockoffs generated using
the proposed method have better power and FDR control on synthetic datasets
sampled from Gaussian and non-Gaussian distributional settings. While the
proposed approach is better than several existing methods for knockoff
generation, it is still computationally expensive compared to MMD. An
important future work is to characterize in a precise manner the behaviour of
sRMMD w.r.t. the sample size, dimension, and the entropy regularizer.

##  7 Acknowledgement

This research was sponsored by the U.S. Army DEVCOM Soldier Center, and was
accomplished under Cooperative Agreement Number W911QY-19-2-0003. The views
and conclusions contained in this document are those of the authors and should
not be interpreted as representing the official policies, either expressed or
implied, of the U.S. Army DEVCOM Soldier Center, or the U.S. Government. The
U. S. Government is authorized to reproduce and distribute reprints for
Government purposes notwithstanding any copyright notation hereon.

We also acknowledge support from the U.S. National Science Foundation under
award HDR-1934553 for the Tufts T-TRIPODS Institute. Shuchin Aeron is also
supported in part by NSF CCF:1553075, NSF RAISE 1931978, NSF ERC planning
1937057, and AFOSR FA9550-18-1-0465.

## References

  * Anderson [1962] T. W. Anderson.  On the distribution of the two-sample cramer-von mises criterion.  _The Annals of Mathematical Statistics_ , pages 1148–1159, 1962.
  * Baringhaus and Franz [2004] L. Baringhaus and C. Franz.  On a new multivariate two-sample test.  _Journal of multivariate analysis_ , 88(1):190–206, 2004.
  * Bauset et al. [2021] C. Bauset, L. Gisbert-Ferrándiz, and J. Cosín-Roger.  Metabolomics as a promising resource identifying potential biomarkers for inflammatory bowel disease.  _Journal of Clinical Medicine_ , 10(4):622, 2021\.
  * Benjamini and Hochberg [1995] Y. Benjamini and Y. Hochberg.  Controlling the false discovery rate: a practical and powerful approach to multiple testing.  _Journal of the Royal statistical society: series B (Methodological)_ , 57(1):289–300, 1995.
  * Bińkowski et al. [2018] M. Bińkowski, D. J. Sutherland, M. Arbel, and A. Gretton.  Demystifying mmd gans.  _arXiv preprint arXiv:1801.01401_ , 2018.
  * Blondel et al. [2018] M. Blondel, V. Seguy, and A. Rolet.  Smooth and sparse optimal transport.  In _International Conference on Artificial Intelligence and Statistics_ , pages 880–889. PMLR, 2018.
  * Candes et al. [2016] E. Candes, Y. Fan, L. Janson, and J. Lv.  Panning for gold: Model-x knockoffs for high-dimensional controlled variable selection.  _arXiv preprint arXiv:1610.02351_ , 2016.
  * Caprilli et al. [2009] R. Caprilli, M. Cesarini, E. Angelucci, and G. Frieri.  The long journey of salicylates in ulcerative colitis: The past and the future.  _Journal of Crohn’s and Colitis_ , 3(3):149–156, 2009.
  * Carlier et al. [2017] G. Carlier, V. Duval, G. Peyré, and B. Schmitzer.  Convergence of entropic schemes for optimal transport and gradient flows.  _SIAM Journal on Mathematical Analysis_ , 49(2):1385–1418, 2017.
  * Chernozhukov et al. [2017] V. Chernozhukov, A. Galichon, M. Hallin, M. Henry, et al.  Monge–kantorovich depth, quantiles, ranks and signs.  _Annals of Statistics_ , 45(1):223–256, 2017\.
  * Clason et al. [2021] C. Clason, D. A. Lorenz, H. Mahler, and B. Wirth.  Entropic regularization of continuous optimal transport problems.  _Journal of Mathematical Analysis and Applications_ , 494(1):124432, 2021.
  * Cuturi et al. [2019] M. Cuturi, O. Teboul, and J.-P. Vert.  Differentiable ranks and sorting using optimal transport.  _arXiv preprint arXiv:1905.11885_ , 2019.
  * Deb and Sen [2019] N. Deb and B. Sen.  Multivariate rank-based distribution-free nonparametric testing using measure transportation.  _arXiv preprint arXiv:1909.08733_ , 2019.
  * Deb et al. [2021] N. Deb, P. Ghosal, and B. Sen.  Rates of estimation of optimal transport maps using plug-in estimators via barycentric projections.  _arXiv preprint arXiv:2107.01718_ , 2021.
  * Feydy et al. [2019] J. Feydy, T. Séjourné, F.-X. Vialard, S.-i. Amari, A. Trouvé, and G. Peyré.  Interpolating between optimal transport and mmd using sinkhorn divergences.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pages 2681–2690. PMLR, 2019.
  * Fretland et al. [1990] D. Fretland, D. Widomski, S. Levin, and T. Gaginella.  Colonic inflammation in the rabbit induced by phorbol-12-myristate-13-acetate.  _Inflammation_ , 14(2):143–150, 1990.
  * Friedman et al. [2010] J. Friedman, T. Hastie, and R. Tibshirani.  Regularization paths for generalized linear models via coordinate descent.  _Journal of statistical software_ , 33(1):1, 2010\.
  * Frommlet et al. [2012] F. Frommlet, F. Ruhaltinger, P. Twaróg, and M. Bogdan.  Modified versions of bayesian information criterion for genome-wide association studies.  _Computational Statistics & Data Analysis_, 56(5):1038–1051, 2012.
  * Gavrilov et al. [2009] Y. Gavrilov, Y. Benjamini, and S. K. Sarkar.  An adaptive step-down procedure with proven fdr control under independence.  _The Annals of Statistics_ , 37(2):619–629, 2009\.
  * Genevay et al. [2019] A. Genevay, L. Chizat, F. Bach, M. Cuturi, and G. Peyré.  Sample complexity of sinkhorn divergences.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pages 1574–1583. PMLR, 2019.
  * Gimenez et al. [2019] J. R. Gimenez, A. Ghorbani, and J. Zou.  Knockoffs for the mass: new feature importance statistics with false discovery guarantees.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pages 2125–2133. PMLR, 2019.
  * Gretton et al. [2012] A. Gretton, K. M. Borgwardt, M. J. Rasch, B. Schölkopf, and A. Smola.  A kernel two-sample test.  _The Journal of Machine Learning Research_ , 13(1):723–773, 2012.
  * Gromski et al. [2014] P. S. Gromski, Y. Xu, H. L. Kotze, E. Correa, D. I. Ellis, E. G. Armitage, M. L. Turner, and R. Goodacre.  Influence of missing values substitutes on multivariate analysis of metabolomics data.  _Metabolites_ , 4(2):433–452, 2014.
  * Hallin et al. [2020] M. Hallin, D. Hlubinka, and Š. Hudecová.  Fully distribution-free center-outward rank tests for multiple-output regression and manova.  _arXiv preprint arXiv:2007.15496_ , 2020.
  * Hallin et al. [2021] M. Hallin, E. del Barrio, J. Cuesta-Albertos, and C. Matrán.  Distribution and quantile functions, ranks and signs in dimension d: A measure transportation approach.  _The Annals of Statistics_ , 49(2):1139–1165, 2021.
  * Hallin et al. [2017] M. Hallin et al.  On distribution and quantile functions, ranks and signs in rd.  _ECARES Working Papers_ , 2017.
  * Heinemann [2019] J. Heinemann.  Machine learning in untargeted metabolomics experiments.  In _Microbial Metabolomics_ , pages 287–299. Springer, 2019.
  * Jordon et al. [2018] J. Jordon, J. Yoon, and M. van der Schaar.  Knockoffgan: Generating knockoffs for feature selection using generative adversarial networks.  In _International Conference on Learning Representations_ , 2018.
  * Jović et al. [2015] A. Jović, K. Brkić, and N. Bogunović.  A review of feature selection methods with applications.  In _2015 38th international convention on information and communication technology, electronics and microelectronics (MIPRO)_ , pages 1200–1205. Ieee, 2015.
  * Kingma and Ba [2014] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kingma and Welling [2013] D. P. Kingma and M. Welling.  Auto-encoding variational bayes.  _arXiv preprint arXiv:1312.6114_ , 2013.
  * Kingma and Welling [2019] D. P. Kingma and M. Welling.  An introduction to variational autoencoders.  _arXiv preprint arXiv:1906.02691_ , 2019.
  * Kuroki et al. [1997] F. Kuroki, M. Iida, T. Matsumoto, K. Aoyagi, K. Kanamoto, and M. Fujishima.  Serum n3 polyunsaturated fatty acids are depleted in crohn’s disease.  _Digestive diseases and sciences_ , 42(6):1137–1141, 1997.
  * Lavelle and Sokol [2020] A. Lavelle and H. Sokol.  Gut microbiota-derived metabolites as key actors in inflammatory bowel disease.  _Nature Reviews Gastroenterology & Hepatology_, 17(4):223–237, 2020.
  * Le Gall et al. [2011] G. Le Gall, S. O. Noor, K. Ridgway, L. Scovell, C. Jamieson, I. T. Johnson, I. J. Colquhoun, E. K. Kemsley, and A. Narbad.  Metabolomics of fecal extracts detects altered metabolic activity of gut microbiota in ulcerative colitis and irritable bowel syndrome.  _Journal of proteome research_ , 10(9):4208–4218, 2011.
  * Lee et al. [2017] T. Lee, T. Clavel, K. Smirnov, A. Schmidt, I. Lagkouvardos, A. Walker, M. Lucio, B. Michalke, P. Schmitt-Kopplin, R. Fedorak, et al.  Oral versus intravenous iron replacement therapy distinctly alters the gut microbiota and metabolome in patients with ibd.  _Gut_ , 66(5):863–871, 2017.
  * Levan et al. [2019] S. R. Levan, K. A. Stamnes, D. L. Lin, A. R. Panzer, E. Fukui, K. McCauley, K. E. Fujimura, M. McKean, D. R. Ownby, E. M. Zoratti, et al.  Elevated faecal 12, 13-dihome concentration in neonates at high risk for asthma is produced by gut bacteria and impedes immune tolerance.  _Nature microbiology_ , 4(11):1851–1861, 2019\.
  * Li et al. [2015] Y. Li, K. Swersky, and R. Zemel.  Generative moment matching networks.  In _International Conference on Machine Learning_ , pages 1718–1727. PMLR, 2015.
  * Liu and Rigollet [2019] J. Liu and P. Rigollet.  Power analysis of knockoff filters for correlated designs.  _arXiv preprint arXiv:1910.12428_ , 2019.
  * Liu and Zheng [2018] Y. Liu and C. Zheng.  Auto-encoding knockoff generator for fdr controlled variable selection.  _arXiv preprint arXiv:1809.10765_ , 2018.
  * Lloyd-Price et al. [2019] J. Lloyd-Price, C. Arze, A. N. Ananthakrishnan, M. Schirmer, J. Avila-Pacheco, T. W. Poon, E. Andrews, N. J. Ajami, K. S. Bonham, C. J. Brislawn, et al.  Multi-omics of the gut microbial ecosystem in inflammatory bowel diseases.  _Nature_ , 569(7758):655–662, 2019.
  * McCann et al. [1995] R. J. McCann et al.  Existence and uniqueness of monotone measure-preserving maps.  _Duke Mathematical Journal_ , 80(2):309–324, 1995\.
  * Mena and Niles-Weed [2019] G. Mena and J. Niles-Weed.  Statistical bounds for entropic optimal transport: Sample complexity and the central limit theorem.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Peyré et al. [2019] G. Peyré, M. Cuturi, et al.  Computational optimal transport: With applications to data science.  _Foundations and Trends® in Machine Learning_ , 11(5-6):355–607, 2019.
  * Price et al. [2017] N. D. Price, A. T. Magis, J. C. Earls, G. Glusman, R. Levy, C. Lausted, D. T. McDonald, U. Kusebauch, C. L. Moss, Y. Zhou, et al.  A wellness study of 108 individuals using personal, dense, dynamic data clouds.  _Nature biotechnology_ , 35(8):747, 2017.
  * Qin [2012] X. Qin.  Etiology of inflammatory bowel disease: a unified hypothesis.  _World journal of gastroenterology: WJG_ , 18(15):1708, 2012.
  * Ramdas et al. [2015] A. Ramdas, S. J. Reddi, B. Póczos, A. Singh, and L. A. Wasserman.  On the decreasing power of kernel and distance based nonparametric hypothesis tests in high dimensions.  In _AAAI_ , 2015.
  * Rietschel et al. [2018] C. Rietschel, J. Yoon, and M. van der Schaar.  Feature selection for survival analysis with competing risks using deep learning.  _arXiv preprint arXiv:1811.09317_ , 2018.
  * Romano et al. [2020] Y. Romano, M. Sesia, and E. Candès.  Deep knockoffs.  _Journal of the American Statistical Association_ , 115(532):1861–1872, 2020.
  * Saber et al. [2019] S. Saber, R. M. Khalil, W. S. Abdo, D. Nassif, and E. El-Ahwany.  Olmesartan ameliorates chemically-induced ulcerative colitis in rats via modulating nf$\kappa$b and nrf-2/ho-1 signaling crosstalk.  _Toxicology and applied pharmacology_ , 364:120–132, 2019\.
  * Salimans et al. [2016] T. Salimans, I. Goodfellow, W. Zaremba, V. Cheung, A. Radford, and X. Chen.  Improved techniques for training gans.  _arXiv preprint arXiv:1606.03498_ , 2016.
  * Santambrogio [2015] F. Santambrogio.  Optimal transport for applied mathematicians.  _Birkäuser, NY_ , 55(58-63):94, 2015.
  * Seguy et al. [2017] V. Seguy, B. B. Damodaran, R. Flamary, N. Courty, A. Rolet, and M. Blondel.  Large-scale optimal transport and mapping estimation.  _arXiv preprint arXiv:1711.02283_ , 2017.
  * Sesia et al. [2019] M. Sesia, C. Sabatti, and E. J. Candès.  Gene hunting with hidden markov model knockoffs.  _Biometrika_ , 106(1):1–18, 2019.
  * Shi et al. [2020a] H. Shi, M. Drton, and F. Han.  Distribution-free consistent independence tests via center-outward ranks and signs.  _Journal of the American Statistical Association_ , pages 1–16, 2020a.
  * Shi et al. [2020b] H. Shi, M. Hallin, M. Drton, and F. Han.  Rate-optimality of consistent distribution-free tests of independence based on center-outward ranks and signs.  _arXiv preprint arXiv:2007.02186_ , 2020b.
  * Sitkin and Pokrotnieks [2018] S. Sitkin and J. Pokrotnieks.  Alterations in polyunsaturated fatty acid metabolism and reduced serum eicosadienoic acid level in ulcerative colitis: is there a place for metabolomic fatty acid biomarkers in ibd?  _Digestive diseases and sciences_ , 63(9):2480–2481, 2018.
  * Smirnov [1939] N. V. Smirnov.  On the estimation of the discrepancy between empirical curves of distribution for two independent samples.  _Bull. Math. Univ. Moscou_ , 2(2):3–14, 1939\.
  * Solakivi et al. [2011] T. Solakivi, K. Kaukinen, T. Kunnas, T. Lehtimäki, M. Mäki, and S. T. Nikkari.  Serum fatty acid profile in subjects with irritable bowel syndrome.  _Scandinavian journal of gastroenterology_ , 46(3):299–303, 2011.
  * Sudarshan et al. [2020] M. Sudarshan, W. Tansey, and R. Ranganath.  Deep direct likelihood knockoffs.  _arXiv preprint arXiv:2007.15835_ , 2020.
  * Székely and Rizzo [2013] G. J. Székely and M. L. Rizzo.  Energy statistics: A class of statistics based on distances.  _Journal of statistical planning and inference_ , 143(8):1249–1272, 2013.
  * Trainor et al. [2017] P. J. Trainor, A. P. DeFilippis, and S. N. Rai.  Evaluation of classifier performance for multiclass phenotype discrimination in untargeted metabolomics.  _Metabolites_ , 7(2):30, 2017.
  * Volkova and Ruggles [2021] A. Volkova and K. V. Ruggles.  Predictive metagenomic analysis of autoimmune disease identifies robust autoimmunity and disease specific microbial signatures.  _Frontiers in microbiology_ , 12:418, 2021.
  * Wald and Wolfowitz [1940] A. Wald and J. Wolfowitz.  On a test whether two samples are from the same population.  _The Annals of Mathematical Statistics_ , 11(2):147–162, 1940.
  * Wilcoxon [1947] F. Wilcoxon.  Probability tables for individual comparisons by ranking methods.  _Biometrics_ , 3(3):119–122, 1947.

##  Appendix A PROOFS

###  A.1 Proof of Lemma 1

###### Proof.

To prove this lemma, using slight abuse of notations, we assume that
$\mu_{X},\mu_{Y}$ are supported on $[a,b]$, and $\mu_{Y}^{s}$ is the shifted
measure defined through $\mu_{Y}^{s}[T]=\mu_{Y}[T-s]$, where $s\in\mathbb{R}$,
$T$ is any measurable set and $T-s=\\{t-s:t\in T\\}$. Let
$\mu_{H}^{s}=\lambda\mu_{X}+(1-\lambda)\mu^{s}_{Y}$ for any $\lambda\in(0,1)$,
denote the mixture measure. Also let $\mathsf{F}_{X},\mathsf{F}_{Y}$ and
$\mathsf{F}_{Y}^{s}$ denote the cumulative distribution functions of
$\mu_{X},\mu_{Y}$ and $\mu_{Y}^{s}$, respectively, where for any
$r\in\mathbb{R}$, by construction $\mathsf{F}_{Y}^{s}(r)=\mathsf{F}_{Y}(r-s)$.
Now suppose $s\in(-\infty,a-b]\cup[b-a,\infty)=S$. Then we must have,

$\displaystyle[a,b]\cap[a+s,b+s]\subset\\{a,b\\}.$

Now, for all $s\in S$, we have,

$\displaystyle\frac{1}{2}\mathop{}\\!\mathrm{RE}_{\lambda}(X,Y^{s})$
$\displaystyle=\int_{-\infty}^{\infty}\big{(}\mathsf{F}_{X}(t)-\mathsf{F}^{s}_{Y}(t)\big{)}^{2}\text{d}\mathsf{H}_{\lambda}^{s}(t)$
$\displaystyle=\int_{-\infty}^{\text{min}(a,a+s)}\\!\\!\\!\\!\\!\\!\\!\\!\big{(}\mathsf{F}_{X}(t)-\mathsf{F}_{Y}^{s}(t)\big{)}^{2}d\text{H}_{\lambda}^{s}(t)+\int_{\text{min}(a,a+s)}^{\text{max}(b,b+s)}\\!\\!\\!\\!\\!\\!\big{(}\mathsf{F}_{X}(t)-\mathsf{F}^{s}_{Y}(t)\big{)}^{2}\text{d}\mathsf{H}_{\lambda}^{s}(t)+\int_{\text{max}(b,b+s)}^{\infty}\\!\\!\\!\\!\\!\\!\\!\\!\\!\big{(}\mathsf{F}_{X}(t)-\mathsf{F}^{s}_{Y}(t)\big{)}^{2}\text{d}\mathsf{H}_{\lambda}^{s}(t)$
$\displaystyle=\int_{\text{min}(a,a+s)}^{\text{max}(b,b+s)}\big{(}\mathsf{F}_{X}(t)-\mathsf{F}^{s}_{Y}(t)\big{)}^{2}\text{d}\mathsf{H}_{\lambda}^{s}(t).$

Case I: $s\leq a-b$,

$\displaystyle\frac{1}{2}\mathop{}\\!\mathrm{RE}_{\lambda}(X,Y^{s})$
$\displaystyle=\int_{a+s}^{b+s}(\mathsf{F}^{s}_{Y}(t)-\mathsf{F}_{X}(t))^{2}\text{d}\mathsf{H}^{s}_{\lambda}(t)+\int_{b+s}^{a}(\mathsf{F}_{Y}^{s}(t)-\mathsf{F}_{X}(t))^{2}\text{d}\mathsf{H}^{s}_{\lambda}(t)+\int_{a}^{b}(\mathsf{F}_{Y}^{s}(t)-\mathsf{F}_{X}(t))^{2}\text{d}\mathsf{H}^{s}_{\lambda}(t)$
$\displaystyle=(1-\lambda)\int_{a+s}^{b+s}(\mathsf{F}_{Y}^{s}(t))^{2}\text{d}\mu_{Y}^{s}(t)+0+\lambda\int_{a}^{b}(1-\mathsf{F}_{X}(t))^{2}\text{d}\mu_{X}(t)$
$\displaystyle=(1-\lambda)\int_{a+s}^{b+s}\\!\\!\\!\\!(\mathsf{F}_{Y}(t-s))^{2}\text{d}\mu_{Y}(t-s)+\lambda\int_{a}^{b}\\!\\!\\!\\!(1-\mathsf{F}_{X}(t))^{2}\text{d}\mu_{X}(t)[\text{assuming
absolute continuity}]$
$\displaystyle=(1-\lambda)\int_{a}^{b}(\mathsf{F}_{Y}(t))^{2}\text{d}\mu_{Y}(t)+\lambda\int_{a}^{b}(1-\mathsf{F}_{X}(t))^{2}\text{d}\mu_{X}(t),$
which is independent of $s$.

Case II: $a-b\leq s\leq 0$;

$\displaystyle\frac{1}{2}\mathop{}\\!\mathrm{RE}_{\lambda}(X,Y^{s})$
$\displaystyle=\int_{a+s}^{a}(\mathsf{F}^{s}_{Y}(t))^{2}\text{d}\mathsf{H}^{s}_{\lambda}(t)+\int_{a}^{b+s}\Big{(}\mathsf{F}^{s}_{Y}(t)-\mathsf{F}_{X}(t)\Big{)}^{2}\text{d}\mathsf{H}^{s}_{\lambda}(t)+\int_{b+s}^{b}\Big{(}1-\mathsf{F}_{X}(t)\Big{)}^{2}\text{d}\mathsf{H}^{s}_{\lambda}(t)$
$\displaystyle=(1-\lambda)\int_{a+s}^{a}\\!\\!\\!\\!\\!(\mathsf{F}^{s}_{Y}(t))^{2}\text{d}\mu^{s}_{Y}(t)+\int_{a}^{b+s}\\!\\!\\!\\!\\!\\!\Big{(}\mathsf{F}^{s}_{Y}(t)-\mathsf{F}_{X}(t)\Big{)}^{2}\\!\\!\Big{(}\lambda\text{d}\mu_{X}(t)+(1-\lambda)\text{d}\mu_{Y}^{s}(t)\Big{)}\\!+\\!\lambda\int_{b+s}^{b}\\!\\!\\!\\!\\!\Big{(}1-\mathsf{F}_{X}(t)\Big{)}^{2}\text{d}\mu_{X}(t)$
this decreases as $s$ moves towards 0.

Case III: $0\leq s\leq b-a$;

$\displaystyle\frac{1}{2}\mathop{}\\!\mathrm{RE}_{\lambda}(X,Y^{s})$
$\displaystyle=\int_{a}^{a+s}(\mathsf{F}_{X}(t))^{2}\text{d}\mathsf{H}^{s}_{\lambda}(t)+\int_{a+s}^{b}\Big{(}\mathsf{F}_{X}(t)-\mathsf{F}^{s}_{Y}(t)\Big{)}^{2}\text{d}\mathsf{H}^{s}_{\lambda}(t)+\int_{b}^{b+s}\Big{(}1-\mathsf{F}^{s}_{Y}(t)\Big{)}^{2}\text{d}\mathsf{H}^{s}_{\lambda}(t)$
$\displaystyle=\lambda\int_{a}^{a+s}\\!\\!\\!\\!\\!\\!\\!(\mathsf{F}_{X}(t))^{2}\text{d}\mu_{X}(t)\\!\\!+\\!\\!\int_{a+s}^{b}\\!\\!\\!\\!\\!\Big{(}\mathsf{F}_{X}(t)-\mathsf{F}^{s}_{Y}(t)\Big{)}^{2}\Big{(}\lambda\text{d}\mu_{X}(t)+(1-\lambda)\text{d}\mu^{s}_{Y}(t)\Big{)}\\!\\!+\\!\\!(1-\lambda)\int_{b}^{b+s}\\!\\!\\!\\!\\!\Big{(}1-\mathsf{F}^{s}_{Y}(t)\Big{)}^{2}\text{d}\mu^{s}_{Y}(t)$
this increases as $s$ moves towards $b-a$.

Case IV : $s\geq b-a$;

$\displaystyle\frac{1}{2}\mathop{}\\!\mathrm{RE}_{\lambda}(X,Y^{s})$
$\displaystyle=\int_{a}^{b}(\mathsf{F}_{X}(t)-\mathsf{F}_{Y}^{s}(t))^{2}\text{d}\mathsf{H}^{s}_{\lambda}(t)+\int_{b}^{a+s}(\mathsf{F}_{X}(t)-\mathsf{F}_{Y}^{s}(t))^{2}\text{d}\mathsf{H}^{s}_{\lambda}(t)+\int_{a+s}^{b+s}(\mathsf{F}_{X}(t)-\mathsf{F}_{Y}^{s}(t))^{2}\text{d}\mathsf{H}^{s}_{\lambda}(t)$
$\displaystyle=\lambda\int_{a}^{b}(\mathsf{F}_{X}(t))^{2}\text{d}\mu_{X}(t)+0+(1-\lambda)\int_{a+s}^{b+s}(1-\mathsf{F}_{Y}^{s}(t))^{2}\text{d}\mu_{Y}^{s}(t)$
$\displaystyle=\lambda\int_{a}^{b}\\!\\!\\!\\!(\mathsf{F}_{X}(t))^{2}\text{d}\mu_{X}(t)+(1-\lambda)\int_{a+s}^{b+s}\\!\\!\\!\\!(1-\mathsf{F}_{Y}^{s}(t-s))^{2}\text{d}\mu_{Y}(t-s)[\text{assuming
absolute continuity}]$
$\displaystyle=\lambda\int_{a}^{b}(\mathsf{F}_{X}(t))^{2}\text{d}\mu_{X}(t)+(1-\lambda)\int_{a}^{b}(1-\mathsf{F}_{Y}(t))^{2}\text{d}\mu_{Y}(t),$
which is independent of $s$.

∎

###  A.2 Proof of Lemma 2

###### Proof.

  * (a)

Using Lemmas 2.2, 2.3 from [Baringhaus and Franz, 2004], for
$\bm{X}_{1},\bm{X}_{2}\overset{i.i.d.}{\sim}\mu_{\bm{X}}$, and
$\bm{Y}_{1},\bm{Y}_{2}\overset{i.i.d.}{\sim}\mu_{\bm{Y}}$, we have,

$\displaystyle 2\mathbb{E}\|\bm{X}_{1}-\bm{Y}_{1}\|-$
$\displaystyle\mathbb{E}\|\bm{X}_{1}-\bm{X}_{2}\|-\mathbb{E}\|\bm{Y}_{1}-\bm{Y}_{2}\|$
$\displaystyle=\gamma_{d}\int_{\mathcal{S}^{d-1}}\int_{\mathbb{R}}\Big{(}\mathbb{P}(\bm{a}^{\top}\bm{X}_{1}\leq
t)-P(\bm{a}^{\top}\bm{Y}_{1}\leq
t)\Big{)}^{2}\text{d}\kappa(\bm{a})\text{d}t,$ (24)

where $\gamma=(2\Gamma(d/2))^{-1}\sqrt{\pi}(d-1)\Gamma\big{(}(d-1)/2\big{)}$
for $d>1$, $\mathcal{S}^{d-1}=\\{\bm{x}\in\mathbb{R}^{d}:\|\bm{x}\|=1\\}$,
$\bm{a}\in\mathcal{S}^{d-1}$, and $\kappa(\cdot)$ denotes the uniform measure
on $\mathcal{S}^{d-1}$. Since,
$\mathbf{R}_{\lambda}^{\varepsilon}(\bm{X}_{1})$,
$\mathbf{R}_{\lambda}^{\varepsilon}(\bm{X}_{2})$,
$\mathbf{R}_{\lambda}^{\varepsilon}(\bm{Y}_{1})$, and
$\mathbf{R}_{\lambda}^{\varepsilon}(\bm{Y}_{2})$ denote the transformed
multivariate random variables corresponds to
$\bm{X}_{1},\bm{X}_{2},\bm{Y}_{1}$ and $\bm{Y}_{2}$, respectively, we can
write the following,

$\displaystyle
2\mathbb{E}\|\mathbf{R}_{\lambda}^{\varepsilon}(\bm{X}_{1})-\mathbf{R}_{\lambda}^{\varepsilon}(\bm{Y}_{1})\|-$
$\displaystyle\mathbb{E}\|\mathbf{R}_{\lambda}^{\varepsilon}(\bm{X}_{1})-\mathbf{R}_{\lambda}^{\varepsilon}(\bm{X}_{2})\|-\mathbb{E}\|\mathbf{R}_{\lambda}^{\varepsilon}(\bm{Y}_{1})-\mathbf{R}_{\lambda}^{\varepsilon}(\bm{Y}_{2})\|$
$\displaystyle=\gamma_{d}\int_{\mathcal{S}^{d-1}}\\!\\!\int_{\mathbb{R}}\\!\\!\
\Big{(}\\!\\!\mathbb{P}(\bm{a}^{\top}\mathbf{R}^{\varepsilon}_{\lambda}(\bm{X}_{1})\leq
t)-\mathbb{P}(\bm{a}^{\top}\mathbf{R}^{\varepsilon}_{\lambda}(\bm{Y}_{1})\leq
t)\Big{)}^{2}\text{d}\kappa(a)\text{d}t.$ (25)

  * (b)

Following the definition of
$\mathop{}\\!\mathrm{sRE}_{\lambda}^{\varepsilon}(\cdot,\cdot)$,

$\displaystyle\mathop{}\\!\mathrm{sRE}^{\varepsilon}_{\lambda}(\bm{X}_{1},\bm{Y}_{1})$
$\displaystyle=\gamma_{d}\int_{\mathbb{R}}\int_{\mathcal{S}^{d-1}}\Big{(}\mathbb{P}(\bm{a}^{\top}\mathbf{R}_{\lambda}^{\varepsilon}(\bm{X}_{1})\leq
t)-\mathbb{P}(\bm{a}^{\top}\mathbf{R}_{\lambda}^{\varepsilon}(\bm{Y}_{1})\leq
t)\Big{)}^{2}d\kappa(\mathbf{a})dt$
$\displaystyle=\gamma_{d}\int_{\mathbb{R}}\int_{\mathcal{S}^{d-1}}\Big{(}\mathbb{P}(\bm{a}^{\top}\mathbf{R}_{\lambda}^{\varepsilon}(\bm{Y}_{1})\leq
t)-\mathbb{P}(\bm{a}^{\top}\mathbf{R}_{\lambda}^{\varepsilon}(\bm{X}_{1})\leq
t)\Big{)}^{2}d\kappa(\mathbf{a})dt$
$\displaystyle=\mathop{}\\!\mathrm{sRE}^{\varepsilon}_{\lambda}(\bm{Y}_{1},\bm{X}_{1}).$
(26)

  * (c)

Assuming $\bm{X}_{1}\overset{d}{=}\bm{Y}_{1}$, we have
$\mathbb{P}(\bm{a}^{\top}\bm{X}_{1}\leq
t)=\mathbb{P}(\bm{a}^{\top}\bm{Y}_{1}\leq t)$ for all
$\bm{a}\in\mathcal{S}^{d-1}$ and $t\in\mathbb{R}$. Given the rank transformed
random variables $\mathbf{R}^{\varepsilon}_{\lambda}(\bm{X}_{1})$,
$\mathbf{R}^{\varepsilon}_{\lambda}(\bm{Y}_{1})$ corresponding to $\bm{X}_{1}$
and $\bm{Y}_{1}$, we have
$\bm{a}^{\top}\mathbf{R}^{\varepsilon}_{\lambda}(\bm{X}_{1})\overset{d}{=}\bm{a}^{\top}\mathbf{R}^{\varepsilon}_{\lambda}(\bm{Y}_{1})$.
Therefore,
$\mathop{}\\!\mathrm{sRE}_{\lambda}^{\varepsilon}(\bm{X}_{1},\bm{Y}_{1})=0$,
if $\bm{X}_{1}\overset{d}{=}\bm{Y}_{1}$.

  * (d)

The unique minimizer $\bm{\pi}^{\varepsilon}$ converges to the optimal Monge-
Kantorovich solution $\bm{\pi}$ with cost function
$c(\cdot)=\|\bm{x}-\bm{y}\|^{p}$ [Carlier et al., 2017],

$\displaystyle\lim_{\varepsilon\rightarrow
0}\mathbf{\bm{\pi}^{\varepsilon}}\rightharpoonup\bm{\pi},$ (27)

where $\rightharpoonup$ denotes convergence w.r.t. the weak topology.
According to the proposition $2.18$ in [Carlier et al., 2017], we have,

$\displaystyle\lim_{\varepsilon\rightarrow
0}\int_{\bm{y}}\bm{y}\text{d}\bm{\pi}^{\varepsilon}(\bm{y}|\bm{x})\rightharpoonup\int_{\bm{y}}\bm{y}\text{d}\bm{\pi}(\bm{y}|\bm{x}),$
(28)

which implies that as $\varepsilon\rightarrow 0$,
$\mathbf{R}^{\varepsilon}_{\lambda}(\cdot)$ will converge to
$\mathbf{R}_{\lambda}(\cdot)$ w.r.t. the weak topology. Consequently,
$\mathop{}\\!\mathrm{sRE}^{\varepsilon}_{\lambda}(\bm{X}_{1},\bm{Y}_{1})$ will
converge to $\mathop{}\\!\mathrm{RE}_{\lambda}(\bm{X}_{1},\bm{Y}_{1})$ as
$\varepsilon\rightarrow 0$.

∎

##  Appendix B List of metabolites discovered using the proposed method

Here, we list the metabolites that are found at least 70 times out of 100
instances of knockoffs generated from the real data using the proposed method.
We add the references (if found) for the selected metabolites.

Table 1: Metabolites selected at an FDR level $\leq 0.05$. “-” denotes either metabolites are not related to IBD or not found in literature as potential biomarkers. METABOLITES | REFERENCES
---|---
1.2.3.4-tetrahydro- -beta-carboline- -1.3-dicarboxylate | [Volkova and Ruggles, 2021]
Urobilin | [Lloyd-Price et al., 2019]
12.13-diHOME | [Levan et al., 2019]
Saccharin | [Qin, 2012]
Salicylate | [Caprilli et al., 2009]
Eicosadienoate | [Sitkin and Pokrotnieks, 2018]
Adrenate | [Lloyd-Price et al., 2019]
Hydrocinnamate | [Price et al., 2017]
Caproate | [Le Gall et al., 2011]
Docosapentaenoate | [Solakivi et al., 2011]
Docosahexaenoate | [Solakivi et al., 2011]
Dodecanedioate | [Lee et al., 2017]
Phenyllactate | [Lavelle and Sokol, 2020]
Eicosatrienoate | [Kuroki et al., 1997]
Olmesartan | [Saber et al., 2019]
Taurolithocholate | [Bauset et al., 2021]
Eicosapentaenoate | [Kuroki et al., 1997]
Arachidonate | [Lloyd-Price et al., 2019]
Myristate | [Fretland et al., 1990]
Tetradecanedioate | -
Porphobilinogen | -
Masilinate | -

##  Appendix C Additional Experiments

###  C.1 Effect of hyperparameters when using sRMMD for knockoff generation

We investigate the effect of the hyperparameters e.g., $\lambda$ (second-order
penalty) and $\delta$ (decorrelation penalty) on the performance of our
proposed knockoff generation technique. We generate knockoffs using different
values of $\lambda$ and $\delta$, and compute the average FDR and power in a
similar way as described in the main paper for each distributional setting. We
then observe the FDR vs. power tradeoff to measure the effect of each
hyperparameter.

#### Multivariate Gaussian.

Figure 4 shows that for all combinations of $\lambda$ and $\delta$, generated
knockoffs can control the FDR at the prespecified level. Using only sRMMD
statistic in the loss function ($\lambda=0,\delta=0$) results into low power
of detection but better FDR control over the small amplitude region. Adding
decorrelation penalty ($\lambda=0,\delta=1$) increases both power and FDR
level. Knockoffs showcase similar power and FDR control if only the second-
order term ($\lambda=1,\delta=0$) is added to the loss function. Addition of
second-order term and decorrelation penalty ($\lambda=1,\delta=1$) marginally
improves the power in the lower amplitude region.

Figure 4: FDR vs amplitude and power vs. amplitude. FDR level is set to 0.1
(black dotted line) for Multivariate Gaussian.

#### Gaussian Mixture Model.

Figure 5 demonstrates that the performance of the proposed method depends on
$\lambda$ and $\delta$ in almost similar way as in multivariate Gaussian case
(Figure 4). sRMMD statistic ($\lambda=0,\delta=0$) can produce knockoffs which
have low power but better control of FDR in the smaller amplitude region.
Decorrelation penalty ($\lambda=0,\delta=1$) increases the power but provides
less control of FDR. A better FDR vs. power tradeoff can be achieved by either
adding only the second-order term ($\lambda=1,\delta=0$) or both
($\lambda=1,\delta=1$).

Figure 5: FDR vs amplitude and power vs. amplitude. FDR level is set to 0.1
(black dotted line) for Gaussian Mixture model.

#### Multivariate Student’s t.

Figure 6 indicates that without the second-order term ($\lambda=0$), knockoffs
gain very low power of detection over the entire amplitude region. With the
increasing $\lambda$, the power increases. However, adding both the second-
order term and the decorrelation penalty $(\lambda=1,\delta=1)$ leads to
higher power but no control of FDR. Reducing decorrelation penalty
($\lambda=1,\delta=0.1$) increases the control of FDR and consequently
achieves better FDR vs. power tradeoff.

Figure 6: FDR vs amplitude and power vs. amplitude. FDR level is set to 0.1
(black dotted line) for multivariate Student’s t distribution.

#### Sparse Gaussian.

From Figure (7), we observe that without the second-order term ($\lambda=0$),
proposed method achieves low power and better control of FDR over the smaller
amplitude region. Increasing decorrelation penalty ($\delta$) marginally
increases the control of FDR.

Figure 7: FDR vs amplitude and power vs. amplitude. FDR level is set to 0.1
(black dotted line) for sparse Gaussian distribution.

In summary, it may therefore be said that proposed method of knockoff
generation is robust to hyperparameters in terms of FDR control in most cases.
However, adding a second-order term or a decorrelation penalty or both can
significantly improve the FDR vs. power tradeoff.

###  C.2 Behavior of the sRE statistic

In this section we investigate the behavior of the sRE statistic w.r.t. the
sample size, dimension and entropic regularization parameter. Similar to
Section 5.1, we consider two Uniform distributions having bounded supports,
$\mathsf{U}[0,1]^{d}$ and $\mathsf{U}[s,s+1]^{d}$ and measure sRE for
$s\in[-10,10]$ between them.

Figure 8: Scaled $\text{RE}(\varepsilon=0)$ and sRE (y-axis) between two
uniform distributions with bounded supports. $n$ and $d$ denote the number of
samples and dimension. RE achieves a constant maximum in the region, where the
two supports are completely disjoint, for all $n,d$. sRE do not saturate for
increased $d$ if $\varepsilon$ and $n$ is larger.

Figure 8 shows the scaled sRE statistic for different combinations of sample
size $n$, dimension $d$ and entropy regularizer $\varepsilon$. sRE with
$\varepsilon=0$ (dubbed as RE) saturates to the constant maximum for all
combinations of $n$ and $d$ over the region where the supports of the two
distributions are completely disjoint. With $\varepsilon=1$, sRE does not
saturate when $d$ is small, but with increasing $d$, it shows similar behavior
as RE. As one can see from the plots that, for increasing $d$, this saturation
can be prevented by using larger $\varepsilon$ and larger sample size $n$.

Comparing the plot of sRMMD statistic (Figure 1, main body), we observe that
both sRE and sRMMD statistics behave in a very similar way w.r.t. $n,d$ and
$\varepsilon$. It may therefore be said that, the behavior of these statistics
depends on how the ranks and the soft ranks are defined.

###  C.3 Use of sRE in learning a generative model

In this section, we validate the use of sRE statistic as a loss function. We
use a similar generative model architecture as described in [Li et al., 2015]
and instead of using maximum mean discrepancy (MMD) [Gretton et al., 2012], we
use sRE as the loss. As recommended, we employ an Autoencoder [Kingma and
Welling, 2013] to project the data into a latent space of reduced dimension.
We set the latent space dimension equal to $d=8$. We train this model with
Adam optimizer [Kingma and Ba, 2014] with the batch size $256$ and the
learning rate $0.001$.

(a) $\varepsilon=0.1$

(b) $\varepsilon=1$

(c) $\varepsilon=5$

Figure 9: MNIST image generation using sRE based generative model.

Figure 9 shows the MNIST image generated using sRE statistic with
$\varepsilon=0.1,1,5$. Though Figure 8 shows the saturation behavior of sRE
statistic with smaller $\varepsilon$ when $d$ is large, which is undesirable
in learning a generative model, the use of Autoencoder facilitates the
generation of valid digits with smaller $\varepsilon$. In addition, this sRE
based generative model can produce unambiguous digits with $\varepsilon=1,5$
with less computational and sample complexity [Mena and Niles-Weed, 2019]
compared to using smaller $\varepsilon=0.1$.