Regularized Tapered Sample Covariance Matrix

Esa Ollila, Senior member, IEEE and Arnaud Breloy Abstract-Covariance matrix tapers have a long history in signal processing and related fields. Examples of applications include autoregressive models (promoting a banded structure) or beamforming (widening the spectral null width associated with an interferer). In this paper, the focus is on high-dimensional setting where the dimension p is high, while the data aspect ratio n/p is low. We propose an estimator called TABASCO (TApered or BAnded Shrinkage COvariance matrix) that shrinks the tapered sample covariance matrix towards a scaled identity matrix. We derive optimal and estimated (data adaptive) regularization parameters that are designed to minimize the mean squared error (MSE) between the proposed shrinkage estimator and the true covariance matrix. These parameters are derived under the general assumption that the data is sampled from an unspecified elliptically symmetric distribution with finite 4th order moments (both real-and complex-valued cases are addressed). Simulation studies show that the proposed TABASCO outperforms all competing tapering covariance matrix estimators in diverse setups. An application to space-time adaptive processing (STAP) also illustrates the benefit of the proposed estimator in a practical signal processing setup.

Index Terms-sample covariance matrix, shrinkage, regularization, elliptically symmetric distributions, tapering, banding, sphericity.

I. INTRODUCTION

Consider a set of p-dimensional (real-valued) vectors {x i } n i=1 sampled from a distribution of a random vector x with unknown mean vector µ = E[x] and unknown positive definite symmetric p × p covariance matrix Σ ≡ cov(x) = E[(x -µ)(xµ) ⊤ ]. In the high-dimensional case and when the sample size n is of the same order as p (p = O(n)) or p ≫ n, one is required to use regularization (shrinkage) in order to improve the estimation accuracy of the sample covariance matrix (SCM) and to obtain a positive definite matrix estimate. A popular estimate of Σ in such a setting is the regularized SCM (RSCM), defined by

S β = βS + (1 -β) tr(S) p I, (1)
where β ∈ [0, 1] is the regularization (or shrinkage) parameter, and where

S = 1 n -1 n i=1 (x i -x)(x i -x) ⊤ , (2)
denotes the unbiased sample covariance matrix (SCM), i.e., E[S] = Σ. Above tr(⋅) denotes the matrix trace, defined as tr(A) = ∑ i a ii for all square matrices A = (a ij ).

E. Ollila is with the Department of Signal Processing and Acoustics, Aalto University, P.O. Box 15400, FI-0007 Aalto, Finland.

Arnaud Breloy is with the LEME (EA4416), University Paris Nanterrre, 92410 Ville-d'Avray, France.

The code used in the presented experiments is available in the Matlab toolbox https://github.com/esollila/Tabasco. Note also that in (2), x = 1 n ∑ n i=1 x i denotes the sample mean vector. Automatic data-adaptive computation of optimal (oracle) parameter β for which S β in (1) attains the minimum mean squared error (MMSE) in Frobenius norm has been an active area of research. See for example [START_REF] Du | Fully automatic computation of diagonal loading levels for robust adaptive beamforming[END_REF], [START_REF] Ledoit | A well-conditioned estimator for largedimensional covariance matrices[END_REF], [START_REF] Chen | Shrinkage algorithms for MMSE covariance estimation[END_REF], [START_REF] Ollila | Optimal shrinkage covariance matrix estimation under random sampling from elliptical distributions[END_REF] to name only a few.

In many applications, the estimation accuracy (or another performance criterion) can alternatively be improved by using a so-called tapered SCM. Such estimate is defined as W • S, where • denotes the Hadamard (or Schur) elementwise product (i.e., (W • S) ij = w ij s ij for (W) ij = w ij and (S) ij = s ij ), and where W is a tapering matrix (also referred to as covariance matrix taper), i.e., a template that imposes some additional structure to the SCM.

Covariance matrix tapers have been used in many applications in diverse fields. A first main example in statistics is related to covariance matrices with a diagonally dominant structure (e.g., in autoregressive models). This means that the variables have a natural order in the sense that |i -j| large implies that the correlation between the ith and the jth variables is close to zero. In this settings, popular estimation approaches are to use a banding-type tapering matrices such as thresholding [START_REF] Bickel | Regularized estimation of large covariance matrices[END_REF], [START_REF]Covariance regularization by thresholding[END_REF]:

(W) ij = 1, |i -j| < k 0, |i -j| ≥ k (3)
for some integer k ∈ [ [1, p]] (called the bandwidth parameter), or softer thresholding variants. Notably, the strong theoretical merits of a linear decay of the form

(W) ij = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 1, |i -j| ≤ k/2 2 -2 |i -j| k , k/2 < |i -j| < k 0, |i -j| ≥ k (4)
were studied in [START_REF] Cai | Optimal rates of convergence for covariance matrix estimation[END_REF]. A second major example concerns the signal processing literature, in which tapering matrices have been developed in order to improve several spectral properties of adaptive beamformers, or to compensate subspace leakage and calibration issues [START_REF] Guerci | Principal components, covariance matrix tapers, and the subspace leakage problem[END_REF]. Most notably, the tapering matrices of the form (W) ij = sinc((ij)∆/π)

where ∆ ∈ R + , attracted interest as a null broadening technique for fluctuating interference [START_REF] Mailloux | Covariance matrix augmentation to produce adaptive array pattern troughs[END_REF], [START_REF] Zatman | Production of adaptive array troughs by dispersion synthesis[END_REF], [START_REF] Guerci | Theory and application of covariance matrix tapers for robust adaptive beamforming[END_REF], [START_REF] Song | Null broadening with snapshot-deficient covariance matrices in passive sonar[END_REF], [START_REF] Rugini | Regularized MMSE multiuser detection using covariance matrix tapering[END_REF].

A first approach to combine regularization with tapering was proposed in [START_REF] Chen | Shrinkage-to-tapering estimation of large covariance matrices[END_REF] with the shrinkage to tapering (ST) estimator, defined as the convex combination of the SCM and the tapered SCM:

S ST,β = βS + (1 -β)(W • S), (6)
where β ∈ [0, 1] is a shrinkage parameter. The authors then derived the optimal oracle parameter β o minimizing the MSE E[∥S ST,β -Σ∥ 2 F ], and proposed a shrinkage to tapering oracle approximating (STOA) estimator βo of β o under the assumption of Gaussian data. Authors in [START_REF] Li | Estimation of large covariance matrices by shrinking to structured target in normal and non-normal distributions[END_REF] also studied the ST estimator and derived an alternative oracle estimator of the shrinkage parameter both under Gaussian and non-Gaussian data. Data adaptive selection of the bandwidth k in (3) was also addressed with cross validation [START_REF] Chen | Shrinkage-to-tapering estimation of large covariance matrices[END_REF] or oracle estimation [START_REF] Li | Estimation of large covariance matrices by shrinking to structured target in normal and non-normal distributions[END_REF]. A possible issue with the ST estimate is that it inherently destroys the tapering template structure (e.g., sparsity for banded matrices) since it can be expressed as the modified tapered SCM S ST,β = (β11

⊤ + (1 -β)W) • S.
Hence, shrinkage is applied to the tapering matrix itself rather than to the SCM. In the high dimensional case, it should also be noted that both W•S and S ST,β are not necessarily positive semidefinite matrices, i.e., they can have negative or null eigenvalues. A possible solution for this problem is to compute their eigenvalue decomposition (EVD) and then replacing the invalid eigenvalues by small positive constants. However, such a post-processing step further deteriorates the template pattern of the covariance matrix estimator, and is computationally restrictive when dealing with high-dimensional data.

In this paper we provide a solution to the aforementioned problems by jointly leveraging shrinkage to identity and tapering: Let W = {W(k)} K k=1 be a finite set of possible tapering matrices 1 satisfying

W(k) ∈ W + ∀k ∈ [[1, K]], with W + = {W ∈ R p×p Sym ∶ w ii = 1, w ij ≥ 0 ∀i, j ∈ [[1, p]]} (7)
and with R p×p Sym denoting the set of all symmetric p×p matrices and [[1, p]] = {1, . . . , p}. We propose an estimator, referred to as TABASCO (TApered or BAnded Shrinkage COvariance matrix), defined as

Σβ,k = β(W(k) • S) + (1 -β) tr(S) p I, (8)
which benefits both from shrinkage (as the classic estimator in (1)) and exploitation of structure via tapering. Note that it also preserves the original scale of the SCM since tr(W • S) = tr(S) ∀W ∈ W + . Obviously, the success of banding and/or tapering depends on one's ability to choose the parameters β and k correctly. In this scope, we derive a fully automatic data-adaptive evaluation of the optimal parameters that jointly minimize the mean squared error

E[∥ Σβ,k -Σ∥ 2 F ]
under the general assumption that the data is sampled from an unspecified elliptically symmetric (ES) distribution with finite 4th order moments. A main interest to consider the general ES model is that it encompasses the standard Gaussian one while still accounting for possibly heavy-tailed distributions. Thus this assumption yields robustness to a large class of possible underlying data distributions. Our empirical experiments evidence that the proposed approach offers a near-tooptimal regularization parameter selection which outperform cross-validation schemes (especially at low sample support). 1 In this paper, we mostly focus on k implying a notion of bandwidth (or model order), for which W can be constructed from (3) or (4) with k ∈ [ [1, p]]. However, the proposed methodology applies to the general setting where W corresponds to any finite collection of possibly envisioned templates. Notably, we will also consider an application where k indexes a set of possible {∆ k } K k=1 used for the template model in [START_REF] Bickel | Regularized estimation of large covariance matrices[END_REF].

