On Weighted Multivariate Sign Functions

TL;DR

This paper introduces data-dependent weights to multivariate sign functions, enhancing robustness and efficiency in estimation and inference, with applications in location estimation, PCA, dimension reduction, and outlier detection.

Contribution

It proposes a novel weighted sign function approach that improves robustness and efficiency over traditional methods in multivariate analysis.

Findings

Weighted signs retain robustness properties.

Significant efficiency improvements demonstrated.

Effective in various multivariate applications.

Abstract

Multivariate sign functions are often used for robust estimation and inference. We propose using data dependent weights in association with such functions. The proposed weighted sign functions retain desirable robustness properties, while significantly improving efficiency in estimation and inference compared to unweighted multivariate sign-based methods. Using weighted signs, we demonstrate methods of robust location estimation and robust principal component analysis. We extend the scope of using robust multivariate methods to include robust sufficient dimension reduction and functional outlier detection. Several numerical studies and real data applications demonstrate the efficacy of the proposed methodology.

Tables

Table 1. Table 1: Finite sample efficiencies of estimators of the first eigenvector based on several scatter matrices in dimension p = 4 𝑝 4 p=4 . The notation H, M or P after Σ ~ ~ Σ \tilde{\Sigma} or Σ ∗ subscript Σ {\Sigma}_{*} indicates the depth function used for the weights: H = halfspace depth, M = Mahalanobis depth, P = projection depth.

4-variate $t_5$	SCM	Tyler	$\tilde{\Sigma }$-H	$\tilde{\Sigma }$-M	$\tilde{\Sigma }$-P	${\Sigma }_{\ast }$-H	${\Sigma }_{\ast }$-M	${\Sigma }_{\ast }$-P
$n$=20	1.04	1.02	1.10	1.07	1.02	1.09	1.07	0.98
$n$=50	1.08	1.08	1.16	1.16	1.13	1.19	1.19	1.13
$n$=100	1.31	1.31	1.42	1.38	1.36	1.46	1.44	1.36
$n$=300	1.46	1.54	1.81	1.76	1.95	1.88	1.88	1.95
$n$=500	1.92	1.93	2.23	2.03	2.31	2.35	2.19	2.39
4-variate $t_6$	SCM	Tyler	$\tilde{\Sigma }$-H	$\tilde{\Sigma }$-M	$\tilde{\Sigma }$-P	${\Sigma }_{\ast }$-H	${\Sigma }_{\ast }$-M	${\Sigma }_{\ast }$-P
$n$=20	1.00	1.05	1.03	1.05	1.00	1.04	1.04	0.95
$n$=50	1.03	1.01	1.13	1.12	1.11	1.19	1.17	1.10
$n$=100	1.08	1.12	1.25	1.23	1.27	1.24	1.25	1.22
$n$=300	1.34	1.36	1.64	1.52	1.60	1.67	1.61	1.68
$n$=500	1.26	1.34	1.55	1.49	1.60	1.65	1.61	1.69
4-variate $t_{10}$	SCM	Tyler	$\tilde{\Sigma }$-H	$\tilde{\Sigma }$-M	$\tilde{\Sigma }$-P	${\Sigma }_{\ast }$-H	${\Sigma }_{\ast }$-M	${\Sigma }_{\ast }$-P
$n$=20	0.90	0.89	0.95	0.98	0.98	0.96	1.01	0.95
$n$=50	0.90	0.91	1.01	0.98	0.98	1.03	1.04	0.99
$n$=100	0.87	0.87	0.93	0.95	1.01	0.99	1.01	1.05
$n$=300	0.87	0.87	1.09	1.09	1.17	1.14	1.16	1.23
$n$=500	0.88	0.92	1.10	1.10	1.23	1.19	1.22	1.29
4-variate $t_{15}$	SCM	Tyler	$\tilde{\Sigma }$-H	$\tilde{\Sigma }$-M	$\tilde{\Sigma }$-P	${\Sigma }_{\ast }$-H	${\Sigma }_{\ast }$-M	${\Sigma }_{\ast }$-P
$n$=20	0.92	0.90	0.94	0.94	0.96	0.95	0.97	0.89
$n$=50	0.82	0.83	0.88	0.91	0.93	0.88	0.93	0.93
$n$=100	0.84	0.87	0.92	0.95	1.00	0.93	0.96	1.00
$n$=300	0.73	0.75	0.96	0.99	1.10	1.00	1.06	1.12
$n$=500	0.73	0.76	0.95	0.96	1.06	0.94	0.97	1.06
4-variate $t_{25}$	SCM	Tyler	$\tilde{\Sigma }$-H	$\tilde{\Sigma }$-M	$\tilde{\Sigma }$-P	${\Sigma }_{\ast }$-H	${\Sigma }_{\ast }$-M	${\Sigma }_{\ast }$-P
$n$=20	0.89	0.92	0.92	0.92	0.90	0.96	0.95	0.89
$n$=50	0.82	0.84	0.89	0.90	0.91	0.93	0.96	0.92
$n$=100	0.77	0.76	0.90	0.90	0.96	0.94	0.98	1.04
$n$=300	0.73	0.77	0.93	0.91	0.98	1.00	0.98	1.03
$n$=500	0.67	0.71	0.83	0.83	0.96	0.88	0.90	1.00
4-variate Normal	SCM	Tyler	$\tilde{\Sigma }$-H	$\tilde{\Sigma }$-M	$\tilde{\Sigma }$-P	${\Sigma }_{\ast }$-H	${\Sigma }_{\ast }$-M	${\Sigma }_{\ast }$-P
$n$=20	0.82	0.84	0.87	0.90	0.91	0.89	0.93	0.89
$n$=50	0.80	0.81	0.87	0.88	0.88	0.88	0.92	0.88
$n$=100	0.68	0.71	0.80	0.85	0.91	0.82	0.86	0.92
$n$=300	0.61	0.63	0.82	0.85	0.93	0.86	0.91	0.96
$n$=500	0.60	0.64	0.77	0.80	0.90	0.82	0.86	0.96

Table 2. Table 2: Table of AREs of the estimator for the first eigenvector estimation using Σ ∗ subscript Σ \Sigma_{*} , relative to using the sample covariance matrix, for different choices of dimension p 𝑝 p . The data-generating distributions are the multivariate Normal (MVN), and multivariate t 𝑡 t -distributions with degrees of freedom 5, 6, 10, 15 and 25. Weights for Σ ∗ subscript Σ \Sigma_{*} are based on either the projection depth (PD) or the half-space depth (HSD).

	PD				HSD
Distribution	$p=2$	$p=5$	$p=10$	$p=20$	$p=2$	$p=5$	$p=10$	$p=20$
$t_5$	4.73	3.99	3.46	3.26	4.18	3.63	3.36	3.15
$t_6$	2.97	3.28	2.49	2.36	2.59	2.45	2.37	2.32
$t_{10}$	1.45	1.47	1.49	1.52	1.30	1.37	1.43	1.49
$t_{15}$	1.15	1.19	1.23	1.27	1.01	1.10	1.17	1.24
$t_{25}$	0.97	1.02	1.07	1.11	0.85	0.94	1.02	1.08
MVN	0.77	0.84	0.89	0.93	0.68	0.77	0.84	0.91

Equations

$$\displaystyle S(x;\mu)=\left\{\begin{array}[]{ll}|x-\mu|^{-1}(x-\mu)&\text{ if }x\neq\mu,\\ 0&\text{ if }x=\mu.\end{array}\right.$$

$$R(X_i;\mu ,\mathbb{F})=S(X_i,\mu )W(X_i,\mathbb{F}).$$

$$\displaystyle\Psi(q;u,{\mathbb{F}})={\mathbb{E}}\Bigl{[}W(X,{\mathbb{F}})\bigl{\{}|X-q|+\langle u,X-q\rangle\bigr{\}}\Bigr{]}={\mathbb{E}}\Phi(q;X,u,{\mathbb{F}}).$$

$${\Psi }_{1W}$$

$${\Psi }_{2W}$$

$$n^{1/2}({\hat{q}}_{nW}-q_0)$$

$$n^{1/2}({\hat{q}}_{nW}-q_0)=-\frac{{\Psi }_{2W}}{\sqrt{n}}\sum _{i=1}^n\nabla \Phi (q;X_i,u,\mathbb{F})+o_P(1),$$

$$\displaystyle ARE(\hat{q}_{nW},\hat{q}_{n})=\Bigl{\{}\frac{det(V_{1})}{det(V_{W})}\Bigr{\}}^{1/p}.$$

$$ARE({\hat{q}}_{nW},{\hat{q}}_n)\ge \frac{{\lambda }_{\min }({\Psi }_1){\lambda }_{\min }^2({\Psi }_{2W})}{W_{\max }{\lambda }_{\max }({\Psi }_{1W}){\lambda }_{\max }^2({\Psi }_2)}.$$

$$\frac{det(V_1)}{det(V_W)}$$

$$=det({\Psi }_1)det({\Psi }_{1W}^{-1})[det({\Psi }_2^{-1})det({\Psi }_{2W}){]}^2.$$

$$\frac{W_{\max }}{[\mathbb{E}|X-q_0{|}^{-1}W(X,\mathbb{F}){]}^2}<\frac{{\lambda }_{\min }({\Psi }_1)}{{\lambda }_{\max }({\Psi }_1)[\mathbb{E}|X-q_0{|}^{-1}{]}^2]}.$$

$$D({\mu }_{\mathbb{F}};\mathbb{F})\ge D({\mu }_{\mathbb{F}}+t(\mathbf{x}-{\mu }_{\mathbb{F}});\mathbb{F}).$$

$$X=\mu +R\Gamma {\Lambda }^{1/2}U,$$

$$S(X;\mu )$$

$$=|{\Lambda }^{1/2}U{|}^{-1}\Gamma {\Lambda }^{1/2}U=|{\Lambda }^{1/2}Z{|}^{-1}\Gamma {\Lambda }^{1/2}Z.$$

