本文作者：Zlin，北京大学数学科学学院应用数学专业博士在读。2022年初，他成为宽德实习生，主要工作集中在变量选择和分布生成模型，以及数值计算的分析等。本篇文章由他向我们介绍Knockoffs，即仿制变量synthetic variables的思想以及其在Variable Selection和False Discovery Rate Control里的应用。
Introduction
仿制变量的思想最早是由斯坦福大学数学系与统计学系的大牛Emmanuel Candès提出并发展的，主要应用于生物统计领域，其思想的发展历程、文章、软件以及源码等可以参考其官网：
https://web.stanford.edu/group/candes/knockoffs
在此，我们对仿制变量的思想做一个Review和Brief introduction。
Knockoffs的产生源于Variable Selection。
比如在线性回归模型$y=X\beta + \epsilon$中，我们有标签(label) $y\in \mathbb{R}^{n}$和特征(feature) $X\in \mathbb{R}^{n\times p}$，其中$X$只有一部分维度对于预测$y$有作用。也就是需要找到集合$\hat{S} = \left\{i\in[p]: \beta_i \neq 0\right\}$，从而可以把冗余的特征去掉，起到降低特征维数的作用。
自然，我们可以直接做最小二乘估计，得到估计值$\hat{\beta}$，从而把$\hat{\beta}$接近为0的特征摒弃，留下$\hat{\beta}$不显著为0的特征。传统的假设检验可以定量地回答$\beta_i=0$是否显著成立这个问题，但这仍然太粗糙了。我们希望对于变量选择的质量能有一个更精细的刻画，并且甚至能对所选变量的错误率(Type I error)做一定的控制。严格地来说，我们希望控制False Discovery Rate：
\begin{equation}\mathrm{FDR} = \mathbb{E}\left[\frac{\#\{j: \beta_j=0~\mathrm{and}~j \in \hat{S}\}}{\#\{j: j \in \hat{S}\}}\right],\end{equation}
（1.1）
也就是需要对变量选择结果中的无用变量的比例有一个定量刻画。构造Knockoffs(或称synthetic variables)的目的就是在变量选择的过程中对FDR做估计。
Literature Review
在下文中，我们始终以$X$表示原始变量，以$\tilde{X}$表示仿制变量。
在[1]中，作者在多元正态分布的假设下，给出了已知$X$，求解$\tilde{X}$的一个实例的解析表达式。在[4]中，直接考虑求解$\tilde{X}$的条件分布$\mathcal{L}(\,\cdot\,|X)$，并依据该分布实现任意多次采样，但文章只能对多元正态分布的情况实现算法。在[2]中，利用Monte-Carlo方法中的metropolized sampling，实现了对任意的已知解析式的分布进行采样的方法，是上一篇文章的进一步推广。
以上文章都必须知悉$X$的分布类型、模型甚至精确的表达式，但结合机器学习算法，我们可以不对$X$做模型假设而直接进行分布的学习。[7]通过深度神经网络直接学习$\tilde{X}$以$X$为条件的分布；[5]则使用GAN学习Knockoffs的条件分布。在[6,9,11]中，作者还采用了VAE等其他可以刻画分布的算法对分布进行学习。在研究Knockoffs的生成方式的同时，[8,10]等文章在降低FDR的Type II error和提升Knockoffs质量等方面有一定的理论分析。
Fixed-X Knockoffs
我们在线性回归模型中引入Knockoffs的概念，并介绍其在回归分析中是如何发挥作用的。
考虑线性模型$y=X\beta + \epsilon$，其中$y\in \mathbb{R}^n$, $X \in \mathbb{R}^{n\times p}$是已知的$n$个样本的标签和特征，$\epsilon \sim \mathcal{N}(0,\sigma^2 I_n)$为标准正态分布白噪声。
在构造仿制变量$\tilde{X}$时，我们希望它通过模仿原变量的相关性结构(correlation structure)，从而在估计FDR时，能发挥有效的作用。我们不妨假设对$X$已经进行了归一化(normalization)，然后估计covariance矩阵$\Sigma = {1\over n}X^\top X$，所构造的$\tilde{X}$需要满足
\begin{equation}\begin{aligned}{1\over n}\tilde{X}^\top \tilde{X} &= \Sigma \,,\\{1\over n}X^\top \tilde{X} &=\Sigma – diag\{s\}.\end{aligned}\end{equation}
（3.1）
其中$s$是一个$p$维非负向量。
$\tilde{X}$保留了$X$的covariance结构，同时$\tilde{X}$与$X$之间也保持$X$的covariance结构，
只有
$${1\over n}X_j^\top\tilde{X}_j=\Sigma_{jj}-s_j = 1-s_j,$$
为了使FDR估计有较好的统计学效力，需要让$s$越大越好，使得$X_j$与$\tilde{X}_j$尽量不相似。否则，极限情况下$s=0$，那么此时$X_j=\tilde{X}_j$，$\tilde{X}_j$与$X_j$对$y$有相同的解释力，这样的构造是没有意义的。
注意到此时
\begin{equation}cov(X,\tilde{X}) = \left[\begin{array}{cc}\Sigma & \Sigma-diag\{s\} \\\Sigma-diag\{s\} & \Sigma\end{array}\right],\end{equation}
（3.2）
对任意的$i,j\in \{1,…,n\}$，交换对应的行和列，covariance矩阵是不变的。也就是说上述构造使得$(X,\tilde{X})$在covariance结构上满足一定的“交换性”。
当$n\geq 2p$时，一种构造方式是选取$s\in \mathbb{R}^p$使得$0\leq diag\{s\} \leq 2\Sigma$，以保证$cov(X,\tilde{X})$的半正定性；再选取$\hat{U}$为$n \times p$矩阵，其$p$个列向量为$n$维空间中的一组标准正交向量，且与$X$张成的子空间正交；再计算Cholesky分解$C^\top C=2diag\{s\}-diag\{s\}\Sigma^{-1} diag\{s\}$；最后取
\begin{equation}\label{FX}\tilde{X} = X(I-\Sigma^{-1}diag\{s\}) + \hat{U}C\end{equation}
（3.3）
构造完毕后，我们将其用于回归。比如考虑回归问题LASSO
\begin{equation*}\hat{\beta} = \mathrm{argmin}_{b\in \mathbb{R}^{2p}} {1\over 2}\|y-[X,\tilde{X}]b\|_2^2 + \lambda \|b\|_1\end{equation*}
并定义$Z_j = \sup\{\lambda: \hat{\beta}_j \neq 0\},~j\in{1,…,2p}$与$\tilde{Z}_j = Z_{j+p}$。注意到$Z_j$其实是一个反映特征重要性的统计量：$Z_j$越大，表明第$j$个变量越重要。
令
$$W_j = f(Z_j,\tilde{Z}_j) = \left\{\begin{array}{cc}Z_j, & Z_j>\tilde{Z}_j,\\-\tilde{Z}_j, & Z_j <\tilde{Z}_j,\\0, &Z_j=\tilde{Z}_j.\end{array}\right.$$
其中$f$关于$Z_j,\tilde{Z}_j$是反对称的。
给定一个threshold: $\lambda$，被选出来的变量即为$S(\lambda) = \{j:W_j\geq \lambda\}$，那么此时对$S(\lambda)$中包含的“零变量”的数量估计为$\#\{j:W_j\leq -\lambda\}$，从而FDR估计为
$$\mathrm{FDR}(\lambda)=\frac{\#\{j:W_j\leq -\lambda\}}{\#\{j:W_j\geq \lambda\}}$$
通过调整threshold的大小，我们就可以通过$\mathrm{FDR}(\lambda)$控制所选择变量的错误率。
图1：将$(Z_j,\tilde{Z}_j)$在坐标中标记
Intuition
我们把$(Z_j, \tilde{Z}_j)$在坐标轴上如图标记出来，红色点是非零变量，黑色点是零变量。那么红色点通常会分布在$y=x$的下方，即$Z_j>\tilde{Z}_j$；黑色点通常会在$y=x$的两侧对称分布。使用仿制变量的intuition就是利用零变量和非零变量的这个分布差异，估计我们所选择的变量里可能混入了多少个零变量。
至此，我们启发性地介绍了仿制变量从构造到使用，再到FDR估计的整个流程，读者对于仿制变量的作用应该有了一个第一印象。由于所述内容几乎不涉及到原理性和推导性的内容，很多细节可能对于读者来说尚且不够清晰，在下一小节中会以更加严谨的形式进行呈现。
Fixed-X Knockoffs仅仅是仿制变量初次被提出时的问题形式。在给定一组样本$X$的情况下，Fixed-X Knockoffs可以给出一组模仿原变量covariance结构的仿制样本$\tilde{X}$。但是这种只关注二阶构造的方式是相对粗糙的，尤其当$X$并不是正态分布的时候，这种构造的作用可能相当有限。下一小节介绍的Model-X Knockoffs直接把$X$作为随机变量，将$X$和$\tilde{X}$的交换性从Cov矩阵延拓到整个联合分布上，并尝试直接刻画条件分布$\tilde{X}\sim \mathcal{L}(\,\cdot\,|X)$。这种方式对分布的刻画更加精细，在理论上更加严谨，拥有更好的统计学效力，我们重点关注Model-X的模型。
Model-X Knockoffs
从名字上看，Model-X就是对$X$进行建模，关注的是$X$的分布。在定义问题的框架时，把$X$视作随机变量，其自身存在一个联合分布，并把样本视作$X$的一组实现(realization)或者实例(instance)，用样本推断$X$的原分布。特别地，在Model-X Knockoffs中，我们用样本$X$推断$\tilde{X}$的条件分布$\mathcal{L}(\,\cdot\,|X)$.
简化起见，我们依然考虑线性模型$y=X\beta + \epsilon$，其中$y\in \mathbb{R}^n$, $X \in \mathbb{R}^{n\times p}$是已知的$n$个样本的标签和特征，$\epsilon \sim \mathcal{N}(0,\sigma^2 I_n)$为标准正态分布白噪声。
在做回归时，我们还引入$\tilde{X}\in \mathbb{R}^{n\times p}$，其满足以下性质：
Definition 4.1 (Model-X Knockoffs) Model-X Knockoffs for the family of random variables $X = (X_1,…,X_p)$ are a new family of random variables $\tilde{X} = (\tilde{X}1,…,\tilde{X}_p)$ constructed with the following two properties:
(1) for any subset $S \subset {1,…,p}$, $$ (X, \tilde{X} )_{swap(S)} \overset{d}{=} (X, \tilde{X} );$$
(2) $\tilde{X} \perp \!\!\!\! \perp Y \,|\, X$ if there is a response $y$. (2) is guaranteed if $\tilde{X}$ is constructed without looking at $y$ .
换言之，构造出来的$\tilde{X}$是一组对$y$没有预测作用的零变量(nulls)。但$\tilde{X}$能起到对照$X$中的零变量的作用，其有如下交换性质
Lemma 4.1 Take any subset $S\subset {1,…,p}$ of nulls of $X$. Then
$$[X,\tilde{X}]~|~y\overset{d}{=} [X,\tilde{X}]_{swap(S)}~|~y.$$
证明可以参考 [4]，在此均略过。
为了简化问题，我们不妨假设$X$满足多元正态分布$\mathcal{N}(0,\Sigma)$，并且假设$\Sigma$是已知的。由性质(1)我们知道$\tilde{X}$也满足多元正态分布。
注1：对于“假设$\Sigma$是已知的”，如果数据中$X$的无标签样本数量远远大于有标签样本数量，或者$X$的样本量本身就是极其丰富的，可以利用丰富的数据对$\Sigma$做一个非常准确的估计。
注2：上述正态分布假设给出了问题最简单的一个形式。在实际情况中，$X$通常不是多元正态分布，因为该假设至少蕴含了：(1) $X$各个分量的边缘分布本身需要是正态分布;(2) $X$各个分量之间是简单的线性关系。换言之，可以在线性模型里使用$X_{-j}$预测$X_j$的值，且残差具有同方差性。但在实际问题中这两个条件通常不仅不满足，而且相去甚远，而$X$的联合分布也不是一个简单的协方差矩阵就可以准确地描述完毕的。因此$\tilde{X}$的分布非常适合使用深度神经网络等容量更大的拟合模型去学习。
正态分布的描述是简单的，只需要知道均值$\mu$和协方差矩阵$\Sigma$，我们就可以完全刻画一个多元正态分布。不妨假设均值为0，
为了满足定义(4.1)中的性质(1)，我们只需要让$cov(X,\tilde{X})$交换任意$i,j\in{1,…,n}$的行列之后保持不变即可。
令$G = cov(X,\tilde{X})$，那么$G$形如
\begin{equation}\label{opt}G = \left[\begin{array}{cc}\Sigma & \Sigma-diag\{s\}\\\Sigma-diag{s} & \Sigma\end{array}\right],\end{equation}
(4.1)
其中$diag{s}$是由$p$维向量构成的对角矩阵，需要我们解一个凸优化问题(semidefinite programming)：
\begin{equation}\label{optpro}\begin{aligned}\mathrm{min}&~\sum_{j=1}^{p} |cov(X_j,\tilde{X}_j)|=\sum_{j=1}^{p} |1-s_j|\\ \mathrm{s.t.}&~s\geq 0\\&~diag\{s\}\leq 2\Sigma\end{aligned}\end{equation}
(4.2)
其中，$X$都预先经过normalization使得$cov(X_j,X_j)=1$；对于$s$的两个限制条件保证$G$仍是一个covariance矩阵，且保证$G$的半正定性质，可以利用矩阵的Schur Complement推导得到。
解上述优化问题的目的是使$X_j$和$\tilde{X}_j$尽量“不相似”。如果不对$s$进行限制，那么在$s=0$时，$\tilde{X}=X$也是一组满足交换性条件的trivial的Knockoffs样本，然而非零变量在Knockoffs里仍然是非零变量，对FDR的估计没有任何作用，这也是定义中设置第(2)条要求的原因。
至此，我们求出了$(X,\tilde{X})$的联合分布$\mathcal{N}(0,G)$，从而给定一个样本$x_i\in \mathbb{R}^p$，我们可以得知仿制变量$\tilde{x}_i$的条件分布
$$p(\tilde{x}|x) \propto exp(-{1\over 2}(x^\top, \tilde{x}^\top)G^{-1}(x^\top, \tilde{x}^\top)^\top)$$
其中假定$G$是正定矩阵，若$G$不是严格正定的矩阵，那么在实操中可以对$G$加一个极小的标量阵$10^{-14}*I_{2p}$，或者直接对$G$求伪逆pinv进行替代。正态分布的采样在计算机中是很容易实现的。
经过一定的代数运算，上述分布密度展开为
\begin{equation}\label{gs}\tilde{X}\sim \mathcal{N}(\,X(I-\Sigma^{-1}diag\{s\})\,,~2diag\{s\}-diag\{s\}\Sigma^{-1}diag\{s\}\,).\end{equation}
(4.3)
读者可以将这个结果与解析式(3.3)进行对比，体会一下Model-X方法的intuition。与此同时，该结果实际上就是对正态分布假设里提到的注(2)的详细解释:
\begin{equation}\begin{aligned}\tilde{\mu}&=X(I-\Sigma^{-1}diag\{s\})\\&=X[(X^\top X)^{-1}X^\top\tilde{X}]\\\tilde{\sigma}^2&=2diag\{s\}-diag\{s\}\Sigma^{-1}diag\{s\}\\&= \tilde{X}^\top (I-X(X^\top X)^{-1}X^\top) \tilde{X}\end{aligned}\end{equation}
(4.4)
如果$\Sigma$就是由样本$X$估计的结果，那么$\tilde{X}$的期望如同使用样本$X$对$\tilde{X}$进行线性回归的预测值，方差如同由回归的残差平方和求均值后得到的。由此可见，如果$X$是多元正态分布，那么$\tilde{X}$也是多元正态分布，并且严格符合上述均值和方差。因此，在现实中，如果$X$的真实分布是一个多元正态分布，那么对Model-X Knockoffs的sampling可以是非常精确和正确的
特别地，如果在实际中，我们有$X$的海量样本，或者是无标签样本数量远大于有标签样本的数量，那么可以利用这些无标签样本对$\Sigma$所一个足够精确的估计，这个估计可以视作对$X$分布的准确的先验知识，求解上述$\tilde{\mu}$和$\tilde{\sigma}^2$。掌握$X$的海量知识，是Model-X方法发挥长处的主要适用场景。
如果$X$的真实分布并不是一个严格的多元正态分布，那么上述过程可以视作构造一套正态分布近似的仿制变量$\tilde{X}$。具体而言做了两个假设：对$\tilde{X}$做多元正态分布的假设，并仅对$\tilde{X}$和$X$的covariance结构做可交换性的条件限制。当然，我们无法奢求这样得到的仿制变量和原变量特别地“相似”，因为该过程实际上强制要求了$\tilde{X}$各个分量的峰度、偏度与正态分布一致，然而$X$的各个分量的峰度、偏度，可能往往并不和正态分布一致。对于这种情况，各类生成模型结合机器学习的强大能力，可以帮助我们突破这些限制，获得更好的效果，具体在下一节中介绍。
仿制变量制备完毕之后，我们接着用其对 FDR 做统计推断。
首先我们需要构造一个特征统计量(feature statistics)
$$W = w([X, \tilde{X}],y) \in \mathbb{R}^p,$$
而$W$的构造又可以分为两步。第一步，构造一个可以直接反映特征的重要性(variable importance)的$2p$维统计量
$$Z=(Z_1,…,Z_p,\tilde{Z}_1,…,\tilde{Z}_p) = \mathbf{t}([X,\tilde{X}],y) \in \mathbb{R}^{2p};$$
第二步，构造反对称函数$f(z,\tilde{z}) = -f(\tilde{z},z)$并定义$W_j = f(Z_j,\tilde{Z}_j)$即可。
举个例子，考虑回归问题 LASSO
\begin{equation}\hat{\beta} = \mathrm{argmin}_{b\in \mathbb{R}^{2p}} {1\over 2}\|y-[X,\tilde{X}]b\|_2^2 + \lambda\|b\|_1\end{equation}
并定义$Z_j = \sup\{\lambda: \hat{\beta}_j \neq 0\},~j\in{1,…,2p}$与$\tilde{Z}_j = Z_{j+p}$.
加入$L^1$正则项之后，若特征$j$的信号强度$|\beta_j|$低于某个函数$\Lambda(\lambda)$，其$\hat{\beta}_j$会被推向为0。$\Lambda(\lambda)$关于$\lambda$是单调增加的，而被选中的变量数目关于$\lambda$是单调减小的。因此，如果$Z_j$比较高，意味着在$\lambda$比较大的时候，变量$j$就已经被选中了，说明变量$j$是重要的。因此$Z_j$是一个反映特征重要性的统计量。
令
$$W_j = f(Z_j,\tilde{Z}_j) = \left\{\begin{array}{cc}Z_j, & Z_j>\tilde{Z}_j,\\-\tilde{Z}_j, & Z_j <\tilde{Z}_j,\\0, &Z_j=\tilde{Z}_j.\end{array}\right.$$
那么我们就得到了一个特征统计量$W=w([X,\tilde{X}],y)$，这里$w$是隐式的。
如果$X_j$是零变量，由引理(4.1)，我们知道$Z_j$和$\tilde{Z}_j$是条件同分布的，从而$P(Z_j>\tilde{Z}_j\,|\,y)=P(\tilde{Z}_j>Z_j\,|\,y)$，即$P(W_j>0\,|\,y)=P(W_j<0\,|\,y)$，故$W_j$的分布在0左右是对称的。
如果$X_j$不是零变量，那么$Z_j>>0$，但$\tilde{X}_j$一定是零变量，$\tilde{Z}_j\approx 0$，从而$W_j$几乎一定是正数（严格来说，是以很高的概率为正数）。
回顾下图，将$(Z_j,\tilde{Z}_j)$点对在坐标中呈现，红色点是非零变量，黑色点是零变量。红色点倾向于满足$Z_j>\tilde{Z}_j$且$Z_j>>0$，黑色点倾向于$(Z_j, \tilde{Z}_j) \approx \mathbf{0}$且在虚线两侧对称分布。
于是我们可以定义FDP (False discovery proportion)
\begin{equation}\mathrm{FDP}(\lambda)=\frac{{j: W_j\leq -\lambda}}{{j: W_j\geq\lambda}},\end{equation}
(4.5)
变量选择的结果为$S(\lambda) = \{j: W_j \geq \lambda\}$，那么$\mathrm{FDP}(\lambda)$是对当前集合中选到的零变量的比例的估计。
取初值$\lambda=\lambda_0$使得第一个变量被选中，然后令$\lambda$单调减小趋近于0，随着正则化系数$\lambda$的减小，被选入的特征数量越来越多，从而零变量的比例越来越高，$\mathrm{FDP}(\lambda)$不断趋近于某个接近1的上界（未必严格单调）。给定一个Type I error的限制$q$ (如q=10%)，我们可以取$\Lambda = \min{\lambda: \mathrm{FDP(\lambda)}\leq q}$和$S(\Lambda)$作为变量选择的结果。
$Z$和$f$的选择还有很多种方式，比如还可以固定$\lambda$，令$Z_j=|\hat{\beta}_j(\lambda)|$, $f(Z_j,\tilde{Z}_j)=Z_j-\tilde{Z}_j$。更多讨论可见[4]。统计量$Z$和反对称函数$f$的构造具有灵活性，读者可以依据实际问题的特点和对统计量解释性与直观性的需要进行大胆设计。该方法的使用还可以推广到GLM等更多模型上，具有广泛性。
上面的叙述均假设$X$满足多元正态分布。对于一般的分布$\mathcal{L}(X)$，同样也可以构造仿制变量符合$\mathcal{L}(\tilde{X}|X)$，使得$\mathcal{L}(X,\tilde{X})=\mathcal{L}(\tilde{X}|X)\mathcal{L}(X)$具有可交换性值。只不过$\mathcal{L}(\tilde{X}|X)$未必存在像(4.3)一样优美简洁的解析解。
最后，读者可能还会关心采样（sampling）的方法。对最一般的情况，[4]提出了Sequential Conditional Independent Pairs，该方法表明只要能对分布$\tilde{X}_j\sim \mathcal{L}(X_j|X_{-j},\tilde{X}_{1:j-1})$采样，然后sequential 地进行下去，就可以得到一组 $(X,\tilde{X})$的实例，其联合分布(的真值)满足Knockoffs定义的可交换性质。[2]提出了Sequential Conditional Exchangeable Pairs的采样过程， 并用Metropolis-Hastings算法（一种Monte-Carlo方法）实现对$\mathcal{L}(X_j|X_{-j},\tilde{X}_{1:j-1})$的采样。
不过上述采样方法均需要我们预先知道$\mathcal{L}(X_j|X_{-j},\tilde{X}_{1:j-1})$的值，在我们面临的更复杂的应用场景中似乎暂时难以派上用场， 因为求解准确的$\mathcal{L}(X_j|X_{-j},\tilde{X}_{1:j-1})$并不是特别容易。
Generating Knockoffs Using Machine Learning
在上文中，我们在线性模型和正态分布的假设下给出了一个knockoff filter of variable selection的应用例子。在这一节里，我们仍然基于线性回归模型，并重点关注仿制变量的生成。上述例子有两个局限，其一是要求$X$必须符合正态分布，其二是要求对分布模型的参数（如例子中的$\Sigma$）有非常丰富的先验知识。
实际情况中，$X$的真实分布一般并不满足严格的正态分布，并且由于维数过高，其真实分布究竟是什么类型都是难以刻画和描述的。
但我们通常拥有的是$X$的丰富样本，于是我们可以考虑用MLP（Multi-layer Perceptron）去表示一个生成函数$$f(x,z),~x\in\mathbb{R}^p,~z\sim \mathcal{N}(0,I_p).$$
从而给定某个样本$x_0$，$f(x_0,z)$把正态分布随机变量$z\sim \mathcal{N}(0,I_p)$变形为$\tilde{x}_0$的条件分布$$\tilde{x}_0 = f(x_0,z)\sim \mathcal{L}(\tilde{x}_0|x_0).$$
在这个问题的框架下，$X$的分布可以是任意的，我们不必知道$X$的分布是什么样的，只需要拥有大量的样本。机器学习强大的表示能力也保证了可以拟合到更奇异的分布类型。
目前，已有[7]用Deep Neural Network实现了上述思想，而[5]使用GAN实现了上述思想，还有作者使用了其他生成模型，但我们在此重点关注和比较DNN和GAN两类网络。
两个算法的Generator均为上述$f(x,z)$的形式，没有太本质的区别。主要是通过设置不同的Loss function来塑造生成的仿制变量的不同性质。
在研究这两个算法前，我们首先需要思考清楚设计算法的需求：我们想要的仿制变量应该是什么样的，应当满足哪些empirical的条件？最核心的有以下三点。
(1) 可交换性：
$$(X, \tilde{X} )_{swap(S)} \overset{d}{=} (X, \tilde{X} )$$
(2) Covariance结构:
\begin{equation}cov(X,\tilde{X}) = \left[\begin{array}{cc}\Sigma & \Sigma-diag\{s\} \\\Sigma-diag\{s\} & \Sigma\end{array}\right],\end{equation}
(3) $X_j$与$\tilde{X}_j$不能太“相似”。
从数学上来说，如果$y$不参与生成$\tilde{X}$的过程，那么满足条件(1)就够了，并且条件(2)其实是(1)的必要不充分条件，似乎显得有些多余。但在线性回归中，变量的covariance结构对结果的影响比较大，因此为了获得更加可靠的结果，在设计loss function时需要对(2)加以强调。
(3)相对前两者，叙述并不那么严谨。其出发点是希望$X_j$和$\tilde{X}_j$所包含的共同的信息越少越好，避免前文中提到的$X_j=\tilde{X}_j$的trivial的情况的产生。最终是希望提高变量选择的Power（即降低Type I error），也就是应当选中的变量尽量都能够被算法选中。
针对这三个条件，在DNN和GAN的两篇文章中，分别提出了不同的两个loss function。我们首先参考一下这两个算法分别是如何设置loss functions，以在数学上描述上述分布的性质，这些是比较有启发性的内容。
【DeepKnockoffs】(DNN) 的 loss function 包含三项:
\begin{equation}J(X,\widetilde{X})=\gamma J_{MMD}(X,\widetilde{X})+\lambda J_{second-order}(X,\widetilde{X}) + \delta J_{decorrelation}(X,\widetilde{X}).\end{equation}
（5.1）
针对条件(1):\begin{equation} J_{MMD}(X,\widetilde{X})=\hat{\mathcal{D}}_{MMD}\left[(X_1,\widetilde{X}_1),(\widetilde{X}_2,X_2)\right] +\hat{\mathcal{D}}_{MMD}\left[(X_1,\widetilde{X}1),(X_2,\widetilde{X}_2)_{swap(S)}\right]\,,\end{equation}
其中$\mathcal{D}_{MMD}$是Maximum Mean Discrepancy （MMD距离），是衡量两个分布的距离的一个度量，后文会简要介绍MMD的定义。$S$则是$\{1,…,p\}$的一个随机子集，每次sample $S$时，任意$i\in\{1,…,p\}$以50%的概率被选进$S$。
针对条件(2):\begin{equation}J_{second-order}(X,\widetilde{X})=\lambda_1 \frac{\|\frac{1}{n}\sum_{i=1}^n(X^i-\widetilde{X}^i)\|^2_2}{p}+ \lambda_2 \frac{\|G_{XX}-G_{\widetilde{X}\widetilde{X}}\|_F^2}{\|G_{XX}\|_F^2}+ \lambda_3 \frac{\|M\circ (G_{XX}-G_{X\widetilde{X}})\|_F^2}{\|G_{XX}\|_F^2}\,,\end{equation}
其中$M=E-I$，$E$是元素全为1的$p$阶方阵，$\circ$为矩阵的逐元素相乘，$\|\cdot\|_F$是矩阵的Frobenius范数，以及
\begin{equation}cov(X,\tilde{X}) = \left[\begin{array}{cc}G_{XX} & G_{X\tilde{X}}\\G_{\tilde{X}X} & G_{\tilde{X}\tilde{X}}\end{array}\right].\end{equation}
针对条件(3):\begin{equation} J_{decorrelation}(X,\widetilde{X})=\sum_{j=1}^p corr(X_j,\widetilde{X}_j). \end{equation}
【KnockoffGAN】(GAN) 的 loss function 也包含三项:
\begin{equation} \mathcal{L}=\min_G\left(\max_D (\mathcal{L}_D) + \lambda\max_M(\mathcal{L}_M) + \mu \max_f (\mathcal{L}_f)\right).\end{equation}
（5.2）
针对条件(1)的可交换性: GAN discriminator
\begin{equation}\mathcal{L}D=\sum_{j=1}^{p} S_j\cdot\log(D_j((X,\tilde{X})_{swap(S)})+(1-S_j)\cdot\log(1-D_j((X,\tilde{X})_{swap(S)})\end{equation}
其中$S$是$p$维的Bernoulli(0.5)随机向量，$S_j=1$表明交换$X$和$\tilde{X}$的第$j$列。$D$是传统GAN中的discriminator: $\mathbb{R}^{2p} \rightarrow [0,1]^{p}$，$D_j=1$表明判别器完全认为$S_j=1$。
针对条件(1)，进一步强调$X\overset{d}{=}\tilde{X}$: WGAN discriminator with gradient penalty
\begin{equation}\mathcal{L}_f = f(X) – f(\tilde{X})-\eta(\| \nabla_{\hat{X}}f(\hat{X})\|_2 – 1)^2\,,\end{equation}
其中函数$f$对应的是Wasserstein度量对偶定义中的函数（即$W(X,\tilde{X})=\sup_{|f|_{Lips}\leq 1} \mathbb{E}_{X}[f(X)]-\mathbb{E}_{\tilde{X}}[f(\tilde{X})]$），在算法中是一个待训练的函数: $\mathbb{R}^{p} \rightarrow \mathbb{R}$，$\hat{X}=\epsilon X + (1-\epsilon)\tilde{X}$。
针对条件(3):使用MINE（Mutual Infomation Neural Estimation互信息神经网络估计），计算互信息$I_{X_j,\widetilde{X}_j}$。互信息越少意味着$X_j$与$\widetilde{X}_j$关系越不密切.
\begin{equation}\begin{aligned}\mathcal{L}_M& =\sum_{j=1}^{p}\hat{I}_{X_j,\tilde{X}_j}\\&=\sum_{j=1}^{p}\left(\sum_{i=1}^{n}T^j(x_j^{(i)},\tilde{x}_j^{(i)}) – \log\left(\sum_{i=1}^{n}\exp(T^j(x_j^{\kappa(i)},\tilde{x}_j^{(i)}))\right)\right)\,,\end{aligned}\end{equation}
其中$T^j$是互信息$I_{X_j,\tilde{X}_j}$的对偶等价表达式中的函数（即$I_{X_j,\tilde{X}_j}=\sup_{T^j}\mathbb{E}_{P_{X,\tilde{X}}}[T^j(X,\tilde{X})]-\log(\mathbb{E}_{P_X\otimes P_{\tilde{X}}}[e^{T^j(X,\tilde{X})}])$，
更详细的内容可见[3]），$T^j: \mathbb{R}^{2p}\rightarrow \mathbb{R}$，也需要通过训练神经网络获得。$n$是样本量，${\kappa(i)}$是一个$n$-permutation。
由于条件(1)蕴含条件(2)，GAN没有另设loss直接约束条件(2)，而是使用两个discriminator强化对条件(1)的约束。原文章的作者为了构造一个形式上统一的min-max问题，loss中的三小项都使用了与原定义对偶的等价表达式，从而需要都需要先求解一个max的优化问题，整体上形成generator $G$与三个函数$D,f,T$相对抗的形式。关于这三项对偶定义的数学阐述并不是本文关心的重点，在此适当略去，感兴趣的读者可自行补充相关知识。
在比较两个算法的表现之前，我们首先需要建立一个评判的标准。评判一组Knockoffs的生成质量如何，是一个比较复杂的问题。在此笔者提出四条标准，仅供参考（注意下文中$X$和$\tilde{X}$都是$n\times p$样本矩阵，从而分布均指经验分布，covariance矩阵均指样本估计得到的矩阵）：
(1) 边缘分布是否足够相似，即$X_j\overset{d}{=}\tilde{X}_j?$。
· 直观地：看分布图像。
· 量化地：看$X_j$和$\tilde{X}_j$的MMD距离（本文具体指RBF kernel MMD）和1-NN分类器的正确率，比较$X$和$\tilde{X}_j$的均值、方差、峰度、偏度，甚至其他有必要的高阶矩。
注1：1-NN分类器指first nearest neighbor classifier。具体而言，将$X$和$\tilde{X}$两组样本混合在一起，对每个样本，以距离其最近的另一个样本所在的类别作为其分类结果，计算其平均的分类正确率。最理想的情况下，正确率应该在50%左右，高于此表明欠拟合，低于此表明过拟合。
注2：这里简要介绍一下MMD。取RBF核$k(x,y)=e^{-{|x-y|_2^2\over 2\xi }}$，那么
$\mathcal{D}_{MMD}(X,Z) = \mathbb{E}_{X,X’}[k(X,X’)] – 2\mathbb{E}_{X,Z}[k(X,Z)] + \mathbb{E}_{Z,Z’}[k(Z,Z’)]，$
其中，$X\overset{i.i.d.}{\sim}X’,~Z\overset{i.i.d.}{\sim}Z’$。直观上，如果$k(\,\cdot\,,\,\cdot\,)$对应的核映射为$x \mapsto \phi(x)$，那么$\mathcal{D}_{MMD}(X,Z)=[\mathbb{E}[\phi(X)]-\mathbb{E}[\phi(Z)]]^2$，注意此处是先求期望，再求差的平方，而不是先求差的平方再求期望。特别地，RBF核映射将$\mathbb{R}^p$中的向量映射到$\mathbb{R}^{\infty}$中，其各个分量$\phi_i(x)$为形如$\phi_i(x_1,x_2,…,x_p)=c_ix_1^{q_i,1}\cdots x_p^{q_i,p}$的高次单项式，因此MMD距离实际上是对$X$和$Z$两个分布的所有“高阶矩”和“高次期望”的差距的综合衡量。
(2) 联合分布是否满足定义中的可交换性质，即$(X, \tilde{X} )_{swap(S)} \overset{d}{=} (X, \tilde{X} )$?
· 量化地：看$(X, \tilde{X} )$和$(X, \tilde{X} )_{swap(S)}$的MMD距离和1-NN分类结果。
(3) $|cov(X_j,\tilde{X}_j)|$需尽量小。
(4) $|cov(X_i,\tilde{X}_j)-cov(X_i,X_j)|$和$|cov(\tilde{X}_i,\tilde{X}_j)-cov(X_i,X_j)|$需尽量小。
特别地，如果要将$\tilde{X}$应用于线性回归中，笔者认为上述标准的优先级是(4)>(3)>(2)>(1)。
需要强调的是，这其中有一个重要的不可忽视的trade-off，即$$\mathrm{minimize} \sum_{j}|cov(X_j,\tilde{X}_j)|$$
与$$\mathrm{minimize} \sum_{i,j}|cov(X_i,\tilde{X}_j)-cov(X_i,X_j)|$$
在训练接近收敛时是明显相互对抗的。在实际经验中，$|cov(\tilde{X}_i,\tilde{X}_j)-cov(X_i,X_j)|$更容易达到一个很小的值，是非常理想的。而 $\sum_{j}|cov(X_j,\tilde{X}_j)|$和 $\sum_{i,j}|cov(X_i,\tilde{X}_j)-cov(X_i,X_j)|$则总是倾向于在某一者结果理想时另一者偏大。原因是对于
\begin{equation}G=\left[\begin{array}{cc}\Sigma & \Sigma-diag\{s\}\\\Sigma-diag\{s\} & \Sigma\end{array}\right] ,\end{equation}
$G$的半正定性是固有的性质，优化问题（4.1）的解$\hat{s}$实际上给出了$s$的一个上界。因此，在实际训练过程中有两种情况会发生：
(1) ${|cov(X_i,\tilde{X}_j)-cov(X_i,X_j)|}_{i,j}$不断趋近于0，而$|cov(X_j,\tilde{X}_j)|$收敛于非负的下界$|1-\hat{s}_j|$。
(2) ${|cov(X_i,\tilde{X}_j)-cov(X_i,X_j)|}_{i,j}$偏差总是很大，难以趋近于0，而$|cov(X_j,\tilde{X}_j)|$可能会低于理论上的下界$|1-\hat{s}_j|$。
换句话说，(2)在最小化$\sum_{j}|cov(X_j,\tilde{X}_j)|$的过程中突破了其正常的下界，而扭曲了$X^\top\tilde{X}$上非对角元的值，而其不能同时趋近于0是由$G$作为半正定矩阵的固有性质所导致的。简而言之，为了满足$G$在行列上的可交换性质，(1)是健康的，(2)是病态的。
这个问题可以通过调整各个loss的权重比值来解决。一般来说，训练过程中，各项loss会逐渐收敛到一个equilibrium point，而这个equilibrium point由各项loss的权重所决定。以DeepKnockoffs为例，当$\lambda_2:\delta$充分大时，训练结果会向形态(1)收敛，$\lambda_2:\delta$极小时，训练结果就可能收敛到(2)这种不正确的情况。
接下来我们来比较一下实际操作中的两个生成模型的质量。
使用两篇文章所公开的原始代码，训练Knockoffs的生成模型，我们比较结果如下。
图2：diag: $|cov(X_j,\tilde{X}_j)|, 1\leq j \leq p = 474$
图3：Distribution of (3)(4) of KnockoffGAN
图 4: Distribution of (3)(4) of DeepKnockoffs
根据图像和上述分析，我们得知KnockoffGAN收敛到了结果(2)，而DeepKnockoffs收敛到了结果(1)。DeepKnockoffs的表现是比KnockoffGAN更好的。并且尽管我们没有刻意地去解优化问题(4.1)，但训练完毕的DeepKnockoffs模型实际上也给出了(4.1)中$s$的一个非常近似的解。
然而，上述结果并不意味着 KnockoffGAN 就一定比 DeepKnockoffs 更差。实际上，我认为 GAN 在生成仿制变量时应该有更大的潜力。
在标准(1)的评判下，GAN的表现显著优于DeepKnockoffs。在检查边缘分布时，DNN生成的仿制变量的各项指标仍然和正态分布太过相似，似乎DNN模型有退化为仅生成正态分布Knockoffs的倾向。既然如此，那我又何必大费周章去训练一个神经网络，直接利用上一节的结果进行正态分布的sampling岂不是更加便捷？可能是因为$\mathcal{D}_{MMD}$刻画分布的能力太过有限，总之，DeepKnockoffs似乎并没有完全发挥出神经网络的强大学习能力，在学习原变量的分布时仍然太过粗糙。
相比之下，KnockoffGAN中引入的第二项loss（即WGAN Discriminator with gradient penalty），在捕捉分布的精细特征上展现出了极其强大的能力。换言之，这项loss使得我们生成的变量不再是一个普普通通的正态分布，而是真正具有了一些奇异随机变量的特征（比如特定的偏度、峰度）。我们可以根据下图直观感受一下两者的差距，其中第一、第四列为同一变量的两个独立同分布样本的经验分布，第二、三列分为DeepKnockoffs和KnockoffGAN生成的仿制变量样本的经验分布。
目前为止，我们已经大致了解了两个生成模型各自的长处和短处。那么一个很自然的想法就是取长补短，以获得一个更好的生成模型。
笔者在此以KnockoffGAN的模型为基础，优化模型的结构。
首先观察到一个现象：当迭代次数充分大时，$|cov(\tilde{X}_i,\tilde{X}_j)-cov(X_i,X_j)|$会收敛到与DeepKnockoffs差不多的水平，因此这一项我们是不必太多担心的。
第二个现象：KnockoffGAN中的三个loss分项，其敏感性（或竞争力）是有差异的。标准(4)中$|cov(X_i,\tilde{X}_j)-cov(X_i,X_j)|$原本 应由Discriminator保证，但MINE太过强势，以至于在$\sum_{j}|cov(X_j,\tilde{X}_j)|$和$|cov(X_i,\tilde{X}_j)-cov(X_i,X_j)|$的相互对抗中，前者获得了胜利。
因此，第一个朴素的想法是降低MINE项loss的权重$\lambda$。然而，笔者经过尝试，发现原算法的Discriminator在引导模型收敛方向时的作用相当有限，大有“摆烂”的嫌疑。
给定一个生成器$G$，一个充分收敛的判别器$D$，其loss(形式上是一个cross-entropy)应当有一个大致的下界$-\log 2$。如果$\mathcal{L}_D<-\log2$，那么意味着这个判别器是无效的，因为判别器的表现太差了，以至于还不如以50%的概率随机给出一个答案。
因此，在对抗生成的过程中，一条健康的$\mathcal{L}_D$的loss曲线，应当至少是先单调增加至严格大于$-\log 2$的某个值，然后再逐渐单调递减趋近于$-\log 2\approx -0.693$。而原模型在实际的训练过程中，$\mathcal{L}_D$从未大于$-0.71$，这显然是不正确的。
对此，笔者做了两个改进，有效地矫正了模型的训练过程。(1)周期性地对$D$加强训练。比如通常在进行一次生成器$G$的更新前，会进行五次判别器$D$的更新；在此基础上，每进行十次$G$的更新之后，额外多进行两百次$D$的更新。(2)适当增加$D$的隐藏层数，提升$D$的拟合容量。
经过上述修改之后，$|cov(X_i,\tilde{X}_j)-cov(X_i,X_j)|$的收敛情况仍然不甚理想，MINE在收敛过程中施加的影响仍然太过显著。
此后，笔者考虑强制加入一项新的loss：$$\mathcal{L}_C = {1\over p^2-p}\left\|{1\over n}(X^\top X-X^\top \tilde{X})\circ M \right\|_F^2.$$
其中$M=E-I$，$E$是元素全为1的$p$阶方阵，$\circ$为矩阵的逐元素相乘。
加入这项损失函数之后，$\mathcal{L}_C$与$\mathcal{L}_M$的竞争博弈更加显著，其loss曲线基本上形如:
此时$|cov(X_i,\tilde{X}_j)-cov(X_i,X_j)|$的收敛情况显著改善，但$|cov(X_j,\tilde{X}_j)|$收敛结果一般，没有像笔者期望的一样收敛于$|1-\hat{s}_j|$。经过多次实验与观察，笔者考虑再加入一项loss:$$\mathcal{L}_d = {1\over p}\sum_{j=1}^{p} cov(X_j,\tilde{X}_j)^2.$$
并直接把$\mathcal{L}_M$摒弃掉，但仍然保留$\mathcal{L}_M$的计算以观察$X_j$与$\tilde{X}_j$的互信息变化情况。
整体上，最终设置的loss为
\begin{equation}\mathcal{L}=\min_G\left(\max_D (\mathcal{L}_D) + \mu \max_f(\mathcal{L}_f) + \lambda\max_M (\mathcal{L}_M) + \gamma\mathcal{L}_C +\zeta \mathcal{L}_d\right)\end{equation}
（5.3）
并用参数来决定保留loss中的哪一些项。其中$\max_f (\mathcal{L}_f)$和$\mathcal{L}_d$可以相互替代，其作用都是为了使生成模型满足条件(3)，故至少有一项的权重不为0即可。
笔者采用\begin{equation}\mathcal{L}=\min_G\left(\max_D (\mathcal{L}_D) + \mu \max_f (\mathcal{L}_f) + \gamma\mathcal{L}_C + \zeta \mathcal{L}_d\right).
\end{equation}
（5.4）
进行参数调试，最终获得了一个在四条标准上结果都相当令人满意的模型。
下文展示一些数值结果。
Loss曲线:
下图中，迭代次数在1500以前时，Discriminator的loss曲线符合我们的预期：在-0.693以上先增后减。尽管iter>1500后突破了-0.693的底线，但实际表现中$(X,\tilde{X})$的可交换性质仍然在继续增强。
因此，empirically speaking，$\mathcal{L}_D$至少在迭代的前期要符合预期，迭代后期即使低于$-\log 2$，可能也不会对生成的仿制变量造成质的影响。
标准(3)和标准(4):
可以看到$cov(X_j,\tilde{X}_j)$的分布和图2中DeepKnockoffs的结果是趋同的，尤其是尖峰的分布位置是相似的，这是两者所隐含的共同的优化问题(4.1)所决定的。
对于标准(1)，参考联合分布的MMD距离$\mathcal{D}_{MMD}(X,\tilde{X})$(越接近0越好):
其中X_origin和X_origin2是原变量的两组独立同分布样本，可以为分布相同时的MMD度量的值作参考。
对于标准(1)，还可以$X$与$\tilde{X}$混合后的1-NN classifier值(越接近0.5越好):
标准(2)的可交换性:
随机采样100组$S$，查看$(X, \tilde{X} )$和$(X, \tilde{X} )_{swap(S)}$的MMD距离和1-NN分类结果。其中第一组$S$特别设置为swap all。
其中MMD距离在$10^{-3}$数量级，是一个可以接受的范围（以$\mathcal{D}_{MMD}(\,\cdot\,,\,\cdot\,)$表格的结果为参照）。而1-NN分类结果表现极其出色。值得一提的是，由于1-NN无法作为loss成为直接的训练目标，在实操中，训练得到的模型的1-NN通常大于0.65，偏离值达0.15以上。当前模型的1-NN偏离值仅在0.02左右，是一个相当优秀的结果。而右图中DeepKnockoffs模型可能存在过拟合或者其他问题。
文章的最后，我们还想简要地回答一个重要的问题：我们如何保证机器学习模型生成的Knockoffs 的质量？达到什么标准的 Knockoffs 生成模型可以放心地被用来进行统计分析？
要严谨地回答这个问题是相当困难的。目前为止，我们关注的重点更多在于“如何制备一套仿制变量”，对于如何评价一套仿制变量的优劣，尚且没有一套完整的体系理论，更别提诸如confidence interval或者confidence bound等推断理论的存在了。
但我们至少有一个 Baseline 可以参考，即上一小节中提到的正态分布近似
\begin{equation}\label{gs2}\tilde{X}\sim \mathcal{N}(\,X(I-\Sigma^{-1}diag\{s\})\,,~2diag\{s\}-diag\{s\}\Sigma^{-1}diag\{s\}\,).\end{equation}
（5.5）
如果依据上述提到的四个标准，机器学习模型都胜过正态分布模型，那么一定程度上我们可以认为机器学习的模型是可靠的。实际中(5.5)的采样是简单的，困难主要在于$s$的计算。笔者尝试了一些已有的求解半定规划的python库，效果均不佳，故考虑直接采用图2中DeepKnockoffs生成的分布的covariance矩阵里的相关信息作为$s$的近似解，会更为靠谱。
经过笔者尝试，若只考虑线性回归模型，那么(5.5)近似的Knockoffs已经可以发挥很不错的统计学效力，尽管在四条标准上表现均显著弱于KnockoffGAN。总的来说，正态分布的采样便捷，但是求解$s$相对困难；KnockoffGAN的训练和参数调整较为困难，但优点是不需要再额外求解一个优化问题，并且无论是最后的采样以及隐含的$s$的解，结果都比较可靠。
References
[1] Rina Foygel Barber and Emmanuel J Candès. Controlling the false discovery rate via knockoffs. The Annals of Statistics, 43(5):2055–2085, 2015.
[2] Stephen Bates, Emmanuel Candès, Lucas Janson, and Wenshuo Wang. Metropolized knockoff sampling. Journal of the American Statistical Association, 116(535):1413–1427, 2021.
[3] Mohamed Ishmael Belghazi, Aristide Baratin, Sai Rajeshwar, Sherjil Ozair, Yoshua Bengio, Aaron Courville, and Devon Hjelm. Mutual information neural estimation. In International conference on machine learning, pages 531–540. PMLR, 2018.
[4] Emmanuel Candès, Yingying Fan, Lucas Janson, and Jinchi Lv. Panning for gold: ’modelx’ knockoffs for high dimensional controlled variable selection. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 80(3):551–577, 2018.
[5] James Jordon, Jinsung Yoon, and Mihaela van der Schaar. Knockoffgan: Generating knockoffs for feature selection using generative adversarial networks. In International Conference on Learning Representations, 2018.
[6] Ying Liu and Cheng Zheng. Deep latent variable models for generating knockoffs. Stat,8(1):e260, 2019.
[7] Yaniv Romano, Matteo Sesia, and Emmanuel Candès. Deep knockoffs. Journal of the American Statistical Association, 115(532):1861–1872, 2020.
[8] Asher Spector and Lucas Janson. Powerful knockoffs via minimizing reconstructability. The Annals of Statistics, 50(1):252–276, 2022.
[9] Mukund Sudarshan, Wesley Tansey, and Rajesh Ranganath. Deep direct likelihood knockoffs. Advances in neural information processing systems, 33:5036–5046, 2020.
[10] Wenshuo Wang and Lucas Janson. A power analysis of the conditional randomization test and knockoffs. arXiv preprint arXiv:2010.02304, 2020.
[11] Zifan Zhu, Yingying Fan, Yinfei Kong, Jinchi Lv, and Fengzhu Sun. Deeplink: Deep learning inference using knockoffs with applications to genomics. Proceedings of the National Academy of Sciences, 118(36):e2104683118, 202