Since both the RSCM in (1) (if W = 11 ⊤ ∈ W) and the tapered SCM (β = 1) appear as special cases of (8), TABASCO performs never worse than these two estimators in terms of MSE independently of the underlying structure of the true covariance matrix.

The paper is structured as follows. In Section II expressions for the oracle regularization parameters β and k that minimize the MSE are derived in the general case of sampling from an unspecified distribution with finite 4th-order moments. In Section III a practical closed-form expression for the optimal regularization parameters are derived when sampling from an unspecified ES distribution, and an adaptive fully automatic procedure for their computation is proposed. It is shown that the optimal parameters depend on the sphericity of the tapered covariance matrix W • Σ, and we address the estimation of this quantity in Section IV. Section V extends our results for complex-valued observations. The special cases of known location (µ = 0) is also briefly discussed in Appendix A. Section VI provides simulation studies while in Section VII the estimator is applied to STAP on a real dataset. Finally, Section VIII concludes. The Appendix contains more technical proofs.

II. ORACLE TABASCO PARAMETERS β AND k

Recall that the TABASCO estimator Σβ,k is defined by [START_REF] Guerci | Principal components, covariance matrix tapers, and the subspace leakage problem[END_REF] for a set W = {W(k)} K k=1 of envisioned tapering matrices (cf. footnote 1 for examples) and a regularization parameter β ∈ [0, 1]. In this section, we derive the expression of the oracle parameters β and k that minimize the MSE in the general case of sampling from an unspecified p-variate distribution with finite 4th-order moments.

A. Oracle shrinkage parameter β for fixed k

First, let us introduce some notations and statistical parameters that are elemental in the proposed method. The scale is defined as the mean of the eigenvalues, and denoted by

η = tr(Σ) p (9)
while the sphericity of Σ [START_REF] Ledoit | Some hypothesis tests for the covariance matrix when the dimension is large compared to the sample size[END_REF], [START_REF] Srivastava | Some tests concerning the covariance matrix in high dimensional data[END_REF] is defined as

γ ≡ γ(Σ) = p tr(Σ 2 ) tr(Σ) 2 = p∥Σ∥ 2 F tr(Σ) 2 , ( 10
)
where ∥⋅∥ F denotes the Frobenius matrix norm, i.e., ∥A∥ 2 F = tr(A ⊤ A). The sphericity measures how close Σ is to a scaled identity matrix: γ ∈ [1, p], where γ = 1 if and only if Σ ∝ I and γ = p if and only if Σ has its rank equal to 1. For any W ∈ W + as in [START_REF] Cai | Optimal rates of convergence for covariance matrix estimation[END_REF], the matrix W • Σ, is called the tapered covariance matrix and we denote by

γ W ≡ γ(W • Σ) = p∥W • Σ∥ 2 F tr(Σ) 2 , (11)
the sphericity parameter of the tapered covariance matrix. Note that in [START_REF] Guerci | Theory and application of covariance matrix tapers for robust adaptive beamforming[END_REF] we utilised the fact that for any p × p matrix A and any W ∈ W + , it holds that tr(W • A) = tr(A). When W = 11

⊤ , we write γ 11 ⊤ ≡ γ for brevity.

We start by assuming that the index k is fixed. This allows us to simply denote the fixed tapering matrix W ≡ W(k) and TABASCO as Σβ ≡ Σβ,k . To find the oracle MMSE shrinkage parameter β ∈ [0, 1] of Σβ , the aim is thus to solve

β o = arg min β∈[0,1] E ∥ Σβ -Σ∥ 2 F . (12)
Notice that the MSE of the tapered SCM is

MSE(W • S) = E ∥W • S -Σ∥ 2 F = E ∥W • S∥ 2 F + ∥Σ∥ 2 F -2∥V • Σ∥ 2 F , (13)
where

V = (v ij ) p×p with v ij = w ij for W ∈ W + . ( 14
)
By normalized MSE (NMSE) we refer to the quantity

