Robust Estimation and Generative Adversarial Nets

Introduction

In the setting of Huber’s $\epsilon$-contamination model , one has i.i.d observations

and the goal is to estimate the model parameter $\theta$. Under the data generating process (1), each observation has a $1-\epsilon$ probability to be drawn from $P_{\theta}$ and the other $\epsilon$ probability to be drawn from the contamination distribution $Q$. The presence of an unknown contamination distribution poses both statistical and computational challenges. For example, consider a normal mean estimation problem with $P_{\theta}=N(\theta,I_{p})$. Due to the contamination of data, the sample average, which is optimal when $\epsilon=0$, can be arbitrarily far away from the true mean if $Q$ charges a positive probability at infinity. Moreover, even robust estimators such as coordinatewise median and geometric median are proved to be suboptimal under the setting of (1) . The search for both statistically optimal and computationally feasible procedures has become a fundamental problem in areas including robust statistics and theoretical computer science.

the maximizer of Tukey’s halfspace depth. Despite the statistical optimality of Tukey’s median, computation of (2) is not tractable. In fact, even an approximate algorithm takes $O(e^{Cp})$ in time .

Recent developments in theoretical computer science are focused on the search of computationally tractable algorithms for estimating $\theta$ under Huber’s $\epsilon$-contamination model (1). The success of the efforts started from two fundamental papers , where two different but related computational strategies “iterative filtering” and “dimension halving” were proposed to robustly estimate the normal mean. These algorithms can provably achieve the minimax rate $\frac{p}{n}\vee\epsilon^{2}$ up to a poly-logarithmic factor in polynomial time. The main idea behind the two methods is the fact that a good robust moment estimator can be certified efficiently by higher moments. This idea was later further extended to develop robust and computable procedures for various other problems.

Compared with these computationally feasible procedures proposed in the recent literature for robust estimation, Tukey’s median (2) and other depth-based estimators have some indispensable advantages in terms of their statistical properties. First, the depth-based estimators have clear objective functions that can be interpreted from the perspective of projection pursuit . Second, the depth-based procedures are adaptive to nuisance parameters in the models such as covariance structures, contamination proportion, and error distributions . In comparison, many of the computationally feasible procedures for robust mean estimation in the literature rely on the knowledge of covariance matrix, and sometimes the order of the contamination proportion as well. Even though these assumptions can be relaxed, nontrivial modifications of the algorithms are required for such extensions and sometimes statistical error rates will be affected. Last but not least, Tukey’s depth and other depth functions are mostly designed for robust quantile estimation, while the recent advancements in the theoretical computer science literature are all focused on robust moments estimation. Although this is not an issue when it comes to the problem of normal mean estimation, the difference becomes fundamental for robust location estimation under general settings such as elliptical distributions where moments do not necessarily exist. For a thorough overview of statistical properties of depth-based estimators, we refer the readers to .

Given the desirable statistical properties discussed above, this paper is focused on the development of computational strategies of depth-like procedures. Our key observation is that robust estimators that are maximizers of depth functions, including halfspace depth, regression depth and covariance matrix depth, can all be derived under the framework of $f$-GAN . As a result, these depth-based estimators can be viewed as minimizers of variational lower bounds of the total variation distance between the empirical measure and the model distribution. This observation allows us to leverage the recent developments in the deep learning literature to compute these variational lower bounds through neural network approximations. Our theoretical results give insights on how to choose appropriate neural network classes that lead to minimax optimal robust estimation under Huber’s $\epsilon$-contamination model. The main contributions of the paper are listed below.

We identify an important subclass of $f$-GAN, called $f$-Learning (Section 2.1), which helps us to unify the understandings of various depth-based estimators, GANs, and MLE in a single framework. The connection between depth functions and $f$-GAN allows us to develop depth-like estimators that not only share good statistical properties of (2), but can also be trained by stochastic gradient ascent/descent algorithms.

In order to choose an appropriate discriminator class for robust estimation, we establish a relation between (JS)-GAN optimization and feature matching (Proposition 3.1). This implies the necessity of hidden layers of neural network structures used in the GAN training. A neural network class without hidden layer is equivalent to matching linear features, and is thus not suitable for robust estimation.

We prove that rate-optimal robust location estimation for both Gaussian distribution (Theorem 3.1 for TV-GAN and Theorem 3.2 for JS-GAN with bounded activations, and Theorem 4.1 for deep ReLU networks) and the general family of elliptical distributions (Theorem 5.1) can be achieved by GANs that use neural network discriminator classes with appropriate structures and regularizations. Extensive numerical experiments are conducted to verify our theoretical findings and show that these procedures can be computed in practice.

Our work is also related to the recent literature on the investigation of statistical properties of GAN. For example, nonparametric density estimation using GAN is studied by . Provable guarantees of learning Gaussian distributions with quadratic discriminators are established by . Theoretical guarantees of learning Gaussian mixtures, exponential families and injective neural network generators are obtained by . The result we obtain in this paper is the first theoretical guarantee of GAN in robust estimation under Huber’s $\epsilon$-contamination model.

The rest of the paper is organized as follows. In Section 2, we introduce an $f$-Learning framework and discuss the connection between robust estimation and $f$-GAN. The theoretical results of robust Gaussian mean estimation using $f$-GAN are given in Section 3. Results for deep ReLU networks are given in Section 4. An extension to robust location estimation for the family of Elliptical distributions is presented in Section 5 that includes both Gaussian distribution and Cauchy distribution whose moments do not exist. In Section 6, we present extensive numerical studies of the proposed procedures. Section 7 collects some discussions on the results of the paper and several possible extensions of the work. Finally, all the technical proofs are given in Section 8.

