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.

I would like to complete a correlation matrix with missing entries, such that the resulting matrix is positive semi-definite. Say I want to find some x and some y in the next matrix that satisfy positive semi-definiteness:

matrix = {{1, y, 0.8}, {y, 1, x}, {0.8, x, 1}};

I have tried FindInstance, Maximize and FindMaximum. They all give good results. First a define a function that is used to feed the optimisation functions:

testmatrix[input_, xtest_, ytest_] :=  input /. {x -> xtest, y -> ytest}

The solutions now are:

sol1 = Maximize[Min[Eigenvalues[testmatrix[matrix, x, y]]], {x, y}]

sol2 = FindInstance[Min[Eigenvalues[testmatrix[matrix, x, y]]] > 0 && 0 < x < 1 && 0 < y < 1, {x, y}, Reals, 1]

sol3 = FindMaximum[Min[Eigenvalues[testmatrix[matrix, x, y]]], {x, y}]

All results give a positive semi-definite correlation matrix:

PositiveDefiniteMatrixQ[testmatrix[matrix, x /. sol1[[2]], y /. sol1[[2]]]]
PositiveDefiniteMatrixQ[testmatrix[matrix, x /. sol2[[2]], y /. sol2[[2]]]]
PositiveDefiniteMatrixQ[testmatrix[matrix, x /. sol3[[1]], y /. sol3[[1]]]]

However, when the size of the matrix and the number of unknowns increase, the calculation time explodes using all methods. Note that I don't require anything else than positive semi-definiteness. I would assume that FindInstance would be the fastest candidate since it does not require an optimal value. Who knows a way to make the procedure more efficient? Your help is much appreciated.

share|improve this question
1  
There are at least 2 fast methods to correct a correlation matrix (in case one of the eigenvalues is negative and so on); the correction will attempt to minimally change the entries and "fix" the matrix. Would such solution to acceptable ? If so, you could just fill the missing entries with random numbers or your best estimate or other criterion. –  b.gatessucks Jan 13 at 12:58
    
How big are the matrices you're trying to generate/fill-in? –  Ymareth Jan 13 at 13:25
    
I know about methods that fix a correlation matrix that is not positive semi-definite. The thing is that I don't want some correlation values to be adjusted whereas I don't care about others. The matrix that I want to fix in the end is of size 14x14. –  Frits Jan 14 at 10:41

3 Answers 3

Let we have a $n\times n$ matrix $M$ with unknowns $x_1,x_2,\dots,x_m$. The measure of the positive definiteness we define as

$$ f(x_1,x_2,\dots,x_m) = \sum_k\varepsilon_k\theta(-\varepsilon_k) $$

where $\varepsilon_k$ is $k$-th eigenvalue of the matrix $M$ and $\theta(x)$ is the Heaviside step function. For a positive definite matrix $f$ is equal to $0$, otherwise it is negative. We should maximize $f$ with respect to $x_1,x_2,\dots,x_m$.

n = 200;
M = #.Transpose[#] &@RandomReal[NormalDistribution[], {n, n}];

Here $M$ is the random positive definite $200\times 200$ matrix (so-called Wishart ensemble in the random matrix theory). Let it have 10% unknown values (4000!)

