{'fraction_non_alphanumeric': 0.0518572336754155, 'fraction_numerical': 0.016377561832107288, 'mean_word_length': 4.2981555733761025, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Regularized Estimation of Kronecker-Structured Covariance Matrix

This paper investigates regularized estimation of Kronecker-structured covariance matrices (CM) for complex elliptically symmetric (CES) data. To obtain a well-conditioned estimate of the CM, we add penalty terms of Kullback-Leibler divergence to the negative log-likelihood function of the associated complex angular Gaussian (CAG) distribution. This is shown to be equivalent to regularizing Tyler’s fixed-point equations by shrinkage. A sufficient condition that the solution exists is discussed. An iterative algorithm is applied to solve the resulting fixed-point iterations and its convergence is proved. In order to solve the critical problem of tuning the shrinkage factors, we then introduce three methods by exploiting oracle approximating shrinkage (OAS) and cross-validation (CV). When the training samples are limited, the proposed estimator, referred to as the robust shrinkage Kronecker estimator (RSKE), has better performance compared with several existing methods. Simulations are conducted for validating the proposed estimator and demonstrating its high performance.


I. INTRODUCTION
C OVARIANCE matrix (CM) estimation is a fundamental problem in many fields, such as adaptive detection and remote sensing [1]- [6]. The most common CM estimator is the sample covariance matrix (SCM), which is the maximum likelihood estimator (MLE) of the CM for Gaussian data. However, the SCM suffers poor performance for data with outliers or heavily-tailed distributions due to the lack of robustness. Such data are common in many applications [7]- [9] and often they can be described by complex elliptically symmetric (CES) distributions [10]. For instance, a subclass of CES distributions known as the compound-Gaussian (CG) distributions have been widely used in modeling the clutter returns in radar applications [11]- [14]. To tackle the heavilytailed data, one class of approaches is to censor the training samples with the aim to exclude outliers from the CM estimation [15]- [22]. Another class of methods is based on robustification. In particular, for CES distributions, various robust CM estimators based on the M-estimator have been developed and characterized [23]- [29]. With such estimators, outlying training samples are usually given small weights when an estimate of the CM is produced.
The SCM also requires an abundant number of samples to achieve satisfactory performance. Many modern applications involve high-dimensional variables whose statistical characteristics remain stationary over a short observation period, where the large sample support assumption does not hold. Regularization provides an effective strategy to improve the CM estimation for addressing the challenge of training shortage. In particular, a class of linear shrinkage algorithms have been introduced [30]- [33] and their integration into robust CM estimators for CES-distributed data have been thoroughly investigated in the recent works [34]- [37]. These algorithms estimate the CM by shrinking an estimate of the CM Σ toward a better-conditioned target matrix T. There can be various choices for Σ and T. For example, one can choose Σ as the SCM and Tyler's estimator [25] for the Gaussian and non-Gaussian data, respectively. Moreover, different types of target matrices T can be used, including the identity and diagonal targets. The linear shrinkage estimators can reduce the requirement of samples and provide positive-definite CM estimates. The choice of shrinkage factors is a fundamental problem for shrinkage estimators. Various criteria and methods have been studied. In particular, Ledoit and Wolf (LW) propose an approach that asymptotically minimizes the mean squared error (MSE) [31]. Then [33] improves the LW approach using the Rao-Blackwell theorem and designs the Rao-Blackwell Ledoit and Wolf (RBLW) estimator. The oracle approximating shrinkage (OAS) method is proposed in [33]. Both estimators have close-form expressions and are easily computed. The problem of determining the shrinkage factors can also be cast as a model selection problem and thus generic model selection techniques such as cross-validation (CV) [38] can be applied. The main challenges faced by CV include the choice of the cost function and the heavy computational cost in its direct implementation. Some efforts are made in [39], [40] to address these challenges for linear shrinkage estimators with unstructured CM.
Exploiting the structural knowledge about the CM can also significantly reduce the number of unknown parameters and improve its high estimation accuracy under limited training data [41]- [51]. Kronecker-structured CM is widely seen in many scenarios such as modeling multiple-input multipleoutput (MIMO) wireless communication channels [52]- [54] and the clutter in polarization-space-time adaptive processing (PSTAP) [55]- [58]. Particularly, [41] proposes a robust estimator for Kronecker-structured CM and proves that a globally optimal solution can be found, [42] proposes a majorization minimization (MM) solution to the Kronecker maximum likelihood estimator (KMLE), and [59] introduces the maximum likelihood (ML) estimation of Kronecker-structured CM with the presence of Gaussian clutter. An extension of KMLE is also studied for compound Gaussian clutter with inverse Gamma-distributed texture and Kronecker normalized sample covariance matrix (KNSCM) is proposed in [60] to estimate the CM. Both KMLE and KNSCM provide considerable performance with sufficient samples but they can still noticeably suffer performance degradation with low sample supports.

A. Contributions
In this paper, we consider the estimation of Kroneckerstructured CM for CES data under low sample supports. We investigate a robust shrinkage Kronecker estimator (RSKE) that aims to achieve well-conditioned 1 and highly accurate CM estimates. With RSKE, the structural knowledge is exploited together with robustfication and regularization techniques. Based on the findings of the previous studies in [10], [34], [36], [39], [62], [63] and others, we investigate the existence of RSKE, its iterative solver and convergence, and also the choice of the shrinkage factors. The contributions of this paper can be summarized as follows: 1) We show that the RSKE can be interpreted as the minimizer of a negative log-likelihood function penalized by the Kullback-Leibler divergence. Based on this, the condition for the existence of RSKE is established under some mild assumptions, which provides insights to the relationship between the dimensionality, sample size and shrinkage factors. 2) We study an iterative solver involving two fixed-point equations to find RSKE and prove its convergence. Following the majorization-minimization framework, we prove the monotonic decrease of the penalized loglikelihood function over iterations. We show that, with fixed shrinkage factors and arbitrary positive-definite initial matrices, the iterative solver converges. 3) We address the critical challenge of shrinkage factor choice in order to exploit the potential of RSKE. We introduce data-driven methods that automatically tune the linear shrinkage factors, based on oracle approximating shrinkage (OAS) and cross-validation (CV). The OAS method adopts a minimum MSE (MMSE) criterion and plug-in estimates of the oracle shrinkage factors. For the CV methods, we start with a quadratic loss for leave-one-out CV (LOOCV) and derive analytical solutions of the shrinkage factors which can approach the performance of the oracle solutions that minimize the MSE of CM estimation. The complexities of these different methods are analyzed. We show that they exhibit different complexities, performance, and require different levels of a prior knowledge and thus may be useful in different scenarios. Note that one main challenge of the CV is its computational cost. However, the proposed CV-based methods have close-form expressions for the shrinkage factors which can address this challenge. Owing to such expressions, the complexity of the CV-based RSKE is about the same as that of the KMLE.