Robust Estimation and f𝑓f-GAN

We start with the definition of $f$-divergence . Given a strictly convex function $f$ that satisfies $f(1)=0$, the $f$-divergence between two probability distributions $P$ and $Q$ is defined by

Here, we use $p(\cdot)$ and $q(\cdot)$ to stand for the density functions of $P$ and $Q$ with respect to some common dominating measure. For a fully rigorous definition, see . Let $f^{*}$ be the convex conjugate of $f$. That is, $f^{*}(t)=\sup_{u\in\text{dom}_{f}}(ut-f(u))$. A variational lower bound of (3) is

Note that the inequality (4) becomes an equality whenever the class $\mathcal{T}$ contains the function $f^{\prime}\left(p/q\right)$ . For notational simplicity, we also use $f^{\prime}$ for an arbitrary element of the subdifferential when the derivative does not exist. With i.i.d. observations $X_{1},...,X_{n}\sim P$, the variational lower bound (4) naturally leads to the following learning method

The formula (5) is a powerful and general way to learn the distribution $P$ from its i.i.d. observations. It is known as $f$-GAN , an extension of GAN , which stands for generative adversarial nets. The idea is to find a $\widehat{P}$ so that the best discriminator $T$ in the class $\mathcal{T}$ cannot tell the difference between $\widehat{P}$ and the empirical distribution $\frac{1}{n}\sum_{i=1}^{n}\delta_{X_{i}}$.

Our $f$-Learning framework is based on a special case of the variational lower bound (4). That is,

where $\widetilde{q}(\cdot)$ stands for the density function of $\widetilde{Q}$. Note that here we allow the class $\widetilde{\mathcal{Q}}_{Q}$ to depend on the distribution $Q$ in the second argument of $D_{f}(P\|Q)$. Compare (6) with (4), and it is easy to realize that (6) is a special case of (4) with

Moreover, the inequality (6) becomes an equality as long as $P\in\widetilde{\mathcal{Q}}_{Q}$. The sample version of (6) leads to the following learning method

The learning method (8) will be referred to as $f$-Learning in the sequel. It is a very general framework that covers many important learning procedures as special cases. For example, consider the special case where $\widetilde{\mathcal{Q}}_{Q}=\widetilde{\mathcal{Q}}$ independent of $Q$, $\mathcal{Q}=\widetilde{\mathcal{Q}}$, and $f(x)=x\log x$. Direct calculations give $f^{\prime}(x)=\log x+1$ and $f^{*}(t)=e^{t-1}$. Therefore, (8) becomes

which is the maximum likelihood estimator (MLE).

The $f$-Learning (8) is related to but is different from the rho-estimation framework . The unpenalized version of the rho-estimator is defined by

where $\psi:[0,+\infty]\rightarrow$ is a non-decreasing function that satisfies $\psi(x)=-\psi(1/x)$. The rho-estimation framework has a different motivation. The function $\psi$ is designed to generalize the logarithmic function (which leads to the MLE) so that the induced procedure is robust to a Hellinger model misspecification. On the other hand, the $f$-Learning (8) is directly derived from a variational lower bound of the $f$-divergence.

2 TV-Learning and Depth-Based Estimators

The TV-Learning (9) is a very useful tool in robust estimation. A closely related idea was previously explored by . We illustrate its applications by several examples of depth-based estimators.

Letting $r\rightarrow 0$, we obtain (2), the exact formula of Tukey’s median. A traditional understanding of Tukey’s median is that (2) maximizes the halfspace depth so that $\widehat{\theta}$ is close to the center of all one-dimensional projections of the data. In the $f$-Learning framework, $N(\widehat{\theta},I_{p})$ is understood to be the minimizer of a variational lower bound of the total variation distance.

The next example is a linear model $y|X\sim N(X^{T}\theta,1)$. Consider the following classes

Here, $P_{y,X}$ stands for the joint distribution of $y$ and $X$. The two classes $\mathcal{Q}$ and $\widetilde{\mathcal{Q}}_{\eta}$ share the same marginal distribution $P_{X}$ and the conditional distributions are specified by $N(X^{T}\eta,1)$ and $N(X^{T}\widetilde{\eta},1)$, respectively. Follow the same derivation of Tukey’s median, let $r\rightarrow 0$, and we obtain

which is the estimator that maximizes the regression depth proposed by . It is worth noting that the derivation of (11) does not depend on the marginal distribution $P_{X}$.

The last example is on covariance matrix estimation. For this task, we set $\mathcal{Q}=\{N(0,\Gamma):\Gamma\in\mathcal{E}_{p}\}$, where $\mathcal{E}_{p}$ is the class of all $p\times p$ covariance matrices. Inspired by the derivations of Tukey depth and regression depth, it is tempting to choose $\widetilde{\mathcal{Q}}_{\Gamma}$ in the neighborhood of $N(0,\Gamma)$. However, a naive choice would lead to a definition that is not even Fisher consistent. We propose a rank-one neighborhood, given by

