2
$\begingroup$

A famous theorem JohnEllipsoids of Fritz John informs us that associated with a convex body are circumscribed and inscribed ellipsoids of minimal and maximal volumes.

Now, a body--argued to be convex in the answer of Nathaniel Johnston to SpectraConvexity--is the set of ordered spectra of absolutely separable two-qubit states. This set is defined by the constraint

1 > x && x >= y && y >= z && z >= 1 - x - y - z >= 0 && x - z < 2 Sqrt[y (1 - x - y - z)

Barring an explicit construction of the associated John ellipsoids--and possibly aiding in the search for them--graphical explorations (using the many tools of Mathematica, including RegionPlot3D, Ellipsoid and RegionMeasure) to find/approximate them would seem of interest.

There are two other sets of associated interest, also containing and contained within the convex body under examination. These are given by the constraints

 1 > x && x >= y && y >= z && z >= 1 - x - y - z >= 0 && x^2 + y^2 + (1 - x - y - z)^2 + z^2 < 3/8]

and

 1 > x && x >= y && y >= z && z >= 1 - x - y - z >= 0 && x^2 + y^2 + (1 - x - y - z)^2 + z^2 < 1/3]

Could these be the ellipsoids in question, and, if not, what geometric shapes might they be?

Here is a plot using RegionPlot3D of the three sets associated with the three constraints given above. I will try to incorporate the Ellipsoid command into these graphics also, as well as the use of RegionMeasure in seeking volumes.

ThreeSetPlot

Here is a very rough, preliminary effort to use Mathematica to explore the question posed. The plot shows an ellipsoid "close" to circumscribing the convex set of ordered spectra of absolutely separable two-qubit states. It seems very challenging, though, to construct that circumscribing ellipsoid of minimal volume--the existence of which is given by the Fritz John theorem. ("John ellipsoids are hard to compute" M-ellipsoids .) What is the objective function to minimize? Also, of course, there is the "dual" inscribed ellipsoid of maximal volume problem.

EllipsoidPlot

It's not clear to me if the Ellipsoid and RegionMeasure commands (among others) can be exploited in this quest.

The (Euclidean) volume of the convex set (of principal interest here) of ordered spectra of absolutely separable two-qubit states is $\approx 0.00227243$ (I should be able to obtain an exact value), while the volume of the ellipsoid shown in the last plot is $\frac{\pi }{150 \sqrt{15}} \approx 0.0054077$.

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ Are you sure that this is a Mathematica question? I only see a math problem. $\endgroup$ Commented Nov 13, 2020 at 15:19
  • 1
    $\begingroup$ Well, MarcoB, I was thinking in terms of the use of RegionPlot3D and the numerous other Mathematica graphics capabilities--an area in which I am not by any means expert. But I thought the question might intrigue Mathematica users with such expertise. I speculated that plots of fitting ellipsoids might be suggestive of further analyses. To make a long story short, I thought this might intrigue the many Mathematica/mathematics to explore the matter--using the abundant tools of Mathematica. Maybe I was excessive in doiing so. $\endgroup$ Commented Nov 13, 2020 at 18:56
  • 1
    $\begingroup$ You can improve the RegionPlot3D by adding PlotPoints->100 ignoring for the moment imaginary results your obtain with the code you referenced. $\endgroup$ Commented Nov 14, 2020 at 18:02
  • $\begingroup$ A related question. $\endgroup$ Commented Dec 24, 2020 at 13:23

6 Answers 6

Reset to default
3
$\begingroup$

More of an extended comment, but in case you weren't aware of the BoundingRegion functionality:

rm=RegionMember[ImplicitRegion[conditionABS,{x,y,z}]]; pts=RandomVariate[UniformDistribution[N[{{1/4,1/8 (2+Sqrt[6])},{1/24 (6-Sqrt[6]),1/8 (2+Sqrt[2])},{1/8 (2-Sqrt[2]),1/3}}]],10^5]; insidePts=Select[pts,rm]; fastEllipsoid=BoundingRegion[insidePts,"FastEllipsoid"] RegionMeasure[fastEllipsoid] Graphics3D[{{Opacity[0.5],fastEllipsoid},Point[insidePts]}]

