Optimization problem with matrix positivity constraints

Question

Is it possible for Mathematica to parametrically find a local minimum of a given function (say of some variable $x$ where $(a,b,c)$ serve as parameters) subject to the constraint that a given symmetric matrix is positive semi-definite? The entries of the matrix are monomials of $(x,a,b,c)$. Namely, the function is: \begin{equation} f(q)=\log\bigg[\frac{(1-q)(1-r^2t^2)(1-x)}{1-q-t^2+qr^2t^2-x+2t\sqrt{qx}-r^2(1-x-t^2+xt^2-xt^4+2t^3\sqrt{qx})}\bigg] \end{equation} (where the logarithm is taken to the base of 2)

       f=Log2[((-1+q)(1-r^2*t^2)(1-x))/(-1+q+t^2-qr^2*t^2+x-2*t*Sqrt[qx]-r^2(-1+t^2+x-t^2*x+t^4*x-2*t^3*Sqrt[q*x]))]

and the matrix is: \begin{equation} S=\left( \begin{array}{cccc} \frac{1}{r^2} & 1 & 1 & t\sqrt{x}\\ 1 & \frac{1}{t^2} & 1 & \frac{\sqrt{x}}{t}\\ 1 & 1 & 1 & \sqrt{q} \\ t\sqrt{x} & \frac{\sqrt{x}}{t} & \sqrt{q} & 1 \end{array} \right). \end{equation}

S={{1/r^2, 1, 1, t Sqrt[x]}, {1, 1/t^2, 1, Sqrt[x]/t}, 
   {1, 1, 1, Sqrt[q]}, {t Sqrt[x], Sqrt[x]/t, Sqrt[q], 1}};

the parameters $(x,r,t)$ are arbitrary value in $(0,1]$, which is also the domain of the variable $q$. Thus I need to find the minimum of $f(q)$ subject to $S\succeq0$ (in the sense of positive semi-definiteness) and $q\geq0$. I have no clue how to do it. I've tried using the FindMimimumcommand but didn't manage to make it work.
If is it indeed possible, it would very much help me to know how to formulate such a command. Thanks!

Answer 1

There are some places in the code above that you have neglected to make qr into q*r and similarly with qx. Addressing that, one might proceed as below.

ff[q_?NumberQ, x_?NumberQ, r_?NumberQ, 
  t_?NumberQ] := ((-1 + q) (1 - r^2*t^2) (1 - x))/(-1 + q + t^2 - 
    q*r^2*t^2 + x - 2*t*Sqrt[q*x] - 
    r^2 (-1 + t^2 + x - t^2*x + t^4*x - 2*t^3*Sqrt[q*x]))

detsmatrix[q_?NumberQ, x_?NumberQ, r_?NumberQ, t_?NumberQ] := 
 Det[{{1/r^2, 1, 1, t Sqrt[x]}, {1, 1/t^2, 1, Sqrt[x]/t}, {1, 1, 1, 
    Sqrt[q]}, {t Sqrt[x], Sqrt[x]/t, Sqrt[q], 1}}]

At this point I tried various things. Results were not great. I will remark that I had to restrict the variables to stay slightly away from the 0, 1 boundaries.

I attempted to enforce the semidefiniteness constraint by having the determinant of the matrix nonnegative. While this is not in general a sufficient condition, it often works if one begins inside the region where it is positive definite holds, as any change in sign of an eigenvalue would make the determinant negative.

{min, vals} = 
 FindMinimum[{Log[
    ff[q, x, r, t]], {.01 <= q <= .99, .01 <= x <= .99, .01 <= 
     r <= .99, .01 <= t <= .99, detsmatrix[q, x, r, t] >= 0}}, {q, x, 
   r, t}]