under the limit $r\rightarrow 0$. Even though the definition of (12) is given by a rank-one neighborhood of the inverse covariance matrix, the formula (2.2) can also be derived with $\widetilde{\Gamma}^{-1}=\Gamma^{-1}+\widetilde{r}uu^{T}$ in (12) replaced by $\widetilde{\Gamma}=\Gamma+\widetilde{r}uu^{T}$ by applying the Sherman-Morrison formula. A similar formula to (2.2) in the literature is given by

3 From f𝑓f-Learning to f𝑓f-GAN

The depth-based estimators (2), (11) and (15) are all proved to be statistically optimal under Huber’s contamination model . This shows the importance of TV-Learning in robust estimation. However, it is well-known that depth-based estimators are very hard to compute , which limits their applications only for very low-dimensional problems. On the other hand, the general $f$-GAN framework (5) has been successfully applied to learn complex distributions and images in practice . The major difference that gives the computational advantage to $f$-GAN is its flexibility in terms of designing the discriminator class $\mathcal{T}$ using neural networks compared with the pre-specified choice (7) in $f$-Learning. While $f$-Learning provides a unified perspective in understanding various depth-based procedures in robust estimation, we can step back into the more general $f$-GAN for its computational advantages, and to design efficient computational strategies. However, there are at least two questions that are unclear:

How to choose the function $f$ that leads to robust learning procedures which are easy to optimize?

How to specify the discriminator class to learn the parameter of interest with minimax rate under Huber’s $\epsilon$-contamination model?

In the rest of the paper, we will study a robust mean estimation problem in detail to answer these questions and illustrate the power of $f$-GAN in robust estimation.

Robust Mean Estimation via GAN

We start with the total variation GAN (TV-GAN) with $f(x)=(x-1)_{+}$ in (5). For the Gaussian location family, (5) can be written as

with $T(x)=D(x)$ in (5). Now we need to specify the class of discriminators $\mathcal{D}$ to solve the classification problem between $N(\eta,I_{p})$ and the empirical distribution $\frac{1}{n}\sum_{i=1}^{n}\delta_{X_{i}}$. One of the simplest discriminator classes is the logistic regression,

With $D(x)={\sf sigmoid}(w^{T}x+b)$ in (17), the procedure (16) can be viewed as a smoothed version of TV-Learning (9). To be specific, the sigmoid function ${\sf sigmoid}(w^{T}x+b)$ tends to an indicator function as $\|w\|\rightarrow\infty$, which leads to a procedure very similar to (10). In fact, the class (17) is richer than the one used in (10), and thus (16) can be understood as the minimizer of a sharper variational lower bound than that of (10).

Assume $\frac{p}{n}+\epsilon^{2}\leq c$ for some sufficiently small constant $c>0$. With i.i.d. observations $X_{1},...,X_{n}\sim(1-\epsilon)N(\theta,I_{p})+\epsilon Q$, the estimator $\widehat{\theta}$ defined by (16) satisfies

Though TV-GAN can achieve the minimax rate, it may suffer from optimization difficulties especially when the distributions $Q$ and $N(\theta,I_{p})$ are far away from each other. The main obstacle is, with optimization based on gradient, the discriminator may be stuck in a local maximum which will then pass wrong signals to the generator. We illustrate this point with a simple one-dimensional example in Figure 1, where samples are drawn from $(1-\epsilon)N(1,1)+\epsilon N(10,1)$ with $\epsilon=0.2$, and we optimize (16) via alternative gradient ascent and descent shown in Algorithm 1. Even with a good initialization, TV-GAN in the form of (16) will continuously increase the value of $\eta$ (from the light area to the dark area in the heatmap) if $w$ cannot achieve its global maximum, and thus fails to learn the saddle point. However, it is almost impossible for ${w}$ to correct its way from $w\to\infty$ to $w\to-\infty$ simply by the information of its local gradient. In comparison, the landscape becomes better when $Q$ and $N(\theta,I_{p})$ are close, where the signal passed to the generator becomes weak before being stuck in the local maximum, as shown in Figure 2.

2 Results for JS-GAN

Given the intractable optimization property of TV-GAN, we next turn to Jensen-Shannon GAN (JS-GAN) with

with $T(x)=\log D(x)$ in (5). This is exactly the original GAN specialized to the normal mean estimation problem. The advantages of JS-GAN over other forms of GAN have been studied extensively in the literature .

Before presenting theoretical properties of (18), we first show a simple numerical result that implies important consequences on the choice of the discriminator class $\mathcal{D}$. Consider i.i.d. observations drawn from the one-dimensional contamination model $(1-\epsilon)N(\theta,1)+\epsilon N(t,1)$ with $\theta=1$ and $\epsilon=0.2$. We consider two estimators in the form of (18) that use different discriminator classes. The first one is the same logistic regression class defined in (17), and the second one is the class of neural networks with one hidden layer. Then, the values of the two estimators are plotted against $t$ in Figure 3. It is clear that the two estimators have completely different behaviors. For the estimator trained by JS-GAN using a logistic regression discriminator class, it is always close to $0.2+0.8t$, which is the grand mean of the entire distribution $(1-\epsilon)N(\theta,1)+\epsilon N(t,1)$. Thus, the estimator is not robust, and its deviation from $\theta$ will become arbitrarily large when the value of $t$ is increased. On the other hand, with an extra hidden layer built into the neural nets, the second estimator is always close to the mean $\theta$ that we want to learn, regardless of the value of $t$. The green curve in Figure 3 first increases as $t$ increases, but it eventually converges to $\theta=1$ as $t$ further increases. The hardest contamination distribution $N(t,1)$ is the one with a $t$ that is not far away from $\theta$, which is well predicted by the minimax theory of robust estimation .