NMSE(W • S) = MSE(W • S)/∥Σ∥
β o = ∥V • Σ -ηI∥ 2 F E ∥W • S -ηI∥ 2 F (15) = p(γ V -1)η 2 E [∥ W • S ∥ 2 F -p -1 E[tr(S) 2 ] (16)
= (γ V -1) γ ⋅ NMSE(W • S) + 2γ V -γ -E[η 2 ]/η 2 (17)
where γ V is defined via [START_REF] Guerci | Theory and application of covariance matrix tapers for robust adaptive beamforming[END_REF] and η = tr(S)/p. Furthermore, the value of the MSE at the optimum is

MSE( Σβ o ) = E tr(S) 2 -tr(Σ) 2 p + ∥Σ∥ 2 F -∥V • Σ∥ 2 F + (1 -β 0 ) ∥V • Σ -ηI∥ 2 F . ( 18
)
Proof. The proof is postponed to Appendix B.

Notice that Theorem 1 also provides the MMSE shrinkage parameter β o for the RSCM S β in (1) since Σβ = S β when W = 11

⊤ . For the RSCM the optimal parameter is

β o = (γ -1) γ ⋅ NMSE(S) + γ -E[η 2 ]/η 2 , ( 19
)
where we used [START_REF] Srivastava | Some tests concerning the covariance matrix in high dimensional data[END_REF] and the facts that γ = γ V and W • S = S for W = 11 ⊤ . The MMSE of the RSCM utilizing the optimal shrinkage parameter in [START_REF] Ollila | Complex elliptically symmetric distributions: survey, new results and applications[END_REF] is

MSE(S β o ) = E tr(S) 2 -tr(Σ) 2 p + (1 -β 0 ) ∥Σ -ηI∥ 2 F ,
where we used [START_REF] Fang | Symmetric Multivariate and Related Distributions[END_REF] and that

V • Σ = Σ for V = 11 ⊤ .

B. Oracle index k

Notice that MSE( Σβ 0 ) in ( 18) implicitly depends on k through W ≡ W(k) and V defined in [START_REF] Chen | Shrinkage-to-tapering estimation of large covariance matrices[END_REF]. We further have the relation

NMSE( Σβ o ) = C - ∥V • Σ∥ 2 F ∥Σ∥ 2 F + (1 -β 0 ) ∥V • Σ -ηI∥ 2 F ∥Σ∥ 2 F = C - γ V γ + (1 -β 0 ) γ V -1 γ = C - 1 γ + β 0 (1 -γ V ) γ , ( 20
)
where C is a constant that is not dependent on k. Equation [START_REF] Muirhead | Aspects of Multivariate Statistical Theory[END_REF] then implies that minimizing the MSE with respect to k is equivalent to set

k o = arg min k β 0 (k)(1 -γ V (k)), (21)
where β o (k) is given by any of the expressions in ( 15)-( 17) and γ V (k) is defined via [START_REF] Guerci | Theory and application of covariance matrix tapers for robust adaptive beamforming[END_REF]. Note that we have made explicit the dependence of β 0 and γ V on k in [START_REF] Horn | Matrix Analysis[END_REF] for clarity of exposition. We use a brute force strategy to determine k 0 , i.e., compute

β 0 (k)(1 -γ V(k)
) for all indices k, and choose k 0 as the index that resulted in minimum value of the objective.

III. DATA ADAPTIVE TABASCO UNDER ES DISTRIBUTIONS

The oracles parameters found in the previous section depend on the true underlying data distribution and covariance matrix through various unknown quantities. A practical implementation of TABASCO thus requires their adaptive evaluation. Rather than resorting to a costly (and potentially inaccurate) cross-validation scheme, we will consider the general case where the data is sampled from an unspecified ES distribution [START_REF] Fang | Symmetric Multivariate and Related Distributions[END_REF], [START_REF] Ollila | Complex elliptically symmetric distributions: survey, new results and applications[END_REF] to refine the result of Theorem 1. In this setting, we show that the oracle parameter β eventually depend on only few scalar-valued statistics that can be accurately estimated.

A. ES distributions

Before stating the main results, we first recall some definitions and key results regarding ES distribution [START_REF] Fang | Symmetric Multivariate and Related Distributions[END_REF], [START_REF] Ollila | Complex elliptically symmetric distributions: survey, new results and applications[END_REF]. The probability density function of an elliptically distributed random vector, denoted by x ∼ E p (µ, Σ, g), is given by

f (x) = C p,g |Σ| -1/2 g (x -µ) ⊤ Σ -1 (x -µ) , (22)
where Σ denotes the positive definite symmetric covariance matrix parameter, µ is the mean vector, g ∶ R ≥0 → R >0 is the density generator, which is a fixed function that is independent of x, µ and Σ, and C p,g is a normalizing constant ensuring that f (x) integrates to 1. Note that here we define g such that "scatter matrix" parameter Σ coincides with the covariance matrix. This can always be assumed (under assumption of finite 2nd order moments) without any loss of generality [START_REF] Fang | Symmetric Multivariate and Related Distributions[END_REF], [START_REF] Ollila | Complex elliptically symmetric distributions: survey, new results and applications[END_REF]. For example, the multivariate normal (MVN) distribution, denoted by N p (µ, Σ), is obtained when g(t) = exp(-t/2).

The flexibility regarding the density generator g allows for modeling a large class of distributions, including heavy-tailed ones such as the multivariate t-distribution (MVT) with ν > 2 degrees of freedom (d.o.f.), denoted by x ∼ t ν (µ, Σ), where ν > 2 needs to be assumed for finite 2nd-order moments and ν > 4 for finite 4th-order moments. The density generator in this case is

g(t) = (1 + t/(ν -2)) -ν+p
2 . The elliptical kurtosis [START_REF] Muirhead | Aspects of Multivariate Statistical Theory[END_REF] parameter κ is defined as

κ = E[r 4 ] p(p + 2) -1 = 1 3 kurt(x i ), (23)
where the expectation is over the distribution of the random variable r = ∥Σ -1/2 (x-µ)∥ and kurt(x i ) denotes the excess kurtosis of any (e.g., ith) marginal variable of x, defined as

kurt(x i ) = E[(x i -µ i ) 4 ] (E[(x i -µ i ) 2 ]) 2 -3,
where

µ i = E[x i ]. Furthermore, observe that E[r 2 ] = p.
The elliptical kurtosis parameter vanishes (so κ = 0) when x has a MVN distribution. We also recall from [4, Lemma 2] that

E ∥S∥ 2 F = (1 + τ 1 + τ 2 ) ∥Σ∥ 2 F + τ 1 tr(Σ) 2 , (24)
E tr(S)

2 = 2τ 1 ∥Σ∥ 2 F + 1 + τ 2 tr(Σ) 2 , ( 25
)
where the scalars

τ 1 = 1 n -1 + κ n and τ 2 = κ n (26)
are dependent on the elliptical distribution (and hence on the density generator g) only via its kurtosis parameter.

B. Oracle shrinkage parameter β under ES distributions

As in the previous derivations, we assume that the index k is fixed and simply denote the fixed tapering matrix W ≡ W(k). This section refines the result of Theorem 1 when assuming that the data is sampled from an unspecified ES distribution. We first remark that the optimal β in ( 16) involves the quantity E ∥W • S∥ 2 F , which can be specified thanks to derivation of the following lemma. Lemma 1. Let {x i } n i=1 be an i.i.d. random sample from E p (µ, Σ, g) with finite 4th order moments. Then for any W ∈ W + , it holds that

E ∥W • S∥ 2 F = (1 + τ 1 + τ 2 )∥W • Σ∥ 2 F + τ 1 tr((D Σ W) 2 )
and

E tr((D S W) 2 ) = 2τ 1 ∥W • Σ∥ 2 F + (1 + τ 2 ) tr((D Σ W) 2 ),
where

D Σ = diag(Σ) and D S = diag(S). 2 Proof. The proof is postponed to Appendix C. Note that if W = 11 ⊤ , then tr((D Σ W) 2 ) = tr(Σ) 2 and
W • S = S so the expectations in Lemma 1 coincide with those of [4, Lemma 2] (i.e., [START_REF] Zou | Multivariate sign-based highdimensional tests for sphericity[END_REF] and ( 25)).

Using Lemma 1 we may now derive a simpler closed form expression of the optimal shrinkage parameter β o given in Theorem 1 that depends only on few summary (scalar-valued) statistics. Let us denote

θ W = d ⊤ Σ (W • W)d Σ p 2 = tr((D Σ W) 2 ) p 2 , ( 27
)
where d Σ = (σ

β o = t t + (n/(n -1))(pθ W /η 2 + γ W -2γ/p) + κ ⋅ A , (28)
where t = n(γ V -1) and

A = pθ W /η 2 -1 + 2γ W -2γ/p. (29)
Proof. Follows from Theorem 1 after substituting the values of E ∥W • S∥ 2 F given in Lemma 1 and of E tr(S) 2 given in ( 25) into the denominator of β o in ( 16) and simplifying the expression.

C. Data adaptive TABASCO implementation

Following from Theorem 2, the proposed data-adaptive implementation of TABASCO consists in applying the oracle procedure of Section II by replacing each of the unknown parameters {η, θ W , κ, γ, γ W , γ V } in the optimal β 0 from (28) by carefully chosen estimates:

• For η and θ W , we use the empirical estimates:

η = tr(S)/p and θW = tr((D S W) 2 )/p 2 (30)
• The elliptical kurtosis κ can be estimated using an estimator κ detailed in [4, Sect. IV] as (bias-corrected) average sample excess kurtosis of the marginal variables scaled by 1/3. Also note that if the data is assumed to follow the MVN distribution, we can set κ = 0, and the last term κ ⋅ A can be ignored in the denominator. • The estimation of the three sphericity statistics: γ, γ W , and γ V is more intricate. It will be addressed in detail in Section IV, which will propose two estimation methods (denoted Ell1 and Ell2) for these quantities. Also notice that V = ( w ij ) p×p , so if W is a selection matrix (i.e., that has only 0-s or 1-s as its off-diagonal elements), as for example in (3), then W = V so only γ W needs to be estimated. Using these plug-in values yields an estimate βo (k) for each template in the set W = {W(k)} K k=1 . Similarly, the index k is estimated based on ( 21) by replacing the unknown β 0 (k) and γ V (k) by their estimates and solving

ko = arg min k β0 (k)(1 -γV (k)), (31)
which allows then to compute a fully data-adaptive TABASCO estimator Σ β0 , k0 as in [START_REF] Guerci | Principal components, covariance matrix tapers, and the subspace leakage problem[END_REF]. The pseudocode of the proposed estimation algorithm is summarized in Algorithm 1.

Algorithm 1: TABASCO

Input : Data {x i } n i=1 , templates set {W(k)} K k=1 1 Compute SCM S in (2). 2 Compute η from (30) 3 Compute κ from [4, Sect. IV] 4 Compute γ (options in Section IV) 5 for k ∈ [[1, K]] do 6 Set W = W(k) and V = V(k) = ( w ij (k)) p×p 7
Compute θW from (30)

8
Compute γW (k) and γV (k) (options in Section IV)

9
Compute βo (k) from ( 28) using plug-in estimates 10 Select optimal k 0 as in [START_REF] Ollila | Shrinking the eigenvalues of mestimators of covariance matrix[END_REF] with

{ βo (k), γV (k)} K k=1 11 Set W = W( ko ) and β = βo (k) Output : Σ = β ⋅ (W • S) + (1 -β)ηI
IV. ESTIMATORS OF SPHERICITY In this section, we detail two new alternative estimators of the sphericity of the tapered covariance matrix W •Σ in order to compute the parameters {γ, γ W , γ V } in TABASCO. These estimators are nontrivial extensions of the (non-tapered) Ell1and Ell2-sphericity estimators in [START_REF] Ollila | Optimal shrinkage covariance matrix estimation under random sampling from elliptical distributions[END_REF]. Notation "Ell" simply refers to the fact that both estimators assume that the data is drawn from an unspecified elliptically symmetric distribution.

A. Ell1-estimator of sphericity

The Ell1-estimator is based on the spatial sign covariance matrix (SSCM) [START_REF] Marden | Some robust estimates of principal components[END_REF], [START_REF] Visuri | Sign and rank covariance matrices[END_REF], which has been popular for constructing robust estimates of sphericity [START_REF] Zou | Multivariate sign-based highdimensional tests for sphericity[END_REF], [START_REF] Zhang | Automatic diagonal loading for Tyler's robust covariance estimator[END_REF]. The robust properties of SSCM comes from the fact that it is distributionfree under elliptical models and has the highest possible breakdown point [START_REF] Magyar | The asymptotic inadmissibility of the spatial sign covariance matrix for elliptically symmetric distributions[END_REF], [START_REF] Croux | The k-step spatial sign covariance matrix[END_REF]. The Ell1-estimator was theoretically studied in [START_REF] Raninen | Linear pooling of sample covariance matrices[END_REF] and we propose here its generalization to the sphericity of the tapered covariance matrix W • Σ.

First, define the shape matrix (or normalized covariance matrix) as Λ = p Σ tr(Σ) and note that tr(Λ) = p. The sphericity measures γ and γ W for any W ∈ W + can then be expressed simply in terms of Λ via the formulas:

γ = ∥Λ∥ 2 F p and γ W = ∥W • Λ∥ 2 F p .
The (scaled) SSCM is defined by

Λ = p n n i=1 (x i -μ)(x i -μ) ⊤ ∥x i -μ∥ 2 , ( 32
)
where μ = arg min µ ∑ n i=1 ∥x i -µ∥ is the sample spatial median [START_REF] Brown | Statistical Uses of the Spatial Median[END_REF]. The sample mean could also be used, but [START_REF] Dürre | The spatial sign covariance matrix with unknown location[END_REF] concludes that the use of spatial median as the location estimator for SSCM is clearly preferable. When µ is known (and without loss of generality assuming µ = 0), the SSCM is defined as

Λ = p n ∑ n i=1 x i x ⊤ i ∥x i ∥ 2 .
Recently, it was shown in [START_REF] Raninen | Linear pooling of sample covariance matrices[END_REF] that the following estimate of sphericity based on the SSCM (when µ is known),

γ = n n -1 ∥ Λ∥ 2 F p - p n , (33)
is asymptotically (as p → ∞) unbiased when sampling from elliptical distributions under the following assumption (A) The sequence of covariance matrix structures being considered with increasing p satisfies γ = o(p) as p → ∞. In other words, E[γ] → γ as p → ∞ when (A) holds. We note that Assumption (A) is sufficiently general and holds for many covariance matrix models as shown in [START_REF] Raninen | Linear pooling of sample covariance matrices[END_REF]Prop. 3]. The following Theorem presents a modification of the Ell1estimator [START_REF] Ollila | Optimal shrinkage covariance matrix estimation under random sampling from elliptical distributions[END_REF] for the sphericity of W • Σ with equivalent asymptotic guarantees. Theorem 3. Let {x i } n i=1 be an i.i.d. random sample from an ES distribution E p (µ, Σ, g) with known µ = 0. Then, for any W ∈ W + and under Assumption (A), the following statistic

γW = n n -1 ⎛ ⎜ ⎝ ∥W • Λ∥ 2 F p - tr (D ΛW) 2 np ⎞ ⎟ ⎠ , (34)
where

D Λ = diag( Λ), is asymptotically, as p → ∞, unbiased estimator of γ W = γ(W • Σ) in (11), i.e., E[γ W ] → γ W as p → ∞, for any fixed n.
Proof. Proof is postponed to the Appendix D.

Also observe that when W = 11 ⊤ , then γW reduces to γ in [START_REF] Wagaman | Discovering sparse covariance structures with the isomap[END_REF].

When the mean is unknown and estimated by spatial median μ, then the results in [START_REF] Zou | Multivariate sign-based highdimensional tests for sphericity[END_REF] derived in the case of W = 11 ⊤ show that when p = O(n 2 ), the estimator of sphericity γ needs to be corrected for bias. Thus the reader should bear in mind that γ (and γW ) may have a non-negligible bias when µ is estimated by μ and p is orders of magnitude larger than n.

B. Ell2-estimator of sphericity

The Ell2-estimator of sphericity was proposed in [START_REF] Ollila | Optimal shrinkage covariance matrix estimation under random sampling from elliptical distributions[END_REF] and we derive here its adaptation to the sphericity of the tapered covariance matrix W • Σ.

The derivation starts by noticing that the obvious plug-in estimate ∥W • S∥ 2 F /p for the parameter

ϑ W = ∥W • Σ∥ 2 F p (35)
is biased by resorting to for E[∥W • S∥

W ∈ W + is θW = b n ⎛ ⎜ ⎝ ∥W • S∥ 2 F p -a n tr (D S W) 2 p ⎞ ⎟ ⎠ ,
where

a n = 1 n + κ n n -1 + κ (36) b n = (κ + n)(n -1) 2 (n -2)(3κ(n -1) + n(n + 1)) . ( 37
)
Proof. Note that a n in (36) can be written as a n = τ 1 /(1+τ 2 ), where definitions of τ 1 and τ 2 are given by ( 26) while b n in (37) can be expressed as b n = (τ 1 +τ 2 -2τ 1 a n -1 . Then using Lemma 1, we notice that b -1

n pE[ θW ] = (τ 1 -a n (1 + τ 2 )) tr((D Σ W) 2 ) + (1 + τ 1 + τ 2 -2τ 1 a n ) ∥W • Σ∥ 2 F = b -1 n ∥W • Σ∥ 2 F
The expressions [START_REF] Kraut | The cfar adaptive subspace detector is a scale-invariant glrt[END_REF] and ( 37) are obtained when replacing the values of τ 1 and τ 2 given in ( 26) into a n ≡ a n (τ 1 , τ 2 ) and b n ≡ b n (τ 1 , τ 2 ) and simplifying the obtained expressions.

Then note that the sphericity of the tapered covariance matrix can also be written as

γ W = ϑ W /η 2 ,
where ϑ W and η are defined in ( 35) and ( 10) respectively. Using this expression, we consider the estimate where θW is computed from Theorem 4, and η2 is obtained from [START_REF] Dürre | The spatial sign covariance matrix with unknown location[END_REF]. This yields the estimator

γW = p bn ⎛ ⎜ ⎝ ∥W • S∥ 2 F tr(S) 2 -ân tr (D S W) 2 tr(S) 2 ⎞ ⎟ ⎠ , (38)
where ân ≡ a n (κ) and bn ≡ b n (κ) are obtained by replacing the unknown κ in ( 36) and ( 37) by its estimate κ [4, Sect. IV].

We refer to [START_REF] Magnus | Matrix Differential Calculus with Applications in Statistics and Econometrics[END_REF] as Ell2-estimator of sphericity γ W . Also note that, if n is reasonably large, then bn ≈ 1 and n/(n + κ) ≈ 1, its expression can be simplified to

γW ≈ p ∥W • S∥ 2 F tr(S) 2 -(1 + κ) p n tr((D S W) 2 )
tr(S) 2 .

In the non-tapered case (W = 11 ⊤ ), the estimator in (38) reduces to the Ell2-estimator of sphericity in [START_REF] Ollila | Optimal shrinkage covariance matrix estimation under random sampling from elliptical distributions[END_REF].

C. Remarks

Although Ell2-estimator of sphericity does not require knowledge of the underlying elliptically symmetric distribution of the data, it is not a robust estimator. Thus we overall favour Ell1-estimator due to robustness of SSCM, and recommend usage of Ell2-estimator when dealing with data that is not heavy-tailed, i.e., which can be approximated by a Gaussian distribution. The non-robustness of Ell2-estimator is due to its usage of 4th-order moments. Namely, tr(S

2 ) = ∑ p i=1 ∑ p j=1 s 2 ij where s ij = (S) ij = 1 n ∑ n =1
x i x j -xi xj , can be written as a sum of mixed 4th-order sample moments. By central limit theorem, any 4th-order sample moment has a limiting normal distribution if the moments of x exists up to 8th order. Without this condition, the estimator tr(S 2 )/p can be highly-variable. Ell1-estimator is highly robust and performs well for heavier-tailed data. Yet, the assumption γ = o(p) needed by Ell1-estimator means that γ ≪ p when dimension is large. Implicitly, this often means that the eigenvalues of Σ tend to be similar as p grows. In [START_REF] Raninen | Linear pooling of sample covariance matrices[END_REF]Theorem 2] it was shown that the bias of the SSCM Λ is of the order of the sphericity, γ = tr(Λ 2 )/p, and becomes negligible when γ = o(p).

Thus for distinctively non-spherical covariance matrices, Ell1estimator has a visible bias, but the bias vanishes for large p (since the estimator is asymptotically unbiased as shown in Theorem 3 when γ = o(p)). Thus Ell2-estimator can be favoured when the data is approximately Gaussian, e.g., when the estimated kurtosis is smaller than 1. For example, for MVT distribution with ν = 10, kurt(x i ) = 1.0. Finally, we mention that in practice, we also always use the thresholding

γ = min(p, max(1, γ)) (39)
for any option in order to guarantee that the final estimator remain in the valid interval, 1 ≤ γ ≤ p.

V. EXTENSION TO THE COMPLEX-VALUED CASE

A. Complex elliptically symmetric (CES) distribution

First we recall some definitions and notations specific to complex-valued case. By ∥x∥ We now assume that the data {x i } n i=1 is a random sample from a circular complex elliptically symmetric (CES) distribution, denoted x ∼ CE p (µ, Σ, g) (cf. [START_REF] Ollila | Complex elliptically symmetric distributions: survey, new results and applications[END_REF] for a detailed review). Similarly to the real-valued case, the probability density function of a CES distributed random vector x ∈ C p is given by

f (x) = C p,g |Σ| -1 g((x -µ) H Σ -1 (x -µ)),
where Σ denotes the positive definite Hermitian covariance matrix, µ = E[x] is the mean vector, g ∶ R ≥0 → R >0 is the density generator, and C p,g is a normalizing constant. Again, we also normalize g so that

Σ = E[(x -µ)(x -µ) H ].
The definitions of the scale and sphericity parameters in [START_REF] Zatman | Production of adaptive array troughs by dispersion synthesis[END_REF] and [START_REF] Guerci | Theory and application of covariance matrix tapers for robust adaptive beamforming[END_REF] remain unchanged. The elliptical kurtosis is however redefined as

κ = E[r 4 ] p(p + 1) -1 = 1 2 kurt(x i ).
where the expectation is over r = ∥Σ -1/2 (xµ)∥ and kurt(x i ) denotes the excess kurtosis of any (e.g., ith) marginal variable of x, defined by kurt

(x i ) = E[|x i -µ i | 4 ] σ 4 i -2,
where

µ i = E[x i ] and σ 2 i = E[|x i -µ i | 2 ]
denote the mean and variance of x i . The theoretical lower bound of the kurtosis in the complex-valued case is κ LB = -1/(p + 1) [START_REF] Ollila | Complex elliptically symmetric distributions: survey, new results and applications[END_REF]. Again κ = 0 if x has a circular complex multivariate normal distribution (x ∼ CN p (µ, Σ)).

B. Oracle β and k under CES distributions

The SCM (2) of complex-valued observations is defined by

S = 1 n -1 n i=1 (x i -x)(x i -x) H (40)
and the complex-valued counterpart TABASCO Σβ is still defined as in [START_REF] Guerci | Principal components, covariance matrix tapers, and the subspace leakage problem[END_REF]. Extending the previous results to this estimator will require some minor adaptations.

We start by noticing that Theorem 1 still holds for complex-valued observations. Then, the next result provides the complex-valued extension of Lemma 1. Lemma 2. Let {x i } n i=1 be an i.i.d. random sample from CE p (µ, Σ, g) with finite 4th order moments. Then for any W ∈ W + , and for the SCM as in (40), it holds that

E ∥W • S∥ 2 F = (1 + τ 2 )∥W • Σ∥ 2 F + τ 1 tr((D Σ W) 2 )
and

E tr((D S W) 2 ) = 2τ 1 ∥W • Σ∥ 2 F + (1 + τ 2 ) tr((D Σ W) 2 ),
where D Σ = diag(Σ), D S = diag(S) and τ 1 and τ 2 are defined in [START_REF] Magyar | The asymptotic inadmissibility of the spatial sign covariance matrix for elliptically symmetric distributions[END_REF].

Proof. The proof is postponed to Appendix E.

This result allows us to derive the complex-valued counterpart of Theorem 2 for the optimal shrinkage parameter β o .

Theorem 5. Let {x i } n i=1 be an i.i.d. random sample from a complex elliptical distribution CE p (µ, Σ, g) with finite 4th order moments. Then the oracle parameter β o in (12) is

β o = t t + (n/(n -1))(pθ W /η 2 -2γ/p) + κ ⋅ A ,
where t = n(γ V -1), and

A = pθ W /η 2 -1 + γ W -2γ/p.
Form this result, the optimal index k can be obtained as in subsection II-B.

C. Data adaptive adaptations for the complex-valued case In the complex-valued case, the TABASCO procedure of Algorithm 1 is kept identical. However, the expression of β o is changed in accordance with Theorem 5. The plug in estimates for the side parameters {η, θ W , κ, γ, γ W , γ V } are then as follows:

• For η and θ W , we keep the empirical estimates as in [START_REF] Dürre | The spatial sign covariance matrix with unknown location[END_REF].

• An estimate of elliptical kurtosis κ is calculated as in [START_REF] Ollila | Shrinking the eigenvalues of mestimators of covariance matrix[END_REF] as average sample excess kurtosis of the marginal variables scaled by 1/2. • The three sphericity statistics γ, γ W , and γ V can still be estimated with the Ell1 and Ell2 estimators detailed in Section IV. The Ell1 estimator is identical, except that the definition of the SSCM in (32) now involves the Hermitian transpose. The Ell2 estimator is defined as earlier, with changes only in expression of a n and b n . Indeed, Theorem 4 holds with a n as in [START_REF] Kraut | The cfar adaptive subspace detector is a scale-invariant glrt[END_REF] and b n given by

b n = n(n -1) 2 (κ + n) 2κn(n 2 -4n + 3) -κ 2 (n -1) 2 + n 2 (n 2 -2n -1)
.

VI. SIMULATION STUDIES

We generate samples from (real-valued) ES distributions with a scatter matrix Σ having a diagonally dominant structure (model 1 and model 2 detailed below). The mean µ is generated randomly as N p (10 ⋅ 1, I) and kept fixed for all trials, and the number of Monte-Carlo trials is 5000.

The estimators included in the study are: i) The Ledoit-Wolf estimator (LWE) [START_REF] Ledoit | A well-conditioned estimator for largedimensional covariance matrices[END_REF] defined by [START_REF] Du | Fully automatic computation of diagonal loading levels for robust adaptive beamforming[END_REF] where β is an estimate of an (oracle) MMSE parameter β o . ii) The shrinkage to tapering oracle approximate (STOA) estimator [START_REF] Chen | Shrinkage-to-tapering estimation of large covariance matrices[END_REF] defined by [START_REF]Covariance regularization by thresholding[END_REF] where β is an estimate of the oracle parameter computed using an iterative procedure. The bandwidth k is selected using a cross-validation scheme with 60%-to-40% split for training and testing. iii) The shrinkage to tapering (ST-)estimators in [START_REF] Li | Estimation of large covariance matrices by shrinking to structured target in normal and non-normal distributions[END_REF] defined by [START_REF]Covariance regularization by thresholding[END_REF] where both β and k are estimates of the oracle MMSE parameters. The estimator ST-gaus assumes Gaussian data, while ST-nong assumes non-Gaussian data. iv) TABASCO (computed via Algorithm 1) using the Ell1estimator of sphericity.

