Solution of 20-variable Constrained NMaximize with Method->"DifferentialEvolution" does not fully satisfy the constraints

Question

I want to find the maximum of

d2=    (Sqrt[Abs[a[1]^2 b[1]^2]] + Sqrt[Abs[a[2]^2 b[2]^2]] + Sqrt[Abs[a[4]^2 b[4]^2]] + Sqrt[Abs[a[5]^2 b[5]^2]] + Sqrt[Abs[a[6]^2 b[6]^2]] + Sqrt[Abs[a[7]^2 b[7]^2]] + (1/Sqrt[2])(\[Sqrt]Abs[a[3]^2 b[3]^2 + a[8]^2 b[8]^2 + a[9]^2 b[9]^2 + a[10]^2 b[10]^2 - \[Sqrt](((a[3] b[3] + a[8] b[8])^2 + (a[9] b[9] - a[10] b[10])^2) ((a[3] b[3] - a[8] b[8])^2 + (a[9] b[9] +a[10] b[10])^2))]) + (1/Sqrt[2])(\[Sqrt]Abs[a[3]^2 b[3]^2 + a[8]^2 b[8]^2 + a[9]^2 b[9]^2 + a[10]^2 b[10]^2 + \[Sqrt](((a[3] b[3] + a[8] b[8])^2 + (a[9] b[9] - a[10] b[10])^2) ((a[3] b[3] - a[8] b[8])^2 + (a[9] b[9] +a[10] b[10])^2))]))^2

subject to a very large set (LeafCount of 27947) of constraints designated "vcf" in the notebook

https://www.wolframcloud.com/obj/slater/Published/MaximizeNorms

The command

G1 = NMaximize[{d2, vcf}, Join[Array[a, 10], Array[b, 10]]]

yields the result

{0.1565, {a[1] -> -0.217171, a[2] -> 0.228215, a[3] -> 0.242746, a[4] -> 0.185082, a[5] -> -0.194495, a[6] -> -0.239223, a[7] -> -0.25139, a[8] -> -0.236464, a[9] -> 0.157249,  a[10] -> 0.152502, b[1] -> -0.187483, b[2] -> 0.197017, b[3] -> 0.217313, b[4] -> 0.201103, b[5] -> 0.19137,b[6] -> 0.269263, b[7] -> -0.256231, b[8] -> 0.197663, b[9] -> 0.0729522, b[10] -> 0.17673}}

which does in fact satisfy the constraint vcf.

However, the command

G2 = NMaximize[{d2, vcf}, Join[Array[a, 10], Array[b, 10]], Method -> "DifferentialEvolution"]

yields the larger (and more plausible/attractive) result

{0.186769, {a[1] -> -0.268552, a[2] -> -0.255555, a[3] -> 0.227544,a[4] -> -0.18653, a[5] -> 0.196016, a[6] -> -0.203221, a[7] -> 0.193385, a[8] -> -0.316546, a[9] -> -0.162042, a[10] -> -0.153076, b[1] -> -0.237133, b[2] -> -0.249411,  b[3] -> 0.215763, b[4] -> 0.20219, b[5] -> -0.192236, b[6] -> -0.193805, b[7] -> 0.203839, b[8] -> -0.31299, b[9] -> -0.026066, b[10] -> -0.201624}}

but does not fully satisfy the constraint vcf.

