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

Question

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

Answer 1

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}} *)

Answer 2

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.