A. Model 1

In Model 1, Σ possesses an auto-regressive AR(1) structure:

(Σ) ij = η |i-j| , (41)
where | | ∈ [0, 1). When ↓ 0, then Σ is close to an identity matrix scaled by η, and when ↑ 1, Σ tends to a singular matrix of rank 1. As illustrated in Figure 1, banding matrices allow for a good approximation, so all tapering-type estimators are computed with W(k) as ( 3) in this subsection. The optimal bandwidth ko is chosen by consider the set of tapering matrices [START_REF] Dürre | The spatial sign covariance matrix with unknown location[END_REF]] ∪ [[p -30, p]]} (this restriction is made to lower the computational cost, but identical results were obtained with k ∈ [[1, p]]).

W = {W(k) ∶ k ∈ [[1,
Figure 2 provides a validation of the theoretical results: it displays the theoretical normalized MSE (NMSE) curves,

L(β) = E[∥ Σβ -Σ∥ 2 F ]/∥Σ∥ 2 F
as a function of shrinkage parameter β for TABASCO estimators using a fixed bandwidths k ∈ [ [START_REF] Du | Fully automatic computation of diagonal loading levels for robust adaptive beamforming[END_REF][START_REF] Bickel | Regularized estimation of large covariance matrices[END_REF]] and k = p (i.e., W = 11 ⊤ ). In this setup, the data is generated from MVN distribution N p (µ, Σ) with p = 100 and n = 50 (similar results were obtained for other ES distributions and dimension setups). The black bullet (•) displays the theoretical minimum NMSE in [START_REF] Fang | Symmetric Multivariate and Related Distributions[END_REF] β o ≡ β o (k) for each bandwidth k. The empirical average NMSE for TABASCO using estimated β0 for each fixed k is displayed using red triangle (▲), where the location on β axis correspond to empirical average β0 . As can be noted from Figure 2, TABASCO estimates the oracle shrinkage parameter β o very accurately since the black bullets and red triangles are mostly overlapping for each bandwidth. The dashed horizontal line shows the average NMSE obtained by TABASCO when also using an estimated optimal bandwidth k0 . One can notice that the optimal bandwidth selection using ( 31) is also accurate. For example, in the case of = 0.4, the optimal bandwidth is k = 3 and TABASCO estimator attains an average NMSE that is very close to the theoretical minimum NMSE. Figure 3 displays the NMSE as a function of bandwidth k. As can be noted, the empirical average NMSE of TABASCO using fixed bandwidth displayed using triangle (e.g., ▲) coincides very accurately with the theoretical NMSE. The NMSE of TABASCO with estimated bandwidth is shown using colored star (e.g., ☆) wherein the location on k-axis correspond to empirical average k0 . As can be noted, the optimal bandwidth k o is also very accurately estimated since ko ≈ k o for each .

β k = 1 k = 2 k = 3 k = 4 k = 5 k = p
Figure 4 compares the performance of TABASCO with the state of the art in various setups. The upper panel displays the NMSE curves as a function of the sample size n for four choices of correlation parameter when the data follows a MVN distribution. The lower panel displays the same results when the data follows a MVT distribution with ν = 5, which is heavy-tailed with marginal kurtosis kurt(x i ) = 6 and elliptical kurtosis κ = kurt(x i )/3 = 2. In the Gaussian case, all banding-type estimators outperform LWE thanks to the exploitation of the diagonally dominant structure of the covariance matrix. In the heavy-tailed case, this is no longer true for STOA and ST-gaus, while ST-nong and TABASCO remain robust. In all scenarios, TABASCO offers the lowest NMSE, and especially improves the performance when n ≪ p. Figure 5 displays the obtained (average) estimated shrinkage parameter βo of TABASCO and LWE as a function of n. The average shrinkage parameter of TABASCO is generally much larger than that of LWE. This means that it assign overall more weight on the banded SCM W • S compared to LWE, which uses W(p) = 11

⊤ . This behavior is expected since banding the SCM should naturally improve the MSE when the true covariance matrix has a diagonally dominant structure.

Figure 6 presents a comparison similar to Figure 4 when the variables are permuted at random for each Monte Carlo trial, thus destroying the diagonally dominant structure of the AR(1) covariance matrix 3 . The hypothesis is that any banding estimator with optimal bandwidth selection should be able to select the bandwidth k = p accordingly. Note that LWE is invariant to variable permutations, and hence its results stays the same for both of these scenarios. In this setup, TABASCO performs better that LWE for n ≪ p and equally well as LWE for n large enough. This result implies that bandwidth selection of TABASCO is consistent: it chooses k = p since the true covariance matrix does not have a diagonally dominant structure. The improvement brought at low sample support can be explained by the fact that an ES distribution is assumed by TABASCO, which allows for a better estimation of the oracle parameter (LWE only assumes finite 4th order moments). This example confirms that TABASCO always benefits from banding and bandwidth selection: it offers significantly improved NMSE compared to RSCM when banding structure is present in the covariance matrix, while it does not perform worse when such structure does not exist, thanks to its robust and efficient bandwidth selection.

B. Model 2

In Model 2 [START_REF] Cai | Optimal rates of convergence for covariance matrix estimation[END_REF], Σ is defined by

(Σ) ij = 1 , i = j ρ|i -j| -(α+1) , i ≠ j, ( 42
)
where α is a decay parameter and ρ is a correlation parameter.

As in the study of [START_REF] Cai | Optimal rates of convergence for covariance matrix estimation[END_REF], we set ρ = 0.6, and Figure 1 illustrates the effect of decay parameter α in the case of p = 100. (oracle) bandwidth [START_REF] Cai | Optimal rates of convergence for covariance matrix estimation[END_REF]Section 6]. The dimension is p = 250. It should be noted that MnMx-Taper has advantage over the other estimators since it uses the true decay parameter α, which is unknown in practice. TABASCO also uses tapering matrices W(k) as in ( 4), but ST-gaus and ST-nong are restricted to tapering matrices whose off-diagonal elements are 0-s or 1-s.

Hence, these are still computed with banding matrices W(k) as in [START_REF] Chen | Shrinkage algorithms for MMSE covariance estimation[END_REF]. In either case, the optimal bandwidth ko is chosen by consider the set of tapering matrices

W = {W(k) ∶ k ∈ [[1, 30]] ∪ [[p -30, p]]}.
As can be noted, TABASCO again outperforms other estimators for all values of n and α and for both sampling distributions. In the MVN case (top panel), TABASCO outperforms MnMx-Taper with a clear margin when n is very small. This can be attributed to its ability to optimally shrink the tapered SCM towards a scaled identity matrix when n/p < 1. However for n ≥ p, TABASCO and MnMx-Taper estimator have similar performance, especially when α = 0.3.

In the MVT case (lower panel of Figure 7), the performance differences are more clear. TABASCO outperforms MnMxtaper by a large margin. ST-gaus estimator completely fails due to the impulsive nature of the underlying sampling distributions. The results also illustrate that the performance of tapered SCM estimator is dependent on the underlying sampling distribution more heavily than TABASCO. This is illustrated further in Figure 8 where we compare the true theoretical NMSE curves of tapered SCM W•S and TABASCO estimator Σβ 0 as a function of bandwidth k in the case where n = 100 and when sampling from a MVN distribution (left panel) and MVT distribution (right panel) with ν = 5 d.o.f. following model 2 with α = 0.1. Figure 8 shows two important points. First, the performance differences between the tapered SCM and TABASCO are larger when the distribution is heavier tailed, which was already evident in Figure 7. Second, TABASCO with optimal bandwidth selection is able to estimate the optimal bandwidth rather accurately since the average (empirical) NMSE value seen in Figure 7 at n = 100 is close to the minimum true (theoretical) NMSE value.

VII. APPLICATION TO SPACE-TIME ADAPTIVE PROCESSING

Space time adaptive processing (STAP) is a technique used in airborne phased array radar to detect moving target embedded in an interference background such as jamming or strong clutter [START_REF] Ward | Space time adaptive processing for airborne radar[END_REF]. The radar receiver consists in an array of Q antenna elements processing P pulses in a coherent processing interval. Within the tested sample x 0 ∈ C p with p = P ⋅ Q, the received signal is composed of i) possible unknown targets responses; ii) unknown interferences (ground clutter) plus thermal noise. A detection problem for a given steering vector p is classically formalized as a binary hypothesis test: under H 0 , x 0 only contains the interference plus noise, or under H 1 , x 0 additionally contains a scaled observation of p, i.e.:

