{'fraction_non_alphanumeric': 0.03759165054572377, 'fraction_numerical': 0.011487664493324735, 'mean_word_length': 4.709299154622307, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Linear Shrinkage Estimation of Covariance Matrices Using Low-Complexity Cross-Validation

Shrinkage can effectively improve the condition number and accuracy of covariance matrix estimation, especially for low-sample-support applications with the number of training samples smaller than the dimensionality. This paper investigates parameter choice for linear shrinkage estimators. We propose data-driven, leave-one-out cross-validation (LOOCV) methods for automatically choosing the shrinkage coefficients, aiming to minimize the Frobenius norm of the estimation error. A quadratic loss is used as the prediction error for LOOCV. The resulting solutions can be found analytically or by solving optimization problems of small sizes and thus have low complexities. Our proposed methods are compared with various existing techniques. We show that the LOOCV method achieves near-oracle performance for shrinkage designs using sample covariance matrix (SCM) and several typical shrinkage targets. Furthermore, the LOOCV method provides low-complexity solutions for estimators that use general shrinkage targets, multiple targets, and/or ordinary least squares (OLS)-based covariance matrix estimation. We also show applications of our proposed techniques to several different problems in array signal processing.


I. INTRODUCTION
In statistical signal processing, one critical problem is to estimate the covariance matrix, which has extensive applications in correlation analysis, portfolio optimization, and various signal processing tasks in radar and communication systems [1]- [5]. One key challenge is that when the dimensionality is large but the sample support is relatively low, the estimated covariance matrix R, which may be obtained using a general method such as sample covariance matrix (SCM) or ordinary least squares (OLS), becomes ill-conditioned or even singular, and suffers from significant errors relative to the true covariance matrix Σ. Consequently, signal processing tasks that rely on covariance matrix estimation may perform poorly or fail to apply. Regularization techniques have attracted tremendous attention recently for covariance matrix estimation. By imposing structural assumptions of the true covariance matrix Σ, techniques such as banding [6], thresholding [7], and shrinkage [8]- [18] have demonstrated great potential for improving the performance of covariance matrix estimation. See [19]- [21] for recent surveys. This paper is concerned with the linear shrinkage estimation of covariance matrices. Given an estimate R of the covariance matrix, a linear shrinkage estimate is constructed as where T 0 is the shrinkage target and ρ and τ are nonnegative shrinkage coefficients. In general, the shrinkage target T 0 is better-conditioned, more parsimonious or more structured, with lower variance but higher bias compared to the original estimate R [11]. The coefficients ρ and τ are chosen to provide a good tradeoff between bias and variance, such that an estimate outperforming both R and T 0 is achieved and a better approximation to the true covariance matrix Σ can be obtained. Compared to other regularized estimators such as banding and thresholding, linear shrinkage estimators can be easily designed to guarantee positive-definiteness. Such shrinkage designs have been employed in various applications which utilize covariance matrices and have demonstrated significant performance improvements. The linear shrinkage approach has also been generalized to nonlinear shrinkage estimation of covariance matrices [22], [23], and is closely related to several unitarily invariant covariance matrix estimators that shrink the eigenvalues of the SCM, such as those imposing condition number constraints on the estimate [24], [25]. There are also a body of studies on shrinkage estimation of precision matrix (the inverse of covariance matrix) [26]- [30] and on application-oriented design of shrinkage estimators. See [31]- [36] for example applications in array signal processing.
Shrinkage has a Bayes interpretation [2], [9]. The true covariance matrix Σ can be assumed to be within the neighborhoods of the shrinkage target T 0 . There can be various different approaches for constructing R and T 0 . For example, when a generative model about the observation exists, one may first estimate the model parameters and then construct R [20]. A typical example of this is linear models seen in communication systems. Furthermore, different types of shrinkage targets, not necessarily limited to identity or diagonal targets, can be used to better utilize prior knowledge. For example, knowledgeaided space-time signal processing (KA-STAP) may set T 0 using knowledge about the environment [3] or past covariance matrix estimates [37]. Even multiple shrinkage targets can be applied when distinct guesses about the true covariance matrix are available [17].
The choice of shrinkage coefficients significantly influences the performance of linear shrinkage estimators. Various criteria and methods have been studied. Under the mean squared error (MSE) criterion, Ledoit and Wolf (LW) [2] derived closed-form solutions based on asymptotic estimates of the statistics needed for finding the optimal shrinkage coefficients, where R and T 0 are assumed as the SCM and identity matrix, respectively. Later the LW solution was extended for more general shrinkage targets [3], [17]. Chen et al [4] assumed Gaussian distribution and proposed an oracle approximating shrinkage (OAS) estimator, which achieves near-optimal parameter choice for Gaussian data even with very low sample supports. The shrinkage coefficients determination can also be cast as a model selection problem and thus generic model selection techniques such as cross-validation (CV) [38]- [40] can be applied. In general, CV splits the training samples for multiple times into disjoint subsets and then fits and assesses the models under different splits based on a properly chosen prediction loss. This has been explored, e.g., in [10], [13], where the Gaussian likelihood is used as the prediction loss.
All these data-driven techniques achieve near-optimal parameter choice when the underlying assumptions hold. However, there are also limitations to their applications: almost all existing analytical solutions to shrinkage coefficients [2]- [4], [17] were derived under the assumption of SCM and certain special forms of shrinkage targets. They need to be re-designed when applied to other cases, which is generally nontrivial. The asymptotic analysis-based methods [2], [3] may not perform well when the sample support is very low. Although the existing CV approaches [10], [13] have broader applications, they assume Gaussian distribution and employ grid search to determine the shrinkage coefficients. The likelihood cost of [10], [13] must be computed for multiple data splits and multiple candidates of shrinkage coefficients, which can be time-consuming.
In this paper, we further investigate data-driven techniques that automatically tune the linear shrinkage coefficients using leave-one-out cross-validation (LOOCV). We choose a simple quadratic loss as the prediction loss for LOOCV, and derive analytical and computationally efficient solutions. The solutions do not need to specify the distribution of the data. Furthermore, the LOOCV treatment is applicable to different covariance matrix estimators including the SCM-and ordinary least squares (OLS)-based schemes. It can be used together with general shrinkage targets and can also be easily extended to incorporate multiple shrinkage targets. The numerical examples show that the proposed method can achieve oracleapproximating performance for covariance matrix estimation and can improve the performance of several array signal processing schemes.
The remainder of the paper is organized as follows. In Section 2, we present computationally efficient LOOCV methods for choosing the linear shrinkage coefficients for both SCM-and OLS-based covariance matrix estimators and also compare the proposed LOOCV methods with several existing methods which have attracted considerable attentions recently. In Section 3, we extend our results for multi-target shrinkage. Section 4 reports numerical examples, and finally Section 5 gives conclusions.
II. LOOCV CHOICE OF LINEAR SHRINKAGE COEFFICIENTS This paper deals with the estimation of covariance matrices of zero-mean signals whose fourth-order moments exist. We study the LOOCV choice of the shrinkage coefficients for the linear shrinkage covariance matrix estimator (1), i.e., Σ ρ,τ = ρR + τ T 0 . The following assumptions are made: 1) The true covariance matrix Σ, the estimated covariance matrix R, and the shrinkage target T 0 are all Hermitian and positive-semidefinite (PSD). 2) T independent, identically distributed (i.i.d.) samples {y t } of the signal are available.
3) The shrinkage coefficients are nonnegative, i.e., Assumption 3 follows the treatments in [2]- [4] and is sufficient but not necessary to guarantee that the shrinkage estimate Σ ρ,τ is PSD when Assumption 1 holds 1 . Two classes of shrinkage targets will be considered in this paper. One is constructed independent of the training samples {y t } for generating R, similarly to the knowledge-aided targets as considered in [3]. The other is constructed from {y t }, but is highly structured with significantly fewer free parameters as compared to R. Examples of the second class include those constructed using only the diagonal entries of R [4], [20] and the Toeplitz approximations of R [17].

