Repair non-Positive Definite Correlation Matrix

Hi, I have a correlation matrix that is not positive definite. Does anyone know how to convert it into a positive definite one with minimal impact on the original matrix?

[1.0000  0.7426  0.1601  -0.7000  0.5500;
0.7426  1.0000  -0.2133  -0.5818  0.5000;
0.1601  -0.2133  1.0000  -0.1121  0.1000;
-0.7000  -0.5818  -0.1121  1.0000  0.4500;
0.5500  0.5000  0.1000  0.4500  1.0000]

Thanks for your help.

Stephen

Answer 1

[John D'Errico]

Your matrix is not that terribly close to being positive definite. As you can see, the negative eigenvalue is relatively large in context.

C =
            1       0.7426       0.1601         -0.7         0.55
       0.7426            1      -0.2133      -0.5818          0.5
       0.1601      -0.2133            1      -0.1121          0.1
         -0.7      -0.5818      -0.1121            1         0.45
         0.55          0.5          0.1         0.45            1

>> [V,D] = eig(C)
V =
       0.4365     -0.63792     -0.14229     -0.02851      0.61763
      0.29085      0.70108      0.28578    -0.064675      0.58141
      0.10029      0.31383     -0.94338     0.012435      0.03649
      0.62481      0.02315     0.048747     -0.64529     -0.43622
     -0.56958    -0.050216    -0.075752     -0.76056      0.29812
D =
     -0.18807            0            0            0            0
            0       0.1738            0            0            0
            0            0       1.1026            0            0
            0            0            0       1.4433            0
            0            0            0            0       2.4684

We can choose what should be a reasonable rank 1 update to C that will make it positive definite.

>> V1 = V(:,1);
>> C2 = C + V1*V1'*(eps(D(1,1))-D(1,1))
C2 =
       1.0358      0.76648      0.16833     -0.64871      0.50324
      0.76648       1.0159     -0.20781     -0.54762      0.46884
      0.16833     -0.20781       1.0019     -0.10031     0.089257
     -0.64871     -0.54762     -0.10031       1.0734      0.38307
      0.50324      0.46884     0.089257      0.38307        1.061

As you can see, it is now numerically positive semi-definite. Any more of a perturbation in that direction, and it would truly be positive definite.

>> eig(C2)
ans =
   2.5276e-16
       0.1738
       1.1026
       1.4433
       2.4684

Edit: The above comments apply to a covariance matrix. For a correlation matrix, the best solution is to return to the actual data from which the matrix was built. Then I would use an svd to make the data minimally non-singular.

Instead, your problem is strongly non-positive definite. It does not result from singular data. I would solve this by returning the solution I originally posted into one with unit diagonals.

>> s = diag(1./sqrt(diag(C2)))
s =
      0.98255            0            0            0            0
            0      0.99214            0            0            0
            0            0      0.99906            0            0
            0            0            0      0.96519            0
            0            0            0            0      0.97082
>> C2 = s*C2*s
C2 =
            1      0.74718      0.16524      -0.6152      0.48003
      0.74718            1     -0.20599     -0.52441      0.45159
      0.16524     -0.20599            1    -0.096732     0.086571
      -0.6152     -0.52441    -0.096732            1      0.35895
      0.48003      0.45159     0.086571      0.35895            1

As you can see, this matrix now has unit diagonals.

[stephen]

Thanks for the solution, John.

My concern though is the new correlation matrix does not appear to be valid, as the numbers in the main diagonal are now all above 1. Is there any way to create a new correlation matrix that is positive and definite but also valid? Thanks!

[Mads Svendsen-Tune]

I have problem similar to this one. I eventually just took absolute values of all eigenvalues. Could you comment a bit on why you do it this way and maybe on if my method makes any sense at all?

The problem with having a very small eigenvalue is that when the matrix is inverted some components become very large. Hence, standard errors become very large.

[John D'Errico]

Taking the absolute values of the eigenvalues is NOT going to yield a minimal perturbation of any sort.

[John D'Errico]

Stephen - true, I forgot that you were asking for a correlation matrix, not a covariance matrix. I've reformulated the solution.

[John D'Errico]

Mads - Simply taking the absolute values is a ridiculous thing to do. If you are computing standard errors from a covariance matrix that is numerically singular, this effectively pretends that the standard error is small, when in fact, those errors are indeed infinitely large!!!!!!!You are cooking the books. Why not simply define the error bars to be of width 1e-16? Wow, a nearly perfect fit!

[Hany Ferdinando]

John, my covariance matrix also has very small eigen values and due to rounding they turned to negative. I implemented you code above but the eigen values were still the same.

Answer 2

[Teja Muppirala]

I guess it really depends on what you mean by "minimal impact" to the original matrix.

Idea 1: This is probably not optimal in any sense, but it's very easy. Shift the eigenvalues up and then renormalize.

A0 = [1.0000  0.7426  0.1601  -0.7000  0.5500;
0.7426  1.0000  -0.2133  -0.5818  0.5000;
0.1601  -0.2133  1.0000  -0.1121  0.1000;
-0.7000  -0.5818  -0.1121  1.0000  0.4500;
0.5500  0.5000  0.1000  0.4500  1.0000];

k = min(eig(A0));

A = A0 - k*eye(size(A0));
Anew = A/A(1,1)

Idea 2: Treat it as a optimization problem. This code uses FMINCON to find a minimal perturbation (by percentage) that yields a matrix that has all ones on the diagonal, all elements between [-1 1], and no negative eigenvalues.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Anew = fixcorrmatrix

Abad = [1.0000  0.7426  0.1601  -0.7000  0.5500;
        0.7426  1.0000  -0.2133  -0.5818  0.5000;
        0.1601  -0.2133  1.0000  -0.1121  0.1000;
        -0.7000  -0.5818  -0.1121  1.0000  0.4500;
        0.5500  0.5000  0.1000  0.4500  1.0000];

x0 = zeros(10,1);

M = zeros(5);
indices = find(triu(ones(5),1));

opts = optimset('Display','iter','maxfunevals',1e4,'TolCon',1e-14,'algorithm','active-set');
x = fmincon(@(x) objfun(x,Abad,indices,M), x0,[],[],[],[],-2,2,...
    @(x) confun(x,Abad,indices,M),opts);

M(indices) = x;
Anew = Abad + M + M';

function E = objfun(x,Abad,indices,M)

M(indices) = x;
Anew = Abad + M + M';

% Set your error condition here
ERR = abs((Anew - Abad)./Abad);
E = max(ERR(:));

function [c,ceq] = confun(x,Abad,indices,M)
ceq = [];

M(indices) = x;
Anew = Abad + M + M';

% Positive definite and every element is between -1 and 1
c = [-min(eig(Anew));
    -1 + max(abs(Anew(:)))];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This yields the following answer:

     [1.0000    0.8345    0.1798   -0.6133    0.4819
      0.8345    1.0000   -0.1869   -0.5098    0.4381
      0.1798   -0.1869    1.0000   -0.0984    0.0876
     -0.6133   -0.5098   -0.0984    1.0000    0.3943
      0.4819    0.4381    0.0876    0.3943    1.0000]

[stephen]

Thanks!

If I knew part of the correlation is positive definite, e.g. the following correlation is positive definite

[1.0000 0.7426 0.1601 -0.7000;

0.7426 1.0000 -0.2133 -0.5818;

0.1601 -0.2133 1.0000 -0.1121;

-0.7000 -0.5818 -0.1121 1.0000]

It is when I added the fifth variable the correlation matrix became non-positive definite. Could I just fix the correlations with the fifth variable while keeping other correlations intact?

Thanks!

[Katerina Rippi]

Idea 2 also worked in my case! However, in case that we have more than 5 parameters, for example 6 arrows and columns then we say:

x0 = zeros(12,1);

M = zeros(6); indices = find(triu(ones(6),1));

or we also change something else?

[C Sung]

I'm trying to use this same Idea 2, but on a 48x48 correlation matrix. What do I need to edit in the initial script to have it run for my size matrix?

Answer 3

[Oi]

Dear Teja,

Thanks for your code, it almost worked to me. I'm also working with a covariance matrix that needs to be positive definite (for factor analysis). Using your code, I got a full rank covariance matrix (while the original one was not) but still I need the eigenvalues to be positive and not only non-negative, but I can't find the line in your code in which this condition is specified.

I'm totally new to optimization problems, so I would really appreciate any tip on that issue.