# NMinimize keeps returning Indeterminate for well defined function

I am attempting to numerically optimize a function of only 2 parameters k, \[Theta].

The function is well-defined and the constraints are simple. However, the optimisation keeps returning Indeterminate or not a real number comment over the parameter space, i.e. function value not a real number at {k,\[Theta]} = {0.918621,0.716689}.

If I instead initialize {k,[Theta]} to these values and run the function, it does provide a value. I am unsure why the optimization fails. The MWE is provided below. Some things I have noticed is that the function value val can return some really small imaginary components that should be zero. If I force the value to be real via either Chop or Re, I keep getting an Indeterminate value.

How can I get around this?

MWE:

s[0 | 0.] = 0;
s[x_] = x Log2[x];
SetAttributes[s, Listable]

A = {{0.0632, 0, 0, 0}, {0, 0.3065, 0, 0}, {0, 0, 0.0632, 0}, {0, 0, 0, 0.3065}};
B = {{0.3065, 0, 0, 0}, {0, 0.0632, 0, 0}, {0, 0, 0.3065, 0}, {0, 0, 0, 0.0632}};
T = {{2 k Cos[\[Theta]]^2 + k + 1, 0, 0, 6 k Cos[\[Theta]]^2 - k - 1}, {0, 2 k Sin[\[Theta]]^2 + 3 - 3 k, 0, 0}, {0, 0, 2 k Sin[\[Theta]]^2 + 3 - 3 k, 0}, {6 k Cos[\[Theta]]^2 - k - 1, 0, 0, 2 k Cos[\[Theta]]^2 + k + 1}};

C1 = A.T.A;
C2 = A.T.B;
C3 = B.T.A;
C4 = B.T.B;
C5 = KroneckerProduct[{{1, 0}, {0, 0}}, C1] + KroneckerProduct[{{0, 1}, {0, 0}}, C2] + KroneckerProduct[{{0, 0}, {1, 0}}, C3] + KroneckerProduct[{{0, 0}, {0, 1}}, C4];

val = Total[s[Eigenvalues[C5]]] - Total[s[Eigenvalues[C1]]]  - Total[s[Eigenvalues[C4]]];


I then try to optimize val using:

NMinimize[{val, 0<=k<=1, 0<\[Theta]<=2\[Pi]}, {k, \[Theta]}]


which does not work. Even val /. {k->0.5, \[Theta]->0.1} fails. What is the workaround?

Problem. OP constructs a symmetric matrix C5 of size $$8 \times 8$$ that in theory has rank $$4$$. Therefore four eigenvalues are exactly zero in theory, but not in floating point arithmetic. Here is a simple example to illustrate:

SeedRandom[1];
M=RandomReal[{-1,1},{4,8}];
Eigenvalues[Transpose[M].M] (* 8x8 matrix that in theory has rank 4)
(*
{4.29967,
3.09331,
1.37365,
0.280941,
3.84658*10^-16,
-1.84698*10^-16,
-6.69403*10^-17,
-6.00174*10^-17} *)


Instead of four zero eigenvalues, we get four eigenvalues close to $$0$$ including slightly negative ones. Therefore, complex values are obtained when OP's function s[...] is applied. One approach is to use the fact that the nonzero eigenvalues of $$M^TM$$ are the eigenvalues of $$MM^T$$:

Eigenvalues[M.Transpose[M]]
(* {4.29967,
3.09331,
1.37365,
0.280941} *)


This yields directly the four eigenvalues of interest.

Solution. A modification of OP's code in this spirit is (please restart Mathematica before using this to avoid conflicts with previous definitions):

s[x_]:=Total[x*Log2[x]];

A=DiagonalMatrix[{0.0632,0.3065,0.0632,0.3065}];
B=DiagonalMatrix[{0.3065,0.0632,0.3065,0.0632}];
AB=MatrixPower[A.A+B.B,1/2];

f[k_?NumericQ,θ_?NumericQ]:=With[{T={
{2 k Cos[θ]^2+k+1,0,0,6 k Cos[θ]^2-k-1},
{0,2 k Sin[θ]^2+3-3 k,0,0},
{0,0,2 k Sin[θ]^2+3-3 k,0},
{6 k Cos[θ]^2-k-1,0,0,2 k Cos[θ]^2+k+1}}},
s[Eigenvalues[AB.T.AB]]-s[Eigenvalues[A.T.A]]-s[Eigenvalues[B.T.B]]
];