A. Oracle Choice
Different criteria may be used for evaluating the covariance matrix estimators. In this paper, we use the squared Frobenius norm of the estimation error as the performance measure. Given Σ, R and T 0 , the oracle shrinkage coefficients minimize where || · || F denotes the Frobenius norm. The cost function in (3) can then be rewritten as a quadratic function of the shrinkage coefficients: where tr(·) denotes the trace of a matrix. As A O is positivedefinite, we can find the minimizer of J O (ρ, τ ) by solving the above bivariate convex optimization problem. We can also apply the Karush-Kuhn-Tucker (KKT) conditions to find the solution analytically. From (4), letting JO(ρ,τ ) The oracle shrinkage coefficients can be obtained by solving (7) and (8): Note that (9) may produce negative shrinkage coefficients, which may not lead to a positive-definite estimate of the covariance matrix. In this case, we clip the negative coefficient to zero and then find the other coefficient using (7) or (8) to guarantee the positive definiteness, for τ = 0 or ρ = 0, respectively. This treatment is similar to [2]- [5] and provides a suboptimal yet simple solution. The oracle estimator requires knowledge of Σ, which is unavailable in real applications, but the result serves as an upper bound of the performance given the linear shrinkage structure.

B. LOOCV Choice for General Cases
Let Σ denote a positive-definite, Hermitian matrix. It can be easily verified that the following cost is minimized when Σ = Σ, where the expectation is taken over y. In this paper, we apply LOOCV [38] to produce an estimate of J S ( Σ) as the proxy for measuring the accuracy of Σ, based on which the shrinkage coefficients can be selected.
With the LOOCV method, the length-T training data Y = [y 1 , y 2 , · · · , y T ] is repeatedly split into two sets with respect to time. For the t-th split, where 1 ≤ t ≤ T , T − 1 samples in Y t (with the t-th column y t omitted from Y) are used for producing a covariance matrix estimate R t and the remaining sample y t is spared for parameter validation. In total, T splits of the training data Y are used and all the training samples are used for validation once. Assuming shrinkage estimation with given shrinkage coefficients (ρ, τ ), we construct from each Y t a shrinkage covariance matrix estimator as We propose to use the following LOOCV cost function to approximate the cost in (10) when Σ is chosen as Σ t,ρ,τ . For notational simplicity, define After some manipulations, the above cost function can be written similarly to (4) as where The shrinkage coefficients can then be found by solving the above bivariate, constant-coefficient quadratic program. Analytical solutions can be obtained under different conditions, as shown below. 1) Unconstrained shrinkage: For unconstrained (ρ, τ ), setting the partial derivatives ∂JCV(ρ,τ ) Solving (18) and (19) produces the unconstrained solution We choose (20) as the optimal shrinkage coefficients if both ρ ⋆ CV and τ ⋆ CV are nonnegative. Otherwise, we consider the optimal choices on the boundary of ρ ≥ 0, τ ≥ 0 specified by (18) or (19) for τ = 0 or ρ = 0 as or 2) Constrained shrinkage: For the more parsimonious design using convex linear combination, the following constraint is imposed: By plugging (23) into the cost function (12) and taking the minimizer, we can also easily find the optimal shrinkage coefficients using 4 In case a negative shrinkage coefficient is produced, we set it to zero and let the other be one according to (23). Note that although the closed-form solution involves multiple matrix operations, the quantities involved need to be computed only once. Furthermore, the computational complexity may be greatly reduced given a specific method of covariance matrix estimation. In the following two subsections, we will show the simplified solutions for SCM-and OLS-based covariance matrix estimation.

