Error bounds of Median-of-means estimators with VC-dimension
Yuxuan Wang Yiming Chen Hanchao Wang Lixin Zhang


Article Title]Error bounds of Median-of-means estimators with VC-dimension










1]Yuxuan Wangwangyuxuan@zju.edu.cn

2]Yiming Chenchenyiming960212@mail.sdu.edu.cn


2]Hanchao Wangwanghanchao@sdu.edu.cn


[1,3]Lixin Zhangstazlx@zju.edu.cn


*[1]School of Mathematical Sciences, Zhejiang University, Hangzhou, 310027, China

[2]Institute for Financial Studies, Shandong University, Jinan, 250100, China

[3]School of Statistics and Mathematics, Zhejiang Gongshang University,  Hangzhou, 310018, China





We obtain the upper error bounds of robust estimators for mean vector, using the median-of-means (MOM) method. The method is designed to handle data with heavy tails and contamination, with only a finite second moment, which is weaker than many others, relying on the VC dimension rather than the Rademacher complexity to measure statistical complexity. This allows us to implement MOM in covariance estimation, without imposing conditions such as L-sub-Gaussian or L_4-L_2 norm equivalence. In particular, we derive a new robust estimator, the MOM version of the halfspace depth, along with error bounds for mean estimation in any norm.



