H 0 ∶ x 0 = n 0 ; x i = n i , ∀ i ∈ [[1, n]] H 1 ∶ x 0 = αp + n 0 ; x i = n i , ∀ i ∈ [[1, n]]
where x i ∈ C p , i = 1, . . . , n is a secondary data set, assumed to contain i.i.d. and target-free realizations of the interference plus noise. Usually, this disturbance n i is modeled as centered complex Gaussian (or elliptically) distributed with covariance matrix Σ. In this context, efficient adaptive detection statistics can be built from the expression of the adaptive coherence estimator (ACE) detector [START_REF] Kraut | The adaptive coherence estimator: a uniformly most-powerful-invariant adaptive detection statistic[END_REF]:

Λ( Σ) = |p H Σ-1 x 0 | 2 |p H Σ-1 p||x H 0 Σ-1 x 0 | H 1 ≷ H 0 δ Σ, ( 43
)
where Σ is a plug-in estimate of Σ computed from {x i } n i=1 . More specifically in STAP, the target p follows the steering vector model of [START_REF] Ward | Space time adaptive processing for airborne radar[END_REF], which is function of the target angle of arrival (AoA) θ and velocity v. The statistic (43) can thus be computed for a dictionary of steering vectors covering a 2D-grid on θ and v, yielding an adaptive detection map.