(I believe that there is an exact value--possibly rational in nature [not too much larger than the 0.186769]--for this [quantum-information-theoretic] problem, but it certainly seems to be too demanding a task to obtain it. So, I'd like to "believe" the 0.186769 result, but the evidence for it is clearly less than 100% convincing.)

There is a second expression

d1=a[1]^2 a[2]^2 a[4]^2 a[5]^2 a[6]^2 a[7]^2 b[1]^2 b[2]^2 b[4]^2 b[5]^2 b[6]^2 b[7]^2 (a[3] a[8] b[3] b[8] - a[9] a[10] b[9] b[10])^2

which it is also of interest to maximize with respect to the same constraints.

Answer 1

These constraints are so complicated that it is very easy to get caught in local maxima. Maybe RandomSearch can find a way out.

I got a much higher maximum than 0.186, namely 3.7354. (First I did some simplifications). Verification yields "True".

d2 = d2 // 
  Simplify[#, Thread[Join[Array[a, 10], Array[b, 10]] \[Element] Reals]] &

vcfl = vcf[[1]] // Rationalize // 
    Simplify[#, Thread[Join[Array[a, 10], Array[b, 10]] \[Element] Reals]] &;

(H1 = NMaximize[{d2, vcfl}, Join[Array[a, 10], Array[b, 10]], 
 Method -> "RandomSearch"]) // Timing

(*   {209.609, {3.7354, {a[1.] -> -0.240179, a[2.] -> -0.228551, 
a[3.] -> 0.196117, a[4.] -> 0.222086, a[5.] -> 0.211336, 
a[6.] -> 0.211402, a[7.] -> 0.222153, a[8.] -> 0.145827, 
a[9.] -> -0.195623, a[10.] -> -0.281385, b[1.] -> -0.836558, 
b[2.] -> -1.1664, b[3.] -> 0.971168, b[4.] -> -1.37882, 
b[5.] -> 0.96628, b[6.] -> -1.12641, b[7.] -> -1.2148, 
b[8.] -> 0.811886, b[9.] -> -0.481526, b[10.] -> 0.812585}}}   *)

vcfl /. H1[[2]]

(*   True   *)

May be you can play with "MaxIterations" to get an even higher maximum.

Answer 2

Well, here's something of a (better/greater) "answer" that I arrived at a number of hours before the quite surprising one of Akku14. If his/her answer does stand up, it makes this one look rather "irrelevant". (But a series of comments of mine to his/her answer now show that by employing vcf[[1]], and not vcf, the ten b[i] variables wound up being unconstrained.)

I tried reducing the problem from a 20-variable one to a 10-variable one, by the command Do[b[i] = a[i], {i, 1, 10}]. (It also works "symmetrically" if one uses Do[a[i] = b[i], {i, 1, 10}].) Then, the use of "DifferentialEvolution" gave me a solution (NOW satisfying the constraints) of

{0.184342, {a[1] -> -0.251841, a[2] -> 0.239653, a[3] -> -0.241901, a[4] -> -0.22769, a[5] -> 0.21667, a[6] -> -0.185068, a[7] -> 0.194481, a[8] -> -0.281275, a[9] -> -0.0274505, a[10] -> -0.0274505}},

but still inferior to 0.186769.

However, the command NMaximize[{d2, vcf}, Array[a, 10]] gave

{0.194815, {a[1] -> -0.255072, a[2] -> -0.268045, a[3] -> 0.228979, a[4] -> 0.186145, a[5] -> 0.195611, a[6] -> -0.194521, a[7] -> 0.204414, a[8] -> -0.315437, a[9] -> -0.149244, a[10] -> -0.149244}}

with 0.194815 > 0.186769, also satisfying the constraints. (If I explicitly use Method->"NelderMead", I get the same result, so it seems that this was the default Method employed.)

This line of reasoning does not explain why the original use of "DifferentialEvolution" appeared flawed (in not fully satisfying the constraints), but in giving a greater-valued answer, seemed to make it a less pressing question.

Also, note that in moving from the 20-variable problem to these 10-variable problems, the [now somewhat notorious] set ("vcf") of constraints can be preliminary simplified.

So, can a maximum greater than 0.194815 be found?

Also, the maximization of the indicated

d1=a[1]^2 a[2]^2 a[4]^2 a[5]^2 a[6]^2 a[7]^2 b[1]^2 b[2]^2 b[4]^2 b[5]^2 b[6]^2 b[7]^2 (a[3] a[8] b[3] b[8] - a[9] a[10] b[9] b[10])^2

is of interest. The maximum appears (from other analyses) to be greater than $2 \cdot 10^{-21}$.