A large covariance matrix estimator under intermediate spikiness regimes

The present paper concerns large covariance matrix estimation via composite minimization under the assumption of low rank plus sparse structure. In this approach, the low rank plus sparse decomposition of the covariance matrix is recovered by least squares minimization under nuclear norm plus $l_1$ norm penalization. This paper proposes a new estimator of that family based on an additional least-squares re-optimization step aimed at un-shrinking the eigenvalues of the low rank component estimated at the first step. We prove that such un-shrinkage causes the final estimate to approach the target as closely as possible in Frobenius norm while recovering exactly the underlying low rank and sparsity pattern. Consistency is guaranteed when $n$ is at least $O(p^{\frac{3}{2}\delta})$, provided that the maximum number of non-zeros per row in the sparse component is $O(p^{\delta})$ with $\delta \leq \frac{1}{2}$. Consistent recovery is ensured if the latent eigenvalues scale to $p^{\alpha}$, $\alpha \in[0,1]$, while rank consistency is ensured if $\delta \leq \alpha$. The resulting estimator is called UNALCE (UNshrunk ALgebraic Covariance Estimator) and is shown to outperform state of the art estimators, especially for what concerns fitting properties and sparsity pattern detection. The effectiveness of UNALCE is highlighted on a real example regarding ECB banking supervisory data.