In other words, $P$ and $Q$ are distinguished by a logistic regression classifier that uses the feature $g(X)$. It is easy to see that ${\sf JS}_{g}(P,Q)$ is a variational lower bound of the original Jensen-Shannon divergence. The key property of ${\sf JS}_{g}(P,Q)$ is given by the following proposition.

Assume $\mathcal{W}$ is a convex set that contains an open neighborhood of . Then, ${\sf JS}_{g}(P,Q)=0$ if and only if $E_{P}g(X)=E_{Q}g(X)$.

Define $F(w)=E_{P}\log{\sf sigmoid}(w^{T}g(X))+E_{Q}\log(1-{\sf sigmoid}(w^{T}g(X)))+\log 4$, so that ${\sf JS}_{g}(P,Q)=\sup_{w\in\mathcal{W}}F(w)$. The gradient and Hessian of $F(w)$ are given by

Therefore, $F(w)$ is concave in $w$, and $\sup_{w\in\mathcal{W}}F(w)$ is a convex optimization with a convex $\mathcal{W}$. Suppose ${\sf JS}_{g}(P,Q)=0$. Then $\sup_{w\in\mathcal{W}}F(w)=0=F(0)$, which implies $\nabla F(0)=0$, and thus we have $E_{P}g(X)=E_{Q}g(X)$. Now suppose $E_{P}g(X)=E_{Q}g(X)$, which is equivalent to $\nabla F(0)=0$. Therefore, $w=0$ is a stationary point of a concave function, and we have ${\sf JS}_{g}(P,Q)=\sup_{w\in\mathcal{W}}F(w)=F(0)=0$. ∎

The proposition asserts that ${\sf JS}_{g}(\cdot,\cdot)$ cannot distinguish $P$ and $Q$ if the feature $g(X)$ has the same expected value under the two distributions. This generalized moment matching effect has also been studied by for general $f$-GANs. However, the linear discriminator class considered in is parameterized in a different way compared with the discriminator class here.

We will show rigorously that a neural net with one hidden layer is sufficient to make (18) robust and optimal. Consider the following class of discriminators,

Consider the estimator $\widehat{\theta}$ defined by (18) with $\mathcal{D}$ specified by (19). Assume $\frac{p}{n}+\epsilon^{2}\leq c$ for some sufficiently small constant $c>0$, and set $\kappa=O\left(\sqrt{\frac{p}{n}}+\epsilon\right)$. With i.i.d. observations $X_{1},...,X_{n}\sim(1-\epsilon)N(\theta,I_{p})+\epsilon Q$, we have

Deep ReLU Networks

Combining with the last sigmoid layer, we obtain the following discriminator class,

Note that all the activation functions are ${\sf ReLU}(\cdot)$ except that we use ${\sf sigmoid}(\cdot)$ in the last layer of feature map $g(\cdot)$. A theoretical guarantees of the class defined above is given by the following theorem.

Assume $\frac{p\log p}{n}\vee\epsilon^{2}\leq c$ for some sufficiently small constant $c>0$. Consider i.i.d. observations $X_{1},...,X_{n}\sim(1-\epsilon)N(\theta,I_{p})+\epsilon Q$ and the estimator $\widehat{\theta}$ defined by (18) with $\mathcal{D}=\mathcal{F}_{L}^{H}(\kappa,\tau,B)$ with $H\geq 2p$, $2\leq L=O(1)$, $2\leq B=O(1)$, and $\tau=\sqrt{p\log p}$. We set $\kappa=O\left(\sqrt{\frac{p\log p}{n}}+\epsilon\right)$. Then, for the estimator $\widehat{\theta}$ defined by (18) with $\mathcal{D}={\mathcal{F}}_{L}^{H}(\kappa,\tau,B)$, we have

with high probability. Hence, we can use $\widehat{\theta}+\widetilde{\theta}$ as the final estimator to achieve the same rate in Theorem 4.1.

On the other hand, our experiments show that this preprocessing step is not needed. We believe that the assumption $\|\theta\|_{\infty}\leq\sqrt{\log p}$ is a technical artifact in the analysis of the Rademacher complexity. It can probably be dropped by a more careful analysis.

Elliptical Distributions

For a unit vector $v$, let the density function of $\xi v^{T}U$ be $h$. Note that $h$ is independent of $v$ because of the symmetry of $U$. Then, there is a one-to-one relation between the distribution of $\xi$ and $h$, and thus the triplet $(\theta,\Sigma,h)$ fully parametrizes an elliptical distribution.

Note that $h$ and $\Sigma=AA^{T}$ are not identifiable, because $\xi A=(c\xi)(c^{-1}A)$ for any $c>0$. Therefore, without loss of generality, we can restrict $h$ to be a member of the following class

This makes the parametrization $(\theta,\Sigma,h)$ of an elliptical distribution fully identifiable, and we use $EC(\theta,\Sigma,h)$ to denote an elliptical distribution parametrized in this way.

where $\mathcal{E}_{p}(M)$ is the set of all positive semi-definite matrix with spectral norm bounded by $M$.

Consider the estimator $\widehat{\theta}$ defined above with $\mathcal{D}$ specified by (19). Assume $M=O(1)$, $\frac{p}{n}+\epsilon^{2}\leq c$ for some sufficiently small constant $c>0$, and set $\kappa=O\left(\sqrt{\frac{p}{n}}+\epsilon\right)$. With i.i.d. observations $X_{1},...,X_{n}\sim(1-\epsilon)EC(\theta,\Sigma,h)+\epsilon Q$, we have

