"Very" high dimensional convex hulls

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

[Oliver]

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.

[Oliver]

One idea that I've come across is using the separating axis theorem. My naive approach would be: for each test point (1) project it onto all N*(N-1)/2 vectors given by the differences of the N potential points defining the hull; (2) test if the projection of the test point overlaps with the projection of the N hull points; (3) if it does check the next vector, else it is outside of the hull so break. The only problem with this is that it seems like for each of the points that fall inside the hull I will need to test all N*(N-1)/2 vectors (potential separating axes). That means that if I want 1e3 points inside the hull I will have to do at least 1e3*N*(N-1)/2 tests, so with N~1000, that makes for O(5e8) tests. If the projection and testing procedure is cheap enough maybe this can be practical, but I'm not sure yet how to do this efficiently. Ideas? Or am I barking up the wrong tree?

[Oliver]

Maybe there's an intelligent way to pick which vectors to test (all you really need to test are the normals to the facets of the hull, but since I don't want to construct the hull explicitly...)?

[Oliver]

Edit: to use the separating axis theorem you have to test the normals to the facets of the hull, but since I don't want to construct the hull, I don't have the facets. You can't project onto the edges as I suggested above, at least it will incorrectly classify points that are close to the hull but outside of it. I'm stuck again...

[Oliver]

I found a post here http://stackoverflow.com/a/14142517/1181402 that is relevant. The answer about finding a separating hyperplane using SVM seems promising, but I don't understand how to use it, or how it could be any faster than what Matt J proposed below.

Here's what was said "The problem of finding a separating hyperplane can be transformed to a convex quadratic programming problem (see: SVM). I did this in python using cvxopt with less then 170 lines of code (including I/O)."

But how is the problem transformed into a convex QP problem?

[Matt J]

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?

[Oliver]

No, no, that was just what someone else posted at stackoverflow (I was quoting their post). I have no idea how they did it. I'm still just as stuck as ever. It just seemed that the method they suggested was promising perhaps.

Answer 1

[Matt J]

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.

[Oliver]

To be honest, I'm still trying to figure out what a good value of N is. It needs to be at least 1e3, and could be as high as 1e6. But somewhere between 1e3 and 1e4 is probably about right.

[Oliver]

I tried your idea on 3D data and it works like a charm! However, when I applied it to my actual 34-dimensional data lsqlin was too slow (and the number of points one has to try before finding a point inside of the convex hull is likely enormous in 34 dimensions). In my application it took ~10sec per point, and in the test case I ran 0/100 of the points were successful

I'm not terribly familiar with the optimization toolbox or QP in general, is there a chance that another solver or some different options from the defaults would be faster?

If more details would be useful, in my application V.' is a 34-by-1000 matrix, so w is 1000-by-1 (obviously a lot of work for a solver).

[Matt J]

