# Graph/Construct (John) ellipsoids circumscribing and inscribing a certain 3D convex set

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$$.

• Are you sure that this is a Mathematica question? I only see a math problem. Nov 13, 2020 at 15:19
• 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. Nov 13, 2020 at 18:56
• You can improve the RegionPlot3D by adding PlotPoints->100 ignoring for the moment imaginary results your obtain with the code you referenced. Nov 14, 2020 at 18:02
• A related question. Dec 24, 2020 at 13:23

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]}]


Note the documentation itself warns:

"FastEllipsoid" gives a bounding Ellipsoid, not necessarily with the minimal volume

• Interesting--certainly did not know of this command! Nov 16, 2020 at 15:02
• I assume there are 4 points touching the ellipsoid. Didn't see any option to find them however. Nov 16, 2020 at 15:26
• 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]]]}]  Nov 16, 2020 at 15:45

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:

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
]


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]

• 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? Nov 19, 2020 at 15:55
• 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. Nov 19, 2020 at 21:57
• 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. Nov 19, 2020 at 23:12

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}]


• 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). Nov 20, 2020 at 13:02
• 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. Nov 20, 2020 at 13:43
• 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.) Nov 21, 2020 at 15:12
• 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. Nov 21, 2020 at 16:57
• Good--glad you're enjoying this whole endeavor. Thanks again for your interest. Nov 21, 2020 at 23:44

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]


• 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.) Nov 14, 2020 at 19:00
• 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 Nov 14, 2020 at 19:59
• Thanks, Dominic! I will study this python code. Nov 14, 2020 at 23:10
• 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. Nov 15, 2020 at 10:33
• 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. :) Nov 15, 2020 at 11:03

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.

• 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. Nov 17, 2020 at 13:53
• Thanks, will check this out. Nov 17, 2020 at 18:01
• 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. Nov 17, 2020 at 18:45

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


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