B. Organization
The remainder of this paper is organized as follows. Section II introduces the RSKE, its existence and iterative solution. Section III gives the choices of the shrinkage factors. Section IV presents simulation results to show the performance of CM estimation. Finally, Section V gives the conclusions.

II. ROBUST SHRINKAGE KRONECKER ESTIMATOR (RSKE)
In this section, we introduce the robust shrinkage estimator for Kronecker-structured covariance matrices. We first discuss the motivation, then discuss the condition for its existence, and finally introduce the iterative solver and its convergence property.

A. Motivation
Let y be an N -dimensional zero-mean random vector following a CES distribution, whose probability density function (p.d.f.) is of the form [10], [29], where (·) H is the conjugate transpose, g(·) denotes the density generator, C N,g a normalizing constant, and R the normalized covariance matrix with its trace Tr(R) = N . Note that R is also known as the scatter matrix [10], [29]. We consider the case that R can be expressed as the Kronecker product of two matrices, i.e., R R A ⊗ R B , where R A ∈ S NA ++ , R B ∈ S NB ++ , ⊗ denotes the Kronecker product, and S n ++ denotes the set of Hermitian, positive definite matrices of size n × n. In order to define R A and R B uniquely, we assume Tr(R A ) = N A and Tr(R B ) = N B , which is consistent with the assumption Tr(R) = N.
This paper considers the estimation of R from Y = {y l } L l=1 which is a set of independent and identically distributed (i.i.d.) samples of y. The normalized samples {x l = y l y l } L l=1 , which belong to a complex unit N -dimensional sphere, follows the complex angular Gaussian (CAG) distribution [10], [29]. The joint distribution function of {x l } L l=1 is expressed as [10] p where det(·) denotes the determinant. After omitting some additive constants and scaling, the negative log-likelihood function of such a joint distribution is given by where R A ∈ S NA ++ , R B ∈ S NB ++ and we have used the fact that in the likelihood function. The above cost function L 0 ( R A , R B ) is non-convex in the classical definitions but is jointly g-convex (geodesic-convex) with respect to R A and R B [41]. Minimizing this cost function produces the KMLE [42], [63]. In the low-sample-support cases, the solution of KMLE can suffer from significant errors and illconditioning. For many applications such as beamforming and spectral estimation [64]- [69], the inverse of the CM estimate is required. Inverting an erroneous, ill-conditioned CM estimate can bring enormous errors. This motivates the design of accurate, well-conditioned CM estimators.

