On the Adversarial Robustness of Multivariate Robust Estimation
Erhan Bayraktar Lifeng Lai



	
	On the Adversarial Robustness of Multivariate Robust Estimation
	
	
		
		
		
		[a]Department of Mathematics, University of Michigan, Ann Arbor, MI 48104.
			e1
		
		[b]Department of Electrical and Computer Engineering, University of California, Davis, CA 95616.
			e2,
		
		
		E. Bayraktar and L. Lai
		
		University of Michigan and University of California, Davis
		
	
	
	
		In this paper, we investigate the adversarial robustness of multivariate M-Estimators. In the considered model, after observing the whole dataset, an adversary can modify all data points with the goal of maximizing inference errors. We use adversarial influence function (AIF) to measure the asymptotic rate at which the adversary can change the inference result. We first characterize the adversary's optimal modification strategy and its corresponding AIF. From the defender's perspective, we would like to design an estimator that has a small AIF. For the case of joint location and scale estimation problem, we characterize the optimal M-estimator that has the smallest AIF. We further identify a tradeoff between robustness against adversarial modifications and robustness against outliers, and derive the optimal M-estimator that achieves the best tradeoff.
	
	
	
		Robustness
		adversarial attack
		M-estimator
		adversarial influence function
	
	




§ ORDINARY TEXT


Most of the existing work on robust statistical inference mainly address distributional robustness issues such as outliers or model uncertainties <cit.>. As machine learning and statistical inference algorithms are being increasingly used in safety critical applications and security related applications <cit.>, there is a growing interest in investigating the robustness of statistical inference algorithms in adversarial environments. In these adversarial environments, we are facing more severe situations than those addressed in the classic robust statistical inference problems. One such scenario is where an adversary can observe the whole dataset and then devise its attack vector to modify all entries in the data point hoping to cause the maximum inference error or to control the inference results. For example, the adversarial example phenomenon in the deep neural network <cit.> where an adversary can observe the whole picture and then carefully modify the pixels in the picture with the goal of fooling the classifier. As another example, it was shown in <cit.> that an adversary can modify the training data so that the model produced by the linear regression is controlled by the adversary. The existence of such powerful adversaries calls for new models and methodologies for adversarially robust inference.


In a canonical statistical inference problem, one infers parameters of interest from given data points {𝐱_1,⋯,𝐱_N}. The classic robust inference mainly deals with distributional robustness, i.e., the shape of the true underlying distribution deviates from the assumed model <cit.>. More specifically, let F_ be the cumulative density function (CDF) of the assumed model with  being parameter, then the classic robust inference deals with the situation where the data points 𝐱_n are independently and identically generated by an unknown CDF F in the η-neighborhood 𝒫_η(F_) of the assumed model F_. The goal of the classic robust inference is to design inference algorithms that perform well for any F∈𝒫_η(F_). For example, 𝒫_η(F_) can be a Lévy neighborhood 𝒫_η(F_)={F|∀ t,  F_(t-η)≤ F(t) ≤ F_(t+η)+η} or contamination neighborhood 𝒫_η(F_)={F| (1-η) F_+η H} in which H can be any probability measure. The contamination neighborhood model can be viewed as having η fraction of the data as outliers, while the Lévy neighborhood (or other related neighborhood) is useful in scenarios with model uncertainties. Various concepts such as influence function (IF), breakdown point, and change of variance etc were developed to quantify the robustness of estimators against the presence of outliers. 

In this paper, we consider a setup with more powerful adversaries than those considered in the classic robust inference and investigate adversarial robustness. In particular, in the considered setup, after data points {𝐱_1,⋯,𝐱_N} are generated, the adversary can observe the whole dataset and then modify all data points to {𝐱_1+Δ𝐱_1,⋯,𝐱_N+Δ𝐱_N} where each Δ𝐱_n is carefully designed and depends on the whole data set. It is easy to see that the adversary in the adversarially robust model is more powerful. In particular, for any F∈𝒫_η(F_) in the classic model, the adversary in the adverarially robust model can mimic the behavior of F by simply replacing the dataset with i.i.d samples generated from F. Clearly, this is not optimal strategy that the adversary in our model will adopt, as the adversary is not limited to this type of i.i.d attacks after observing the whole dataset. It can construct correlated attack signals that are based on the whole dataset. As the result, it is important to understand the following questions: 1) What is the attacker's optimal attack strategy in choosing Δ𝐱?; 2) What are the impacts of these attacks?; 3) How shall we design inference algorithms to minimize the impact?

In our recent work <cit.>, we made some progress in addressing these problems for the case of scalar parameter estimation, in which the parameter to be estimated is a scalar and each sample x_n is also a scalar. In particular, given a data set 𝐱={x_1,⋯,x_N} with x_n being i.i.d realizations of random variable X that has CDF F_θ(x) with unknown scalar parameter θ, we would like to estimate the unknown parameter θ. There is an adversary who can observe the whole dataset and can modify the dataset to 𝐱^Δ=𝐱+Δ𝐱:={x_1+Δ x_1, ⋯,x_N+Δ x_N}, in which Δ𝐱={Δ x_1, ⋯,Δ x_N} is the attack vector chosen by the adversary after observing 𝐱. Certain restrictions need to be put on Δ𝐱, otherwise the estimation problem will not be meaningful. In <cit.>, we assume that 

Δ𝐱∈𝒫_η:={1/N||Δ𝐱||_p^p≤δ^p},

in which ||·||_p is the ℓ_p norm. This type of constraints are reasonable and are motivated by real life examples. For example, in generating adversary examples in images <cit.>, the total distortion should be limited, otherwise human eyes will be able to detect such changes. 
The classic setup with contamination model <cit.> can be viewed as a special case of our formulation by letting p→ 0, i.e., the classic setup has a constraint on the total number of data points that the attacker can modify. For a given estimator, we would like to characterize how sensitive the estimator is with respect to the adversarial attack. In <cit.>, we considered a scenario where the goal of the attacker is to maximize the estimation error caused by the attack. We introduced a concept named “adversarial influence function” (AIF) to quantify the asymptotic rate at which the attacker can introduce estimation error through its optimal attack. From the defender's perspective, the smaller AIF is, the more adversarially robust the estimator is. In <cit.>, building on the characterization of AIF, we characterized the optimal estimator, among a certain class of estimators, that minimizes AIF. From this characterization, we show that there is a tradeoff between the robustness against adversarial attacks and robustness against outliers. In <cit.>, we further designed optimal estimator that achieve the optimal tradeoff among these quantities for the scalar case.

In this paper, we extend our work in <cit.> to multivariate setup, in which the goal is to jointly estimate multiple parameters from vector observations. The multivariate setup includes many important cases such as the joint location-scale estimation and robust linear regression etc. In this multivariate setup, we have N data points {𝐱_n, n=1,⋯,N} with each data point 𝐱_n∈ℝ^m being a vector. These data points are realizations of a random variable that has CDF F_() with unknown parameter vector ∈ℝ^q. We use m× N matrix :=[_1,⋯,_N] to denote the given data matrix. From this given data set, we would like to estimate the unknown parameter . The adversary will modify the data to ^Δ=+Δ𝐗, in which Δ is the attack matrix chosen by the adversary after observing . Similar to the situation in the classic robust inference problem <cit.>, the multivariate adversarial robustness setup is significantly more challenging than the scalar case.


Firstly, the characterization of the optimal attack strategy is much more difficult. There are many more degrees of freedom for the attacker to choose from, as the dimension of Δ𝐗 is m× N. Furthermore, each modification will affect all q components of the estimated vector  in a different but coupled manner. In this paper, we focus on the class of M-estimators specified by q-dimension functions . For this class of estimators, we characterize the optimal attack vector Δ𝐗 and the corresponding AIF. We further simplify this general formula for robust linear regression and evaluate the adversarial robustness of various existing robust algorithms. 

Secondly, the characterization of the optimal defense strategy is also much harder. For example, in the M-estimator case, now  is a q-dimension function, and the corresponding optimization problem of maximizing AIF becomes a coupled multi-dimension calculus of variation problem, which is in general very challenging. In this paper, for the important case of joint location-scale estimation problem, we show that we can decouple the characterization of optimal defense problem into two scalar problems. Building on this, we identify the optimal M-estimator that minimizes AIF.  In addition, similar to the scalar problem, we show that there exist a tradeoff between the robustness against adversarial attack and robustness again outliers. We further characterize the optimal M-estimator that achieves the optimal tradeoff between these robustness metrics.



The remainder of the paper is organized as follows. In Section <ref>, we introduce the considered model. In Section <ref>, we derive general AIF results and simplify the results for robust linear regression problems. In Section <ref>, we focus on the special case of joint location-scale estimation problem and characterize the optimal estimators that achieve the optimal AIF for both cases with and without constraints on the robustness against outliers. In Section <ref>, we use several numerical examples to illustrate results derived in this paper. Section <ref> provides concluding remarks. 



§ MODEL
 

In this section, we first introduce our problem formulation. We will then briefly review results from classic robust estimation that are directly related to our study.





 §.§ Problem Formulation