Using NMinimize yields a point where 6 k Cos[θ]^2-k-1 vanishes, that is, where $T$ is diagonal. I am guessing (without having analyzed in detail) that all points on that curve are minima, as indicated in this plot:

Show[{
Plot3D[f[k,θ],{k,0,1},{θ,0,Pi},AxesLabel->{"k","θ"}],
ParametricPlot3D[{1/(-1+6 Cos[θ]^2),θ,f[1/(-1+6 Cos[θ]^2),θ]},
{θ,0,Pi},PlotStyle->{Thick,Green}]
}]


The value along this curve is

With[{a=0.0632,b=0.3065},
8/Log[2]* ((a^2+b^2) Log[a^2+b^2]-a^2 Log[a^2]-b^2 Log[b^2])
]
(* 0.192641 *)


Explanation for Eigenvalues[AB.T.AB] ignoring details: OP's matrix C5 can be defined as follows:

X = ArrayFlatten[{{A,B}}];
C5 = Transpose[X].T.X;


Let S be a matrix square root of T, where S is also symmetric. Then equivalently C5 = Transpose[X].S.S.X. Denoting M = S.X this gives C5 = Transpose[M].M. By the example at the beginning of this answer, the nonzero eigenvalues of C5 are equal to the eigenvalues of M.Transpose[M] = S.X.Transpose[X].S = S.(A.A+B.B).S = S.AB.AB.S and they are equal (see this) to the eigenvalues of AB.S.S.AB = AB.T.AB.

• This gives a good indication of what is causing the issue. However, I'm not sure I follow why C5 in the OP is equivalent to AB.T.AB with AB defined in this answer.
– Sid
Dec 2, 2022 at 10:09
• They are not equal as matrices, but the nonzero eigenvalues are equal, see the edit. The numerical advantage is that AB.T.AB has size $4\times 4$, like in the example at the beginning. Dec 2, 2022 at 10:37
• I asked a follow-up question in mathematica.stackexchange.com/questions/276974/…, in case you are able to look at it?
– Sid
Dec 6, 2022 at 10:38
\$Version

(* "13.1.0 for Mac OS X x86 (64-bit) (June 16, 2022)" *)

Clear["Global*"]

s[0 | 0.] = 0;
s[x_] = x Log2[x];
SetAttributes[s, Listable]


Use exact values for A andB

A = {{0.0632, 0, 0, 0}, {0, 0.3065, 0, 0}, {0, 0, 0.0632, 0}, {0, 0, 0,
0.3065}} // Rationalize;
B = {{0.3065, 0, 0, 0}, {0, 0.0632, 0, 0}, {0, 0, 0.3065, 0}, {0, 0, 0,
0.0632}} // Rationalize;
T = {{2 k Cos[θ]^2 + k + 1, 0, 0, 6 k Cos[θ]^2 - k - 1}, {0,
2 k Sin[θ]^2 + 3 - 3 k, 0, 0}, {0, 0,
2 k Sin[θ]^2 + 3 - 3 k, 0}, {6 k Cos[θ]^2 - k - 1, 0, 0,
2 k Cos[θ]^2 + k + 1}};

C1 = A . T . A;
C2 = A . T . B;
C3 = B . T . A;
C4 = B . T . B;
C5 = KroneckerProduct[{{1, 0}, {0, 0}}, C1] +
KroneckerProduct[{{0, 1}, {0, 0}}, C2] +
KroneckerProduct[{{0, 0}, {1, 0}}, C3] +
KroneckerProduct[{{0, 0}, {0, 1}}, C4];

val = Total[s[Eigenvalues[C5]]] - Total[s[Eigenvalues[C1]]] -
Total[s[Eigenvalues[C4]]];

{min, arg} =
NMinimize[{val, 0 <= k <= 1, 0 < θ <= 2 π}, {k, θ}]

(* {0.192641, {k -> 0.218967, θ -> 3.41361}} *)
`