$$\tilde{X}=W(X,\mathbb{F})S(X;\mu )\equiv W(Z,{\mathbb{F}}_Z)|{\Lambda }^{1/2}Z{|}^{-1}\Gamma {\Lambda }^{1/2}Z.$$

$$\tilde{\Sigma }=\mathbb{E}\tilde{X}{\tilde{X}}^T=\mathbb{E}{W}^2(X,\mathbb{F})\mathbb{S}(X;\mu ).$$

$$\displaystyle\tilde{U}_{j}=\left\{\begin{array}[]{ll}U_{j}&\text{ if }j\neq i,\\ -U_{i}&\text{ if }j=i.\end{array}\right.$$

$$\mathbb{E}\frac{U_i{U}_j}{U^T\Lambda U}$$

$$\tilde{\Sigma }$$

$$=\mathbb{E}{W}^2(X)\mathbb{E}S(X;\mu )S(X;\mu {)}^T=\Gamma {\Lambda }_W{\Gamma }^T,$$

$$\underset{x\in \mathcal{X}}{\sup }|W(x,{\mathbb{F}}_n)-W(x,\mathbb{F})|\to 0$$

$$\displaystyle\sup_{x\in\mathcal{X}}\Bigl{|}W\bigl{(}x,{\mathbb{F}}\bigr{)}-W\bigl{(}x,(1-\epsilon){\mathbb{F}}+\epsilon\delta_{x}\bigr{)}\Bigr{|}\leq\epsilon.$$

$$\frac{1}{n}\sum _{i=1}^n{W}_n^2(X_i,{\mathbb{F}}_n)\mathbb{S}(X_i;{\hat{\mu }}_n)=\frac{1}{n}\sum _{i=1}^n{W}^2(X_i,\mathbb{F})\mathbb{S}(X_i;\mu )+R_n,$$

$$\frac{1}{n}\sum _{i=1}^n{W}_n^2(X_i,{\mathbb{F}}_n)\mathbb{S}(X_i;{\hat{\mu }}_n)$$

$$=$$

$$\displaystyle+\frac{1}{n}\sum_{i=1}^{n}W^{2}(X_{i},{\mathbb{F}})\bigl{\{}{\mathbb{S}}(X_{i};\hat{\mu}_{n})-{\mathbb{S}}(X_{i};{\mu})\bigr{\}}$$

$$\displaystyle+\frac{1}{n}\sum_{i=1}^{n}\bigl{\{}W_{n}^{2}(X_{i},{\mathbb{F}}_{n})-W^{2}(X_{i},{\mathbb{F}})\bigr{\}}\bigl{\{}{\mathbb{S}}(X_{i};\hat{\mu}_{n})-{\mathbb{S}}(X_{i};{\mu})\bigr{\}}$$

$$=\frac{1}{n}\sum _{i=1}^n{W}^2(X_i,\mathbb{F})\mathbb{S}(X_i;\mu )+T_2+T_3+T_4.$$

$$\displaystyle n^{-1}\sum_{i=1}^{n}|X_{i}-\mu|^{-2}\bigl{\{}W_{n}^{2}(X_{i},{\mathbb{F}}_{n})-W^{2}(X_{i},{\mathbb{F}})\bigr{\}}(X_{i,j}-\mu_{j})(X_{i,k}-\mu_{k}),$$

$$c^T{T}_{2}c=\sum _{j,k}{c}_j{c}_k{T}_{2,j,k}$$

$$\displaystyle=n^{-1}\sum_{i=1}^{n}|X_{i}-\mu|^{-2}\bigl{\{}W_{n}^{2}(X_{i},{\mathbb{F}}_{n})-W^{2}(X_{i},{\mathbb{F}})\bigr{\}}\bigl{(}\sum_{j}c_{j}(X_{i,j}-\mu_{j})\bigr{)}^{2}$$

$$\displaystyle\leq Mn^{-1}\sum_{i=1}^{n}|X_{i}-\mu|^{-2}\bigl{\{}|W_{n}(X_{i},{\mathbb{F}}_{n})-W(X_{i},{\mathbb{F}})|\bigr{\}}\bigl{(}c^{T}(X_{i}-\mu)\bigr{)}^{2}$$

$$\displaystyle\leq Mn^{-1}\sum_{i=1}^{n}|X_{i}-\mu|^{-2}\bigl{\{}|W_{n}(X_{i},{\mathbb{F}}_{n})-W_{n}(X_{i},{\mathbb{F}}_{n,-i})|\bigr{\}}\bigl{(}c^{T}(X_{i}-\mu)\bigr{)}^{2}$$

$$\displaystyle\hskip 28.45274pt+Mn^{-1}\sum_{i=1}^{n}|X_{i}-\mu|^{-2}\bigl{\{}|W_{n}(X_{i},{\mathbb{F}}_{n,-i})-W(X_{i},{\mathbb{F}})|\bigr{\}}\bigl{(}c^{T}(X_{i}-\mu)\bigr{)}^{2}$$

Full text

On Weighted Multivariate Sign Functions

Subhabrata Majumdar

and

Snigdhansu Chatterjee

University of Minnesota, Minneapolis, MN, USA Currently at AT&T Labs Research. Email: subho@research.att.comCorresponding author. Email: chatt019@umn.edu

Abstract: Multivariate sign functions are often used for robust estimation and inference. We propose using data dependent weights in association with such functions. The proposed weighted sign functions retain desirable robustness properties, while significantly improving efficiency in estimation and inference compared to unweighted multivariate sign-based methods. Using weighted signs, we demonstrate methods of robust location estimation and robust principal component analysis. We extend the scope of using robust multivariate methods to include robust sufficient dimension reduction and functional outlier detection. Several numerical studies and real data applications demonstrate the efficacy of the proposed methodology.

Keywords: Multivariate sign, Principal component analysis, Data depth, Sufficient dimension reduction

1 Introduction

Given a point $\mu$ in a normed linear space $\mathcal{X}$ with norm denoted by $|\cdot|$, the generalized sign function $S:\mathcal{X}\times\mathcal{X}\mapsto\mathcal{X}$ with center $\mu$ is defined as

[TABLE]

This is a functional and multivariate generalization of the real-valued sign function, that takes the values one, negative one or zero if the point $x\in{\mathbb{R}}$ is to the right, left or equal $\mu\in{\mathbb{R}}$ respectively. This generalized sign function was introduced by Möttönen and Oja, (1995) for $\mathcal{X}={\mathbb{R}}^{p}$, the $p$-dimensional real Euclidean space.

The function $S$ maps $\mu$ to the origin and all other points of $\mathcal{X}$ to the unit sphere ${\mathcal{S}}_{0;1}=\{x\in\mathcal{X}:|x|=1\}$. Given a dataset $\{X_{i}\in{\mathbb{R}}^{p}:i=1,\ldots,n\}$, that we collect together in the $n\times p$ matrix ${\mathbf{X}}=(X_{1}:\ldots:X_{n})^{T}$, an approach for robust estimation and inference in multivariate data starts by evaluating (3) on each observation, thus defining $S_{i}=S(X_{i};\mu_{i})$ with respect to some suitable center $\mu_{i}\in{\mathbb{R}}^{p}$, and then using these for robust location and scale estimation and inference, including inference for $\mu_{i}$ (Locantore et al.,, 1999; Oja,, 2010; Wang et al.,, 2015). Suppose ${\mathbf{S}}=(S_{1}:\ldots:S_{n})^{T}\in{\mathbb{R}}^{n\times p}$. If the data $\{X_{i}\in{\mathbb{R}}^{p}:i=1,\ldots,n\}$ are independent, identically distributed (hereafter, i.i.d.) from an elliptically symmetric distribution, then the eigenvectors of ${\mathbb{E}}(X_{1}-\tilde{\mu})(X_{1}-\tilde{\mu})^{T}$ and of ${\mathbb{E}}S_{1}S_{1}^{T}$ are the same for suitable centering parameter $\tilde{\mu}\equiv\mu_{i}$, that is, the population principal components from the original data and from its sign transformations are the same (Taskinen et al.,, 2012). However, valuable information is lost in the form of magnitudes of sample points. As a result, spatial sign-based procedures suffer from low efficiency. For example, eigenvector estimates obtained from the covariance matrix of $S$ are asymptotically inadmissible (Magyar and Tyler,, 2014) and Tyler’s M-estimate of scatter (Tyler,, 1987) has uniformly lower asymptotic risk.

In this paper, we propose to alleviate this low efficiency problem, by associating a data-driven weight $W_{i}$ with the generalized sign $S_{i}$, that can be used to adaptively trade-off between efficiency and robustness considerations in any given application. We demonstrate the utility of using the proposed weighted generalized sign functions in a number of problems of current interest, including robust estimation of location and scatter.

Specifically, we propose using product of the generalized sign function and a weight function derived as a transformation of a data-depth function (Serfling,, 2006; Zuo and Serfling,, 2000). Like data-depth functions, the weight functions used in this paper are non-negative reals defined over $\mathcal{X}\times\mathcal{F}$, where $\mathcal{F}$ is a fixed family of probability measures. For every choice of parameters $\mu\in\mathcal{X}$ and ${\mathbb{F}}\in\mathcal{F}$, in this paper

[TABLE]

is used as a robust surrogate for observation $X_{i}$. Notice that for the trivial choice $W(x,{\mathbb{F}})=|x-\mu_{{\mathbb{F}}}|$, $\mu=\mu_{{\mathbb{F}}}$, we get $R(X_{i};\mu,{\mathbb{F}})=X_{i}-\mu_{{\mathbb{F}}}$, the original centered observations. With the other trivial choice of $W(x,{\mathbb{F}})\equiv 1$, we get the generalized sign $R(X_{i};\mu,{\mathbb{F}})=S(X_{i},\mu)=S_{i}$. However, in this paper we illustrate how using other weight functions can lead to interesting robustness and efficiency trade-offs in a variety of situations. The various technical conditions and assumptions that we impose on the weight function $W(x,{\mathbb{F}})$ are valid for weights derived from three well-known data depth functions: the half-space depth, the Mahalanobis depth, and the projection depth. Note that for i.i.d. data, there are three parameters involved here: the location parameter $\mu$ in (3), the distribution ${\mathbb{F}}$ used for the weight function, and the distribution ${\mathbb{F}}_{X}$ of $X_{1}$. For clarity and to mirror the contexts of how data depth has been used in the literature (Liu et al.,, 1999; Serfling,, 2006), we fix ${\mathbb{F}}={\mathbb{F}}_{X}$ for this paper, although in Section 2 we briefly remark on the case when these distributions are different. Also note that in many problems of interest, ${\mathbb{F}}$ is unknown and data-depths are computed using ${\mathbb{F}}_{n}$, however, under very standard regularity conditions (for example, see assumptions B1-B3 below), the properties of $W(x,{\mathbb{F}})$ and $W(x,{\mathbb{F}}_{n})$ are close enough for both asymptotic theory and practical applicability.

We primarily focus on the task of robust dispersion/scatter estimation and robust principal component analysis in this paper. Figure 1 presents an illustrative example of bivariate data with outliers in the top left panel, where the outliers are marked with red points. In the other panels, the generalized sign values of the same data are presented as black points on the unit circle, with the outliers again marked with red points. Notice that the black points from either the top right or bottom panels have very similar eigenvector structure as the original data without the outliers. The green and blue triangles are examples of the proposed weighted sign values: the top right (respectively, bottom left, bottom right) panels depict these values where the weights have been generated using Mahalanobis depth (respectively Tukey’s half-space depth and the projection depth). The blue triangles are the weighted sign values of the outliers. Notice that the eigenvectors from the weighted signs also capture the pattern from the original data without the outliers.

There are two unknown quantities in the generalized sign function as defined in (3): $\mu$ and ${\mathbb{F}}$, To estimate dispersion and its eigen-structure robustly, we must start with a robust estimator for $\mu$. In Section 2 we present the case for weighted spatial quantiles, which can be defined and studied in very general spaces $\mathcal{X}$. One special case of this is the weighted spatial median. As a location estimator, it has several interesting robustness properties and can be shown to be more efficient that some existing robust location estimators, thus making it a perfect candidate to estimate $\mu$.

Following that, we restrict to the ${\mathbb{R}}^{p}$ for fixed $p$, and present detailed discussions on our primary proposal for a robust measure of dispersion in Section 3, followed by a proposed * affine equivariant version* of it in Section 4, robust estimation of eigenvalues and a third robust estimator for dispersion in Section 5, and a thorough study of robustness and efficiency using influence functions in Section 6. We then report multiple simulation-based numeric studies in Section 7, present several real data examples in Section 8, and concluding remarks in Section 9.

In the rest of this paper, all finite-dimensional vectors are column vectors, and for a vector or matrix $a$, the notation $a^{T}$ stands for its transpose. The Gaussian distribution with mean $\mu$ and variance $\Sigma$ is denoted by $N(\mu,\Sigma)$. The identity matrix is denoted by ${\mathbb{I}}$, with or without a subscript to denote its dimension. The notations $A^{-1},det(A),\lambda_{\min}(A),\lambda_{\max}(A)$ respectively stand for the inverse, determinant, minimum and maximum eigenvalues of a matrix $A$, whenever these are well-defined. For a scalar or vector valued random variable $Y$, ${\mathbb{E}}Y$ denotes its expected value, while ${\mathbb{V}}Y$ denotes its variance or covariance matrix.

2 The Weighted Spatial Median

Suppose the open unit sphere in $\mathcal{X}$ is given by $\textrm{int}{\mathcal{X}}_{0;1}=\{x\in\mathcal{X}:|x|<1\}$, and let $u\in\textrm{int}{\mathcal{X}}_{0;1}$. We also fix the set of probability measures $\mathcal{M}$, and select ${\mathbb{F}}\in\mathcal{M}$. Consider a random element $X\in\mathcal{X}$, and define the function \Phi(q;X,u,{\mathbb{F}})=W(X,{\mathbb{F}})\bigl{\{}|X-q|+\langle u,X-q\rangle\bigr{\}}. We define the $(u,{\mathbb{F}})$-th weighted spatial quantile of $\mathcal{X}$ as the minimizer $q(u,{\mathbb{F}})\in\mathcal{X}$ of the expectation of $\Phi(q;X,u,{\mathbb{F}})$, that is

[TABLE]

This is a natural generalization of the spatial median (Chaudhuri,, 1996; Haldane,, 1948; Koltchinskii,, 1997) ($W(X,{\mathbb{F}})\equiv 1$ and $u={\bf 0}_{p}$), or more general spatial quantiles (Cardot et al.,, 2017; Chakraborty and Chaudhuri,, 2014; Minsker,, 2015) ($W(X,{\mathbb{F}})\equiv 1$). We assume that $\Phi(q;X,u,{\mathbb{F}})$ is convex in $q$ for ${\mathbb{F}}$-almost all values of $x\in\mathcal{X}$.

In what follows, for brevity we elaborate only the case of the weighted spatial median (thus \Psi(q;0,{\mathbb{F}})={\mathbb{E}}\bigl{[}W(X,{\mathbb{F}})|X-q|\bigr{]}) for the case of finite dimensional $\mathcal{X}$. Specifically, we demonstrate the utility of using the weighted spatial median as opposed to using the traditional, unweighted versions found in literature. The analysis and computation of the general weighted spatial quantiles will largely follow by extending the results of the above cited references, and we postpone details to a future study.

Assume that we have independent and identically distributed observations $X_{1},\ldots,X_{n}\in\mathcal{X}$, and the sample weighted spatial median is computed by minimizing $\Psi_{n}(q;0,{\mathbb{F}})=\sum_{i=1}^{n}W(X_{i},{\mathbb{F}})|X_{i}-q|$, and is denoted by $\hat{q}_{nW}$, the second subscript is in acknowledgement that the weight function is used. We denote the unweighted version of this estimator, ie, the case where $W(X,{\mathbb{F}})\equiv 1$ as $\hat{q}_{n}$. Assume the following technical conditions:

(A1)

$\Psi(q;0,{\mathbb{F}})$ is finite for all $q\in\mathcal{X}\subseteq{\mathbb{R}}^{p}$ and has a unique minimizer $q_{0}$.

(A2)

$\Psi(q;0,{\mathbb{F}})$ is twice differentiable at $q_{0}$ and the second derivative is positive definite.

(A3)

${\frac{\partial^{2}}{\partial q^{2}}}\Psi(q;0,{\mathbb{F}})$ exists for all $q$ in a neighborhood of $q_{0}$, and we use the notations

[TABLE]

These assumptions are very broad-based and general. The first one essentially requires the existence of a population parameter, the second one requires that the minimization approach is meaningful in the population, and the third one essentially requires that the weight function has a finite variance. No further restrictions are placed on the tuning parameter ${\mathbb{F}}$ or the choice of the weight function.

Theorem 2.1**.**

Under assumptions [A1]-[A3], we have

[TABLE]

Thus, under very standard regularity conditions, the sample weighted spatial median is consistent and is asymptotically normal. Theorem 2.1 can be proved in several different ways. Here we use techniques following Haberman, (1989); Niemiro, (1992). Specifically, following Theorem 4 in Niemiro, (1992), which traces back to Haberman, (1989) with slightly relaxed conditions, we get

[TABLE]

where $\nabla\Phi(q;X_{i},u,{\mathbb{F}})$ is any measurable subgradient of $\Phi(q;X_{i},u,{\mathbb{F}})$. Theorem 2.1 follows by applying the central limit theorem.

Remark.

Note that the technical conditions for the result presented in Theorem 2.1 is one of several alternatives that can conceived, and the scope of this result is broader than what is presented above. First, note that if ${\mathbb{F}}$ and ${\mathbb{F}}_{X}$ are different and the weights are not a function of $q$, a situation that may arise in hypothesis testing problems where the weights are based on the null distribution, the convexity of $\Phi(q;X_{i},u,{\mathbb{F}})$ follows automatically and is not an assumption. Second, even if $\Phi(q;X_{i},u,{\mathbb{F}})$ is not convex but sufficiently smooth, we can have a central limit theorem, for example, by using techniques similar to Chatterjee and Bose, (2005). Choices of ${\mathbb{F}}$ other than ${\mathbb{F}}_{X}$, e.g. Majumdar and Chatterjee, (2018), may lead to interesting interpretations of $W(.,{\mathbb{F}})$ and the resulting location estimator and will be explored further in future.

2.1 Asymptotic efficiency of weighted spatial median

Let $V_{W}=\Psi_{2W}^{-1}\Psi_{1W}\Psi_{2W}^{-1}$ be the asymptotic variance of $\hat{q}_{n,W}$ from Theorem 2.1, where we use the subscript “W” to denote that this depends on the weight function. We use the notation $V_{1}$ for the case where $W(x,{\mathbb{F}})\equiv 1$, that is, all weights are one. The asymptotic relative efficiency of two statistics is the $p$-th root of the reciprocals of their determinants. That is,

[TABLE]

One easy result from Theorem 2.1 is that under reasonable conditions the asymptotic relative efficiency of the weighted spatial median over the unweighted version is always greater than one. We document this in the following corollary:

Corollary 2.2**.**

Assume that the weight function $W(X,{\mathbb{F}})$ is bounded above by some $W_{\max}>0$, and the matrices $\Psi_{1}={\mathbb{E}}S(X;q_{0})S^{T}(X;q_{0})$ and $\Psi_{1W}$ are positive definite. Then

[TABLE]

Consequently, if $W_{\max}/\lambda_{\min}^{2}(\Psi_{2W})<\lambda_{\min}(\Psi_{1})/(\lambda_{\max}(\Psi_{1W})\lambda_{\max}^{2}(\Psi_{2}))$ then this asymptotic relative efficiency is larger than 1.

Proof of Corollary 2.2.

Using the facts that $det(AB)=det(A)det(B)$ for square matrices $A,B$ and $det(A^{-1})=1/det(A)$ for non-singular $A$, we write

[TABLE]

The result follows, using the facts that $\det(A)\geq\lambda_{\min}(A)$ and $\det(A^{-1})\geq 1/\lambda_{\max}(A)$, and the upper bound on $W$. ∎

We may wish to further explore the conditions of Corollary 2.2. Let us concentrate on the case of where the distribution of $X$ is spherically symmetric. Following Fang et al., (1990):

Definition 2.3**.**

A $p$-dimensional random vector $X$ is said to elliptically distributed if there exist a vector $\mu\in{\mathbb{R}}^{p}$, a positive semi-definite matrix $\Sigma\in{\mathbb{R}}^{p\times p}$ and a function $\phi:(0,\infty)\rightarrow{\mathbb{R}}$ such that the characteristic function of $X$ is $\exp\{it^{T}\mu\}\phi({\bf{t}}^{T}\Sigma{\bf{t}})$ for ${\bf{t}}\in{\mathbb{R}}^{p}$.

There are several alternative formulations, see the above reference and citations within it for details.

Corollary 2.4**.**

Assume that ${\mathbb{F}}\equiv{\mathbb{F}}_{X}$ is an elliptically symmetric distribution, with location parameter $\mu=q_{0}$ and the covariance matrix $\Sigma$ satisfies the following conditions:

$\exp\left\{-\frac{Tr^{2}(\Sigma)}{256\lambda_{\max}^{2}(\Sigma)}\right\}=o\left(\min\left(\frac{\lambda_{\max}(\Sigma)}{Tr(\Sigma)},\frac{\lambda_{\min}(\Sigma)}{\lambda_{\max}(\Sigma)}\right)\right)$, 2. 2.

$\lambda_{\max}(\Sigma)=o(Tr(\Sigma))$.

Also assume that the weight function is (a) affine invariant, i.e. $W(Ax+b,A{\mathbb{F}}+b)=W(x,{\mathbb{F}})$ for $A\in{\mathbb{R}}^{p\times p},b\in{\mathbb{R}}^{p}$, and (b) satisfies the following:

[TABLE]

Then we have $ARE(\hat{q}_{n,W},\hat{q}_{n})>1$.

The conditions 1 and 2 in Corollary 2.4 are due to Wang et al., (2015), and ensure that as $p\equiv p_{n}$ diverges with $n$, the eigenvalues of $\Sigma\in{\mathbb{R}}^{p\times p}$ are bounded away from 0 and $\infty$. Corollary 2.4 can be proved using Corollary 2.2, the fact that for elliptical distributions $\Psi_{1}$ is non-singular (Taskinen et al.,, 2012), then using and expanding upon Lemma A.5 in Wang et al., (2015). We give the technical details in the supplementary material.

Obtaining affine invariant weight functions is not challenging. Most functions arising in the context of data depths are affine invariant. Corollary 2.4 implies that there is a wide frame of distributions and choices of weight functions where there is a benefit to considering weighted spatial medians. In fact, Corollary 2.4 can be established under just location invariance condition on weights. We use the stronger affine invariance condition here only because the subsequent sections use weights based on affine invariant data depth functions. Also note that in case (6) does not hold for the original form of a weight function, one can scale all weights by an appropriate constant to reduce the value in the left hand side of (6) to satisfy the condition.

In actual computations, in place of $W(x,{\mathbb{F}})$ we propose using $W(x,{\mathbb{F}}_{n})$ where ${\mathbb{F}}_{n}(z)=n^{-1}\sum_{i=1}^{n}\prod_{j=1}^{p}\mathcal{I}(X_{ij}\leq z_{j}),z\in{\mathbb{R}}^{p}$ is the empirical distribution function. Up to first order asymptotics, the analysis remain unchanged from the above for this modified weight function.

2.2 Examples of affine invariant weights

We now illustrate some specific choices of weight functions that are compatible with the conditions of Corollary 2.4 and the results in the rest of this paper. These arise as easy transformations of data-depth functions. A data depth function is defined on $\mathcal{X}\times\mathcal{F}$, where $\mathcal{F}$ is a fixed set of probability measures. The main property of a data-depth function is that for every probability measure ${\mathbb{F}}\in\mathcal{F}$, there exists a constant $\mu_{{\mathbb{F}}}\in\mathcal{X}$ such that for any $t\in[0,1]$ and $x\in\mathcal{X}$,

[TABLE]

That is, for every fixed ${\mathbb{F}}$, the data-depth function achieves a supremum at $\mu_{F}$, and is non-decreasing in every direction away from $\mu_{F}$, thus providing a center-outward partial ordering of points in $\mathcal{X}$. There are generally several algebraic and analytic properties assumed for data-depth functions to elicit interesting mathematical properties, see for example Serfling, (2006); Zuo and Serfling, (2000) for details.

The spherically symmetric case of an elliptical distribution is realized with $\Sigma=\sigma^{2}{\mathbb{I}}_{p}$ for some $\sigma^{2}>0$. We fix the notation $Z=\Sigma^{-1/2}(X-\mu)$, and let $Z\sim{\mathbb{F}}_{Z}$. Note that ${\mathbb{F}}_{Z}$ is a spherically symmetric distribution and hence depends only on $|z|$, and ${\mathbb{E}}Z=\mathbf{0}_{p}\in{\mathbb{R}}^{p}$ and ${\mathbb{V}}Z={\mathbb{I}}_{p}$. Taking affine invariant data depth functions as weights ensures that $W(X,{\mathbb{F}})=W(Z,{\mathbb{F}}_{Z})$. It is now easy to show that results in this paper are valid for the weight functions (with appropriate scaling to satisfy (6)) ($i$) $W_{HSD}(X)={\mathbb{F}}_{Z_{1}}(|Z|)$ derived from the half-space depth, ($ii$) $W_{MhD}(X)=|Z|^{2}/(1+|Z|^{2})$ derived from the Mahalanobis depth, and ($iii$) $W_{PD}(X)=|Z|/(1+|Z|/MAD(Z_{1}))$, where $MAD$ stands for median absolute deviation, derived from the projection depth. We omit the technical details. These three weight functions give a center-inward partial ordering, thus essentially quantifying peripherality instead of depth. Note however, that our results below are of much more general form, and these three special choices of weights only serve as important illustrative examples to achieve desirable robustness and efficiency balance in data analysis.

3 A robust measure of dispersion

From this section onwards, we assume that $\mathcal{X}={\mathbb{R}}^{p}$, that is, the support of the random variable under study is the $p$-dimensional Euclidean plane, and the data $X_{1},\ldots,X_{n}$ are independent and identically distributed from an elliptical distribution ${\mathbb{F}}$ with parameters $\mu$ and $\Sigma$. We also assume that $X_{1}$ is absolutely continuous, with ${\mathbb{P}}[|X_{1}|=0]=0$, and that $\Sigma$ is positive definite. This eliminates technicalities arising from rank deficient cases, and makes the weight functions affine invariant. Thus, we essentially restrict the rest of this paper to the same framework as in Corollary 2.4. Most of the results below generalize to the case where $\mathcal{X}$ is a separable Hilbert space, however additional technicalities are involved, as in Bali et al., (2011), and will be considered in a future project.

3.1 The Weighted Sign Covariance Matrix

Consider the spectral decomposition of $\Sigma$ given by $\Sigma=\Gamma\Lambda\Gamma^{T}$, where $\Gamma$ is an orthogonal matrix and $\Lambda$ is diagonal with positive diagonal elements $\lambda_{1}\geq\ldots\geq\lambda_{p}$. Also denote the $i$-th eigenvector of $\Sigma$ by ${\boldsymbol{\gamma}}_{i}=(\gamma_{i,1},\ldots,\gamma_{i,p})^{T}$ for $1\leq i\leq p$. Thus, the $i$-th column of $\Gamma$ is ${\boldsymbol{\gamma}}_{i}$. In the rest of this paper we use the notation $\Sigma^{-1/2}=\Lambda^{-1/2}\Gamma^{T}$, and hence $Z=\Lambda^{-1/2}\Gamma^{T}(X-\mu)$. Recall from Section 2 that we use the notation ${\mathbb{F}}_{Z}$ for the distribution of $Z$, and that ${\mathbb{F}}_{Z}$ is a spherically symmetric distribution and hence depends only on $|z|$. Additionally, to simplify notations, for any random variable $X\sim{\mathbb{F}}$, we occasionally use the abbreviated notation $W(X)\equiv W(X,{\mathbb{F}})$. Note that $W(X)$ is a random weight, and takes the same value as $W(Z,{\mathbb{F}}_{Z})$. Also, as noted in Corollary 2.4, $W(Z)$ is a function of $|Z|$ only. We additionally assume that ${\mathbb{E}}W^{2}(X)<\infty$.

It is convenient to write

[TABLE]

where $U$ is a random variable uniformly distributed on the unit sphere ${\mathcal{S}}_{0;1}=\{x\in\mathcal{X}:|x|=1\}$ and $R$ is another random variable independent of $U$ satisfying ${\mathbb{E}}R^{2}=p$. Note that $Z=RU$, and $|Z|=R$, $Z/|Z|=U$. Then we have

[TABLE]

As a robust surrogate for $X-\mu$, we consider the following random variable

[TABLE]

In samples, the equivalent for $\tilde{X}$ is $\hat{\tilde{X}}=W(X,{\mathbb{F}}_{n})S(X;\hat{\mu})$ for a suitable location estimator $\hat{\mu}$, for example, the weighted spatial median. We fix the notation ${\mathbb{S}}(X;\mu)=S(X;\mu)S(X;\mu)^{T}$, and define the following dispersion parameter:

[TABLE]

In the following Theorem, we establish that the eigenvectors of ${\Sigma}$ and $\tilde{\Sigma}$ are identical, although their eigenvalues may be different.

Theorem 3.1**.**

Under the conditions listed above, we have $\tilde{\Sigma}=\Gamma\tilde{\Lambda}\Gamma^{T}$, where $\tilde{\Lambda}=\Lambda^{1/2}{\mathbb{E}}W^{2}(X){\mathbb{E}}[UU^{T}/(U^{T}\Lambda U)]\Lambda^{1/2}$ is a diagonal matrix. Thus, the eigenvectors of ${\Sigma}$ and $\tilde{\Sigma}$ are identical.

Proof of Theorem 3.1.

Fix any index $i\in\{1,\ldots,p\}$. Consider the vector $\tilde{U}$ such that

[TABLE]

Then $\tilde{U}$ and $U$ have the same distribution, and note that $U^{T}\Lambda U=\tilde{U}^{T}\Lambda\tilde{U}$ almost surely. Consequently, for any $j\neq i$ we have

[TABLE]

Consequently, ${\mathbb{E}}S(X;\mu)S(X;\mu)^{T}=\Gamma\Lambda_{S}\Gamma^{T}$, as established in Theorem 1 of Taskinen et al., (2012).

Also, since the weight $W(X)$ is a function of $|Z|=R$, we have that $W(X)$ is independent of $S(X;\mu)$. Consequently, we have

[TABLE]

where $\Lambda_{W}$ is a diagonal matrix. ∎

3.2 Sample version of $\tilde{\Sigma}$

We now discuss the properties of the sample version of $\tilde{\Sigma}$, say $\widehat{\tilde{\Sigma}}$ computed from ${\mathbf{X}}$. In practice, we cannot obtain $W(x)\equiv W(x,{\mathbb{F}})$, and consequently use $W(x,{\mathbb{F}}_{n})$ instead. We assume the following conditions:

(B1)

Bounded weights: The weights $W(\cdot,\cdot)$ are bounded functions.

(B2)

Uniform convergence:

[TABLE]

almost surely as $n\rightarrow\infty$.

(B3)

Smoothness under perturbation: For all ${\mathbb{F}}\in\mathcal{F}$, there exists a $\delta>0$, possibly depending on ${\mathbb{F}}$, such that for any $\epsilon\in(0,\delta)$

[TABLE]

In the above, $\delta_{x}$ denotes point mass at $x$. These properties are easily satisfied for weight functions derived from standard depth functions, for example, $W_{HSD}(\cdot)$, $W_{MhD}(\cdot)$ and $W_{PD}(\cdot)$ discussed earlier.

The following result allows us to use the empirical, plug-in weights and an estimated location parameter in the weighted dispersion estimator. A natural choice for the location parameter estimator is the solution to $\sum_{i=1}^{n}\tilde{X}_{i}=0$, which is the same as the sample version of the weighted spatial median discussed in Section 2.

Lemma 3.2**.**

Assume that ${\mathbb{E}}|X-\mu|^{-4}<\infty$. Also assume that we have a location estimator $\hat{\mu}_{n}$ satisfying ${\mathbb{E}}|\hat{\mu}_{n}-\mu|^{4}=O(n^{-2})$. Then

[TABLE]

where for any $c\in{\mathbb{R}}^{p}$ such that for $|c|=1$, we have ${\mathbb{E}}c^{T}R_{n}c=o(n^{-1})$.

Proof of Lemma 3.2.

Note that the moment condition above is satisfied by default since $X$ has an elliptical distribution (Dürre et al.,, 2014). The proof is mostly algebra, and we provide a sketch of the main arguments. We have

[TABLE]

Using the stated technical conditions, we can now show that ${\mathbb{E}}c^{T}T_{j}c=o(n^{-1})$ for $j=2,3,4$. For illutration, we present the case for $T_{2}$ below.

Notice that the $(j,k)$-th element of $T_{2}$ is given by

[TABLE]

and hence

[TABLE]

Let H(X_{i})=|X_{i}-\mu|^{-2}\bigl{(}c^{T}(X_{i}-\mu)\bigr{)}^{2}, and notice that $H(X_{i})\leq 1$ almost surely for $|c|=1$. Now notice that conditional on $X_{i}$ except for a null set $A_{i}$ (possibly depending on $X_{i}$) we have $T_{21i}\leq n^{-1}H(X_{i})$. Thus, except for a null set $A_{1}\cap\ldots\cap A_{n}$, $T_{21}\leq Mn^{-2}H(X_{i})$ and the conclusion follows for this part.

The argument for $T_{22}$ follows a similar argument. ∎

Let $vec({\mathbb{S}}(X;\mu))$ be the vectorized version of ${\mathbb{S}}(X;\mu)$. We are now in a position to state the result for consistency of the sample version of $\tilde{\Sigma}$,

Theorem 3.3**.**

Assume the conditions of Lemma 3.2. Then

[TABLE]

where $V_{W}={\mathbb{V}}[W^{2}(X,{\mathbb{F}})vec({\mathbb{S}}(X;{\mu}))]$.

The asymptotic normality follows from our assumptions and as a direct consequence of Lemma 3.2. Incidentally, an expression for $V_{W}$ can be explicitly obtained in terms of $\Gamma$, $\Lambda$ and ${\mathbb{F}}$, but is algebraic in nature. We present it in the supplementary material.

We now use Theorem 3.3 to obtain consistency results for the eigenvectors obtained from $\widehat{\tilde{\Sigma}}=n^{-1}\sum_{i=1}^{n}W^{2}(X_{i},{\mathbb{F}}){\mathbb{S}}(X_{i};\hat{\mu}_{n})$. Suppose that $\tilde{\Lambda}_{1}>\tilde{\Lambda}_{2}>\ldots>\tilde{\Lambda}_{p}$ are the eigenvalues of $\tilde{\Sigma}$, which we assume are all distinct values.

Theorem 3.4**.**

Suppose the spectral decomposition of $\widehat{\tilde{\Sigma}}$ is given by $\widehat{\tilde{\Sigma}}=\widehat{\Gamma}\widehat{\tilde{\Lambda}}\widehat{\Gamma}^{T}$. Then the matrix of centered and scaled eigenvectors $G_{n}=n^{1/2}(\widehat{\Gamma}-\Gamma)$ and the vector of centered and scaled eigenvalues $L_{n}=n^{1/2}(\widehat{\tilde{\Lambda}}-{\tilde{\Lambda}})$ have independent distributions. The distribution of the random variable $vec~{}(G_{n})$ converges weakly to a $p^{2}$-variate normal distribution with mean ${\bf 0}_{p^{2}}$ and the variance matrix whose $(i,j)$-th block of $p\times p$ entries are given by

[TABLE]

The distribution of $L_{n}$ converges weakly to a $p$-dimensional normal distribution with mean ${\bf 0}_{p}$ and the variance-covariance matrix whose $(i,j)$-the element is

[TABLE]

The proof of this result follows from using Theorem 3.3 and using techniques similar to a corresponding result in Taskinen et al., (2012). We omit the algebraic details here and put it in the supplementary material.

Recall that the asymptotic variance of the $i$-th eigenvector of the sample covariance matrix, say $\hat{{\boldsymbol{\gamma}}}_{i}$ is (Anderson,, 2009):

[TABLE]

Suppose $\widehat{\tilde{{\boldsymbol{\gamma}}}}_{i}$ is the $i$-th eigenvector of $\widehat{\tilde{\Sigma}}$, whose asymptotic behavior is presented above in Theorem 3.4.

This leads to the following useful result:

Corollary 3.5**.**

The asymptotic relative efficiency of $\widehat{\tilde{{\boldsymbol{\gamma}}}}_{i}$, relative to $\hat{{\boldsymbol{\gamma}}}_{i}$, is given by

[TABLE]

The proof of this Corollary is immediate.

4 An affine equivariant robust measure of dispersion

A desirable invariance property of any dispersion parameter $T_{X}$ corresponding to a random variable $X$ is that under affine transformation $Y=AX+b$ the dispersion parameter scales to $T_{Y}=AT_{X}A^{T}$. It is clear that $\tilde{\Sigma}$ does not possess this property, since it remains unchanged for $X$ and $Y=cX$ for any scalar $c>0$.

We follow the general framework of M-estimation with data-dependent weights (Huber,, 1981) to construct an affine equivariant counterpart of the $\tilde{\Sigma}$. Specifically, we implicitly define

[TABLE]

To ensure existence and uniqueness of $\Sigma_{*}$, consider the class of dispersion parameters $\Sigma_{M}$ that are obtained as solutions of the following equation:

[TABLE]

with $Z_{M}=\Sigma_{M}^{-1/2}(X-\mu)$. Under the following assumptions on the scalar valued functions $u$ and $v$, the above equation produces a unique solution (Huber,, 1981):

(C1)

The function $u(r)/r^{2}$ is monotone decreasing, and $u(r)>0$ for $r>0$;

(C2)

The function $v(r)$ is monotone decreasing, and $v(r)>0$ for $r>0$;

(C3)

Both $u(r)$ and $v(r)$ are bounded and continuous,

(C4)

$u(0)/v(0)<p$,

(C5)

For any hyperplane in the sample space $\mathcal{X}$, (i) $P(H)={\mathbb{E}}\{\mathcal{I}_{\{X\in H\}}\}<1-pv(\infty)/u(\infty)$ and (ii) $P(H)\leq 1/p$.

Putting things into context, in our case we have $v(\cdot)=p^{-1}{\mathbb{V}}W(X)$, $u(\cdot)=W^{2}(X)$. We proceed to verify the other conditions for the weight functions $W_{HSD}(\cdot)$, $W_{MhD}(\cdot)$ and $W_{PD}(\cdot)$ discussed earlier.

It is easy to verify that the resulting $u(\cdot)$ from the above choices satisfy (C1) and (C3). Note that $v(\cdot)$ is a finite positive constant, and (C2) and (C3) are also easily satisfied. Since $u(0)=0$ in all the above cases, (C4) is also easy to check. Since $X$ is absolutely continuous, (C5) holds trivially.

To compute the sample version of $\Sigma_{*}$, we solve (13) iteratively by obtaining a sequence of positive definite matrices $\hat{\Sigma}^{(k)}_{*}$ until convergence. Thus, using the location estimator $\hat{\mu}_{n}$, we may iterate

[TABLE]

The asymptotic properties of $\hat{\Sigma}_{*}$ can be obtained using methods similar to those of Section 3, and techniques presented in Croux and Haesbroeck, (2000) and elsewhere. We state the following result and omit its proof.

Theorem 4.1**.**

The asymptotic covariance matrix of an eigenvector of the sample affine equivariant scatter functional $\hat{\Sigma}_{*}$ is given by

[TABLE]

where $V_{12}$ is the asymptotic variance of an off-diagonal element of $\hat{\Sigma}_{*}$ when the underlying distribution is ${\mathbb{F}}_{Z}$. It follows that if $\hat{{\boldsymbol{\gamma}}}_{*,i}$ is the $i$-th eigenvector of $\hat{\Sigma}_{*}$,

[TABLE]

5 Robust estimation of eigenvalues, and a plug-in estimator of $\Sigma$

As seen in Theorem 3.1, eigenvalues of the $\tilde{\Sigma}$ are not same as the population eigenvalues. In this section, we discuss on robust estimation of $\lambda_{i}$’s using $\tilde{\Sigma}$. Assume the data is centered, the robust estimator from Section 2 suffices. We start by computing the sample version $\widehat{\tilde{\Sigma}}$ and its spectral decomposition:

[TABLE]

We then use the following steps:

Randomly divide the sample indices $\{1,2,\ldots,n\}$ into $k$ disjoint groups $\{G_{1},\ldots,G_{k}\}$ of size $\lfloor n/k\rfloor$ each. 2. 2.

Transform the data matrix: ${\mathbf{S}}=\widehat{{\Gamma}}^{T}{\mathbf{X}}$. 3. 3.

Calculate coordinate-wise variances for each group of indices $G_{j}$:

[TABLE]

where $\bar{{\mathbf{S}}}_{G_{j}}=(\bar{S}_{G_{j},1},\ldots,\bar{S}_{G_{j},p})^{T}$ is the vector of column-wise means of ${\mathbf{S}}_{G_{j}}$, the submatrix of ${\mathbf{S}}$ with row indices in $G_{j}$. 4. 4.

Obtain estimates of eigenvalues by taking coordinate-wise medians of these variances:

[TABLE]

We collect ${\lambda}^{\dagger}_{i}$, $i=1,\ldots,p$ in the diagonal matrix ${\Lambda}^{\dagger}=\operatorname*{diag}({\lambda}^{\dagger}_{1},\ldots,{\lambda}^{\dagger}_{p})$. The number of subgroups used to calculate this median-of-small-variances estimator can be determined following Minsker, (2015). There can be other ways of estimating the eigenvalues of $\Sigma$ using ${\mathbf{S}}$ also, we will pursue such methods elsewhere. We construct a consistent plug-in estimator of the population covariance matrix $\Sigma$ as

[TABLE]

Let $|A|_{F}$ denote the Frobenius norm of a matrix $A$, in other words, $|A|_{F}=(\text{trace}A^{T}A)^{1/2}$. The following result establishes that this is a consistent estimator of $\Sigma$:

Theorem 5.1**.**

Suppose that as $n\rightarrow\infty$, $k\rightarrow\infty$ and $n/k\rightarrow\infty$. Then we have

[TABLE]

Proof of Theorem 5.1.

This proof has many algebraic steps, and we sketch the main arguments below.

Suppose $\hat{A}=\widehat{{\Gamma}}^{T}\Sigma\widehat{{\Gamma}}$. Owing to the fact that the Frobenius norm is invariant under rotations and that $p$ is finite and fixed, it suffices to show that the off-diagonal elements of $\hat{A}$ converge in probability to zero, and that the difference between the $i$-th diagonal element of $\hat{A}$ and ${\lambda}^{\dagger}_{i}$ converges to zero for any $i=1,\ldots,p$.

Now notice that from Theorem 3.4 we have that $\widehat{{\Gamma}}={{\Gamma}}+R_{n1}$, where the $(i,j)$-th element of the remainder $R_{n1,i,j}$ satisfies ${\mathbb{E}}R_{n1,i,j}^{2}=O(n^{-1})$. We can show, using standard algebra, that

[TABLE]

where the $(i,j)$-th element of the remainder $R_{n2,i,j}$ satisfies ${\mathbb{E}}R_{n2,i,j}^{2}=O(n^{-1})$. This follows immediately from above, the fact that $p$ is finite and fixed, and all elements of $\Lambda$ are constants. This immediately establishes the case for the off-diagonal elements.

For the diagonal elements, notice that since $k\rightarrow\infty$, each coordinate-wise variance $\lambda_{i,j}^{\dagger}$ for each group of indices $G_{j}$ is a consistent estimator of $\lambda_{i}$. The result follows.

∎

6 Influence Functions of Dispersion Measures

We retain the framework adopted in Section 3, and discuss in this section the robustness and efficiency properties associated with $\tilde{\Sigma}$ and $\Sigma_{*}$, and principal components derived therefrom. We do not discuss $\Sigma^{\dagger}$ here, since the properties of that approach follow from those of $\tilde{\Sigma}$. We additionally assume that the eigenvalues of $\Sigma$ are distinct, and given by $\lambda_{1}>\lambda_{2}>\ldots>\lambda_{p}$, to avoid several additional technical conditions for the theoretical results to follow. The case where the eigenvalues of $\Sigma$ can have multiplicity greater than one requires no additional conceptual development, but does require considerable algebraic manipulations.

For studying the robustness aspect, we first present some results relating to influence functions in the current context. Influence functions quantify how much influence a sample point, especially an infinitesimal contamination, has on any functional of a probability distribution (Hampel et al.,, 1986). Given any probability distribution ${\mathbb{H}}\in\mathcal{M}$, the influence function of any point $x_{0}\in\mathcal{X}$ for some functional $T({\mathbb{H}})$ on the distribution is defined as

[TABLE]

where ${\mathbb{H}}_{\epsilon}=(1-\epsilon){\mathbb{H}}+\epsilon\delta_{x_{0}}$; $\delta_{x_{0}}$ being the distribution with point mass at $x_{0}$. When $T({\mathbb{H}})=E_{\mathbb{H}}f$ for some ${\mathbb{H}}$-integrable function $f$, $IF(x_{0};T,{\mathbb{H}})=f(x_{0})-T({\mathbb{H}})$.

It now follows that

[TABLE]

Recall that $\tilde{\lambda}_{1}>\tilde{\lambda}_{2}>\ldots>\tilde{\lambda}_{p}$ are the eigenvalues of $\tilde{\Sigma}$, which we assume are all distinct values.

Proposition 6.1**.**

The influence function of ${\boldsymbol{\gamma}}_{i}$ as follows:

[TABLE]

The proof of Proposition 6.1 follows from Sibson, (1979) and Croux and Haesbroeck, (2000), we omit the details.

If the weight function $W(\cdot)$ is a bounded function, as is the case of $W_{HSD}$, $W_{MhD}$, and $W_{PD}$, the influence function given in Proposition 6.1 is bounded, indicating good robustness properties of the principal component analysis.

We now derive the influence function for $\Sigma_{*}$.

Proposition 6.2**.**

The influence function of $\Sigma_{*}$ is given by

[TABLE]

Proof of Proposition 6.2.

Let $z_{0}=\Lambda^{-1/2}\Gamma^{T}(x_{0}-\mu)=(z_{01},\ldots,z_{0p})^{T}$. As a first step, since $\Sigma_{*}$ is affine equivariant, we obtain from Croux and Haesbroeck, (2000) that

[TABLE]

From Lemma 1 of (Hampel et al.,, 1986), page 276, we obtain that there exist scalar valued functions $\alpha_{\Sigma_{*}}(|x_{0}|;{\mathbb{F}}_{Z})$ and $\beta_{\Sigma_{*}}(|x_{0}|;{\mathbb{F}}_{Z})$ such that

[TABLE]

consequently we obtain

[TABLE]

∎

Suppose ${\lambda}_{*1}>{\lambda}_{*2}>\ldots>{\lambda}_{*p}$ are the eigenvalues of ${\Sigma}_{*}$, which we assume are all distinct values. Also denote the $i$-th eigenvector of $\Sigma_{*}$ by ${\boldsymbol{\gamma}}_{*i}=(\gamma_{*i1},\ldots,\gamma_{*ip})^{T}$ for $1\leq i\leq p$.

Proposition 6.3**.**

The influence function of ${\boldsymbol{\gamma}}_{*i}$ may be obtained as

[TABLE]

We omit the proof of Proposition 6.3, which follows along similar lines to to rest of the computations of this section. It can be shown that when $W(\cdot)$ is a bounded function, $\alpha_{\Sigma_{*}}(|x_{0}|;{\mathbb{F}}_{Z})$ is also bounded, along the lines of Huber, (1981), which in turn implies that the influence function for a principal component based on $\Sigma_{*}$ is also bounded.

7 Simulation Studies

We report a number of numerical simulation studies on several properties relating to $\tilde{\Sigma}$ and $\Sigma_{*}$, and their eigenvalues and eigenvectors, on datasets with or without influential points, to illustrate the finite sample efficiency and robustness properties of the proposed weighted estimators. We compare these proposed estimators with techniques that exists in literature, specifically, the Sign Covariance Matrix (SCM) and Tyler’s shape matrix (Tyler,, 1987).

7.1 Efficiency of different robust estimators

We compare the performance of $\tilde{\Sigma}$ and $\Sigma_{*}$ with that of the SCM and Tyler’s scatter matrix. For this study, we fix the dimension $p=4$. We consider six elliptical distributions, and from every distribution draw 10000 samples each for sample sizes $n=20,50,100,300$ and $500$. All distributions are centered at ${\bf 0}_{p}$, and have covariance matrix $\Sigma=\operatorname*{diag}(4,3,2,1)$.

We use the concept of principal angles (Miao and Ben-Israel,, 1992) to find out error estimates for the first eigenvector of a scatter matrix. In our case, the first eigenvector is

[TABLE]

We measure the prediction error for an eigenvector estimator (say, $\tilde{\boldsymbol{\gamma}}_{1}$), using the smallest angle between the true and predicted vectors, i.e. $\cos^{-1}|\tilde{\boldsymbol{\gamma}}_{1}^{T}\hat{\boldsymbol{\gamma}}_{1}|$. A small absolute value of this angle means to better prediction. We repeat this 10,000 times and calculate the Mean Squared Prediction Angle:

[TABLE]

where $\tilde{\boldsymbol{\gamma}}^{(m)}_{1}$ is the value of $\tilde{\boldsymbol{\gamma}}_{1}$ in the $m$-th replication, $m=1,\ldots,10,000$. The finite sample efficiency of $\tilde{\boldsymbol{\gamma}}_{1}$ relative to that from the sample covariance matrix, i.e. $\hat{\boldsymbol{\gamma}}_{1}$ is obtained as:

[TABLE]

The results from this simulation exercise are presented in Table 1. It can be seem that $\tilde{\Sigma}$-based estimators (columns 3-5) outperform SCM and Tyler’s $M$-estimator of scatter. Among the 3 depth functions considered, projection depth gives the best results. Its finite sample performances are better than Tyler’s and Huber’s M-estimators of scatter, as well as their symmetrized counterparts that are much more computationally intensive (see Table 4 in Sirkiä et al., (2007)). The affine equivariant ${\Sigma}_{*}$-based estimators (columns 6-8) are even more efficient.

7.2 Influence function comparison

In Figure 2 we consider first eigenvectors of $\tilde{\Sigma}$, the Sign Covariance Matrix (SCM) and Tyler’s shape matrix (Tyler,, 1987). We generate data from and set ${\mathbb{F}}\equiv\mathcal{N}_{2}({0},\operatorname*{diag}(2,1))$ and plot norms of the eigenvector influence functions for different values of $x_{0}$. Let us denote the $i$-th eigenvector of the Sign Covariance Matrix and Tyler’s shape matrix by $\gamma_{S,i}$ and $\gamma_{T,i}$, respectively. Their influence functions are given as follows:

[TABLE]

Panels (b) and (c) in Figure 2, corresponding to Sign Covariance Matrix and Tyler’s shape matrix respectively, exhibit an ‘inlier effect’, that is, points close to the center having high influence, which results in loss of efficiency. On the other hand, the influence function for eigenvector estimates of the sample covariance matrix (panel (a)) is unbounded and makes the corresponding estimates non-robust. In comparison, the $\tilde{\Sigma}$ corresponding to weights derived from projection depth, half-space depth and Mahalanobis depth have bounded influence functions and small values of the influence function at ‘deep’ points.

7.3 Efficiency of affine equivariant robust estimator

We study the finite sample efficiency properties of $\Sigma_{*}$ using a simulation exercise. We consider 6 different elliptic distributions, namely, the $p$-variate multivariate Normal distribution and the multivariate $t$ distributions corresponding to degrees of freedom 5, 6, 10, 15 and 25. We compute the ARE of the estimator for the first eigenvector using $\Sigma_{*}$, using weights based on the projection depth (PD) and the halfspace depth (HSD), thus this simulation is an illustration of how different choices of weights affect the results in the context of Theorem 4.1. We consider using the sample covariance matrix as the baseline method for this study. The ARE values are computed by using Monte-Carlo simulation of $10^{6}$ samples and subsequent numerical integration. We report the results of this exercise in Table 2. Based on these results, we notice that $\Sigma_{*}$ is particularly efficient in lower dimensions for distributions with heavier tails ($t_{5}$ and $t_{6}$), while for distributions with lighter tails, the AREs increase with data dimension. At higher values of $p$, note that $\Sigma_{*}$ is almost as efficient as the sample covarnace matrix even when the data comes from multivariate normal distribution.

7.4 Robust sufficient dimension reduction and supervised learning

One of the main usages of obtaining dispersion estimators and their eigenvalues and eigenvectors is in dimension reduction techniques. Examples of such uses are in principal component regression, partial least squares and envelope methods. We illustrate below the latter technique, in the context of sufficient dimension reduction (SDR). For details on envelope methods and other uses of robust estimators of dispersion and eigen-structures, see Cook et al., (2010); Adragni and Cook, (2009); Cook and Zhang, (2015) and references and citations of these articles. In the context of multivariate-response ($Y_{i}\in{\mathbb{R}}^{q}$) linear regression, the envelope method proposes the model $Y_{i}=\alpha+\Gamma_{1}\eta x_{i}+e_{i}$, where $e_{i}$ are independent mean zero Gaussian noise terms with covariance matrix $\Sigma$ whose spectral representation can be written as

[TABLE]

Thus, the eigenvectors of $\Sigma$ are partitioned into two blocks: $\Gamma_{1}\in{\mathbb{R}}^{q}\times{\mathbb{R}}^{d}$ and $\Gamma_{0}\in{\mathbb{R}}^{q}\times{\mathbb{R}}^{q-d}$, and the regression coefficient of $Y_{i}$ on $x_{i}$ is given by $\Gamma_{1}\eta$ for some $\eta\in{\mathbb{R}}^{d}\times{\mathbb{R}}^{p}$. Dimension reduction is achieved when $d\ll p$, typically without extraneous assumptions like sparsity. The envelope model for generalized linear models is discussed in Adragni and Cook, (2009), and may be used for supervised learning. Nonlinear regression models may also be handled similarly.

Given a set of examples $\{(Y_{i},X_{i}),i=1,\ldots,n\}$, an envelope-based prediction for the response $Y$ for any $X$ may be obtained from

[TABLE]

The above assumes that the covariates come from the Gaussian distribution $N_{p}({\bf 0}_{p},\sigma^{2}{\mathbb{I}}_{p})$, and appropriate changes may be made for other distributions.

We design a robust version of the above, by using weighted spatial medians for location parameters corresponding to the distributions of $X$ and $X|Y$, and using the first $d$ eigenvectors of $\tilde{\Sigma}$ as $\hat{\Gamma}_{1}$. A robust location estimator for the distribution of $X|Y$ is required for the estimation of $\sigma^{2}$. Details are available in Adragni and Cook, (2009).

In a non-linear regression model, we compare the performance of the robust version of SDR with the original method of Adragni and Cook, (2009) with or without the presence of bad leverage points in $\Sigma$. For any given choice of covariate dimension $p$, we take $n=200$ and $d=1$, and generate the responses $Y_{1},\ldots,Y_{n}$ as independent standard normal, and $X|Y$ as Normal with mean $Y+Y^{2}+Y^{3}$ in each of the $p$ coordinates, and variance $25{\mathbb{I}}_{p}$. We measure performances of the SDR models by their mean squared prediction error on another set of 200 observations generated similarly, and taking the average of these errors on 100 such training-test pairs of datasets. The above steps are repeated for the choices of $p=5,10,25,50,75,100,125,150$.

Panel (a) of figure 3 compares prediction errors using both robust and maximum likelihood SDR estimates when the covariates contain no outliers: here the two methods are virtually indistinguishable. We then introduce outliers in each of the 100 datasets by adding 100 to first $p/5$ coordinates of the first 10 observed covariate values, and repeat the analysis. Panel (b) of the figure shows that the robust SDR method remains more accurate in predicting out of sample observations for all values of $p$ than the standard SDR.

8 Real Data Applications

We now present an application of our proposed approach to some real data problems.

Robust techniques are useful when in identifying outlying observations, and we illustrate below how to make use of our fixed-dimensional methods presented earlier for functional (and hence infinite-dimensional) data.

We follow the approach of Boente and Salibian-Barrera, (2015) for performing robust principal component analysis on functional data using the estimated eigenvectors from $\widehat{\tilde{\Sigma}}$. Suppose the data consistents of $n$ curves, say $\mathcal{F}=\{f_{1},\ldots,f_{n}\}\in L^{2}[0,1]$, each observed at a set of common design points $\{t_{1},...,t_{m}\}$. We model each of these functions as a linear combination of $p$ mutually orthogonal B-spline basis functions $\mathcal{D}=\{\delta_{1},...,\delta_{p}\}$. We map data for each of the functions onto the coordinate system formed by the spline basis:

[TABLE]

We then model the $i$-th row of the $n\times p$ matrix $T(\mathcal{F},\mathcal{D})\equiv T$, denoted by ${\mathbf{T}}_{i}$ as follows:

[TABLE]

where $\mu$ is a location parameter, $P$ is a $p\times q$ loading matrix, $s_{i}$ is a $q\times 1$ score vector, and $e_{i}$ is the error term. We obtain robust estimators of $\mu$, $P$ and consequently $s_{i}$ using $\widehat{\tilde{\Sigma}}$. Define $\widehat{{\mathbf{T}}}_{i}=\hat{\mu}+\hat{P}\hat{s}_{i}$. The orthogonal distance (OD) corresponding to this projection is defined as

[TABLE]

Analogously, the score distance (SD) is defined as

[TABLE]

where $\hat{\lambda}_{1},\ldots,\hat{\lambda}_{q}$ are the top eigenvalues from $\widehat{\tilde{\Sigma}}$. For outlier detection, following Hubert et al., (2005) we set the upper cutoff values for score distances at $(\chi^{2}_{2,.975})^{1/2}$ and orthogonal distances at

[TABLE]

where $\Phi(\cdot)$ is the standard normal cumulative distribution function.

We apply the above outlier detection method on two data sets. First, we consider the monthly average sea surface temperature anomaly data from June 1970 to May 2004, available from http://www.cpc.ncep.noaa.gov/data, depicted in top-left panel of Figure 4)). Second, we consider the octane data, which consists of 226 variables and 39 observations (Esbensen et al.,, 1994). Each sample is a gasoline compound with a certain octane number, and has its NIR absorbance spectra measured in 2 nm intervals between 1100 - 1550 nm. There are 6 outliers here: compounds 25, 26 and 36-39, which contain alcohol. This data is presented in top-right panel of Figure 4)).

In the sea surface temperature data, using a cubic spline basis with knots at alternate months starting in June, we get a close approximation as depicted in middle-left panel of Figure 4. Using our proposed methodology with $q=1$ results in two points having their SD and OD larger than cutoff, depicted in bottom-left panel of Figure 4. These points correspond to the time periods June 1982 to May 1983 and June 1997 to May 1998 are marked by black curves in panels a and c, and pinpoint the two seasons with strongest El-Niño events.

On the octane data, we use the same methodology, and again the top robust PC turns out to be sufficient in identifying all 6 outliers. Details are available in the right side panels of Figure 4.

9 Conclusions

We propose the use of a weighted multivariate sign transformation for robust estimation and inference, and as demonstrated by theoretical results and several simulation studies and data examples, in many situations using a data-depth driven weight function leads to considerable efficiency gain without compromising on robustness properties. Our methodology seems to suggest new ways of identifying El-Niño or La-Niña events from the sea-surface temperature anomaly data, which will be studied further later.

Several of our results stated above are for data from the Euclidean space ${\mathbb{R}}^{p}$, where $p$ is fixed. The cases where $p$ increases with sample size and may be higher than sample size, and where data are from a separable Hilbert space, will be considered in a future work. There are only few conceptual challenges to such extensions, however, there are several technical and algebraic challenges.

Acknowledgements

The research of SC is partially supported by the National Science Foundation (NSF) under grants # DMS-1622483, # DMS-1737918.

Supplementary material

Appendix A Form of $V_{W}$

First observe that for $F$ having covariance matrix $\Sigma=\Gamma\Lambda\Gamma^{T}$,

[TABLE]

where $V_{W,\Lambda}$ is the covariance matrix of ${\mathbb{F}}_{\Lambda}$, the elliptic distribution with mean $\mu$ and covariance matrix $\Lambda$. Now,

[TABLE]

The matrix $\operatorname*{vec}(\tilde{\Lambda})\operatorname*{vec}^{T}(\tilde{\Lambda})$ consists of elements $\lambda_{i}\lambda_{j}$ at $(i,j)^{\text{th}}$ position of the $(i,j)^{\text{th}}$ block, and 0 otherwise. These positions correspond to variance and covariance components of on-diagonal elements. For the expectation matrix, all its elements are of the form ${\mathbb{E}}[\sqrt{\lambda_{a}\lambda_{b}\lambda_{c}\lambda_{d}}Z_{a}Z_{b}Z_{c}Z_{d}.W^{4}(Z,{\mathbb{F}}_{Z})/(Z^{T}\Lambda Z)^{2}]$, with $1\leq a,b,c,d\leq p$. Since $W^{4}(Z,{\mathbb{F}}_{Z})/(Z^{T}\Lambda Z)^{2}$ is even in $Z$, which has a spherically symmetric distribution, all such expectations will be 0 unless $a,b,c,d$ are all equal or pairwise equal. Following a similar derivation for spatial sign covariance matrices in Magyar and Tyler, (2014), we collect the non-zero elements and write the matrix of expectations:

[TABLE]

where ${\mathbb{I}}_{k}=(e_{1},...,e_{k}),{\mathbb{K}}_{m,n}=\sum_{i=1}^{m}\sum_{j=1}^{n}{\mathbb{J}}_{ij}\otimes{\mathbb{J}}_{ij}^{T}$ with ${\mathbb{J}}_{ij}\in{\mathbb{R}}^{m\times n}$ having 1 as $(i,j)$-th element and 0 elsewhere, and $\tilde{\gamma}_{mn}={\mathbb{E}}[\lambda_{m}\lambda_{n}Z_{m}^{2}Z_{n}^{2}.W^{4}(Z,{\mathbb{F}}_{Z})/(Z^{T}\Lambda Z)^{2}];1\leq m,n\leq p$.

Putting everything together, denote by $\hat{\tilde{\Lambda}}$ the sample version of $\tilde{\Lambda}$, the weighted covariance matrix obtained from ${\mathbb{F}}_{\Lambda}$, i.e. $\hat{\tilde{\Lambda}}=\sum_{i=1}^{n}W_{n}^{2}(Z_{i},{\mathbb{F}}_{Z,n}){\mathbb{S}}(\Lambda^{1/2}Z_{i};\hat{\mu}_{n})/n$. Then the different types of elements in the matrix $\hat{\tilde{\Lambda}}$ are as given below ($1\leq a,b,c,d\leq p$):

• •

Variance of on-diagonal elements

[TABLE]

• •

Variance of off-diagonal elements ($a\neq b$)

[TABLE]

• •

Covariance of two on-diagonal elements ($a\neq b$)

[TABLE]

• •

Covariance of two off-diagonal elements ($a\neq b,c\neq d$)

[TABLE]

• •

Covariance of one off-diagonal and one on-diagonal element ($a\neq b\neq c$)

[TABLE]

The above give all the elements of $V_{W,\Lambda}$. We plug these in (17) to recover $V_{W}$.

Appendix B Proofs

Proof sketch of Corollary 2.4.

We start with slightly modified versions of Lemmas A.4 and A.5 in Wang et al., (2015):

Lemma B.1**.**

Given that condition 1 in Corollary 2.4 holds, we have $\lambda_{\max}(\Psi_{1W})\leq 2W_{\max}\frac{\lambda_{\max}(\Sigma)}{\text{Tr}(\Sigma)}(1+o(1))$.

Lemma B.2**.**

Define $\Psi_{3W}={\mathbb{E}}\left[\frac{W^{2}(\epsilon,{\mathbb{F}}_{\epsilon})}{|\epsilon|^{2}}({\mathbb{I}}_{p}-S(\epsilon)S(\epsilon)^{T})\right]$. Then $\lambda_{\max}(\Psi_{2W})\leq W_{\max}{\mathbb{E}}(|\epsilon|^{-1})$ and $\lambda_{\max}(\Psi_{3W})\leq(W_{\max})^{2}{\mathbb{E}}(|\epsilon|^{-2})$. Further, if conditions 1 and 2 of Corollary 2.4 hold then $\lambda_{\min}(\Psi_{2W})\geq{\mathbb{E}}(W(\epsilon,{\mathbb{F}}_{\epsilon})/|\epsilon|)(1+o(1))/\sqrt{3}$.

The lemmas can be proved using similar steps as the proofs of the original lemmas in Wang et al., (2015) and using the upper bound on the weight function. Corollary 2.4 is now proved by applying Corollary 2.2, plugging in upper bound of $\lambda_{\max}(\Psi_{1W})$ from Lemma B.1 and lower bound of $\lambda_{\max}(\Psi_{1W})$ from Lemma B.2 into the ARE expression. ∎

Proof of Theorem 3.4.

We suppose $G_{n}=(g_{1},\ldots,g_{p}),L_{n}=\operatorname*{diag}(l_{1},\ldots,l_{p})$. In spirit, this corollary is similar to Theorem 13.5.1 in Anderson, (2009). We start with the following result, due to (Taskinen et al.,, 2012), allows us to obtain asymptotic joint distributions of eigenvectors and eigenvalues of $\hat{\tilde{\Sigma}}$, provided we know the limiting distribution of $\hat{\tilde{\Sigma}}$ itself:

Theorem B.3**.**

Let ${\mathbb{F}}_{\Lambda}$ be defined as before, and $\hat{C}$ be any positive definite symmetric $p\times p$ matrix such that at $F_{\Lambda}$ the limiting distribution of $\sqrt{n}\operatorname*{vec}(\hat{C}-\Lambda)$ is a $p^{2}$-variate (singular) normal distribution with mean zero. Write the spectral decomposition of $\hat{C}$ as $\hat{C}=\hat{P}\hat{\Lambda}\hat{P}^{T}$. Then the limiting distributions of $\sqrt{n}\operatorname*{vec}(\hat{P}-{\mathbb{I}}_{p})$ and $\sqrt{n}\operatorname*{vec}(\hat{\Lambda}-\Lambda)$ are multivariate (singular) normal and

[TABLE]

The first matrix picks only off-diagonal elements of the left-hand side and the second one only diagonal elements. We shall now use this as well as the form of the asymptotic covariance matrix of the vectorized $\hat{\tilde{\Sigma}}$, i.e. $V_{W}$ to obtain limiting variance and covariances of eigenvalues and eigenvectors.

Due to the decomposition (18) we have, for ${\mathbb{F}}_{\Lambda}$, the following relation between any off-diagonal element of $\hat{\tilde{\Lambda}}$ and the corresponding element in the estimate of eigenvectors, say $\hat{\tilde{\Gamma}}_{\Lambda}$:

[TABLE]

So that for eigenvector estimates of the original ${\mathbb{F}}$ we have

[TABLE]

Now $\sqrt{n}(\hat{\tilde{\gamma}}_{\Lambda,i,i}-1)=o_{P}(1)$ and $A{\mathbb{V}}(\sqrt{n}\hat{\tilde{\Lambda}}(i,k),\sqrt{n}\hat{\tilde{\Lambda}}(i,l))=0$ for $k\neq l$, so the above equation implies

[TABLE]

For the covariance terms, from (19) we get, for $i\neq j$,

[TABLE]

The exact forms given in the statement of the corollary now follows from the Form of $V_{W}$ in Section A.

For the on-diagonal elements of $\hat{\tilde{\Lambda}}$, using Theorem B.3 we have for the $i^{\text{th}}$ eigenvalue of $\hat{\tilde{\Lambda}}$, say $\lambda_{\Lambda,i}$,

[TABLE]

for $i=1,...,p$. Hence

[TABLE]

A similar derivation gives the expression for $A{\mathbb{V}}(l_{i},l_{j});i\neq j$. Finally, since the asymptotic covariance between an on-diagonal and an off-diagonal element of $\hat{\tilde{\Lambda}}$, it follows that the elements of $G_{n}$ and diagonal elements of $L_{n}$ are independent. ∎