# Generate constraints that ensure positive definiteness

What is a good way to generate algebraic constraints that ensure matrix be positive definite? Ideally, I'd be able to do something like below

Solve[# \[Element] Reals & /@ Eigenvalues[A]]

However, this doesn't directly work. Practical example below uses this to find the norm of a positive linear operator (related issue). It works, but requires AposDefiniteConstraints to be specified manually which I'd like to avoid.

(also tried Thread[Eigenvalues[X] > 0] suggestion from Find minimum with matrix positive-definiteness constraint but I get Maximize returning unevaluated)

(* Find norm of a positive transformation of a positive definite \
d-by-d matrix *)
SeedRandom[1];
d = 2;
symmetricMatrix[d_] := Array[a[Min[#1, #2], Max[#1, #2]] &, {d, d}];
extractVars[mat_] := DeleteDuplicates@Cases[Flatten@A, _a];

(* using built-in Norm/Simplify too slow, use this helper instead *)

norm[A_] :=
Max[x /. # & /@ Solve[CharacteristicPolynomial[A, x] == 0, x]];

A = symmetricMatrix[d];
Avars = extractVars[A];

B = Mean[#\[Transpose].A.# & /@
Table[RandomReal[{-1, 1}, {d,
d}], {d^2}]]; (* random positive transformation of A *)
normA =
norm[A];
normB = norm[B];
AposDefiniteConstraints =
a[1, 1]^2 + 4 a[1, 2]^2 - 2 a[1, 1] a[2, 2] + a[2, 2]^2 >= 0 &&
a[1, 1]^2 + 4 a[1, 2]^2 - 2 a[1, 1] a[2, 2] + a[2, 2]^2 >= 0;
Maximize[{normB, normA < 1,
AposDefiniteConstraints}, Avars] (* => {0.7853700810760375,{a[1,1]\
\[Rule]0.999855037823971,a[1,2]\[Rule]0.00017274783320670866,a[2,2]\
\[Rule]0.9997941436806035}} *)

$$$$

• If M is Hermitian, then look at the trace inequalities (Tr[M]^2/Tr[MatrixPower[M, 2]] > n - 1 ) && Tr[M] > 0 – flinty Sep 9 '20 at 16:46
• Your norm looks redundant - it's just the same asMax[Eigenvalues[A]]. – flinty Sep 9 '20 at 16:54

Instead of using constraints, you could use a penalty in the objective. Whenever the constraints are violated it subtracts a large penalty with the hope of pushing NMaximize away from bad values:

(** Given random matrix X, find max eigenvalue of (Transpose[X].A.X)
where A is posdef and max eigenvalue of A is < 1 **)
penalty = 10^20;
d = 2;

(* this is a hack to shut up Max when complex numbers appear *)
norm[m_] := Max[If[Not[Element[#, Reals]],-penalty,#] & /@ Eigenvalues[m]]

normtest[m_] := AllTrue[Eigenvalues[m], Element[#, Reals]&]

(* refer to the trace inequalities *)
positivedef[m_] :=
Tr[m]^2/Tr[MatrixPower[m, 2]] > First[Dimensions[m]] - 1 && Tr[m] > 0

A = Array[a[Min[#1, #2], Max[#1, #2]] &, {d, d}];

f[B_] := NMaximize[
norm[B] - penalty*Boole[Not[positivedef[A]]] -
penalty *Boole[Not[normtest[A] && Max[Eigenvalues[A] < 1]]],
Variables[A], Method -> "RandomSearch"]

SeedRandom[1];
(* random positive transformation of A *)
b = Mean[Transpose[#].A.# & /@ Table[RandomReal[{-1, 1}, {d, d}], {d^2}]];

{maxn, asub} = f[b]
Eigenvalues[A /. asub]
PositiveDefiniteMatrixQ[A /. asub]

(** results:
{0.738925, {a[1, 1] -> 0.799585, a[1, 2] -> 0.176808, a[2, 2] -> 0.815972}}
{0.984776, 0.630781}
True **)


There are problems for d > 2 so we need another approach. One idea I had was to use the CholeskyDecomposition. If matrix $$A$$ is positive-definite and Hermitian, then it has a decomposition $$U^\top U$$ where $$U$$ is upper triangular, and real valued with a positive diagonal. We then need only find the entries $$u_i$$ of $$U$$ to form $$A$$ with the constraint that $$\mathrm{diag}(U)\succeq \mathbf{0}$$.

This eliminates the need for the first penalty, but there are issues with convergence for d > 2 and the result may not be close enough to optimal:

penalty = 10^20;
d = 3;

(*this is a hack to shut up Max when complex numbers appear*)
norm[m_] := Max[If[Not[Element[#,Reals]],-penalty,#]& /@ Eigenvalues[m]]

normtest[m_] := AllTrue[Eigenvalues[m], Element[#, Reals] &]

A = Transpose[U].U;

f[B_] := NMaximize[{
norm[B] - penalty*Boole[Not[normtest[A] && Max[Eigenvalues[A] < 1]]],
Method -> "RandomSearch"]

SeedRandom[1];
(*random positive transformation of A*)
b = Mean[Transpose[#].A.# & /@ Table[RandomReal[{-1,1}, {d,d}], {d^2}]];

{maxn, asub} = f[b]
Eigenvalues[A /. asub]
PositiveDefiniteMatrixQ[A /. asub]

(** NMaximize::cvmit: Failed to converge to the requested accuracy or precision within 100 iterations. **)

(** results:
{0.491483, {u[1] -> 0.159054, u[2] -> 0.619449, u[3] -> -0.0776527, u[4] ->
0.595631, u[5] -> 0.0898834, u[6] -> 0.588458}}
{0.751889, 0.360839, 0.0114554}
True **)

• Neat trick! Is this supposed to work for more dimensions? If I replace d=2 with d=3 in that example, it gives a solution with large negative eigenvalues – Yaroslav Bulatov Sep 10 '20 at 13:14
• Yeah it seems to perform badly for d > 2. I have a suspicion this is because of the positivedef. I'm looking at another way to generate positive definite matrices through a CholeskyDecomposition. – flinty Sep 10 '20 at 14:34
• ^ this approach gave me a result for d == 3, but it didn't go as high for d == 2. – flinty Sep 10 '20 at 15:19