How are you initializing lsqlin? I would recommend you find the closest V(i,:) to y,

 [~,imin]=min(sum(bsxfun(@minus,V.',y).^2))

Then set w(imin)=1 and all other w(i) to 0.

[Oliver]

Unfortunately it took just as long with this modification. Here's the code, with fabricated data points in the columns of "c_temp":

nsamples = 1e2;
npts = 1e3;
ndims = 34;
c_temp = rand(ndims,npts);

%---set up QP problem---%
Aeq = ones(1,npts);
beq = 1;
lb = zeros(npts,1);
ndims = size(c_temp,1);
options.LargeScale = 'off'; %system is underdetermined, so a warning is thrown and it defaults to medium scale anyways

temp_samples = zeros(ndims,0);
N = 0;
while N < nsamples
    %---sample hypercube---%
    samples = rand(ndims,nsamples);

      %---keep the ones inside of the hull without explicitly computing it---%
      d = zeros(nsamples,1);
      for i = 1:nsamples
          y = samples(:,i);
          [~,imin] = min(sum(bsxfun(@minus,c_temp,y).^2));
          x0 = zeros(npts,1);
          x0(imin) = 1;
          tic
          x = lsqlin(c_temp,y,[],[],Aeq,beq,lb,[],x0,options);
          toc
          d(i) = norm(c_temp*x-y).^2; %distance from closest point inside of the hull, C*x, to the test point, y.
      end
      in = (d == 0);
      temp_samples = [temp_samples,samples(:,in)];

      %---how many samples have we kept so far---%
      N = size(temp_samples,2);
  end
  temp_samples = temp_samples(:,1:nsamples);

Answer 2

[Matt J]

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.

[Oliver]

The first constraint says that all of the projections of vectors from a test point to each hull point need to be non negative (explicitly you've said the negative of these vectors need to be non positive). I'm on board so far.

But I'm unsure what the meaning of the second constraint "f*x<0" means? What does your variable "f" represent physically?

[Matt J]

But I'm unsure what the meaning of the second constraint "f*x<0" means?

It's mainly to ensure that you find a solution x which is not zero. The point x=0 is always a trivial solution to A*x<=0 and does not tell you anything about separability.

By choosing f=sum(A,1), it ensures that not only will x satisfy A*x<=0, but also that A(i,:)*x<0 (strict inequality) for at least one i, which should always be satisfied by any separating hyperplane as long as the convex hull is solid in R^D.

I'm assuming you know your convex hull to be solid, i.e., not confined to a hyperplane in R^D. Otherwise, the probability of landing in it when sampling from the larger hypercube will be zero.

[Oliver]

Ahhh, I see. Brilliant!

Yes, I believe it is solid (though your comment made me double check, am I correct in assuming that if the matrix containing the hull points has full rank then it is solid?).

[Oliver]

Your suggestion works brilliantly in 3D (validated it against John D'Errico's "inhull" function from the FEX), and the speed appears to be adequate to handle 34D. I tested 1e6 points in 34D in about 10 min using a parfor loop and 4 cores. However, 0 out of the 1e6 points was found to be contained in the hull. I suppose it is possible that the convex hull of the points is very small, or very close to being a plane or line (nearly coplanar or colinear points), but is there any way to test for that? (something like being "almost" not full rank?)

[Matt J]

am I correct in assuming that if the matrix containing the hull points has full rank then it is solid?

Not exactly. You need to test that the rank of this matrix is D

 Q=bsxfun(@minus,V(2:end,:), V(1,:));

The point being subtracted V(1,:) is arbitrary. It can be any of the hull points. BTW, there are occasionally issues with Matlab's RANK function, so you're taking your chances by using it. It uses funny thresholding strategies. I use a different rank calculation approach in VERT2LCON.

[Matt J]

I suppose it is possible that the convex hull of the points is very small, or very close to being a plane or line (nearly coplanar or colinear points), but is there any way to test for that?

I think you need to work in a transformed space where your hypercube fits as snugly as possible around the hull. How about this. Manipulate your hull samples as follows

   X=bsxfun(@minus,V(2:end,:), V(1,:)).';

   [r,idx,T]=lindep(X); %a sub-function in lcon2vert.m

   if r<D, error 'Convex hull is not solid'; end

   X=(T\X).';

   cube_bounds=[min(X).',max(X).'];

Now, apply the sampling method to the transformed hull data, X, and using the hypercube bounded by the data in cube_bounds. You can untransform your samples at the end if you want,

 samples=bsxfun(@plus, samples*T.', V(1,:));

[Oliver]

Why do you remove V(1,:), rather than just subtracting it off? I.e. why use:

X=bsxfun(@minus,V(2:end,:), V(1,:)).';

instead of:

X=bsxfun(@minus,V, V(1,:)).';

Presumably, one point won't change the hull very much...but it might if you happen to pick the right one.

[Matt J]

You're right. You should leave V(1,:) in there.

I'm starting to have doubts that this hypercube idea is going to work, thoug :-(

The problem is that the ratio between the volume of a polyhedron and its bounding hypercube can get exponentially smaller in higher dimensions. You could need a vast amount of samples to have a probability of landing in the polyhedron in 34D. You can test this by trying a solid unit simplex

 V=eye(D+1,D)

for larger and larger D and with the unit hypercube bounding it. I'll bet you'll see your percentage of hits fall off exponentially fast. Here's more evidence

    s=@(D) eye(D+1,D);

    v=1;
    for D=2:20;

        [~,v(D)]=convhulln(s(D));

    end

     plot(v)

[Oliver]

I'm beginning to suspect the same thing, unfortunately. I tried it with and without the transformation, and both result in 0 out of 1e6 trials landing in the convex hull. Hmmm...

[Oliver]

There is one strange thing though. Just for kicks I tried to see if the method would work on points that I know apriori are inside the hull, 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? I tried the same thing on 3D and found that, while most of the points that define the hull were correctly identified as belonging to it, not all of them were!

%---generate random points---%
ndims = 3;
npts = 1000;
V = rand(npts,ndims);

%---fake samples---%
nsamples = npts;
samples = V;

%---LP algorithm for finding separating hyperplane axis---%
in = false(nsamples,1);
options = optimset('algorithm','active-set','largescale','off','display','off');
C = V;
x0 = zeros(ndims,1);
for i = 1:nsamples
    y = samples(i,:);
    A = -bsxfun(@minus,y,C);
    b = zeros(npts,1);
    f = sum(A,1).';
    x = linprog(f,A,b,[],[],[],[],x0,options);
    if all(x == x0)
        in(i) = true;
    end
end

[Matt J]

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.

Answer 3

[Matt J]

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.

[Oliver]

Hi Matt. No, taking random convex combinations does not work. The samples are all concentrated tightly around the "centroid" of the convex hull. I think fundamentally that the reason it doesn't work is because the convex hull is a convex polytope, but not a simplex. Here's some code to show what I mean

%---generate random points---%
V = rand(20,3);

%---get convex hull---%
k = convhulln(V);
ktess = k;
k = unique(k(:));

%---sample random convex combinations---%
nsamples = 1e3;
nverts = numel(k);
u = rand(nsamples,nverts); %uniform samples
p = bsxfun(@rdivide,u,sum(u,2));

%---map to xyz coordinates---%
samples = p*V(k,:);

%---plot---%
figure
h = trisurf(ktess,V(:,1),V(:,2),V(:,3));
set(h,'edgecolor',[0.9 0 0.1],'edgealpha',0.5,'facecolor','r','facealpha',0.5)
hold on
plot3(samples(:,1),samples(:,2),samples(:,3),'k.')
axis equal tight vis3d
box off
grid off
xlabel('x')
ylabel('y')
zlabel('z')

[Matt J]

Can you elaborate on why it is important that the sampling be uniformly distributed? In R^30, it would take on the order of 1e30 samples to uniformly cover a solid convex hull. You are taking much fewer samples than that.

[Oliver]

The application has to do with mapping from one space to another. Each of the dimensions of my convex hull corresponds to the coefficient of a generalized fourier series. Call these coefficients the "c" coefficients. In reality it is infinite dimensional, but I'm looking at just the first ~30 independent, non-degenerate, lowest order terms. The allowable values of the coefficients define a convex hull in the infinite dimensional fourier space (and all orthogonal projections are also convex, which is why the first 3 or 30 dimensions taken alone also define a convex hull).

I have another set of fourier coefficients, call them the "t" coefficients, which, under certain assumptions, can be expressed as functions, of the "c" coefficients. I am trying to look at the unrestricted convex hull of the "t" coefficients in their fourier space, and compare it to the convex hull of the "t" coefficients that are computed from the "c" coefficients. This convex hull is smaller and fully contained within the unrestricted hull. I want to quantify how much smaller and what it "looks" like.

My idea for doing this was sampling uniformly over the convex hull of the "c" coefficients in their fourier space, then calculating the corresponding "t" coefficients for each of these points, then compute the convex hull and compare. (For comparison I'm looking at low dimensional orthogonal projections, but the mapping from c -> t requires many of the c coefficients, i.e. a single t coefficient may depend on the values of 30 c coefficients, which is why I need to sample over a 30 dimensional projection of the full infinite dimensional c hull).

So, to answer your question, in reality I don't need statistical uniformity for this to work. I just need to make sure that I've sampled from the entire c-hull. Taking random convex combinations does not accomplish this because all of the samples are near the center so I don't get anything near the boundary. A uniform distribution would just be the most efficient way to fill the hull with points, but any other kind of distribution that has some samples from all over the hull would be acceptable.

[Matt J]

This convex hull is smaller and fully contained within the unrestricted hull. I want to quantify how much smaller and what it "looks" like.

Even if you could find the t-space convex hull, how would you quantify its size and shape? Is it of lower dimension than the c-space hull? CONVHULLN could normally be used to measures its volume, but if it is high-dimensional like the c-space hull, then you will have the same problem running it through CONVHULLN as in c-space.

On the other hand, if the set in t-space lies in a significantly lower-dimensional manifold, then the best approach might be to sample t-space instead. You could solve an inverse problem to see if a sample point in t-space maps from the c-space convex hull and discard it otherwise. To evaluate this approach, it would help to see the form of the c-->t mapping. Is it linear, for example?

[Oliver]

The mapping is non-linear, and is of the form:

t(m) = sum(alpha(i,j,k).*c(i).*c(j).*c(k))

Where t is a vector of the t-coefficients, and c is a vector of the c-coefficients, and alpha is a constant dependent upon the indices i,j,k. In other words, t is a linear combination of products of 3 c coefficients.

It is true, quantifying the size and shape is difficult. I was hoping that one might get a qualitative feel for their relative sizes by looking at the ratios of the volumes of the convex hulls of t and c, in various 3D orthogonal projections.

From your suggestions, and our discussion it appears that the best solution to this whole issue is to find a way to avoid dealing with 30+ dimensions in the first place. I think I've found a way to do that. The first 5 t coefficients depend only on 6 c coefficients, so I can sample from an orthogonal projection of c-space that is of dimension 6, and follow the procedure above. This makes explicit construction of the convex hull feasible, and allows me to use your vert2lcon function and cprnd to sample without doing rejection sampling. 6D still takes a while (2 days to do 1e6 samples on my machine), but its definitely doable.

Thank you so much for your time and suggestions, they really helped me attack the problem from different angles, which ultimately I think was what I needed to do. I wish I could just accept our whole conversation as an answer, but being unable to, I just selected one of your (many) helpful answers.

[Matt J]

OK, glad you found a solution!