B. Regularization via KL Divergence Penalty
In this subsection, we introduce a penalized estimator that promotes well-conditioned estimates of the sub-CMs R A and R B . We adopt penalty terms of the Kullback-Leibler divergence for Gaussian distributions [70], i.e., where X, Y ∈ S N ++ . As shown in [71], the KL divergence D KL (X, I N ) can effectively constrain the condition number of X. We thus add the penalty terms to the negative log-likelihood function in (3) to promote well-conditioned estimates R A and R B , where α A = NB ρA 1−ρA and α B = NAρB 1−ρB with ρ A ∈ [0, 1) and ρ B ∈ [0, 1). Ignoring some additive constants which are irrelevant to R A and R B , the penalized negative log-likelihood function is obtained as By adding the penalty terms which are convex, the obtained objective function is also g-convex w.r.t. R A and R B . This guarantees that all local minimizers of L( R A , R B ) are also globally optimal, following [41, Proposition 1].
Minimizing the penalized log-likelihood function by setting (5b) In the above, we have defined where y (i) l denotes the ith entry of y l and unvec NBNA (·) reshapes a vector into a N B × N A matrix as shown above. Therefore, the solution to (5), if exists, can be interpreted as the minimizer of the penalized negative log-likelihood function (4). These fixed-point equations interestingly have the same form as the linear shrinkage estimators for unstructured CM [31], [33]- [37]. Following these work, we refer to the resultant CM estimator as the robust shrinkage Kronecker estimator (RSKE), with shrinkage factors ρ A and ρ B . The KMLE [42] can be obtained as a special case of RSKE by letting ρ A = ρ B = 0.
It should be noted that in [63], estimators that exploit robustification and shrinkage for the unstructured CM and robust estimators for the Kronecker-structured CM have been thoroughly studied via the geodesic convexity. The KL divergence penalty has also been exploited in [36] for robust estimation of unstructured CM. We here extend these studies to the estimation of Kronecker-structured CM by simultaneously exploiting robustification and shrinkage.

C. Existence of RSKE
In this subsection, we examine the conditions under which the RSKE exists. When ρ A and ρ B are small, it is possible that the cost function (4) tends to −∞ on the boundary of the set S NA ++ and S NB ++ , i.e., (4) becomes unbounded below and there is no solution to the fix-point equations of (5). The existence of the shrinkage Tyler's estimator for unstructured CM has been studied in [36], where the relationship between the shrinkage factors, sample size, and dimensionality is revealed. By establishing the condition under which the cost function tends to +∞ on the boundary of the set of positivedefinite, Hermitian matrix, the minimum shrinkage factor for the existence of the CM estimator is obtained [36]. This result, however, can not directly determine the conditions of the two shrinkage factors affecting each other. In this work, we follow [36, Theorem 3] and its proof to study the RSKE. We first construct auxiliary functions by which the penalized negative log-likelihood function (4) can be lowerbounded. The two auxiliary functions have a similar form as (15) in [36]. Thus, using the same treatment of [36], we can examine the conditions for the auxiliary functions tending to +∞ at the boundary. Based on the results, we can obtain the following sufficient condition for the existence of a solution to the RSKE: Proposition 1: The cost function (4) has a finite lower bound over the set of positive-definite R A and R B , i.e., a solution to (5) exists if the following conditions are satisfied: (1) None of r j,l and c i,l is an all-zero vector, where r j,l ∈ C NA×1 denotes the jth row of Y l and c i,l ∈ C NB ×1 denotes the ith column of Y l ; such that for any proper subspace S A ⊂ C NA×1 and S B ⊂ C NB ×1 in the space of length-N A and -N B vectors, respectively, where , 1 x denotes the indicator function.
Proof : See Appendix A.
In general, the above conditions require that the number of samples to be sufficiently large, and the samples are evenly spread out in the whole space.
Corollary 1: If the samples are evenly spread out in the whole space, such that Rearranging (9), one has which is exactly (8a). Similarily, we have (8b). Remark 1: Condition (2) in Corollary 1 shows the relationship between the shrinkage factors, the number of samples L, and the dimension of the sub-CMs N A and N B . In general, a larger shrinkage factor ρ A is required when L decreases or N A increases. Moreover, Condition (2) can be easily checked. For example, when β 1 = 1 and β 2 = 0, , which agrees with the result in [35], [36] for the case of unstructured CM.