Note that Theorem 5.1 guarantees the same convergence rate as in the Gaussian case for all elliptical distributions. This even includes multivariate Cauchy where mean does not exist. Therefore, the location estimator (20) is fundamentally different from , which is only designed for robust mean estimation.

To achieve rate-optimality for robust location estimation under general elliptical distributions, the estimator (20) is different from (18) only in the generator class. They share the same discriminator class (19). This underlines an important principle for designing GAN estimators: the overall statistical complexity of the estimator is only determined by the discriminator class.

The estimator (20) also outputs $(\widehat{\Sigma},\widehat{h})$, but we do not claim any theoretical property for $(\widehat{\Sigma},\widehat{h})$ in this paper.

Numerical Experiments

In this section, we give extensive numerical studies of robust mean estimation via GAN. After introducing the implementation details in Section 6.1, we verify our theoretical results on minimax estimation with both TV-GAN and JS-GAN in Section 6.2. Comparison with other methods on robust mean estimation in the literature is given in Section 6.3. The effects of various network structures are studied in Section 6.4. Finally, adaptation to unknown covariance structure and elliptical distributions are investigated in Section 6.5 and Section 6.6.

The implementation for JS-GAN is given in Algorithm 1, and a simple modification of the objective function leads to that of TV-GAN. A PyTorch implementation is available at https://github.com/zhuwzh/Robust-GAN-Center or https://github.com/yao-lab/Robust-GAN-Center. Several important implementation details are listed below.

How to tune parameters? The choice of learning rates is crucial to the convergence rate, but the minimax game is hard to evaluate. We propose a simple strategy to tune hyper-parameters including the learning rates. Suppose we have estimators $\widehat{\theta}_{1},\ldots,\widehat{\theta}_{M}$ with corresponding discriminator networks $D_{\widehat{w}_{1}}$,…, $D_{\widehat{w}_{M}}$. Fixing $\eta=\widehat{\theta}$, we further apply gradient descent to $D_{w}$ with a few more epochs (but not many in order to prevent overfitting, for example 10 epochs) and select the $\widehat{\theta}$ with the smallest value of the objective function (18) (JS-GAN) or (16) (TV-GAN). We note that training discriminator and generator alternatively usually will not suffer from overfitting since the objective function for either the discriminator or the generator is always changing. However, we must be careful about the overfitting issue when training the discriminator alone with a fixed $\eta$, and that is why we apply an early stopping strategy here. Fortunately, the experiments show that if the structures of networks are same (then of course, the dimensions of the inputs are same), the choices of hyper-parameters are robust to different models.

When to stop training? Judging convergence is a difficult task in GAN trainings, since sometimes oscillation may occur. In computer vision, people often use a task related measure and stop training once the requirement based on the measure is achieved. In our experiments below, we simply use a sufficiently large $T$ (see below), which works well in practice. It is interesting to explore an efficient early stopping rule in the future work.

How to design the network structure? Although Theorem 3.1 and Theorem 3.2 guarantee the minimax rates of TV-GAN without hidden layer and JS-GAN with one hidden layer, one may wonder whether deeper network structures will perform better. From our experiments, TV-GAN with one hidden layer is better than TV-GAN without any hidden layer. Moreover, JS-GAN with deep network structures can significantly improve over shallow networks especially when the dimension is large (e.g. $p\geq 200$). For a network with one hidden layer, the choice of width may depend on the sample size. If we only have 5,000 samples of 100 dimensions, two hidden units performs better than five hidden units, which performs better than twenty hidden units. If we have 50,000 samples, networks with twenty hidden units perform the best.

How to stabilize and accelerate TV-GAN? As we have discussed in Section 3.1, TV-GAN has a bad landscape when $N(\theta,I_{p})$ and the contamination distribution $Q$ are linearly separable (see Figure 1). An outlier removal step before training TV-GAN may be helpful. Besides, spectral normalization is also worth trying since it can prevent the weight from going to infinity and thus can increase the chance to escape from bad saddle points. To accelerate the optimization of TV-GAN, in all the numerical experiments below, we adopt a regularized version of TV-GAN inspired by Proposition 3.1. Since a good feature extractor should match nonlinear moments of $P=(1-\epsilon)N(\theta,I_{p})+\epsilon Q$ and $N(\eta,I_{p})$, we use an additional regularization term that can accelerate training and sometimes even leads to better performances. Specifically, let $D(x)={\sf sigmoid}(w^{T}\Phi(x))$ be the discriminator network with $w$ being the weights of the output layer and $\Phi_{D}(x)$ be the corresponding network after removing the output layer from $D(x)$. The quantity $\Phi_{D}(x)$ is usually viewed as a feature extractor, which naturally leads to the following regularization term , defined as

2 Numerical Supports for the Minimax Rates

In this section, we verify the minimax rates achieved by TV-GAN (Theorem 3.1) and JS-GAN (Theorem 3.2) via numerical experiments. The TV-GAN has no hidden layer, while the JS-GAN has one hidden layer with five hidden units in our experiments. All activation functions are sigmoid. Two main scenarios we consider here are $\sqrt{p/n}<\epsilon$ and $\sqrt{p/n}>\epsilon$, where in both cases, various types of contamination distributions $Q$, are considered.