We have N data points {𝐱_n, n=1,⋯,N} with 𝐱_n:=[_n,1,⋯,_n,m]^T∈ℝ^m. These data points are i.i.d realizations of a random variable that has CDF F_() with unknown parameter ∈⊂ℝ^q. Here,  is a compact set. We will use f_() to denote the corresponding probability density function (pdf). We use m× N matrix :=[_1,⋯,_N] to denote the given data matrix. From this given data set, we would like to estimate the unknown parameter . However, as the adversary has access to the whole dataset, it will modify the data to ^Δ=+Δ𝐗:=[_1+Δ_1, ⋯,_N+Δ_N], in which Δ:=[Δ_1, ⋯,Δ_N] is the attack matrix chosen by the adversary after observing . We will discuss the attacker's optimal attack strategy in choosing Δ in the sequel. In this work, we consider the case where the attacker can modify all data points, which is a more suitable setup for recent data analytical applications. However, certain restrictions need to be put on Δ, otherwise the estimation problem will not be meaningful. In this paper, we assume that 

    1/mN||Δ||_p^p≤δ^p,

in which ||𝐀||_p is the entry-wise matrix p-norm:

    ||𝐀||_p=||Vec(𝐀)||_p=(∑∑ |𝐀_i,j|^p)^1/p,

with Vec(𝐀) being the vectorization of matrix 𝐀. In (<ref>), we have the normalization term 1/mN as the matrix Δ is of size m× N. The normalization factor implies that the per-dimension change (on average) is upper-bound by δ^p. As mentioned in the introduction, this type of constraints are reasonable and are motivated by real life examples. 


Following notation used in robust statistics <cit.>, we will use _N(), a q dimensional vector, to denote an estimator. For a given estimator _N, we would like to characterize how sensitive the estimator is with respect to the adversarial attack. In this paper, we consider a scenario where the goal of the attacker is to maximize the estimation error caused by the attack. In particular, the attacker aims to choose Δ by solving the following optimization problem

    max_Δ       ||_N(+Δ)-_N()||_1,
    s.t.      1/mN||Δ||_p^p≤δ^p,

in which ||·||_1 is the ℓ_1 norm.

We use f(_N,,p,δ) to denote the optimal value obtained from the optimization problem (<ref>), and define the adversarial influence function (AIF) of estimator _N at  under ℓ_p norm constraint as

    AIF(_N,,p)=lim_δ↓ 0f(_N,,p,δ)/δ.

This quantity, a generalization of the concept of influence function (IF) used in classic robust estimation (a brief review of IF will be provided in Section <ref>), quantifies the asymptotic rate at which the attacker can introduce estimation error through its attack. 

From the defender's perspective, the smaller AIF is, the more robust the estimator is. In this paper, building on the characterization of AIF(_N,,p), we will characterize the optimal estimator _N, among a certain class of estimators 𝒯, that minimizes AIF(_N,,p). In particular, we will investigate 

    min__N∈𝒯AIF(_N,,p).



Note that AIF(_N,𝐗,p) depends on the data matrix 𝐗. Based on the characterization of AIF for a given data realization matrix 𝐗 with N columns (each column representing one data point), we will then study the population version of AIF where each column of 𝐗 is i.i.d generated by F_. We will examine the behavior of AIF(_N,𝐗,p) as N increases. We will see that for a large class of estimators AIF(_N,𝐗,p) has a well-defined limit as N→∞. We will use AIF(,F_,p) to denote this limit when it exists.




From the defense's perspective, we would like to design an estimator that is least sensitive to the adversarial attack. Again, we will characterize the optimal estimator , among a certain class of estimators 𝒯, that minimizes AIF(,F_,p). That is, for a certain class of estimators 𝒯, we will solve

    min_∈𝒯AIF(,F_,p).


It will be clear in the sequel that the solution to the optimization problem (<ref>), even though is robust against adversarial attacks, has poor performance in guarding against outliers. This motivates us to design estimators that strike a desirable tradeoff between these two robustness measures. In particular, we will solve (<ref>) with an additional constraint on IF. After the corresponding quantities are introduced in later sections, precise formulation of this optimization problem with additional IF constraint will be given in (<ref>). 



 §.§ M-Estimator and Influence Function (IF)


In this paper, we will mainly focus on a class of commonly used estimator in robust statistic: M-estimator <cit.>, in which one obtains an estimate _N() of  by solving 

    ∑_n=1^N (_n,_N)=0.


Here (,):ℝ^m×ℝ^q→ℝ^q is a vector function of data point  and parameter  to be estimated. We use ψ_i, i=1,⋯,q, to denote each component of . Different choices of  lead to different robust estimators. For example, the most likely estimator (MLE) can be obtained by setting =-f_^'/f_. 

As the form of  determines _N, in the remainder of the paper, we will use  and _N interchangeably. For example, we will denote AIF(_N,,p) as AIF(,,p). 

It is typically assumed that (,) is continuous and almost everywhere differentiable. This assumption is valid for all 's that are commonly used. It is also typically required that the estimator is Fisher consistent <cit.>:

    𝔼_F_[(,)]=0,

in which 𝔼_F_ means expectation under F_. Intuitively speaking, this implies that the true parameter  is the solution of the M-estimator if there are increasingly more i.i.d. data points generated from F_.

In the contamination model of the classic robust estimation setup, it is assumed that a fraction δ of data points are outliers, while the remainder of data points are generated from the true distribution F_. For a given estimator , the concept of IF introduced by Hamper <cit.> is defined

    IF(,,F_)=lim_δ↓ 0((1-δ)F_+δ i_)-(F_)/δ.

Here, IF(,,F_) is a q dimensional vector. In this definition, i_ is a distribution that puts mass 1 at point . In addition, (F_) is the obtained estimate when all data points are generated i.i.d from F_, and ((1-δ)F_+δ i_) is the obtained estimate when 1-δ fraction of data points are generated i.i.d from F_ while δ fraction of the data points are at . Hence, IF(,,F_) measures the asymptotic influence of having outliers at point  as δ↓ 0. Similar as above, as  is determined by  in M-estimator, in the following, we will also denote IF(,,F_) as IF(,,F_).

Furthermore, to characterize the impact of the worst outliers, Hamper <cit.> introduced the (unstandardized) gross-error sensitivity:

    γ^*_u(,F_)=sup_{||IF(,,F_)||_2},

in which ||·||_2 is the ℓ_2 norm.

For M-estimator, IF(,, F_) was shown to be <cit.> 

    IF(,, F_)=(,F_)^-1(,(F_)),

with the q× q matrix  given by

    (,F_)=-∫.∂/∂(𝐱,)|_(F_)dF_(),

see (4.2.9) of <cit.>. 



§ CHARACTERIZING AIF


In this section, for a given data matrix , we analyze the AIF for any given M-estimator  as specified in (<ref>). As  is q-dimension vector, there are q equations. To simplify the presentation, we write each equation as 
    G_i(,_N):=∑_n=1^N ψ_i(_n,_N)=0,i=1,⋯,q
 and denote 𝐆=[G_1,⋯,G_q]^T. Using this notation, (<ref>) can be written as

    𝐆(,_N)=0.




 §.§ General Case


To proceed further, we write

    𝐆^'_(,_N)   =   .∂𝐆/∂|_,_N,
    𝐆^'_(,_N)   =   .∂𝐆/∂Vec(𝐗)|_,_N,

and

    ^'_N()=.∂𝐭_N/∂Vec(𝐗)|_,_N.



