10
$\begingroup$

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.

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Jan 13, 2014 at 12:58
  • $\begingroup$ How big are the matrices you're trying to generate/fill-in? $\endgroup$ Commented Jan 13, 2014 at 13:25
  • $\begingroup$ 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. $\endgroup$ Commented Jan 14, 2014 at 10:41

3 Answers 3

Reset to default
6
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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? $\endgroup$ Commented Jan 14, 2014 at 10:44
  • $\begingroup$ 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. $\endgroup$ Commented Jan 14, 2014 at 13:08
  • $\begingroup$ @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. $\endgroup$ Commented Jan 14, 2014 at 16:19
1
$\begingroup$

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[%]
Improve this answer
$\endgroup$
2
  • $\begingroup$ pos = Position[UpperTriangularize@matrix, _Symbol, Heads -> False] is a bit more convenient. $\endgroup$ Commented Jan 14, 2014 at 16:25
  • $\begingroup$ ah, yes that's much nicer. I knew it could have been more elegant. thanks! $\endgroup$ Commented Jan 15, 2014 at 16:07
0
$\begingroup$

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.

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.