10
$\begingroup$

I recently needed to deal with a large data set of 600,000 points in three dimensions. My task was to find a convex hull for this data. I have Mathematica 10, so I could use the function ConvexHullMesh; I obtained this:

screenshot convex hull

I was wondering if there is some way to find a smooth convex hull (maybe an ellipsoid) for my data.

$\endgroup$
  • 4
    $\begingroup$ What exactly do you mean by 'smooth convex hull', or, how is an ellipsoid a convex hull for your data? Could you explain that a bit. $\endgroup$ – user21 Aug 21 '14 at 7:39
  • 4
    $\begingroup$ I second user21's question; there is no standard notion of a "smooth convex hull". If you really really want an ellipsoid, you could try finding the minimum-volume ellipsoid enclosing your set of points. $\endgroup$ – Rahul Aug 21 '14 at 9:11
  • 1
    $\begingroup$ Yes, I mean, smooth convex hull as differential function. So, the other comment is probably that I need. I'll read the paper, and If I success, I write hear my trie. $\endgroup$ – jonaprieto Aug 21 '14 at 15:49
  • $\begingroup$ There are routines to calculate MVE ellipsoids in MATLAB. You could just translate one to Mathematica. $\endgroup$ – user24615 Jan 9 '15 at 18:17
10
$\begingroup$

Minimum Volume Ellipsoid

Translated from here, this uses the Khachiyan algorithm, and should work for any dimension.

    MinVolEllipse[P_, tolerance_] := 
     Module[{d, n, Q, count, err, u, X, M, maximum, j, stepSize, newu, U, A, c},

      {n, d} = Dimensions[P];
      Q = Append[1] /@ P;

      count = 1;
      err = 1;
      u = ConstantArray[1./n, n];

      While[err > tolerance,
       X = Q\[Transpose].DiagonalMatrix[u].Q;
       M = Diagonal[Q.Inverse[X].Q\[Transpose]];
       maximum = Max[M];
       j = Position[M, maximum][[1, 1]];
       stepSize = (maximum - d - 1)/((d + 1) (maximum - 1));
       newu = (1 - stepSize) u;
       newu[[j]] += stepSize;
       count += 1;
       err = Norm[newu - u];
       u = newu;
       ];

      U = DiagonalMatrix[u];

      A = (1/d) Inverse[P\[Transpose].U.P - Outer[Times, u.P, u.P]];
      c = u.P;

      {c, A}
      ]

Usage:

    pts = RandomVariate[
       MultinormalDistribution[RandomReal[{-1, 1}, {2}], 
        With[{m = RandomReal[{0, 1}, {2, 2}]}, m.m\[Transpose]]], 500];

    P = MeshCoordinates[ConvexHullMesh[pts]];
    tolerance = 0.0001;

    {c, A} = MinVolEllipse[P, tolerance];

    X = {x, y};
    Show[
     ConvexHullMesh[pts],
     Graphics[{
       Point[pts],
       {Red, Point[P]}
       }],
     ContourPlot[(X - c).A.(X - c) == 1, {x, -4, 4}, {y, -4, 4}]
     ]

enter image description here

In 3D:

    pts = RandomVariate[
       MultinormalDistribution[RandomReal[{-1, 1}, {3}], 
        With[{m = RandomReal[{0, 1}, {3, 3}]}, m.m\[Transpose]]], 100];

    P = MeshCoordinates[ConvexHullMesh[pts]];
    tolerance = 0.0001;

    {c, A} = MinVolEllipse[P, tolerance];

    X = {x, y, z};
    Show[
     ConvexHullMesh[pts, MeshCellStyle -> {{2, All} -> Opacity[0.5, Green]}],
     Graphics3D[{
       Point[pts],
       {Red, Point[P]}
       }],
     ContourPlot3D[(X - c).A.(X - c) == 1, {x, -3, 3}, {y, -3, 3}, {z, -3, 3},
      ContourStyle -> Directive[Red, Opacity[0.2]]
      ]
     ]

enter image description here

$\endgroup$
8
$\begingroup$

Here is a faster implementation of the Khachiyan minimum-volume ellipsoid algorithm, using updating formulae presented in this paper:

mve[pts_?MatrixQ, tol_: 1.*^-8] := 
    Module[{prec = Precision[pts], bet, bm, c, d, del, dp, h, j, kap, n,
            qj, qm, sc, sig, u, zero},
           zero = SetPrecision[0, prec]; {n, d} = Dimensions[pts]; dp = d + 1;
           qm = PadRight[pts, {n, dp}, N[1, prec]]; 
           u = N[ConstantArray[1/n, n], prec];
           bm = Transpose[qm].DiagonalMatrix[u].qm; 
           kap = Diagonal[qm.LinearSolve[bm, Transpose[qm]]];
           c = FixedPoint[(j = Ordering[kap, -1]; qj = Extract[qm, j];
                           sc = 1/Extract[kap, j]; sig = (1 - sc dp)/(1 - sc);
                           bet = sig/dp; h = 1 - bet; del = dp (1 - sc)/d;
                           kap = del (kap - sc sig (qm.LinearSolve[bm, qj])^2);
                           bm = h bm + bet Outer[Times, qj, qj];
                           h # + SparseArray[{j -> bet}, n, zero]) &, u, 
                          SameTest -> (SquaredEuclideanDistance[##] <= tol &)].pts;
           Ellipsoid[c, d (Drop[bm, -1, -1] - Outer[Times, c, c])]]

The routine returns an Ellipsoid[] object, suitable for display or any other further manipulations (e.g. one can apply an eigendecomposition to the matrix in the second argument of the Ellipsoid[] object to get the axes).

BlockRandom[SeedRandom[42, Method -> "MersenneTwister"]; (* for reproducibility *)
            pts = RandomVariate[MultinormalDistribution[RandomReal[{-1, 1}, {2}], 
                                With[{m = RandomReal[1, {2, 2}]}, m.Transpose[m]]], 1024];]

Graphics[{Red, mve[pts], Blue, MeshPrimitives[ConvexHullMesh[pts], 2], Yellow, Point[pts]}]

minimum-volume ellipse and 2D convex hull

BlockRandom[SeedRandom[42, Method -> "Rule30CA"]; 
            pts = RandomVariate[MultinormalDistribution[RandomReal[{-1, 1}, {3}], 
                                With[{m = RandomReal[1, {3, 3}]}, Transpose[m].m]], 1024];]

Graphics3D[{Opacity[1/2, Red], mve[pts], Blue, MeshPrimitives[ConvexHullMesh[pts], 2]}]

minimum-volume ellipsoid and 3D convex hull

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.