m = Round[n^2/10];
pos = RandomInteger[{1, n}, {m, 2}];
vars = Array[Symbol["x" <> ToString[#]] &, m];
mat[x_] := ReplacePart[M, Thread[{#, #[[{2, 1}]]} & /@ pos -> x]]

The measure is

f[x_ /; And @@ NumericQ /@ x] := Total[# UnitStep[-#]] &@Eigenvalues@mat[x];

The pattern x_ /; And @@ NumericQ /@ x is necessary to forbid a symbolic simplification (there is no chance to calculate eigenvalues symbolically). Also we can obtain the gradient of the measure from the perturbation theory (it will produce 25x speedup!)

$$ \frac{\partial f}{\partial x_l} = \sum_k(2-\delta_{i,j})(v_k)_i(v_k)_j\theta(-\varepsilon_k) $$

where $v_k$ is $k$-th eigenvector, $(i,j)$ is the position of $x_l$ in $M$.

g[x_ /; And @@ NumericQ /@ x] := Table[(2 - KroneckerDelta @@ pos[[i]]) Total[
       UnitStep[-#] #2[[All, pos[[i, 1]]]] #2[[All, pos[[i, 2]]]]], 
           {i, m}] & @@ Eigensystem@mat[x];

Let's go!

x = FindArgMax[f[#], Transpose@{#, 0. #}, Gradient -> g[#]] &@vars;

It takes about 4 seconds on my old laptop.

Verification:

f[x]

0.

Histogram@Eigenvalues@mat[x]

enter image description here

PositiveDefiniteMatrixQ@mat[x]

True

Update: I assumed that $M$ is a general positive definite matrix. The correlation matrix has the additional properties

$$ M_{ii} = 1, \quad -1\le M_{ij}\le 1. $$

However, it does not require any constraints in the above algorithm since for any positive definite matrix $|M_{ij}|\le M_{ii}$. So $x_1,x_2,\dots,x_m$ will be in the interval $[-1,1]$ automatically.

share|improve this answer
    
Thanks a lot for your extensive reply. I'm impressed! I see that the resulting matrix mat[x]//MatrixForm does not have ones at the diagonal. How should I adjust your code for that? And how should I enter the matrix to start with? –  Frits Jan 14 at 10:44
    
I've adjusted your code a bit and created a module that does exactly what I was aiming for. I will post it as an answer to my question. Thanks a lot for your help. It is really appreciated. –  Frits Jan 14 at 13:08
    
@Frits I read "correlation" as "covariance" and wrote a general algorithm (the covariance matrix is more general than the correlation matrix). See my notes in the update. –  ybeltukov Jan 14 at 16:19
up vote 1 down vote accepted

With the brilliant answer of ybeltukov below, I was able to create a module that solves exactly my problem.

So we have a matrix with unknown correlations:

matrix = {{1, y, 0.8}, {y, 1, z}, {0.8, z, 1}};

The method of ybeltukov implemented in a Module:

completeMatrix[inputMatrix_] := Module[{pos, nrUnknowns, vars, mat, f, g},  
pos = Position[matrix, _?(! NumericQ[#] &)];
pos = Cases[pos, {_, _}];
pos = Cases[pos, {_?(# > 0 &), _?(# > 0 &)}];
nrUnknowns = Length[pos];
vars = Array[Symbol["x" <> ToString[#]] &, nrUnknowns];
mat[x_] := ReplacePart[matrix, Thread[{#, #[[{2, 1}]]} & /@ pos -> x]];
f[x_ /; And @@ NumericQ /@ x] := Total[# UnitStep[-#]] &@Eigenvalues@mat[x];
g[x_ /; And @@ NumericQ /@ x] := Table[(2 - KroneckerDelta @@ pos[[i]]) Total[UnitStep[-#] #2[[All, pos[[i, 1]]]] #2[[All, pos[[i, 2]]]]], {i, nrUnknowns}] & @@Eigensystem@mat[x];
corrs = FindArgMax[f[#], Transpose@{#, 0. #}, Gradient -> g[#]] &@vars;
mat[corrs]
]

The resulting matrix is found in a split second:

completeMatrix[matrix] // MatrixForm
PositiveDefiniteMatrixQ[%]
share|improve this answer
    
pos = Position[UpperTriangularize@matrix, _Symbol, Heads -> False] is a bit more convenient. –  ybeltukov Jan 14 at 16:25
    
ah, yes that's much nicer. I knew it could have been more elegant. thanks! –  Frits Jan 15 at 16:07

This could be done by semidefinite programming: The condition that the $(i,j)$ entry of $X$ has a given value $b$ is of the form $\langle E_{ij}, X \rangle = b$ where $E_{ij}$ has a $1$ in position $(i,j)$ and a zero everywhere else; you want to find an $X$ which solves various conditions like this and is positive semi-definite. This matches the description in , the Wikipedia article I linked except that the set up there also has a linear functional $\langle C,X \rangle$ which you are trying to optimize. Since you don't care how the missing entries are chosen, you can choose $C$ any way you like. (I might try minimizing the sum of the entries of $X$.)

Unfortunately, Mathematica doesn't seem to have native handling for semidefinite programming; see this question.

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.