We introduce the contamination distributions $Q$ used in the experiments. We first consider $Q=N(\mu,I_{p})$ with $\mu$ ranges in $\{0.2,0.5,1,5\}$. Note that the total variation distance between $N(0_{p},I_{p})$ and $N(\mu,I_{p})$ is of order $\|0_{p}-\mu\|=\|\mu\|$. We hope to use different levels of $\|\mu\|$ to test the algorithm and verify the error rate in the worst case. Second, we consider $Q=N(1.5*1_{p},\Sigma)$ to be a Gaussian distribution with a non-trivial covariance matrix $\Sigma$. The covariance matrix is generated according to the following steps. First generate a sparse precision matrix $\Gamma=(\gamma_{ij})$ with each entry $\gamma_{ij}=z_{ij}*\tau_{ij},i\leq j$, where $z_{ij}$ and $\tau_{ij}$ are independently generated from Uniform$(0.4,0.8)$ and Bernoulli$(0.1)$. We then define $\gamma_{ij}=\gamma_{ji}$ for all $i>j$ and $\bar{\Gamma}=\Gamma+(|\min\textnormal{eig}(\Gamma)|+0.05)I_{p}$ to make the precision matrix symmetric and positive definite, where $\min\textnormal{eig}(\Gamma)$ is the smallest eigenvalue of $\Gamma$. The covariance matrix is $\Sigma=\bar{\Gamma}^{-1}$. Finally, we consider $Q$ to be a Cauchy distribution with independent component, and the $j$th component takes a standard Cauchy distribution with location parameter $\tau_{j}=0.5$.

3 Comparisons with Other Methods

We perform additional experiments to compare with other methods including dimension halving and iterative filtering under various settings.

Dimension Halving. Experiments conducted are based on the code from https://github.com/kal2000/AgnosticMeanAndCovarianceCode. The only hyper-parameter is the threshold in the outlier removal step, and we take $C=2$ as suggested in the file outRemSperical.m.

Iterative Filtering. Experiments conducted are based on the code from https://github.com/hoonose/robust-filter. We assume $\epsilon$ is known and take other hyper-parameters as suggested in the file filterGaussianMean.m.

We emphasize that our method does not require any prior knowledge the nuisance parameters such as the contamination proportion $\epsilon$. Tuning GAN is only a matter of optimization and one can tune parameters based on the objective function only.

Table 4 shows the performances of JS-GAN, TV-GAN, dimension halving, and iterative filtering with i.i.d. observations sampled from $(1-\epsilon)N(0_{p},I_{p})+\epsilon Q$. The network structure, for both JS-GAN and TV-GAN, has one hidden layer with 20 hidden units when the sample size is 50,000 and 2 hidden units when sample size is 5,000. With fixed network structure, the hyper parameters are robust to various sampling distributions. For the network with 20 hidden units, the critical parameters to reproduce the results in the table are $\gamma_{g}=0.02$, $\gamma_{d}=0.2$, $K=5$, $T=150$ ($p=100$), $T=250$ ($p=200$), $T_{0}=25$ for JS-GAN and $\gamma_{g}=0.0001$, $\gamma_{d}=0.3$, $K=2$, $T=150$ ($p=100$), $T=250$ ($p=200$), $T_{0}=1$, $\lambda=0.1$ for TV-GAN, where $\lambda$ is the penalty factor of the additional regularization term (21). For the network with 2 hidden units, the critical parameters to reproduce the results below are $\gamma_{g}=0.01$, $\gamma_{d}=0.2$, $K=5$, $T=150$ ($p=100$), $T_{0}=25$ for JS-GAN and $\gamma_{g}=0.01$, $\gamma_{d}=0.1$, $K=5$, $T=150$ ($p=100$), $T_{0}=1$ for TV-GAN. We use Xavier initialization for both JS-GAN and TV-GAN trainings.

To summarize, our method outperforms other algorithms in most cases. TV-GAN is good at cases when $Q$ and $N(0_{p},I_{p})$ are non-separable but fails when $Q$ is far away from $N(0_{p},I_{p})$ due to optimization issues discussed in Section 3.1 (Figure 1). On the other hand, JS-GAN stably achieves the lowest error in separable cases and also shows competitive performances for non-separable ones.

4 Network Structures

In this section, we study the performances of TV-GAN and JS-GAN with various structures of neural networks. The experiments are conducted with i.i.d. observations drawn from $(1-\epsilon)N(0_{p},I_{p})+\epsilon N(0.5*1_{p},I_{p})$ with $\epsilon=0.2$. Table 5 summarizes results for $p=100$, $n\in\{5000,50000\}$ and various network structures. We observe that TV-GAN that uses neural nets with one hidden layer improves over the performance of that without any hidden layer. This indicates that the landscape of TV-GAN is improved by a more complicated network structure. However, adding one more layer does not improve the results. For JS-GAN, we omit the results without hidden layer because of its lack of robustness (Proposition 3.1). Deeper networks sometimes improve over shallow networks, but this is not always true. Table 6 illustrates the improvements of network with more than one hidden layers over that with only one hidden layer for JS-GAN when $p\in\{200,400\}$. We also observe that the optimal choice of the width of the hidden layer depends on the sample size.

5 Adaptation to Unknown Covariance

The robust mean estimator constructed through JS-GAN can be easily made adaptive to unknown covariance structure, which is a special case of (20). We define