[
    *
    Received Month dd, yyyy; accepted Month dd, yyyy
====================================================





§ INTRODUCTION


Inspired by applications in machine learning and data science, there has been a growing interest in constructing μ mean estimators in recent years. As the most basic method of estimation, the sample mean μ̅_N=1/N∑_i=1^n X_i on a sample (X_1, …, X_N) of N independent and identically distributed random variables possesses favorable statistical properties established by the central limit theorem. However, the asymptotic properties often require a large sample size in practical applications, significantly increasing the difficulty level. Alternatively, non-asymptotic estimators with faster convergence rates have emerged, implying the need for fewer samples.

Simultaneously, in situations where the distribution exhibits heavy-tailed characteristics or outliers present in the data, the empirical mean may no longer be sufficient to meet the requirements. There is an urgent need to enhance the quality of mean estimation, particularly in a non-asymptotic context. It is worth mentioning that in mean estimation, statistical optimality, including convergence rates and computational complexity, are crucial factors. We will primarily focus on estimators that maintain high accuracy while providing substantial confidence in mean estimation.




Unlike asymptotic estimators, we consider a non-asymptotic estimator known as the "L-sub-Gaussian" estimator as below. For a mean estimator μ_N and any δ∈(0,1), let σ^2 be the variance, then there exists a constant L > 0, for any sufficiently large sample size N, the following inequality holds with at least 1 - δ probability:

    |μ_N-μ| ≤L σ√(log (2 / δ))/√(N).

This definition comes from the accuracy of empirical mean estimation in the context of sub-Gaussian distributions. Many well-known robust mean estimators exhibit this property, with the most common being the median-of-means (MOM) estimator. MOM method has experienced rapid development, as evidenced by works such as  <cit.>, <cit.> and <cit.>.

The estimation of vector means and real-valued sample means are fundamentally different, as is illustrated in <cit.>. The former involves transforming into a problem of concentration inequalities for the upper bounds of a stochastic process indexed by vectors in ℝ^d. This represents an essential distinction and also marks a departure from previous work.

For the estimation of the mean in the multivariate sub-Gaussian case, the empirical mean of independent and identically distributed samples with mean μ and covariance matrix Σ satisfies, with at least 1-δ probability:


    ∑_i=1^NX_iN-μ≤√(Tr(Σ)N)+ √(Σlog(1/δ)N),

where Tr(Σ) is the trace of Σ and Σ is the operator norm of Σ. A mean estimator is considered sub-Gaussian, as defined in <cit.>, if it satisfies an inequality of the form above (with possibly different constant factors). <cit.> used the MOM estimator and obtained near-optimal confidence bounds for mean estimation under general heavy-tailed conditions, specifically when only the second moment is finite. This means that for δ∈ (0,1), the optimal confidence upper bound holds with at least a probability of 1-δ:
    μ_N-μ≤c/√(N)(max{𝔼Y_N, 𝔼G+R √(log (2 / δ))}),

where c is an absolute constant and

    R=sup _x^* ∈ℬ^∘(𝔼(x^*(X-μ))^2)^1 / 2,


    Y_N=1/√(N)∑_i=1^N ε_i(X_i-μ),

where ℬ^∘ is the unit ball of the dual space to (ℝ^d,), (ε_i)_i=1^N are i.i.d. symmetric {-1,1}-valued random variables that are also independent of (X_i)_i=1^N, and 𝔼Y_N/√(N) is called the Rademacher complexity. Also, observe that by the central limit theorem, Y_N tends, in distribution, to the centered Gaussian random vector G that has the same covariance as X.

 <cit.> furthered this line of research by constructing estimators for the covariance matrix of random vectors that are robust to heavy tails and outliers. However, their approach imposes stricter conditions, such as the L_4-L_2 norm equivalence.
Some other covariance estimators (<cit.>) adopt a two-step approach: firstly, estimating μ by MOM or other methods, and then employing truncation techniques to estimate the covariance, aiming to mitigate the influence of heavy-tails. This two-step method is often used to improve the robustness and accuracy of covariance matrix estimation in the presence of heavy-tailed distributions or outliers.


    As one of the earliest studies on sub-Gaussian mean estimator, <cit.> introduced a sharp example for distributions with known variances and distributions with finite fourth moments and known kurtosis upper bounds. Utilizing a specific function ϕ, the M-estimator, also known as the Catoni estimator, provides a well-performed confidence interval, naturally extended to mean estimation of heavy-tailed random vectors (<cit.>) and covariance matrix estimation (<cit.>). However, the construction of Catoni estimator inevitably requires the distribution variance σ^2 or Σ to be known.  
 
In computer science and logic design, Boolean functions are basic, representing functions with binary outputs. In machine learning, they are used to model simple classification problems. The VC dimension, first introduced by Vapnik and Chervonenkis (<cit.>), measures the maximum complexity a model's hypothesis space can handle, especially in classification tasks. While it is primarily used for assessing classifier complexity, directly computing the VC dimension can be challenging, especially in high-dimensional and complex models. It serves as a theoretical guide for understanding a model's learning and generalization abilities. 

<cit.> introduced a novel general approach to constrain the estimation error of MOM estimators. The author applied VC dimension instead of Rademacher complexity to measure statistical complexity, which does not take into account
the unknown structure of the covariance matrix, but is related only to the dimension
of the dual space. 

In the context of multivariate analysis, the varying definitions of the median give rise to distinct MOM estimators. Among those, geometric median is proved, by <cit.>, to be available in constructing robust MOM estimators.
 <cit.> then discussed the construction of sub-Gaussian estimators of a mean vector by VC dimension, using the MOM versions of the Stahel-Donoho outlyingness (SDO) and Median Absolute Deviation (MAD) functions. 
 
Inspired by the groundbreaking work of Depersin, our objective is to establish the median-of-means (MOM) procedure to attain a sub-Gaussian rate in estimating mean vectors. This MOM approach holds the promise of being implemented for covariance estimation directly, bypassing the need for L-sub-Gaussian or L_4-L_2 norm equivalence conditions. Moreover, we aim to extend its utility to practical applications such as Principal Component Analysis (PCA).

Furthermore, ever since <cit.> introduced the concept of data depth (also known as halfspace depth), it has emerged as a fundamental tool for assessing the centrality of data points in multivariate datasets. Consequently, in this paper, we endeavor to delve into the error bounds achievable by a novel estimator: the MOM adaptation of Tukey's median. We anticipate exploring the bounds within certain confidence levels, leveraging the VC dimension as a guiding framework. Through this exploration, we aim to shed light on the robustness and efficacy of MOM estimators in practical data analysis scenarios.


The structure of this paper is as follows: In Section <ref>, we will provide necessary symbol explanations and introduce the definitions and lemmas; in Section <ref> and Section <ref>, we  present error bounds of mean estimation introduced by <cit.> and covariance estimation, respectively; finally, in Section <ref>, we will give a Tukey MOM estimator and discuss the applications of this technique in regression problems.



§ PRELIMINARY



 §.§ Notation

In this paper, we assume that the covariance matrix of interest is non-degenerate. We use ||·|| to represent a norm on ℝ^d, and assume the existence of an inner product ⟨·, ·⟩ that induces this norm. ||·||_* denotes its dual norm. ℬ represents the unit ball of the norm ||·||, and ℬ^* represents the unit ball of the norm ||·||_*. ℬ_0^* is defined as the set of extreme vectors of ℬ^*. We also introduce an operator norm ||A||=sup_u ∈ℬ^* ||Au||_2, where ||·||_2 is the Euclidean norm on ℝ^d. In particular, for a vector u = (u_i), ||u||_2 = √(∑_i u_i^2) represents the ℓ^2 norm. The set S^d-1 = {u ∈ℝ^d : ||u|| = 1} is the unit sphere in ℝ^d. For a matrix A = (A_ij), when A = A^T ∈ℝ^p× p is symmetric, we use λ_j(A) to denote its j-th largest singular value. The operator norm of A is represented as ||A||_op = λ_1(A), and the Frobenius norm is denoted as ||A||_F = √(∑_ij A^2_ij). 

Given an integer d and a, b ∈ℝ, we use [d] to denote the set {1, 2, …, d} and write a ∨ b = max(a, b) and a∧ b = min(a, b). For two non-negative sequences {a_n}, {b_n}, for some constant C > 0 independent of n, a_n≲ b_n means a_n≤ C b_n, and a_n≳ b_n means a_n≥ C b_n.  Throughout the entire paper, C, c and their variations, whose specific values may vary, represent universal constants independent of n. Additionally, the indicator function 1_B(·) is defined as
    1_B(u)= 1,     if  u∈ B, 
     0,     if  u ∉ B .





 §.§ Median of mean

Recall the definition of the classic median-of-means (MOM). First, we randomly divide the data into K equally sized blocks B_1, …, B_K (if K does not divide N evenly, we discard some data). Then, we calculate the empirical mean within each block. For k=1, …, K,

    X̅_k=1/m∑_i ∈ B_k X_i.


In the one-dimensional case, for x_1,…,x_n∈ℝ, Med(x_k)=x_i, such that

    #{j ∈[n]: x_j ≤ x_i}≥n/2 and #{j ∈[n]: x_j ≥ x_i}≥n/2,

where #(·) denotes the cardinality of the set, and if there are multiple i satisfying this condition, the median is defined as the smallest among them.
Let the MOM estimator be μ̃_0Med(X̅_k), then it can be shown that under suitable second-moment conditions, μ̃_0 is a sub-Gaussian estimator.

Now, let Y_1, …, Y_N denote N independent and identically distributed random vectors in ℝ^d. Our goal is to estimate 𝔼Y_1 = μ, assuming that Y_1 has finite second moments. Define Σ = 𝔼((Y_1-μ)(Y_1-μ)^T), sometimes also denoted as 𝔼((Y_1-μ) ⊗(Y_1-μ)), to represent the unknown covariance matrix of Y_1.



 §.§ Contamination model

As in practice, sometimes, we cannot directly observe the vectors Y_1, …, Y_N. Instead, this dataset may already be contaminated or corrupted. One of the most famous examples is the so-called Huber's contamination model. In this setting, instead of observing samples directly from the true distribution P, we observe samples drawn from P_ε, which for an arbitrary distribution Q is deﬁned as a mixture model,
    P_ε = (1-ε)P +ε Q.
   This setting is called the ε-contamination model, first proposed in a path-breaking paper by <cit.>.

More generally, our problem is that the contamination may be adversarial (<cit.>). This means that when an ε fraction of all observed values is maliciously tampered with by an adversary, who is aware of both the "clean" samples and our estimators, there exists a (possibly random) set 𝒪 such that for any i ∈𝒪^c, X_i=Y_i. Here, the size of 𝒪 satisfies |𝒪| ≤⌊ε N ⌋. Thus, the dataset we observe is {X_i: i=1, …, N}, and this model is commonly referred to as a strong contamination model. The contaminated samples {X_i: i=1, …, N} will be called ε-contaminated samples. Furthermore, our task is to recover μ and Σ.



 §.§ VC dimension

Boolean classes ℱ arise in the problem of classification, where ℱ can be taken to consist of all functions
f of the form 1_{g(X) Y} for mappings g. VC dimension was first studied by Vapnik and Cervonenkis
in the 1970s, and let us recall the classical definitions.

	Every f ∈ℱ, taking values in {0, 1}, is called a boolean function. And ℱ is called a boolean class of functions.



	Let 𝒞 be a class of subsets of any set 𝒳. We say that 𝒞 picks out a certain subset from {x_1, . . . , x_n} if this can be formed as a set of the form C ∩{x_1, . . . , x_n} for some C∈𝒞. The collection 𝒞 is said to shatter {x_1, . . . , x_n} if each of its 2^n subsets can be picked out by 𝒞. The VC dimension  VC(𝒞) is the largest cardinality of a set shattered by 𝒞, more formally, 
    VC (𝒞) = sup{ n : max_x_1 ,...,x_n∈𝒳#{C∩{x_1, . . . , x_n}:C∈𝒞}= 2^n} ,
 and in particular, VC(𝒞)=-1 if 𝒞 is empty.

The definition of VC dimension can be easily extended to a function class ℱ in which every
function f is binary-valued, taking the values within {0, 1}.  In this case, we define
    VC (ℱ) = sup{ n : max_x_1, . . . , x_n∈𝒳#{( f(x_1) ,...,f(x_n)) :f∈ℱ}= 2^n}.

In particular, we derive the equivalent definition of the VC dimension as for a set C be a subset of Euclidean space E, the VC dimension of the set of half-spaces generated by the vectors of C,
    VC(C)=VC( { x∈ E→1_⟨ x,v⟩≥ 0: v∈ C}) .

There are some basic facts about VC dimension from Section 7 in <cit.>.
	
		
  * VC(ℝ^d)=d+1. If F is a set of real-valued functions in a k-dimensional linear space, then Pos(F):={x →1_f(x) ≥ 0, f ∈ F} has VC dimension k+1 .
		
		
  * For a function g: 𝒴→𝒳, if we note ℱ∘ g={f ∘ g | f ∈ℱ}, then we have VC(ℱ∘ g) ≤VC(ℱ).
		
		
  * For any r>0, VC({x ∈ E →1_⟨ x, v⟩≥ r, v ∈ C}) ≤VC(C-C) ≲VC(C). 
	


The following lemma can be found in <cit.>:

	If VC(𝒞_i)=V_i , i=1,…,m, 
	let V ≡∑_j=1^m V_j, and let
	
    ⊔_j=1^m 𝒞_j    ≡{∪_j=1^m C_j: C_j ∈𝒞_j, j=1, …, m}, 
    ⊓_j=1^m 𝒞_j    ≡{∩_j=1^m C_j: C_j ∈𝒞_j, j=1, …, m},
    ⊠_j=1^m 𝒞_j    ≡{C_1 ×…× C_m: C_j ∈𝒞_j, j=1, …, m}.

	Then the following bounds hold:
	
    {[ V(⊔_j=1^m 𝒞_j); V(⊓_j=1^m 𝒞_j);   V(⊠_1^m 𝒞_j) ]}≤ c_1 V log(c_2 m),

	where c_1 = e/(e-1) log (2)≈ 2.28231.

The following is a classical result of Vapnik-Chervonenkis theory (one can see more details in <cit.>), which shows the connection between the VC dimension and empirical processes.

	Let ℱ be a class of Boolean functions on a probability space (Ω, 𝒜, ℙ) with finite VC dimension VC(ℱ) ≥ 1. Let X, X_1, X_2, …, X_N be independent random points in Ω distributed according to the law ℙ. Then
	
    𝔼sup _f ∈ℱ|1/N∑_i=1^N f(X_i)-𝔼 f(X)| ≤ C √(VC(ℱ)/N) .






§ MEAN ESTIMATION 


For the mean estimation of a multi-dimensional random vector, we have the following class of median-of-means (MOM) estimators:

For ϵ>0, the sample (X_i)_i=1^N can be partitioned into K blocks B_k, each of size m=N / K. Let X̅_k=1/m∑_i ∈ B_k X_i.
For each x^* ∈ℬ^*, we obtain the set

    S_x^*={y ∈ℝ^d: | Med(x^*(X̅_k) : k∈ [K])-x^*(y) | ≤ϵ} .



Let 𝕊(ϵ)=⋂_x^* ∈ℬ_0^* S_x^*, and μ_K(ϵ, δ) can be taken as any point in 𝕊(ϵ).
It can be shown that the proposed estimator satisfies the following.

	For any δ∈[e^-c N, 1 / 2], there exists an estimator μ_δ such that, with probability at least 1-δ,	
    μ_δ-μ≲ R(√(VC(ℬ_0^*)/N)+√(128log (1 / δ)/N)+√(ε)) ,

	where R^2=sup _v ∈ℬ_0^*𝔼(⟨ Y_1-μ, v⟩^2).


We will show that, in fact, we obtain the concentration inequality about the MOM estimator 
    ℙ(μ-μ≥ 8R√(K/N)) ≤exp(-K/128),
whenever K≥ C(VC(ℬ_0^*)∨|𝒪|), where C is a universal constant. Hence, by taking K≥ 128log (1/δ), we get the upper bound of the error.


	Let ℱ={(x_i)_i ≤ m→1_⟨1/m∑_i x_i -μ, v⟩≥ r_K, v ∈ℬ_0^*}, where
    r_K=4 sup _v ∈ℬ_0^*𝔼(⟨ Y_1-μ, v⟩^2)^1 / 2√(K/N).
 
	For any k ∈ [K], let 𝐗_k:=(X_i)_i ∈ B_k and 𝐘_k:=(Y_i)_i ∈ B_k. The functions f ∈ℱ are compositions of the function x →1/m∑_i x_i-μ and of the functions x →1_⟨ x, v⟩≥ r_K for v ∈ℬ_0^*. The VC-dimension of the set of these compositions is smaller than the VC-dimension of the set of indicator functions indexed by ℬ_0^*. We just get VC(ℱ) ≤ c_0 VC(ℬ_0^*) for some constant c_0.
	
	Notice that, the definition in (<ref>) is equivalent to
    S_x^*={y ∈ℝ^d:|x^*(X̅_k)-x^*(y)| ≤ϵ for more than K/2 blocks } .
For any f∈ℱ, or, for the corresponding v ∈ℬ_0^*, there exist at least (K-∑_k=1^K f(𝐗_k)) blocks B_k, where
	
    |⟨X̅_k-μ,v⟩| ≤ r_K.

	So it is sufficient to compute the sum of f(𝐗_k). Now we write
	
	
    sup _f ∈ℱ∑_k=1^K f(𝐘_k)=   [ sup _f ∈ℱ∑_k=1^K f(𝐘_k)-𝔼(sup _f ∈ℱ∑_k=1^K f(𝐘_k))]
       +𝔼(sup _f ∈ℱ∑_k=1^K f(𝐘_k)) .

	Let f_j(𝐘)=sup _f ∈ℱ( ∑_k j f(𝐘_k)+f(𝐘)), since f is binary-valued, we have | f_j(𝐘)-f_j(𝐘')|≤ 1 , for any j∈ [K] and 𝐘, 𝐘'∈ℝ^d× m. By the bounded difference inequality (see <cit.>),
    ℙ(sup _f ∈ℱ∑_k=1^K f(𝐘_k)-𝔼(sup _f ∈ℱ∑_k=1^K f(𝐘_k)) ≥ t)  ≤exp(-2 t^2/K).
 Therefore, by taking t=K/16, we can derive that, with probability at least 1-exp (-K/128) , the first term in (<ref>) is bounded above by K/16.
	
	For the second term, we further write
	
    𝔼(sup _f ∈ℱ∑_k=1^K f(𝐘_k)) ≤𝔼(sup _f ∈ℱ∑_k=1^K f(𝐘_k)-K 𝔼(f(𝐘_k)))+sup _f ∈ℱ K 𝔼(f(𝐘_k)) .

	By Markov's inequality, for any v ∈ℬ_0^*,
	
    ℙ(|⟨1/m∑_i ∈ B_1 Y_i-μ, v⟩| ≥ r_K) ≤𝔼(∑_i ∈ B_1⟨ Y_i-μ, v⟩^2)/m^2 r_K^2≤1/16.
Then, sup _f ∈ℱ K 𝔼(f(𝐘_k)) ≤ K/16. And by Lemma <ref>, there exists a universal constant C^' such that
	
    𝔼(sup _f ∈ℱ1/K∑_k=1^K f(𝐘_k)- 𝔼(f(𝐘_k))) ≤ C^'√(VC(ℱ)/K) .

	Hence, if K ≥ 256 C^' 2VC(ℱ),
	
    𝔼(sup _f ∈ℱ∑_k=1^K f(𝐘_k)-K 𝔼(f(𝐘_k))) ≤K/16 .

	Thus,  the second term in (<ref>) is no more than K/16+K/16=K/8.  
	
	Back to the possibly contamination model, if C ≥ 16, then K ≥ 16|𝒪| and one can deduce that 
    sup _f ∈ℱ∑_k=1^K f(𝐗_k) ≤sup _f ∈ℱ∑_k=1^K f(𝐘_k)+ K /16.

	
	Putting everything together, we derive that, if C ≥ 256 C^' 2∨ 16, the following event ℰ has probability ℙ(ℰ) ≥ 1-exp (-K / 128): for all f ∈ℱ,
	
    ∑_k=1^K f(𝐗_k) ≤K /16+K /8+K /16=K /4.

	Hence, we have proved that, for any x^* ∈ℬ_0^*,
	
    Med|x^*(X̅_k-μ)  | ≤ r_K,

	whenever K≥ C(VC(ℬ_0^*)∨|𝒪|). 	
	We conclude that taking ϵ=r_K, 𝕊(ϵ) is nonempty as it contains μ, at least on ℰ. By definition, μ_K∈𝕊(ϵ) for this choice of ϵ. Observe that for every x^* ∈ℬ_0^* there is some index j such that
	
    |x^*(X̅_j)-x^*(μ_K)| ≤ r_K  and |x^*(X̅_j)-x^*(μ)| ≤ r_K,

	because both conditions hold for more than half of the indices j. Thus,
	
    |x^*(μ_K)-x^*(μ)| ≤|x^*(X̅_j)-x^*(μ_K)|+|x^*(X̅_j)-x^*(μ)| ≤ 2 r_K .

	
	Finally, recalling that v=sup _x^* ∈ℬ_0^* x^*(v), one has that
	
    μ_K-μ=sup _x^* ∈ℬ_0^*|x^*(μ)-x^*(μ)| ≤ 2 r_K.

	Thus, for μ_K∈𝕊(r_K), and K≥ C(VC(ℬ_0^*)∨|𝒪|∨ 128log (1 / δ)), we obtain that with probability at least 1-δ,
    μ_K-μ≤ 8sup _v ∈ℬ_0^*𝔼(⟨ Y_1-μ, v⟩^2)^1 / 2√(K/N).



	Note that the construction of our estimator μ_K depends on the number of splitting blocks K, which is relative to the level δ.   Besides, the drawback is that this estimator is more theoretical than practical, for it is not convenient to construct a set like 𝕊(ϵ). As is shown in <cit.>, for every ϵ>0, the sets 𝕊(ϵ) are compact and nested set and nonempty for a sufficiently large ϵ. Therefore, the set
	
    𝕊=⋂_ϵ>0: 𝕊(ϵ)≠∅𝕊(ϵ)

	is not empty. We can define the mean estimator as any element μ∈𝕊.

Since the result is true for any norm equipped with inner product, if we consider the Euclidean norm, it shows that with probability no more than 1-δ, there exists an estimator such that

    μ-μ_δ_2≲Σ^1 / 2(√(VC(ℬ_0^*)/N)+√(128log (1 / δ)/N)+√(ε)) .

 We refer to the first two terms on the right-hand side of the above inequality as the weak and strong terms, respectively: the strong term  is a global component with the order of √(λ_1(Σ)d)/√(N), compared with the error of sub-Gaussian empirical mean in (<ref>), where the corresponding term is √(Tr(Σ))/√(N), it matches only when Σ≃λI_d; and the weak term with directional information, corresponds to the largest variance of a one dimensional marginal of Y_1, that is,  sup _v ∈B_0^*σ(v):=sup _v ∈ℬ_0^*√(𝔼(⟨ Y_1-μ, v⟩^2)); while the third terms, reflects the corruption level. Strong-weak inequalities are an important notion in high-dimensional probability, and for more improvements, one can see <cit.>, which constructs an estimator that, up to the optimal strong term, preforms robust in every direction.








§ COVARIANCE ESTIMATION


Let 0<δ<1 and consider the given sample X_1, …, X_N. Again, the sample (X_i)_i=1^N can be partitioned into K blocks B_k, each of size m=N / K.
Set M_k=1/m∑_i ∈ B_kX_i ⊗X_i. Recall the well-known fact that the dual norm to the operator norm is the nuclear norm. And, since a linear functional z acts on the matrix x via trace duality—that is, z(x)=[z, x]:=Tr(z^T x). It follows that
T={u⊗ u| u ∈ℬ_0^*(ℝ^d)} is the set of extreme points of the corresponding dual unit ball B^∘. For ϵ>0 and a fixed u∈ℬ_0^*, let U=u⊗ u, and 

    S_u(ϵ)={Y ∈ℝ^d × d:|[ M_k-Y, U] | ≤ε for more than  K / 2  blocks } .

Set

    S(ϵ)=⋂_U ∈ T S_u(ϵ) .

The estimator Σ_δ is taken as any points in S(ϵ). Again, we derive the bound of robust covariance estimator.

	For any δ∈[e^-c N, 1 / 2], there exists an estimator Σ_δ∈ S(ϵ) such that	
    Σ-Σ_δ≲σ(√(VC(ℬ_0^*)/N)+√(128log (1 / δ)/N)+√(ε)) ,

	where σ^2=sup _u ∈ℬ_2𝔼(⟨ u,(Σ-Y_1 Y_1^T) u⟩^2)<∞.


	In comparison with the similar estimator in <cit.>, our estimator does not require a two-step estimation for the trace and truncation level. Additionally, it imposes fewer assumptions. For example, L-sub-Gaussian or L_4-L_2 norm equivalence for the sample distribution are not necessary any more.

The proof of Theorem <ref> is similar to that of Theorem <ref>, and we just give the first part where we adjust the dual space and some coefficients.

	Let 
    ℱ={(x_i)_i ≤ m→1_[1/m∑_i x_i⊗ x_i -Σ, U] ≥ r_K, u ∈ℬ_0^*},
where
	r_K=4 σ√(K/N). The function f ∈ℱ are compositions of the function x →1/m∑_i x_i⊗ x_i -Σ and of the functions x →1_[ x, U] ≥ r_K for U ∈ T. The VC-dimension of the set of these compositions is smaller than the VC-dimension of the set of indicator functions indexed by T. We just get VC(ℱ)≤VC(T) ≤VC(ℬ_0^*⊗ℬ_0^*) , which is, by Lemma <ref>, bounded by c_0 VC(ℬ_0^*) for some constant c_0.
	
	By Markov's inequality, for any u ∈ℬ_0^*, that is, U∈ T,
	
    ℙ(|[1/m∑_i Y_i⊗ Y_i -Σ, U]| ≥ r_K) ≤𝔼(∑_i ∈ B_1[Y_i⊗ Y_i -Σ, U]^2)/m^2 r_K^2≤1/16 .
	
	The rest of this proof is the same as that of the preceding theorem.
	There exists an absolute constant C such that, if K ≥ C(d ∨|𝒪|), then, with probability larger than 1-exp (-K / 128),
	
    Σ_δ-Σ≤ 8 σ√(K/N) .


Similarly, we can immediately obtain expressions regarding the Frobenius norm as follows:

	For any δ∈[e^-c N, 1 / 2], there exists an estimator Σ_δ such that	
    Σ-Σ_δ_F≲σ(√(VC(ℬ_0^*)/N)+√(128log (1 / δ)/N)+√(ε)) ,

	where σ is the same as that in Theorem <ref>.

Another immediate corollary of Theorem <ref> is the quantitative result for the performance of PCA based on the estimator Σ_δ. Let Proj_k be the orthogonal projector on a subspace corresponding to the k largest positive eigenvalues λ_1, …, λ_k of Σ (here, we assume for simplicity that all the eigenvalues are distinct), and Proj_k be the orthogonal projector of the same rank as Proj_k corresponding to the k largest eigenvalues of Σ_δ. The following bound follows from the Davis-Kahan perturbation theorem in <cit.>.


	Let Δ_k=λ_k-λ_k+1, and assume that Δ_k ≥ 16 σ√(K/N). Then
	
    Proj_k-Proj_k≤8/Δ_kσ√(K/N),

	with probability at least (1-exp (-K / 128)).





§ TUKEY'S MEDIAN

We consider a special type of median of mean, where we take the Tukey's median (see <cit.>), as a robust estimator. First, we need to introduce halfspace depth function. For any η∈ℝ^d and a distribution ℙ on ℝ^d, the halfspace depth of η with respect to ℙ is defined as

    𝒟(η, ℙ)=inf _u ∈ S^d-1ℙ{u^T X ≤ u^T η}   where  X ∼ℙ .


Given i.i.d. observations {X_i}_i=1^N, the halfspace depth of η with respect to the observations {X_i}_i=1^N is defined as

    𝒟(η,{X_i}_i=1^N)=𝒟(η, ℙ_N)=min _u ∈ S^p-11/N∑_i=1^N 1_{u^T X_i ≤ u^T η},

where ℙ_N=1/N∑_i=1^N δ_X_i is the empirical distribution. Then Tukey's median is defined to be the deepest point with respect to the observations, that is,

    θ̂=max _η∈ℝ^d𝒟(η,{X_i}_i=1^N) .



 §.§ Halfspace depth


In the matter of the coherence with the preceding expressions, we consider the use of its equivalent definition
    θ̂=min _η∈ℝ^dsup_v ∈ℬ_0^*∑_i=1^N1_⟨ X_i-η,v⟩ > 0,

when (<ref>) has multiple maxima, θ̂ is understood as any vector that attains the deepest level. As is known from <cit.>, the maximum depth N𝒟(θ̂,{X_i}_i=1^N) is bounded below by ⌈ N /(d+1)⌉. Cause the natural form of boolean functions, it suggests us to use the VC-dimension technique. From now on, we take our estimator, the Tukey median of mean, as
    μ̂=min _η∈ℝ^dsup_v ∈ℬ_0^*∑_k=1^K1_⟨X̅_k-η,v⟩ > 0.
  Then we obtain the following dimension dependent bound:

	There exists an universival constant C such that if K ≥ C(VC(ℬ_0^*) ∨|𝒪|), then, with probability larger than 1-exp (-K/32d^2),
	
    μ̂_K-μ≤√(8d)sup _v ∈ℬ_0^*𝔼(⟨ Y_1-μ, v⟩^2)^1 / 2√(K/N) .


In particular, for any δ∈[e^-c N, 1 / 2], and if we use the Euclidean distance, there exists an estimator μ̂ such that

    μ̂-μ≲√(d)Σ^1 / 2(√(VC(ℬ_0^*)/N)+d√(log (1 / δ)/N)+√(ϵ)).

Before providing the proof of the theorem, for ease of manipulation, we first provide a lemma in the general form without proof, which can in fact be derived from the proof process of Theorem <ref>.


	For any k ∈ [K], let 𝐗_k:=(X_i)_i ∈ B_k and 𝐘_k:=(Y_i)_i ∈ B_k. Let ℱ be a Boolean class of functions satisfying the following two assumptions:
	
		
  * For all f ∈ℱ, ℙ(f(Y_1)=1) ≤1 4α,
		
  * K ≥ C(VC(ℱ) ∨|𝒪|) where C is a universal constant,
	
	where α>1 can be any constant. Then, with probability at least 1-exp (-K/8α^2), for all f ∈ℱ, there is at least α-1α K blocks B_k on which f(X_k)=0.
	
 

	Let K ≥ C(VC(ℱ) ∨|𝒪|) with C the universal constant from Lemma  <ref>, and let  ℱ={(x_i)_i ≤ m→1_⟨1/m∑_i x_i -μ, v⟩≥ r_K,d, v ∈ℬ_0^*}, where
    r_K,d=√(8d)sup _v ∈ℬ_0^*𝔼(⟨ Y_1-μ, v⟩^2)^1 / 2√(K/N).
 Based on the analysis in the previous Section <ref>, VC(ℱ) ≤ c_0 VC(ℬ_0^*) for some constant c_0.
	By Markov's inequality, for any v ∈ℬ_0^*,
	
    ℙ(|⟨1/m∑_i ∈ B_1 Y_i-μ, v⟩| ≥ r_K,d) ≤𝔼(∑_i ∈ B_1⟨ Y_i-μ, v⟩^2)/m^2 r_K,d^2≤1/8d .

	
	By Lemma <ref>, applied with ℱ and α=2d, we derive that, with probability ≥ 1-exp (-K/32d^2), there is at least 2d-1 2d K blocks B_k on which the following event happens:
	
    sup _v ∈ℬ_0^*|⟨1/m∑_i ∈ B_k X_i-μ, v⟩| ≤ r_K,d.

	We claim that, for any v ∈ℬ_0^*,
	
    ⟨μ-μ̂, v⟩≤  r_K,d.

	That is because, for any η∈ℝ^d attains the deepest level, if there exists v^* ∈ℬ_0^* such that ⟨μ-η, v^*⟩> r_K,d, then, on the above event,
    ⟨X̅_k-η, v^*⟩   =⟨X̅_k-μ, v^*⟩+⟨μ-η, v^*⟩
       >⟨X̅_k-μ, v^*⟩+r_K,d
       ≥ 0

	holds for at least 2d-1 2d K blocks B_k, which means that, 
    sup_v ∈ℬ_0^*∑_k=1^K1_⟨X̅_k-η,v⟩ > 0≥∑_k=1^K1_⟨X̅_k-η,v^*⟩ > 0≥ d K / (d+1).

	Therefore η≠μ̂. Hence, we have proved the claim, ⟨μ-μ̂, v⟩≤  r_K,d. Take the supremum over v ∈ℬ_0^* and since sup _v ∈ℬ_0^*⟨μ-μ̂, v⟩=μ-μ̂, we conclude the proof.
	

The results obtained from Theorem <ref> in inequality (<ref>) are not entirely satisfactory due to the influence of the dimension d. In comparison to the traditional empirical mean, it provides a better approximation rate only for the contamination model. However, this upper bound is not ideal for high-dimensional heavy-tailed data, especially when d≥ N.

While <cit.> point out that the depth of our Tukey's median lies between ⌈ N /(d+1)⌉ and ⌈ N / 2⌉, and if the dataset is nearly symmetric (about some point x_0), the maximum depth will be much larger than N /(d+1), in fact approximately N / 2. According to the above proof, it can be inferred that, under the assumption of a symmetric enough distribution and the level 1-δ, Tukey MOM estimator can achieve the same error, as our initially proposed estimator, of Σ^1 / 2(√(VC(ℬ_0^*)/N)+√(log (1 / δ)/N)+√(ε)) .



 §.§ Discussion on regression

The problem of regression function estimation involves estimating conditional expectations, making it a natural extension of the mean estimation concepts discussed in this paper. This section explores recent advancements in regression problems driven by uniform median-of-means estimators.

The standard framework for regression function estimation is as follows. Consider a pair of random variables (Y, V), where Y takes values in the set 𝒳 and V is real-valued. In a class ℱ comprising real-valued functions defined on 𝒳, the goal is to identify f ∈ℱ such that f(Y) serves as a reliable prediction of V. The efficacy of a predictor f ∈ℱ is measured through the mean-squared error 𝔼(f(Y)-V)^2, known as the risk. The optimal performance within the class is achieved by the risk minimizer

    f^*=f ∈ℱargmin𝔼(f(Y)-V)^2.


The joint distribution of (Y, V) is usually unknown. Instead, an i.i.d. sample 𝒟_N=(Y_i, V_i)_i=1^N is provided, distributed according to the joint distribution of Y and V. Given a sample size N, a learning procedure is a mapping assigning to each sample 𝒟_N a function in ℱ, denoted as f.

The effectiveness of f is assessed based on the trade-off between accuracy ϵ and confidence δ in which f achieves that accuracy. In other words, one seeks f that satisfies the condition

    ℙ(𝔼((f(Y)-V)^2 |𝒟_N) ≤inf _f ∈ℱ𝔼(f(Y)-V)^2+ϵ) ≥ 1-δ

for values of ϵ and δ as small as possible. The exploration of this accuracy/confidence trade-off has been the focus of extensive research(see <cit.>, <cit.>, and <cit.>). 

Consider the standard linear regression setting, where, ℱ={⟨β, Y⟩: β∈ℝ^d}, and

    β^*=β∈ℝ^dargmin l(β)=β∈ℝ^dargmin𝔼(V_1-⟨β, Y_1⟩)^2 .


Since for all β,
l(β)-l(β^*)=2 𝔼(ξ_1⟨β-β^*, Y_1⟩)+(β-β^*) Σ(β-β^*) ≤(β-β^*) Σ(β-β^*),
the key to control excess risk is to bound the β-β^*_Σ= (β-β^*) Σ(β-β^*) . One can see the natural decomposition of the quadratic component and the multiplier one (see <cit.>) leads to applying the small-ball method to learning problems. It can thus be suggested that we could introduce Tukey MOM along with VC-dimension.



The advantage of using halfspace depth is that it offers rich geometric properties related to depth, which are valuable for exploration. Additionally, there are existing algorithms that simplify the computation process. We expect that there will be favorable properties for depth functions related to linear regression problems. These properties are crucial for obtaining robust estimates of regression coefficients β, thereby constraining the excess risk in regression (learning) problems. 



§ CONCLUDING REMARKS
 
As noted earlier, the method we present may be implemented in mean estimation, covariance estimation, and other learning problems. There has been a substantial amount of prior work on MOM estimation related to depth. In this context, we present an estimation method of Tukey MOM. Besides, <cit.> offers a theoretical framework for random approximations of Tukey depth, establishing consistency and convergence rates for their consistency with theoretical depths under general probability measures, the computation of Tukey depth constitutes a scholarly pursuit, as elucidated in  <cit.>. 

The introduction of VC dimension has provided new ideas for estimating the upper bounds of empirical processes. It must be reiterated that we have circumvented Rademacher complexity and, correspondingly, replaced it with the VC dimension term. The critical issue in the proof of bound errors is to find the proper boolean functions class and the use of Lemma <ref>.  Consequently, in the estimation of the error bound, the focus is more on the dimensional structure of the set/space, thereby mitigating the impact of heavy tails and contamination in the samples. At the same time, we have solely addressed the statistical convergence rates in the estimation problem, without considering efficiency concerns in computation. 

<cit.> shows that, using the generic chaining method, Catoni estimator can also be linked to the upper bounds of random processes, which can be used to analyze and control the upper bounds of risks in the empirical risk minimization process. Moreover, in the multivariate setting, different definitions of the median lead to different MOM estimators.  Apart from Tukey median,
many other types of median have been developed, such as the coordinate-wise median, the geometric (or spatial) median, the Oja median, and the Liu median, among others; see <cit.> for a survey.  Therefore, one can also consider similar results of the Catoni estimator and other types of MOM estimators with VC dimension.



§ ACKNOWLEDGEMENT
Hanchao Wang  was supported by  the National Natural Science Foundation of China (No. 12071257 and No. 11971267); National Key R&D Program of China (No. 2018YFA0703900 and No. 2022YFA1006104); Shandong Provincial Natural Science Foundation (No. ZR2019ZD41).

Lixin Zhang was supported by grants from the NSF of China (Grant Nos. U23A2064 and 12031005).












