Improved Shrinkage Estimator of Large-Dimensional Covariance Matrix under the Complex Gaussian Distribution

Estimating the covariance matrix of a random vector is essential and challenging in large dimension and small sample size scenarios.*e purpose of this paper is to produce an outperformed large-dimensional covariance matrix estimator in the complex domain via the linear shrinkage regularization. Firstly, we develop a necessary moment property of the complex Wishart distribution. Secondly, by minimizing the mean squared error between the real covariance matrix and its shrinkage estimator, we obtain the optimal shrinkage intensity in a closed form for the spherical target matrix under the complex Gaussian distribution. *irdly, we propose a newly available shrinkage estimator by unbiasedly estimating the unknown scalars involved in the optimal shrinkage intensity. Both the numerical simulations and an example application to array signal processing reveal that the proposed covariance matrix estimator performs well in large dimension and small sample size scenarios.


Introduction
e problem of estimating the covariance matrix of a random vector arises in both multivariate statistical theory and various applications [1,2]. In large sample setting, where the dimension of a random vector is small and the sample size is large enough, the sample covariance matrix (SCM) is a reliable estimator of the real covariance matrix and is widely employed in many scenarios. However, suffering from the curse of dimensionality, the SCM becomes ill-conditioned or even singular in large dimension scenarios [3]. en, severe consequences may appear if the SCM remained as the covariance matrix estimator [4,5].
During the last two decades, scientists have proposed many regularization strategies to generate outperformed covariance matrix estimators in large dimension scenarios [6][7][8][9][10]. Among these, the linear shrinkage estimation is an effective strategy to inspire a well-conditioned covariance matrix estimator when the dimension is large compared to the sample size [11,12]. When the prior information of the covariance structure is available, the linear shrinkage estimator is modeled as a linear combination between the SCM and a proper target matrix. In the existing literature, the target matrices, which are usually formed through structuring the SCM in line with the prior information, include the spherical target and others, such as the diagonal target, the Toeplitz rectified target, and the tapered SCM [13].
With the aid of prior information, the linear shrinkage estimator can always outperform the SCM when the involved tuning parameter is carefully selected [14]. erefore, one of the crucial difficulties in linear shrinkage estimation is to determine the optimal tuning parameter which is also called as the shrinkage intensity. By minimizing the mean squared error (MSE) between the shrinkage estimator and the real covariance matrix, the optimal tuning parameter can be expressed in a closed form for an arbitrary target. However, it is comprised of unknown scalars which involve the expectation operator and the real covariance matrix, leading to a chief difficulty in generating an available shrinkage estimator. In particular, when the data follow a specific distribution such as Gaussian distribution or elliptical distribution, the expectation can be further calculated [11]. What is noteworthy is that the optimal tuning parameter owns different expressions under different distributions, even for the same target [15]. erefore, it is necessary to discover the specific properties of shrinkage intensity under some typical distributions. In many applications, such as array signal processing, the data come from the complex domain. ere have been some related research studies on it. In [16], a linear shrinkage estimator for Toeplitz rectified target is developed under the complex Gaussian distribution, whereas the involved unknown scalars in the optimal tuning parameter are not unbiasedly estimated, resulting in a suboptimal covariance matrix estimator [17]. In [18], a linear shrinkage estimator is proposed via lowcomplexity crossvalidation under an arbitrary complex distribution.
erefore, when the data follow a specific distribution, the linear shrinkage estimator could be further improved by making full use of the distribution information.
In this paper, we further research the linear shrinkage estimator under the complex Gaussian distribution. e target matrix is chosen as the spherical target which has been widely studied under the real number field [6,11,14]. e optimal tuning parameter is obtained by minimizing the MSE. We remind that the above optimal tuning parameter involves both the expectation operator and the real covariance matrix. By developing a novel moment property of the complex Wishart distribution, we can calculate the expectation operator. en, the optimal tuning parameter turns to be only related to some unknown scalars concerning the real covariance matrix. A popular approach is adopted by replacing these unknown scalars with their estimates to obtain an available tuning parameter. Furthermore, good estimates of unknown scalars can benefit the corresponding available tuning parameter and the corresponding shrinkage estimator [11,19]. e main contributions of this paper are summarized as three-fold: (1) A necessary moment property of the complex Wishart distribution is developed. On this basis, the optimal tuning parameter for the spherical target is analytically expressed under the complex Gaussian distribution. (2) All the unknown scalars involved in the optimal tuning parameter are unbiasedly estimated. en, the corresponding available linear shrinkage estimator under the complex Gaussian distribution is proposed.
(3) e performance of the proposed covariance matrix estimator is verified with comparison to the existing estimators in numerical simulations and an example application to adaptive beamforming. e rest of this paper is organized as follows: Section 2 formulates the linear shrinkage estimation under the complex Gaussian distribution as a quadratic programming problem. e optimal solution is analytically obtained. Section 3 unbiasedly estimates the relative unknown scalars and subsequently proposes a new shrinkage estimator for the spherical target. Section 4 provides some numerical simulations and an example application for verifying the performance of the proposed covariance matrix estimator. Section 5 concludes.