D. Iterative Solver and Its Convergence
Similarly to [34]- [37], we solve (5) by applying the process below, which involves two fixed-point iterations: where and R

(k)
A and R (k) B denote the estimates of the sub-CMs at the kth iteration. In this paper, we choose the initial CM estimates as R (0) It is useful to examine the convergence property of the above iterative estimator which generalizes Tyler's estimator [25] and its shrinkage extension [34], [36], [37] to the case of Kronecker-structured CM. The works [25], [34], [36], [37] assume unstructured CM and thus their solutions can be characterized by a single fixed-point equation. The convergence of the iterative process for Tyler's estimator is proved in [25] by examining the fixed-point iterations. For the shrinkage extension of Tyler's estimator, the convergence is proved in [34] by applying the concave Perron-Frobenius theory, in [36] by applying the majorization-minimization theorem, and in [37] by applying the monotone bounded convergence theorem. For the Kronecker-structured CM, though the case of the KMLE has been studied in [41], in this work we incorporate shrinkage into the estimator and the convergence has not been analyzed earlier to the authors' best knowledge. Exploiting the majorization-minimization framework [72], we have the following proposition that establishes the converging property of the fixed-point iterations in (10).
Proposition 2: The fixed-point iterations in (10) converge to the solution of (5) for arbitrary positive-definite initial matrices R Proof : See Appendix B. Remark 2: The iterations in (10) can be terminated by using a distance metric where B and · denotes the Frobenius norm. This metric measures the variation of the solution over iterations. Then a stopping criterion can be set to terminate the iterations when D( R (k+1) , R (k) ) < δ (13) or k > K max is met, where δ denotes a preset threshold and K max the maximum number of iterations allowed.

III. CHOICE OF THE SHRINKAGE FACTORS
The performance of the RSKE can be optimized by properly choosing the shrinkage factors ρ A and ρ B . In this section, we propose two different classes of choices, based on oracle approximating shrinkage (OAS) and leave-one-out cross validation (LOOCV), respectively, to provide solutions with different performance and complexity.