enter image description here

Note the documentation itself warns:

"FastEllipsoid" gives a bounding Ellipsoid, not necessarily with the minimal volume
Improve this answer
$\endgroup$
3
  • $\begingroup$ Interesting--certainly did not know of this command! $\endgroup$ Commented Nov 16, 2020 at 15:02
  • $\begingroup$ I assume there are 4 points touching the ellipsoid. Didn't see any option to find them however. $\endgroup$ Commented Nov 16, 2020 at 15:26
  • 1
    $\begingroup$ Unfortunately, I think it might only guarantee to have 1 point touching.. In any case here's a (slow) way to find the closest points to the bounding ellipsoid: rd=RegionDistance[RegionBoundary[fastEllipsoid]]; dists=rd[insidePts]; Graphics3D[{{Opacity[0.5],fastEllipsoid},Point[insidePts],Red,Point[MinimalBy[Thread[{insidePts,dists}],Last,10][[All,1]]]}] $\endgroup$ Commented Nov 16, 2020 at 15:45
1
$\begingroup$

Here is a method to circumscribe an ellipsoid, given 4 points:

Form the conditions:

cond = 1 > x && x >= y && y >= z && z >= 1 - x - y - z >= 0 && x - z <= 2 Sqrt[y (1 - x - y - z)];

,where we changed < to <=, we first determine the four extremal points using Minimize and Maximize: E.g. Maximize[{y, cond}, {x, y, z}]. This gives the 4 points:

pts={{1/3, 1/3, 1/3}, {1/4, 1/4, 1/4}, {1/2, 1/6, 1/ 6}, {1/8 (2 + Sqrt[2]), 1/8 (2 + Sqrt[2]), 1/2 (1 + 1/4 (-2 - Sqrt[2]))}}//N;

Next we determine the two points that a farthest apart. In our case this are pts[[2]] and pts[[4]]. We choose the midpoint of the line pts[[2]] to pts[[4]] as the center of our ellipsoid: com (center of mass). And half of the distance will be the largest half axis: a3 of our ellipsoid:

com = (pts[[2]] + pts[[4]])/2 // N; a3 = Norm[pts[[2]] - pts[[4]]]/2 // N;

To make the following calculations easier, we translate the points so that com lies at the origin. And then, we rotate the coordinate system so that the half axis a3 points in the z direction:

pts1 = (# - com) & /@ pts // N; pts2 = (r2 = RotationMatrix[{pts1[[2]] - pts1[[4]], {0, 0, 1}}]).# & /@ pts1;

Now we determine which of the points 1 or 3 (in our case point 3) is further from the origin and rotate around the z-axis so that this points lays in the y-z plane:

pts3 = (r3 = RotationMatrix[ ArcTan[pts2[[3, 1]], pts2[[3, 2]]], {0, 0, 1}]).# & /@ pts2;

Next we determine the half axis along the y axis so that point 3 lays on the ellipse, the y-z plane cuts out of the ellipsoid:

a2 = Sqrt[pts3[[3, 2]]^2/(1 - (pts3[[3, 3]]/a3)^2)]

Now we determine the half axis a1 in direction of the x coordinates so that the last point 1 lays on the ellipsoid:

a1 = Sqrt[ pts3[[1, 1]]^2/(1 - (pts3[[1, 2]]/a2)^2 - (pts3[[1, 3]]/a3)^2)]

We now have all the data to plot the ellipsoid and the transformed points in the new coordinate system:

enter image description here

Finally, we need to transform the ellipsoid back to the original coordinates by writing the ellipsoid formula in the old coordinates:

fun[{x_, y_, z_}] = Total[((r3.r2.({x, y, z} - com))/{a1, a2, a3})^2];

With this we can now plot the ellipsoid in the original coordinates:

Show[ ContourPlot3D[ fun[{x, y, z}] == 1, {x, .1, .6}, {y, .1, .55}, {z, -.1, .4}, AxesLabel -> {"x", "y", "z"}, ContourStyle -> Opacity[0.5], Mesh -> None] , Graphics3D[{PointSize[0.03], Point[pts]}, Axes -> True] , reg ]

enter image description here

And for convenience, all the code in one piece:

cond = 1 > x && x >= y && y >= z && z >= 1 - x - y - z >= 0 && x - z <= 2 Sqrt[y (1 - x - y - z)]; pts = {{1/3, 1/3, 1/3}, {1/4, 1/4, 1/4}, {1/2, 1/6, 1/6}, {1/8 (2 + Sqrt[2]), 1/8 (2 + Sqrt[2]), 1/2 (1 + 1/4 (-2 - Sqrt[2]))}} // N; com = (pts[[2]] + pts[[4]])/2 // N; a3 = Norm[pts[[2]] - pts[[4]]]/2 // N; pts1 = (# - com) & /@ pts // N; pts2 = (r2 = RotationMatrix[{pts1[[2]] - pts1[[4]], {0, 0, 1}}]).# & /@ pts1; pts3 = (r3 = RotationMatrix[ ArcTan[pts2[[3, 1]], pts2[[3, 2]]], {0, 0, 1}]).# & /@ pts2; a2 = Sqrt[pts3[[3, 2]]^2/(1 - (pts3[[3, 3]]/a3)^2)] a1 = Sqrt[ pts3[[1, 1]]^2/(1 - (pts3[[1, 2]]/a2)^2 - (pts3[[1, 3]]/a3)^2)] Show[ContourPlot3D[ Total[({x, y, z}/{a1, a2, a3})^2] == 1, {x, -.2, .2}, {y, -.25, .25}, {z, -.2, .21}, AxesLabel -> {"x", "y", "z"}], Graphics3D[{PointSize[0.03], Point[pts3], Line[{pts3[[2]], pts3[[4]]}], Line[{{0, 0, 0}, pts3[[3]]}]}, Axes -> True]] fun[{x_, y_, z_}] = Total[((r3.r2.({x, y, z} - com))/{a1, a2, a3})^2]; reg = RegionPlot3D[ cond, {x, 1/4, 1/8 (2 + Sqrt[6])}, {y, 1/24 (6 - Sqrt[6]), 1/8 (2 + Sqrt[2])}, {z, 1/3, 1/8 (2 - Sqrt[2])}, PlotPoints -> 100]; Show[ContourPlot3D[ fun[{x, y, z}] == 1, {x, .1, .6}, {y, .1, .55}, {z, -.1, .4}, AxesLabel -> {"x", "y", "z"}, ContourStyle -> Opacity[0.5], Mesh -> None], Graphics3D[{PointSize[0.03], Point[pts]}, Axes -> True], reg]
Improve this answer
$\endgroup$
3
  • $\begingroup$ This is all great--and very impressive. Thanks! I tried calculating the volume of the ellipsoid, and got 0.0283059176...Right? It's known that the volume of the convex set is 1/576 (8 - 6 Sqrt[2] + 15 Sqrt[2] [Pi] - 48 Sqrt[2] ArcTan[Sqrt[2]])=0.00227243 mathematica.stackexchange.com/questions/234792/…, with 0.0283059 being consirably larger. Can your approach be extended to get the minimal volume (John) ellipsoid here? Would it have to go through these particular four points? How about the inscribed case? $\endgroup$ Commented Nov 19, 2020 at 15:55
  • $\begingroup$ You are right, the volume of the ellipsoid is 4/3 Pi a1 a2 a3. I also completed the code at the end, so that it is self contained. Unfortunately I do not think that my method can be used for the inscribed ellipsoid. The reason is, that I started with the four extreme points that are easily determined. For the inscribed problem, these points would correspond to the points where your body touches the ellipse. The body and the ellipse share there a tangent plane. But how to get these points is another problem. $\endgroup$ Commented Nov 19, 2020 at 21:57
  • $\begingroup$ So--using the ellipsoid volume formula,and your a1 a2 a3, the volume computes to 1/32 Sqrt[1/553 (29 + 12 Sqrt[2])] [Pi] \approx 0.0283059. (Amazingly, it seems also 1/32 Sqrt[1/ (29 - 12 Sqrt[2])] [Pi]. What happened to the 553 ?) So what are the 0.236089 and 0.186964 appearing in your new output? So, is this the requested (unique) circumscribing John ellipsoid? If so, further congratulations. $\endgroup$ Commented Nov 19, 2020 at 23:12
1
$\begingroup$

Here's my solution to finding inscribed ellipsoid:

  1. First generate one million points in a bounding region and select points inside conditionABS.
  2. Generate ConvexHull of points,
  3. Use LinearOptimization to construct polytopes,
  4. Run ConicOptimization to find ellipsoid.

However, had to change sign of the translation vector,d. Volume of ellipsoid (machine precision) is 0.001442. See Link to PF for largest ellipse in polygon for more information.

conditionABS = 1 > x && x >= y && y >= z && z >= 1 - x - y - z >= 0 && x - z < 2 Sqrt[y (1 - x - y - z)]; (* generate one million points in bounding region and select points \ inside conditionABS *) rm = RegionMember[ImplicitRegion[conditionABS, {x, y, z}]]; pts = RandomVariate[ UniformDistribution[ N[{{1/4, 1/8 (2 + Sqrt[6])}, {1/24 (6 - Sqrt[6]), 1/8 (2 + Sqrt[2])}, {1/8 (2 - Sqrt[2]), 1/3}}]], 10^6]; insidePts = Select[pts, rm]; (* generate a convex hull for the points *) mesh = ConvexHullMesh[insidePts]; meshP = Show[Graphics3D@{Opacity[0.02, Blue], mesh}, Axes -> True] (* Obtain polytope inequalities to represent the region *) {A, b} = LinearOptimization[0, {}, x \[Element] mesh, "LinearInequalityConstraints"]; Length@A (* use ConicOptimization to find max ellipsoid *) polyA = A; polyB = b; constraints = Table[Norm[polyA[[i]].c] + polyA[[i]].d <= polyB[[i]], {i, Length[polyA]}]; {cEllipse, dEllipse} = {c, d} /. ConicOptimization[-Tr[c], constraints, {c \[Element] Matrices[{3, 3}], d}] (* compute volume *) eVolume = 4 Pi/3 (Norm[cEllipse[[All, 1]]] Norm[cEllipse[[All, 2]]] Norm[cEllipse[[All, 3]]]) (* construct affine paramaterization for ellipsoid *) aFine[d_, m_, \[Theta]_, \[Phi]_] := d + m[[All, 1]] Cos[\[Theta]] Cos[\[Phi]] + m[[All, 2]] Cos[\[Theta]] Sin[\[Phi]] + m[[All, 3]] Sin[\[Theta]]; (* generate plots *) pp1 = ParametricPlot3D[ aFine[-dEllipse, cEllipse, t, p], {t, -Pi/2, Pi/2}, {p, 0, 2 Pi}] Show[{meshP, pp1}, Axes -> True, BoxRatios -> {1, 1, 1}]

enter image description here

Improve this answer
$\endgroup$
5
  • $\begingroup$ Quite skillful--thanks! More a fitting--rather than an exact analysis as that of Daniel Huber in his answer to the circumscribed case. Puzzled by the 0.11037 volume reported. The convex set has (Euclidean/flat) volume (as indicated at end of question) 1/576 (8 - 6 Sqrt[2] - 9 Sqrt[2] [Pi] + 24 Sqrt[2] ArcCos[1/3]) = 0.00227243 and the circumscribed ellipsoid of DH has volume 1/32 Sqrt[1/ (29 - 12 Sqrt[2])] [Pi] = 0.0283059 (reported in my comment to answer of DH). $\endgroup$ Commented Nov 20, 2020 at 13:02
  • $\begingroup$ Thanks for that. I see I made a mistake with the volume calculation. Volume of ellipsoid is pi/6(abc) of the semi-axes and not the sum. I updated my code above to reflect this and obtain as the volume, 0.000179. Would like to know if anyone gets a different value. The calculations (and the coding) are very new to me. $\endgroup$ Commented Nov 20, 2020 at 13:43
  • $\begingroup$ In a comment to his answer, Daniel Huber gives 4 Pi a b c/3 as the volume formula. (I thought I had posted a comment to this effect yesterday, but apparently not.) $\endgroup$ Commented Nov 21, 2020 at 15:12
  • $\begingroup$ Thanks. Apparently I don't know the difference between principal and semi-major axes. Volume is 4/3 Pi (abc) where a,b,c are semi-major axes which are specified nicely as the norms of column vectors of the matrix C above. Will update above. This gives volume ratio ellipsoid/convex set of 0.00144/0.00227. Would be nice though to have another way to compute it to compare the answer. Haven't worked with convex optimization before. $\endgroup$ Commented Nov 21, 2020 at 16:57
  • $\begingroup$ Good--glad you're enjoying this whole endeavor. Thanks again for your interest. $\endgroup$ Commented Nov 21, 2020 at 23:44
0
$\begingroup$

So for example, you're trying to find the circumscribed and inscribed elipsoids of minimal and maximal volumes for the region below (edited your code a bit)?

 conditionABS = 1 > x && x >= y && y >= z && z >= 1 - x - y - z >= 0 && x - z < 2 Sqrt[y (1 - x - y - z)]; RegionPlot3D[conditionABS, {x, 1/4, 1/8 (2 + Sqrt[6])}, {y, 1/24 (6 - Sqrt[6]), 1/8 (2 + Sqrt[2])}, {z, 1/3, 1/8 (2 - Sqrt[2])}, AxesLabel -> {Style["x", 16, Bold, Black], Style[ "y", 16, Bold, Black], Style[ "z", 16, Bold, Black]}, PlotPoints -> 100]

enter image description here

Improve this answer
$\endgroup$
6
  • $\begingroup$ OK! The term "conditionABS" fully defines the convex set of interest--for which we seek the circumscribed and inscribed (John) ellipsoids of minimal and maximal volumes, respectively. The ranges of the x, y, z coordinates employed in the RegionPlot3D command are only judgmental and certainly not fixed--and not key, in any case, to the central computational issue at hand.. (There is a mix of > and >= signs in the constraint, of perhaps no critical importance. But maybe it should just be simply of one or another type.) $\endgroup$ Commented Nov 14, 2020 at 19:00
  • 1
    $\begingroup$ Don't find any reference using Mathematica to compute John Ellipsoids. Here's a python example that may be helpful: docs.mosek.com/9.2/pythonfusion/case-studies-ellipsoids.html $\endgroup$ Commented Nov 14, 2020 at 19:59
  • $\begingroup$ Thanks, Dominic! I will study this python code. $\endgroup$ Commented Nov 14, 2020 at 23:10
  • $\begingroup$ Maybe post the problem in the math section of Stackexchange but suggest doing so from the perspective of the 3D mathematics entirely and not the code. $\endgroup$ Commented Nov 15, 2020 at 10:33
  • $\begingroup$ Also one more comment: They do a lot of geometry, including 3D, in (high school Math) AoPS. If it were my problem, I would take the time to write up a nicely worded (appropriate for the forum) question about John ellipsoids there. :) $\endgroup$ Commented Nov 15, 2020 at 11:03