C. LOOCV Choice for SCM-Based Estimation
We consider in this subsection that R is the SCM estimate of Σ. In this case, which is a sufficient statistic for Gaussian-distributed data when the mean vector is the zero vector. For the t-th split, the SCM constructed from all the samples except the t-th is We can then verify the following expressions for quickly computing the relevant quantities in (16) and (17): Plugging these into (16) and (17) and after some manipulations, we can rewrite the LOOCV cost function (15) as The optimal shrinkage coefficients can then be obtained analytically from the SCM R, the shrinkage target T 0 , and the training samples {y t }, as discussed below. 1) Unconstrained shrinkage: It can be verified from (19) and (30) that the optimal shrinkage coefficients (ignoring the nonnegative constraint ρ ≥ 0, τ ≥ 0) satisfy The closed-form solution to ρ is given by In case ρ ⋆ CV,SCM > 1 or ρ ⋆ CV,SCM < 0, we apply (21) or (22), respectively, to determine the solution, using the expressions in (27)- (30).
Note that for the typical shrinkage target T 0 = tr(R) This provides another justification for the convex linear combination design with an identity target, which has been widely adopted in the literature, e.g., [4]. This also shows that for such a special target the unconstrained solution is equivalent to the constrained solution, which does not hold for more general shrinkage targets.
2) Constrained shrinkage: For the widely considered convex linear combination with constraint ρ + τ = 1, the optimal ρ (ignoring the nonnegative constraint) is computed as (34) Similarly, in case a negative shrinkage coefficient is obtained, we set it to zero and let the other be one.
The above results show that the optimal shrinkage coefficients for the covariance matrix estimate (1) can be computed directly from the samples and shrinkage target, without the need of specifying any user parameters. The constrained shrinkage design may lead to certain performance loss as compared to the unconstrained one.

