text:
\begin{document}

\maketitle

\begin{abstract}
Many modern datasets are collected automatically and are thus easily contaminated by outliers. This led to a regain of interest in robust estimation, including new notions of robustness such as robustness to adversarial contamination of the data. However, most robust estimation methods are designed for a specific model. Notably, many methods were proposed recently to obtain robust estimators in linear models (or generalized linear models), and a few were developed for very specific settings, for example beta regression or sample selection models.
In this paper we develop a new approach for robust estimation in arbitrary regression models, based on Maximum Mean Discrepancy minimization. We build two estimators. Both are proven to be robust to Huber-type contamination. One of them satisfies a strong non-asymptotic error bound and is also robust to adversarial contamination, at the cost of being more computationally expensive.
\end{abstract}

\section{Introduction}
\label{section:intro}

Robustness  is a fundamental problem in statistics, which aims at using statistical procedures that remain stable in presence of outliers. Historically, outliers were mistakes in data collection, or observations of individuals belonging to a different population than the population of interest. Robustness became even more important in the modern context where automatically collected datasets are often heterogeneous. Moreover, some strategic datasets are susceptible of malevolent manipulations.

In the statistical literature, the development of robust estimation methods for regression models generally focusses on the construction of $Z$-estimators   \cite[Chapter 5]{Vaart2000} for which each individual  observation can only have a bounded impact on the  estimating equations, and which therefore have a bounded influence function \cite{hampelrobust,hampelrobust2}. This   strategy has been successfully applied for the  robust estimation of the regression coefficients in generalized linear models    \cite{Kunsch1989,cantoni2001,cantoni2006}  as well as e.g.\ robust inference in the negative binomial regression model with unknown overdispersion parameter \cite{Aeberhard2014} or in the Heckman sample selection model \cite{Zhelonkin2016}. Based on other approaches, robust estimators for    mixtures of linear regression models  (see e.g.~\cite{bai2012}), for the Beta regression model with unknown precision parameter  \cite{Ghosh2019} or for  robust linear least square regression \cite{audibert2011robust} have been developed.

In the past ten years there was a renewed interest for robust methods in the machine learning community. Catoni developed a loss function whose minimization leads to robust estimators of the expectation of a random variable~\cite{catoni2012challenging}, a technique which was then adapted to many situations including linear regression~\cite{catoni2017dimension}. More generally, Lipschitz loss functions such as the absolute loss, or Huber's loss~\cite{huber1964robust}, lead to robustness of the empirical risk minimization procedure, a fact that was used in~\cite{chinot2018statistical,alquier2019estimation,chinot2019robust,holland2019pac} to study robust procedures of classification and regression. The Median-of-Means  approach of~\cite{nemi1983problem,devroye2016sub} was also adapted to regression~\cite{lugosi2016risk,lecueSPA,lugosi2019regularization,depersin} and classification~\cite{lecueML}. Minimax rates for regression in terms of the sample size and the contamination rate were derived in~\cite{diako2018,dalalyan2019outlier}.

In the discussion by Sture Holm in \cite{bickel1976another}, as well as in \cite{parr1980minimum}, minimum distance estimation is identified as a way to obtain robust estimators.  Building on this idea, \cite{Basu1998} introduced a density power divergence minimization  approach for robust inference in parametric models for i.i.d.\ observations. This procedure  is extended to  regression models in \cite{Ghosh2013} but suffers from two limitations. Firstly, the optimization of the objective function is, in general, a computationally challenging problem. Secondly, there is no general result which guarantees that the resulting $M$-estimator is robust. Its influence function is however known to be bounded   for the Gaussian linear regression model  \cite{Ghosh2013},  for the Poisson  and Logistic regression models \cite{Ghosh2016} and  for the Beta regression model with unknown precision parameter \cite{Ghosh2013}.

In this paper we introduce a new minimum distance estimation strategy for parameter inference in regression models which (a) is proven to be robust to outliers under general conditions on the statistical model and (b) only requires to be able to sample from the model and to compute the gradient of its log-likelihood function to be  applicable. In this sense, the approach proposed in this work defines a \textit{universal robust regression method}.

More specifically, following an idea recently introduced in \cite{barp2019minimum}, the    minimum distance estimation procedure presented in this paper relies on the Maximum Mean Discrepancy (MMD) distance. The MMD metric already turned out to be very useful in statistics and machine learning problems, e.g.\ for two sample test~\cite{smola2007hilbert,gretton2012kernel}, change-point detection~\cite{arlot2012kernel}, goodness-of-fit tests~\cite{jitkrittum2017linear}, or training GANs~\cite{li2015generative,dziugaite2015training}. The robustness of minimum MMD estimation with bounded kernels when the distribution of the data is completely specified is studied in~\cite{barp2019minimum,cherief2019finite}, see also~\cite{cherief2019mmd} for a Bayesian type estimator. (Note that unbounded kernels are used in~\cite{lerasle2019monk}, but the ``automatic'' robustness induced by bounded kernels is then lost, and the authors have to use a Median-of-Mean procedure to robustify their MMD estimator.) The main novelty in this paper is to extend the MMD method of~\cite{barp2019minimum,cherief2019finite} to the regression setting, where we only want to estimate the distribution $Y|X$, and not the distribution of the pair $(Y,X)$. This step turns out to be non-trivial, especially in the random design case.

 We prove that the proposed estimator $\hat{\theta}_n$ of the model parameter is (i) universally consistent, in the sense that it will always converge to the best approximation of the truth in the model, without \textit{any} assumption on the distribution generating the observations, and (ii) robust to adversarial contaminations. However, the computation time of $\hat{\theta}_n$ is quadratic with respect to the sample size $n$ and, for this reason, we introduce an alternative estimator $\tilde{\theta}_n$  which, as argued below, is expected to have a similar behaviour as  $\hat{\theta}_n$ in practice. The computation time of $\tilde{\theta}_n$ is linear in $n$ and we establish   that this estimator  is itself   robust to outliers,  in the sense that its  influence function is bounded.

The rest of this paper is organized as follows. In Section~\ref{section:notations}  we define the two estimators $\hat{\theta}_n$ and $\tilde{\theta}_n$, and discuss algorithms to compute them. In Section~\ref{section:theory}  we provide their theoretical analysis, both in the deterministic and random design case for $\hat{\theta}_n$ and in the random design case only for $\tilde{\theta}_n$. In  the random design case, our theoretical result for $\hat{\theta}_n$ requires assumptions on the reproducing kernel  that are quite involved. Section~\ref{section:kernels} provides examples of kernel satisfying these assumptions while Section~\ref{section:simulation} is devoted to the simulation study.  Section \ref{section:conclusion} concludes, and all the proofs are gathered in the Appendix \ref{app:B}.

\section{MMD-based regression}
\label{section:notations}

\subsection{ Set-up and  first assumptions}

Let $\mathcal{X}$ and  $\mathcal{Y}$ be two topological   spaces, equipped respectively with the Borel $\sigma$-algebra $\mathfrak{S}_\mathcal{X}$ and $\mathfrak{S}_\mathcal{Y}$, and let $\mathcal{Z}= \mathcal{X}\times \mathcal{Y}$, $\mathfrak{S}_\mathcal{Z}=\mathfrak{S}_\mathcal{X}\otimes \mathfrak{S}_\mathcal{Y}$ and $\mathcal{P}(\mathcal{Z})$ be the set of  probability distributions on $(\mathcal{Z},\mathfrak{S}_\mathcal{Z})$. All the random variables considered in this work are assumed to be defined on the same  probability space $(\Omega,\mathcal{F},\mathbb{P})$.

Let  $\{P_{\lambda},\lambda \in \Lambda\}$ be a set of probability distributions on $\mathcal{Y}$, $\Theta$ be a Polish space  and $g:\Theta\times \mathcal{X} \rightarrow \Lambda$ be such that the mapping $x\mapsto P_{g(\theta,x)}(A)$ is $\mathfrak{S}_\mathcal{X}$-measurable for all $A\in \mathfrak{S}_\mathcal{Y}$ and all $\theta\in\Theta$. Then, given an $\mathcal{Z}$-valued random variable $(Y,X)$, with $Y$ taking values in $\setY$,  we consider the statistical model $\{ ( P_{g(\theta,x)})_{x\in\setX},\,\theta\in\Theta\}$ for   the conditional distribution of $Y$ given $X$. For example, the Gaussian linear regression model with known variance is obtained by taking $P_\lambda = \mathcal{N}_1(\lambda,\sigma^2)$ and $g(\theta,x)=\theta^\top x$, the logistic regression model by taking $P_\lambda = \mathcal{B}er(\lambda)$  and $g(\theta,x)=1/(1+\exp(-\theta^\top x))$, and the Poisson regression model  by taking $P_\lambda=\mathcal{P}ois(\lambda)$ and $ g(\theta,x) = \exp(\theta^\top x)$. Other classical examples include binomial, exponential, gamma and inverse-Gaussian regression models.

