Generate constraints that ensure positive definiteness

Question

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`}} *)

```

Answer 1

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] &]

U = PadLeft@Internal`PartitionRagged[Array[u,d(d+1)/2], Range[d,1,-1]];
A = Transpose[U].U;

f[B_] := NMaximize[{
  norm[B] - penalty*Boole[Not[normtest[A] && Max[Eigenvalues[A] < 1]]],
  Splice[Thread[Diagonal[U] > 0]]}, 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]

(** 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 **)