10
$\begingroup$

Are there any known ways of getting an unbiased estimate of the condition number of the true covariance matrix being estimated, or at least correct within a small number of orders of magnitude? For concreteness, you can consider the 2-norm condition number; but the one-norm (or equivalently infinity norm for a symmetric matrix) would also be of interest, and might be easier to do. Or an (almost) unbiased estimate of any reasonable approximate condition number or reciprocal condition number in any norm, such as a reasonable estimate of MATLAB's RCOND (which estimates the reciprocal of the 1-norm condition number) of the true but unknown matrix?

For some perspective, the sample covariance matrix is singular, and therefore has an infinite condition number unless the number of samples exceeds the dimension of the matrix. If the true covariance matrix is not singular, the condition number of the sample covariance exhibits a rather extreme bias as an estimate of the condition number of the true covariance matrix when the number of samples does not exceed the dimension. Aware of the dangers of using an estimated covariance matrix with an excessively high condition number, LeDoit and Wolf wrote the serious, but humorously titled paper "Honey, I Shrunk the Sample Covariance Matrix But their proposed method essentially assumes you know what you want to shrink the covariance to, not dealing with the situation in which the true condition number might be very high.

If every element of a matrix is estimated with an unbiased estimate of the individual element, the condition number of the matrix consisting of these elements could be a severely biased estimate of the condition number of the true matrix being estimated. That is kind of the point of the LeDoit and Wold paper.

The 2-norm condition number of a real matrix is the ratio of the largest to smallest singular value, or equivalently if the matrix is symmetric, the ratio of the largest to smallest eigenvalue. So we can already see a ratio thing going on there, which usually is not wonderful for unbiased estimation.

Another manifestation is calculating Hessian (which are symmetric matrices) estimates, possibly via Quasi-Newton update rules, whereby errors in the gradients being differenced can eventually cause runaway condition number growth of the Hessian estimate. The trick in that situation is distinguishing whether the condition number growth reflects the condition number of the true Hessian (as for the famous Rosenbrock banana function), or is due to errors in the gradient cascading over the course of several updates resulting in runaway condition numbers. In the case of Hessian calculation and estimation, as an alternative to condition number, it would be useful to get a ballpark (within s mall number of orders of magnitude) estimate of the largest eigenvalue of the true matrix being estimated.

$\endgroup$
5
  • $\begingroup$ Just to be clear, I'm talking about estimating the condition number of the true matrix being estimated, not the condition number of the estimate of that matrix, which is a straightforward Linear Algebra calculation. $\endgroup$ Commented Nov 1, 2022 at 20:36
  • 3
    $\begingroup$ I appreciate from your comment that you want to estimate the condition number of the true covariance matrix $\text{cond}(\Sigma)$, not $\text{cond}(\hat\Sigma)$ the condition number of an estimator of $\Sigma$. Perhaps you could edit the title and body of your post to make that clear? As it stands, the title seems a tad misleading, since you don't want to estimate $\text{cond}(\hat\Sigma)$, but you use the term "estimated covariance matrix." $\endgroup$
    – Sycorax
    Commented Nov 1, 2022 at 23:03
  • $\begingroup$ Precondition what matrix? I don't know the true matrix, otherwise I could just calculate the condition number. I know how to constrain the condition number of an estimated covariance matrix. However, that doesn't tell me anything about what the true condition number of the true covariance matrix is. That's the intended point of my question. $\endgroup$ Commented Nov 5, 2022 at 2:40
  • $\begingroup$ (1) I feel it would be more fruitful to pursue an unbiased estimate of $log(cond(A)) = log(s_{max}) - log(s_{min})$. (2) I feel it would have something to do with the Marchenko-Pastur distribution. (3) The condition number is typically only a heuristic for indicating that “results are not going to be good”, so the motivation for an unbiased estimate is not clear to me — sorry if I’m offending the OP. $\endgroup$
    – Number
    Commented Nov 7, 2022 at 14:20
  • $\begingroup$ @Number Hey great. I'll take an (approximately) unbiased estimate of the log of condition number in any norm. I'm trying to get in the ballpark. Within a 2 or 3 orders of magnitude is good enough. Your 250 points await if you can deliver!! I am not offended, "at least ballpark" is in the question title, and "within a small number of orders of magnitude" is in the first sentence. $\endgroup$ Commented Nov 7, 2022 at 15:07

1 Answer 1

4
+250
$\begingroup$

I would like to preface my answer by clarifying that I fully understand that the question requests estimating the condition number of the true unknown matrix, and not the condition number of the estimate of the matrix which is available to us (from the available samples).

My answer does not directly address this requirement; instead, I propose 2 things:

  1. An answer to the last part of the question: an estimate of the largest eigenvalue of the true matrix being estimated. My answer uses the available matrix as a proxy for the unknown true matrix, and my assumption is that this estimate will apply to the true matrix.
  2. Links to 3 papers which directly deal with estimating large covariance matrices using shrinkage, including a recent paper by Ledoit and Wolf: Quadratic Shrinkage for Large Covariance Matrices (published 2022).

First, an answer for the last part of the question:

as an alternative to condition number, it would be useful to get a ballpark (within s mall number of orders of magnitude) estimate of the largest eigenvalue of the true matrix being estimated.

I would like to emphasize that my answer uses the available matrix as a proxy for the unknown true matrix, and my assumption is that this estimate will apply to the true matrix. A justification for this assumption (although not proof) has to do with the radii in the Gershgorin circle theorem, which are directly affected by the "sample variation" between the samples which are available to us, from which the available sample covariance matrix is created, and the true unknown covariance matrix. I'm contemplating a proof of this, but do not yet have one, meanwhile I'd like to present my unproven answer for now.

The Gershgorin circle theorem (https://en.wikipedia.org/wiki/Gershgorin_circle_theorem) can provide estimates of the locations of eigenvalues of a given square matrix. To summarize, this theorem says: the eigenvalues of any matrix are located within discs whose centers are the diagonal elements, and radii are row or column sums of absolute-values of the non-diagonal elements corresponding to each diagonal element (the minimum of the row or column sum can be used because the same theorem applies also to the transpose of the given matrix).

A nice property of this method is that it gives an $O(n^2)$ algorithm for an $n$-by-$n$ matrix, which is much faster than a general eigensolver $O(n^3)$, when the covariance matrix $C$ does not have any known structure to exploit.

In the current case of a covariance matrix $C$, it is symmetric, so the eigenvalues are all real and non-negative. Symmetry also means row sums are identical to column sums, of course.

Therefore, a reasonable estimate for the largest eigenvalue $\lambda_{max} = \lambda_1$ can be obtained as: $$\hat{\lambda_{max}} := \max_i (c_i + r_i)$$ where $c_i = C_{ii}$ are the diagonal elements, and $r_i = \sum_{j \neq i}{|C_{ij}|}$ are the radii corresponding to each diagonal element.

A different method for estimating the largest eigenvalue is using a few iterations of the Power iteration (https://en.wikipedia.org/wiki/Power_iteration) from a random initialization, where each iteration requires only a matrix-by-vector product, so is $O(n^2)$, but the number of iterations for a reasonable estimate is unknown in advance, since it is determined by the ratio of the unknown 2nd-largest to largest eigenvalues $\lambda_2 / \lambda_1$.

Instead of several iterations of the Power iteration, one could try several different random initializations, doing a single iteration in each, and picking the maximum estimated eigenvalue, with an appropriate scaling. This has been proposed at:

End of chapter 2 (p. 58) in Numerical Recipes in C, 2nd edition, by Press, Teukolsky, Vetterling, and Flannery.


Returning to the original problem of estimating the condition number.

It would be tempting to use the Gershgorin circle theorem to also estimate the smallest eigenvalue $\lambda_{min} = \lambda_n$ as: $$\hat{\lambda_{min}} := \max \{ \min_i (c_i - r_i), 0 \} \geq 0$$

However, in a few simple cases I tried, this is an extremely disappointing estimate.


Since the motivation of the OP was estimation of covariance matrices, not just their condition numbers, then additional material which may be useful are the following pair of papers, which deal with a covariance matrix having a so-called "spiked" model, which roughly means a signal subspace of low rank within a full-rank noise subspace.

Donoho, Gavish, and Johnstone: Optimal Shrinkage of Eigenvalues in the Spiked Covariance Model, https://arxiv.org/abs/1311.0851

Donoho and Ghorbani: Optimal Covariance Estimation for Condition Number Loss in the Spiked Model, https://arxiv.org/abs/1810.07403

The OP already mentioned the 2004 paper by Ledoit and Wolf, which is also referenced in the papers above. Olivier Ledoit's home-page (http://www.ledoit.net) presents updated research in this area, including the most recent paper which claims an improvement over Stein's shrinkage:

Ledoit and Wolf: Quadratic Shrinkage for Large Covariance Matrices, submitted 2021: http://www.ledoit.net/BEJ1911-021R1A0.pdf, published 2022: https://www.econstor.eu/bitstream/10419/228874/1/1743676301.pdf

$\endgroup$
3
  • $\begingroup$ The bulk of your text is way off the mark, because my interest is not in fast low accuracy computation of eigenvalues of sample covariance matrix, it is for eigenvalues of population covariance. Accepted answer and awarded bounty points based on the reference links, which are great. $\endgroup$ Commented Nov 11, 2022 at 1:44
  • $\begingroup$ @MarkL.Stone I strongly believe that the estimation of the largest population eigenvalue using the ordinary sample covariance matrix with the estimator I proposed, based on Gershgorin circles, is statistically superior to the exact computation of the largest eigenvalue from the the ordinary sample covariance matrix (using a standard library function); However, regretfully I cannot yet prove this. The computational advantage is just a nice by-product which I didn't want to ignore to mention. Nonetheless, I did explicitly admit my answer was not exactly what you require, sorry. Thanks $\endgroup$
    – Number
    Commented Nov 11, 2022 at 2:01
  • $\begingroup$ have you made any progress on proofs? $\endgroup$ Commented Dec 7, 2023 at 11:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.