D. LOOCV Choice for OLS-Based Covariance Estimation
One advantage of the LOOCV method is that it can be applied to different covariance matrix estimators. In this subsection, we discuss the LOOCV method for OLS-based covariance matrix estimation. Note that most existing analytical solutions for choosing the shrinkage coefficients assume SCM and specific shrinkage targets and need to be re-derived for general cases. Also, in contrast to general applications of LOOCV which require a grid search of the parameters and thus a high computational complexity, we have shown that for SCM, fast analytical solutions can be obtained for choosing the shrinkage coefficients. This will also be the case for the OLS-based covariance matrix estimation.
Consider the case with observation y ∈ C N modeled as where H ∈ C N ×M is a deterministic channel matrix and z ∈ C N a zero-mean, white noise with covariance matrix σ 2 I, which is uncorrelated with the zero-mean input signal x ∈ C M with covariance matrix I. If both training samples of x and y are known, we may first estimate the channel matrix H and the covariance matrix of the noise z using the ordinary least squares (OLS) approach. Let the block of training data be (X, Y), where the input signal X can be designed to have certain properties such as being orthogonal. The OLS estimates of the channel matrix and noise variance are then obtained as where (·) denotes the estimate of a quantity. In this case, the covariance matrix of y can be estimated as Such OLS-based covariance matrix estimation may be useful for designing signal estimation schemes in wireless communications. We can apply the linear shrinkage design (1) to enhance its accuracy and apply the LOOCV method (12) to choose the shrinkage coefficients. Note that in this case, in the t-th split, we generate the covariance matrix estimate R t by applying the OLS estimate to the leave-one-out samples The LOOCV cost is the same as (15). In this case, the leave-one-out estimate of the covariance matrix for the t-th data split is where H t and σ 2 t denote the channel matrix and noise variance estimated from (X t , Y t ), respectively. A direct computation of (16) and (17) for evaluating the LOOCV cost performs OLS estimation for T times, which incurs significant complexity. The complexity can be greatly reduced by observing that the leave-one-out OLS estimate of the channel matrix is related to the OLS channel matrix estimate H in (36) by a rank-one update: where In the above, is the t-th diagonal entry of Similarly, the leave-one-out estimate of the noise variance can be updated as where Note that both updates can be achieved with low complexity when a few matrices are computed in advance and reused. In this way, the covariance matrix estimate can be computed as where φ t and ψ t are defined as This shows that the leave-one-out OLS covariance matrix estimate can be obtained from R by corrections involving a scaled identity matrix and two rank-one updates. Eqn. (47) can be exploited to compute the closed-form LOOCV solution quickly. From (47), the most involved computation for finding the solution of the optimization problem (15) can be implemented as where R(·) denotes the real part of a scalar. When R is already computed, the right-hand side of (50) can be evaluated using inner products and matrix-vector products. The terms tr(T 2 0 ) and tr(T 0 S t ) are the same as those for SCM. For the other two terms, we have Note that the computational complexities of (50)-(52) are low because the major operations are matrix-vector products and inner products.