Using the SCM as estimate in (43) yields a generalized likelihood ratio test (GLRT) [START_REF] Kraut | The cfar adaptive subspace detector is a scale-invariant glrt[END_REF]. However, plug-in detectors can benefit from refined estimation processes in order to improve robustness, or to deal with limited sample support issues. For example shrinkage to identity (also referred to as diagonal loading or robust beamforming [START_REF] Li | On robust Capon beamforming and diagonal loading[END_REF]) is a common procedure to improve several properties of the detector's output. In the context of interference cancellation, tapering templates have been considered as a spectrum notch-widening technique [START_REF] Guerci | Theory and application of covariance matrix tapers for robust adaptive beamforming[END_REF], or to deal with modulation effects [START_REF] Guerci | Principal components, covariance matrix tapers, and the subspace leakage problem[END_REF].

This section presents an experimental validation of TABASCO to illustrate the interest of both approach on real data. The STAP data is provided by the French agency DGA/MI: the clutter is real but the targets were synthetically added in the dataset conception. The number of sensors is Q = 4 and the number of coherent pulses is P = 64, the size of the data is then p = QP = 256. The center frequency and the bandwidth are respectively equal to f 0 = 10GHz and the bandwidth B = 5MHz. The radar celerity is V = 100m/s. The inter-element spacing is d = 0.3m and the pulse repetition frequency is f r = 1kHz. The clutter to noise ratio is evaluated around 20dB. We consider two different scenarios: one where the tested cell is under H 1 with a single target at (θ = 0, v = 4m/s), and one where the tested cell is under H 1 with 10 targets at various speed/angle. The signal to clutter ratio of each target was estimated to be around -5dB. In both scenarios, n = 397 (all available) target-free secondary data are used to estimate the interference covariance matrix.

The tapering matrix is constructed as proposed in [START_REF] Guerci | Theory and application of covariance matrix tapers for robust adaptive beamforming[END_REF] 4 , i.e.,

W(k) = T f ⊗ T θ T f ij = (1 + sinc((i -j)k/π))/2 ∈ R P ×P [T θ ] ij = (1 + sinc((i -j)k/π))/2 ∈ R Q×Q (44)
Note that index k is here a "null-spectrum width" parameter in R + and not a bandwidth parameter in [ [1, p]] as in (3) or (4). We also point out that the banding-type tapering matrices (even involving a Kronecker-product structure) were tested but appeared not well suited to the data, nor beneficial to the detection process. Thus they will not be discussed in the following.

Figure 9 presents the detection map of Λ( Σ) constructed with: i) the SCM; ii) the tapered SCM W(k) • S using bandwidth k = 0.05 (selected manually to obtain the best visual results); iii) TABASCO with the proposed adaptive selection of β for k = 0 (equivalent to RSCM, yielding β = 0.9324); iv) TABASCO with the proposed adaptive estimation of β and k allowing k ∈ [10 -3 , 10

-1 ] (TABASCO, yielding k = 0.0143 and β = 0.9929). ST-type estimators presented results visually identical to the tapered-SCM so they are omitted. This result can be explained because the oracle shrinkage coefficient from ST tends, in practice, to push the estimate towards the tapered matrix only: for example, the oracle coefficient (1-β) in [START_REF] Chen | Shrinkage-to-tapering estimation of large covariance matrices[END_REF] is always greater than 0.9 and is most likely close to 1.

[T θ ] ij = sinc((ij)k/π), which performs a sliding window average on the estimated signal spectrum. The one considered here performs a linear combination of the original spectrum with such average. This modification was made so that the tapering matrix always conforms to the theoretical requirements w ii = 1 and w ij ≥ 0, but did not significantly impacted the output of the tested detectors.

First we can notice that the SCM provides an unreliable detection maps in both scenarios, which is due to insufficient sample support in this configuration (n < 2p). As observed in [START_REF] Guerci | Theory and application of covariance matrix tapers for robust adaptive beamforming[END_REF] on another dataset, the covariance matrix tapering can widen the clutter notch, i.e., properly cancel the response of the detector on the anti-diagonal of the detection map. This permits to clearly distinguish several targets when compared to the detection map of the SCM. However, this improvement is at the cost of canceling the response of slower targets in the second scenario. This was to be expected since slow targets are hard to distinguish from the ground response in practice. Within the considered framework, hard cancellation of the ground clutter generally implies to also cancel these targets, which implies a trade-off when selecting the nullspectrum width of the tapering matrix. The shrinkage to identity of RSCM also greatly improves the detection process, as it allows us to detect all the 10 targets in the second scenario. However, it still presents some high false alarms on the clutter ridge. Finally, TABASCO appears as an interesting trade-off by combining the two effects, and illustrates that the proposed NMSE-driven method still allows for a reasonable regularization parameters (both β and k) selection in this detection application. More precisely, TABASCO yields the best detection map in the first scenario (no null-spectrum pattern and no false alarms around the clutter ridge). In the second scenario, slower targets are still canceled but TABASCO allows for recovering some targets that were initially canceled by a "tapering only" approach.

