# Correcting a correlation matrix to be positive semidefinite

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.

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[#]&]]


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


• Why not have the test be ! PositiveDefiniteMatrixQ[mat]? Aug 31, 2012 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.
– kglr
Aug 31, 2012 at 8:15
• Does this ensure that the diagonal has only ones ? Aug 31, 2012 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.
– kglr
Aug 31, 2012 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, 2012 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.

• 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... Aug 31, 2012 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. Aug 31, 2012 at 13:53

I like the approach in this paper: The most general methodology to create a valid correlation matrix for risk management and option pricing purposes. Below is a starting point for an implementation.

MakeGoodCorrelationMatrix[badCorrMat_,flag_]:=
Which[
flag==1,
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}};

• @Jagra I think the above works but it's far from being "production grade". The first approach (flag==1) was quite slow for matrices above 5x5 (maybe use NMinimize ?), the second approach shouldn't have a hard-coded cutoff and I did not check it will give reasonable results for large matrices. Do post your own version when ready please. Mar 6, 2018 at 8:28