E. Comparisons with Alternative Choices of Linear Shrinkage Coefficients
In the above, we have introduced LOOCV methods with analytical solutions for choosing the coefficients for linear shrinkage covariance matrix estimators. We now discuss several alternative techniques which have received considerable attentions recently and compare them with the LOOCV methods proposed in this paper.
In 2004, Ledoit and Wolf (LW) [2] studied estimators that shrink SCM toward an identity target, i.e., T 0 = I. Such estimators do not alter the eigenvectors but shrink eigenvalues of the SCM, which is well supported by the fact that sample eigenvalues tend to be more spread than population eigenvalues. The optimal shrinkage coefficients under the MMSE criterion (3) can be written as where the parameters µ tr(Σ) F depend on the true covariance matrix Σ and other unknown statistics. Ref. [2] shows that δ 2 = α 2 +β 2 and proposes to approximate these quantities by their asymptotic estimates under T → ∞, which can all be computed from the training samples. By substituting these into (53), estimators that significantly outperform SCM are obtained, which also approach the oracle estimators when the training length is large enough.
The above LW estimator is extended by Stoica et al in 2008 [3] for complex-valued signals with general shrinkage targets T 0 , with applications to knowledge-aided space-time adaptive processing (KA-STAP) in radar applications. Several estimators with similar performance are derived there. For the general linear combination (GLC) design of [3], it is shown that the oracle shrinkage coefficients for (1) satisfy where The quantity β 2 is estimated in the same way as (55), and a computationally efficient expression for β 2 is given by Furthermore, ν and E[||R − νT 0 || 2 F ] are estimated as ν = tr(T0R) ||T0|| 2 F and ||R− νT 0 || 2 F , respectively. This leads to the result given by Eqns. (34) and (35) of [3], which can recover the LW estimator [2] when the identity shrinkage target T 0 = I is assumed.
More recently, Chen et al [4] derived the oracle approximating shrinkage (OAS) estimator, which assumes SCM, real-valued Gaussian samples, and scaled identity target with T 0 = tr(R) N I and ρ = 1 − τ . They first derive the oracle shrinkage coefficients for SCM obtained from i.i.d. Gaussian samples, which is determined by N, T, tr(Σ) and tr(Σ 2 ). Then, they propose an iterative procedure to approach the oracle estimator. In the iterations, tr(Σ 2 ) and tr(Σ) are estimated by tr( Σ j R) and tr( Σ j ), respectively, where Σ j is the covariance matrix estimate at the j-th iteration. It is further proved that Σ j converges to the OAS estimator with the following analytical expression for τ : This approach achieves superior performance for (scaled) identity target and Gaussian data and dominates the LW estimator [2] when T is small. It was later generalized by Senneret et al [20] to a shrinkage target chosen as the diagonal entries of the SCM. Other related techniques include [14], which also assumes SCM, Gaussian data, and identity/diagonal shrinkage targets.
All the above techniques provide analytical solutions and achieve near-oracle performance when the underlying assumptions (e.g., large dimensionality, large size of training data, identity/diagonal shrinkage targets) hold. However, they also have limitations. A common restriction is that all these analytical solutions assume SCM and are not optimized for other types of covariance matrix estimators such as modelbased estimators. In particular, the LW and GLC methods [2], [3], which employ asymptotic approximations, may exhibit a noticeable gap to the oracle choice when the sample support is low, which may be relevant in some applications. The OAS method [4] assumed identity target, but its extensions to more general cases, e.g., with multiple/general shrinkage targets, are not trivial. By contrast, the LOOCV method proposed in this paper allows different designs and achieves near-oracle performance in general.
Cross-validation has also been applied previously for choosing shrinkage coefficients for covariance matrix estimation.
The key issues for applying this generic tool include finding appropriate predictive metrics for scoring the different estimators and fast computation schemes. In [10], [13], the Gaussian likelihood was chosen as such a proxy. The computations with likelihood are generally involved as multiple matrix inverses/determinants are required, and a grid search is required for finding the optimal parameters. In this paper, we use the distribution-free, Frobenius norm loss in (12) as the metric, which leads to analytical solutions and is computationally more tractable.

III. MULTI-TARGET SHRINKAGE
In Section 2, we have considered linear shrinkage designs with a single target. Multiple shrinkage targets may be used to further enhance performance, which may be obtained from a priori knowledge, e.g., a past covariance matrix estimate from older training samples or from neighboring frequencies. We can easily extend our proposed LOOCV method to multiple targets.

A. Oracle choice of shrinkage coefficients
Consider the multi-target shrinkage design where all the shrinkage coefficients are nonnegative to guarantee PSD covariance matrix estimates, i.e., The oracle multi-target shrinkage minimizes the squared Frobenius norm of the estimation error which can be rewritten as The oracle shrinkage coefficients can then be obtained by solving the problem of minimizing the cost function J O,MT (ρ, τ ) of (63), which is a strictly convex quadratic program (SCQP) with K + 1 variables.

B. LOOCV choice of shrinkage coefficients
We now extend the LOOCV method in Section 2 to the multi-target shrinkage here. Following the same treatment as in Section II-B, in each split of the training data, R t and S t are constructed to generate and validate the covariance matrix estimate, respectively. The multiple shrinkage coefficients are chosen to minimize the LOOCV cost The above cost function can be rewritten in a form similar to (15) as . . .
The constant entries of A CV,MT and b CV,MT can be computed in the same way as for the single-target case. When K is small, which is typically the case, the solution that minimizes the LOOCV cost can be found quickly using standard optimization tools. Alternatively, we may find first the global optimizer that ignores the nonnegative constraint by and check if the nonnegative condition is satisfied. If a negative shrinkage coefficient is produced, we then consider the boundaries of ρ ≥ 0, τ k ≥ 0, k = 1, 2, · · · , K, which are equivalent to removing a certain number of shrinkage targets from the shrinkage design. The solution can be found in exactly the same way as (70) but with fewer targets.
Similarly to the single-target case, we may also consider a constrained case, where the shrinkage targets {T k } have the same trace as the estimated covariance matrix R, and Then the LOOCV cost function can be rewritten as where The optimal shrinkage coefficients can be found similarly as for the unconstrained case by minimizing where the entries of A ′ CV,MT and b ′ CV,MT are defined by These entries may also be evaluated quickly. For example, with SCM, can be computed using (27). The solution to τ can be found as if the nonnegative condition is satisfied. Otherwise, find the solution in a similar way as for the unconstrained case. Note that for multi-target shrinkage, Lancewicki and Aladjem [17] recently introduced another method for finding the shrinkage coefficients. They assume SCM and shrinkage targets which belong to a set that can be characterized by Eqn. (21) of [17]. Then, they follow the Ledoit-Wolf (LW) framework [2] to derive unbiased estimates of the unknown coefficients needed for minimizing the expectation of the cost in (62), based on which {ρ, τ } can be optimized. By contrast, our approach resorts to a LOOCV estimate of the cost in (10), which does not rely on the aforementioned assumptions in [17]. As will be shown later, the LOOCV method can achieve similar performance as [17] for the shrinkage targets considered there. However, it can be applied to general estimators other than SCM and shrinkage targets which are not covered by Eqn. (21) of [17], offering wider applicability.