During evaluation of In[397]:= FindMinimum::eit: The algorithm does not converge to the tolerance of 4.806217383937354`*^-6 in 500 iterations. The best estimated solution, with feasibility residual, KKT residual, or complementary residual of {0.00149150538821,0.101362614558,0.0000278753724068}, is returned. >>

(* Out[397]= {0.0000743601741057, {q -> 0.956875379004, 
  x -> 0.98962743842, r -> 0.0333929679494, t -> 0.973088803016}} *)

One might also try adding a small multiple of the log of the determinant of the "s" matrix. This is, I believe, how interior point methods for LMI (semidefinite cone) programming proceed.

In[398]:= {min, vals} = 
 FindMinimum[{Log[ff[q, x, r, t]] + 
    Log[detsmatrix[q, x, r, t]]/10000, {.01 <= q <= .99, .01 <= 
     x <= .99, .01 <= r <= .99, .01 <= t <= .99, 
    detsmatrix[q, x, r, t] >= 0}}, {q, x, r, t}]

During evaluation of In[398]:= FindMinimum::eit: The algorithm does not converge to the tolerance of 4.806217383937354`*^-6 in 500 iterations. The best estimated solution, with feasibility residual, KKT residual, or complementary residual of {0.000510944299476,0.0386687828749,1.60849373579*10^-6}, is returned. >>

(* Out[398]= {0.0000881135619057, {q -> 0.921804018018, 
  x -> 0.989538131275, r -> 0.032681009984, t -> 0.955085543776}} *)

One can check that all eigenvalues of the s matrix are nonnegative at these approximated optima.

I think a problem is that, as one approaches certain vertices of the unit hypercibe in {q,x,r,t}, an eigenvalue goes to infinity and at least one other goes to zero. Possibly there are worse issues, such as one or more becoming negative. The upshot being, all I can do is outline plausible approaches. You'll need to play around with this to get results you think are usable.

--- edit ---

(1) I did a replacement incorrectly. The function ff should be defined as below (I had a q*r^2 that was intended to be (q*r)^2).

ff[q_?NumberQ, x_?NumberQ, r_?NumberQ, t_?NumberQ] := ((-1 + q) (1 - r^2*t^2) (1 - x))/(-1 + q + t^2 - (q*r)^2*t^2 + x - 2*t*Sqrt[q*x] - r^2 (-1 + t^2 + x - t^2*x + t^4*x - 2*t^3*Sqrt[q*x]))

(2) What I mean by solving for the derivative of the exponentiated objective being zero is achieved as follows. Compute the appropriate derivative, clear denominators, separate radicals (can do this since all factors inside radicals are assumed nonnegative).

f = ((-1 + q) (1 - r^2*t^2) (1 - x))/(-1 + q + 
     t^2 - (q*r)^2*t^2 + x - 2*t*Sqrt[q*x] - 
     r^2 (-1 + t^2 + x - t^2*x + t^4*x - 2*t^3*Sqrt[q*x]));

numer = Numerator[Together[D[f, q]]] /. (a_*b_)^c_ :> a^c*b^c

(* Out[30]= (-1 + r^2 t^2) (-1 + x) (Sqrt[q] r^2 Sqrt[x] + 
   Sqrt[q] t^2 Sqrt[x] - Sqrt[q] r^2 t^2 Sqrt[x] - 
   2 q^(3/2) r^2 t^2 Sqrt[x] + q^(5/2) r^2 t^2 Sqrt[x] - t x - q t x +
    r^2 t^3 x + q r^2 t^3 x + Sqrt[q] x^(3/2) - Sqrt[q] r^2 x^(3/2) + 
   Sqrt[q] r^2 t^2 x^(3/2) - Sqrt[q] r^2 t^4 x^(3/2)) *)

Now set to zero and solve for q in terms of the others.

In[28]:= sol = q /. Solve[numer == 0, q]