Introduction
Estimation of population covariance matrices from samples of multivariate data is of interest in many high-dimensional inference problems -principal components analysis, classification by discriminant analysis, inferring a graphical model structure, and others. Depending on the different goal the interest is sometimes in inferring the eigenstructure of the covariance matrix (as in PCA) and sometimes in estimating its inverse (as in discriminant analysis or in graphical models). Examples of application areas where these problems arise include gene arrays, fMRI, text retrieval, image classification, spectroscopy, climate studies, finance and macro-economic analysis.
The theory of multivariate analysis for normal variables has been well worked out (see, for example, Anderson (1984)). However, it became soon apparent that exact expressions were cumbersome, and that multivariate data were rarely Gaussian. The remedy was asymptotic theory for large samples and fixed, relatively small, dimensions. However, in recent years, datasets that do not fit into this framework have become very common, since nowadays the data can be very high-dimensional and sample sizes can be very small relative to dimension.
The most traditional covariance estimator, the sample covariance matrix, is known to be dramatically ill-conditioned in a large dimensional context, where the process dimension p is larger than or close to the sample size n, even when the population covariance matrix is well-conditioned. Two key properties of the matrix estimation process assume a particular relevance in large dimensions: well conditioning (i.e. numerical stability) and identifiability.
Both properties are crucial for the theoretical recovery and the practical use of the estimate. A bad conditioned estimate suffers from collinearity and causes its inverse, the precision matrix, to dramatically amplify any error in the data. A large dimension may cause the impossibility to identify the unknown covariance structure thus hampering the interpretation of the results.
Regularization approaches to large covariance matrices estimation have therefore started to be presented in the literature, both from theoretical and practical points of view (see Fan et al. (2016) for an exhaustive overview). Eigenvalue regularization approaches include linear (Ledoit and Wolf, 2004) and nonlinear shrinkage (Ledoit and Wolf (2015), Lam et al. (2016)).
A different approach is based on the assumption of a low rank plus sparse structure for the covariance matrix: where L * is low rank with rank r < p, S * is positive definite and sparse with at most s nonzero off-diagonal elements, and Σ * is a positive definite matrix. The generic covariance estimator Σ can be written asΣ where W is an error term. The error matrix W may be deterministic or stochastic, as explained in Agarwal et al. (2012). If the data are Gaussian andΣ is the unbiased sample covariance matrix Σ n , W is distributed as a re-centered Wishart.
In Fan et al. (2013), a large covariance matrix estimator, called POET (Principal Orthogonal complEment Thresholding), is derived under this assumption. POET combines Principal Component Analysis for the recovery of the low rank component and a thresholding algorithm for the recovery of the sparse component. The underlying model assumptions prescribe an approximate factor model with spiked eigenvalues (i.e. growing with p) for the data, thus allowing to reasonably use the truncated PCA of the sample covariance matrix. Furthermore, at the same time, sparsity in the sense of Bickel and Levina (2008a) is imposed to the residual matrix. The latent rank r is chosen by the information criteria of Bai and Ng (2002).
Indeed, rank selection represents a relevant issue: if p is large, setting a large rank would cause the estimateΣ to be non-positive definite, while setting a small rank would cause a too relevant variance loss. In the discussion of Fan et al. (2013), Yu and Samworth point out that the probability to underestimate the latent rank does not asymptotically vanish if the eigenvalues are not really spiked at rate O(p). In addition, we note that POET systematically overestimates the proportion of variance explained by the factors (given the true rank) because the eigenvalues of Σ n are more spiky than the true ones (as showed in Ledoit and Wolf (2004)).
POET asymptotic consistency holds given that a number of assumptions is satisfied. The key assumption is the pervasiveness of latent factors, which causes the PCA of Σ n to asymptotically identify the eigenvalues and the eigenvectors of Σ * as p diverges. The results of Fan et al. (2013) provide the convergence rates of the relative norm ofΣ P OET − Σ * (defined as ||Σ P OET −Σ * || Σ = p −1/2 ||Σ * − 1 2Σ P OET Σ * − 1 2 − I p || F ro ), the maximum norm ofΣ P OET −Σ * and the spectral norm ofŜ P OET −S * . Under stricter conditions,Ŝ P OET andΣ P OET are proved to be non-singular with probability approaching 1.
At the same time, a number of non-asymptotic methods has been presented. In Chandrasekaran et al. (2011) the exact recovery of the covariance matrix in the noiseless context is first proved. The result is achieved minimizing a specific convex non-smooth objective, which is the sum of the nuclear norm of the low rank component and the l 1 norm of the sparse component. In Chandrasekaran et al. (2012), which is an extension of Chandrasekaran et al. (2011), the exact recovery of the inverse covariance matrix by the same numerical problem in the noisy graphical model setting is provided. The authors prove that, in the worst case, the number of necessary samples in order to ensure consistency is n = O p 3 /r 2 , even if the required condition for the positive definiteness of the estimate is p ≤ 2n.
An approximate solution to the recovery and identifiability of the covariance matrix in the noisy context is described in Agarwal et al. (2012). Even there, the condition p ≤ n is unavoidable, for standard results on large deviations and non-asymptotic random matrix theory. An exact solution to the same problem, based on the results in Chandrasekaran et al. (2012), is then shown in Luo (2011b). The resulting estimator is called LOREC (LOw Rank and sparsE Covariance estimator) and is proved to be both algebraically and parametrically consistent in the sense of Chandrasekaran et al. (2012).
In Chandrasekaran et al. (2012) algebraic consistency is defined as follows Definition 1.1 A pair of symmetric matrices (S, L) with S, L ∈ R p×p is an algebraically consistent estimate of the low rank plus sparse model (2) for the covariance matrix Σ * if the following conditions hold: 1. The sign pattern of S is the same of S * : sign(S ij ) = sign((S * ) i,j ), ∀i, j. Here we assume that sign(0) = 0.
2. The rank of L is the same as the rank of L * .
3. Matrices L + S, S and L are such that L + S and S are positive definite and L is positive semidefinite.
Parametric consistency holds if the estimates of (S, L) are close to (S * , L * ) in some norm with high probability. In Chandrasekaran et al. (2012) such norm is LOREC shows several advantages respect to POET. The most important is that the estimates are both algebraically and parametrically consistent, while POET provides only parametric consistency. In spite of that, LOREC suffers from some drawbacks, especially concerning fitting properties. What is more, the strict condition p ≤ n is required, while POET allows for p log(p) ≫ n.
For these reasons, we propose a new estimator, UNALCE (UNshrunk ALgebraic Covariance Estimator), based on the unshrinkage of the estimated eigenvalues of the low rank component, which allows to improve the fitting properties of LOREC systematically. We assume that the non-zero eigenvalues of L * and Σ * are proportional to p α , α ∈ [0, 1] (the so called generalized spikiness context). Under the assumption that the maximum number of non-zeros per row in S * , called "maximum degree", is O(p δ ) (with δ ≤ 1 2 ), we prove that our estimator possesses a non-asymptotic error bound admitting that n is as small as O(p 3δ ).
We derive absolute bounds depending on α for the low rank, the sparse component, and the overall estimate, as well as the conditions for rank consistency, positive definiteness and invertibility. In this way we provide a unique framework for covariance estimation via composite minimization under the low rank plus sparse assumption.
The remainder of the paper is organized as follows. In Section 2 we first define ALCE (AL-gebraic Covariance Estimator) with the necessary assumptions for algebraic and parametric consistency, and then define UNALCE, proving that the unshrinkage of thresholded eigenvalues of the low rank component is the key to improve fitting properties as much as possible given a finite sample, preserving algebraic consistency. In Section 3 we propose a new model selection criterion specifically tailored to our model setting. In Section 4 we provide a real Euro Area banking data example which clarifies the effectiveness of our approach. Finally, in Section 5 we draw the conclusions and discuss the most relevant findings.
2 Numerical estimation and spiked eigenvalues: the ALCE approach 2.1 The model First of all, we recall the definitions of the matrix norms used throughout the paper. Let us define a p × p symmetric positive-definite matrix M . We denote by λ i (M ), i = 1, . . . , p, the eigenvalues of M in descending order. Then we recall the following norms definitions: 1. element-wise: (a) L 0 norm: ||M || 0 = p i=1 p j=1 ½(m ij = 0), which is the total number of non-zeros.
2. induced by vector: (a) ||M || 0,v = max i≤p j≤p ½(m ij = 0), which is the maximum number of non-zeros per column, defined as the maximum "degree" of M ; (c) spectral norm: ||M || 2 = λ 1 (M ); 3. Schatten: (a) nuclear norm of M , here defined as the sum of the eigenvalues of M : Let us suppose the population covariance matrix of our data is the sum of a low rank and a sparse component. A p-dimensional random vector x is said to have a low rank plus sparse structure if its covariance matrix Σ * satisfies the following relationship: where: 1. L * is a positive semidefinite symmetric p × p matrix with at most rank r ≪ p; 2. S * is a positive definite p × p sparse matrix with at most s ≪ p(p − 1)/2 nonzero off-diagonal elements and maximum degree s ′ .
According to the spectral theorem, we can write Let us suppose that the p × 1 random vector x is generated according to the following model: where f is a r × 1 random vector with E(f ) = 0 r , V (f ) = I r and ǫ is p × 1 random vector with E(ǫ) = 0 p ,V (ǫ) = S * . The random vector x is thus assumed to be zero mean, without loss of generality. Given a sample x i , i = 1, . . . , n, Σ n = 1 It is easy to observe that x follows a low rank plus sparse structure: Assuming p ≤ n, it is also useful to recall that for n → ∞ If we assume a normal distribution for f and ǫ, the above equality is true for any fixed n and the matrix W := Σ n − (L * + S * ) is distributed as a re-centered Wishart noise. In any case, the normality assumption is not essential for our setting.

