Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Does Mathematica have a way to "fix" a correlation matrix that is not positive semi-definite?

I looked through the documentation and search the internet but could not find anything.

share|improve this question

3 Answers 3

up vote 12 down vote accepted

Here is simple (unweighted) Mma version of the Matlab implementation of Covariance Bending.

ClearAll[covBending];
covBending[mat_, tol_: 1/10000]:=If[PositiveDefiniteMatrixQ[mat], mat, 
NestWhile[(Eigensystem[#][[2]].DiagonalMatrix[
  Max[#, tol] & /@(Eigensystem[#][[1]])].Transpose[
  Eigensystem[#][[2]]]) &, N@mat,  Min[Eigensystem[#][[1]]] < tol &]]

examples:

 matrices = ((1/2) (# + Transpose[#]) & /@  RandomReal[{0, 1}, {10, 4, 4}]);
{MatrixForm[#], PositiveDefiniteMatrixQ[#], MatrixForm[covBending[#]],
 PositiveDefiniteMatrixQ[covBending[#]]} & /@  matrices // TableForm

enter image description here

Update: the previous method continues the iteration until all eigenvalues are at or above the threshold, that is, it does not check and stop if the current matrix is PD. The following is a version that tests for positive definiteness using PositiveDefiniteMatrixQ at every iteration:

ClearAll[covBending2];
covBending2[mat_, tol_: 1/10000] := If[PositiveDefiniteMatrixQ[mat], mat,
   NestWhile[ (Eigensystem[#][[2]].DiagonalMatrix[
     Max[#, tol] & /@ (Eigensystem[#][[1]])].Transpose[
     Eigensystem[#][[2]]]) &,  N@mat, !PositiveDefiniteMatrixQ[#]&]]
share|improve this answer
    
Why not have the test be ! PositiveDefiniteMatrixQ[mat]? –  J. M. Aug 31 '12 at 7:59
    
@J.M. yes indeed, checking if mat is positive definite is actually much faster than checking if minimum eigenvalue of mat less than a threshold. –  kguler Aug 31 '12 at 8:15
    
Does this ensure that the diagonal has only ones ? –  Faysal Aberkane Aug 31 '12 at 12:20
    
@Faysal, no. This is for a generic square matrix. For a covariance matrix symmetry is preserved, I think, unless you use a an arbitrary weight vector/matrix. For correlation matrices you need to add renormalization. –  kguler Aug 31 '12 at 12:29
    
@kguler Maybe I'm wrong, but isn't this a fundamentally different problem than "fixing" non psd correlation/covariance matrices? This looks like it belongs to the class of shrinkage estimators and is converting a general square symmetric matrix to a psd matrix. One could argue that if a matrix is demonstrably far from psd (in some sense) other than due to minor numerical errors, then it was never a covariance matrix to begin with. Of course, OP was very vague about what they meant by "fixing" and I do recognize that yours could've very well been their implied meaning... –  rm -rf Aug 31 '12 at 21:04

Generally, the reason why matrices that were supposed to be positive semi-definite but are not, is because the constraint of working in a finite precision world often introduces a wee bit of perturbation in the lowest eigenvalues of the matrix, making it either negative or complex. These errors are generally of the order of machine precision, but is enough to return False if you use standard algorithms to check for positive semi-definiteness. This can be easily fixed with Chop. For example, something like the following:

psdMat = #1.Chop@#2.#3\[ConjugateTranspose] & @@ SingularValueDecomposition[mat];

Now if this isn't sufficient to make your matrix positive semi-definite, you should go back and take a closer look at your problem to see if there are other reasons to not expect it to be positive semi-definite.

share|improve this answer
4  
The last sentence cannot be emphasized enough; in fact, before applying any "correction" method, you should try to figure out why things are not coming out to be positive (semi)definite in the first place. "Correction" is more of a Band-Aid than anything... –  J. M. Aug 31 '12 at 8:22
    
When calculating the correlation of multiple entities, there can be missing values between entity 1 and entity 3, but not between 2 and 3. If there are many of these, you can produce a correlation matrix that is not positive semi-definite. The easy answer is to just use the data where it exists for all 3 entities, but if you want to use all the available data this can be a problem. Another example of this, is that you may have some of the correlation matrix be known but other values are not and the user needs to enter what they think is the correct correlation amount. –  Stephen Lien Aug 31 '12 at 13:53

I like the approach in this paper. Below is a starting point for an implementation.

MakeGoodCorrelationMatrix[badCorrMat_,flag_]:=
Which[
    flag==1,
        Module[{LocalN=Length[badCorrMat],LocalTheta,LocalB,LocalOutput,LocalSolution},
            LocalTheta=Outer[Unique[] &, Range[LocalN], Range[LocalN-1], 1, 1];
            LocalB = Outer[If[#2 < LocalN, Cos[LocalTheta[[#1, #2]]] (Times @@ Sin[LocalTheta[[#1, 1 ;; #2 - 1]]]),
                Times @@ Sin[LocalTheta[[#1, 1 ;; #2 - 1]]]] &, Range[LocalN], Range[LocalN], 1, 1] //Chop;
            LocalOutput = Dot[LocalB, Transpose[LocalB]] //Chop;
            LocalSolution = Minimize[Total[Flatten[(LocalOutput - badCorrMat)^2]], Flatten[LocalTheta]]//Chop;
            LocalOutput /. LocalSolution[[2]] 

            ],  
    True,
        Module[{LocalLambdap = (Max[0.000001, #] & /@ Eigenvalues[badCorrMat]) , LocalS = Eigenvectors[badCorrMat] // Transpose, LocalT, LocalB, LocalOutput}, 
            LocalT = DiagonalMatrix[Dot[#^2,LocalLambdap]^(-1) & /@ LocalS] ; 
            LocalB = Dot[Sqrt[LocalT],LocalS,Sqrt[DiagonalMatrix[LocalLambdap]]]; 

            LocalOutput = Dot[LocalB,Transpose[LocalB]] ;
            LocalOutput = 1/2 (LocalOutput + Transpose[LocalOutput]) (* just to make sure the output is symmetric *)

            ]
]

Example :

badMat = {{1, 0.05, 0.9}, {0.05, 1, 0.8}, {0.9, 0.8, 1}};
Eigenvalues[badMat]
(* {2.22926, 0.950346, -0.179603} *)

goodMat1 = MakeGoodCorrelationMatrix[badMat, 1];
Eigenvalues[goodMat1]
(* {2.11798, 0.882022, 8.32667*10^-17} *)

goodMat2 = MakeGoodCorrelationMatrix[badMat, 2];
Eigenvalues[goodMat2]
(* {2.08911, 0.910888, 9.35756*10^-7} *)
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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