1
$\begingroup$

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?

$\endgroup$

2 Answers 2

2
$\begingroup$

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}]
}]

enter image description here

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.

$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – Sid
    Dec 2, 2022 at 10:09
  • $\begingroup$ 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. $\endgroup$
    – user293787
    Dec 2, 2022 at 10:37
  • $\begingroup$ I asked a follow-up question in mathematica.stackexchange.com/questions/276974/…, in case you are able to look at it? $\endgroup$
    – Sid
    Dec 6, 2022 at 10:38
2
$\begingroup$
$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}} *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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