We have the following characterization.


	
	
	
	
	
	
	
	
	For p≥ 1, suppose 𝐆^'_(,_N) is invertible, we have
	
    (,,p)=(mN)^1/pmax_∈||^T[𝐆^'_(,_N)]^-1𝐆^'_(,_N)||_p/p-1,

	in which  is the set of length q vectors with each entry being either 1 or -1.


	
	First, from (<ref>), we have 
	
    _N^'()=-[𝐆^'_(,_N)]^-1𝐆^'_(,_N).

	
	Using Taylor expansion, we have
	
    _N(+Δ)-_N()=_N^'()Vec(Δ)+higher order terms.

	
	
	When δ is small, the adversary can focus on the following optimization problem to obtain an o(δ) optimal solution
	
    (P1):max_Vec(Δ)      ||_N^'()Vec(Δ)||_1,
    s.t.       ||Δ||_p^p≤ Nmδ^p.

	
	Let f^* be the optimal value obtained in the optimization problem (P1). Let g^* be the optimal value of the following optimization problem
	
    (P2):max_∈max_Vec(Δ)      ^T_N^'()Vec(Δ),
    s.t.       ||Δ||_p^p≤ Nmδ^p.

	
	In Appendix <ref>, we show that f^*=g^*. Hence, we can focus on problem (P2). 
	
	The inner maximization problem in (P2) is the same as
	
    min_Δ      -^T_N^'()Vec(Δ),
    s.t.       ||Δ||_p^p≤ Nmδ^p.

	
	Using (<ref>), we have 
	
    -^T_N^'()Vec(Δ)=^T[𝐆^'_(,_N)]^-1𝐆^'_(,_N)Vec(Δ).
	
	
	To simplify the notation, we denote
	
    𝐚=^T[𝐆^'_(,_N)]^-1𝐆^'_(,_N),

	which is a row vector with mN entries. Even though 𝐚 is only a row vector, we denote these elements as a_i,j for i=1,⋯,m and j=1,⋯,N to better connect with each elements of Δ. Hence, a_i,j corresponds to Δ_i,j. Using this notation, the optimization problem can be written as 
	
    min_Δ      𝐚Vec(Δ),
    s.t.       ||Δ||_p^p≤ mNδ^p.

	
	For p=1, this is a linear programing problem, whose solution is simple. In particular, let a^*=max_i,j | a_i,j|, and (i^*,j^*) be the corresponding index, it is easy to check that we have
	
    Δ_i^*,j^*^*=-sign{ a_i^*,j^*}mNδ,

	and Δ_i,j^*=0 for other i,js. Hence,
	
    (,,p)=max_∈mNa^*=mNmax_∈||𝐚||_∞.

	
	
	For p> 1, (<ref>) is a convex optimization problem. To solve this, we form Lagrange
	
    ℒ(Δ,λ)= 𝐚Vec(Δ)+ λ(||Δ||_p^p-mNδ^p).

	
	The corresponding optimality conditions are:
	
    a_i,j+λ^* p sign (Δ_i,j^*) |Δ_i,j^*|^p-1   =   0, ∀ n
    λ^*   ≥   0,
    λ^*(||Δ||_p^p-mNδ^p)   =   0.

	
	From (<ref>), we know that λ^*≠ 0, hence 
	
    ||Δ^*||_p^p=mNδ^p,

	and
	
    sign (Δ_i,j^*) |Δ_i,j^*|^p-1=-a_i,j/λ^* p.

	
	From (<ref>) and the fact that λ^* p is positive, we know sign (Δ_i,j^*)=-sign(a_i,j), and hence we have
	
    |Δ_i,j^*|^p-1=|a_i,j|/λ^* p,

	which can be simplified further to
	
    Δ_i,j^*=-(|a_i,j|/λ^* p)^1/(p-1)sign (a_i,j).

	
	Combining these with (<ref>), we obtain the value of λ^*:
	
    λ^*=1/p(∑_n=1^N |a_i,j|^p/(p-1)/mNδ^p)^(p-1)/p.

	
	As the result, we have
	
    Δ_i,j^*=-|a_i,j|^1/(p-1)(mN)^1/p/(∑|a_i,j|^p/(p-1))^1/psign(a_i,j)δ.

	
	Hence, the optimal value of the inner maximization of (P2) is 
	
    ∑_i,j a_i,j|a_i,j|^1/(p-1)(mN)^1/p/(∑|a_i,j|^p/(p-1))^1/psign(a_i,j)	=(mN)^1/p||𝐚||_p/(p-1),

	which finishes the proof.

By setting m=1 and q=1, (<ref>) recovers the result on the scalar case presented in <cit.>.




 §.§ Robust Regression


In this section, we use robust linear regression, an important multivariate parameter estimation problem, as an example to illustrate the result derived in Section <ref>. In linear regression problems, the data points are _n=([ 𝐱_n; y_n ]), n=1,⋯,N with 𝐱_n∈ℝ^q and y_n∈ℝ. Hence, _n∈ℝ^q+1. In the following, we let =[_1,⋯,_N]∈ℝ^(q+1)× N, and still denote =[_1,⋯,_N]∈ℝ^q× N. From the data, we would like to fit a linear model, i.e., we would like to find =[θ_1,⋯,θ_q]^T∈ℝ^q such that 𝐱_n^T is a good approximation of y_n. Hence, the parameters to be estimated are ∈ℝ^q. Furthermore, each data point _n∈ℝ^q+1, hence m=q+1. We denote 
    r_n=y_n-𝐱_n^T
 as the residual error.  

The commonly used ordinary least square (OLS) approach finds  by solving

    min_∑_n=1^N r_n^2,

which is equivalent to solving

    ∑_n=1^N r_n𝐱_n=0.

The solution is well known _OLS=(^T)^-1𝐲. In the subsequent discussion, we will use a related quantity named hat matrix 
    =^T(^T)^-1.


It is known that OLS solution is not robust to outliers <cit.>. Various robust linear regression schemes were proposed <cit.>. They generally set (_n,_N) in the form η(r_n v_n)w_n_n 
with function η: ℝ→ℝ and weights w_n and v_n. That is, for robust linear regression, one obtains the estimate _N of  by solving

    ∑_n=1^N(_n,_N)=∑_n=1^N η(r_n v_n)w_n_n=0.


The weights w_n and v_n can be chosen to not only depend on _n but also the whole data matrix . For example, it is common <cit.> to use w_n=√(1-h_nn), in which h_nn=_n^T(^T)^-1_n is the nth diagonal element of the hat matrix (<ref>). It is known <cit.> that 0≤ h_nn≤ 1. Comparing (<ref>) with (<ref>), we can see that one replaces _n in (<ref>) with w_n_n and replaces r_n in (<ref>) with η(r_n v_n). The main idea is to use w_n to limit the impact of outliers in _n and use η(r_n v_n) to limit the impact of outliers in the residual r_n while taking the location of _n into consideration. From (<ref>), we have 

    ψ_i(_n,_N)   =   η(r_n v_n)w_n_n,i=η((y_n-𝐱_n^T_N)v_n)w_n_n,i,
    
    G_i(,_N)   =   ∑_n=1^Nψ_i(_n,_N)=∑_n=1^Nη((y_n-𝐱_n^T_N)v_n)w_n_n,i.


Different choices of functions lead to different classes of robust linear regression methods. For example:

	
  * w_n=1, v_n=1 leads to Huber's proposal <cit.>, in which the idea was to replace r_n in (<ref>) with η(r_n) so as to limit the influence of large residuals (similar to the M-estimator in a single variable case). 
	
  * v_n=1 leads to Mallows's proposal.
	
	
  * v_n=1/w_n is Schewppe's approach <cit.>. 


In the following, we calculate 𝐆^'_(,_N) and 𝐆^'_(,_N).

First, we compute 𝐆^'_(,_N). For i=1,⋯,q; j=1,⋯,N; and k=1,⋯,q, we have

    .∂ G_i/∂𝐱_j,k|_,_N   =   ∂(∑_n=1^N w_nη(r_n v_n)𝐱_n,i)/∂𝐱_j,k
       =   ∑_n=1^N [ η(r_n v_n) 𝐱_n,i∂  w_n/∂𝐱_j,k+w_n 𝐱_n,iη^'(r_n v_n)(-_N(k)v_n𝕀(j=n)+r_n∂ v_n /∂𝐱_j,k).
          .+w_nη(r_n v_n)𝕀(j=n,k=i)],

in which 𝕀(·) is the indicator function and _N(k) is the kth element of . In addition,

    .∂ G_i/∂ y_j|_,_N=w_jη^'(r_jv_j)v_j𝐱_j,i:=c_j𝐱_j,i,

in which we denote 
    c_j=w_jη^'(r_j v_j)v_j.


Furthermore, each entry of 𝐆^'_(,_N) can be computed as

    .∂ G_i/∂θ_j|_,_N=-∑_n=1^N w_nη^'(r_n v_n)v_n𝐱_n,i𝐱_n,j:=-∑_n=1^N c_n𝐱_n,i𝐱_n,j.

From this, we know that

    𝐆^'_(,_N)=-diag[c_1,⋯,c_N]^T.


Using the result in Theorem <ref>, assuming 𝐆^'_(,_N) is invertible, we have the following characterization of AIF of robust linear regression.

	
	
	
	
	
	
	
	
	For robust linear regression,
	
    (,,p)=(mN)^1/pmax_∈||^T[diag[c_1,⋯,c_N]^T]^-1𝐆^'_(,_N)||_p/p-1,

	with c_j defined in (<ref>).


We now apply Proposition <ref> to various specific (robust) linear regression approaches.

OLS: For OLS, we have w_n=1, v_n=1 and η(r)=r. In this case, we have

    .∂ G_i/∂𝐱_j,k|_,_N=-𝐱_j,i_N(k)+r_j𝕀(k=i),

and
.∂ G_i/∂ y_j|_,_N=𝐱_j,i.
Furthermore, c_j=1 and hence
𝐆^'_(,_N)=-^T.

Huber's proposal <cit.>: In Huber's proposal <cit.>, we have w_n=1, v_n=1. In this case, we have

    .∂ G_i/∂𝐱_j,k|_,_N=-𝐱_j,i_N(k)η^'(r_j)+η(r_j)𝕀(k=i),

and
.∂ G_i/∂ y_j|_,_N=η^'(r_j)𝐱_j,i. Furthermore
𝐆^'_(,_N)=-diag[η^'(r_1),⋯,η^'(r_N)]^T.


Mallow's proposal: In Mallow's proposal, we have v_n=1. In this case, we have

    .∂ G_i/∂𝐱_j,k|_,_N   =   ∑_n=1^N [ η(r_n ) 𝐱_n,i∂  w_n/∂𝐱_j,k]-w_j𝐱_j,i_N(k)η^'(r_j)+w_jη(r_j)𝕀(k=i),

and
.∂ G_i/∂ y_j|_,_N=w_jη^'(r_j)𝐱_j,i. Furthermore,
𝐆^'_(,_N)=-diag[c_1,⋯,c_N]^T
with
c_j=w_jη^'(r_j).

Schewppe's approach <cit.>: In Schewppe's approach, we have v_n=1/w_n. In this case, we have

    .∂ G_i/∂𝐱_j,k|_,_N   =   ∑_n=1^N [ η(r_n v_n) 𝐱_n,i∂  w_n/∂𝐱_j,k-w_n 𝐱_n,iη^'(r_n v_n)(_N(k)v_n𝕀(j=n)+r_n/w_n^2∂ w_n /∂𝐱_j,k).
          .+w_nη(r_n v_n)𝕀(j=n,k=i)].

and
.∂ G_i/∂ y_j|_,_N=η^'(r_jv_j)𝐱_j,i. Furthermore, c_j=η^'(r_jv_j) and
𝐆^'_(,_N)=-diag[c_1,⋯,c_N]^T.

We will compare these methods numerically in Section <ref>.



§ OPTIMAL ADVERSARIAL ROBUSTNESS VS OUTLIER ROBUSTNESS TRADEOFF


In this section, we specialize the results to the joint estimation of location and scale. Building on these results, we will design M-estimators that minimizes AIF or achieves the optimal tradeoff between AIF, i.e., adversarial robustness, and IF, i.e., outlier robustness.


In the joint location-scale estimation, given {x_n,n=1,⋯,N} with x_n∈ℝ, 
the goal is to jointly estimate the location parameter θ_1 and the scale parameter θ_2. We will that assume θ_1 is bounded and there is a constant c such that |θ_2|>c. Hence, the dimension of each data point m=1 and the dimension of parameter q=2. We focus on a large class of model named the location-scale model <cit.>. In the location-scale model, we have 

    F_(X)=ℒ(θ_1+θ_2Z),

in which Z is a random variable with symmetric pdf f_0 and CDF F_0, and ℒ(θ_1+θ_2Z) means the distribution of θ_1+θ_2Z. This class of model includes many important models in statistics. For example the joint estimation of mean and variance of Gaussian random variable belongs to this model. Another important example is the linear regression discussed in Section <ref>. 

Following the convention, we use T_N to denote the estimation of location θ_1 and S_N to denote the estimation of the scale θ_2. In the location-scale model, one typically obtains (T_N,S_N) by solving the following equations <cit.>:

    ∑_n=1^Nψ_1(x_n-T_N/S_N)   =   0,
    ∑_n=1^Nψ_2(x_n-T_N/S_N)   =   0,

with properly chosen ψ_1 and ψ_2. In the following, for presentation and notation convenience, we denote 
    z_n:=x_n-T_N/S_N.

Hence, 
    (x,)=([ ψ_1(z); ψ_2(z) ]),
 
with z=x-θ_1/θ_2.

For this class of model, one typically focuses on equivariant M-estimator <cit.> with 
    ψ_1(-z)=-ψ_1(z)  and ψ_2(-z)=ψ_2(z).

This implies that the first component of  is an odd function, while the second component of  is an even function. Furthermore, ψ_1 and ψ_2 are assumed to be monotone functions in z≥ 0. Without loss of generality, we will focus on monotone increasing functions, hence ψ_1^'(z)≥ 0 and ψ_2^'(z)≥ 0 for z≥ 0. 



 §.§ Given Sample Case

We first use the results derived in Section <ref> to characterize the AIF for a given data matrix . We note that in this joint location-scale estimation problem, m=1, q=2, and

    G_i=∑_n=1^Nψ_i(z_n).


We have

    .∂ G_i/∂ x_n|_,_N   =   1/S_Nψ_i^'(z_n),.∂ G_i/∂θ_1|_,_N=-1/S_N∑_n=1^Nψ_i^'(z_n),
    .∂ G_i/∂θ_2|_,_N   =   -1/S_N∑_n=1^Nz_nψ_i^'(z_n).


Hence,

    𝐆^'_(,_N)   =   1/S_N[[ ψ_1^'(z_1),         ,⋯,  ψ_1^'(z_N); ψ_2^'(z_1),          ⋯,  ψ_2^'(z_N) ]],
    𝐆^'_(,_N)   =   -1/S_N[[   ∑_n=1^Nψ_1^'(z_n), ∑_n=1^Nz_nψ_1^'(z_n);   ∑_n=1^Nψ_2^'(z_n), ∑_n=1^Nz_nψ_2^'(z_n) ]]:=-1/S_N[[ a,  b; c,  d ]].

Since, q=2, ={(1,1)^T,(-1,-1)^T,(1,-1)^T,(-1,1)^T}. Due to symmetry, we only need to consider  being either (1,1)^T or (1,-1)^T. When ∈{(1,1)^T,(1,-1)^T}, we have

    ^T[𝐆^'_(,_N)]^-1𝐆^'_(,_N)
          =-[(d-(2)c)ψ_1^'(z_1)+((2)a-b)ψ_2^'(z_1),⋯,(d-(2)c)ψ_1^'(z_N)+((2)a-b)ψ_2^'(z_N)]/ad-bc.












Hence, using Theorem <ref>, for p>1, we have

    (,,p)
       =   max_N^1/p||^T[𝐆^'_(,_N)]^-1𝐆^'_(,_N)||_p/p-1
       =   N^1/p/|ad-bc|max_(2)∈{1,-1}[∑_n=1^N|(d-(2)c)ψ_1^'(z_n)+((2)a-b)ψ_2^'(z_n)|^p/(p-1)]^(p-1)/p
       =   N/|ad-bc|max_(2)∈{1,-1}[1/N∑_n=1^N|(d-(2)c)ψ_1^'(z_n)+((2)a-b)ψ_2^'(z_n)|^p/(p-1)]^(p-1)/p
       ≥   N/|ad-bc|max_(2)∈{1,-1}1/N∑_n=1^N|(d-(2)c)ψ_1^'(z_n)+((2)a-b)ψ_2^'(z_n)|
       ≥   1/|ad-bc|max_(2)∈{1,-1}|∑_n=1^N((d-(2)c)ψ_1^'(z_n)+((2)a-b)ψ_2^'(z_n))|
       =   1.




Summarizing the discussion, we have the following proposition.

	For joint location-scale estimation, 
	
    (,,p)=N^1/p/|ad-bc|max_(2)∈{1,-1}[∑_n=1^N|(d-(2)c)ψ_1^'(z_n)+((2)a-b)ψ_2^'(z_n)|^p/(p-1)]^(p-1)/p.

	
	
	Furthermore, we have
	
    (,,p)≥ 1.

	


In the following, we provide several examples to analyze AIF of existing schemes.


	The coupled mean and sample standard deviation estimator (page 233 of <cit.>) is specified by
	ψ_1(z)=z and ψ_2(z)=z^2-1. For this estimator, we have ψ_1^'(z_n)=1, ψ_2^'(z_n)=2z_n, and
	
    ∑_n=1^Nz_n=0,   and  ∑_n=1^N(z_n^2-1)=0.

	As the result, 
	
    a   =   ∑_n=1^Nψ_1^'(z_n)=N,
    
    	b   =   ∑_n=1^Nz_nψ_1^'(z_n)=∑_n=1^Nz_n=0,
    
    	c   =   ∑_n=1^Nψ_2^'(z_n)=∑_n=1^N2z_n=0,
    
    	d   =   ∑_n=1^Nz_nψ_2^'(z_n)=∑_n=1^N2z_n^2=2N.

	Hence, 
	
    (,,p)   =   N^1/p/|ad-bc|max_(2)∈{1,-1}[∑_n=1^N|(d-(2)c)ψ_1^'(z_n)+((2)a-b)ψ_2^'(z_n)|^p/(p-1)]^(p-1)/p
       =   N^1/p/2N^2max_(2)∈{1,-1}[∑_n=1^N|2N+2N(2)z_n|^p/(p-1)]^(p-1)/p
       =   max_(2)∈{1,-1}[1/N∑_n=1^N|1+(2)z_n|^p/(p-1)]^(p-1)/p.

	In particular, when p=2, we have
	
    (,,2)   =   max_(2)∈{1,-1}[1/N∑_n=1^N(1+(2)z_n)^2]^1/2=√(2).





	Combination of the asymptotic minimax estimates of location and scale:
	
    ψ_1(z)=max[-K,min(K,z)],  and  ψ_2(z)=min(α^2,z^2)-β,

	where 0<β<α^2. The parameters K, α and β are properly chosen so that the corresponding location and scale estimators are minimax estimators respectively. If one further sets α=K, this estimator corresponds to Huber's Proposal 2 <cit.>.	For this joint estimator, we have
	
    ∑_n=1^Nψ_1(z_n)   =   ∑_n∈𝒜_-:={z_n≤ -K}-K+∑_n∈𝒜:={|z_n|<K}z_n+∑_n∈𝒜_+:={z_n≥ K}K=0,
    ∑_n=1^Nψ_2(z_n)   =   ∑_n∈ℬ:={z_n^2≤ c^2}(z_n^2-β)+∑_n∈ℬ̅(α^2-β)=0,

	ψ_1^'(z_n)=𝕀(n∈𝒜) and ψ_2^'(z_n)=2z_n𝕀(n∈ℬ).
	
	
	As the result,
	
    a   =   ∑_n=1^Nψ_1^'(z_n)=|𝒜|,
    
    	b   =   ∑_n=1^Nz_nψ_1^'(z_n)=∑_n∈𝒜z_n=K[|𝒜_-|-|𝒜_+|],
    
    	c   =   ∑_n=1^Nψ_2^'(z_n)=2∑_n∈ℬz_n,
    
    	d   =   ∑_n=1^Nz_nψ_2^'(z_n)=∑_n∈ℬ2z_n^2=2(Nβ-|ℬ̅|c^2).

	
	Plugging these into Proposition <ref>, we can obtain the formulator for AIF.
	


In Section <ref>, we will see that these existing schemes are not optimal in terms of minimizing AIF with or without additional IF constraints. We will derive optimal schemes that minimizes AIF with or without IF constraint in Section <ref>.



 §.§ Population Case


In this section, we analyze the behavior of AIF as N→∞. Based on this analysis, we will characterize the optimal  that minimizes AIF. We will further identify a tradeoff between AIF and IF, and design  that achieves this optimal tradeoff.


Using Proposition <ref>, we have

    (,,2)
       =   N/|ad-bc|max_(2)∈{1,-1}[1/N∑_n=1^N|(d-(2)c)ψ_1^'(z_n)+((2)a-b)ψ_2^'(z_n)|^2]^1/2
       =   1/|ãd̃-b̃c̃|max_(2)∈{1,-1}[∑_n=1^N((d̃-c̃)^2ψ_1^'(z_n)^2+2(ã-b̃)((d̃-c̃))ψ_1^'(z_n)ψ_2^'(z_n)+(ã-b̃)^2ψ_2^'(z_n)^2)/N]^1/2,

in which ã=(2)a/N,b̃=b/N,c̃=(2)c/N,d̃=d/N. 

It has been shown in Theorem 2.4 of <cit.> that, under certain mild regularity conditions, T_Na.s.→θ_1 and S_Na.s.→θ_2. In the following, we will need the following additional regularity conditions:

	
  * ψ_1^'(z) and ψ_2^'(z) are continuous functions. 
	
  * There exist a function K(z) such that |ψ_1^'(z)|≤ K(z), |zψ_1^'(z)|≤ K(z), |ψ_2^'(z)|≤ K(z), |zψ_2^'(z)|≤ K(z), and 𝔼[K(Z)]<∞.
 
Recall that the CDF of Z is F_0, hence 𝔼{·} here (and the subsequent discussion) denotes expectation under F_0. The conditions here are slightly stronger than those conditions needed for the strong law of large numbers, as we will need to use the uniform strong law of large numbers (see Theorem 16 (a) <cit.>). Under these regularity assumptions, using the uniform strong law of large numbers, Slutsky Theorem (see Chapter 6 of <cit.>) and the fact that T_Na.s.→θ_1 and S_Na.s.→θ_2, we have 

    ã   =   (2)a/Na.s.→(2)𝔼{ψ_1^'(Z)},
    d̃   =   d/Na.s.→𝔼{Zψ_2^'(Z)}.

Furthermore, 

    b̃=b/Na.s.→𝔼{Zψ_1^'(Z)}=0,

as ψ_1^'(z) is an even function (since ψ_1(z) is an odd function) and f_0 is symmetric. Similarly, c̃=(2)c/Na.s.→𝔼{ψ_2^'(Z)}=0 since ψ_2(Z) is an even function. 



As a result, 

    (,,2)a.s.→√(𝔼{ψ_1^'(Z)^2}/(𝔼{ψ_1^'(Z)})^2+𝔼{ψ_2^'(Z)^2}/(𝔼{Zψ_2^'(Z)})^2):=(,F_0,2).

In deriving this equation, we use the fact that ψ_1^'(Z)ψ_2^'(Z) is an odd function, hence 
    1/N∑_n=1^Nψ_1^'(z_n)ψ_2^'(z_n)a.s.→𝔼{ψ_1^'(Z)ψ_2^'(Z)}=  0.




  §.§.§ Minimizing (,F_0,2) 


From the defender's perspective, one would like to design  such that (,F_0,2) is small. We now characterize ψ_1 and ψ_2 that minimize (,F_0,2). As (<ref>) can be decomposed into two independent terms, we can minimize ψ_1(z) and ψ_2(z) separately. Since ψ_1(z) is an odd function and ψ_2(z) is an even function, we only need to characterize the functions for z≥ 0.

To obtain the optimal ψ_1(z), due to Jensen's inequality, we have that

    𝔼{ψ_1^'(Z)^2}/(𝔼{ψ_1^'(Z)})^2≥ 1

for which the equality holds when ψ_1^'(z) is a constant. By setting ψ_1(0)=0, this function satisfies the requirements on ψ_1(z) (i.e., ψ_1(z) is Fisher consistent and is an odd function).

For ψ_2(z), we need to solve

    min      𝔼{ψ_2^'(Z)^2}/(𝔼{Zψ_2^'(Z)})^2=∫_0^∞ψ_2^'(z)^2f_0(z) z/(∫_0^∞zψ_2^'(z)f_0(z) z)^2, 
    s.t.      𝔼{ψ_2}=2∫_0^∞ψ_2(z)f_0(z) z=0,
          ψ_2^'(z)≥ 0.

