2
$\begingroup$

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

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

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 **)
$\endgroup$
3
  • $\begingroup$ 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 $\endgroup$ Sep 10 '20 at 13:14
  • $\begingroup$ 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. $\endgroup$
    – flinty
    Sep 10 '20 at 14:34
  • $\begingroup$ ^ this approach gave me a result for d == 3, but it didn't go as high for d == 2. $\endgroup$
    – flinty
    Sep 10 '20 at 15:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.