3
$\begingroup$

There is this highly upvoted but unanswered post at Math.SE. I would like to try some ideas. This requires numerical experiments with Hermitian and positive definite matrices $A$, $B$, and $C$ such that $$ A+B+C=I. $$

Following this answer, I tried to generate such matrices with a translation of their Matlab algorithm:

n = 3; Do[ a = DiagonalMatrix[RandomReal[{0, 1}, n]]; R = MatrixPower[IdentityMatrix[3] - a, 1/2]; (* {U,S,V}=SingularValueDecomposition[ RandomReal[{0,1},{n,n}]+I*RandomReal[{0,1},{n,n}]]; *) U = RandomVariate[GaussianUnitaryMatrixDistribution[n]]; b = R . U . DiagonalMatrix[RandomReal[{0, 1}, n]] . ConjugateTranspose[U] . R; c = IdentityMatrix[3] - a - b; Print[PositiveDefiniteMatrixQ[#] & /@ {a, b, c}]; , {k, 10}]

Unfortunately, matrix $C$ is almost never PSD. How to fix the code? Additionally, any ideas how to speed up the calculation or steer it towards matrices with a given norm of deviation from $I$ are welcome.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
1
$\begingroup$

Here is a brute force method:

Find three matrices with positive eigenvalues that satisfy $A + B + C = I$.

a = DiagonalMatrix[RandomReal[{0, .1}, 3]] b = DiagonalMatrix[RandomReal[{0, .1}, 3]] c = IdentityMatrix[3] - a - b While[And @@ # < 0 & /@ Diagonal[c], a = DiagonalMatrix[RandomReal[{0, .1}, 3]]; b = DiagonalMatrix[RandomReal[{0, .1}, 3]]; c = IdentityMatrix[3] - a - b; ]

That will be slow for large n, but could be made faster by picking random diagonal entries sequentially.

A random rotation matrix and its inverse

rm = RotationMatrix[RandomReal[{0, 2 Pi}], RandomReal[{0, 1}, 3]] rmi = Inverse[rm]

Three Hermitian matrices:

{aH, bH, cH} = rmi . # . rm & /@ {a, b, c}

Check:

PositiveDefiniteMatrixQ /@ {aH, bH, cH} (*{True, True, True}*) Chop[aH + bH + cH] == N[IdentityMatrix[3]] (*True*)

post comment edit:

One can make the choices of the diagonal matrices more general and faster narrowing the range of random number selection:

sumtoOne[n_] := Module[{list = {RandomReal[]}}, Do[ AppendTo[list, RandomReal[{0, 1 - Total[list]}]], n - 2 ]; Append[list, 1 - Total[list]] ]

The constraint being that each entry in the diagonal must be less than 1 and the combined entries sum to 1. This can be used to generate positive definite Hermitian matrices, but I don't think they are sampling from all possible choices.

Then as before:

{{a, b, c} = Table[DiagonalMatrix[sumtoOne[3]], 3]

with a more general similarity transformation

mt = RandomComplex[{-10 - 10 I , 10 + 10 I}, {3, 3}] mti = Inverse[mt]

Less limited space of Hermitian matrices

{aH, bH, cH} = Chop[mti . # . mt & /@ {a, b, c}] Chop[Eigenvalues /@ {aH, bH, cH}] (*all positive*)
Improve this answer
$\endgroup$
2
  • $\begingroup$ Sorry, my previous comment was not relevant, but I still have a question. By generating eigenvalues in the [0, 1/10] interval, don't we miss some large space of possibilities? I have a feeling that by using the same transformation matrix, we essentially reduce everything to a scalar case. Isn't it? $\endgroup$ Commented May 10 at 15:14
  • $\begingroup$ Yes, It does limit the space. It sounded like you are after a method to create instances for numerical experiments? I've edited my response to make the space a bit bigger. One could make the space bigger yet by employing more general simularity transformations. $\endgroup$ Commented May 10 at 20:18
1
$\begingroup$

Here is an approach that (I hope) should work every time.

Define a random positive definite matrix (I don't guarantee that this is a good distribution - you can replace it with a better!).

rnd[n_] := Module[{m}, m = RandomVariate[NormalDistribution[0, 1], {n, n}]; Transpose[m] . m]

Now scale the matrices so that they remain positive definite, but have the correct sum. You might eliminate asymmetry due to rounding errors (average with half the transpose).

rnd3[n_] := Module[{a, b, c, L}, a = rnd[n]; b = rnd[n]; c = rnd[n]; L = Inverse[CholeskyDecomposition[a + b + c]]; Map[Transpose[L] . # . L &, {a, b, c}]]

Generate a test example

 {a, b, c} = rnd3[3];

Check the identity

 a + b + c - IdentityMatrix[3] // Chop (* {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}} *)

And positive definite

PositiveDefiniteMatrixQ /@ {a, b, c} (* {True, True, True} *)
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.