Notations.
e notation C m is the set of all m-dimensional complex column vectors, and H n is the set of all n × n Hermitian matrices. e symbol E denotes the mathematical expectation. e bold symbols 0 and 1 respectively denote the column vectors having all entries 0 and 1 with an appropriate dimension. e symbol I n denotes the n × n identity matrix. For a matrix A, A H and ‖A‖ respectively denote its conjugate transpose and Frobenius matrix norm. For a squared matrix A, A − 1 and tr(A) respectively denote its inverse and trace. For two real numbers a and b, a∧b and a∨b respectively mean the maximum and minimum of a and b.

Formulation and the Optimal Solution
Assume a p-dimensional random vector x ∈ C p follows the complex Gaussian distribution CN(0, Σ), where Σ is the unknown covariance matrix. Let x 1 , x 2 , . . . , x n ∈ C p be an independent and identically distributed (i.i.d.) sample, then the sample covariance matrix S is defined by For an arbitrary prespecified target matrix T ∈ H p which represents an aspect prior information of the real covariance matrix structure [20], the linear shrinkage estimator of covariance matrix Σ is modeled as where w ∈ [0, 1] is the tuning parameter which is also called shrinkage intensity [21]. Because S and T are Hermitian, we have Σ∈ H p for an arbitrary w ∈ [0, 1]. To find the optimal shrinkage intensity, we employ the MSE criterion: 2 ] is a constant. erefore, the optimal shrinkage intensity can be obtained through solving the following optimization problem: It is worth noticing that the objective function in optimization problem (5) is a convex quadratic function of w, and the optimal shrinkage intensity can be expressed in a closed form as follows: 2 Mathematical Problems in Engineering Furthermore, for the spherical target T � (tr(S)/p)I p , the optimal shrinkage intensity w * 0 given by (6) becomes Denote the matrix We can obtain the following moment property of the complex Gaussian distribution.
Furthermore, when a random matrix W � (w ij ) p×p follows complex Wishart distribution CW(I p , n) with degree of freedom n, we have By taking expectation on both sides, we have where b 1 � n 2 and b 2 � n. For a random matrix W which follows complex Wishart distribution CW(Σ, n) with degree of freedom n, let Σ � GG H ; then, we have G − 1 WG − H ∼ CW(I p , n). In the same manner, we can obtain E tr 2 (W) � n 2 tr 2 (Σ) + ntr Σ 2 .
erefore, we have By (10) and (16), equality (9) holds. □ Theorem 1. When the target matrix is T � (tr(S)/p)I p , the optimal shrinkage intensity under the complex Gaussian distribution is Proof. By plugging equalities (10) and (16) into (7), we have erefore, we can obtain By Cauchy-Schwarz inequality, we have w * ∈ [0, 1]. Hence, we have By eorem 1, the corresponding optimal linear shrinkage estimator is We remind that the optimal shrinkage estimator concerns with the real covariance matrix. us, it is unavailable in practical applications. Despite this, it provides a theoretical optimal value for evaluating the available ones.

Available Linear Shrinkage Estimator
Theorem 2. Under the complex Gaussian distribution, the unbiased estimates of tr(Σ 2 ) and tr 2 (Σ) are, respectively, given by Proof. Because the inverse matrix of E n is we have erefore, we can obtain revealing that α and β are unbiased estimates of tr(Σ 2 ) and tr 2 (Σ), respectively. rough plugging the unbiased estimates given by (22) into the optimal shrinkage intensity w * , we can obtain the available shrinkage intensity: erefore, the available linear shrinkage estimator is e linear shrinkage estimator given by (27) is positive definite even when the dimension exceeds the sample size, except that Σ degenerates into the SCM.

Numerical Simulations and Adaptive Beamforming
In this section, we provide some numerical simulations and an example application to adaptive beamforming for verifying the performance of the proposed covariance matrix estimator. e proposed linear shrinkage estimator is denoted as T1cg.
e linear shrinkage estimator corresponding to the spherical target matrix in [18] is denoted as T1cv.