Finally, we let  $D_n=\{(X_i,Y_i)\}_{i=1}^n$  be a   set of $n$ random variables   taking values in $\mathcal{Z}$. Throughout this paper we assume that a realization of $D_n$ is used to fit the regression model $\{ ( P_{g(\theta,x)})_{x\in\setX}$  and, for that reason, in what follows we will refer to $D_n$ as the observations.

\subsection{Background on RKHS and first assumptions}

We let $k:\mathcal{Z}^2 \rightarrow \mathbb{R}$ be a kernel on $\mathcal{Z}$ (i.e.\ $k$ is
symmetric and positive semi-definite) and  denote by $(\setH, <\cdot,\cdot>_{\setH})$   the reproducing kernel Hilbert space (RKHS) over $\mathcal{Z}$ having $k$ as reproducing kernel (see \cite{muandet2017kernel} for a comprehensive introduction to RKHS). We assume  throughout this work that the following condition on $k$ holds:
\begin{asm}
 \label{asm:k}
 The kernel $k$ is $\mathfrak{S}_\mathcal{Z}$-measurable and such that $|k|\leq 1$.
\end{asm}
Under Assumption \ref{asm:k},  for any probability distribution $P\in\mathcal{P}(\mathcal{Z})$ the quantity  $\mu(P) := \mathbb{E}_{Z\sim P}[k(Z,\cdot)] $ is well defined in $\mathcal{H}$. If in addition $k$ is such that the mapping $P\mapsto \mu(P)$ is one-to-one,   $k$ is said to be characteristic. In this case the MMD, that is the mapping  $\mathbb{D}_k:\mathcal{P}(\mathcal{Z})^2\rightarrow [0,2]$  defined by
$$
\mathbb{D}_k(P,Q) = \left\| \mu(P) - \mu(Q) \right\|_{\mathcal{H}}, \quad P,Q\in \mathcal{P}(\mathcal{Z})^2,
$$
is a metric on $\mathcal{P}(\mathcal{Z})$. We stress that none of the results presented in this work   requires $k$ to be characteristic. However, it is only under this assumption on $k$ that they provide useful and interpretable convergence guarantees for the proposed estimators.

Next,  we let $k_X$ be a kernel on $\mathcal{X}$, $k_Y$ be a kernel  on $\setY$ and we denote by $k_X\otimes k_Y$  the product kernel on $\mathcal{Z}$ such that $k_X\otimes k_Y((x,y),(x',y')) = k_X(x,x') k_Y (y,y')$ for all $(x,y),(x',y')\in \mathcal{Z}$.
 With this notation in place we can state a second assumption on $k$  that will be required for some of the results presented below.

\begin{asm}
 \label{asm:kxy}
 There is a continuous kernel $k_X$ on $\mathcal{X}$ and a kernel $k_Y$ on $\mathcal{Y}$ such that $|k_X|\leq 1$, $|k_Y|\leq 1$ and $k=k_X\otimes k_Y$.
\end{asm}

From Theorems 3-4 in~\cite{szabo2018characteristic}  it follows that, under Assumption \ref{asm:kxy}, $k$ is characteristic if  $k_X$ and $k_Y$   are continuous, translation invariant, characteristic and bounded. When $\setX=\R^d$ and $\setY=\R^{d_y}$ examples of such kernels $k_X$ and $k_Y$ include the Gaussian kernel, the exponential kernel and the Mat\'ern kernels.

\subsection{Definition of the estimators  \texorpdfstring{$\hat{\theta}_n$}{Lg} and \texorpdfstring{$\tilde{\theta}_n$}{Lg}}\label{sec:est}

Let   $\hat{P}^{n} = (1/n)\sum_{i=1}^{n} \delta_{(X_i,Y_i)}$ be the empirical distribution of the observations  and,   for every $\theta\in\Theta$,   let $\hat{P}_\theta^{n}$ be the (random) probability distribution on $\mathcal{Z}$ defined by
\begin{align}\label{eq:model}
\hat{P}_\theta^{n}\big(A\times B\big)=\frac{1}{n} \sum_{i=1}^{n} \delta_{X_i }(A)P_{g(\theta,X_i)}(B),\quad  A\times B\in \mathfrak{S}_\mathcal{X}\otimes \mathfrak{S}_\mathcal{X}.
\end{align}
In other words, if $(X,Y)\sim  \hat{P}_\theta^{n}$ then $X$ is uniformly distributed on the set $\{X_1,\dots,X_n\}$ and $Y|(X=x)\sim  P_{g(\theta,x)}$.

The first  estimator introduced in this work, $\hat{\theta}_n$, is  defined through the minimization of the MMD distance  between the probability  distributions $\hat{P}_\theta^n$ and $\hat{P}^{n}$, that is\footnote{When such a minimizer does not exist, we can use an $\epsilon$-minimizer instead and that follows can be trivially adapted.  In addition, we implicitly assume that $\hat{\theta}_n$ and $\tilde{\theta}_n$ are measurable, for all $n\geq 1$.}
\begin{equation}\label{dfn:estim}
\begin{split}
\hat{\theta}_n(D_n) &\in  \argmin_{\theta\in\Theta} \mathbb{D}_k^2\big( \hat{P}^n_\theta, \hat{P}^n)= \argmin_{\theta\in\Theta} \hat{F}_n(\theta),\quad   \hat{F}_n(\theta):=\sum_{ i,j=1}^n \hat{\ell}(\theta,X_i,X_j, Y_j)
\end{split}
\end{equation}
where,  for all $\theta\in\Theta$, $(x,x')\in\setX^2$ and $y\in\setY$,
\begin{equation}\label{dfn:estim_l}
\hat{\ell}(\theta, x,x',y)=\E_{  Y\sim P_{g(\theta,x)},  \,\,Y'\sim P_{g(\theta,x') }} \big[k\big((x,Y),(x',Y')\big) - 2  k\big((x,Y),(x',y)\big)\big].
\end{equation}

Remark that the number of terms in the definition of the function $\hat{F}_n$ to minimize is $\mathcal{O}(n^2)$. Although  below we propose an approach that makes possible to efficiently compute  $\hat{\theta}_n$ for moderate values of $n$ (i.e.\ for $n$ equals to a few thousands), this feature of   $\hat{F}_n$ limits the applicability of  $\hat{\theta}_n$  in large datasets. For large scale problems we propose the alternative  estimator $\tilde{\theta}_n$, defined by
\begin{equation}
\label{eq:estim:approx}
 \tilde{\theta}_n(D_n) \in  \argmin_{\theta\in\Theta}\tilde{F}_n(\theta),\quad   \tilde{F}_n(\theta):=\sum_{i=1}^n \tilde{\ell}(\theta,X_i,Y_i)
\end{equation}
where,  for all $\theta\in\Theta$, $x\in\setX$ and $y\in\setY$,
\begin{equation}\label{eq:estim:approx_l}
\tilde{\ell}(\theta, x,y)=\E_{\substack{ Y,Y'\,\,\iid P_{g(\theta,x)}}} \big[k_Y(Y,Y')  - 2 k_Y(Y,y)\big].
\end{equation}

The function $\tilde{F}_n$ defined in \eqref{eq:estim:approx} has the advantage to  involve  only $n$ terms but, on the other hand,     $\tilde{\theta}_n$  cannot be interpreted as the minimizer of a measure of discrepancy between  $\hat{P}_\theta^n$ and $\hat{P}^{n}$. For this reason,  our general theory for $\hat{\theta}_n$ does not apply to this estimator. Although   theoretical results regarding the robustness properties of  $\tilde{\theta}_n$ are provided in the next section,   they are weaker than those obtained for $\hat{\theta}_n$.

Henceforth we  use the shorthand  $\hat{\theta}_n = \hat{\theta}_n(D_n)$ and $\tilde{\theta}_n = \tilde{\theta}_n(D_n)$, which is standard in statistics.

\subsection{Link between the two estimators}\label{sub:link}

Assume that  $k=k_{\gamma_X}\otimes k_Y$, with $k_{\gamma_X}$ a strictly positive kernel on $\setX$ such  that $k_{\gamma_X}(x,x)=1$, and let   $\ell(\theta,x,x',y) =\hat{\ell}(\theta, x,x',y)/k_{\gamma_X}(x,x')$ if $k_{\gamma_X}(x,x')\neq 0$ and $\ell(\theta,x,x',y)=0$ otherwise. Remark that under the above assumption on $k$ the quantity $\ell(\theta,x,x',y)$ does not depend on   $k_{\gamma_X}$ and is such that $\ell(\theta,x,x,y)=\tilde{\ell}(\theta,x,y)$. Therefore,  letting
\begin{align}\label{eq:Hn}
h_n\big(\gamma_X,\theta,D_n\big)=2\sum_{i=1}^{n-1}\sum_{j=i+1}^n k_{\gamma_X}(X_i,X_j)\ell(\theta, X_i, X_j,Y_j),
\end{align}
it follows that the  estimators $\hat{\theta}_n$ and $\tilde{\theta}_n$ are such that
\begin{align*}
&\hat{\theta}_n\in\argmin_{\theta\in\Theta}\bigg\{ \sum_{i=1}^n \tilde{\ell}(\theta, X_i, Y_i) + h_n\big(\gamma_X,\theta,D_n\big)\bigg\},\quad \tilde{\theta}_n\in\argmin_{\theta\in\Theta} \sum_{i=1}^n\tilde{\ell}(\theta,  X_i,Y_i).
\end{align*}
Consequently, using $\tilde{\theta}_n$ in place of $\hat{\theta}_n$ amounts to discarding,  in the definition of this latter estimator, the  term $h_n(\gamma_X,\theta, D_n)$ whose computation  requires $\mathcal{O}(n^2)$ operations.

For standard choices of  $k_{\gamma_X}$ we have  $\lim_{\gamma_X\rightarrow 0}k_{\gamma_X}(x,x')=0$ for all $x'\neq x$. This property holds e.g.\ if $k_{\gamma_X}$  is the exponential or the Gaussian kernel with bandwidth parameter $\gamma_X>0$, or if $k_{\gamma_X}$  is the kernel defined in Section \ref{section:kernels} (with $\gamma_X>0$ as in that section). In this case, assuming that the $X_i$'s are $\P$-a.s.\ distinct, for all $\theta\in\Theta$ we have    $ \lim_{\gamma_X\rightarrow 0}h_n(\gamma_X,\theta, D_n)=0$, $\P$-a.s.  Consequently, under   suitable continuity assumptions, $\hat{\theta}_n \rightarrow \tilde{\theta}_n$ as $\gamma_X \rightarrow 0$, $\P$-a.s. For this reason,  and as illustrated in Section \ref{section:simulation}, for a small value of $\gamma_X$ we expect the two estimators $\hat{\theta}_n$ and $\tilde{\theta}_n$   to have a very similar behaviour in practice.

\subsection{Computation of the two estimators}
\label{subsec:algo}

Computing the estimator  $\hat{\theta}_n$ (resp.\ $\tilde{\theta}_n$) requires to optimize the function $\hat{F}_n$ (resp.\ $\tilde{F}_n$), which is defined through an expectation. In some models, the expectation   appearing in \eqref{dfn:estim_l} (resp.\ in  \eqref{eq:estim:approx_l}) can be computed explicitly, in which case the function  $\hat{F}_n$ (resp.\ $\tilde{F}_n$) has a known expression and standard optimization algorithms can be used to    optimize it. This is for example the case in logistic or multinomial regression, since for these two models the expectation in \eqref{dfn:estim_l} (and in  \eqref{eq:estim:approx_l}) is a finite sum.

In general, the functions  $\hat{F}_n$ and  $\tilde{F}_n$ are however intractable, and as a general strategy for computing $\hat{\theta}_n$ and $\tilde{\theta}_n$ we propose the use of a stochastic gradient   algorithm. As shown in the  Appendix \ref{app:A}, under suitable regularity conditions we have
\begin{align*}
&\nabla_\theta   \hat{\ell} (\theta, x,x',y)\\
 &=2 \E_{Y\sim P_{g(\theta,x)}, \,\,Y'\sim P_{g(\theta,x') }}   \bigg[\Big(k\big((x,Y),(x',Y')\big) -  k\big((x,Y),(x',y)\big)\Big) \nabla_{\theta}  \log p_{g(\theta,x)}(Y)  \bigg]
 \end{align*}
 and
 $$
\nabla_\theta \tilde{\ell}(\theta, x,y)=2\E_{\substack{ Y,Y'\,\,\iid P_{g(\theta,x)}}}\Big[  \Big(k_Y(Y,Y') - k_Y(Y,y) \Big) \nabla_{\theta}  \log p_{g(\theta,x)}(Y)   \Big].
$$
In the general case, we can easily compute an unbiased estimate of $\nabla \hat{F}_n(\theta)$ and of $\nabla \tilde{F}_n(\theta)$ under the mild conditions that (i) we can sample from  $P_\lambda$ for all $\lambda\in\Lambda$ and  (ii) that we can compute the gradient of the log-likelihood function of a single observation. The stability of this procedure obviously depends on the model, but it holds generally at least for compact parameter spaces $\Theta$, as discussed in Appendix \ref{app:A}

In both cases, it is important to stress that the functions $\hat{F}_n$ and $\tilde{F}_n$ being typically non-convex and multi-modal (in particular in the presence of outliers), the starting value of the stochastic gradient algorithm must be chosen with care. When the maximum likelihood estimator  of the model parameter can be efficiently computed we recommand to use its value to initialize the optimization procedure. Otherwise, the (gradient-free) algorithm introduced in \cite{gerber2020global}, designed to compute the global optimum of a function  defined through an expectation, can be used to find a good starting value for the stochastic gradient algorithm.

It is important to mention at this stage that the computation of $\hat{\theta}_n$ can be greatly facilitated by taking $k=k_{\gamma_X}\otimes k_Y$ for some small $\gamma_X>0$, with $k_{\gamma_X}$ as in Section \ref{sub:link}. Indeed, for such a kernel $k$  it often true that with high probability we have $k_{\gamma_X}(X_i,X_j)\approx 0$ for all $i\neq j$, and thus that $h_n\big(\gamma_X,\theta,D_n\big)\approx 0$ (with $h_n\big(\gamma_X,\theta,D_n\big)$ as defined in \eqref{eq:Hn}).  In this case, we can efficiently compute $\hat{\theta}_n$ with a stochastic gradient algorithm whose cost per iteration  is linear in the sample size $n$, as explained  in   Appendix \ref{subsec:SGD_hat}.

\section{Theoretical analysis }
\label{section:theory}

\subsection{Set-up and summary of the main results}

In this section we introduce theoretical results concerning the robustness of the   estimators $\hat{\theta}_n$ and $\tilde{\theta}_n$ when outliers are present in the dataset  used by the statistician   to fit the regression model. To this aim, we need to interpret $D_n$ as a contaminated version of a dataset $D^0_n:=\{(X^0_i,Y^0_i)\}_{i=1}^n$, and below   we consider two contamination models, namely the Huber contamination model \cite{huber1964robust} and the adversarial  contamination model, defined in Definition \ref{def:adv} and in Definition \ref{def:Huber}, respectively. The various possible contamination models are discussed in Section 2 in \cite{bateni2020confidence}.

\begin{definition}\label{def:adv}
We say that the observed dataset $D_n$ is an $\epsilon$-adversarial contamination of $D_n^0$ if there exists a set  $I\subset\{1,\dots,n\}$ such that  $|I|/n\leq \epsilon$ and such that $(X_i,Y_i)=(X^0_i,Y^0_i)$ for all $i\not\in I$.
\end{definition}

\begin{definition}\label{def:Huber}
We say that the observed dataset $D_n$ is an $(\epsilon,Q)$-Huber contamination of $D_n^0$ if there exists a $Q\in\mathcal{P}(\mathcal{Z})$ such that $(X_i,Y_i)\sim \epsilon Q+(1-\epsilon)\delta_{(X_i^0,Y_i^0)}$  for all $i\in\{1,\dots,n\}$.
\end{definition}

Below we   derive   non-asymptotic  bounds for  $\hat{\theta}_n$ under both  the fixed design and the  random design  scenarios, which prove the robustness of this estimator to adversarial contaminations of the data. For $\tilde{\theta}_n$, we derive an asymptotic result in the random design case which establishes the robustness of this estimator to Huber contaminations of the sample.

 We recall the reader that results in the fixed design case only provide guarantees on the estimation of the distribution of $Y$ when $X$ is equal to one of the observed $X_i$'s, while in practice regression is often used for out-of-sample predictions. In this case, and assuming that the pairs $(Y^0_i,X^0_i)$'s are i.i.d, this means that we want   guarantees on the estimation of the distribution of $Y$ when $X$ is drawn from the same (unknown) distribution as the observed $X_i$'s, and independently from them. This is precisely what  theoretical results in the random design case provide.

In what follows we let $(P^0_{Y|x})_{x\in\setX}$ be a regular conditional probability\footnote{We assume throughout  this work that $(P^0_{Y|x})_{x\in\mathcal{X}}$ exists, which is for instance the case if $\setX$ and $\setY$ are two Polish spaces.} of $Y$ given $X$, so $Y^0_i|(X^0_i=x)\sim P^0_{Y|x}$ for all $x\in\setX$ and all $i=1,\dots,n$.

\subsection{Non-asymptotic bounds for the estimator \texorpdfstring{$\hat{\theta}_n$}{Lg}}

\subsubsection{Fixed design case}

In the fixed design case   the $X^0_i$'s are deterministic and therefore it does not really make sense to assume that they can be manipulated by an adversary. For this reason, in this subsection we  assume that the contamination of the sample occurs only on the $Y^0_i$'s, so   that $X_i = X^0_i$ for all $i\in\{1,\dots,n\}$.

Then, letting
$$
\bar{P}^{0}_n(A\times B) = \frac{1}{n}\sum_{i=1}^n \delta_{X^0_i}(A) P_{Y|X^0_i}^0(B),\quad \forall (A\times B)\in \mathfrak{S}_\mathcal{Z},
$$
we set up our objective as the reconstruction   of   $\bar{P}^{0}_n\in\mathcal{P}(\mathcal{Z})$ by a distribution in the set $\{\hat{P}_\theta^n, \theta\in\Theta\}$.

The following lemma  gives  a non-asymptotic bound on the performances of the estimator  $\hat{\theta}_n$ for this task, under   an adversarial contamination of the sample.

\begin{lemma}\label{thm:fixed1}
Assume that $D_n$ is an $\epsilon$-adversarial contamination of $D_n^0$ such that $X_i=X_i^0$ for all $i\in\{1,\dots,n\}$. Then, under Assumption~\ref{asm:k},
$$
\mathbb{E} [  \mathbb{D}_k(\hat{P}^n_{\hat{\theta}_n} ,\bar{P}^{0}_n ) ]  \leq 4\epsilon + \inf_{\theta\in\Theta}  \mathbb{D}_k(\hat{P}^n_{\theta} ,\bar{P}^{0}_n ) + \frac{2}{\sqrt{n}}
$$
and, for all $\eta\in(0,1)$,
 \begin{align}
 \label{eq:tthm:fixed1:p}
 \mathbb{P}\bigg\{\mathbb{D}_k(\hat{P}^n_{\hat{\theta}_n} ,\bar{P}^{0}_n )  & < 4\epsilon + \inf_{\theta\in\Theta}  \mathbb{D}_k(\hat{P}^n_{\theta} ,\bar{P}^{0}_n) +  \frac{1}{\sqrt{n}}\big(2+\sqrt{ 2\log(1/\eta)}\big) \bigg\} \geq 1-\eta.
 \end{align}
\end{lemma}
\begin{remark}
Lemma \ref{thm:fixed1}   does not require \textit{any} assumption on the distribution of the data $\{(X_i^0,Y_i^0)\}_{i=1}^n$.
\end{remark}

In statistical theory  we often   assume that the ``truth is in the model'', that is, that there is a $\theta_0\in\Theta$ such that $\hat{P}^n_{\theta_0} = \bar{P}^0_n$. In this case, Lemma \ref{thm:fixed1} shows that
\begin{align}\label{eq:cont1}
\mathbb{E}\big[  \mathbb{D}_k(\hat{P}^n_{\hat{\theta}_n} ,\hat{P}^n_{\theta_0} )  \big]\leq 4\epsilon + 2/\sqrt{n}
\end{align}
while, in the non-contaminated case where $D_n = D_n^0$,
\begin{align}\label{eq:cont2}
\mathbb{E} [  \mathbb{D}_k(\hat{P}^n_{\hat{\theta}_n} ,\bar{P}^{0}_n ) ]  \leq  2/\sqrt{n}.
\end{align}
In  words, when  computed from the uncontaminated dataset $D_n^0$ the  estimator $\hat{\theta}_n$ is consistent for estimating $\theta_0$ (with respect to the MMD distance) and,   provided that  $\epsilon$ is small, an $\epsilon$-adversarial contamination of the sample will have only a negligible  impact on the estimated parameter value. Similar consequences can    be derived from the inequality in probability given in \eqref{eq:tthm:fixed1:p}.

Lemma~\ref{thm:fixed1}  implies  the convergence of $\hat{\theta}_n$ with respect to the MMD distance. However, under additional assumptions, it is possible to relate this form of convergence to the convergence with respect to a more usual criterion.

 \begin{thm}
\label{cor:fixed}
Under the assumptions of Lemma~\ref{thm:fixed1}, let $\|\cdot\|$  be a  semi-norm on $\Theta$ and $\theta_0\in\Theta$. Assume  that there exist  a neighbourhood $U$ of $\theta_0$ and a constant $\mu>0$ such that
$\mathbb{D}_k(\hat{P}^n_{\theta},\hat{P}^n_{\theta_0} )  \geq \mu\|\theta-\theta_0\|$ for all $\theta\in U$.
Let $\alpha = \inf_{\theta\in U^c}  \mathbb{D}_k(\hat{P}^n_{\theta}, \hat{P}^n_{\theta_0}) \in (0,2]$ and assume that $\epsilon \in[0, \alpha /32)$ and     that   $P^0_{Y|X^0_i}=P_{g(\theta_0, X^0_i)} $ for all $i\in \{1,\dots,n\}$ (well-specified model).
Then,
$$
\mathbb{P}\big(\lim\sup_{n\rightarrow\infty } \|\hat{\theta}_n-\theta_0\| \leq 4\epsilon/\mu\big)=1
$$
 and, for all $n\geq 64/\alpha^2$ and all $\eta \in [ 2e^{-n \alpha^2/38},1)$,
$$
\mathbb{P}\bigg\{  \|\hat{\theta}_n-\theta_0\| < \frac{4\epsilon}{\mu}+ \frac{1}{\mu \sqrt{n}} \big(2+\sqrt{2\log(2/\eta)}\big) \bigg\} \geq 1-\eta.
$$
 \end{thm}

\subsubsection{Random design case}\label{subsec:random_hat}

We   assume now that the pairs $(Y^0_i,X^0_i)$'s are i.i.d.\ from some probability distribution $P^0\in\mathcal{P}(\mathcal{Z})$. Let $P_X^0$ be the marginal distribution of the $X^0_i$'s and, for every $\theta\in\Theta$, let $P_\theta\in\mathcal{P}(\mathcal{Z})$ be defined by
\begin{align*}
P_\theta(A\times B)=\E_{X\sim P^0_X}\big[\ind_A(X) P_{g(\theta,X)}(B)],\quad  A\times B\in \mathfrak{S}_\mathcal{Z}.
\end{align*}
Then, we set up our objective as the reconstruction of $P^0$ by a distribution in $\{P_\theta,\,\theta\in\Theta\}$, and our natural candidate is $P_{\hat{\theta}_n}$. Since the approximating set $\{P_\theta,\,\theta\in\Theta\}$ is unknown (because it depends on $P_X^0$)   achieving this objective requires  more care  than in the fixed design setting. Notably, the results presented below are limited to the case where $k=k_X\otimes k_Y$ and impose some conditions on the RKHSs  $\setH_{X}$ and $\setH_Y$.
\begin{asm}
 \label{asm:X}
$\mathcal{X} \subseteq \mathbb{R}^{d}$ for some integer $d$. Moreover, $\mathcal{X}$ is path-wise connected and such that $\Lambda_d(\mathcal{X})>0$, with $\Lambda_d$  the Lebesgue measure on $\mathbb{R}^{d}$. (Remark that these latter two conditions impose no loss of generality as we can just replace $\mathcal{X} \subseteq \mathbb{R}^{d}$ by any set that contains $\mathcal{X}$ and that satisfies them.)
\end{asm}
\begin{asm}
 \label{asm:constants}
 Either the function $f=0$ is the only constant function in $\mathcal{H}_X$, either both $\mathcal{H}_X$ and $\mathcal{H}_Y$ contain non-zero constant functions.
\end{asm}
\begin{asm}
 \label{asm:cond:expect}
 For all $f\in\mathcal{H}_Y$ and $\theta\in\Theta$, $g\in\mathcal{H}_X$ where $g(x)=\mathbb{E}_{Y\sim P_{g(\theta,x)} }\left[f(Y)\right] $.
\end{asm}

Under the assumption that $k=k_X\otimes k_Y$, Assumption \ref{asm:constants}    holds when $k_X$ and $k_Y$ are such that $\inf_{(x,x')\in\setX^2}k_X(x,x')>0$ and $\inf_{(y,y')\in\setX^2}k_Y(y,y')>0$. Assumption \ref{asm:cond:expect} is harder to fulfil but we provide in  Section \ref{section:kernels} a class of characteristic kernels $k=k_X\otimes k_Y$ that verify this assumption (as well as Assumption  \ref{asm:constants}).

\begin{lemma}\label{thm:random1}
Assume that $D_n$ is an $\epsilon$-adversarial contamination of $D_n^0$. Then, under Assumptions~\ref{asm:k}-\ref{asm:cond:expect},  there exists a constant $C_{k_X}\in(0,\infty)$ that depends only on $k_X$  such that
$$
   \mathbb{E} [ \mathbb{D}_k  ( P_{\hat{\theta}_n} , P^0  )  ]  \leq 8\epsilon +
 \inf_{\theta\in\Theta}  \mathbb{D}_k  ( P_{\theta} , P^0  )
 +
 (\sqrt{2}C_{k_X} + 3)/\sqrt{n}
 $$
   and such that, for all $\eta\in (0,1)$,
 \begin{align*}
 \mathbb{P} & \bigg\{
  \mathbb{D}_k( P_{\hat{\theta}_n} , P^0 )
  <
  8\epsilon +
  \inf_{\theta\in\Theta}  \mathbb{D}_k( P_{\theta} , P^0)
  + \frac{(3+C_{k_X})(1+\sqrt{2\log(4/\eta) } }{\sqrt{n}}
 \bigg\} \geq 1-\eta.
 \end{align*}
\end{lemma}

From Lemma \ref{thm:random1} we can readily obtain the random design counterpart of the inequalities \eqref{eq:cont1}-\eqref{eq:cont2}, obtained in the fixed design setting, to prove the consistency of $\hat{\theta}_n$ in the well-specified case and in the absence of contamination.

The following theorem is the main result of this subsection.
 \begin{thm}
 \label{cor:random}

Assume that there is a unique $\theta_0\in\Theta$ such that $\mathbb{D}_k(P_{\theta_0},P^0 )=\inf_{\theta\in\Theta}\mathbb{D}_k(P_{\theta},P^0 )$ and  that there exist a neighbourhood $U$ of $\theta_0$ and a constant $\mu>0$ such that
\begin{align}\label{eq:loweDk}
\mathbb{D}_k(P_{\theta},P^0  ) \geq \mathbb{D}_k(P_{\theta_0},P^0  ) +\mu\|\theta-\theta_0\|,\quad\forall \theta\in U.
\end{align}
Let  $\alpha = \inf_{\theta\in U^c}  \mathbb{D}_k(P_{\theta},P^0)- \mathbb{D}_k(P_{\theta_0},P^0  ) \in (0,2]$ and  assume that $D_n$ is an $\epsilon$-adversarial contamination of $D_n^0$ for some $\epsilon\in[0, \alpha/64)$.  Assume also that Assumptions~\ref{asm:k}-\ref{asm:cond:expect} hold and let $C_{k_X}\in(0,\infty)$ be as in Lemma \ref{thm:random1}. Then,
$$
\mathbb{P}\big(\limsup_{n\rightarrow\infty } \|\hat{\theta}_n-\theta_0\| \leq 8 \epsilon/\mu\big)=1
$$
and there exist     constants $(C_1,C_2)\in(0,\infty)^2$, that depend only on $\alpha$ and on $C_{k_X}$, such that for all $ n\geq C_1$ and all  $\eta \in [ 8e^{- C_2 n},1)$ we have
$$
\mathbb{P}\bigg\{  \|\hat{\theta}_n-\theta_0\| < \frac{8\epsilon}{\mu}+ \frac{(C_{k_X}+3)}{\mu \sqrt{n}} \Big(1+\sqrt{2\log(8/\eta)}\Big)  \bigg\} \geq 1-\eta.
$$

\end{thm}

\begin{remark}
If $P^0=P_{\theta_0}$  for some $\theta_0\in\Theta$ (i.e.\ if the model is well-specified) then a sufficient condition for $\theta_0$ to be the unique global minimizer of the mapping $\theta\mapsto \mathbb{D}_k(P_{\theta}, P^0  )$ is that $k$ is characteristic and the model $\{P_\theta,\,\,\theta\in\Theta\}$ is identifiable, in the sense that $\theta_1\neq\theta_2\Rightarrow P_{\theta_1}\neq P_{\theta_2}$.
\end{remark}

In Theorem \ref{cor:random} the distribution $P_{\theta_0}$ should be interpreted as the best approximation of $P^0$ in the sense of the MMD distance $\mathbb{D}_k$ and thus, if the model contains  $P^0$, we have $P_{\theta_0}=P^0$ in the theorem. If a small fraction $\epsilon$ of the data is corrupted then the theorem ensures that $\hat{\theta}_n$  still estimates well $\theta_0$.  In particular, the first part of Theorem \ref{cor:random} implies that the influence function of the estimator $\hat{\theta}_n$ is bounded. It is also worth noting that since Theorem \ref{cor:random} imposes no conditions on $\Theta$, taking $\epsilon=0$ in this result establishes  the almost sure convergence of $\hat{\theta}_n$ towards $\theta_0$ for non-compact parameter spaces. This is an important difference with the popular maximum likelihood estimator  which, typically, requires   $\Theta$ to be  compact in order to converge with probability one to $\theta_0$.

Condition \eqref{eq:loweDk} of Theorem \ref{cor:random} requires that the function
$
\theta\mapsto \mathbb{D}_k(P_{\theta}, P^0   )- \mathbb{D}_k(P_{\theta_0}, P^0  )
$
is strongly convex in a neighbourhood $U$ of $\theta_0$, a condition which is rather weak, as shown in the next proposition.

\begin{prop}\label{prop:Taylor}
Assume that $\theta_0\in\argmin_{\theta\in\Theta}\mathbb{D}_k(P_{\theta},\tilde{P}^0 )$ is unique  and that the function   $\theta\mapsto \mathbb{D}_k(P_{\theta},\tilde{P}^0)$ is twice continuously differentiable at $\theta_0$. Then, condition \eqref{eq:loweDk} of Theorem \ref{cor:random} holds.
\end{prop}

\subsection{Asymptotic  guarantees for the estimator \texorpdfstring{$\tilde{\theta}_n$}{Lg}}

For this estimator we set up our objective as the reconstruction of the regular conditional probability $(P^0_{Y|x})_{x\in\mathcal{X}}$ by a distribution in the set $\{ (P_{g(\theta,x)})_{x\in\setX},\,\theta\in\Theta\}$. As in Section \ref{subsec:random_hat}, in what follows  we denote by $P^0$   the   distribution of the $(X_i^0, Y_i^0)$'s. In addition,  for $Q\in\mathcal{P}(\mathcal{Z})$ we let $Q_X$  and $(Q_{Y|x})_{x\in\setX}$ be the distribution of $X$ and the regular probability of $Y$ given $X$  when $(X,Y)\sim Q$, respectively.

It is direct to see that $\tilde{\theta}_n$ is an $M$-estimator,  and therefore   sufficient conditions on $k_Y$ and on the statistical model to ensure that $\tilde{\theta}_n$ converges to some value $\tilde{\theta}_{0}\in\Theta$ as $n\rightarrow\infty$ can be obtained from the general theory on $M$-estimators (\cite{Vaart2000}, Chapter 5). The focus of the present paper being on robust estimation, the following proposition only establishes this result under the strong assumption that $\Theta$ is a compact set.

\begin{prop}\label{prop:random2}
Let $k_Y$ satisfy $|k_Y|\leq 1$,   $D_n=D_n^0$  for all $n$, and assume  that the following conditions hold:
\begin{itemize}
\item   $\Theta$ is  compact,
\item    the mapping  $
\theta\mapsto \mathbb{E}_{Y, Y'\iid P_{g(\theta,x)} } \big[k_Y(Y,Y') -2 k_Y(Y,y)\big]
$ is continuous on $\Theta$  for all $(x,y)\in\setX\times\setY$,
\item    the mapping $\theta\mapsto \E\big[ \mathbb{D}_{k_Y}\big(P_{g(\theta,X_1)},  P^0_{Y|X_1}\big)^2\big]$ has a unique global minimum at $\tilde{\theta}_{0}\in\Theta$.
\end{itemize} Then, $\tilde{\theta}_n\rightarrow \tilde{\theta}_{0}$ in $\mathbb{P}$-probability.
\end{prop}

\begin{remark}
If the model is well-specified, that is if $(P^0_{Y|x})_{x\in\mathcal{X}}\in\{ (P_{g(\theta,x)})_{x\in\setX},\,\theta\in\Theta\}$, then a sufficient condition for $\tilde{\theta}_0$ to be well-defined is that $k_Y$ is a characteristic kernel and  the model $\{ (P_{g(\theta,x)})_{x\in\setX},\,\theta\in\Theta\}$ is identifiable, in the sense that $\mathbb{P}\big(\theta_1\neq\theta_2\Rightarrow P_{g(\theta_1,X)}\neq P_{g(\theta_2,X)}\big)=1$.
\end{remark}
\begin{remark}

Using the general theory on $M$-estimators (\cite{Vaart2000}, Chapter 5)  one can obtain sufficient conditions on $k_Y$ and on the statistical model which ensure that $\tilde{\theta}_n$ is $\sqrt{n}$-consistent and asymptotically Gaussian.
\end{remark}

\begin{remark} Unless the model is well-specified,    the parameter value $\tilde{\theta}_{0}$ defined in  Proposition \ref{prop:random2}  may be different from the parameter value $\theta_0$  the estimator $\hat{\theta}_n$ converges to (see Theorem \ref{cor:random}). \end{remark}

Recalling  that for $\tilde{\theta}_n$ we focus on Huber's type contaminations of the sample, in Proposition \ref{prop:random2} there is no loss of generality to assume that  $D_n=D_n^0$ for all $n$ (no contamination). Indeed, in the random design setting, under an  $(\epsilon,Q)$-Huber contamination of the sample we have that  $(X_i,Y_i)\iid \tilde{P}^0:=\epsilon Q+(1-\epsilon) P^0$ and thus, since   Proposition \ref{prop:random2} requires no assumption on   $P^0$, it results remains valid of $P^0$ is replaced by $\tilde{P}^0$.

The following theorem provides an asymptotic guarantee regarding  the robustness of  $\tilde{\theta}_n$  to Huber type contamination of the data.  Notably, a direct implication of this theorem is that the influence function   of $\tilde{\theta}_n$ is bounded.

\begin{thm}\label{thm:random2}
Let $k_Y$ satisfy $|k_Y|\leq 1$, $Q\in\mathcal{P}(\mathcal{Z})$ and assume that following two conditions hold:
\begin{itemize}
\item  the mapping  $\theta\mapsto \E\big[ \mathbb{D}_{k_Y}\big(P_{g(\theta,X^0_1)},  P^0_{Y|X^0_1}\big)^2\big]$ has a unique global minimum at $\tilde{\theta}_0$,

\item  there exist a neighbourhood $U$ of $\tilde{\theta}_0$ and a constant $\mu>0$ such that, for all $\theta\in U$,
\begin{align}\label{eq:cond_tildeU}
\E_{X \sim P_X^0}\big[ \mathbb{D}_{k_Y}(P_{g(\theta,X)}, P_{Y|X}^0)^2\big]\geq \E_{X \sim P_X^0}\big[ \mathbb{D}_{k_Y}(P_{g(\tilde{\theta}_0 ,X)}, P_{Y|X}^0)^2\big]+\mu\|\theta-\tilde{\theta}_0\|.
\end{align}
\end{itemize}
Let
$$
\alpha= \inf_{\theta\in U^c} \E_{X \sim P_X^0}[ \mathbb{D}_{k_Y}(P_{g(\theta,X)}, P_{Y|X}^0)^2]-\E_{X \sim P_X^0}[ \mathbb{D}_{k_Y}(P_{g(\tilde{\theta}_0,X)}, P_{Y|X}^0)^2]\in(0,4]
$$
and assume that $\epsilon\in [0, \alpha/(52+\alpha))$. Then,  for all
$$
\tilde{\theta}_{Q,\epsilon}\in \argmin_{\theta\in\Theta}\E_{X_1\sim\epsilon Q_X+(1-\epsilon) P^0_X}\big[ \mathbb{D}_{k_Y}\big(P_{g(\theta,X_1)}, \epsilon Q_{Y|X_1}+(1-\epsilon)P^0_{Y|X_1}\big)^2\big]
$$
 we have $\|\tilde{\theta}_{Q,\epsilon}-\tilde{\theta}_0\|\leq   52 \epsilon/(\mu -\epsilon\mu)$.
\end{thm}

Following similar steps as in the proof of Proposition \ref{prop:Taylor} it is easily checked  that condition \eqref{eq:cond_tildeU} holds when the function $\theta\mapsto \E_{X \sim P_X^0}[ \mathbb{D}_{k_Y}(P_{g(\theta,X)},P_{Y|X}^0)^2]$ has $\tilde{\theta}_0$ as unique global minimizer and is twice continuously differentiable around  this parameter value.

\section{A class of characteristic kernels satisfying the assumptions of Lemma \ref{thm:random1}}
\label{section:kernels}

Let $\psi:\R\rightarrow (0,1)$ be such that $\psi(0)=1/2$ and $\psi(v) = 1/2 + (\sqrt{v^2 +1}-2)/v$ for all $v\neq 0$. For any $d\in\mathbb{N}$ we let $ \psi_{(d)}(v) = (\psi(v_1),\dots,\psi(v_d)) $, noting that $\psi_{(d)}:\R^d\rightarrow(0,1)^d$ is a $C^1$-diffeomorphism. Then,  we let $k_{\alpha,\gamma_X}$ be the kernel on $\mathcal{X}$ defined by $k_{\alpha, \gamma_X}(x,x') = K_{\alpha,\gamma_X}(\|\psi_{(d)}(x)-\psi_{(d)}(x')\|)$, with $K_{\alpha,\gamma_X}$ the  Mat\'ern kernel with smoothness parameter $\alpha>0$ and bandwidth parameter $\gamma_X>0$. We refer the reader to  Example 2.2 in \cite{kanagawa2018gaussian} for the definition $K_{\alpha,\gamma_X}$ and note that  $K_{\alpha,\gamma_X}$ reduces to the exponential kernel when $\alpha=1/2$, i.e.\ $K_{1/2,\gamma_X}(\|x-x'\|)=\exp(-\|x-x'\|/\gamma_X)$.

The following proposition  shows that the class of kernels $\{k_{\alpha,\gamma_X},\,(\alpha,\gamma_X)\in(0,\infty)^2\}$  can be used to construct a characteristic kernel $k=k_X\otimes k_Y$ on $\setX\times\setY$ that verifies Assumptions \ref{asm:k}, \ref{asm:kxy} and \ref{asm:constants}.
\begin{prop}\label{prop:generalK}
Let $\setX=\R^d$,  $k_X=c \,k_{\frac{m}{2},\gamma_X}$ for some constants $(m,\gamma_X,c)\in\mathbb{N}\times(0,\infty)^2$ such that  $|k_X|\leq 1$, and let   $k_Y=\beta K+(1-\beta)$ for some $\beta\in(0,1)$ and for some continuous translation invariant characteristic kernel $K$ on $\setY$ satisfying $|K|\leq 1$. Then, the  kernel $k:=k_X\otimes k_Y$ on $\setX\times\setY$ is characteristic and Assumptions \ref{asm:k}, \ref{asm:kxy} and \ref{asm:constants} hold.
\end{prop}

The next result notably implies that, under some regularity conditions on the model $\{ (P_{g(\theta,x)})_{x\in\setX},\,\theta\in\Theta\}$, the kernel $k$   defined in Proposition \ref{prop:generalK}   satisfies  Assumption  \ref{asm:cond:expect}.

\begin{thm}\label{thm:kernel1}
Let $\setX=\R^d$, $k_X=c \,k_{\frac{m}{2},\gamma_X}$ for some constants $(m,\gamma_X,c)\in\mathbb{N}\times(0,\infty)^2$ such that  $|k_X|\leq 1$, let  $k_Y$ be a continuous   kernel  on $\setY$ satisfying $|k_y|\leq 1$, and let $s=(d+m)/2$ if $(d+m)$ is even and $s=(d+m+1)/2$ if $(d+m)$ is odd. Let  $A_s= \big\{a\in\mathbb{N}_0^{d}:\,\sum_{i=1}^{d}a_i\leq s\big\}$ and assume that  each distribution $P_\lambda$ has a density $p_\lambda$ with respect to a measure $\mu$.

For all  $(\theta,x,y,a)\in\Theta\times\setX\times\setY\times A_s$  let $I_a=\{\lambda\in A_s:\, \lambda_i\leq a_i,\,\forall i\in\{1,\dots,d\}\}$ and
$$
h_{a,\theta}(y,x)=\bigg(\max_{\lambda\in I_a}\prod_{i=1}^d\big(1+  x_i^2 \big)^{\lambda_i+1}\bigg)\,\max_{\lambda\in I_a}\bigg| \frac{\partial^{\sum_{i=1}^d \lambda_i}}{\partial x_1^{\lambda_1}\dots \partial x_d^{\lambda_d}}p_{g(\theta,x)}(y)\bigg|,
$$
 with the convention $h_{a,\theta}(y,x)=p_{g(\theta,x)}(y)$ when $a=(0,\dots,0)$.

 Assume that the following conditions hold for all $(\theta,a)\in\Theta\times A_s$:
\begin{enumerate}[label=(C\arabic*)]
\item\label{t1} $\int_{\setX} h_{a,\theta}(y,x)^2\Lambda_d(\dd x)<\infty$ for all $y\in\setY$,
\item\label{t2}  $\int_{\setY\times \setX} h_{a,\theta}(y,x)^2 \Lambda_d\otimes \mu(\dd(x ,y))<\infty$,
\item\label{t3} one of the following two conditions hold:
\begin{enumerate}
\item\label{t3a} The set $\setY$ is countable,
\item \label{t3b}

\begin{enumerate}
\item[(i)] For some $M\in\mathbb{N}$ there exist separable sets $\{\setY_m\}_{m=1}^M$ such that     $\cup_{m=1}^M \setY_m=\setY$ and such that,  for all $m\in\{1,\dots,M\}$, $\setY_m\in \mathfrak{S}_\setY$ and the function $y\mapsto \frac{\partial^{\sum_{i=1}^{d}a_i}}{\partial x_1^{a_1}\dots \partial x_{d}^{a_{d}}}p_{g(\theta,x)}(y)$ is  continuous on $\setY_m$,   for all $x\in\setX$.
\item[(ii)]\label{t3b3} $\sup_{x\in\setX} |h_{a,\theta}(y,x)|<\infty$ for all $y\in\setY$.
\end{enumerate}
\end{enumerate}

\end{enumerate}
Then, Assumption  \ref{asm:cond:expect} is  satisfied.
\end{thm}

\begin{rmk}
If the measurable space $(\setX, \mathfrak{S}_\mathcal{X},\Lambda_d)$ is complete then part (ii) of Condition \ref{t3b}) can be  omitted. We conjecture that under additional   technical assumptions on $(\setY, \mathfrak{S}_\mathcal{Y},\mu)$ Condition \ref{t1} may be omitted too.
\end{rmk}

\begin{rmk}
It is worth noting that the estimator $\hat{\theta}_n$ does not depend on the constants $c$ and $\beta$ appearing  in the definition of $k_X$ and $k_Y$. Thus, in practice we can always set $c=\beta=1$.
\end{rmk}

\begin{rmk}
 Other $C^1$-diffeomorphisms  $\psi:\R\rightarrow(0,1)$   than the one considered above can used in the definition of  $k_{\alpha,\gamma_X}$, such as the logistic function or the (rescaled) $\mathrm{arctan}$ function. However, other choices for the mapping $\psi$ may require stronger  assumptions on the model for the conclusion of Theorem \ref{thm:kernel1} to hold.
\end{rmk}

The following corollary summarizes the conclusions of  Proposition \ref{prop:generalK} and of Theorem \ref{thm:kernel1} for some popular regression models.
\begin{corollary}\label{cor:example}
Let $\setX=\R^d$, $k_X$ and $k_Y$ be as in Proposition \ref{prop:generalK} and as in Theorem \ref{thm:kernel1}, and let $k=k_X\otimes k_Y$. Then $k$ is characteristic and Assumptions \ref{asm:k}-\ref{asm:cond:expect} are satisfied for:
\begin{enumerate}
\item The Gaussian linear regression model with unknown variance (see Section \ref{subsec:expe:linear}),

\item For the mixture of Gaussian linear regression models with unknown variances and unknown mixture weights; that is for $P_{g(\theta,x)}=\sum_{m=1}^M\alpha_m \mathcal{N}_1(\beta_m^\top x,\sigma_m^2)$ with $M\in\mathbb{N}$ known and with $\theta=(\beta_1,\dots,\beta_M, \sigma_1,\dots,\sigma_M,\alpha_1,\dots,\alpha_M)\in\Theta$ where
$$
\Theta=\R^{Md}\times(0,\infty)^M\times\Big\{(a_1,\dots,a_M)\in [0,1]^M:\sum_{m=1}^M a_m=1\Big\}.
$$

\item For the Poisson linear regression model;  that is for $P_{g(\theta,x)}=\mathcal{P}ois(\exp(\theta^\top x))$ with $\theta\in\Theta:=\R^d$.

\item For the Logistic regression model; that is for $P_{g(\theta,x)}=\mathcal{B}er\big(1/(1+\exp(\theta^\top x))$ with   $\theta\in\Theta:=\R^d$.

\item For the Heckman sample selection model (see Section \ref{sub:Heck}).

\item For the Gamma regression model with unknown shape parameter (see Section \ref{sub:gamma}).
\end{enumerate}
\end{corollary}

\section{Numerical experiments}\label{section:simulation}

\subsection{Set-up\label{sub:sim_setup}}

All the results presented in this section are obtained with the kernel $k=k_Y \otimes k_X$, where $k_Y$ is the exponential kernel with parameter $\gamma_Y=1$, that is with $k_Y(y,y')=\exp(-\|y-y'\|)$,  and where  $k_X=c \,k_{\frac{m}{2},\gamma_X}$ with $c=m=1$ and $\gamma_X=0.01$.

For each  experiment   the  observations  used to fit the model  are obtained from an uncontaminated dataset $d^0_{N}=\{(x_i^0,y_i^0)\}_{i=1}^N$   that we  contaminate  as follows. We choose an $\epsilon\in [0,1)$ and  randomly select a set $I\subset\{1,\dots,N\}$ such that $|I|=\lfloor \epsilon N \rfloor$. For all $i\in\{1,\dots,N\}$  we then let $(x_i, y_i)=(x^0_i, y^0_i)$ if $i\not\in I$ and  $(x_i, y_i)=(x^c_i, y^c_i)$ for some $(x^c_i, y^c_i)\in\mathcal{Z}$ if $i\in I$. The way we generate the  $(x^c_i, y^c_i)$'s will vary from one example to the next and will therefore be specified in due course. Finally,  for $n\leq N$ we let $d_n:=\{(x_i,y_i)\}_{i=1}^n$ be the sample available   to estimate the model parameters.

Below the value of the estimators   $\hat{\theta}_n$ and $\tilde{\theta}_n$ are obtained using AdaGrad  \cite{Duchi}, an adaptive stochastic gradient algorithm, and the strategy exposed in Appendix \ref{app:A}  for computing the  former estimator is implemented. As  suggested in Section \ref{subsec:algo}, the  algorithms  that compute   $\hat{\theta}_n$ and $\tilde{\theta}_n$ use the maximum likelihood estimator  as starting value.

\subsection{Gaussian linear regression}\label{subsec:expe:linear}

We let $d=8$ and, for every $x\in\R^d$, we let $P_{g(\theta,x)}=\mathcal{N}_1(\beta^\top x,\sigma^2)$ with $\theta=(\beta,\sigma)\in\Theta:=\R^d\times(0,\infty)$. For this example the dataset $d_N^0$ is constructed by simulating  independent observations  using
$$
Y^0_i=\beta_0^\top X^0_i+\epsilon_i,\quad (X^0_i,\epsilon_i)\iid\mathcal{N}_d(0,I_d)\otimes \mathrm{Laplace}(0,\sigma_0)
$$
with $\beta_0=(4,4,3,3,2,2,1,1)$ and $\sigma_0=1$. The model $\big\{ \big(P_{g(\theta,x)}\big)_{x\in\R^d},\,\theta\in\Theta\big\}$ is therefore misspecified and, in what follows, we focus on the estimation of $\beta_0$.

We let $N=5\,000$ and for $n\leq N$ and $\epsilon>0$ the contaminated dataset $d_n$ is generated as explained in Section \ref{sub:sim_setup}, where two types of outliers $(x_i^c,y^c_i)$  are considered for this example. More precisely, we say that the outliers are of type   $\mathsf{X}$ when $y_i^c=y_i$ and $x_i^c$ is such that $x^c_{ij}=x_{ij}$ for all $j>1$ and such that $x_{i1}^c$ is a random draw from the from the $\mathcal{N}(5,1)$ distribution, and of type  $\mathsf{Y}$ when $y_i^c$  is a random draw from the $\mathcal{N}(10,1)$ distribution and $x^c_i=x_i$.

In Table \ref{table1} we report, for different values of $n\leq N$ and of $\epsilon\in[0,0.03]$, the root mean squared error (RMSE)  for the estimation of $\beta_0$ obtained for four estimators, namely the ordinary least square estimator $\beta_{\mathrm{ols},n}$, the least absolute deviation estimator $\beta_{\mathrm{lad},n}$, the robust estimator of $\beta$ in linear Gaussian regression models proposed by \cite{koller2011} (computed using the \texttt{R} package \texttt{robustbase}),  and the two proposed estimators $\hat{\beta}_n$ and $\tilde{\beta}_n$.

When the sample is not contaminated (i.e.\ when $d_n=d_n^0$) we observe that the OLS estimator is the best estimator. However, this estimator is extremely sensitive to the presence of outliers (of both types), a fact that is already well documented in the literature, see for example Chapter 1 in \cite{Rousseeuw}. The estimator $\beta_{\mathrm{lad},n}$ is robust to outliers  but its RMSE is almost always higher than the one obtained with the  two MMD based estimators $\hat{\beta}_n$ and $\tilde{\beta}_n$, especially in the case of outliers of type $\mathsf{X}$. We also observe that $\hat{\beta}_n$ and $\tilde{\beta}_n$ are robust to the two considered  types of outliers, as predicted by our theory, and have a performance which is almost identical to that of the robust estimator $\beta_{\mathrm{rob},n}$ proposed \cite{koller2011}, which is the one which tends to have the lowest RMSE in Table \ref{table1}.

\begin{table}[t]
\begin{center}
\begin{scriptsize}
  \begin{tabular}{|c|c|c|c|c|c|c|c|}
   \hline

   $\tau$ &type &$n$&  $\beta_{\mathrm{ols},n}$& $\beta_{\mathrm{lad},n}$&$\beta_{\mathrm{rob},n}$&  $\hat{\beta}_n$
          &  $\tilde{\beta}_n$
   \\

   \hline \hline
   \multirow{3}{*}{0} & & $100$   &0.372&0.353&0.350& 0.355& 0.334\\
    &&1\,000 & 0.116& 0.092&0.104&0.108&0.107  \\
    &&5\,000  &0.053&0.039&0.046&0.049&0.047  \\
   \hline

    \multirow{3}{*}{1}  & \multirow{3}{*}{$\mathsf{Y}$} &100 &0.464& 0.339&0.385&0.350&0.342\\
    &&1\,000 &0.181& 0.094&  0.106& 0.105& 0.097  \\
    &&5\,000  & 0.103& 0.043&  0.049&0.054&0.051 \\
   \hline
     \multirow{3}{*}{2}  & \multirow{3}{*}{$\mathsf{Y}$} &100&0.647&0.351&0.359&  0.337& 0.333 \\
    &&1\,000 &0.241&0.097&0.110&  0.114 &0.115  \\
    &&5\,000  & 0.175&0.039&  0.047& 0.051&0.052 \\
   \hline
       \multirow{3}{*}{3}  & \multirow{3}{*}{$\mathsf{Y}$} &100&0.724&0.331& 0.343&0.329&0.320 \\
    &&1\,000 & 0.309& 0.100& 0.108&0.113&0.110 \\
    &&5\,000  & 0.250&0.043&  0.048&0.053& 0.055  \\
    \hline

          \multirow{3}{*}{1}  & \multirow{3}{*}{$\mathsf{X}$} &100&0.870&0.356&0.374&0.342&0.338   \\
    &&1\,000 &0.836& 0.111& 0.105&  0.104&0.096    \\
    &&5\,000  & 0.818& 0.065&0.049&0.054&   0.052    \\
    \hline
          \multirow{3}{*}{2}  & \multirow{3}{*}{$\mathsf{X}$} &100& 1.575&0.400&0.347&0.337&0.331 \\
    &&1\,000 &1.467&0.160&0.110&  0.112& 0.115  \\
    &&5\,000  & 1.401& 0.119& 0.046&0.051& 0.052         \\
     \hline
          \multirow{3}{*}{3}  & \multirow{3}{*}{$\mathsf{X}$} &100&1.838&0.442&0.344&0.331& 0.323  \\
    &&1\,000 & 1.805&  0.216&0.108& 0.113& 0.109\\
    &&5\,000  & 1.771& 0.183&  0.048& 0.054&0.056  \\
    \hline
  \end{tabular}
  \hspace{0.1cm}

  \vspace{0.3cm}

 \caption{Results for the Gaussian linear regression model.  For each experimental setting  we report the RMSE over
25 replications.}
 \label{table1}
 \end{scriptsize}
 \end{center}
\end{table}

\begin{table}
\begin{center}
\begin{scriptsize}
  \begin{tabular}{|c|c|c|c|c|c|}
   \hline

   $\epsilon$ &$n$&  $\theta_{\mathrm{mle},n}$ & $\theta_{\mathrm{rob},n}$&  $\hat{\theta}_n$
          &  $\tilde{\theta}_n$
   \\

   \hline \hline
   \multirow{3}{*}{0\%} & $100$  &1.504&1.737 &2.132 &2.012  \\
   &  $ 1\,000$  &0.565&0.608 & 0.733&0.735  \\
   &  $5\,000$  &0.210& 0.296& 0.290&0.289  \\
    \hline
   \multirow{3}{*}{1\%} & $100$  & 1.767&2.020 &2.017&1.802  \\
   &  $ 1\,000$  &1.325&1.164 & 0.706&0.695  \\
   &  $5\,000$  & 1.293&1.083 &0.283 & 0.283  \\
    \hline
     \multirow{3}{*}{2\%} & $100$  &2.218&4.131 &2.238 &2.186   \\
   &  $ 1\,000$  &1.766&1.551 &0.683 &0.681      \\
   &  $5\,000$  & 1.936& 1.669&0.252 &0.257 \\
    \hline
     \multirow{3}{*}{3\%} & $100$  & 2.657&3.048 &1.786 &1.755   \\
   &  $ 1\,000$  &2.496&2.253 & 0.631& 0.642 \\
   &  $5\,000$  &2.404&2.142 & 0.243& 0.242   \\
    \hline

  \end{tabular}
  \hspace{0.5cm}\begin{tabular}{|c|c|c|c|c|c|}
   \hline

   $\epsilon$ &$n$& $\beta_{\mathrm{mle},n}$ & $\beta_{\mathrm{rob},n}$&  $\hat{\beta}_n$
          &  $\tilde{\beta}_n$
   \\

   \hline \hline
   \multirow{3}{*}{0\%} & $100$  &1.451&1.696&2.0638&1.944  \\
   &  $ 1\,000$  &0.536&0.574&0.705&0.708  \\
   &  $5\,000$  &0.202&0.238&0.273&0.269   \\
    \hline
   \multirow{3}{*}{1\%} & $100$  &1.598&1.731&1.947&1.727  \\
   &  $ 1\,000$  &1.009&0.584&0.675&0.669  \\
   &  $5\,000$  &0.878&0.273 &0.267&0.266  \\
    \hline
     \multirow{3}{*}{2\%} & $100$  &1.958&3.684 &2.189&2.147 \\
   &  $ 1\,000$  &1.312&0.638&0.654&0.651 \\
   &  $5\,000$  &1.389&0.341&0.239&0.241 \\
    \hline
     \multirow{3}{*}{3\%} & $100$  &2.200&1.989&1.732&1.697   \\
   &  $ 1\,000$  &1.861&0.707&0.604&0.611 \\
   &  $5\,000$  &1.762&0.418&0.231&0.230   \\
    \hline
  \end{tabular}

  \vspace{0.3cm}

 \end{scriptsize}
   \caption{Results for the Heckman sample selection model (synthetic data). The left table is for the estimation of $\theta=(\beta,\gamma,\sigma,\rho)$  while right table is for he estimation of $\beta$ only. For each experimental setting, we report the RMSE  over 25 replications.}
 \label{table2}
 \end{center}

\end{table}

\subsection{Heckman sample selection model\label{sub:Heck}}

Let $\Lambda=\R^2\times(0,\infty)\times (-1,1)$ and, for every  $\lambda=(\mu_1,\mu_2,\sigma,\rho)\in \Lambda$, let $P_\lambda$ be the distribution of $(Z_1,Z_2)$ where
\begin{align*}
 Z_{2}=\ind_{(0,\infty)}(Z_{2}^*),\quad Z_{1}=Z_{2}Z_{1}^*,\quad
 \begin{pmatrix}
 Z_{1}^*\\
  Z^*_{2}
  \end{pmatrix}
  \sim \mathcal{N}_2
\left(\begin{pmatrix}
\mu_1\\
\mu_2
\end{pmatrix},
\begin{pmatrix}
\sigma^2&\rho\sigma\\
\rho\sigma&1
\end{pmatrix}\right).
\end{align*}
Letting $d\in\mathbb{N}$,  $\Theta\subseteq\R^{d}\times(0,\infty)\times(-1,1)$ and
$$
g(\theta,x)=(\beta^\top x,\gamma^\top x,\sigma,\rho),\quad\forall x\in\R^{d},\quad\forall \theta=(\beta,\gamma,\sigma,\rho)\in\Theta,
$$
the Heckman sample selection model  is defined by  $\big\{ \big(P_{g(\theta,x)}\big)_{x\in\R^d},\,\theta\in\Theta\big\}$.

For this model,  in addition to the two MMD based estimators $\hat{\theta}_n$ and $\tilde{\theta}_n$, results are presented for $\theta_{\mathrm{mle},n}$, the maximum likelihood estimator   of $\theta$ (computed using the \texttt{R} package \texttt{sampleSelection} \cite{Toomet2008}) and for   $\theta_{\mathrm{rob},n}$, the robust two-step estimator of $\theta$ proposed by \cite{Zhelonkin2016}   (computed using the \texttt{R} package \texttt{ssmrob} \cite{Zhelonkin2013}).  We stress that the  estimator  $\theta_{\mathrm{rob},n}$ is designed specifically for robust estimation in Heckman sample selection models.

\subsubsection{Synthetic  data}
For this example we let   $d=8$,
$$
\Theta=\Big\{(\beta,\gamma,\sigma,\rho)\in \R^{d}\times(0,\infty)\times(-1,1) \text{ such that }\beta_{d-i+1}=\gamma_i=0,\,\forall i\in\{1,\dots,d/2\}\Big\}
$$
 and construct the dataset $d_N^0$ by simulating $N=5\,000$ independent   observations  using
$$
Y^0_i|X^0_i\sim P_{g(\theta_0,X^0_i)},\quad X^0_i \iid \mathcal{N}_{d}(0,I_{d})
$$
 where $\theta_0=(\beta_0 ,\gamma_0 ,\sigma_0,\rho_0)$ with
 $$
 \beta_0=(4,3,2,1,0,0,0,0),\quad \gamma_0=(0,0,0,0,4,3,2,1),\quad \sigma_0=1.5,\quad \rho_0=0.5.
 $$
Remark that the model we consider here assumes that the outcome equation depends only on the first $d/2$ components of $X_i^0$ while the selection equation depends only on the last $d/2$ components of $X_i^0$, so that the total number of parameter to estimate is $d+2=10$.

For $n\leq N$ and $\epsilon>0$ the contaminated dataset $d_n$ is constructed as explained in Section \ref{sub:sim_setup}, with $y_i^c=y_i$ and $x_i^c$ such that $x^c_{ij}=x_{ij}$ for all $j>1$ and such that $x_{i1}^c$ is a random draw from the from the $\mathcal{N}_1(5,1)$ distribution.

Table \ref{table2} shows, for different values of $n\leq N$ and of $\epsilon\in[0,0.03]$, the  RMSE  obtained with the different estimator for the estimation of $\theta_0$ as well as for the estimation of $\beta_0$, which is often the main parameter of interest in this model. We observe that the maximum likelihood estimator is the best estimator when there are no outliers, as expected from the asymptotic theory. On the other hand, this estimator is sensitive to the presence of outliers. The robust estimator $ \theta_{\mathrm{rob},n}$ of \cite{Zhelonkin2016} improves upon $ \theta_{\mathrm{mle},n}$ when the sample is contaminated. However,  for all the considered values of $\epsilon>0$ and  $n$, this estimator is dominated by $\hat{\theta}_n$ and by $\tilde{\theta}_n$. When restricting our attention to the estimation of $\beta$, we observe in  Table \ref{table2} that $\beta_{\mathrm{rob},n}$ dominates the two proposed estimator when $\epsilon=0.01$ and $n\in\{100,1\,000\}$ and when $\epsilon=0.02$ and $n=1\,000$. For all the other  18 pairs  $(\epsilon,n)$ with $\epsilon>0$ studied in Table \ref{table2}, the MMD based estimators $\hat{\beta}_n$ and $\tilde{\beta}_n$ outperform $\beta_{\mathrm{rob},n}$.

\begin{table}
\begin{center}
\begin{scriptsize}
  \begin{tabular}{|c|c|c|c|c||c|c|c|c||c|c|c|c|}
  \hline
  &\multicolumn{4}{c}{$\epsilon=0\%$}& \multicolumn{4}{c}{$\epsilon=1\%$}&\multicolumn{4}{c}{$\epsilon=3\%$}\\
  \hline
 & $\theta_{\mathrm{mle},n}$&$\theta_{\mathrm{rob},n}$& $\hat{\theta}_n$ &$\tilde{\theta}_n$& $\theta_{\mathrm{mle},n}$&$\theta_{\mathrm{rob},n}$& $\hat{\theta}_n$ &$\tilde{\theta}_n$& $\theta_{\mathrm{mle},n}$&$\theta_{\mathrm{rob},n}$& $\hat{\theta}_n$ &$\tilde{\theta}_n$\\
  \hline
  $\beta_1$&5.037 & 5.401&  5.157 & 5.111&5.445 & 5.441  & 5.183 & 5.090& 5.726 &-0.735 &  5.028 & 5.085   \\
   $\beta_2$& 0.212 &  0.201 & 0.200 &  0.212&0.180 & 0.190  &0.198  &0.213&0.172  & 0.120 & 0.207  & 0.214  \\
   $\beta_3$&0.350 & 0.255 &  0.436  & 0.318 &0.250 & 0.244  & 0.439  & 0.320& 0.193 & 0.665 & 0.442 & 0.319  \\
   $\beta_4$&0.019 & 0.013  & 0.025 &  0.027 &0.013  & 0.015  & 0.024  & 0.026 & 0.006 & 0.064  & 0.030 & 0.027  \\
   $\beta_5$&-0.220 &-0.155 & -0.193 &-0.210 &-0.129 & -0.131  &-0.188 &-0.200 & -0.096 & -0.388 &-0.193 &-0.212 \\
   $\beta_6$& 0.541 & 0.481 &  0.590&  0.526 &0.489  & 0.474   &0.590 & 0.517& 0.475  & 0.753  & 0.596  & 0.519 \\
   $\beta_7$& -0.029 & -0.067&  -0.084&  -0.084&-0.056 & -0.064 &-0.081 & -0.075 &-0.090  & 0.177& -0.075 & -0.078    \\
    $\gamma_1$&-0.724&  -0.749 & -1.067&  -0.827 &-0.613& -0.735 & -1.015 &-0.825&-0.508  & 5.441 & -0.631& -0.819    \\
    $\gamma_2$& 0.098 &  0.105&  0.147 & 0.116 &0.109 &  0.120 & 0.154 & 0.125 &0.106  & 0.190 & 0.124  & 0.120  \\
    $\gamma_3$& 0.644 & 0.687 &  0.661 & 0.736 &0.584 &  0.665 & 0.657 & 0.752 & 0.505 &  0.244 &  0.650 & 0.742  \\
    $\gamma_4$& 0.070  & 0.070 &  0.085 &  0.071&0.060 &  0.064 &  0.080 &  0.069& 0.053 & 0.015  & 0.060 & 0.068 \\
    $\gamma_5$& -0.373 &-0.398& -0.447 &-0.415 &-0.344& -0.388 &-0.453 & -0.417 & -0.302 & -0.131 & -0.465& -0.412  \\
    $\gamma_6$& 0.795 & 0.834 & 1.148 & 0.906 &0.559 &  0.753 & 1.162 &  0.888 & 0.344 & 0.474  & 1.152  & 0.869 \\
    $\gamma_7$& 0.182 &  0.183 & 0.146&  0.212 &0.179 & 0.177 & 0.144  & 0.214& 0.144 & -0.064 & 0.140 & 0.207 \\
    $\sigma$& 1.271 & 1.318 & 1.237 & 1.217 &1.351 &  1.357 & 1.243 & 1.229  & 1.454 & 1.357 & 1.251  & 1.242 \\
    $\rho$&  -0.124 & -0.051& -0.172 &-0.315  &-0.520& -0.047 & -0.184 & -0.307& -0.720 & -0.047 & -0.157& -0.308
\\
\hline
  \end{tabular}

  \vspace{0.3cm}

 \caption{Estimated value of $\theta$ in the Heckman sample selection model for the  2001 Medical Expenditure Panel Survey dataset. \label{table22}}

 \end{scriptsize}
 \end{center}
\end{table}

\subsubsection{Real  data}

We now consider the same dataset, and  model, than the one used in \cite[Table 4]{Zhelonkin2016} (and available from the \texttt{R} package \texttt{ssmrob}). The dataset contains $N=3\,328$ observations $d_N^0:=\{(y_{1i}^0,y_{2i}^0, x_i^0)\}_{i=1}^N$, extracted from the 2001 Medical Expenditure Panel Survey, and the Heckman sample selection model is used to regress $y^0_{i1}$, the  log ambulatory expenses for the $i$-th individual,  on the vector $x_i^0$ containing $d=7$ covariates (including an intercept). Both the selection and the outcome equation is assumed to depend on all the components of $x_i^0$, so that $\Theta=\R^{d}\times(0,\infty)\times(-1,1)$ and the number of parameter to estimate is   $2d+2=16$. We let $n=N$ and, for a given $\epsilon>0$, we contaminate  the dataset $d_N^0$, where    $y_{2i}^0=\ind_{(0,\infty)}(y_{1i}^0)$, by applying the approach described in Section \ref{sub:sim_setup}  with $(y_{1i}^c,y_{2i}^c, x_i^c)=(y_{1i}^0,1-y_{2i}^0, x_i^0)$. Below  we denote by $d_{n,\epsilon}$ the resulting contaminated version of $d_n^0$.

The estimated parameters values, obtained for  $\epsilon\in\{0,0.01,0.03\}$ and the four considered estimators,  are presented in Table \ref{table22}. For $\epsilon=0$, the results obtained for $\theta_{\mathrm{mle},n}$ and for $\theta_{\mathrm{rob},n}$ reproduce those given in  \cite{Zhelonkin2016}, and the two MMD based estimators provide   similar estimated values of the model parameters. In order to assess the sensibility of the different estimators to a contamination of the data, in Table \ref{table:SS_real} we give the value of $\Delta_{\epsilon,n}(\theta_n):=\|\theta_n(d_{n,\epsilon})-\theta_n(d^0_{n})\|$ for  $\theta_n\in\{\theta_{\mathrm{mle},n}, \theta_{\mathrm{rob},n},\hat{\theta}_n,\tilde{\theta}_n\}$. As expected, the maximum likelihood estimator is very sensitive to the presence of outliers and $\theta_{\mathrm{rob},n}$ improves upon $\theta_{\mathrm{mle},n}$. The most striking feature of Table \ref{table:SS_real} is the remarkable performance of $\tilde{\theta}_n$. Notably, the   value  $\Delta_{\epsilon,n}(\tilde{\theta}_n)$ is about 2.675 times smaller than that of $\Delta_{\epsilon,n}(\theta_{\mathrm{rob},n})$ when $\epsilon=0.01$, and about 8.5
times smaller when $\epsilon=0.03$. Finally, we observe that the impact of the contamination of the data on $\hat{\theta}_n$ is similar to that on $\theta_{\mathrm{rob},n}$. We however recall the reader that $\hat{\theta}_n$ depends on a parameter $\gamma_X>0$  which is kept constant throughout this Section \ref{section:simulation}, and thus which is not optimized (in any sense) to the particular problem at hand.

\begin{table}
\begin{center}
\begin{tabular}{c||c|c|c|c}
$\epsilon$&$\theta_{\mathrm{mle},n}$&$\theta_{\mathrm{rob},n}$&$\hat{\theta}_{n}$&$\tilde{\theta}_{n}$\\
\hline
1\%&0.651&0.107&0.063&0.040\\
3\%& 1.090&0.468&0.458&0.055
\end{tabular}
\caption{Value of $\Delta_{\tau,n}(\theta_n)$ for the Heckman sample selection model for the  2001 Medical Expenditure Panel Survey dataset, with $\theta_n\in\{\theta_{\mathrm{mle},n}, \theta_{\mathrm{rob},n},\hat{\theta}_{n},\tilde{\theta}_{n}\}$.\label{table:SS_real}}
\end{center}
\end{table}

\subsection{Gamma regression model\label{sub:gamma}}

We finally consider the Gamma regression model $\big\{ \big(P_{g(\theta,x)}\big)_{x\in\R^d},\,\theta\in\Theta\big\}$ obtained by taking $\Lambda=(0,\infty)^2$, $P_\lambda=Gamma(\lambda_1,\lambda_1/\lambda_2)$ for all  for  $\lambda\in\Lambda$, $\Theta=(0,\infty)\times \R^d$ and $g:\Theta\times\R^d\rightarrow \Lambda$  defined by
$$
g(\theta,x)=\big(\nu, \nu\exp(x^\top \beta)\big),\quad \theta=(\nu,\beta)\in\Theta.
$$

For this example we let $d=8$  and the dataset $d_N^0$ is obtained by simulating $N=5\,000$  independent   observations using $Y^0_i|X^0_i\sim P_{g(\theta_0,X^0_i)}$ and $X^0_i\iid \mathcal{N}_d(0,I_d)$, with $\theta_0=(1,\dots,1)$. Then, for every $n\leq N$ and $\epsilon\in [0,0.03]$, the contaminated dataset $d_n$ is constructed as described in Section \ref{sub:sim_setup},   with $y_i^c=y_i$ and $x_i^c$ such that $x^c_{ij}=x_{ij}$ for all $j>1$ and such that $x_{i1}^c$ is a random draw from the from $\mathcal{N}_1(-0.5,1)$ distribution.

Table \ref{table3} presents the RMSE for the estimation of $\theta_0$ and of $\beta_0$ obtained with     $\hat{\theta}_n$ and $\tilde{\theta}_n$. Results are also given for  $\theta_{\mathrm{mle},n}$, the maximum likelihood estimator of $\theta$, and for  $\theta_{\mathrm{rob},n}$, the estimator proposed by \cite{cantoni2001,cantoni2006} for robust inference in generalized linear model (computed using the \texttt{R} package \texttt{robustbase}). We observe from this table that, as expected, $\theta_{\mathrm{mle},n}$ is not robust to the presence of outliers, while the two proposed MMD based estimators are. The RMSEs obtained with the estimator $\theta_{\mathrm{rob},n}$ proposed by  \cite{cantoni2001,cantoni2006} is however slightly smaller than the ones obtained with $\hat{\theta}_n$ and $\tilde{\theta}_n$.

We stress that, in some sense, the better performance of $\theta_{\mathrm{rob},n}$ is reassuring, since this estimator is precisely  designed  for robust inference in generalized linear model and, in particular, has been motivated in  \cite{cantoni2006} for robust inference in Gamma regression models. By contrast, the applicability of the estimators $\hat{\theta}_n$ and $\tilde{\theta}_n$, and their theoretical guarantee   derived in Section \ref{section:theory} regarding their robustness, hold  for a much broader class of regression models.

\begin{table}[t]
\begin{center}
\begin{scriptsize}
  \begin{tabular}{|c|c|c|c|c|c|}
   \hline

   $\epsilon$ &$n$&  $\theta_{\mathrm{mle}}$ & $\theta_{\mathrm{rob}}$&  $\hat{\theta}_n$
          &  $\tilde{\theta}_n$
   \\

   \hline \hline
   \multirow{3}{*}{0\%} & $100$  & 0.478 &0.425&0.501 &0.498    \\
   &  $ 1\,000$  &0.123&0.120 &0.142&0.144 \\
   &  $5\,000$  &0.055&0.053&0.067&0.069     \\
    \hline
   \multirow{3}{*}{1\%} & $100$  &0.394&0.398 &0.521&0.521   \\
   &  $ 1\,000$  &0.250& 0.108 &0.149 &0.148\\
   &  $5\,000$  &0.350&0.049 &0.071 &0.070  \\
    \hline
     \multirow{3}{*}{2\%} & $100$  &0.395&0.381&0.539  &0.545  \\
   &  $ 1\,000$  &0.353 & 0.113& 0.146& 0.145 \\
   &  $5\,000$  &0.551 &0.056 &0.073&0.073 \\
    \hline
     \multirow{3}{*}{3\%} & $100$  &0.441&0.407&0.471 &0.475     \\
   &  $ 1\,000$  &  0.349 &0.119& 0.157& 0.158 \\
   &  $5\,000$  & 0.631&0.069 & 0.076 & 0.078 \\
    \hline
  \end{tabular} \hspace{0.5cm} \begin{tabular}{|c|c|c|c|c|c|}
   \hline

   $\epsilon$ &$n$&  $\beta_{\mathrm{mle}}$ & $\beta_{\mathrm{rob}}$&  $\hat{\beta}_n$
          &  $\tilde{\beta}_n$
   \\

   \hline \hline
   \multirow{3}{*}{0\%} & $100$  &0.308&0.306 &0.444 &0.443     \\
   &  $ 1\,000$  &0.087&0.093&0.132&0.133     \\
   &  $5\,000$  &0.041&0.043&0.061 &0.063       \\
    \hline
   \multirow{3}{*}{1\%} & $100$  &0.297 &0.313&0.473&0.475    \\
   &  $ 1\,000$  &0.114&0.098& 0.132&0.133  \\
   &  $5\,000$  &  0.070&0.042&0.061&  0.059 \\
    \hline
     \multirow{3}{*}{2\%} & $100$  &0.302 &0.301& 0.490&0.493 \\
   &  $ 1\,000$  & 0.120&0.094&0.127& 0.127  \\
   &  $5\,000$  & 0.105& 0.046&0.066& 0.067\\
    \hline
     \multirow{3}{*}{3\%} & $100$  &0.296 & 0.296&0.403& 0.406     \\
   &  $ 1\,000$  &0.119&  0.091&0.138& 0.136   \\
   &  $5\,000$  & 0.129& 0.044& 0.067& 0.069  \\
    \hline
    \end{tabular}
  \vspace{0.3cm}
 \caption{Results for the Gamma regression model. For each experimental setting, we report the mean square error over
25 replications.}
 \label{table3}
 \end{scriptsize}
 \end{center}
\end{table}

\section{Conclusion}
\label{section:conclusion}

We introduced a new family of parametric regression estimators based on MMD minimization. These estimators are universally $\sqrt{n}$-consistent and have strong robustness properties. Moreover, these estimators can be computed using a (stochastic) gradient algorithm whose cost per iteration is linear in the sample size $n$. Test on simulated data confirmed the good behaviour of these estimators in a wide range of regression models. Some questions remain open, including the dependence of the rate of convergence with respect to the dimension of the parameter space $\Theta$.

\section*{Acknowledgements}
We would like to acknowledge the Heilbronn Institute for Mathematical Research who funded a two weeks visit for P. A. to the University of Bristol, during which this work was started.

\small


TITLE: Homogeneous riemannian manifolds are complete. Trouble understanding proof.
QUESTION [16 upvotes]: I came across this proof while looking for hints on my homework, and I think it's only gotten me more confused. This is from Global Lorentzian Geometry.

Lemma 5.4 If $(H,h)$ is a homogeneous Riemannian manifold, then $(H,h)$ is complete.
Proof. By the Hopf-Rinow Theorem, it suffices to show that $(H,h)$ is geodesically complete. Thus suppose that $c:[a,1) \rightarrow H$ is a unit speed geodesic which is not extendible to $t=1$. Choosing any $p \in H$, we may find a constant $\alpha > 0$ such that any unit speed geodesic starting at $p$ has length $\ell \geq \alpha$. Set $\delta = \min\{\alpha/2, (1-a)/2\} > 0$. Since isometries preserve geodesics, it follows from the homogeneity of $(H,h)$ that any unit speed geodesic starting at $c(1-\delta)$ may be extended to a geodesic of length $\ell \geq 2 \delta$. In particular, $c$ may be extended to a geodesic $c:[a, 1+\delta) \rightarrow H$, in contradiction to the inextendibility of $c$ to $t=1$.

How exactly does it "follow from homogeneity" that we are able to extend unit speed geodesics? Isn't that what we're trying to show in the lemma? And how does considering geodesics based at $c(1-\delta)$ lead us to extending $c$ to $c:[a,1+\delta) \rightarrow H$?

REPLY [2 votes]: Another way to prove this is by Hopf Rinow's Theorem. I'll prove that a homogeneous Riemannian manifold is a complete metric space.
Let $\{x_n\}\subset M$ be a Cauchy sequence and fix $p\in M$. For each $n\in\mathbb{N}$ there is a $g_m\in \operatorname{Iso}(M)$ s.t. $g_m(x_m)=p$. Let $B_\epsilon (p)$ a  normal ball around $p$. There is a $N\in\mathbb{N}$ s.t. for every $n\geq N$,
$$\epsilon>d(x_N,x_n)=d(g_Np,g_np)=d(p,g_N^{-1}g_np)$$
So $d(p,x_n)=d(p,g_np)\in g_N(B_\epsilon(p))$. Then $\{x_n\}$ admits a convergente subsequence. As $x_n$ is a cauchy sequence, $x_n$ converges to the limit of the subsequence.