Is this a correct implementation of the principal-leading-minors test for positive-semidefiniteness?

Question

I am asking this question in response to comments by mikado and Daniel Lichtblau on my question Maximize a six-dimensional function subject to joint positive-semidefiniteness constraints

I gave two matrices there

{{1/3 - a2/2, -((I a1)/2), (I a3)/2}, {(I a1)/2, 1/3 + a2/2, 0}, {-((I a3)/2), 0, 1/3}}

and

{{1/4, 0, b1/2, 0}, {0, 1/4, 1/2 (I b2 - b3), 0},
 {b1/2, 1/2 (-I b2 - b3), 1/4, 0}, {0, 0, 0, 1/4}}

Let's call the first $C1$ and the second $C2$. I want to ensure that these two ($3 \times 3$ and $4 \times 4$ "density") matrices are positive-semidefinite.

The code I used for the "principal leading minors" test in order to implement it was

T = Array[1, 3];
Do[T[[k]] = FullSimplify[ComplexExpand[FullSimplify[Det[Take[C1, {1, k}, {1, k}]]]]{k, 1,3}];
constraint1 = T[[1]] >= 0;
Do[constraint1 = constraint1 && T[[i]] >= 0, {i, 2, 3}];
constraint1 = FullSimplify[constraint1]

giving

9 (a1^2 + a2^2) <= 4 && 18 (a1^2 + a2^2) + 9 (2 + 3 a2) a3^2 <= 8

and

T = Array[1, 4];
Do[T[[k]] = FullSimplify[ComplexExpand[FullSimplify[Det[Take[C2, {1, k}, {1, k}]]]]{k, 1,4}];
constraint2 = T[[1]] >= 0;
Do[constraint2 = constraint2 && T[[i]] >= 0, {i, 2, 4}];
constraint2 = FullSimplify[constraint2]

giving

4 (b1^2 + b2^2 + b3^2) <= 1

The question posed in Maximize a six-dimensional function subject to joint positive-semidefiniteness constraints

was to maximize

Abs[a1 b1] + Abs[a2 b2] + Abs[a3 b3]

subject to the intersection

9 (a1^2 + a2^2) <= 4 && 18 (a1^2 + a2^2) + 9 (2 + 3 a2) a3^2 <= 8 && 4 (b1^2 + b2^2 + b3^2) <= 1

of the two constraints.

All kosher?

(Note: In the indicated prior question, $C2$ had a "typo" of 01/2--as pointed out by mikado--rather than 1/2--I guess due to my unartful cutting-and-pasting.)

Incidentally, "density matrices" are "self-adjoint (or Hermitian), positive semi-definite, of trace one".

Answer 1

An efficient way to check for positive (semi)definiteness is to try to compute the Cholesky decomposition of a Hermitian matrix. To simplify things, here's a modified version of the routine from this answer, which compute the root-free variant for a Hermitian matrix:

LDLH[mat_?SquareMatrixQ] := Module[{n = Length[mat], mt = mat, v, w},
     Do[If[j > 1,
           w = mt[[j, ;; j - 1]]; v = Conjugate[w] Take[Diagonal[mt], j - 1];
           mt[[j, j]] -= w.v;
           If[j < n,
              mt[[j + 1 ;;, j]] -= mt[[j + 1 ;;, ;; j - 1]].v]];
        mt[[j + 1 ;;, j]] /= mt[[j, j]], {j, n}];
     {LowerTriangularize[mt, -1] + IdentityMatrix[n], Diagonal[mt]}]

Then, here is how to generate the required positive semidefiniteness conditions from your matrices:

m1 = {{1/3 - a2/2, -I a1/2, I a3/2}, {I a1/2, 1/3 + a2/2, 0}, {-I a3/2, 0, 1/3}};

Reduce[And @@ Thread[Simplify[ComplexExpand[Last[LDLH[m1]]]] >= 0],
       {a1, a2, a3}, Reals]
   -(2/3) < a1 < 2/3 &&
   -(1/3) Sqrt[4 - 9 a1^2] < a2 < 1/3 Sqrt[4 - 9 a1^2] &&
   -Sqrt[((8 - 18 a1^2 - 18 a2^2)/(18 + 27 a2))] <= a3 <=
   Sqrt[(8 - 18 a1^2 - 18 a2^2)/(18 + 27 a2)]

For the other matrix:

m2 = {{1/4, 0, b1/2, 0}, {0, 1/4, 1/2 (I b2 - b3), 0},
      {b1/2, 1/2 (-I b2 - b3), 1/4, 0}, {0, 0, 0, 1/4}};

And @@ Thread[Simplify[ComplexExpand[Last[LDLH[m2]]]] >= 0]
   1/4 - b1^2 - b2^2 - b3^2 >= 0