in which we use the fact that ψ_2^'(z) is an odd function and f_0(z) is a symmetric function, while we use the requirement that ψ_2(z) is an even function. Here, (<ref>) is the Fisher consistency requirement.

	The optimal ψ_2(z) is given by the following form: 
	
		
  * ψ_2^'(z)=z/∫_0^∞z^2f_0(z)dz.
		
  * ψ_2(0)=-∫_0^∞(∫_0^zψ_2^'(t) t)f_0(z) z.
	


	Please see Appendix <ref>.




  §.§.§ Optimal  vs IF Tradeoff


In Section <ref>, we characterize the optimal ψ_1 and ψ_2 that minimize (,F_0,2). As we will see in the following, this choice will lead to an unbounded IF function. Hence, in this section, we aim to minimize (,F_0,2) while putting an upperbound on IF.

In particular, for the location-scale model, the influence function can be written as (see page 233 of <cit.> and using the Fisher consistent constraint)

    IF(z;T,F_0)=([ B_1^-1ψ_1(z); B_2^-1ψ_2(z) ]),

in which
B_1=𝔼{ψ_1^'(Z)} and 
B_2=𝔼{Zψ_2^'(Z)}.

Hence, the unstandardized gross-error sensitivity defined in (<ref>) for the location-scale model is given by

    γ_u^*=sup_z√(B_1^-2ψ_1^2(z)+B_2^-2ψ_2^2(z)).
 
It is easy to check that, if one uses ψ_1(z) and ψ_2(z) characterized in Section <ref>, γ_u^* is unbounded.

In the following, we aim to characterize ψ_1(z) and ψ_2(z) that minimize (,F_0,2) while making sure that γ_u^* is bounded. In particular, we have the following optimization problem

    min      (,F_0,2),
    s.t.      sup_z ||IF(z;T,F_0)||_2^2≤ξ^2,

for any given positive constraint ξ. As mentioned at the beginning of the section, we also require ψ_1 to be an odd function with nonnegative gradient and to be Fisher consistent, and require ψ_2 to be an even function with nonnegative gradient when z≥ 0 and to be Fisher consistent. Again, we only need to focus on the case when z≥ 0.

Plugging the expression of AIF and IF, we have

    min      √(𝔼{ψ_1^'(Z)^2}/(𝔼{ψ_1^'(Z)})^2+𝔼{ψ_2^'(Z)^2}/(𝔼{Zψ_2^'(Z)})^2)
    s.t.      sup_z[ ψ_1^2(z)/(𝔼{ψ_1^'(Z)})^2+ ψ_2^2(z)/(𝔼{Zψ_2^'(Z)})^2]≤ξ^2.


As both ψ_1^'(z) and ψ_2^'(z) are nonnegative, the sup in (<ref>) are achieved when z is either ∞ or 0. Hence, using the constraint that ψ_1(z) is odd, the IF constraint can be written as

    max{ψ_2^2(0)/𝔼{Zψ_2^'(Z)}^2,ψ_2^2(∞)/(𝔼{Zψ_2^'(Z)})^2+ψ_1^2(∞)/(𝔼{ψ_1^'(Z)})^2}≤ξ^2.


By setting ξ_1^2+ξ_2^2=ξ^2, we can first solve the following two problems.

    P1:min      𝔼{ψ_1^'(Z)^2}/(𝔼{ψ_1^'(Z)})^2
    s.t.
          ψ_1^'(z)≥ 0,
          ψ_1^2(∞)/(𝔼{ψ_1^'(Z)})^2≤ξ_1^2,

and

    P2:min      𝔼{ψ_2^'(Z)^2}/(𝔼{Zψ_2^'(Z)})^2
    s.t.      𝔼{ψ_2}=0,
          ψ_2^'(z)≥ 0,
          ψ_2^2(∞)/(𝔼{Zψ_2^'(Z)})^2≤ξ_2^2,
          ψ_2^2(0)/(𝔼{Zψ_2^'(Z)})^2≤ξ^2.


After solving these problems for a given ξ_1^2 and ξ_2^2, we can then adjust the values ξ_1^2 and ξ_2^2 to obtain the overall solution to the optimization problem (<ref>). In P1, we do not write the Fisher consistent constraint, as an odd function ψ_1(z) will automatically satisfy this constraint for symmetric f_0(z).


	The optimal odd function ψ_1(z) is specified by
	
    ψ_1^'(z)=[ν^*-ϑ_1^*/ f_0(z)]^+,

	in which ν^* and ϑ_1^* are chosen to satisfy
	
    ∫_0^∞ψ_1^'(z)f_0(z)dz=1,
          ϑ^*_1( ∫_0^∞ψ_1^'(z)dz
    	-2 ξ_1  )=0.

	
	The optimal even function ψ_2(z) have the following form:
	
    ψ_2^'(z)=[
    	ν^* z-ϑ_2^*+(ϑ_1^*-ϑ_2^*)F_0(z)/ f_0(z)]^+

	and 
	ψ_2(0)=-∫_0^∞ψ_2^'(z)(1-F_0(z)) z, in which
	the parameters ν^*, ϑ^*_1≥ 0 and ϑ^*_2≥ 0 satisfy
	
    ∫_0^∞zψ_2^'(z)f_0(z)dz=1,
          ϑ^*_1( ∫_0^∞ψ_2^'(z)F_0(z)dz
    	- 2ξ_2  )=0,
          ϑ^*_2(∫_0^∞ψ_2^'(z)[1-F_0(z)]dz-2ξ)=0.

	
	


	Please see Appendix <ref>.	




§ NUMERICAL EXAMPLES


In this section, we provide numerical examples to illustrate the results obtained in this paper.


 §.§ Joint Location-Scale Estimation

In this example, we consider the adversarially robust joint location-scale estimation for Laplace random variables. In this case, f_0 has the following form

    f_0(z)=1/2exp(-|z|).


From Proposition <ref>, we know that the optimal  that minimizes AIF(,F_0,2) has the following form:

	
  * ψ_1(0)=0, ψ_1^'(z)=c with c being a constant.
	
  * ψ_2(0)=-1 and ψ_2^'(z)=2z.


In the following, using Theorem <ref>, we characterize the optimal  that minimizes AIF(,F_0,2) subject to the constraint that the ℓ_2 norm of IF(x;T,F_0) is upper-bounded by ξ. We will only state the form of the functions for z≥ 0 as ψ_1(z) is an odd function and ψ_2(z) is an even function.

First, from (<ref>), we know that the optimal form ψ_1^'(z) has the following form: 

    ψ_1^'(z)=ν^*-2ϑ_1^*exp(z), 0≤ z≤ a_1,

in which a_1 is the solution of 
    ν^*/2ϑ_1^*=exp(z),
 and (ν^*,ϑ_1^*) are chosen to satisfy

    ∫_0^a_1(ν^*exp(-z)-2ϑ_1^*)dz=2,
          ∫_0^a_1(ν^*-2ϑ_1^*exp(z))dz
    = 2ξ_1.


It is easy to check that these conditions can be simplified to 

    ν^*-2(a_1+1)ϑ_1^*=2,
           (a_1-1)ν^*+2ϑ_1^*=2ξ_1.

Using these, we can express the values of ν^* and ϑ_1^* in terms of a_1

    ν^*   =   2+2(a_1+1)ξ_1/a_1^2,
    ϑ_1^*   =   ξ_1-(a_1-1)/a_1^2,

and the value of a_1 is determined by

    e^a_1=1+(a_1+1)ξ_1/ξ_1+1-a_1.

For a given ξ_1, the value of a_1 can be obtained by solving (<ref>) numerically. We can then plug the value of a into (<ref>) and (<ref>) to obtain ν^* and ϑ_1^*. 

Secondly, for a given ξ_2, we determine the form of ψ_2 in Theorem <ref>. For f_0 considered in Laplace random variables, we know that 1-F_0(z)≤ F_0(z) when z≥ 0. Using this fact in (<ref>) and (<ref>) along with the fact that ξ_2≤ξ, we know that ϑ_2^*=0 and ϑ_1^*≠ 0. Hence from Theorem <ref>, we know that ψ_2^'(z) has the following form

    ψ_2^'(z)=
    ν^* z-ϑ_1^*(2exp(z)-1),  0≤ a_2≤ z≤ b,

in which both a_2 and b satisfy 
    ν^*/ϑ_1^* =2exp(z)-1/z.

Furthermore ν^* and ϑ_1^* are determined by the following two equations (simplified from (<ref>) and (<ref>))

    ∫_a_2^b(ν^* z-ϑ_1^*(2exp(z)-1))zexp(-z) z=2,
    ∫_a_2^b(ν^* z-ϑ_1^*(2exp(z)-1))(2-exp(-z)) z=4ξ_2.

After simple integral and using the fact that a_2, b satisfy (<ref>), we can simplify the above two equations to

    ν^*/ϑ_1^*=2(e^b-e^a_2)-2(ξ_2+1)(b-a_2)+ξ_2(a_2^2-b^2+e^-b_2-e^-a)/0.5(b^2-a_2^2)+(2ξ_2+0.5)(e^-b-e^-a_2),
    ϑ_1^*[ν^*/ϑ_1^*2(e^-b-e^-a_2)+b^2-a_2^2-e^-b+e^-a_2+2(b-a_2)]=-2.


To find the value of ν^* and ϑ_1^*, we can first obtain ν^*/ϑ_1^* (and hence the values of a_2 and b) numerically from (<ref>) and (<ref>). After knowing ν^*/ϑ_1^*, we can then use (<ref>) to obtain ϑ_1^*. 





Fig <ref> illustrates the tradeoff between AIF and IF of the optimal M-estimator characterized using the approach outlined above. In the figure, we also plot the curve for Huber proposal 2 discussed in Example <ref>. From the figure, we can see that as IF increases (less robust to outliers), AIF decreases (more robust to adversary modifications) and vice verse. Furthermore, there is a gap between the tradeoff achieved by the Huber proposal 2 and the optimal tradeoff achieved by the estimator characterized above. The tradeoff achieved by Huber proposal 2 in turn is better than that of the coupled mean and standard deviation estimator discussed in Example <ref>, for which achieves AIF=√(2) but IF=∞.



 §.§ Robust Linear Regression


In this section, we compare AIF of various robust regression methods discussed in Section <ref>. In the following, we will adopt the commonly used coefficients w_n=√(1-h_nn), v_n=1/w_n and use Huber function for η, that is 
    η(x)=max[-K,min(K,x)].
 To proceed further, we need to compute 

    ∂ w_n/∂_j,k   =   -1/2√(1-h_nn)∂ h_nn/∂_j,k,
    ∂ v_n/∂_j,k   =   1/2(1-h_nn)^3/2∂ h_nn/∂_j,k,
 
both of which depend on ∂ h_nn/∂_j,k. Recall that 
h_nn=_n^T(^T)^-1_n, which is a complicated function of _j,k. To address this, we use _(-j) to denote the data matrix  but with the jth column removed. To simplify notation, we let =^T and _(-j)=_(-j)_(-j)^T. We have =^T=_(-j)_(-j)^T+_j_j^T=_(-j)+_j_j^T. Using Sherman-Morrison formula, we have

    ^-1=(_(-j)+_j_j^T)^-1=_(-j)^-1-_(-j)^-1_j_j^T_(-j)^-1/1+_j^T_(-j)^-1_j.

Hence,

    h_nn=_n^T^-1_n=_n^T_(-j)^-1_n-(_n^T_(-j)^-1_j)^2/1+_j^T_(-j)^-1_j.


As the result, if j=n, we have

    ∂ h_nn/∂_n,k=2_(-n)^-1(k,:)_n/(1+_n^T_(-n)^-1_n)^2.

If j≠ n, we have

    ∂ h_nn/∂_j,k=-2_(-n)^-1(k,:)_n(1+_j^T_(-j)^-1_j)-2(_n^T_(-j)^-1_j)^2_(-n)^-1(k,:)_j/(1+_j^T_(-j)^-1_j)^2.


Plugging (<ref>) and (<ref>) into the corresponding equations of OLS, Huber's proposal, Mallow's proposal and Schewppe's proposal discussed in Section <ref>, we obtain the corresponding AIF.  



Figure <ref> illustrates the comparison of AIF for these methods for different value of K, the parameter in the Huber function (<ref>). In generating this figure, we set q=5, N=500. We first generate  using Gaussian zero mean and variance 1. After  is generated, it is fixed throughout the simulation. We let K to be from 3.05 to 5. For each position, we run 100 times and obtain the average. We generate each entry of  using i.i.d with zero mean and variance 1. We then obtain y_n by adding zero mean variance 2 noise to ^T_n. From the figure, we can see that for the same value of K, Mallow's proposal has the smallest value of AIF (i.e., it is the most robust again adversary modifications), the AIF value of Scheweppe's proposal is similar to Huber's proposal. Furthermore, all three approaches are more adversary modification resistant than OLS, which does not depend on K. Furthermore, as K increases, the performance of all three methods approach that of OLS. This is expected, as K increases, all three approaches are similar to OLS.



§ CONCLUSION


In this paper, we have investigated the adversarial robustness of multivariate M-Estimators. We have characterized the adversary's optimal modification strategy and its corresponding AIF. Under certain regularity conditions, we have characterized the optimal M-estimator for the case of joint location-scale estimation problem. We have further identified a tradeoff between robustness against adversarial modification and robustness against outliers, and have derived the optimal M-estimator that achieves the best tradeoff.





§ PROOF OF F^*=G^*

First, we show f^*≥ g^*.

To simplify the presentation, we let _i, i=1,⋯,q be the i-th row of _N^'(). For any ∈, we have

    ||_N^'()Vec(Δ)||_1=∑_i=1^q|_i Vec(Δ)|
    (a)=∑_i=1^q|(i)_i Vec(Δ)|
    ≥∑_i=1^q (i)_i Vec(Δ),
	
in which (a) is true as each entry of  is either 1 or -1.
As this holds for any ∈, hence f^*≥ g^*.

Next we show f^*≤ g^*.

||_N^'()Vec(Δ)||_1=∑_i=1^q|_i Vec(Δ)|=∑_i=1^q sign(_i Vec(Δ))_i Vec(Δ). Hence


    f^*   =   max_Vec(||Δ)||_p^p≤ Nmδ^p||_N^'()Vec(Δ)||_1
       =   max_∈max_sign(_i Vec(Δ))=(i),||Δ)||_p^p≤ Nmδ^p∑_i=1^q (i)_i Vec(Δ)
       ≤   max_∈max_||Δ)||_p^p≤ Nmδ^p∑_i=1^q (i)_i Vec(Δ)=g^*.




§ PROOF OF PROPOSITION <REF>

From (<ref>), we have the following variational optimization problem

    min      ∫_0^∞ψ_2^'(z)^2f_0(z)  z/(∫_0^∞zψ_2^'(z)f_0(z) z)^2,
    s.t.      ∫_0^∞ψ_2(z)f_0(z) z=ψ_2(0)+∫_0^∞(∫_0^zψ_2^'(t) t)f_0(z) z=0,
          ψ_2^'≥ 0.


As ψ_2(0) does not appear in the objective function, we can solve (<ref>) without the constraint (<ref>) first. After that, we can simply set 
    ψ_2(0)=-∫_0^∞(∫_0^zψ_2^'(t) t)f_0(z) z
 so that the constraint (<ref>) will be satisfied. 

In the following, to simplify the notation, we will use g(z) to denote ψ_2^'(z). It is clear that the optimization problem is scale invariant in the sense that if g^*(z) is a solution to this problem, then for any positive constant c,  cg^*(z) is also a solution to this problem. As a result, without loss of generality, we can assume ∫_0^∞zg(z)f_0(z)dz=1. Using this, we can further simplify the optimization problem to 

    min      1/2∫_0^∞g^2(z)f_0(z)dz,
    s.t.      ∫_0^∞zg(z)f_0(z)dz=1,
           g(z)≥ 0.

For this convex calculus of variations problem, we form Lagrange function

    ℒ=∫_0^∞1/2g^2(z)f_0(z)dz+ν(-∫_0^∞zg(z)f_0(z)dz+1)-∫_0^∞λ(z)g(z) z.







For any given z, the optimal value of g(z) can be found from

    g^*(z)f_0(z)-ν^* zf_0(z)-λ^*(z)=0,

in which the parameters ν^* and λ^*(z)≥ 0 satisfy <cit.>

    ∫_0^∞g^*(z)f_0(z)dz=1,
          λ^*(z)g(z)= 0.


From (<ref>), for z≥ 0 in the range of f_0(z), we have

    g^*(z)=λ^*(z)+ν^* zf_0(z)/f_0(z).



As z≥ 0 and f_0≥ 0, we know from (<ref>) that λ^*(z)=0, and hence

    g^*(z)=ν^* z.

and the optimal value of ν^* is 

    ν^*=1/∫_0^∞z^2f_0(z)dz.

As the result, for z in the range of f_0(z), the optimal g^*(z) is 

    g^*(z)=z/∫_0^∞z^2f_0(z)dz,

and 
ψ_2(0)=-∫_0^∞(∫_0^zg^*(t) t)f_0(z) z.



§ PROOF OF THEOREM <REF>


We first focus on P1, and rewrite P1 into the following form

    P1a:min      𝔼{ψ_1^'(Z)^2}/(𝔼{ψ_1^'(Z)})^2=∫_0^∞ψ_1^'(z)^2f_0(z) z/2(∫_0^∞ψ_1^'(z)f_0(z) z)^2, 
    s.t.
          ψ_1^'(z)≥ 0,
          ψ_1^2(∞)/(𝔼{ψ_1^'(Z)})^2=(∫_0^∞ψ_1^'(z) x)^2/4(∫_0^∞ψ_1^'(z)f_0(z) z)^2≤ξ_1^2.

Here, we use the constraint that ψ_1(z) is an odd function.

To ease the notation, in the following, we use h(z) to denote ψ_1^'(z). It is clear that the optimization problem is scale invariant, hence we can assume ∫_0^∞h(z)f_0(z) z=1. Hence, P1a can be simplified to

    min      1/2∫_0^∞h^2(z)f_0(z) z, 
    s.t.       h(z)≥ 0,
          ∫_0^∞h(z) z≤ 2ξ_1,
          ∫_0^∞h(z)f_0(z) z=1.


To solve this convex functional minimization problem, we first form the Lagrangian function

    ℒ   =   1/2∫_0^∞h^2(z)f_0(z)dz+ν(-∫_0^∞h(z)f_0(z)dz+1)-∫_0^∞λ(z)h(z)dz
          +ϑ_1(∫_0^∞h(z)dz- 2ξ_1).



For any given z, the optimal value of h(z) can be found from

    h^*(z)f_0(z)-ν^* f_0(z)+ϑ_1^*-λ^*(z)=0,

in which the parameters ν^*, ϑ_1^*≥ 0 and λ^*(z)≥ 0 satisfy <cit.>

    ∫_0^∞h^*(z)f_0(z)dz=1,
          ϑ^*_1( ∫_0^∞h^*(z)dz
    -2 ξ_1  )=0,
          λ^*(z)h(z)= 0.


From (<ref>), for z≥ 0 in the range of f_0(z), we have

    h^*(z)=λ^*(z)+ν^* f_0(z)-ϑ^*_1/f_0(z).


Combining this with condition (<ref>), we know that if ν^* f_0(z)-ϑ^*_1>0, then λ^*(z)=0. On the other hand, if ν^* f_0(z)-ϑ^*_1<0, then h^*(z)=0. As a result, we have

    h^*(z)={[  ν^*-ϑ_1^*/ f_0(z), ν^* f_0(z)> ϑ_1^* ;;                  0,          otherwise. ].

This characterizes the optimal odd function ψ_1(z).


We now focus on ψ_2. We rewrite P2 into

    P2a:min      𝔼{ψ_2^'(Z)^2}/(𝔼{Zψ_2^'(Z)})^2=∫_0^∞ψ_2^'(z)^2f_0(z) z/2(∫_0^∞zψ_2^'(z)f_0(z) z)^2,
    s.t.      𝔼{ψ_2}=2ψ_2(0)+2∫_0^∞[∫_0^zψ_2^'(t) t]f_0(z) z=0,
          ψ_2^'(z)≥ 0,
          ψ_2^2(∞)/(𝔼{Zψ_2^'(Z)})^2=(ψ_2(0)+∫_0^∞ψ_2^'(z) z)^2/4(∫_0^∞zψ_2^'(z)f_0(z) z)^2≤ξ_2^2,
          ψ_2^2(0)/(𝔼{Zψ_2^'(Z)})^2=ψ_2^2(0)/4(∫_0^∞zψ_2^'(z)f_0(z) z)^2≤ξ^2.


We note that in P2a,

    ∫_0^∞[∫_0^zψ_2^'(t) t]f_0(z) z=∫_0^∞ψ_2^'(z)(1-F_0(z)) z.

To satisfy the Fisher consistent constraint (<ref>) in P2a, we need to set 
    ψ_2(0)=-∫_0^∞ψ_2^'(z)(1-F_0(z)) z.

Using these, P2a can be simplified to

    P2b:min      𝔼{ψ_2^'(Z)^2}/(𝔼{Zψ_2^'(Z)})^2=∫_0^∞ψ_2^'(z)^2f_0(z) z/2(∫_0^∞zψ_2^'(z)f_0(z) z)^2, 
    s.t.      ψ_2^'(z)≥ 0,
          (-∫_0^∞ψ_2^'(z)(1-F_0(z)) z+∫_0^∞ψ_2^'(z) z)^2/4(∫_0^∞zψ_2^'(z)f_0(z) z)^2≤ξ_2^2,
          (-∫_0^∞ψ_2^'(z)(1-F_0(z)) z)^2/4(∫_0^∞zψ_2^'(z)f_0(z) z)^2≤ξ^2.


Similar to other cases, it is clear that P2b is scale invariant, and hence we can without loss of generality assume that ∫_0^∞zψ_2^'(z)f_0(z)dz=1. Using this fact and denoting g(z)=ψ_2^'(z), P2b can be simplified to

    P2b:min      1/2∫_0^∞g(z)^2f_0(z) z 
    s.t.      ∫_0^∞g(z)F_0(z) z≤ 2ξ_2,
          ∫_0^∞g(z)(1-F_0(z)) z≤ 2ξ,
          ∫_0^∞zg(z)f_0(z) z=1,
           g(z)≥ 0.

To solve this convex functional minimization problem, we first form the Lagrangian function

    ℒ   =   1/2∫_0^∞g^2(z)f_0(z)dz+ν(-∫_0^∞zg(z)f_0(z)dz+1)-∫_0^∞λ(z)g(z)dz
          +ϑ_1(∫_0^∞g(z)F_0(z)dz
    - 2ξ_2 )+ϑ_2(∫_0^∞g(z)[1-F_0(z)]dz- 2ξ).






For any given z, the optimal value of g(z) can be found from

    g^*(z)f_0(z)-ν^* z f_0(z)+ϑ^*_2+(ϑ^*_1-ϑ^*_2)F_0(z)-λ^*(z)=0,

in which the parameters ν^*, ϑ^*_1≥ 0, ϑ^*_2≥ 0, λ^*(z)≥ 0 satisfy

    ∫_0^∞xg^*(z)f_0(z)dz=1,
          ϑ^*_1( ∫_0^∞g^*(z)F_0(z)dz
    - 2ξ_2  )=0,
          ϑ^*_2(∫_0^∞g^*(z)[1-F_0(z)]dz- 2ξ)=0,
          λ^*(z)g^*(z)= 0.


From (<ref>), for those z≥ 0 with f_0(z)>0, we have

    g^*(z)=λ^*(z)+ν^* z f_0(z)-ϑ^*_2-(ϑ^*_1-ϑ^*_2)F_0(z)/f_0(z).


Combining this with the condition (<ref>), we know that if ν^* z f_0(z)-ϑ^*_2-(ϑ^*_1-ϑ^*_2)F_0(z)>0, then λ^*(z)=0. On the other hand, if ν^* f_0(z)-ϑ^*_2-(ϑ^*_1-ϑ^*_2)F_0(z)<0, then g^*(z)=0. As the result, we have


    g^*(z)={[    ν^* z-ϑ_2^*+(ϑ_1^*-ϑ_2^*)F_0(z)/ f_0(z), ν^* z f_0(z)> ϑ_1^* F_0(z)+ϑ_2^*(1-F_0(z));;                                          0,                                  otherwise. ].

Coupled with ψ_2(0)=-∫_0^∞g^*(z)(1-F_0(z)) z, this characterizes the optimal even ψ_2(z).



§ ACKNOWLEDGEMENTS

The work of E. Bayraktar was supported in part by the National Science Foundation under grant DMS-1613170 and by the Susan M. Smith Professorship. The work of L. Lai was supported by the National Science Foundation under grants CCF-17-17943 and ECCS-17-11468. 





28
	
	
	urlstyle[Buczak and Guven(2016)]mlcyber 
	A. L. Buczak and E. Guven. (2016)
	A survey of data mining and machine learning methods for cyber
	security intrusion detection.
	IEEE Communications Surveys Tutorials, 180
	(2):0 1153–1176.
	
	
	[Carlini and Wagner(2017)]Carlini:ISSP:17
	Nicholas Carlini and David Wagner. (2017)
	Towards evaluating the robustness of neural networks.
	In Proc. IEEE Intl. Symposium on Security and Privacy, San
	Jose, CA.
	
	[Ferguson(1996)]Ferguson:Book:96
	T. Ferguson. (1996)
	A Course in Large Sample Theory.
	Chapman and Hall, London, UK.
	
	[Goodfellow et al.(2015)Goodfellow, Shlens, and
	Szegedy]Goodfellow:ICLR:15
	Ian J. Goodfellow, Jonathon Shlens, and Christian Szegedy. (2015)
	Explaining and harnessing adversarial examples.
	In Proc. International Conference on Learning Representations,
	San Diego, CA.
	
	[Gregory and Lin(1992)]Gregory:Book:92
	J. Gregory and C. Lin. (1992)
	Constrained Optimization in the Calculus of Variations and
		Optimal Control Theory.
	Van Nostrand Reinhold, New York, NY.
	
	[Hampel(1968)]Hampel:PHD:68
	F. Hampel. (1968)
	Contributions to the theory of robust estimation.
	Ph.D. thesis, University of California, Berkeley.
	
	[Hampel et al.(2009)Hampel, Ronchetti, Rousseeuw, and
	Stahel]Hampel:Book:86
	Frank Hampel, E. Ronchetti, Peter Rousseeuw, and W. Stahel. (2009)
	Robust statistics: The Approach Based on Influence Functions.
	Wiley.
	
	[Hoermann et al.(2018)Hoermann, Bach, and Dietmayer]8460874
	S. Hoermann, M. Bach, and K. Dietmayer. (2018)
	Dynamic occupancy grid prediction for urban autonomous driving: A
	deep learning approach with fully automatic labeling.
	In Proc. IEEE International Conference on Robotics and
		Automation, pages 2056–2063, Brisbane, Australia.
	
	[Huber(1964)]Huber:AMS:64
	P. J. Huber. (1964)
	Robust estimation of a location parameter.
	Ann. Math. Statist., 35:0 73–101.
	
	[Huber(1973)]Huber:AoS:73
	P. J. Huber. (1973)
	Robust regression: Asymptotics, conjectures, and Monte Carlo.
	Ann. Statist., pages 799–821.
	
	[Huber and Ronchetti(2009)]Huber:Book:09
	Peter Huber and E. Ronchetti. (2009)
	Robust statistics.
	Wiley.
	
	[Hubert et al.(2008)Hubert, Rousseeuw, and Van Aelst]hubert2008
	Mia Hubert, Peter J. Rousseeuw, and Stefan Van Aelst. (2008)
	High-breakdown robust multivariate methods.
	Statist. Sci., 230 (1):0 92–119.
	
	[Huval et al.(2015)Huval, Wang, Tandon, Kiske, Song, Pazhayampallil,
	Andriluka, Rajpurkar, Migimatsu, Cheng-Yue, et al.]huval2015empirical
	B. Huval, T. Wang, S. Tandon, J. Kiske, W. Song, J. Pazhayampallil,
	M. Andriluka, P. Rajpurkar, T. Migimatsu, R. Cheng-Yue, et al. (2015)
	An empirical evaluation of deep learning on highway driving.
	arXiv preprint arXiv:1504.01716.
	
	[Jagielski et al.(2018)Jagielski, Oprea, Biggio, Liu, Nita-Rotaru, and
	Li]Jagielski:ISSP:18
	Matthew Jagielski, Alina Oprea, Battista Biggio, Chang Liu, Cristina
	Nita-Rotaru, and Bo Li. (2018)
	Manipulating machine learning: Poisoning attacks and
	countermeasures for regression learning.
	In Proc. IEEE Intl. Symposium on Security and Privacy, San
	Francisco, CA.
	
	[Kot(2014)]Kot:Book:14
	Mark Kot. (2014)
	A first course in the calculus of variations.
	American Mathematical Society, Providence, RI.
	
	[Lai and Bayraktar(2018)]Lai:TIT:18
	L. Lai and E. Bayraktar. (2018)
	On the adversarial robustness of robust estimators.
	IEEE Trans. Inform. Theory.
	Submitted.
	
	[Litjens et al.(2016)Litjens, Sánchez, Timofeeva, Hermsen,
	Nagtegaal, Kovacs, Hulsbergen-Van De Kaa, Bult, Van Ginneken, and Van
	Der Laak]litjens2016deep
	G. Litjens, C. I. Sánchez, N. Timofeeva, M. Hermsen, I. Nagtegaal,
	I. Kovacs, C. Hulsbergen-Van De Kaa, P. Bult, B. Van Ginneken, and J. Van
	Der Laak. (2016)
	Deep learning as a tool for increased accuracy and efficiency of
	histopathological diagnosis.
	Scientific reports, 60 (1):0 26286.
	
	[Mallows(1975)]Mallows:BTL:75
	C. L. Mallows. (1975)
	On some topics in robustness.
	Technical Memorandum.
	Bell Telephone Laboratories, Murray Hill.
	
	[Merrill and Schweppe(1971)]Merrill:TP:71
	H. Merrill and F. Schweppe. (1971)
	Bad data suppression in power system static state estimation.
	IEEE Trans. Power App. Syst., pages 2718–2725.
	
	[Nelson et al.(2008)Nelson, Barreno, Chi, Joseph, Rubinstein, Saini,
	Sutton, Tygar, and Xia]Nelson:2008:EML
	Blaine Nelson, Marco Barreno, Fuching Jack Chi, Anthony D. Joseph, Benjamin
	Rubinstein, Udam Saini, Charles Sutton, J. D. Tygar, and Kai Xia. (2008)
	Exploiting machine learning to subvert your spam filter.
	In Proc. Usenix Workshop on Large-Scale Exploits and Emergent
		Threats, pages 7:1–7:9, San Francisco, CA.
	
	[Nguyen et al.(2015)Nguyen, Yosinski, and Clune]easyfooldnn
	A. Nguyen, J. Yosinski, and J. Clune. (2015)
	Deep neural networks are easily fooled: High confidence predictions
	for unrecognizable images.
	In Proc. IEEE Conference on Computer Vision and Pattern
		Recognition, pages 427–436, Boston, MA.
	
	[Pimentel-Alarcon et al.(2017)Pimentel-Alarcon, Biswas, and
	Solis-Lemus]Pimentel:ISIT:17
	Daniel L. Pimentel-Alarcon, Aritra Biswas, and Claudia R. Solis-Lemus. (2017)
	Adversarial principal component analysis.
	pages 2363–2367.
	
	[Soule et al.(2005)Soule, Salamatian, and Taft]Soule:2005:CFS
	Augustin Soule, Kavé Salamatian, and Nina Taft. (2005)
	Combining filtering and statistical methods for anomaly detection.
	In Proc. ACM SIGCOMM Conference on Internet Measurement, pages
	31–31, Berkeley, CA.
	
	[Stamp(2018)]Stamp:Book:18
	M. Stamp. (2018)
	Introduction to Machine Learning with Applications in
		Information Security.
	New York: Chapman and Hall/CRC.
	
	[Suthaharan(2014)]Suthaharan:2014:BDC:2627534.2627557
	Shan Suthaharan. (2014)
	Big data classification: Problems and challenges in network intrusion
	prediction with machine learning.
	SIGMETRICS Perform. Eval. Rev., 410 (4):0
	70–73.
	
	[Szegedy et al.(2013)Szegedy, Zaremba, Sutskever, Bruna, Erhan,
	Goodfellow, and Fergus]szegedy2013intriguing
	C. Szegedy, W. Zaremba, I. Sutskever, J. Bruna, D. Erhan, I. Goodfellow, and
	R. Fergus. (2013)
	Intriguing properties of neural networks.
	arXiv preprint arXiv:1312.6199.
	
	[Vallon et al.(2017)Vallon, Ercan, Carvalho, and Borrelli]7995936
	C. Vallon, Z. Ercan, A. Carvalho, and F. Borrelli. (2017)
	A machine learning approach for personalized autonomous lane change
	initiation and control.
	In Proc. IEEE Intelligent Vehicles Symposium, pages
	1590–1595, Redondo Beach, CA.
	
	[Wilcox(2005)]Wilcox:Book:05
	R. Wilcox. (2005)
	Introduction to Robust Estimation and Hypothesis Testing.
	Elsevier Academic Press, Burlington, MA.
	