A. The KOAS Method
In [34], an OAS strategy for choosing the shrinkage factor for unstructured CM is derived by exploiting the MMSE criterion and plug-in estimates. We can extend this strategy to the RSKE. The choice of the two shrinkage factors will be decoupled into separate problems to enable a low-complexity solution. Following [34], we begin by assuming that the true CM R A and R B are already "known". Then, we choose to choose the shrinkage factors (ρ A , ρ B ) that achieve the MMSE of the covariance matrix estimates as and where E{·} denotes the mathematical expectation and The following proposition extends the OAS solution of [34] to the Kronecker-structured CM. Proposition 3: The shrinkage factors that achieve the MMSE are given as .
where we have used the assumption Tr( Proof : See Appendix C. In practice, R A and R B in (17) are unknown. Similarly to [34], we propose to replace them by their trace-normalized estimates R A and R B , such as the KNSCM [60] and KMLE [42]. We will show the performance of the resulting shrinkage factors (ρ A,KOAS , ρ B,KOAS ), referred to as the Kronecker OAS (KOAS) choice, in Section IV. Note that, if N A = 1 or N B = 1, the Kronecker-structured CM reduces to the unstructured CM and (17) agrees with (17) in [34]. If ρ A,KOAS < 0 is produced, we then truncate it to ρ A,KOAS = 0. If ρ A,KOAS ≥ 1, we simply set the covariance matrix estimate to be the shrinkage target matrix. The treatments are similar for ρ B,KOAS < 0 and ρ B,KOAS ≥ 1 and also the LOOCVbased choices of the shrinkage factors to be introduced in the next subsection.

B. The LOOCV Methods
We next provide alternative methods for choosing the shrinkage factors based on LOOCV. In order to achieve good performance and complexity tradeoff, the cost for LOOCV must be carefully chosen. In this work, we extend the quadratic cost used in [39] to obtain data-driven, analytical solutions. Note that [39] considers unstructured CM for Gaussian data, whereas this paper considers Kronecker structured CM estimation with elliptically distributed data for which iterative solvers are required. Two LOOCV solutions, i.e., CV-I and CV-II, will be introduced in this subsection.
1) CV-I choice: Let Σ A and Σ B be two positive-definite, Hermitian matrices. Define the following cost function where the expectation is with respect to Y = unvec NB NA (y), Proof : See Appendix D. Inspired by Proposition 4, we aim to estimate the cost function in (18) and then minimize it over the shrinkage factors. This may be achieved using different strategies, e.g., [31]. In this paper, we apply the LOOCV strategy [38] to estimate J A (Σ A ) and J B (Σ B ) and minimize them to determine the shrinkage factors. With the standard LOOCV, the samples Y are repeatedly split into two sets. For the lth split, the samples in the training set Y l (with the lth sample y l omitted from Y) are used for producing shrinkage CM estimates {Σ A , Σ B } and the remaining sample y l is used for constructing {S A , S B } to estimate J A (Σ A ) and J B (Σ B ). The standard LOOCV process requires the iterative estimator to be applied for L times for each pair of candidate shrinkage factors (ρ A , ρ B ), which can lead to significant complexity, especially when grid search of (ρ A , ρ B ) is conducted. In order to address this complexity challenge, we propose alternative solutions by using proxy estimators so that closed-form expressions can be found for the optimized shrinkage factors.
Similarly to KOAS, we first assume that the covariance matrices are "known" and consider estimates of the covariance matrices from the samples Y l = {Y j , j = l} as where Following [39], we adopt the quadratic cost functions below: where S (l) Substituting (20a) into (22a), the cost function can be rewritten as We treat J A,CV (ρ A ) as a proxy of J A (Σ A ) and choose the shrinkage factor ρ A as the minimizer of (24) as: Similarly, we choose ρ B as Tr Alternative expressions can be derived for (25) and (26) to reduce the computational costs. Let Recalling (21a) and (23a), we have A and I NA are all Hermitian matrices. By using (28), we have Substituting (29) into (25), we obtain (31a) on the next page to quickly evaluate the shrinkage factors ρ A,CV-I . Similarly, we can obtain (31b) there for ρ B,CV-I , where The shrinkage factors determined by (31) still require the true CM R A and R B to be known to compute (27), (30), and (23). Similarly to KOAS, we propose to substitute them by their trace-normalized estimates R A and R B . We refer to the resultant solutions as the CV-I choice.
2) CV-II choice: As shown later in the simulation, the performance of KOAS and CV-I depend on the choice of the plug-in estimates R A and R B . In order to address this challenge, we propose to choose R A and R B used in the CV-I method as the current CM estimates at each iteration, i.e., R where where R

