The purpose of this notebook is to allow the numerical experiments described in the paper to be reproduced easily. Running this code may take a few hours on a graphical graphical processing unit.
import numpy as np
from DeepKnockoffs import KnockoffMachine
from DeepKnockoffs import GaussianKnockoffs
import data
import parameters
We model $X \in \mathbb{R}^p $ as a multivariate Student's-t distribution, with $p=100$ and the covariance matrix of an auto-regressive process of order one. The default correlation parameter for this distribution is $\rho =0.5$ and the number of degrees of freedom $\nu = 3$.
# Number of features
p = 100
# Load the built-in multivariate Student's-t model and its default parameters
# The currently available built-in models are:
# - gaussian : Multivariate Gaussian distribution
# - gmm : Gaussian mixture model
# - mstudent : Multivariate Student's-t distribution
# - sparse : Multivariate sparse Gaussian distribution
model = "mstudent"
distribution_params = parameters.GetDistributionParams(model, p)
# Initialize the data generator
DataSampler = data.DataSampler(distribution_params)
Let's sample $n=10000$ observations of $X$. This dataset will be used later to train a deep knockoff machine.
# Number of training examples
n = 10000
# Sample training data
X_train = DataSampler.sample(n)
print("Generated a training dataset of size: %d x %d." %(X_train.shape))
After computing the empirical covariance matrix of $X$ in the training dataset, we can initialize a generator of second-order knockoffs. The solution of the SDP determines the pairwise correlations between the original variables and the knockoffs produced by this algorithm.
# Compute the empirical covariance matrix of the training data
SigmaHat = np.cov(X_train, rowvar=False)
# Initialize generator of second-order knockoffs
second_order = GaussianKnockoffs(SigmaHat, mu=np.mean(X_train,0), method="sdp")
# Measure pairwise second-order knockoff correlations
corr_g = (np.diag(SigmaHat) - np.diag(second_order.Ds)) / np.diag(SigmaHat)
print('Average absolute pairwise correlation: %.3f.' %(np.mean(np.abs(corr_g))))
The default parameters of the machine are set below, as most appropriate for the specific built-in model considered. The figures in the paper were obtained by setting the number of epochs to 1000 and the learning rate to 0.001, while in order to reduce the runtime this notebook uses the values 100 and 0.01 respectively.
# Load the default hyperparameters for this model
training_params = parameters.GetTrainingHyperParams(model)
# Set the parameters for training deep knockoffs
pars = dict()
# Number of epochs
pars['epochs'] = 100
# Number of iterations over the full data per epoch
pars['epoch_length'] = 100
# Data type, either "continuous" or "binary"
pars['family'] = "continuous"
# Dimensions of the data
pars['p'] = p
# Size of the test set
pars['test_size'] = 0
# Batch size
pars['batch_size'] = int(0.5*n)
# Learning rate
pars['lr'] = 0.01
# When to decrease learning rate (unused when equal to number of epochs)
pars['lr_milestones'] = [pars['epochs']]
# Width of the network (number of layers is fixed to 6)
pars['dim_h'] = int(10*p)
# Penalty for the MMD distance
pars['GAMMA'] = training_params['GAMMA']
# Penalty encouraging second-order knockoffs
pars['LAMBDA'] = training_params['LAMBDA']
# Decorrelation penalty hyperparameter
pars['DELTA'] = training_params['DELTA']
# Target pairwise correlations between variables and knockoffs
pars['target_corr'] = corr_g
# Kernel widths for the MMD measure (uniform weights)
pars['alphas'] = [1.,2.,4.,8.,16.,32.,64.,128.]
The machine will be stored in the tmp/ subdirectory for later use and continuously updated during training after each epoch.
# Where to store the machine
checkpoint_name = "tmp/" + model
# Where to print progress information
logs_name = "tmp/" + model + "_progress.txt"
# Initialize the machine
machine = KnockoffMachine(pars, checkpoint_name=checkpoint_name, logs_name=logs_name)
Let's fit the machine to the training data. The value of the loss function on the training will be printed after each epoch, along with other diagnostics based on the MMD, the second moments and the pairwise correlations between variables and knockoffs.
# Train the machine
print("Fitting the knockoff machine...")
machine.train(X_train)