Repair non-Positive Definite Correlation Matrix

Question

stephen 2011-4-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/6057-repair-non-positive-definite-correlation-matrix

编辑： Walter Roberson 2017-7-19

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

John D'Errico 2011-4-22

3
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/6057-repair-non-positive-definite-correlation-matrix#answer_8413

在 MATLAB Online 中打开

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.

6 个评论
显示 4更早的评论隐藏 4更早的评论

John D'Errico 2011-4-26

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 2016-6-1

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 2011-4-26

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/6057-repair-non-positive-definite-correlation-matrix#answer_8623

在 MATLAB Online 中打开

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]

3 个评论
显示 1更早的评论隐藏 1更早的评论

Katerina Rippi 2015-6-22

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 2017-7-19

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 2012-4-10

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/6057-repair-non-positive-definite-correlation-matrix#answer_43926

编辑：Walter Roberson 2017-7-19

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.