Out[28]= {Root[-t^2 x + 2 r^2 t^4 x - r^4 t^6 x + (r^4 + 2 r^2 t^2 - 2 r^4 t^2 + t^4 - 2 r^2 t^4 + r^4 t^4 + 2 r^2 x - 2 r^4 x - 4 r^2 t^2 x + 4 r^4 t^2 x + 6 r^2 t^4 x - 4 r^4 t^4 x - 2 r^2 t^6 x + x^2 - 2 r^2 x^2 + r^4 x^2 + 2 r^2 t^2 x^2 - 2 r^4 t^2 x^2 - 2 r^2 t^4 x^2 + 3 r^4 t^4 x^2 - 2 r^4 t^6 x^2 + r^4 t^8 x^2) #1 + (-4 r^4 t^2 - 4 r^2 t^4 + 4 r^4 t^4 - t^2 x - 4 r^2 t^2 x + 4 r^4 t^2 x + 2 r^2 t^4 x - 4 r^4 t^4 x + 3 r^4 t^6 x) #1^2 + (2 r^4 t^2 + 2 r^2 t^4 + 2 r^4 t^4 + 2 r^2 t^2 x - 2 r^4 t^2 x + 2 r^4 t^4 x - 2 r^4 t^6 x) #1^3 - 4 r^4 t^4 #1^4 + r^4 t^4 #1^5 &, 1], Root[-t^2 x + 2 r^2 t^4 x - r^4 t^6 x + (r^4 + 2 r^2 t^2 - 2 r^4 t^2 + t^4 - 2 r^2 t^4 + In[28]:= sol = q /.

Solve[numer == 0, q]

(* Out[28]= {Root[-t^2 x + 2 r^2 t^4 x - 
    r^4 t^6 x + (r^4 + 2 r^2 t^2 - 2 r^4 t^2 + t^4 - 2 r^2 t^4 + 
       r^4 t^4 + 2 r^2 x - 2 r^4 x - 4 r^2 t^2 x + 4 r^4 t^2 x + 
       6 r^2 t^4 x - 4 r^4 t^4 x - 2 r^2 t^6 x + x^2 - 2 r^2 x^2 + 
       r^4 x^2 + 2 r^2 t^2 x^2 - 2 r^4 t^2 x^2 - 2 r^2 t^4 x^2 + 
       3 r^4 t^4 x^2 - 2 r^4 t^6 x^2 + 
       r^4 t^8 x^2) #1 + (-4 r^4 t^2 - 4 r^2 t^4 + 4 r^4 t^4 - 
       t^2 x - 4 r^2 t^2 x + 4 r^4 t^2 x + 2 r^2 t^4 x - 
       4 r^4 t^4 x + 3 r^4 t^6 x) #1^2 + (2 r^4 t^2 + 2 r^2 t^4 + 
       2 r^4 t^4 + 2 r^2 t^2 x - 2 r^4 t^2 x + 2 r^4 t^4 x - 
       2 r^4 t^6 x) #1^3 - 4 r^4 t^4 #1^4 + r^4 t^4 #1^5 &, 1], 
 Root[-t^2 x + 2 r^2 t^4 x - 
    r^4 t^6 x + (r^4 + 2 r^2 t^2 - 2 r^4 t^2 + t^4 - 2 r^2 t^4 + 
       r^4 t^4 + 2 r^2 x - 2 r^4 x - 4 r^2 t^2 x + 4 r^4 t^2 x + 
       6 r^2 t^4 x - 4 r^4 t^4 x - 2 r^2 t^6 x + x^2 - 2 r^2 x^2 + 
       r^4 x^2 + 2 r^2 t^2 x^2 - 2 r^4 t^2 x^2 - 2 r^2 t^4 x^2 + 
       3 r^4 t^4 x^2 - 2 r^4 t^6 x^2 + 
       r^4 t^8 x^2) #1 + (-4 r^4 t^2 - 4 r^2 t^4 + 4 r^4 t^4 - 
       t^2 x - 4 r^2 t^2 x + 4 r^4 t^2 x + 2 r^2 t^4 x - 
       4 r^4 t^4 x + 3 r^4 t^6 x) #1^2 + (2 r^4 t^2 + 2 r^2 t^4 + 
       2 r^4 t^4 + 2 r^2 t^2 x - 2 r^4 t^2 x + 2 r^4 t^4 x - 
       2 r^4 t^6 x) #1^3 - 4 r^4 t^4 #1^4 + r^4 t^4 #1^5 &, 2], 
 Root[-t^2 x + 2 r^2 t^4 x - 
    r^4 t^6 x + (r^4 + 2 r^2 t^2 - 2 r^4 t^2 + t^4 - 2 r^2 t^4 + 
       r^4 t^4 + 2 r^2 x - 2 r^4 x - 4 r^2 t^2 x + 4 r^4 t^2 x + 
       6 r^2 t^4 x - 4 r^4 t^4 x - 2 r^2 t^6 x + x^2 - 2 r^2 x^2 + 
       r^4 x^2 + 2 r^2 t^2 x^2 - 2 r^4 t^2 x^2 - 2 r^2 t^4 x^2 + 
       3 r^4 t^4 x^2 - 2 r^4 t^6 x^2 + 
       r^4 t^8 x^2) #1 + (-4 r^4 t^2 - 4 r^2 t^4 + 4 r^4 t^4 - 
       t^2 x - 4 r^2 t^2 x + 4 r^4 t^2 x + 2 r^2 t^4 x - 
       4 r^4 t^4 x + 3 r^4 t^6 x) #1^2 + (2 r^4 t^2 + 2 r^2 t^4 + 
       2 r^4 t^4 + 2 r^2 t^2 x - 2 r^4 t^2 x + 2 r^4 t^4 x - 
       2 r^4 t^6 x) #1^3 - 4 r^4 t^4 #1^4 + r^4 t^4 #1^5 &, 3], 
 Root[-t^2 x + 2 r^2 t^4 x - 
    r^4 t^6 x + (r^4 + 2 r^2 t^2 - 2 r^4 t^2 + t^4 - 2 r^2 t^4 + 
       r^4 t^4 + 2 r^2 x - 2 r^4 x - 4 r^2 t^2 x + 4 r^4 t^2 x + 
       6 r^2 t^4 x - 4 r^4 t^4 x - 2 r^2 t^6 x + x^2 - 2 r^2 x^2 + 
       r^4 x^2 + 2 r^2 t^2 x^2 - 2 r^4 t^2 x^2 - 2 r^2 t^4 x^2 + 
       3 r^4 t^4 x^2 - 2 r^4 t^6 x^2 + 
       r^4 t^8 x^2) #1 + (-4 r^4 t^2 - 4 r^2 t^4 + 4 r^4 t^4 - 
       t^2 x - 4 r^2 t^2 x + 4 r^4 t^2 x + 2 r^2 t^4 x - 
       4 r^4 t^4 x + 3 r^4 t^6 x) #1^2 + (2 r^4 t^2 + 2 r^2 t^4 + 
       2 r^4 t^4 + 2 r^2 t^2 x - 2 r^4 t^2 x + 2 r^4 t^4 x - 
       2 r^4 t^6 x) #1^3 - 4 r^4 t^4 #1^4 + r^4 t^4 #1^5 &, 4], 
 Root[-t^2 x + 2 r^2 t^4 x - 
    r^4 t^6 x + (r^4 + 2 r^2 t^2 - 2 r^4 t^2 + t^4 - 2 r^2 t^4 + 
       r^4 t^4 + 2 r^2 x - 2 r^4 x - 4 r^2 t^2 x + 4 r^4 t^2 x + 
       6 r^2 t^4 x - 4 r^4 t^4 x - 2 r^2 t^6 x + x^2 - 2 r^2 x^2 + 
       r^4 x^2 + 2 r^2 t^2 x^2 - 2 r^4 t^2 x^2 - 2 r^2 t^4 x^2 + 
       3 r^4 t^4 x^2 - 2 r^4 t^6 x^2 + 
       r^4 t^8 x^2) #1 + (-4 r^4 t^2 - 4 r^2 t^4 + 4 r^4 t^4 - 
       t^2 x - 4 r^2 t^2 x + 4 r^4 t^2 x + 2 r^2 t^4 x - 
       4 r^4 t^4 x + 3 r^4 t^6 x) #1^2 + (2 r^4 t^2 + 2 r^2 t^4 + 
       2 r^4 t^4 + 2 r^2 t^2 x - 2 r^4 t^2 x + 2 r^4 t^4 x - 
       2 r^4 t^6 x) #1^3 - 4 r^4 t^4 #1^4 + r^4 t^4 #1^5 &, 5]}
       r^4 t^4 + 2 r^2 x - 2 r^4 x - 4 r^2 t^2 x + 4 r^4 t^2 x + 
       6 r^2 t^4 x - 4 r^4 t^4 x - 2 r^2 t^6 x + x^2 - 2 r^2 x^2 + 
       r^4 x^2 + 2 r^2 t^2 x^2 - 2 r^4 t^2 x^2 - 2 r^2 t^4 x^2 + 
       3 r^4 t^4 x^2 - 2 r^4 t^6 x^2 + 
       r^4 t^8 x^2) #1 + (-4 r^4 t^2 - 4 r^2 t^4 + 4 r^4 t^4 - 
       t^2 x - 4 r^2 t^2 x + 4 r^4 t^2 x + 2 r^2 t^4 x - 
       4 r^4 t^4 x + 3 r^4 t^6 x) #1^2 + (2 r^4 t^2 + 2 r^2 t^4 + 
       2 r^4 t^4 + 2 r^2 t^2 x - 2 r^4 t^2 x + 2 r^4 t^4 x - 
       2 r^4 t^6 x) #1^3 - 4 r^4 t^4 #1^4 + r^4 t^4 #1^5 &, 2], 
 Root[-t^2 x + 2 r^2 t^4 x - 
    r^4 t^6 x + (r^4 + 2 r^2 t^2 - 2 r^4 t^2 + t^4 - 2 r^2 t^4 + 
       r^4 t^4 + 2 r^2 x - 2 r^4 x - 4 r^2 t^2 x + 4 r^4 t^2 x + 
       6 r^2 t^4 x - 4 r^4 t^4 x - 2 r^2 t^6 x + x^2 - 2 r^2 x^2 + 
       r^4 x^2 + 2 r^2 t^2 x^2 - 2 r^4 t^2 x^2 - 2 r^2 t^4 x^2 + 
       3 r^4 t^4 x^2 - 2 r^4 t^6 x^2 + 
       r^4 t^8 x^2) #1 + (-4 r^4 t^2 - 4 r^2 t^4 + 4 r^4 t^4 - 
       t^2 x - 4 r^2 t^2 x + 4 r^4 t^2 x + 2 r^2 t^4 x - 
       4 r^4 t^4 x + 3 r^4 t^6 x) #1^2 + (2 r^4 t^2 + 2 r^2 t^4 + 
       2 r^4 t^4 + 2 r^2 t^2 x - 2 r^4 t^2 x + 2 r^4 t^4 x - 
       2 r^4 t^6 x) #1^3 - 4 r^4 t^4 #1^4 + r^4 t^4 #1^5 &, 3], 
 Root[-t^2 x + 2 r^2 t^4 x - 
    r^4 t^6 x + (r^4 + 2 r^2 t^2 - 2 r^4 t^2 + t^4 - 2 r^2 t^4 + 
       r^4 t^4 + 2 r^2 x - 2 r^4 x - 4 r^2 t^2 x + 4 r^4 t^2 x + 
       6 r^2 t^4 x - 4 r^4 t^4 x - 2 r^2 t^6 x + x^2 - 2 r^2 x^2 + 
       r^4 x^2 + 2 r^2 t^2 x^2 - 2 r^4 t^2 x^2 - 2 r^2 t^4 x^2 + 
       3 r^4 t^4 x^2 - 2 r^4 t^6 x^2 + 
       r^4 t^8 x^2) #1 + (-4 r^4 t^2 - 4 r^2 t^4 + 4 r^4 t^4 - 
       t^2 x - 4 r^2 t^2 x + 4 r^4 t^2 x + 2 r^2 t^4 x - 
       4 r^4 t^4 x + 3 r^4 t^6 x) #1^2 + (2 r^4 t^2 + 2 r^2 t^4 + 
       2 r^4 t^4 + 2 r^2 t^2 x - 2 r^4 t^2 x + 2 r^4 t^4 x - 
       2 r^4 t^6 x) #1^3 - 4 r^4 t^4 #1^4 + r^4 t^4 #1^5 &, 4], 
 Root[-t^2 x + 2 r^2 t^4 x - 
    r^4 t^6 x + (r^4 + 2 r^2 t^2 - 2 r^4 t^2 + t^4 - 2 r^2 t^4 + 
       r^4 t^4 + 2 r^2 x - 2 r^4 x - 4 r^2 t^2 x + 4 r^4 t^2 x + 
       6 r^2 t^4 x - 4 r^4 t^4 x - 2 r^2 t^6 x + x^2 - 2 r^2 x^2 + 
       r^4 x^2 + 2 r^2 t^2 x^2 - 2 r^4 t^2 x^2 - 2 r^2 t^4 x^2 + 
       3 r^4 t^4 x^2 - 2 r^4 t^6 x^2 + 
       r^4 t^8 x^2) #1 + (-4 r^4 t^2 - 4 r^2 t^4 + 4 r^4 t^4 - 
       t^2 x - 4 r^2 t^2 x + 4 r^4 t^2 x + 2 r^2 t^4 x - 
       4 r^4 t^4 x + 3 r^4 t^6 x) #1^2 + (2 r^4 t^2 + 2 r^2 t^4 + 
       2 r^4 t^4 + 2 r^2 t^2 x - 2 r^4 t^2 x + 2 r^4 t^4 x - 
       2 r^4 t^6 x) #1^3 - 4 r^4 t^4 #1^4 + r^4 t^4 #1^5 &, 5]} *)