0
$\begingroup$

Clearly not a full answer per se, but simply an amplification of two matters.

Firstly, it appears AreaVolumeRatio that the area/volume ratio of the convex set in question is 6. If so, this might assist in the identification of the nature of the set, if among known families of convex sets.

Secondly, the central ("troublesome") inequality constraint

 x - z < 2 Sqrt[y (1 - x - y - z)

is equivalent to the positive semidefiniteness of the $2 \times 2$ matrix,

P= {{2 (1 - x - y - z), -x + z}, {-x + z, 2 y}},

as pointed out by Nathaniel Johnston (citing work of R. Hildebrand) at the end of his answer to PositiveSemidefiniteness .

Now, might this matrix P be the one required (also denoted by P) in the python code "Inner and outer Löwner-John Ellipsoids" PythonCode noted by user Dominic is one of his comments to this question?

If so (I am somewhat skeptical at this point that P in some way represents a polytope, as seems to be required), then attempted implementation of the python code (utilizing ConfigurePythonForExternalEvaluate, it would seem), would clearly be in order.

Perhaps one can construct a matrix for which the positive-semidefiniteness condition is the defining constraint

1 > x && x >= y && y >= z && z >= 1 - x - y - z >= 0 && x - z < 2 Sqrt[y (1 - x - y - z)

given at the outset of the question. Possibly such a matrix would be the appropriate one to input to the python code.

A rather trivial way to obtain such a matrix would be to embed the indicated $2 \times 2$ matrix (yielding the main inequality constraint)

 P= {{2 (1 - x - y - z), -x + z}, {-x + z, 2 y}}

in the upper corner of an originally null $6 \times 6$ matrix and insert entries of 1-x, x-y, y-z and z-(1-x-y-z) into the remaining four diagonal positions.

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ Looked into solving max ellipse in polygon and found the function ConicOptimization: reference.wolfram.com/language/ref/ConicOptimization.html. Perhaps you can adapt your problem (linearize it maybe) using this function. $\endgroup$ Commented Nov 17, 2020 at 13:53
  • $\begingroup$ Thanks, will check this out. $\endgroup$ Commented Nov 17, 2020 at 18:01
  • $\begingroup$ There seem to be several Optimization commands. Not sure which one--if any--is the right one to use for the two (cricumscribing and inscribing) problems at hand. $\endgroup$ Commented Nov 17, 2020 at 18:45
0
$\begingroup$

Again, not by any means, the requested construction of the two ellipsoids in question, but rather an effort to bring two developments of some interest in this matter to attention.

Firstly, user Dominic in a comment here has noted a sophisticated (Mosek--a software package) python code entitled "Inner and outer Löwner-John Ellipsoids" Mosekpythoncode. Neither being a python user--nor much of an expert in the underlying optimization procedures--I posted this question pythonQuestion .

A user Michal Adamaszek--apparently affiliated with Mosek--commented:

"The Mosek code is intended for the ellipsoid inscribed in a polytope P. If P is convex but not a polytope than it may or may not be possible, depending on if you can rewrite the "for all u" part into something more manageable. It seems that your set has an SDP representation, so at the very least you can get an approximation by sampling sufficiently many u and constraining the corresponding Cu+d to lie in P." (By the "SDP representation", I believe is meant the $6 \times 6$ matrix

{{2 (1 - x - y - z), -x + z, 0, 0, 0, 0}, {-x + z, 2 y, 0, 0, 0, 0}, {0, 0, 1 - x, 0, 0, 0}, {0, 0, 0, x - y, 0, 0}, {0, 0, 0, 0, y - z, 0}, {0, 0, 0, 0, 0, -1 + x + y + 2 z}}

constructed in my previous "answer".

I replied:

"Thanks very, very much Michal Adamaszek--exactly the expertise I was hoping to obtain through posing the question. Since I'm not a python user, I may have to struggle some more with implementing your suggested approach. At this point, Ihave no firm knowledge as to whether or not P is a polytope--I would suspect that was "too good to be true". Again, within my limited understanding, whether or nor P is a polytope is itself a challenging question."

As to the second development I wanted to highlight here, it is now known that the area/volume ratio of the convex ("ordered spectra") set is 6. This is an immediate consequence of the obtaining (by user JimB) AreaVolumeRatio of the volume of the set as

1/576 (8 - 6 Sqrt[2] - 9 Sqrt[2] π + 24 Sqrt[2] ArcCos[1/3]) ,

coupled with my previous finding of the area as six times this expression.

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.