The estimator $\widehat{\theta}$, as a result, is rate-optimal even when the true covariance matrix is not necessarily identity and is unknown (see Theorem 5.1). Below, we demonstrate some numerical evidence of the optimality of $\widehat{\theta}$ as well as the error of $\widehat{\Sigma}$ in Table 7.

6 Adaptation to Elliptical Distributions

To illustrate the performance of (20), we conduct a numerical experiment for the estimation of the location parameter $\theta$ with i.i.d. observations $X_{1},...,X_{n}\sim(1-\epsilon)\textnormal{Cauchy}(\theta,I_{p})+\epsilon Q$. The density function of $\textnormal{Cauchy}(\theta,I_{p})$ is given by $p_{\theta}(x)\propto\left(1+\|x-\theta\|\right)^{-(1+p)/2}$.

Table 8 shows the comparison with other methods. Our method still works well under Cauchy distribution, while the performance of other methods that rely on moment conditions deteriorates in this setting.

Discussions

In this paper, we study robust estimation via the technique of generative adversarial nets. We show that the presence of hidden layers are crucial for the estimators trained by JS-GAN to be robust. To better understand the intuition of the results in the paper, we give some further discussion from the perspective of variational lower bounds. In view of (4), we have

for any discriminator class $\mathcal{D}$. Moreover, according to , the optimal discriminator is achieved at

Interestingly, (23) is in the form of logistic regression, and this immediately implies that the variational lower bound (22) is sharp when we take $\mathcal{D}$ to be the class of logistic regression defined in (17). Indeed, when there is no contamination or $\epsilon=0$, the sample version of JS-GAN (18) with the logistic regression discriminator class (17) leads to the estimator $\widehat{\theta}=\frac{1}{n}\sum_{i=1}^{n}X_{i}$ according to Proposition 3.1, and this is obviously a minimax optimal estimator .

In contrast, when there is contamination or $\epsilon>0$, the logistic regression discriminator class (17) does not even lead to a consistent estimator. This is because the population objective function to be minimized is

instead of ${\sf JS}(N(\theta,I_{p}),N(\eta,I_{p}))$. The variational lower bound with the logistic regression discriminator class (17) is not sharp anymore because of the presence of the contamination distribution $Q$. In fact, a discriminator class $\mathcal{D}$ that leads to a sharp variational lower bound has to include the function

However, since there is no assumption on the contamination distribution $Q$, the discriminator function (24) can take an infinite many of forms. As a consequence, a discriminator class $\mathcal{D}$ that includes all possible functions in the form of (24) will certainly overfit the data, and thus is not practical at all. On the other hand, we show that for the purpose of robust mean estimation, we only need to add an extra hidden layer to the logistic regression discriminator class (17). The class (19) of neural nets with one hidden layer does not lead to a sharp variational lower bound, but it is rich enough for the estimator trained by JS-GAN to be robust against any contamination distribution. Moreover, the complexity of the class (19) is well controlled so that overfitting does not happen and thus the estimator achieves the minimax rate of the problem.

Future Projects.

Besides the topic of robust mean estimation, other important problems include robust covariance matrix estimation, robust high-dimensional regression, robust learning of Gaussian mixture models, and robust classification. It will be interesting to investigate what class of discriminators are suitable for these tasks. Another line of research is motivated from the goal to understand the class of divergence functions that are suitable for robust estimation. In addition to JS-GAN and TV-GAN studied in this paper, we would like to know whether it is possible to train robust estimators using GAN derived from other $f$-divergence functions. A further question is whether it is possible to use GAN derived from integral probability metrics including Wasserstein distance and maximum mean discrepancy . Finally, the landscapes and optimization properties of various GANs under robust estimation settings are topics to be explored.

Proofs

In this section, we present proofs of all technical results in the paper. We first establish some useful lemmas in Section 8.1, and the the proofs of main theorems will be given in Section 8.2.

with probability at least $1-\delta$ for some universal constant $C>0$.

with probability at least $1-\delta$. Using a standard symmetrization technique , we obtain the following bound that involves Rademacher complexity,

where $\epsilon_{1},...,\epsilon_{n}$ are independent Rademacher random variables. The Rademacher complexity can be bounded by Dudley’s integral entropy bound, which gives

with probability at least $1-\delta$ for some universal constant $C>0$.

Therefore, by McDiarmid’s inequality , we have

which uses Theorem 12 of . By Hölder’s inequality, we further have

Given i.i.d. observations $X_{1},..,X_{n}\sim N(\theta,I_{p})$ and the function class ${\mathcal{F}}_{L}^{H}(\kappa,\tau,B)$. Assume $\|\theta\|_{\infty}\leq\sqrt{\log p}$ and set $\tau=\sqrt{p\log p}$. We have for any $\delta>0$,

with probability at least $1-\delta$ for some universal constants $C>0$.

This leads to the desired result under the conditions on $\tau$ and $\|\theta\|_{\infty}$. ∎

2 Proofs of Main Theorems

We first introduce some notations. Define $F(P,\eta)=\sup_{w,b}F_{w,b}(P,\eta)$, where

With probability at least $1-\delta$, the above inequalities hold. We will explain each inequality. Since

which implies (27) and (31). The inequalities (28) and (30) are implied by Lemma 8.1 and the fact that

The inequality (29) is a direct consequence of the definition of $\widehat{\theta}$. Finally, it is easy to see that $F(P_{\theta},\theta)=0$, which gives (32). In summary, we have derived that with probability at least $1-\delta$,

