# Find minimum with matrix positive-definiteness constraint

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,
a >= 1, b >= 1, -1 < d < 1},
{a, b, c, {d, 0}}]


• if M is symmetric, you can use FindMinimum[{Det[M], a >= 1, b >= 1, Thread[Eigenvalues[M] > 0]}, {a, b}]?
– kglr
Feb 20, 2020 at 16:33
• For your given matrix M, the eigenvalues are $a$ and $b$ and the Det[M] is $a b$. Since you are constraining $a$ and $b$ to be >=1, the minimum is clearly at $a=b=1$. Feb 20, 2020 at 16:42
• @kglr It is symmetric, so I can use that approach. However, for the found minimum values Thread[...] returns {True,True,False,False} (for a 4x4 matrix) so the constraint is not being respected Feb 20, 2020 at 17:01
• @bills Yes, but that's a simplified case for demonstration. My actual matrix has a more complex (but symmetric) form and a longer list of constraints. Feb 20, 2020 at 17:02
• Cailean, can you update your question with (1) the information that M is symmetric and (2) the 4X4 example you mention in the comment?
– kglr
Feb 20, 2020 at 17:09

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.

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[]] > 0]


{True, True, True, True}

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}


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[]] > 0]
(* {True, True, True, True} *)


• In v13.0, there are some errors in your code min.@MikeY Mar 28 at 3:35
• Interesting, the problem is in the Thread[ ] command, which works fined outside of the min[ ] call now, and worked fined inside it before. Will update answer. Mar 28 at 22:40
• I can't get it to work, short of pasting the results of Thread[lpm>0] directly into the code. Tried every trick I could think of or find. Bug? Mar 28 at 23:43