So for specific values of the parameters {t,r,x} there are five points to check for possible minima. Then check boundary points of the positive cone. As remarked earlier, a superset of these will be where Det[matrix] = 0.

--- end edit ---

Answer 2

Seems like the problem is not solved completely in the previous post. For this problem, it is straightforward to obtain the positive semi-definite condition of S by Cholesky decomposition, which all the diagonal elements should be nonnegative.

Diagonal[CholeskyDecomposition[S]]

After some simplification, it can be shown the diagonal elements are

{1/r, Sqrt[-r^2 + 1/t^2], Sqrt[((1 - r^2) (1 - t^2))/(1 - r^2 t^2)], Sqrt[1 - ((1 -r^2 t^2) (Sqrt[q] - t Sqrt[x])^2)/((1 - r^2) (1 - t^2)) - x]}

With the given conditions, the constraint of positive semi-definiteness is

(1 - r^2) (1 - t^2) (1 - x) > (1 - r^2 t^2) (Sqrt[q] - t Sqrt[x])^2

As an example, at q=0.5, the minimization point is:

f = ((q - 1) (1 - r^2 t^2) (1 - x))/(-1 + q + t^2 - (q r)^2 t^2 + x - 2 t Sqrt[q x] - r^2 (-1 + t^2 + x - t^2 x + t^4 x - 2 t^3*Sqrt[q x]));
sol = Last@FindMinimum[{f, (1 - r^2) (1 - t^2) (1 - x) > (1 - r^2 t^2) (Sqrt[q] - t Sqrt[x])^2, 0 < x <= 1, 0 < r <= 1, 0 < t <= 1} /. q -> 0.5, Thread[{{r, x, t}, 0.5}]]

This naive way would not give the minimization point as a function of q directly since the minimization algorithm would fail if the initial point is not chosen properly. A possible way to avoid the issue is by numerical continuation. Just to provide the very basic idea behind it. One could compute the minimization point at a particular value q=q1 by starting at a known initial minimization point at q=q0 then gradually varying the parameter and at the same time performing the minimization successively until the desired parameter value q=q1 is reached.