Is a seven-dimensional constrained maximization problem too demanding for meaningful analysis?

I am interested in finding the maximum of

a[1]^2 a[2]^2 a[3]^2 a[4]^2 a[5]^2 a[6]^2 a[7]^2


subject to the (entanglement) constraint

9 (a[1]^2 + a[2]^2 + a[3]^2) <= 4 && 9 (-2 a[2]^2 +
a[2] (-6 a[5] a[6] + 6 a[4] a[7]) + (-2 + 3 a[3]) (a[6]^2 +
a[7]^2)) >= (2 + 3 a[3]) (-4 + 6 a[3] + 9 a[4]^2 +
9 a[5]^2) + 18 a[1] (a[1] + 3 a[4] a[6] + 3 a[5] a[7])


I would suspect that this is too demanding a problem to solve symbolically. If so, can a numerical result of some degree of accuracy/precision be obtained for use in subsequent constrained integration problems.

Conceptually, at least, there is a companion 14-dimensional problem (cf. Maximize a six-dimensional function subject to joint positive-semidefiniteness constraints) to maximize

Abs[a[1] b[1]] + Abs[a[2] b[2]] + Abs[a[3] b[3]] + Abs[a[4] b[4]] +Abs[a[5] b[5]] + Abs[a[6] b[6]] + Abs[a[7] b[7]


subject to the intersection of the constraint above and a second version of it in which the $$a$$'a are replaced by $$b$$'s.

As a matter of some background, these problems are in pursuit of an attempt to implement an $$8 \times 8$$ Hadamard extension of an already pursued study based on a $$4 \times 4$$ Hadamard matrix. https://mathoverflow.net/questions/351790/are-n-times-n-special-orthogonal-matrices-all-the-entries-of-which-have-the

It seems to readily solve numerically. Does the answer make sense to you?

obj = a[1]^2 a[2]^2 a[3]^2 a[4]^2 a[5]^2 a[6]^2 a[7]^2

const = 9 (a[1]^2 + a[2]^2 + a[3]^2) <= 4 &&
9 (-2 a[2]^2 + a[2] (-6 a[5] a[6] + 6 a[4] a[7]) + (-2 + 3 a[3]) (a[6]^2 +a[7]^2))
>= (2 + 3 a[3]) (-4 + 6 a[3] + 9 a[4]^2 + 9 a[5]^2) + 18 a[1] (a[1] + 3 a[4] a[6] + 3 a[5] a[7])

res = NMaximize[{obj, const}, {a[1], a[2], a[3], a[4], a[5], a[6], a[7]}]


Instantaneous solve...

{1.15429*10^-10, {a[1] -> -0.100557, a[2] -> 0.367443,
a[3] -> -0.056478, a[4] -> 0.298236, a[5] -> 0.331095,
a[6] -> -0.179873, a[7] -> 0.289865}}


EDIT

Switching over to FindMaximum[ ] using the prior results as an initial condition...

init = Array[ao, 7];
var = Array[a, 7];

init = var /. (res // Last);

FindMaximum[{obj, const}, Transpose@{var, init}]

{7.58829*10^-11, {a[1] -> -0.0950384, a[2] -> 0.347098,
a[3] -> -0.0536012, a[4] -> 0.296238, a[5] -> 0.329195,
a[6] -> -0.176924, a[7] -> 0.28554}}


EDIT BASED ON UPDATED CONSTRAINTS

const2 = 9 (a[1]^2 + a[2]^2 + a[3]^2) <= 8 &&
9 (-4 a[1]^2 + 6 Sqrt[2] a[2] (-a[5] a[6] + a[4] a[7])
- 6 Sqrt[2] a[1] (a[4] a[6] + a[5] a[7]) + (-4 + 3 Sqrt[2] a[3])
(a[6]^2 + a[7]^2)) >= 36 a[3]^2 + 27 Sqrt[2] a[3] (a[4]^2 + a[5]^2)
+  4 (-8 + 9 a[2]^2 + 9 a[4]^2 + 9 a[5]^2)


res = NMaximize[{obj, const2}, res = NMaximize[{obj, const2}, var]]

(* {7.4268*10^-6, {a[1] -> 0.480501, a[2] -> 0.595736, a[3] -> 0.237296,
a[4] -> -0.268952, a[5] -> -0.562672, a[6] -> 0.528649,
a[7] -> 0.501494}} *)


Use this to initialize FindMaximum[ ]

FindMaximum[{obj, const2}, Transpose@{var, var /. res[[2]]}]

(* {0.000236499, {a[1] -> 0.730195, a[2] -> 0.421605, a[3] -> 0.421605,
a[4] -> -0.350482, a[5] -> -0.607013, a[6] -> 0.982126,
a[7] -> 0.567067}} *)

• Since it's rather unexplored territory, can't really vouch for its plausibility, but does seem off-hand not unreasonable. The square of the result ($\approx 1*10^-20$) would be used in subsequent constrained integrations. But is the result stable with the use of greater precision and increased MaxIterations,..--or should one even go in such a direction? As something of a pipe-dream, one might hope for sufficient precision to intuit an underlying exact value. Commented Feb 3, 2020 at 20:53
• Well, if I use WorkingPrecision->16 in the set-up of MikeY, I get 6.996813688626978*10^-7--not too encouraging, I guess. Commented Feb 3, 2020 at 21:22
• I find a much higher solution by using the numerically more stable NMaximize[{obj^(1/7), const}, {a[1], a[2], a[3], a[4], a[5], a[6], a[7]}], giving obj == 0.00015333. Playing with the Method may give even higher results. Commented Feb 3, 2020 at 21:51
• Very interesting, Roman!--and I don't think an implausible result. What led you to such a strategy, may I ask? I guess you're thinking "higher" increases plausibility. Commented Feb 3, 2020 at 22:07
• @PaulB.Slater what led me to this strategy is the thought that your objective function is numerically not well-behaved because even a small change in all parameters (say, dividing all parameters by two) gives a huge change in the objective function. First I tried using Log[obj] as an objective function, but that's not good because singular if any parameter is zero. Any monotonous transformation of obj will do, and obj^(1/7) or obj^(1/14) seem natural because they scale as the parameters. Commented Feb 3, 2020 at 22:11

I think you have to prevent to get into a lokal maximum. Since a[1] to a[3] is spherical, a[1]^2 a[2]^2 a[3]^2 has 8 equivalent maxima. Look at the two-dimensional analagon: Plot3D[a[1]^2*a[2]^2, {a[1], -1, 1}, {a[2], -1, 1}, RegionFunction -> Function[{x, y, z}, 9 (x^2 + y^2) <= 4]] .

Maxima and minima for the a[1],a[2] and a[3] show, they are all +- 2/3 and not restricted by the socond condition.

max = NMaximize[{a[1], const}, {a[1], a[2], a[3], a[4], a[5], a[6],
a[7]}]

(*   {0.666667, {a[1] -> 0.666667, a[2] -> 1.20991*10^-6,
a[3] -> -2.05127*10^-6, a[4] -> 0.0747194, a[5] -> -0.585257,
a[6] -> -0.074698, a[7] -> 0.58523}}    and so on.   *)


Find this octant where obj is maximal by trying all 8.

max = NMaximize[{obj,
const && a[1] > 0 && a[2] > 0 && a[3] < 0}, {a[1], a[2], a[3],
a[4], a[5], a[6], a[7]}]

(*   {1.17934*10^-6, {a[1] -> 0.438758, a[2] -> 0.298584,
a[3] -> -0.395921, a[4] -> -0.726167, a[5] -> 0.411118,
a[6] -> 0.189507, a[7] -> -0.370075}}   *)


and 1.17934*10^-6 as overall maximum.

• Wow!--lots of room for ingenuity/creativity here--but still (it would seem) quite a difference from the Roman result of 0.00015333 (being lesser--maybe not a good sign). Thanks everyone for all the effort--hope to a good purpose. Commented Feb 3, 2020 at 22:32
• Incidentally, I'm letting my machine run in an effort to get an exact solution--but probably foolhardy! Commented Feb 3, 2020 at 22:33
• Sorry, one and all--it seems that my constraint was flawed. I ran a FindInstance[const,Array[a,7],i] command in a loop, and for i = 2, I got {a[1] -> -(41/151), a[2] -> 13/83, a[3] -> Sqrt[489400354]/37599, a[4] -> -(39/10), a[5] -> 17/5, a[6] -> 598377/(10 (-25066 + Sqrt[489400354])), a[7] -> -(23487/(2 (-25066 + Sqrt[489400354])))}, which gave me a whopping 724.052. So, I went back and found that I had not done some scaling in setting up a $3 \times 3$ density matrix correctly. I'll give what I now believe to be the proper constraint in the next comment-but pl. accept my apologies. Commented Feb 4, 2020 at 1:22
• The proper constraint now appears to be const= 9 (a[1]^2 + a[2]^2 + a[3]^2) <= 8 && 9 (-4 a[1]^2 + 6 Sqrt[2] a[2] (-a[5] a[6] + a[4] a[7]) - 6 Sqrt[2] a[1] (a[4] a[6] + a[5] a[7]) + (-4 + 3 Sqrt[2] a[3]) (a[6]^2 + a[7]^2)) >= 36 a[3]^2 + 27 Sqrt[2] a[3] (a[4]^2 + a[5]^2) + 4 (-8 + 9 a[2]^2 + 9 a[4]^2 + 9 a[5]^2) Commented Feb 4, 2020 at 2:31