VIII. CONCLUSIONS AND PERSPECTIVES

We proposed TABASCO: a new covariance matrix estimator that jointly benefits from shrinkage to a scaled identity matrix and tapering of the SCM. By assuming the samples to be generated from an unspecified ES distribution, we also derived an efficient and robust estimation method for the oracle regularization parameters that minimize the MSE. Simulations studies illustrated that TABASCO outperforms existing regularized and tapered estimators in numerous setups. Interestingly, if W = 11 ⊤ belongs to the set of tapering matrices W considered, the estimator can avoid applying tapering if this option does not provide reduction to the MSE. Thus TABASCO performs similarly to the regularized SCM proposed in [START_REF] Ollila | Optimal shrinkage covariance matrix estimation under random sampling from elliptical distributions[END_REF] in this case, while significantly outperforming it when the tapering templates are valid. We also proposed two new novel estimators that measure the sphericity of the tapered covariance matrix.

APPENDIX

A. Known location µ

In some applications, the mean vector µ = E[x] is known and assumed to be µ = 0 without loss of generality. This is the case in STAP application for example. In this case, the covariance matrix Σ = E[xx ⊤ ] is estimated by the SCM, defined by

S = 1 n n i=1 x i x ⊤ i . (45)
This result follows as a special case of [31, Lemma 1] for a Gaussian weight function. Similarly Lemma 1 holds when using τ 1 and τ 2 in (46). The change to the optimal β 0 parameter is also minimal: one may ignore the term (n/(n -1)) that appears as the multiplier of the 2nd last term pθ W /η 2 +γ W -2γ/p in the denominator of β 0 in Theorem 2.

Theorem 4 also holds with

a n = 1 + κ n + κ and b n = n(n + κ) (n -1)(n + 2 + 3κ) . B. Proof of Theorem 1 Write L(β) = MSE( Σβ ) = E[∥ Σβ -Σ∥ 2 F ]. Then note that L(β) = E β(W • S) + (1 -β)p -1 tr(S)I -Σ 2 F = E β(W • S -Σ) + (1 -β) p -1 tr(S)I -Σ) 2 F = β 2 a 1 + (1 -β) 2 a 2 + 2β(1 -β)a 3 ( 47
)
where the constants a i -s are defined by

a 1 = MSE(W • S), a 2 = E p -1 tr(S)I -Σ 2 F , a 3 = E tr((W • S -Σ)(p -1 tr(S)I -Σ)) .
Then define 2 ) while a 2 can be expressed as

ã3 = p -1 E [tr (W • S -Σ) tr(S)] = p -1 E tr(S)
a 2 = ã3 + ∥Σ∥ 2 F -pη 2 . ( 50
)
Note that L(β) is a convex quadratic function in β with a unique minimum given by

β o = a 2 -a 3 (a 1 -a 3 ) + (a 2 -a 3 ) . (51)
Using ( 49) and (50), the numerator is

a 2 -a 3 = ∥V • Σ∥ 2 F -pη 2 (52) = ∥V • Σ -ηI∥ 2 F = p(γ V -1)η 2 ,
where we used that tr(V • Σ) = tr(Σ). Using the expression for a 1 = MSE(W • S) given in (13) together with equations (49) and (48), the first term in the denominator of

β o is a 1 -a 3 = E ∥W • S∥ 2 F -∥V • Σ∥ 2 F -p -1 E tr(S) 2 + η 2 p.
Summing this with term a 2 -a 3 from (52) shows that the denominator of

β o is E[∥W • S∥ 2 F ] -p -1 E[tr(S) 2
]. These results thus yield expression given in ( 15) and ( 16) for β o . The final expression [START_REF] Srivastava | Some tests concerning the covariance matrix in high dimensional data[END_REF] for β o can be deduced from ( 16) by using [START_REF] Rugini | Regularized MMSE multiuser detection using covariance matrix tapering[END_REF] and then simplifying the expression.

The expression for MSE of Σβ o follows by substituting β o into expression for L(β) in (47) and using the relation,

(1 -β o )(a 2 -a 3 ) = β o (a 1 -a 3 ), which follows from (51). This gives L(β o ) = a 2 -β o (a 2 -a 3 ) = a 3 + (1 -β o )(a 2 -a 3 ).
This gives the stated MSE expression since a 2 -a 3 = ∥V • Σ -ηI∥ 2 F and using (49) together with (48) for a 3 .

C. Proof of Lemma 1

Before proceeding with the proof we introduce some definitions and results that are used in the sequel. First, we let K p denote the p 2 × p 2 commutation matrix defined as a block matrix whose ijth block is equal to a p × p matrix that has a 1 at element ji and zeros elsewhere, i.e., K p = ∑ i,j e i e ⊤ j ⊗ e j e ⊤ i . It also has the following important properties [START_REF] Magnus | Matrix Differential Calculus with Applications in Statistics and Econometrics[END_REF]:

K p vec(A) = vec(A ⊤ ) and K p (A ⊗ B)K p = (B ⊗ A)
for any p × p matrices A and B, where vec(A) vectorizes matrix A by stacking the columns of the matrix on top of each other. We then have the following identities. Lemma 3. The following holds:

a) ∥A • B∥ 2 F = tr vec(A)vec(A) ⊤ • vec(B)vec(B) ⊤ for all A, B ∈ R m×n . b) ∥A•B∥ 2 F = tr vec(A)vec(A) ⊤ • K p (B ⊗ B) ∀A ∈ R m×m and ∀B ∈ R m×m Sym . c) d ⊤ B (A • A)d B = tr vec(A)vec(A) ⊤ • (B ⊗ B) ∀A, B ∈ R m×m . d) tr (D B A) 2 = d ⊤ B (A • A)d B for all A ∈ R
where τ 1 and τ 2 are constants defined in [START_REF] Magyar | The asymptotic inadmissibility of the spatial sign covariance matrix for elliptically symmetric distributions[END_REF]. Equations ( 54 which proves the latter claim.

D. Proof of Theorem 3

Let us express the SSCM as

Λ = 1 n n i=1 v i v ⊤ i
, where v i = p

x i ∥x i ∥ .

Hence ∥W • Λ∥ 2 F p = 1 pn 2 tr W • v 1 v ⊤ 1 + ⋯ + W • v n v ⊤ n 2 = n i=1 ∥W • v i v ⊤ i ∥ 2 F pn 2 + n i≠j tr W • v i v ⊤ i W • v j v ⊤ j pn 2 .
Then since v i -s are i.i.d., and

E[ Λ] = E[v i v ⊤ i
] for all i, the expectation of the 2nd term is

i≠j E tr (W • v i v ⊤ i )(W • v j v ⊤ j ) pn 2 = n -1 n ∥W • Λ sgn ∥ 2 F p ,
where Λ sgn = E[ Λ]. The expectation of the 1st terms is Using ( 59) and ( 60) we then obtain that (62)

i E ∥W • v i v ⊤ i ∥ 2 F pn 2 = E ∥W • vv ⊤ ∥ 2 F pn = E d ⊤ (W • W)d
E[γ W ] = ∥W • Λ sgn ∥ 2 F p + 1 n ε, (61)
Next note that Λ sgn = E[ Λ] = Λ + o(∥Λ∥ F ) when (A) holds by [START_REF] Raninen | Linear pooling of sample covariance matrices[END_REF]Theorem 2]. This fact together with (61) and (62) imply that

E[γ W ] → ∥W • Λ∥ 2 F p = γ W
as p → ∞ under assumption (A). Thus we have proven the claim.

E. Proof of Lemma 2: complex case

In our proof we will use the following identities.

where τ 1 and τ 2 are constants defined in [START_REF] Magyar | The asymptotic inadmissibility of the spatial sign covariance matrix for elliptically symmetric distributions[END_REF]. Equations ( 64) and (65) then imply that

E[vec(S)vec(S)

H ] = cov(vec(S)) + vec(Σ)vec(Σ)

H = τ 1 (Σ * ⊗ Σ) + (1 + τ 2 )vec(Σ)vec(Σ) H . (66)
Inserting (66) into (63) yields

E ∥W • S∥ 2 F = (1 + τ 2 )∥W • Σ∥ 2 F + τ 1 tr (D Σ W) 2 .
simply by invoking identities in Lemma 4. This proves the first identity. The proof of latter part E tr((D S W)

2 ) is as earlier in the real-valued case in Appendix C.

2 F

 2 ] in Lemma 1. As a remedy, the following theorem derives a proper unbiased estimator of ϑ W which extends [4,Theorem 4] (provided that the elliptical kurtosis parameter κ is known).