(k)
A and R Similarly to the CV-I choice, we then adopt the following LOOCV cost functions The CV-II choice of the shrinkage factors are then given as the minimizers of (34a) and (34b) and their analytical expressions are given by (36) on the following page. As such, the shrinkage factors in (10) are updated at each iteration as ρ , for updating the CM in the iterative process in (10). As such, computing the shrinkage factors costs O(N it L(N 2 A + N 2 B )). Note that the CV-II choice does not require extra plug-in CM estimates.
It can be seen that, ignoring the cost for finding the plugin CMs, the complexity of finding the shrinkage factors is dominated by that of iteratively updating the CMs in (10). In contrast to CV-II, both KOAS and CV-I require extra cost for finding the plug-in CMs. For example, finding the KNSCM and KMLE requires complexity of O(L( it is the number of iterations for KMLE. In particular, the KMLE has the same order of complexity as the RSKE.
Remark 4: The CV-II choice iteratively updates the shrinkage factors by treating the current estimates of the CMs R (36b)

IV. SIMULATION RESULTS
In this section, we show the performance of the proposed RSKE estimators. We compare the proposed estimators with the following CM estimators: SCM [78], KMLE [42], [60], and KNSCM [60]. For the RSKE, in addition to the KOAS and CV choices of the shrinkage factors, the oracle choice of the shrinkage factors is also considered, which minimizes the NMSE defined in (37) at each iteration under the assumption that the true CM is known.
In order to evaluate the CM estimation accuracy, we use the following normalized mean-square error (NMSE) as the performance metric [79]: Consider an example with the autoregressive (AR) covariance matrices, whose (i, j)th entries are given as which has been widely considered for evaluating CM estimation techniques [34], [39], [60]. Here we first set N A = 6, N B = 6, ǫ A = 0.1 and ǫ B = 0.9. Then the samples are generated according to y l = √ τ l u l , l = 1, 2, · · · , L, where the texture τ l follows a Gamma distribution [75] of shape parameter ν and scale parameter 1/ν, i.e., τ l ∼ Γ(ν, 1/ν), u l ∼ CN (0, R). The generated samples {y l } follow a zeromean CES distribution. The estimated sub-CMs R (k) in (11) are initialized as identity matrices for simplicity but other initilizations can produce similar results.
As CV-II leads to shrinkage factors varying over iterations and a proof of the convergence of the resulting RSKE is unavailable, we resort to numerical studies. 200 Monte-Carlo experiments are performed. Fig. 1 shows the average of D( R (k) , R (k−1) ) varying with the iterations where CV-II is applied. We can see that R A and R B of CV-II tend to stabilize over iterations. This indicates that CV-II also leads to converging solutions. The distance between R (k) and R also decreases over iterations, indicating improved accuracy of the CM estimate over iterations. Furthermore, a small number of iterations is sufficient for approaching the performance limit of the RSKE. In the rest of this section, for terminating the iterations, we choose the threshold δ in (13) as 10 −3 and K max = 15. see that the NMSE of the proposed estimators approaches that of the oracle estimator as N A (and N ) increases. Fig. 5 shows the condition number of the estimated CM of RSKE (with CV-I, CV-II, KOAS), KMLE and KNSCM. Here we set N A = 8, N B = 8, ǫ A = 0.1 and ǫ B = 0.9. We set the plug-in estimator for CV-I and KOAS as KNSCM. One can see that the proposed CV-I, CV-II and KOAS algorithms yield CM estimates which are better-conditioned than those with KNSCM and KMLE, especially when the number of samples is small. As they also improve the NMSE, it is expected that the RSKE with the proposed shrinkage factor choices can improve the performance for applications where the inverse of the CM is required, such as beamforming and spectral estimation applications.

V. CONCLUSIONS
In this paper, we investigate a robust, iterative shrinkage estimator, which is referred to as RSKE, for estimating the CM with the Kronecker product structure. The RSKE can be obtained by minimizing a negative log-likelihood function penalized by Kullback-Leibler divergence and interpreted by integrating linear shrinkage into the fixed-point iterations. The conditions for the existence of the RSKE are investigated and the convergence of the iterative solver is investigated. We also introduce three methods for choosing the shrinkage factors by exploiting oracle approximating shrinkage (OAS) and crossvalidation (CV), respectively. Compared with the state-ofthe-art estimators, the RSKE achieves more accurate CM estimation and improves the condition number by significantly reducing the number of unknown parameters and integrating shrinkage into the robust estimation. APPENDIX A PROOF OF PROPOSITION 1 In this appendix, we examine the conditions under which a solution to (5) exists by constructing two auxiliary functions to lowerbound the cost function in (4). Let λ be the eigenvalues of R A and R B . Then we have where we have utilized Jensen's inequality in the last step. Similarly, we have  Here we have assumed that none of r j,l and c i,l is an all-zero vector, such that r H j,l R −1 A r j,l = 0, c H i,l R −1 B c i,l = 0, ∀i, ∀j, ∀l. Then let us define the following auxiliary functions: A , where β 1 + β 2 = 1 and β 1 , β 2 ∈ [0, 1]. From (39) and (40), we have In the following, we check the conditions under which F 1 ( R A ) → +∞ and Here we consider the case with r A ≥ 1, i.e., there exists at least one eigenvalue diverging, following [36], in order to examine the condition for F 1 ( R A ) → +∞ at the boundary of feasible set for R A .
Define  and A 1 = I NA , and hence a ρ = a ′ ρ = β 1 N A and a 1 = +∞, a ′ 1 = 0 2 . Note also that 2 (aρ, a ′ ρ ) and (a 1 , a ′ 1 ) are respectively defined for ρ(s) and αh 1 (s) according to [36,Definition 2] for any ǫ > 0, where o(·) denotes the higher order infinitesimal and ϕ Then we impose the same condition as the first line 3 of (36) in [36, Appendix A], i.e., Under this condition, G 1 ( R A ) goes to zero , i.e., F 1 ( R A ) → +∞ on the boundary of positive-definite and Hermitian R A [36]. Letting ǫ → 0 and rearranging the terms, one has for arbitrary d = 1, · · · , N A − 1. Intuitively, this requires that the samples are evenly spread in the subspaces spanned by the eigenvectors of R A . The condition (42) is then rewritten in a general form as (7a). Similarily, we have (7b).
In summary, we have obtained conditions (7a) and (7b) under which the cost function (4) tends to positive infinity at the boundary of the set of positive definite and Hermitian matrix. By [36,Lemma 1], these also give a sufficient condition that a solution to (5) exists.

APPENDIX B PROOF OF PROPOSITION 2
In this Appendix, we prove the convergence of the proposed iteration process, following the methodology of [36], [41]. By the concavity of the logarithm function, one has log x ≤ log a + x a − 1, ∀a > 0. The equality holds when x = a. Then we have 3 The second line of (36) in [36, Appendix A] is always met since a 1 = +∞ in this paper.
where the equality holds when R A = R (k) A . We then construct the surrogate function Recalling (43), we have and the equality holds when It is easy to verify that the minimizer of (44) is exactly (10a) by setting the gradient of (44) with respect to R A to zero. It follows that Therefore, Then define Similarly, we can verify that the minimizer of (49) is exactly (10b), and where the equality holds when R B = R (k) It follows that Combining (48) and (52), we have i.e., the penalized log-likelihood function L( R A , R B ) in (4) is decreasing with iterations.
Since L( R A , R B ) is g-convex, its minimizer exists and denote it by B )} is bounded by an infimum. Then according to the monotone convergence theorem [80], the sequence will converge to the infimum as k increases, i.e., ( R

APPENDIX C PROOF OF PROPOSITION 3
We here complete the proof by exploiting results from random matrix theory. Following [81], when the true covariance matrix R A and R B are known, the oracle shrinkage factor ρ ⋆ B , i.e., the solution to (15), is given by where Re(·) denotes the real part and and C B is defined by (30). The resulting optimal shrinkage estimate can be interpreted as the projection of the true CM onto the linear space spanned by C B and I NB .  It is easy to see that z l 2 = 1 and {z l } are independent of each other. Moreover, the whiten vectors {z l } are isotropically distributed [82] and satisfy [33], [34] E z l z H l = 1 N I N , E z H l Λz l 2 = Tr R 2 + Tr 2 (R) N (N + 1) , B . We then reshape z l into a matrix satisfying which can be easily verified by vectorizing both sides of (57). In order to determine the shrinkage factor for the robust shrinkage estimator of unstructured CM, [34] analyzed the feature of Z l where it reduces to a vector. We here extend the analysis to the more general case of matrix-valued Z l by exploiting random matrix theory and properties of Kronecker product. Let z (i) l be the ith entry of z l . From (56), one has This indicates that {z (i) l } N i=1 are i.i.d. with zero mean and variance 1/N . Consequently, we have Note that D −1 y l 2 2 = y H l (R A ⊗ R B ) −1 y l , and we have Note that E{·}, Re(·) and Tr(·) are exchangeable to each other. Substituting (60) into (55), one has From [83], [84], we have Since {z l } L l=1 are i.i.d, we have Therefore, (55) can be rewritten as Utilizing [83, Lemma 1.1] and subsitituting (62), (63) into (64), E 1 is obtained as (65) in the following page. Substituting (65) and (61) into (54), (17b) is obtained. Similarly, we can have the optimal ρ ⋆ A , i.e., (17a). The resulting expressions of ρ ⋆ A and ρ ⋆ B can be used to produce the KOAS choice ρ A,KOAS and ρ B,KOAS by plugging estimates of R A and R B into (17).
APPENDIX D PROOF OF PROPOSITION 4 This proposition can be proven by combining the results in Appendix C. Recalling (19), (59) and (60), we have Moreover, (18) can be rewritten as By setting the derivative of (67a) and (67b) with respect to Σ A and Σ B to zero, we have the minimizer of (18) as Σ A = R A and Σ B = R B .