Nuclear norm plus l 1 norm heuristics
Under model (2), the need rises to develop a method able at the same time to consistently estimate the covariance matrix Σ * as well as to catch the sparsity pattern of S * and the spikiness pattern of the eigenvalues of L * simultaneously. Such estimation problem is stated as min L,S where ||S|| 0,of f = p−1 i=1 p j=i+1 ½(s * ij = 0) (because the diagonal of S is preserved as in Fan et al. (2013)). This is a combinatorial problem, which is known to be NP-hard, since both rank(L) and ||S|| 0,of f are not convex.
The tightest convex relaxation of problem (7), as shown in Fazel (2002), is where ψ and ρ are non-negative threshold parameters, and ||S|| of f, The use of nuclear norm for covariance matrix estimation was introduced in Fazel et al. (2001).
The feasible set of (8) is the set of all p × p positive definite matrices S and all p × p positive semi-definite matrices L.
From a statistical point of view, (8) is a penalized least squares heuristics, composed by a smooth least squares term ( 1 2 ||(L + S) − Σ n || 2 F ro ) and a non-smooth composite penalty (ψ||L|| * + ρ||S|| 1 ). The choice of (8) allows to lower the condition number of the estimates and the parameter space dimensionality simultaneously. The optimization of (8) requires the theory of non-smooth convex optimization provided by Rockafellar (2015) and Clarke (1990) (the solution algorithm is reported in the Supplement).
In principles, different losses could be used, like Stein's one (Dey and Srinivasan, 1985).
However, the classical Frobenius loss does not require normality and is computationally appealing. The study of different fitting terms, including the ones performing eigenvalue regularization, is left to future research.
From an algebraic point of view, (8) is an algebraic matrix variety recovery problem. In the noisy covariance matrix setting described in equation (3), matrices L * and S * are assumed to come from the following sets of matrices: B(r) is the variety of matrices with at most rank r. A(s) is the variety of (element-wise) sparse matrices with at most s nonzero elements, where support(S) is the orthogonal complement of ker(S).
In Chandrasekaran et al. (2011) the notion of rank-sparsity incoherence is developed, which is defined as the uncertainty principle between the sparsity pattern of a matrix and its row/column space. Denoting by T (L) and Ω(S) the tangent spaces to B(r) and A(s) respectively, the following rank-sparsity incoherence measures between Ω(S * ) and T (L * ) are defined: µ(Ω(S * )) = max In order to identify T (L * ) and Ω(S * ), we need quantities ξ(T (L * )) and µ(Ω(S * )) to be as small as possible, because the smaller they are, the better is the decomposition. The product µ(Ω(S * ))ξ(T (L * )) is the rank-sparsity incoherence measure and bounding it controls both for identification and recovery.
The described approach was first used for deriving LOREC estimator in Luo (2011b).
Therein, the reference matrix class imposed to Σ * is which is the class of positive definite matrices having uniformly bounded eigenvalues. In the context so far described, Luo proves that L and S can be identified and recovered with bounded error, and the rank of L as well as the sparsity pattern of S are exactly recovered.
The key model-based results for deriving LOREC consistency bounds are a lemma by Bickel and Levina (2008a) for the sample loss in infinity (element-wise) norm: and a lemma by Davidson and Szarek (2001) for the sample loss in spectral norm: We stress that (15) strictly requires the assumption p ≤ n.
From a theoretical point of view, LOREC approach presents some deficiencies and incongruities. Differently from POET approach, where the sparsity assumption is imposed to the sparse component S * , LOREC approach imposes it directly to the covariance matrix Σ * . As a consequence, the assumption Σ * ∈ Σ * (ǫ 0 ) (see (13)) is necessary and causes, jointly with the identifiability assumptions, uncertainty on the underlying structure of Σ * .
In fact, assuming uniformly bounded eigenvalues may conflict with the main necessary identifiability condition: the transversality between Ω and T . Since the eigenvalue structures of Σ * and S * are somehow linked, requiring class (13) for Σ * may cause S * to be not enough sparse, and simultaneously the row/column space of L * to have high values of incoherence, because we have no spiked eigenvalues. This may result in possible non-identifiability issues.

