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

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


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".

## 1 Answer

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

• Quite impressive! Clearly, the simple (spherical-type) results for the $4\times4$ case agree for the two analyses. As for the $3\times3$ case, I asked FindInstance unsuccessfully to give me an example where your constraint was satisfied and mine wasn't. As to an example, where my constraint was satisfied and yours wasn't I got {{a1 -> Sqrt[799]/48, a2 -> -(5/16), a3 -> 0}}, yielding the matrix {{a1 -> Sqrt[799]/48, a2 -> -(5/16), a3 -> 0}}, the apparent density matrix of a "pure state" (determinant 0). So, has your methodolgy failed in this situation? What's going on here? Is my method kosher? Feb 2, 2020 at 17:08
• If I simplify the conditions generated by LDLH[] from your matrix, I get {3 a2 <= 2, 2 + 3 a2 + (9 a1^2)/(-2 + 3 a2) >= 0, 2 + (9 (2 + 3 a2) a3^2)/(-4 + 9 a1^2 + 9 a2^2) >= 0} Look at the form of the last entry. If a3 = 0, then the last entry becomes 2, which is indeed nonnegative. However, the denominator also becomes 0 for those values of a1 and a2. Thus, you will need to be careful in taking these conditions apart. As an example: CylindricalDecomposition[And @@ Thread[FullSimplify[ComplexExpand[Last[LDLH[m1]]]] >= 0], {a3, a2, a1}] Feb 2, 2020 at 17:19
• Well, I'm employing these constraints in three- and six-dimensional integrations, so singular/lower-dimensional cases such as these should not play a role. Feb 2, 2020 at 17:42
• Show how this algorithm works for a matrix containing symbols
– dtn
Dec 7, 2021 at 17:51