It seems like the ideal solution would be one in which the convex hull did not explicitly need to be computed. One idea is to generate samples uniformly over a hypercube and throw out the ones that are outside of the convex hull without actually computing the hull itself. But, since I don't know which points are on the hull apriori without constructing it I'm not sure how to do this either.
"Very" high dimensional convex hulls
19 views (last 30 days)
Show older comments
I am trying to generate points that are uniformly distributed over a convex hull in a high dimensional space.
What I have as input is a set of N points in D dimensions. I want to sample uniformly over the convex hull of this set points. The best solution that I have found so far is by using Matt Jacobson's "vert2lcon" (file 30892 on the FEX) to convert from vertex representation to a linear constrain representation (i.e. inequalities defining the facets of the hull), and then to use Tim Benham's "cprand" (file 34208 on the FEX) to sample uniformly over the hull.
This works great for a few dimensions, however, the hull's that I am interested in are very high dimensional, somewhere between 30-60 dimensions usually. When I try to use vert2lcon on all D~30 dimensions MATLAB crashes. The reason is that vert2lcon constructs the convex hull of my points by calling "convhulln", which in turn calls "Qhull", and the crash report says that the error is happening during the mex function for "Qhull". My guess is that it can't handle that many dimensions.
Any suggestions?
I have thought of trying to take a subset of the dimensions at a time, but I'm not sure that that would work (thought about the example of a sphere, if you take the one dimensional orthogonal projections, and compute their convex hulls you get three lines parallel to the coordinate axes, the 3D convex hull of these convex hulls is going to be an octahedron so it definitely doesn't reproduce the sphere).
7 Comments
Matt J
on 30 May 2013
I did this in python using cvxopt with less then 170 lines of code (including I/O)
Does this mean you already have a solution that you like? What is the status of the problem?
Accepted Answer
Matt J
on 29 May 2013
Edited: Matt J
on 29 May 2013
One idea is to generate samples uniformly over a hypercube and throw out the ones that are outside of the convex hull without actually computing the hull itself. But, since I don't know which points are on the hull apriori without constructing it I'm not sure how to do this either.
You could do that by solving the QP (e.g., using lsqlin)
min d(w) = norm(V.'*w-y).^2
s.t. sum(w)=1
w>=0
where y is the point you want to test for inclusion in the convex hull. If the optimal d(w) is zero, then y is in the hull
You haven't mentioned how big N is so I can't know in advance how practical this would be for you.
More Answers (2)
Matt J
on 30 May 2013
Here's another possibility based on linear programming and separating plane ideas. Define
A=bsxfun(@minus, V,y(:));
f=sum(A,1);
A separating hyperplane exists if there is a vector, x satisfying
A*x<=0
f*x<0
which can be found by running a linear program solver with objective f*x and constraints A*x<=0 until a feasible solution satisfying f*x<0 is reached. This would probably require just one iteration of LINPROG starting at the initial feasible point x0=0 and using the 'active set' algorithm option.
11 Comments
Matt J
on 1 Jun 2013
Edited: Matt J
on 3 Jun 2013
in particular I used the points that define the hull as my samples, and in 34D it said that none of the points that define the hull are inside of it! What is going on?
Only points strictly in the interior of the polyhedron will fail the separation criterion f*x<0.
Also, one can only ever hope to distinguish between points that are either strictly in the interior of the polyhedron or strictly in the exterior. Boundary points can never be reliably tested by any method because floating point errors can make them look either slightly inside or slightly outside the polyhedron.
BTW, you should also be setting MaxIter=1. You probably also need some upper and lower bounds, e.g. lb=-ones(D,1), ub=-lb.
Matt J
on 29 May 2013
Taking random convex combinations of the vertices V(i,:) like in the following definitely won't work?
W=rand(M,N);
W=bsxfun(@rdivide,W,sum(W,2));
samples=W*V; %rows of V are the vertices.
6 Comments
See Also
Categories
Find more on Bounding Regions in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!