Find minimum with matrix positive-definiteness constraint

Question

Say I wish to find the minimum value of the determinant of a symmetric matrix under the condition that the matrix be positive definite. So I attempt:

M = {{a,0},{0,b}}

FindMinimum[{Det[M],a>=1,b>=1,PositiveDefiniteMatrixQ[M]},{a,b}]

This returns an error that Constraints in {False} are not all equality or inequality constraints..., suggesting to me that the PositiveDefiniteMatrixQ is being evaluated immediately for arbitrary a,b and not evaluated each iteration for a,b values.

Then I might try to delay the evaluation of PositiveDefiniteMatrixQ with Delayed, which returns a similar error Constraints in {Delayed[PositiveDefiniteMatrixQ[M]],a>=1,b>=1} are not all equality or inequality constraints.

How can I impose such a constraint on the FindMinimum function?

---

The specific problem I am trying to minimise involves two coupled matrices, one 2x2 and one 4x4. Simplified as far as possible while still exhibiting an issue with the Thread[Eigenvalues[B] > 0] approach:

A = {{a, 0}, {0, d^2*b + a - 2 d*c*Sign[d]}};
B = {{a, 0, c, 0}, {0, a, 0, -c}, {c, 0, b, 0}, {0, -c, 0, b}};

min = FindMinimum[{
    Det[A],
    a^2 + b^2 - 2 c^2 >= 0,
    Thread[Eigenvalues[B] > 0],
    a >= 1, b >= 1, -1 < d < 1}, 
    {a, b, c, {d, 0}}]

Thread[Eigenvalues[B /. min[[2]]] > 0]

Answer 1

One simple way to do this, is to just introduce a function that returns a huge number when the constraint is not met:

ClearAll[det2];
det2[mat_?PositiveSemidefiniteMatrixQ] := Det[mat];
det2[_?(MatrixQ[#, NumericQ] &)] := 10^100;

FindMinimum[{det2[M], a >= 1, b >= 1}, {a, b}]
(* {1., {a -> 1., b -> 1.}} *)

This works for this toy example, but it's unlikely that this is going to give good results for realistic problems since it will be very difficult for the solver to figure out where exactly the constraint holds and where not and how to optimize the objective function around this constraint.

Instead, it's almost always better to try and find a way to formulate the problem in such a way that the constraint is guaranteed to hold. For example, any Hermitian, positive definite matrix can be written as a Cholesky decomposition L . Transpose[L] (with L lower triangular). So the general way to write a 2 x 2 positive definite matrix is:

M = With[{L = {{a, 0}, {b, c}}}, L. Transpose[L]]
(* {{a^2, a b}, {a b, b^2 + c^2}} *)

The eigenvalues of this matrix are >= 0:

Minimize[#, {a, b, c}] & /@ Eigenvalues[M]
(*{{0, {a -> -1, b -> -1, c -> 0}}, {0, {a -> 0, b -> 0, c -> 0}}}*)

So this matrix you could use for your minimization option without having to worry about the PositiveSemidefiniteMatrixQ constraint.

Answer 2

For a real symmetric matrix X, we can replace the constraint PositiveDefiniteMatrixQ[X] with Thread[Eigenvalues[X] > 0].

To get around the issue caused by Sign in OP's example A we can apply PiecewiseExpand on Det[A] in the objective function:

min = FindMinimum[{PiecewiseExpand[Det[A], Element[d, Reals]], 
   a^2 + b^2 - 2 c^2 >= 0, Thread[Eigenvalues[B] > 0], a >= 1, 
   b >= 1, -1 < d < 1}, {a, b, c, {d, 0}}]

{1., {a -> 1., b -> 1.98128, c -> 0.641898, d -> 0.}}

Thread[Eigenvalues[B /. min[[2]]] > 0]

{True, True, True, True}

Answer 3

Since the OP's matrix B is intended to be positive definite, let's take a look at what its Cholesky decomposition is supposed to look like (while incorporating all the other known constraints on the variables):

Simplify[CholeskyDecomposition[{{a, 0, c, 0}, {0, a, 0, -c}, {c, 0, b, 0}, {0, -c, 0, b}}],
         a >= 1 && b >= 1 && -1 < d < 1 && c ∈ Reals]
   {{Sqrt[a], 0, c/Sqrt[a], 0}, {0, Sqrt[a], 0, -(c/Sqrt[a])},
    {0, 0, Sqrt[b - c^2/a], 0}, {0, 0, 0, Sqrt[b - c^2/a]}}

Remembering that the diagonal elements of this matrix are intended to be real, we instantly pick out the conditions a >= 0 (which can be eliminated since we already have the a >= 1 constraint) and b - c^2/a >= 0. With that, we can easily formulate the required FindMinimum[] call:

FindMinimum[{Det[A], a^2 + b^2 - 2 c^2 >= 0 && a >= 1 && b >= 1 && -1 < d < 1 && b - c^2/a >= 0},
            {a, b, c, {d, 0}}, WorkingPrecision -> 20]
   {1.6799589323874212830*10^-8,
    {a -> 1.0385747471268498288, b -> 4.1403066994506190213,
     c -> 2.0736484553649652263, d -> 0.50084415969450801764}}

We can then verify that

{a^2 + b^2 - 2 c^2 >= 0 /. Last[%], PositiveDefiniteMatrixQ[B /. Last[%]]}
   {True, True}

Answer 4

Following on @kglr, but using the definition of positive definiteness as all of the leading principal minors are > 0. Create a function to compute the leading principal minors. Assumes square matrix.

makeLPM[mat_] := Table[Det@mat[[1 ;; i, 1 ;; i]], {i, 1, Length@mat}]

lpm = makeLPM[B];
(* {a, a^2, a^2 b - a c^2, a^2 b^2 - 2 a b c^2 + c^4} *)

Replicating @kglr's code, substituting principal minors for eigenvalues

min = FindMinimum[{PiecewiseExpand[Det[A], Element[d, Reals]], 
                  a^2 + b^2 - 2 c^2 >= 0, 
                  Thread[lpm > 0], a >= 1, 
                   b >= 1, -1 < d < 1}, 
                  {a, b, c, {d, 0}}]
 (* {1., {a -> 1., b -> 2.16112, c -> 0.941631, d -> 0.}} *)

This is a different answer from @kglr, but...

Thread[Eigenvalues[B /. min[[2]]] > 0]
 (* {True, True, True, True} *)

 Thread[makeLPM[B /. min[[2]]] > 0]
 (* {True, True, True, True} *)