1

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

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Paul B. Slater
  • 2,335
  • 10
  • 16

2 Answers2

3

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}} *)
MikeY
  • 7,153
  • 18
  • 27
  • 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. – Paul B. Slater Feb 03 '20 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. – Paul B. Slater Feb 03 '20 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. – Roman Feb 03 '20 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. – Paul B. Slater Feb 03 '20 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. – Roman Feb 03 '20 at 22:11
  • Making use of NMaximize[{Log[obj], const}, {a[1], a[2], a[3], a[4], a[5], a[6], a[7]}, Method -> "DifferentialEvolution"], I obtained {-13.2005, {a[1] -> -0.298146, a[2] -> 0.51639, a[3] -> 0.298142, a[4] -> 0.247844, a[5] -> 0.42925, a[6] -> -0.400976, a[7] -> 0.694565}} which corresponds to obj == 1.84968*10^-6. – user64494 Feb 04 '20 at 15:20
  • The command of Maple DirectSearch:-GlobalOptima(log(a1^2a2^2a3^2a4^2a5^2a6^2a7^2), [9a1^2 + 9a2^2 + 9a3^2 <= 4, (2 + 3a3)(-4 + 6a3 + 9a4^2 + 9a5^2) + 18a1(a1 + 3a4a6 + 3a5a7) <= -18a2^2 + 9a2(-6a5a6 + 6a4a7) + 9(-2 + 3a3)(a6^2 + a7^2)],variables=[a1,a2,a3,a4,a5,a6,a7],maximize,evaluationlimit=20000,number=200); produces [0.1458275695631476*10^(-5), Vector[column](7, [-0.2225791355399396, 0.5673723704664783, 0.2701530424982285, -0.43802349890515774, 0.2691884463267418, -0.6310872094756825, -0.4756791517960255]), 703] . – user64494 Feb 04 '20 at 15:53
  • const /. {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} performs False. – user64494 Feb 04 '20 at 17:05
0

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.

Akku14
  • 17,287
  • 14
  • 32
  • 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. – Paul B. Slater Feb 03 '20 at 22:32
  • Incidentally, I'm letting my machine run in an effort to get an exact solution--but probably foolhardy! – Paul B. Slater Feb 03 '20 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. – Paul B. Slater Feb 04 '20 at 01: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) – Paul B. Slater Feb 04 '20 at 02:31