1

For a given problem, I am interested in estimating a correlation matrix.

In this problem, I can somewhat easily get estimates of the pairwise correlations. Each of these estimates should be consistent for the true pairwise correlation. From this, I can theoretically construct an estimate of the correlation matrix that is consistent as well.

However, I don't think there's any guarantee that any finite-sized correlation matrix will be non-negative definite! This is similar to the issue that constructing a correlation matrix from pairwise estimates of the correlation with missing data can lead to a non-negative definite correlation matrix.

For various reasons, I would really like a non-negative definite matrix. Is there any established methods for doing so? My first guess would be to just multiply the off diagonals by $\eta$, where $\eta$ was the largest values such that $\alpha \hat C $ is non-negative definite ($\hat C$ is the naive correlation matrix estimated by filling in the off diagonals). Is there better ideas? If this idea is good, is there any justification for it?

Cliff AB
  • 17,741
  • 1
  • 39
  • 84

1 Answers1

2

This can be solved by (linear convex) semidefinite programming. Choose a loss function, such as the frobenius norm of the difference between the original matrix $M$, and the adjusted matrix $C$. Minimize this loss function subject to the constraints that $C$ is positive semi-definite and has all diagonal elements equal to $1$. More generally, a minimum eigenvalue, possibly positive, can be specified for $C$.

Here is a made-up example, showing the solution using minimum eigenvalue for $C$ of $0$, i.e., that $C$ be psd. CVX under MATLAB is used, but this can be done under R, Python, or many other languages.

Note: everything in a line after % is a comment, per MATLAB syntax

>> M = [1 .9 .8;.9 1 .99;.8 .99 1]  % input matrix
M =
   1.000000000000000   0.900000000000000   0.800000000000000
   0.900000000000000   1.000000000000000   0.990000000000000
   0.800000000000000   0.990000000000000   1.000000000000000

>> disp(eig(M)) % calculate and display eigenvalues of M
  -0.007591945123683
   0.212026312718646
   2.795565632405038
>> cvx_begin
>> variable C(3,3) semidefinite % constrains C to be symmetric psd
>> minimize(norm(C-M,'fro')) % objective function is frobenius norm of C - M
>> diag(C) == 1 % constrains all disgonal elements of C to = 1
>> cvx_end % numerically solves the optimization problem
>> disp(C) % display optimal value of C
   1.000000000000000   0.897254763937705   0.802029497299419
   0.897254763937705   1.000000000000000   0.983333732823410
   0.802029497299419   0.983333732823410   1.000000000000000
>> disp(eig(C))
   0.000000000927173
   0.209549435800457
   2.790450563272369

Note that the minimum eigenvalue of C is not exactly zero due to optimization solver tolerance.

Many other variants are possible, including imposition of various constraints and/or changing the objective function. For instance, certain elements of $C$ could be constrained to have the same values as the corresponding elements of $M$.

Mark L. Stone
  • 12,546
  • 1
  • 31
  • 51
  • Very nice! In this situation, we see that the solution obtained by minimizing the Frobenius norm appears to have minimal impact on the eigen values, and presumably the eigen vectors as well. As I think about it more, that is exactly what I want (closest in terms of eigen vectors and values). So another method I could use would be to look at the eigen values and vectors, and merely set all the negative eigen values to 0. – Cliff AB Jan 23 '17 at 00:24
  • 2
    Using your proposed procedure, the diagonal elements will (generally) not = 1. In the example in my answer, using your proposed procedure results in diagonal elements of the adjusted matrix being 1.000433368904208, 1.004627516774850, 1.002531059444626. But fi you're from the liberal wing of probability in which you can give 150% and have probabilities > 1, see http://stats.stackexchange.com/questions/4220/can-a-probability-distribution-value-exceeding-1-be-ok/160979#160979 , then I suppose self-correlations exceeding 1 would also be o.k. :) – Mark L. Stone Jan 23 '17 at 00:41
  • Good point. If I wanted to keep driving, I could use that method to derive a valid (but not necessarily full rank) covariance matrix, from which I could derive a valid correlation matrix. What I'm really hoping for in this problem is estimation of the eigen values/vectors, and it's not clear what method will perform better for this. – Cliff AB Jan 23 '17 at 00:50
  • The most critical item is likely to be the minimum eigenvalue specified for the adjusted matrix. Non-zero minimum eigenvalue, eigmin, can be reduced to the example shown in my answer by constraining $C$ - eigmin * Identity_matrix, rather than $C$, to be symmetric psd. – Mark L. Stone Jan 23 '17 at 00:53
  • As another example, if all you cared about was matching the spectral norm, you could minimize the 2-norm of the difference between the adjusted and unadjusted matrices, subject to psd and diagonal elements = 1 constraints. You could even add in a constraint on the maximum frobenius and/or 1-norm (which is same as infinity norm given symmetry). This is just as easy to solve as the frobenius norm minimization. Directly optimizing eigenvectors is likely to be a more difficult non-convex optimization problem. – Mark L. Stone Jan 23 '17 at 01:04