where $f(t)=\int\frac{1}{1+e^{z+t}}\phi(z)dz$, with $\phi(\cdot)$ being the probability density function of $N(0,1)$. It is not hard to see that as long as $|f(t)-f(0)|\leq c$ for some sufficiently small constant $c>0$, then $|f(t)-f(0)|\geq c^{\prime}|t|$ for some constant $c^{\prime}>0$. This implies

with probability at least $1-\delta$. The proof is complete. ∎

We continue to use $P_{\eta}$ to denote $N(\eta,I_{p})$. Define

with $D(x)={\sf sigmoid}\left(\sum_{j\geq 1}w_{j}\sigma(u_{j}^{T}x+b_{j})\right)$. Then,

The inequalities (33)-(8.2) follow similar arguments for (27)-(31). To be specific, (34) and (36) are implied by Lemma 8.2, and (35) is a direct consequence of the definition of $\widehat{\theta}$. To see (33) and (8.2), note that for any $w$ such that $\|w\|_{1}\leq\kappa$, we have

A similar argument gives the same bound for $|\log(2(1-D(X)))|$. This leads to

which further implies (33) and (8.2). To summarize, we have derived that with probability at least $1-\delta$,

for all $\|w\|_{1}\leq\kappa$, $\|u_{j}\|\leq 1$ and $b_{j}$. Take $w_{1}=\kappa$, $w_{j}=0$ for all $j>1$, $u_{1}=u$ for some unit vector $u$ and $b_{1}=-u^{T}\theta$, and we get

with $Z\sim N(0,1)$. Direct calculations give

we have $\kappa f_{\delta}^{\prime}(0)\leq f_{\delta}(\kappa)+\kappa^{2}/4$. In view of (38), we have

It is easy to see that for the choices of $\sigma(\cdot)$, $\int\sigma(z)\phi(z)dz-\int\sigma(z+t)\phi(z)dz$ is locally linear with respect to $t$. This implies that

Therefore, with a $\kappa\lesssim\sqrt{\frac{p}{n}}+\epsilon$, the proof is complete. ∎

We continue to use $P_{\eta}$ to denote $N(\eta,I_{p})$. Define

Follow the same argument in the proof of Theorem 3.2, use Lemma 8.3, and we have

Therefore, $\max(x_{h},0),\max(-x_{h},0)\in\mathcal{G}_{l+1}^{H}(B)$ as long as $B\geq 2$. Hence, the above construction satisfies $D(x)={\sf sigmoid}(\kappa{\sf sigmoid}(\widetilde{u}^{T}(x-\theta)))\in\mathcal{F}_{L}^{H}(\kappa,\tau,B)$, and we have

where the definition of $f_{\delta}(t)$ is given by (39) with $Z\sim N(0,1)$ and $\sigma(\cdot)$ is taken as ${\sf sigmoid}(\cdot)$. Apply the a similar in the proof of Theorem 3.2, we obtain the desired result. ∎

We use $P_{\theta,\Sigma,h}$ to denote the elliptical distribution $EC(\theta,\Sigma,h)$. Define

with $D(x)={\sf sigmoid}\left(\sum_{j\geq 1}w_{j}\sigma(u_{j}^{T}x+b_{j})\right)$. The same argument in Theorem 3.2 leads to the fact that with probability at least $1-\delta$,

for all $\|w\|_{1}\leq\kappa$, $\|u_{j}\|\leq 1$ and $b_{j}$. Take $w_{1}=\kappa$, $w_{j}=0$ for all $j>1$, $u_{1}=u/\sqrt{u^{T}\widehat{\Sigma}u}$ for some unit vector $u$ and $b_{1}=-u^{T}\theta/\sqrt{u^{T}\widehat{\Sigma}u}$, and we get

where $\delta=\frac{u^{T}(\widehat{\theta}-\theta)}{\sqrt{u^{T}\widehat{\Sigma}u}}$ and $\Delta=\frac{\sqrt{u^{T}\Sigma u}}{\sqrt{u^{T}\widehat{\Sigma}u}}$. A similar argument to the proof of Theorem 3.2 gives

where $H(\delta)=\int\sigma(\delta+s)\widehat{h}(s)ds$. The above bound also holds for $\frac{\kappa}{2}(H(\delta)-H(0))$ by a symmetric argument, and therefore the same bound holds for $\frac{\kappa}{2}|H(\delta)-H(0)|$. Since $H^{\prime}(0)=\int\sigma(s)(1-\sigma(s))\widehat{h}(s)ds=1$, $H(\delta)$ is locally linear at $\delta=0$, which leads to a desired bound for $\delta=\frac{u^{T}(\widehat{\theta}-\theta)}{\sqrt{u^{T}\widehat{\Sigma}u}}$. Finally, since $u^{T}\widehat{\Sigma}u\leq M$, we get the bound for $u^{T}(\widehat{\theta}-\theta)$. The proof is complete by taking supreme of $u$ over the class of all unit vectors. ∎

Acknowledgement

The research of Chao Gao was supported in part by NSF grant DMS-1712957 and NSF Career Award DMS-1847590. The research of Yuan Yao was supported in part by Hong Kong Research Grant Council (HKRGC) grant 16303817, National Basic Research Program of China (No. 2015CB85600), National Natural Science Foundation of China (No. 61370004, 11421110001), as well as awards from Tencent AI Lab, Si Family Foundation, Baidu Big Data Institute, and Microsoft Research-Asia.

References