IV. NUMERICAL EXAMPLES
In this section, we present numerical examples to demonstrate the effectiveness of the proposed shrinkage design and compare it with alternative methods. The quality of covariance matrix estimation is measured by the MSE normalized by the average of the squared Frobenius norm ||Σ|| 2 F , i.e., We show examples of covariance matrix estimation and its applications in array signal processing. We denote by N (µ, σ 2 ) a real-valued Gaussian distribution with mean µ and variance σ 2 . Example 1: Shrinkage toward an identity target: We first consider a real-valued example with an autoregressive (AR) covariance matrix, whose (i, j)-th entry is given by N I. "LW", "GLC" and "OAS" refer to the methods of [2], [3] and [4], respectively, which are also described in Section II-E; "CV" refers to our proposed LOOCV method; "Oracle" refers to the coefficient choice in Section II-A; and "Con" and "Unc" indicate that the constraint ρ + τ = 1 is imposed or not, respectively. which has been widely considered for evaluating covariance matrix estimation techniques [4]- [7]. Let Σ 1/2 be the Cholesky factor of Σ. The training samples are randomly generated as y t = Σ 1/2 n t , where n t consists of i.i.d. entries drawn from N (0, 1). The typical shrinkage target T 0 = tr(R) N I is considered for single-target shrinkage. Our proposed LOOCV method is compared with the widely used alternative methods [2]- [4] for choosing the shrinkage coefficients. The simulation results (averaged over 1000 repetitions for each training length T ) in Fig. 1 confirm that the LOOCV methods with and without the constraint ρ + τ = 1 produce the same results for the scaled identity target and they achieve performance almost identical to the OAS estimator [4], which was derived by assuming Gaussian data and identity target. The LW [2] and GLC [3] methods, which are equivalent for the scaled identity target here, do not perform well for very low sample support, but are able to approximate the oracle choice very well when more samples are available, which is consistent with the observations from [4]. All of these shrinkage designs significantly outperform the SCM, confirming the effectiveness of shrinkage for covariance matrix estimation. Recall that these methods were derived using different strategies and assumptions and have different analytical solutions.
Example 2: Shrinkage toward a nondiagonal target: We then consider an example of the linear model given by (35). For each training length, 1000 random realizations of Σ = HH † +σ 2 I are generated and estimated through training, where σ 2 = 0.1. The entries of H are independently generated from N (0, 1) and then fixed for the whole training process. Given H, T training samples are generated by y = Hx + z, with the entries of x and z generated independently from N (0, 1) and N (0, σ 2 ), respectively. In order to demonstrate the effectiveness of the LOOCV method for general shrinkage targets, we assume a scenario where H is slowly time-varying and the shrinkage target T 0 can be constructed as a well- conditioned estimate of a past covariance matrix where and the entries of ∆ are independently drawn from N (0, 0.2) and are fixed for each repetition. Specifically, we construct T 0 as the shrinkage estimate of Σ past using SCM and the scaled identity target. This construction is similar to the knowledgeaided target considered in [3] and the resulting T 0 is not diagonal. We assume that the numbers of samples used for estimating Σ and T 0 are both equal to T . The simulation results are included in Fig. 2 for the normalized MSE. It can be seen that the LOOCV methods generally achieve near-oracle performance and outperform the GLC method. Also, the nondiagonal shrinkage target achieves better performance than the scaled identity target. Example 3: Shrinkage with multiple targets: A multi-target example is illustrated in Fig. 3. An AR covariance matrix is estimated by shrinking SCM with three targets which can be represented by Eqn. (21) of [17]: T 1 = tr(R) N I, T 2 = Diag(R), and T 3 is a symmetric, Toeplitz matrix which was considered in [17]: where C i is a symmetric, Toeplitz matrix with unit entries on the i-th sub-and super-diagonals and zeros elsewhere. It is seen that multi-target shrinkage can significantly outperform single-target shrinkage with T 0 = tr(R) N I when the number of samples is large enough. For the oracle parameter choices, the unconstrained shrinkage design, which allows a larger set of shrinkage factors to be chosen, can noticeably outperform the design constrained by (71). However, when the proposed Normalized MSE (dB) Figure 3. MSE of covariance matrix estimation with multi-target (MT) shrinkage and LOOCV parameter choices. AR covariance matrix with N = 50, r = 0.9 is assumed. "LA" refers to the method proposed by Lancewicki and Aladjem [17]. Note that the LOOCV methods and the LA method achieve similar performance for this example.
LOOCV methods are used, the gap is significantly reduced. We can show that when the number of samples is small, using a more parsimonious design with constrained shrinkage coefficients or fewer shrinkage targets may achieve better performance. It is seen that the multi-target shrinkage method of [17] (indicated by "MT-LA" in Fig. 3) performs similarly to the LOOCV method for this example. Note that the method of [17] assumes SCM and shrinkage targets satisfying certain structures and does not apply directly to model-based covariance matrix estimation or more general shrinkage targets. Example 4: Application to MMSE estimation of MIMO channels. A potential application of the proposed technique is the design of MMSE estimator of MIMO channels. Consider a point-to-point MIMO system with N t transmitting antennas and N r receiving antennas. Let B be the length of the pilot sequence. The received signal matrix during the training stage is modelled as where Y ∈ C Nr×B is the received signal matrix, H ∈ C Nr×Nt the channel matrix, P ∈ C Nt×B the pilot matrix, and N ∈ C Nr×B the noise which is uncorrelated with H. Vectorizing Y in (86) gives where y = vec(Y), P = P T ⊗ I, h = vec(H), n = vec(N), vec(·) denotes vectorization, and ⊗ denotes Kronecker product. We assume a Rayleigh fading channel and denote by Σ h ∈ C NtNr×NtNr the covariance matrix of the channel vector h. We also assume that the disturbance n is complex Gaussian-distributed with a zero mean and identity covariance matrix. Given Σ h , the MMSE estimate of h from y can be computed as [41] h MMSE = Σ h P † ( PΣ h P † + I) −1 y.
The covariance matrix Σ h , which can be very large, must be estimated in order to compute h MMSE . In communication systems, h is not directly observable and thus the SCM estimator can not be directly applied to generate Σ h . One may estimate Σ h from least squares (LS) estimates of H, i.e., When orthogonal training signal with P = √ P I is applied, where P determines the power for training signals, it can be shown that Denote by h LS the vectorization of H LS . It can be shown that the covariance matrix of h LS is Therefore, if Σ hLS is estimated as Σ hLS , we can then use (91) to estimate Σ h as which can be used in (88). The estimation of Σ hLS can be achieved using the different shrinkage estimators introduced in this paper.
An example is shown in Fig. 4. The covariance matrix is assumed to be where Σ t and Σ r are, respectively, the transmitter side and receiver side covariance matrix, with entries given by r t = 0.7e −j0.9349π and r r = 0.9e −j0.9289π . While applying shrinkage to estimate Σ hLS , two shrinkage targets are tested: the identity matrix and the shrinkage estimate (with a scaled identity target) of a past covariance matrix. The second is considered based on the assumption that Σ h is slowly varying in time and a well-conditioned estimate of a past covariance matrix Σ past h can be available. In our simulations, Σ past h is modeled by randomly perturbing r t and r r in (93) and (94) by δ t and δ r whose real and imaginary parts are both randomly and uniformly generated from − 1 10 √ 2 , 1 10 √ 2 . The normalized MSE of channel estimation is defined as where h MMSE is the MMSE channel estimate obtained from (88) with the true channel covariance matrix replaced by its shrinkage estimate.
From the simulation results in Fig. 4, when the number of samples T of channel estimates is small, the MMSE channel estimator constructed using the SCM estimate of Σ h is poorer than the LS estimator which does not require any knowledge of Σ h . Therefore, an accurate estimate of the covariance matrix is necessary to exploit the potential of the MMSE channel estimator. Shrinkage with LOOCV choice of the shrinkage coefficients improves the performance of the MMSE channel estimator by providing a better estimate of Σ h . Two-target shrinkage can further enhance performance. Note that the multi-target method of [17] is not directly applicable to the shrinkage target used here. Similarly to [42], we do not exploit the Kronecker product structure in (92) and the exponential modeling of (93) and (94) while estimating the covariance matrix and similar trends can be observed when the channel covariance matrix follows different models such as those in [43], [44].
Example 5: Application to LMMSE signal estimation: Another example application is the design of linear minimum mean squared error (LMMSE) estimator [45], [46] for estimating the transmitted signal x in MIMO communications. The received signal is modeled by (35) and the LMMSE estimate of x is obtained as where we have assumed that x has identity covariance matrix and Σ y is the covariance matrix of y. The OLS-based covariance matrix estimation in Section II-D can be used to estimate Σ y in (96). In Fig. 5, we show an example where the shrinkage target T 0 is chosen as the diagonal matrix of the OLS estimate (38) of the covariance matrix. This results in a shrunk LMMSE signal estimator Orthogonal training of length T constructed from the discrete Fourier transform (DFT) matrix is assumed for the OLS channel estimate and finding the shrinkage coefficients using our proposed LOOCV method is achieved at a low complexity. The normalized MSE of signal estimation is defined as (98) Fig. 5 presents the simulation results averaged over 1000 random realizations of H for each T . It can be seen that the shrinkage estimate of the covariance matrix can lead to noticeable improvement of the MSE performance of signal estimation. The resulting performance can approach the oracle choice of (ρ, τ ) that minimizes the MSE of estimating x [35]. Note that in contrast to the cross-validation methods in [35] and [36] which choose shrinkage factors by a grid search for optimizing the signal estimation performance, the method proposed in this paper has an analytical solution and optimizes covariance matrix estimation. It also differs from [47] which targets the design of a signal estimator that shrinks the sample LMMSE filter toward the matched filter. Example 6: Application to MVDR beamforming: Finally, we show an example application to minimum variance distortionless response (MVDR) beamforming [31], [33]. We assume a N = 30-element uniform linear array (ULA) with halfwavelength spacing between neighboring antennas. As in [33], we assume that the desired complex Gaussian signal has an angle of arrival (AoA) of θ 0 = 0 }, all with an average power 10 dB higher than the desired signal. The noise is assumed to be additive white Gaussian noise (AWGN) with an average power 10 dB lower than the desired signal. The MVDR beamformer is given by where s is the steering vector of the desired signal and Σ is the covariance matrix of the received signal. We consider a practical scenario where the desired signal's steering vector suffers from an AoA error uniformly distributed in [−5 • , 5 • ] and Σ is estimated from the training samples by shrinking the SCM R toward the scaled identity matrix tr(R) N I. We Average SINR (dB) Figure 6. Average output SINR for a MVDR beamformer with AoA mismatch and the estimated covariance matrix. The results labeled by "SCM" is obtained by replacing Σ −1 in (99) with the pseudo-inverse of the SCM. Note that the LOOCV and OAS methods achieve almost the same performance, which is slightly better than the GLC method when T is very small.
focus on the low-sample-support case and compare the result with an approach that uses the pseudo-inverse of the SCM for computing w. The output signal-to-interference-and-noise ratio (SINR) averaged over 1000 repetitions are plotted in Fig. 6. It is seen that though the proposed approach targets covariance matrix estimation only and is not optimized for beamformer designs, it still provides noticeable gains as compared to the pseudo-inverse approach in the low-sample-support regime.

V. CONCLUSIONS
In this paper, we have introduced a leave-one-out crossvalidation (LOOCV) method for choosing the coefficients for linear shrinkage covariance matrix estimators. By employing a quadratic loss as the LOOCV prediction error, analytical expressions of the optimal shrinkage coefficients are obtained, which do not require a grid search of the parameters. As a result, the coefficients can be computed at low costs for the SCM-and OLS-based estimation of the covariance matrix. The LOOCV method is generic in the sense that it can be applied to different covariance matrix estimation methods and different shrinkage targets. Numerical examples show that it can approximate the oracle parameter choices in general and have wider applications than several existing analytical methods that have been widely applied.
Zero-mean signals have been assumed in this paper. When nonzero-mean signals are considered, our proposed approach may be applied after subtracting an estimate of the mean from the samples. However, the inaccuracy in the mean vector estimate may introduce extra errors to the covariance matrix estimation. Jointly estimating the mean and covariance matrix in a robust manner may be further explored. Other future work includes theoretical study of the properties of the proposed approach and low-complexity cross-validation schemes for choosing shrinkage factors for specific signal processing applications such as beamforming, space-time adaptive processing, correlation analysis, etc.