Theorem 4 . 2 F

 42 Let {x i } n i=1 be an i.i.d. random sample from a pvariate elliptical distribution E p (µ, Σ, g) with finite 4th order moments. Then, an unbiased estimator of ϑ W = ∥W•Σ∥ /p for any finite n and p and any

  2 = x H x we denote the usual Euclidean norm in complex vector spaces, while ∥B∥ F = tr(B H B) denotes the Frobenius norm of a matrix B ∈ C m×n , where (⋅) H = [(⋅) * ] ⊤ denotes the conjugate transpose (or Hermitian transpose). For any x ∈ C, the notation | ⋅ | refers to modulus, so |x| 2 = xx *

3 Fig. 1 .

 31 Fig. 1. (Σ) ij as a function of |i -j|. Left: Model 1 in (41) with various correlation parameters . Right: Model 2 in (41) with various decay parameters α and ρ = 0.6. The dimension is p = 100.

Fig. 2 .

 2 Fig. 2. Solid line display NMSE of TABASCO estimator using fixed W(k) as in (3) when samples are drawn from a MVN distribution, Σ as in Model 1 in (41) with = 0.2 (top left), = 0.4 (top right), = 0.6 (bottom left), = 0.8 (bottom right); n = 50, p = 100. The horizontal dashed line correspond to the empirical NMSE obtained by TABASCO of Algorithm 1.

Fig. 3 .

 3 Fig. 3. Oracle NMSE of TABASCO estimator with fixed bandwidth k (◯) when Σ has an AR(1) structure and ∈ {0.2, 0.4, 0.6, 0.8}. The average NMSE of estimated TABASCO using fixed k and estimated k over 5000 MC trials is superimposed to the curves using symbols ▲ and ☆, respectively. Samples are drawn from a MVN distribution, n = 50 and p = 100. Banding matrices W(k) as in (3) are used.

Fig. 4 . 8 Fig. 5 .

 485 Fig. 4. Average NMSE curves when samples are from a MVN distribution (upper panel) and t-distribution with ν = 5 d.o.f. (lower panel), Σ has an AR(1) structure with ∈ {0.2, 0.4, 0.6, 0.8} from left to right. Dimension is p = 100 and banding matrices are used in STOA, ST-gaus, ST-nong and TABASCO.

Fig. 6 .

 6 Fig. 6. Average NMSE curves when samples are from a MVN distribution (top row) and MVT distribution (bottom row) with ν = 5 d.o.f., Σ has a permuted AR(1) structure with = 0.2 (left panel) and = 0.4 (right panel), and dimension is p = 100.

Figure 7 Fig. 7 .

 77 Figure 7 presents a comparison similar to Figure 4 where we also included the minimax risk tapering (MnMx-Taper) estimator W(k * )•S, where k * = ⌊n 1/(2(α+1)) ⌋ is the optimal

Fig. 8 .

 8 Fig. 8. The true (theoretical) NMSE curves as a function of bandwidth k for the tapered SCM W • S and TABASCO Σβ when sampling from a MVN distribution (left panel) and MVT distribution (right panel) with ν = 5 d.o.f., Σ follows model 2 with α = 0.1, n = 100 and p = 250.

Fig. 9 . 2 F

 92 Fig. 9. Output of various STAP detectors for the first (top) and second (bottom) scenarios

2 F -∥V • Σ∥ 2 F

 22 that E[S] = Σ, tr(W • S) = tr(S) and E[tr(S)] = tr(E[S]) = ηp with η = tr(Σ)/p denoting the scale. We may write a 3 in the form a 3 = ã3 + E tr((W • S -Σ)Σ) = ã3 + ∥Σ∥ (49) by using E[S] = Σ and E[tr((W•S)Σ)] = tr((W•Σ)Σ) = tr((V • Σ)

B∥ 2 F 2 F

 22 Let A = (a ij ) and B = (b ij ). a) First note that ∥A • = tr vec(A • B)vec(A • B) ⊤ = tr vec(A)vec(A) ⊤ • vec(B)vec(B) ⊤ . b) It is a simple matter to verify that for all B ∈ R m×m Sym it holds that diag(vec(B)vec(B) ⊤ ) = diag(K p (B ⊗ B)). Thus tr vec(A)vec(A) ⊤ • K p (B ⊗ B) = tr vec(A)vec(A) ⊤ • vec(B)vec(B) ⊤ which gives the stated result due to a)-part. c) It is a simple task to verify that the trace of the Hadamard product of vec(A)vec(A) ⊤ with B ⊗ B equals ∑ i,j b ii a 2 ij b jj which is equivalent with d ⊤ B (A • A)d B . d) Follows from [21, Lemma 7.5.2].Write w = vec(W). Using Lemma 3a) we first notice thatE ∥W • S∥ = tr ww ⊤ • E[vec(S)vec(S) ⊤ ] .(53)We then recall that the (variance-)covariance matrix of S when sampling from an elliptical population E p (µ, Σ, g) is given by [4, Theorem 2]: cov(vec(S)) = E[vec(S)vec(S) ⊤ ] -vec(Σ)vec(Σ) ⊤ (54) = τ 1 (I + K p )(Σ ⊗ Σ) + τ 2 vec(Σ)vec(Σ) ⊤ ,

2 F = ( 1 + τ 1 + τ 2 )∥W • Σ∥ 2 FE[s ii s jj ] = 2τ 1 σ 2 ij + ( 1 + τ 2 2 ) = 2τ 1 p= 2τ 1 ∥W • Σ∥ 2 F

 21122212212 ) and (55) then imply thatE[vec(S)vec(S)⊤ ] = cov(vec(S)) + vec(Σ)vec(Σ)⊤ = τ 1 (I + K p )(Σ ⊗ Σ) + (1 + τ 2 )vec(Σ)vec(Σ) ⊤ . (56)Inserting (56) into (53) yieldsE ∥W • S∥ + τ 1 tr (D Σ W) 2 .simply by invoking identities in Lemma 3. This proves the first identity. Next we note that E tr((D S W) inserting (58) into (57) yields E tr((D S W) + (1 + τ 2 ) tr((D Σ W) 2 )

pn where d = (v 2 1 2 = 2 =

 222 , . . . , v 2 p )⊤ contains the diagonal elements of vv ⊤ , where v = d v i and = d reads "has the same distribution as". Furthermore, write D = diag(vv ⊤ ). Thus we have thatn n -1 ⋅ E ∥W • Λ∥ 2 F p = ∥W • Λ sgn ∥ 2 F p + E d ⊤ (W • W)d p(n -1) . (59)Next note that D Λ = diag( Λ) can be written asD Λ = 1 n D 1 + . . . + D n ,whereD i = diag(v i v ⊤ i ). Furthermore, let d i = (v 2 i1 , . . . , v 2 ip )⊤ denote a random vector containing the diagonal elements of v i v ⊤ i . Then we gettr DΛW 2 = 1 n 2 tr D 1 W + . . . + D n W E[d(W • W)d] pn(n -1) + 1 pn tr E[D]W E[d(W • W)d] pn(n -1) + 1 pn E[d] ⊤ (W • W)E[d].(60)

where ε = 1 pw

 1 E[d(W • W)d] -E[d] ⊤ (W • W)E[d] ij cov(v i v j ) → 0 as p → ∞.

Lemma 4 . 2 F 2 F

 422 The following holds: a) ∥A•B∥ = tr vec(A)vec(A)H • vec(B)vec(B) H for all A, B ∈ C m×n . b) d ⊤ B (A • A)d B = tr vec(A)vec(A) ⊤ • (B * ⊗ B) ∀A ∈ C m×m and B ∈ C m×m Sym . c) tr (D B A) 2 = d ⊤ B (A • A)d B for all A ∈ Rm×m Sym and B ∈ C m×m . Proof. a,b) proofs of the identities are as proofs of Lemma 3a),b). c) follows directly from [21, Lemma 7.5.2].Write w = vec(W). Using Lemma 4a) we first notice thatE ∥W • S∥ = tr ww ⊤ • E[vec(S)vec(S) H ] . (63)We then recall that the (variance-)covariance matrix of S when sampling from a complex elliptically symmetric distribution CE p (µ, Σ, g) is [39, Theorem 3]:cov(vec(S)) = E[vec(S)vec(S) H ] -vec(Σ)vec(Σ) H (64) = τ 1 (Σ * ⊗ Σ) + τ 2 vec(Σ)vec(Σ) H ,

  contains the variances of the variables, i.e., the diagonal elements of Σ. Obviously, if W = 11

2 1 , . . . , σ 2 p ) ⊤ ⊤ , then θ W = η 2 . The 2nd equality in (27) follows from [21, Lemma 7.5.2]. Theorem 2. Let {x i } n i=1 be an i.i.d. random sample from an ES distribution E p (µ, Σ, g) with finite 4th order moments. For any W ∈ W + , the oracle parameter β o in (12) is

We denote diag(A) ≡ diag(a 11 , . . . , a pp ) when the operator is applied to any matrix A = (a ij ) p×p . Conversely diag(a) denotes a diagonal matrix with the entries of vector a on the main diagonal.

Prominent algorithms for recovering hidden ordering-structure in the variables are the Best Permutation Analysis (BPA)[START_REF] Rajaratnam | Best permutation analysis[END_REF] or Isoband[START_REF] Wagaman | Discovering sparse covariance structures with the isomap[END_REF]. The perspective of their joint use with TABASCO is left for further studies.

The tapering in[START_REF] Guerci | Theory and application of covariance matrix tapers for robust adaptive beamforming[END_REF] actually uses T f ij = sinc((ij)k/π) and