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

I am interested in finding the maximum of

a^2 a^2 a^2 a^2 a^2 a^2 a^2


subject to the (entanglement) constraint

9 (a^2 + a^2 + a^2) <= 4 && 9 (-2 a^2 +
a (-6 a a + 6 a a) + (-2 + 3 a) (a^2 +
a^2)) >= (2 + 3 a) (-4 + 6 a + 9 a^2 +
9 a^2) + 18 a (a + 3 a a + 3 a a)


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 b] + Abs[a b] + Abs[a b] + Abs[a b] +Abs[a b] + Abs[a b] + Abs[a b


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^2 a^2 a^2 a^2 a^2 a^2 a^2

const = 9 (a^2 + a^2 + a^2) <= 4 &&
9 (-2 a^2 + a (-6 a a + 6 a a) + (-2 + 3 a) (a^2 +a^2))
>= (2 + 3 a) (-4 + 6 a + 9 a^2 + 9 a^2) + 18 a (a + 3 a a + 3 a a)

res = NMaximize[{obj, const}, {a, a, a, a, a, a, a}]


Instantaneous solve...

{1.15429*10^-10, {a -> -0.100557, a -> 0.367443,
a -> -0.056478, a -> 0.298236, a -> 0.331095,
a -> -0.179873, a -> 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 -> -0.0950384, a -> 0.347098,
a -> -0.0536012, a -> 0.296238, a -> 0.329195,
a -> -0.176924, a -> 0.28554}}


EDIT BASED ON UPDATED CONSTRAINTS

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


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

(* {7.4268*10^-6, {a -> 0.480501, a -> 0.595736, a -> 0.237296,
a -> -0.268952, a -> -0.562672, a -> 0.528649,
a -> 0.501494}} *)


Use this to initialize FindMaximum[ ]

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

(* {0.000236499, {a -> 0.730195, a -> 0.421605, a -> 0.421605,
a -> -0.350482, a -> -0.607013, a -> 0.982126,
a -> 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. 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. Feb 3, 2020 at 21:22
• I find a much higher solution by using the numerically more stable NMaximize[{obj^(1/7), const}, {a, a, a, a, a, a, a}], giving obj == 0.00015333. Playing with the Method may give even higher results. 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. 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. Feb 3, 2020 at 22:11

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

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

max = NMaximize[{a, const}, {a, a, a, a, a, a,
a}]

(*   {0.666667, {a -> 0.666667, a -> 1.20991*10^-6,
a -> -2.05127*10^-6, a -> 0.0747194, a -> -0.585257,
a -> -0.074698, a -> 0.58523}}    and so on.   *)


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

max = NMaximize[{obj,
const && a > 0 && a > 0 && a < 0}, {a, a, a,
a, a, a, a}]

(*   {1.17934*10^-6, {a -> 0.438758, a -> 0.298584,
a -> -0.395921, a -> -0.726167, a -> 0.411118,
a -> 0.189507, a -> -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. Feb 3, 2020 at 22:32
• Incidentally, I'm letting my machine run in an effort to get an exact solution--but probably foolhardy! 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 -> -(41/151), a -> 13/83, a -> Sqrt/37599, a -> -(39/10), a -> 17/5, a -> 598377/(10 (-25066 + Sqrt)), a -> -(23487/(2 (-25066 + Sqrt)))}, 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. Feb 4, 2020 at 1:22
• The proper constraint now appears to be const= 9 (a^2 + a^2 + a^2) <= 8 && 9 (-4 a^2 + 6 Sqrt a (-a a + a a) - 6 Sqrt a (a a + a a) + (-4 + 3 Sqrt a) (a^2 + a^2)) >= 36 a^2 + 27 Sqrt a (a^2 + a^2) + 4 (-8 + 9 a^2 + 9 a^2 + 9 a^2) Feb 4, 2020 at 2:31