ALCE estimator
Let us suppose that the eigenvalues of Σ * are intermediately spiked with respect to p. This equals to assume the generalized spikiness of latent eigenvalues in the sense of Yu and Samworth (Fan et al. (2013), p. 656): Assumption 2.1 All the eigenvalues of the r × r matrix p −α B ′ B are bounded away from 0 for all p and α ∈ [0, 1].
If p is finite, Assumption 2.1 is equivalent to state that for some δ α > 0. Hence, we aim to study the properties of the covariance estimates obtained by heuristics (8) under the generalized spikiness assumption in a non-asymptotic context.
In order to do that, we need to study the behaviour of the model-based quantity P (||Σ n − Σ * ||), which is the only probabilistic component. We bound P (||Σ n − Σ * || ∞ ) exploiting the property ||.|| ∞ ≤ ||.|| 2 . Therefore, our aim is to show that which is verified if it holds with very high probability (C 1 , C 2 , C 3 and C are positive constants). Exploiting the consistency norm of Chandrasekaran et al. (2012), which is it follows from (17) that with very high probability (see Luo (2011b) for technical details).
In order to reach this goal, we need to impose that the following assumptions hold in our finite sample context.
Assumption 2.2 is needed to ensure algebraic consistency. In fact, an identifiability condition for problem (8), as shown in Theorem 2.1, is ξ(T (L * ))µ(Ω(S * )) ≤ 1 54 . According to Chandrasekaran et al. (2011), it holds . It descends that ξ(T (L * )) = 1 with δ = 0 in the worst case scenario and ξ(T (L * )) = r p with δ = 1 2 in the best case scenario, under the condition The assumption δ ≤ α is made to prevent the violation of Assumption 2.1 under the condition λ r (L * ) > C 2 ψ ξ 2 (T ) of Theorem 2.1 . Assumption 2.3 is necessary to ensure that the large deviation theory can be applied to f ij , ǫ ij ′ and f ij ǫ ij ′ for all i ≤ n, j ≤ r and j ′ ≤ p. Assumption 2.4 is necessary to apply the results of Bickel and Levina (2008a) on the thresholding of the sparse component, which prescribe that S * must be well conditioned with uniformly bounded diagonal elements. We stress that the maximum degree s ′ must be bounded. This condition is stronger than the corresponding one in Fan et al. (2013), which prescribes max i≤p j≤p |s * ij | q < c 4 , q ∈ [0, 1], c 4 > 0. This is the price to pay for algebraic consistency, because our assumption ensures Assumption 2.5 prescribes that the latent rank is infinitesimal with respect to p and the sample size n is possibly smaller than p, but not smaller than δ 4 p 3 2 δ . The need for this assumption rises throughout the proof of (17), and to ensure consistency with Assumption The inequality (20), under Assumptions 2.2, 2.4 and 2.5, boils down to 1 k L δ 4 p α−2δ < δ 2 p δ , which holds if Assumption 2.6 is respected. As a consequence, we can allow for ||S * || 2 ≤ All outlined propositions must hold for finite values of p, α and n. The following theorem provides a non-asymptotic consistency result particularly useful when p is not that large and α < 1, because the absolute rate of Σ n under POET assumptions, O( p √ n ), may be too strong and prevent consistency.
In addition, suppose that the minimum singular value of L * (λ r (L * )) is greater than C 2 ψ ξ 2 (T ) and the smaller absolute value of the nonzero entries of S * , S min,of f , is greater than C 3 ψ µ(Ω) . Then, with probability greater than 1−C 4 p −C 5 , the pair (L,Ŝ) minimizing (8) recovers the rank of L * and the sparsity pattern of S * exactly: rank(L) = rank(L * ) and sign(Ŝ) = sign(S * ).
Moreover, with probability greater than 1 − C 4 p −C 5 , the matrix losses for each component are bounded as follows: We call the resulting covariance estimator ALCE ( We stress that the conditions λ r (L * ) > C 2 ψ ξ 2 (T ) and S min,of f > C 2 ψ µ(Ω) under Assumptions 2.2 and 2.5 become λ r (L * ) > C 2 p α and S min,of f > C 3 p α−2δ respectively. The latter in turn leads to (20) under Assumption 2.6. Therefore, the resultant model setting is fully consistent with Assumptions 2.1 and 2.4.
Our results are non-asymptotic in nature, thus giving some probabilistic guarantees for finite values of p and n. This fact depends on the algebraic consistency properties, which ensure the exact recovery of the rank and the sparsity pattern, with some parametric guarantees for the estimation error in g γ norm (see (19)). The shape of the probabilistic bound ψ depends on the assumption n ≥ δ 4 p 3δ , which is also necessary in order to ensure asymptotic consistency, as the following Corollary shows.
We stress that the probabilistic bound ψ decreases to 0 as long as n ≥ δ 4 p 3 2 δ . This assumption leads to overcome the restrictive condition p ≤ n, since δ ≤ 1/2. In addition, an immediate consequence of (17) is reported in the following Corollary (the proof is reported in the Supplement).
Corollary 2.2 Letλ r,ALCE be the r−th largest eigenvalue ofL ALCE . Then under the assumptions of Theorem 2.1λ r,ALCE > C 1 p α with probability approaching 1 for some C 1 > 0.
Corollary 2.2 states thatL U N ALCE is rank-consistent as the latent degree of spikiness α is not smaller than the maximum degree of the residual component δ. If α = 1, we fall back to the POET setting. If α = 0, r = log(1) = 0, and we fall back to the pure sparsity estimator of Bickel and Levina (2008a).
A representative selection of the latent eigenvalue and sparsity patterns admitted under described conditions is reported in the Supplementary Material. We emphasize that the algebraic consistency does no longer force the latent eigenvalues to scale to p, provided that the spectral norm of the residual component is scaled accordingly. In general, it is needed that the minimum latent eigenvalue and absolute nonzero residual entry are large enough to ensure consistency, but, unlike POET, they can be scaled to p α , α < 1. The exponent α plays the role of an adaptive spikiness degree.
In particular, if we increase α (ceteris paribus), both λ r (L * ) and S min,of f must be larger to ensure identifiability. The same happens if p increases, because, according to Chandrasekaran et al. (2011), both ξ(T (L * )) and µ(Ω(S * )) depend inversely on p. On the contrary, to ensure consistency, if r increases L * can have less spiked eigenvalues, while if s ′ increases S min,of f can be smaller. This occurs because, according to Chandrasekaran et al. (2011), ξ(T (L * )) and µ(Ω(S * )) directly depend on r and s ′ respectively.
From Theorem 2.1, we can derive with probability larger than 1 − C 1 p −C 2 the following bounds forΣ ALCE : which hold if and only if λ min (Σ * ) > φ. The same bounds hold for the inverse covariance estimateΣ −1 ALCE with the same probability: given that λ min (Σ * ) ≥ 2φ.
Within the same framework, we can complete our analysis with the bounds forŜ. From S is positive definite if and only if λ min (S * ) > φ S .Ŝ −1 has the same bound ofŜ if and only To sum up, by ALCE estimator we offer the chance to recover consistently a relaxed spiked eigen-structure, thus overcoming the condition p ≤ n, even using the sample covariance matrix as estimation input (the ratio p/n directly impacts on the error bound). Our bounds are in absolute norms, and reflect the underlying degree of spikiness α. Our theory relies on the probabilistic convergence of the sample covariance matrix under the assumption that the data follow an approximate factor model with a sparse residual. If p and n are in a proper relationship, both parametric and algebraic consistency are ensured.

UNALCE estimator: a re-optimized ALCE solution
A key aspect of Theorem 2.1 is that the two losses in L * and S * are bounded separately. This fact results in a negative effect on the overall performance ofΣ ALCE , represented by the loss ||∆ Σ || 2 , since ||∆ Σ || 2 is simply derived as a function of ||∆ L || 2 and ||∆ S || 2 according to the triangle inequality ||∆ Σ || 2 ≤ ||∆ L || 2 + ||∆ S || 2 . Therefore, the need rises to correct for this drawback, re-shapinĝ Σ ALCE , as ALCE approach is somehow sub-optimal for the whole covariance matrix.
We approach this problem by a finite-sample analysis, which could be referred to as a re-optimized least squares method. We refer to the usual objective function (8) the l 1 norm of S excluding the diagonal entries, consistently with POET approach. We define Y pre and Z pre the last updates in the gradient step of the minimization algorithm of (8) (see the Supplement for more details). Y pre and Z pre are the two matrices we condition upon in order to derive our finite-sample re-optimized estimates.
Suppose thatB(r) andÂ(ŝ) are the recovered varieties ensuring the algebraic consistency of (8). One might look for the solution (say (L N ew ,Ŝ N ew )) of the problem min L∈B(r),S∈Â(ŝ) where T L(L, S) stands for Total Loss. The sample covariance matrix follows the model Σ n = L * + S * + W , given a sample of p−dimensional data vectors x i , i = 1, . . . , n. Our problem essentially is: which pair L ∈B(r), S ∈Â(ŝ) satisfying algebraic consistency shows the best approximation properties of Σ n ?
We prove the following result.
Theorem 2.2 Suppose thatB(r), andÂ(s) are the recovered matrix varieties, and that such that its off-diagonal elements are the same asŜ ALCE andΣ N ew such that its diagonal elements are the same asΣ ALCE respectively. Then, the minimum conditioning on Y pre and Z pre is achieved if and only if whereψ > 0 is any prescribed threshold parameter.
Theorem 2.2 essentially states that the sample total loss (27) is minimized if we un-shrink the eigenvalues ofL ALCE (re-adding the thresholdψ). We call the resulting overall estimator Σ new =L new +Ŝ new UNALCE (UNshrunk ALgebraic Covariance Estimator). We stress the importance of conditioning on Y pre and Z pre . Since Y pre and Z pre are the matrices minimizing 1 2 ||(L+S)−Σ n || 2 F ro conditioning on the contemporaneous minimization ofψ||L|| * +ρ||S|| 1 , our finite-sample re-optimization step aims to re-compute min ||Σ n − (L + S)|| 2 F ro once removed the effect of the composite penalty.
As proved in the Supplement (which we refer to for the details), problem (27) can be decomposed in two problems: one in L and one in S. The problem in L is solved by the covariance matrix formed by the topr principal components of Y pre , which belongs by construction toB(r) and is equal toÛ ALCE (D ALCE +ψI r )Û ′ ALCE =L U N ALCE . The problem in S collapses to the problem in L under the prescribed assumptions on the off-diagonal elements of S U N ALCE (which causesŜ U N ALCE ∈Â(ŝ)), and on the diagonal elements ofΣ U N ALCE . The new estimate of the diagonal of S * is simply the difference between the diagonal of the original Σ ALCE and the diagonal of the newly computedL U N ALCE . Note that our re-optimization step depends entirely on Σ n , as Y pre and Z pre are Σ n -dependent.
Four consequences of Theorem 2.2 are reported in Corollary 2.3.

Corollary 2.3
The gains in terms of spectral loss forL U N ALCE ,Ŝ U N ALCE in comparison toL ALCE ,Ŝ ALCE respectively are all strictly positive and bounded byψ: The gains in terms of Frobenius norm are all strictly positive and bounded as follows: Two further relevant consequences of Theorem 2.2 are reported in Corollary 2.4.

Corollary 2.4
The gain in terms of spectral sample total loss forΣ U N ALCE respect toΣ ALCE is strictly positive and bounded byψ: The gain in terms of Frobenius sample total loss forΣ U N ALCE respect toΣ ALCE is strictly positive and bounded by √ rψ: The following result compares the losses ofΣ U N ALCE andΣ ALCE from the target Σ * .
Theorem 2.3 Conditioning on Σ n , the gains in terms of spectral loss and Frobenius loss for Σ U N ALCE respect toΣ ALCE are strictly positive and bounded as follows: The rationale of the reported claims is the following. We accept to pay the price of a nonoptimal solution in terms of nuclear norm (we allow to increment ||L|| * by rψ) but we have a best fitting performance for the whole covariance matrix, decrementing the squared Frobenius loss ofΣ by a quantity bounded by rψ 2 . The l 1 norm of S excluding the diagonal, ||Ŝ|| of f , is unvaried, while the norm ||S|| 1 (included the diagonal) is decreased by a quantity bounded by √ rψ.
The following Corollary extends our framework to the performance of (Σ U N ALCE ) −1 .
Corollary 2.5 The gains in terms of spectral loss and Frobenius loss forΣ −1 U N ALCE respect toΣ −1 ALCE are strictly positive and bounded as follows: The outlined results allow to improve the estimation performance given the finite sample.
However, the non-asymptotic bounds forL U N ALCE ,Ŝ U N ALCE andΣ U N ALCE are exactly the ones ofL ALCE ,Ŝ ALCE andΣ ALCE . UNALCE improves systematically the fitting performance of ALCE, inheriting all its algebraic and parametric consistency properties. The proofs of all theorems and corollaries can be found in the Supplement.
Finally, we study how the necessary conditions to ensure the positive definiteness of UN-ALCE estimates evolve respect to the ALCE ones. The following Corollary holds.
We stress that the improvement of the condition forL U N ALCE is numerically much larger than the worsening of the conditions forŜ U N ALCE andΣ U N ALCE .

A new model selection criterion: MC
In empirical applications, the selection of thresholds ψ and ρ in equation (8) requires a model selection criterion consistent with the described estimation method. The motivation rises from the consistency norm g γ used in Luo (2011b) (see (18)). Our aim is to detect the optimal threshold pair (ψ, ρ) in respect to the spikiness/sparsity trade-off. In order to exploit (18) with model selection purposes, we need to make the two terms comparable, i.e., the need of rescaling both arguments of g γ rises.
First of all, we note that if all the estimated latent eigenvalues are equal, we have ||L|| * = r||L||. As the condition number ofL increases, we haver||L|| > ||L|| * . As a consequence, the quantityr||L|| acts as a penalization term against the presence of too small eigenvalues.
In order to compare the magnitude of the two quantities, we divide the former by the trace ofL, estimated byθtrace(Σ n ), and the latter by the trace ofŜ, estimated by (1−θ)trace(Σ n ).
Our maximum criterion M C can be therefore defined as follows: where γ = ρ ψ is the ratio between the sparsity and the spikiness threshold. MC criterion is by definition mainly intended to catch the proportion of variance explained by the factors. For this reason, it tends to choose quite sparse solutions with a small number of non zeros and a small proportion of residual covariance, unless the non-zero entries ofŜ are prominent, as Theorem 2.1 prescribes. The M C method performs considerably better than the usual cross-validation using H-fold Frobenius loss (used in Luo (2011b)). In fact, minimizing a loss based on a sample approximation like the Frobenius one causes the parameterθ to be shrunk too much. The threshold setting which shows a minimum for M C criterion (given that the estimateΣ is positive definite) is the best in terms of composite penalty, taking into account the latent low rank and sparse structure simultaneously.

A Euro Area banking data example
This Section provides a real example on the performance of POET and UNALCE based on a selection of Euro Area banking data. We acknowledge the assistance of the European Central   In order to to obtain a POET estimate, we exploit the algebraic consistency ofΣ U N ALCE setting the rank to 6 and we perform cross-validation for threshold selection. The results are reported in Table 2, where we note that the number of estimated non-zeros is 404 (0.56%).
Apparently, one could argue that POET estimate is better: the estimated proportion of common variance is 0.6123, and the proportion of residual covariance is 0.0161. On the contrary, UNALCE method outputsθ = 0.3247 andρ corr = 0.1687. A relevant question arises: how much is the true proportion of variance explained by the factors? In fact, a so high latent proportion variance, which depends on the use of PCA with 6 components, causes the residual covariance proportion to be very low. Therefore, POET procedure gives a priori a preference for the low rank part. This pattern does not change even if we choose a lower value for the rank.
On the contrary, the UNALCE estimate, which depends on a double-step iterative thresholding procedure, requires a larger magnitude of the non-zero elements in the sparse component. In fact, the proportion of lost covariance during the procedure is here 29.39%. As a consequence, via rank/sparsity detection UNALCE shows better approximation properties respect to POET: its Sample Total Loss is relevantly lower than the one of the competitor (0.337 VS 0.645).
For UNALCE method, the covariance structure appears so complex that a relevant proportion of residual covariance is present. This allows us to explore the importance of variables, that is to explore which variables have the largest systemic power (i.e. the most relevant communality) or the largest idiosyncrasy (i.e. the most relevant residual variance).
In Figure 2  In Figure 3 we report the top 6 variables by estimated number of non-zero residual covariances. They are mainly credit-based variables: financial assets through profit and loss, central banks impaired assets, allowances to credit institutions and non-financial corporations, cash.
These variables are related to the largest number of other variables.

Variable
Maximum degree Financial assets designated at fair value through profit or loss 34 Central banks Impaired assets [gross carrying amount] 25 Credit institutions Collective allowances for incurred but not reported losses 20 Other financial corporations Collective allowances for incurred but not reported losses 19 Cash, cash balances at central banks and other demand deposits 16 Other financial corporations Specific allowances for financial assets, collectively estim. 16 We note that those variables have a marginal power in the explanation of the common covariance structure, and are much less relevant for supervisory analysis than the previous five.
In conclusion, our UNALCE procedure offers a more realistic view of the underlying  In particular, we prove that UNALCE can effectively recover the covariance matrix even in presence of spiked eigenvalues with rate O(p), exactly as POET estimator does, allowing n to be as small as O(p 3δ ), where δ is the maximum degree of the sparse component. In addition, we prove that the recovery is actually effective even if the latent eigenvalues show an intermediate degree of spikiness α ∈ [0, 1]. The resulting loss is bounded accordingly to α and the r−th latent eigenvalue is asymptotically strictly positive under the assumption δ ≤ α.

Idiosyncrasy proportion
In this way we encompass both LOREC and POET theory in a generalized theory of large covariance matrix estimation by low rank plus sparse decomposition.
A real example on a set of Euro Area banking data shows that our tool is particularly useful for mapping the covariance structure among variables even in a large dimensional context.

Supplementary material
The present paper is complemented by an Appendix containing a simulation study and the proofs of stated theorems and corollaries. In addition, the MATLAB functions UNALCE.m and POET.m, performing UNALCE and POET procedures respectively, can be downloaded at Farné and Montanari (2018). Both functions contain the detailed explanation of input and output arguments. Finally, the MATLAB dataset supervisory_data.mat, which contains the covariance matrix, C, and the relative labels of supervisory indicators, Labgood, can also be downloaded at the same link, which we refer to for the details.

A A simulation study A.1 Simulation settings
In order to compare the performance of UNALCE, LOREC and POET, we take into consideration five simulated low rank plus sparse settings reported in Table 3. The key simulation parameters are: 1. the dimension p, the sample size n; 2. the rank r and the condition number c of the low rank component L * ; 3. the trace of L * , τ θp, where τ is a magnitude parameter and θ is the percentage of variance explained by L * ; 4. the (half) number non-zeros s in the sparse component S * ; 5. the proportion of nonzeros prop s ; 6. the proportion of (absolute) residual covariance ρ corr .
The reported settings give an exhaustive idea of the low rank plus sparse settings recoverable under our assumptions. The critical parameters are: 1. the spectral norm of L * , which controls for the degree of spikiness. ||L * || is a direct function of τ and an inverse function of c, which together control for the magnitude of λ r (L * ); 2. the spectral norm of S * , which controls for the degree of sparsity. ||S * || is a direct function of prop s , which control for s and ρ corr . All the norms relative to our simulated settings are reported in Table 4. The data generation algorithm is described in detail in Farné (2016).
In Table 5 we summarize the features of our settings. Settings 1, 2 and 3 vary according to the degree of spikiness and sparsity. Setting 4 and 5 are intermediately spiked and sparse, and vary according to the ratio p/n. In particular, Settings 1,2 and 3 have p/n = 0.1, while Setting 4 has p/n = 1 and Setting 5 p/n = 2. The described features are pointed out in Figures 6 and 7, which show the degree of spikiness and sparsity across settings.
The reported scheme allows to achieve a convergence speed proportional to O(t 2 ). We define t * the number of steps needed for convergence. We set Y pre = Y t * −1 − 1 2 (Y t * −1 + Z t * −1 − Σ n ) and Z pre = Z t * −1 − 1 2 (Y t * −1 +Z t * −1 −Σ n ). The computational cost of the solution algorithm is proportional to p 4 √ ε , where ε is the required precision, while POET has the cost of a full-SVD (proportional to p 3 ). For more details see Farné (2016).
Lots of quantities are computed in order to describe comparatively the performance of the three methods on the same data. We call the low rank estimateL, the sparse estimateŜ, and the covariance matrix estimateΣ =Ŝ +L. The error norms used are: The correct classification rates of (true) non-zeros and zero elements (denoted respectively by sens and spec) are derived, as well as the correct classification rates of positive and negative elements separately considered (denoted respectively by senspos and specpos).

A.2 Simulation results
We start analyzing the performance ofΣ U N ALCE in comparison to the one ofΣ LOREC on our reference setting (Setting 1). In Figure 8 and 9 we report the differences between the Sample Total Losses and the Total Losses of LOREC and UNALCE for a grid of 20 × 20 = 400 threshold pairs. We note that the gain is positive everywhere, with the exception of the threshold pairs which do not return the exact rank (because they do not satisfy the range of Theorem 2.1). This pattern is more remarkable for Sample Total Loss than for Total Loss.
For both losses and each ψ, we note that, as explained, the gain across ρ never overcomes its maximum √ rψ (plotted for each ψ).
In Figure 10 we report the plot of the estimated proportion of latent variance θ across thresholds forΣ U N ALCE (in solid line the true θ = 0.7). In Figure 11 the same plot is reported forΣ LOREC . The shape is exactly the same as forΣ U N ALCE , the only difference is that all patterns are negatively shifted. In particular,θ gets closer to θ forΣ U N ALCE respect tô Σ LOREC in correspondence to all threshold combinations. shows the superior performance of UNALCE concerning the proportion of latent variance, of        Table 7, which contains the same results across different ratios p/n (Settings 2,4,5). Nevertheless, we note that the gap progressively decreases as p/n increases, due to the increased consistency with POET assumptions. We note, for instance, that the proportion of residual covariance is underestimated also by UNALCE for p/n = 2. At the same time, the performance of the proportion of residual variance detected by POET is upper biased, due to the natural bias of sample eigenvalues, and the bias decreases as the degree of spikiness increases. Tables 8 and 9 contain the error measures about the detection of the residual pattern across the degree of spikiness and the ratio p/n respectively. We note that POET, due to the lack of algebraic consistency, is completely unable to classify positive and negative elements.
On the contrary, UNALCE shows a recovery rate around 70% when p/n is small, while the detection capability deteriorates as p/n increases.
Tables 10 and 11 report the Euclidean distance between the vectors of estimated and true eigenvalues (denoted by eig), the condition number of the estimates and the estimated spectral norms. This Table can be compared to Table 4  To sum up, our UNALCE estimator outperforms POET concerning fitting and conditioning properties, detection of sparsity pattern, and eigen-structure recovery. We note that POET does not detect positive and negative elements at all. This is because it only has parametric consistency and not also the algebraic one. In addition, in order to obtain a positive definite estimate, cross validation selects a very high threshold for POET, and this causes the sparse estimate to be almost completely diagonal if p/n is large. On the contrary, the mathematical optimization procedure of UNALCE gets closer to the target and ensures to catch the algebraic spaces behind the two components.
under our assumption setting. Lemma (B.1) in turn relies on the following Lemmas by Fan et al. (2011): f it s jt ≤ C log(p) √ n .
As a consequence, we need to explore how Lemmas B.2, B.3, B.4 change under our assumptions.
The first step of the proof consists in decomposing E n = Σ n − Σ * in its four components: where: where f i and ǫ i are respectively the vectors of factor scores and residuals for each observation.
In order to show how Lemma B.3 changes, we need to recall some key results of Bickel and Levina (2008a). Differently from Luo's approach, in that setting (as in ours and in the POET one) the sparsity assumption is imposed to S * , and not to Σ * . While in Fan et al. (2013) the assumption max i≤p j≤p |s * ij | q = o(p), q ∈ [0, 1] is needed to ensure POET consistency, here Assumption 2.4 prescribes max j≤p ½(s * ij = 0) ≤ δ 2 p δ .
Consider now the uniformity class of sparse matrices in Bickel and Levina (2008a) Under Assumption 2.4 this class is no longer appropriate, because we can no longer write (see Bickel and Levina (2008a), page 2580) since the quantity c 0 (p) no longer scales to p but to p δ . Therefore, we need to replace c 0 (p) by C(p δ−1 ) in the proof which derives the rate of the sample covariance matrixŜ n under class (41) (see Bickel and Levina (2008a), page 2582), thus proving under Assumptions 2.3 and 2.4 that: Using Lemma B.5, we can derive because ||D 2 || ≤ p||D 2 || ∞ . Since log(p) ≪ n, we can write To conclude, we analyze Lemma B.4: Exploiting Assumption 2.5 we obtain r n ≤ C(p −3δ ). Therefore, the bound above becomes C p −2δ log p n .
In other words, the bound (44) means The proof relies on the combined use of proof tools by Fan et al. (2013), Fan et al. (2011), Fan et al. (2008 and Bickel and Levina (2008a).
Exploiting the basic property ||.|| ∞ ≤ ||.|| 2 and the minimum for γ in the range of Theorem 2.1, we can simply write According to Chandrasekaran et al. (2012) and Luo (2011a), the only probabilistic component of the error norm g γ (Ŝ − S * ,L − L * ) is g γ (E n ). Therefore, following the proof of Luo (2011a) and setting ψ = 1 ξ(T ) p α √ n , we can finally prove, under all the assumptions and conditions of Theorem 2.1, the thesis

B.2 Proof of Corollary 2.1
The proof directly descends by bound (44), because p α √ n → 0 if and only if p 2α n = o(1) as min(p, n) → ∞. As expected, the absolute bound vanishes only in the small dimensional case (n ≫ p α log(p)).

B.3 Proof of Corollary 2.2
Defined Σ r as the covariance matrix formed by the first r principal components, we know by dual Lidskii inequality that We start studying the behaviour of λ r (Σ n ). From bound (44) it descends the following Lemma Lemma B.6 Letλ r be the r−th largest eigenvalue of Σ n . If α ≥ δ, thenλ r > C 1 p α √ n with probability approaching 1 for some C 1 > 0.
Finally, assuming α ≥ δ, λ p−r+1 (L U N ALCE − Σ r ) vanishes asymptotically dividing by p α because both the relative errors ofL U N ALCE andL P OET = Σ r vanish (once assumed that r is known a priori or consistently estimated by UNALCE).
Then we can write By Cauchy-Schwartz inequality, it can be shown that Hence min L∈B(r),S∈Â(ŝ),Σ=L+S The problem in L is solved taking out the firstr principal components of Y pre . By The problem in S, assuming that the diagonal ofΣ ALCE is given and the off-diagonal elements ofŜ are invariant, leads to: The question now becomes: which diagonal elements of L ensure the minimum of B ′ + B ′′ ? Term B ′ is fixed respect to L, because we are assuming the invariance of diagonal elements inΣ (diag(Σ U N ALCE ) = diag(Σ ALCE )). The minimization of term B ′′ , given that rank(L) =r, falls back into the previous case, i.e. B ′′ is minimum if and only ifL = Optimality holds over the cartesian product of the set of all positive semi-definite matrices with rank smaller or equal to r,B(r), and the set of all sparse matrices with the same sparsity pattern asŜ ALCE such that diag(S) = diag(Σ ALCE − L), L ∈B(r) (we call this setÂ diag (ŝ)).
Consequently, we can write:

B.5 Proof of Corollary 2.3
We know that ||L U N ALCE −L ALCE || 2 =ψ. We can prove that because Y pre is uniquely determined by the conditioning event. The same inequality holds in spectral norm.
Since it holds given the conditioning event. As a consequence, since ||L U N ALCE −L ALCE || F ro = tr(L U N ALCE −L ALCE ) 2 = rψ 2 , we obtain The analogous triangular inequality for the sparse component is In order to quantify ||Ŝ U N ALCE −Ŝ ALCE || 2 F ro , we need to study the behaviour of the term p i=1 (l U N ALCE,ii −l ALCE,ii ) 2 , which is less or equal to rψ 2 , because it is less or equal to tr(L U N ALCE −L ALCE ) 2 . ||S − Z pre || 2 F ro + ||Z pre − S * || 2 F ro , because Z pre is uniquely determined by the conditioning event.
Given the matrix of eigenvectorsÛ , we can write diag( where K ii is a null matrix except for the i-th diagonal element equal toψ andû ′ i is the i-th row ofÛ . Similarly we can write of f −diag(L U N ALCE −L ALCE ) = p i=1 j =iû ′ iû j K ij where K ij is a null matrix except for the element ij equal toψ. Note that the rows ofÛ , differently from the columns, are not orthogonal.
Therefore, the real gain is terms of approximation of Σ n respect to ALCE measured in squared Frobenius norm is strictly positive and bounded from rψ 2 .
The quality of the estimation input ||Σ n − Σ * || does not depend on the estimation method.
Analogously, it is easy to prove that B.8 Proof of Corollary 2.5 We recall the following expression: We prove in sequence the three claims of the Corollary.
1. We start noting thatL U N ALCE ,L ALCE and U ALCEψ I r U ′ ALCE are r-ranked. We denote the respective spectral decompositions by: As a consequence, by Lidskii dual inequality we note that Therefore we obtain which proves the claim onŜ U N ALCE .