Numerical Simulations.
As mentioned before, an accurate shrinkage intensity estimate can benefit the linear shrinkage estimator. In this section, we compare the proposed available shrinkage intensity and the existing one based on crossvalidation in [18] to reveal the advantage of the proposed shrinkage estimator. In our simulations, the real covariance matrix is Σ � (σ ij ) p×p with σ ij � t |i− j| .
(28) e model parameter t is set to be 0.5, resulting in the real covariance matrix being close to a spherical structured matrix. e data come from the complex Gaussian distri- bution CN(0, Σ).
e MSE of each available shrinkage intensity relative to the optimal intensity given by (17) is computed by averaging 5 × 10 4 Marlo runs. Figure 1 reports the MSEs of available shrinkage intensities versus the sample size and the dimension. We can see that the MSEs of available shrinkage intensities in T1cv and T1cg decrease as the sample size or the dimension gets larger. Because the proposed T1cg by plugging in the unbiased estimates of unknown scalars employs the complex Gaussian distribution information, it outperforms the T1cv based on the nonparameter approach.

Adaptive Beamforming.
In this section, we apply the proposed covariance estimators to array signal processing. Specifically, we consider a uniform linear array (ULA) which consists of p sensors with half-wavelength spacing. At time t � 1, . . . , n, the received signal can be modeled as where θ 0 and θ k are the directions of desired signal s 0 (t) and interference signals s k (t), respectively, a(θ 0 ) and a(θ k ) are the corresponding array responses, and n(t) is the noise [23]. en, the minimum variance distortionless response (MVDR) beamformer is expressed as e covariance matrix Σ in (30) is unknown and suggested to be replaced with its estimate Σ [24]. en, the corresponding output signal-to-interference-plus-noise ratio (SINR) is where w is the estimated beamformer corresponding to Σ. e covariance matrix estimator with larger output SINR is preferred in array signal processing.
In our simulations, we assume the desired signal has an angle of arrival of θ 0 � 5°with power σ 2 0 � 10 dB, and the interference signals come from the directions − 10°, 0°, 10°{ } with power 8 dB. For each covariance estimator, the corresponding output SINR is approximated by averaging 5 × 10 6 repetitions. Figures 2 and 3 report the SINR and corresponding elapsed time of the adaptive beamformers based on different covariance matrix estimators under the complex Gaussian scenario, where the noise n(t) follows the complex Gaussian distribution with power 0 dB. In Figure 2, the dimension is p � 60 and the sample size ranges from 20 to 120. In Figure 3, the sample size is n � 60 and the dimension ranges from 20 to 120. Our observations and analyses are summarized as follows: (1) Even though enjoying the lowest computation cost, the classic covariance estimator SCM has an unsatisfactory performance in small sample size scenarios. erefore, it is not an ideal covariance matrix estimator any more in these scenarios. (2) e SINR based on each covariance matrix estimator increases when the sample size gets larger but   Figures 2 and 3. When the received signal comes from the complex non-Gaussian distribution, the proposed estimator T1cg performs inferior to the existing T1cv in adaptive beamforming. It is worthy noticing that both T1cg and T1cv have analytical expressions, and there are (np 2 + np + 3p 2 ) multiplication in the proposed estimators T1cg and (np 2 + 3np + 5p 2 ) multiplication in T1cv. erefore, the proposed T1cg can always enjoy a lower computation complexity than T1cv.
On the whole, by employing additional distribution information, the proposed estimator T1cg outperforms T1cv in the complex Gaussian scenario and enjoys a comparable performance with T1cv in the complex non-Gaussian scenario. Moreover, the proposed estimator T1cg always enjoys an advantage over T1cv in computation cost.

Conclusion
In this paper, we have proposed a new covariance matrix estimator via linear shrinkage procedure under the complex Gaussian distribution. rough calculating the moment of Wishart distribution, we obtain the optimal shrinkage intensity for the spherical target. Furthermore, the involved unknown scalars are unbiasedly estimated. Subsequently, we propose the corresponding available linear shrinkage estimator. Numerical simulations and application to adaptive beamforming show that the proposed covariance matrix estimator is outperformed compared with the existing estimators. In future work, we will investigate the Cramér-Rao bound for the linear shrinkage estimation and develop nonlinear shrinkage estimation of the large-dimensional covariance matrix.

Data Availability
All data included in this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interests regarding the publication of this paper. Mathematical Problems in Engineering 7