File Exchange

Analyze N-dimensional Polyhedra in terms of Vertices or (In)Equalities

version 1.9.0.0 (14.8 KB) by
Express polyhedra by vertices or (in)equalities. Also, find intersections and unions.

24 Downloads

Updated 04 Oct 2017

View Version History

View License

This submission contains a set of files for analyzing N-dimensional polyhedra. It is intended for fairly low dimensions N -- basically low enough so that vertex and facet enumeration using MATLAB's convhulln() command is tractable. For now, it is also limited to bounded polyhedra (i.e., polytopes). A bounded polyhedron can be represented either as the convex hull of a finite set of vertices V(i,:) or by using a combination of linear constraint equalities and inequalities,
A*x<=b,
Aeq*x=beq

Here, A and Aeq are MxN and PxN matrices while b and beq are Mx1 and Px1 column vectors, respectively. The (in)equality representation expresses the polyhedron as the intersection of two regions. One region is a solid N-dimensional shape, described by the inequalities, while the other is a possibly lower-dimensional sub-space, described by the equalities. The screenshot above illustrates this, showing how a triangle in 3D can be represented as the intersection of a tetrahedron (a solid shape in R^3) and a plane.
The package contains tools for converting between the two representations (see VERT2LCON and LCON2VERT) as well as for taking intersections and unions of polyhedra in either form (see intersectionHull and unionHull).

The package was inspired by Michael Kleder's vert2con and con2vert functions, which were limited to N-dimensional polyhedra possessing non-zero volume in R^N.Thus, for example, they could not handle a triangle floating in R^3 as depicted in the above screenshot. Although a triangle has non-zero volume (i.e., area) in R^2 it has zero volume in R^3. The extension made in this package covers general bounded polyhedra. NOTE: however, when using linear constraint data A, b, Aeq, beq to represent a given polyhedron, it is important that the inequalities A*x<=b still be chosen to correspond with a region of non-zero N-dimensional volume. Zero-volume polyhedra are captured by adding equality constraints Aeq*x=beq.

EXAMPLE:

Consider the 3D polyhedron defined by x+y+z=1, x>=0, y>=0, z>=0. Equivalent constraint in/equalities can be obtained from the known vertices by doing,

>> [A,b,Aeq,beq]=vert2lcon(eye(3))

A =

0.4082 -0.8165 0.4082
0.4082 0.4082 -0.8165
-0.8165 0.4082 0.4082

b =

0.4082
0.4082
0.4082

Aeq =

0.5774 0.5774 0.5774

beq =

0.5774

Conversely, the vertices can be obtained by doing,

>> V=lcon2vert(A,b,Aeq,beq)

V =

1.0000 0.0000 0.0000
0.0000 0.0000 1.0000
-0.0000 1.0000 0.0000

When an interior point of the polyhedron is known a priori, one can instead use QLCON2VERT, a quicker version of LCON2VERT which uses the known point to get around some computationally intensive steps.

Cite As

Matt J (2020). Analyze N-dimensional Polyhedra in terms of Vertices or (In)Equalities (https://www.mathworks.com/matlabcentral/fileexchange/30892-analyze-n-dimensional-polyhedra-in-terms-of-vertices-or-in-equalities), MATLAB Central File Exchange. Retrieved .

Comments and Ratings (111)

Zefan Lin

Hi Matt,
I learned how to cancel this printing. By setting 'Display' to 'Off'. Thank you!

Zefan Lin

Hi Matt,
Thank you! It is really helpful. But I would like to know how to cancel the printing of the following prompt:

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

<stopping criteria details>

They always pop up when I run lsqlin().

Zefan Lin

Matt J

Hi Zefan,
You would use lsqlin() to find the closest point, but you could use vert2lcon to calculate the A,b,Aeq,beq inputs that lsqlin requires.

Zefan Lin

Hi Matt,
Does this code allow one to calculate the minimum distance between a convex polyhedron and a point? The convex polyhedron is composed of five vertices (xi ∈ R^4, i = 1, 2, 3, 4, 5) in a 4D space. The point is outside the convex polyhedron.

Matt J

@SK,
It sure does.

LC

Hi. Does this code allow one to calculate the intersection points of a polyhedron having a set of vertices V in 3D and a straight line segment with end points P0 and P1?

AHMED Al-Ajeli

Hi Matt,
Thank you very much for your quick reply.
Yes, it might be because of the large number of dimensions.

Regards,
Ahmed

Matt J

Hi Ahmed,
I think you may have posted the wrong data because A and Aeq have different numbers of columns. Also, it doesn't really make sense that you have more equality constraints than you do unknowns. In any case, I am not too optimistic that you could apply lcon2vert to problem dimensions this large.

AHMED Al-Ajeli

Hi Matt,
Thanks a lot for the nice piece of code.
I used V=lcon2vert(A,b,Aeq,beq) on the following materices:
Aeq = [0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0;
0 0 -1 0 -1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0;
1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0;
-1 1 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0;
0 -1 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0;
0 0 0 -1 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0;
0 0 0 0 -1 -1 0 1 0 -1 0 0 0 0 0 0 1 0 0 0 0 0;
0 0 0 0 0 1 0 -1 0 1 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 1 -1 0 0 -1 0 0 0 0 0 0 0 0 1 0 0 0;
0 0 0 0 0 -1 1 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0;
0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 0 0 1 0;
0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 0 0 1;
0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
-1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1]

beq =

0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1

A = [ -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1]

b =

0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

Unfortunately, I got this error message "qhull internal warning (main): did not free 176 bytes of long memory (1 pieces)
Something's wrong. We should have found a recession vector (bb<0)."

Could you please help me for this?
Many thanks...
Regards,
Ahmed

Hassan Dehghani

Matt J

@H. Tong,
The FEX package was written with the deliberate intention not to require the Optimization Toolbox. I believe most linear programming packages use the ellipsoid method to find the initial point. I've been meaning to look into finding an implementation of it, but it hasn't gotten to the top of my TODO list yet. If you have the Optimization Toolbox, you could use linprog to get an initial point x0 and then use qlcon2vert - the quicker version of lcon2vert that lets you supply your own initial point.

H. Tong

Is it possible to just use linear programming to find the initial point in con2vert?

Matt J

Thanks, Benjamin. I'm not sure, though, what would be giving you trouble concatenating with []. It should work fine with matrices of any dimension, e.g.,

>> [ [] ; [1,2,3;4,5,6] ]

ans =

1 2 3
4 5 6

Benjamin Bernard

Thank you! Very useful piece of code.

Minor suggestion for improvement: if the intersection is empty, return an empty set of the right dimensions, i.e., return V=zeros(0, size(A, 2)). This would make it a little bit more user friendly if, for example, one wants to append those vertices to some existing set of vertices.

Shachar

Matt J

Hi Aaron,
It sounds like the hulls you have input simply do not intersect, or at least their intersection is not numerically stable. But if you are sure that is not the case, please present the example for study.

Aaron Strohbusch

Hey Matt,

Awesome code. There seems to be an issue with lcon2vert when calling the intersectionHull function, or maybe its not a problem idk. In the Initializer2 sub-function, it requires a large amount of iteration. I've tried playing around with the max iterations and I keep running into a problem where it will return empty vert and lcon fields, or it requires tons of iteration, making the function run very slowly. It also returns empty vert and lcon fields when the intersectionHull is empty, meaning the two hulls are not intersecting at all, which means the function is throwing false negatives. I was really excited to use this code, so hopefully there is a fix.

Looking forward to hearing from you,
Aaron

Luis Almeida

Hi Matt,
Yes, I was talking about vert2lcon, sorry not to made it clear.
That makes sense. Thank you for helping me solve my problem.

Matt J

Hi Luis,
I'm guessing that you are talking about vert2lcon. If so, yes, the workings of the code is largely as you describe. To solve the "distance-to-boundary" problem, just orthogonally project your point onto the planes passing through all the facets. Then check which projections satisfy all the constraints and which among those has the minimum distance from the given interior point.

Luis Almeida

So I have a similar problem to this one:
https://www.mathworks.com/matlabcentral/answers/107595-how-can-i-find-the-minimum-distance-from-convex-boundary
But I want to compute the closest point from the boundary to a point inside of the convex hull (to the origin).
With your answer, the point that I got is the origin itself and I was expecting one point on convex hull boundary.

Also I would love to understand better your code, so roughly speaking what you do is: for each three points that define a facet, you compute the normal to that facet, n = (a,b,c) and then you make an plane inequation. Let's say it's a*x+b*y+c*z <= d, then you put [a b c] in each row of A and "d" in "b".
And all this inequalities describe the convex hull (boundary and the inside).
Am I thinking right?

Thanks a lot for your time

Matt J

Hi Luis,
I don't understand the question, I'm afraid.

Luis Almeida

Hi Matt!
How do I enforce constraints for points to lie only on the boundary of the convex hull?
Thanks in advance

Devin

rubindan

Great work! manage to work with high complexity polyhedra which MPT3 cannot handle.

Raphaël Boudreault

@Matt J
Thank you for your quick answer. I'll try something else to avoid this situation in my algorithm.

Matt J

Hi Raphael,
Unfortunately, there is no numerical algorithm that will correctly find the vertex in your example. The reason is that the solutions to A*x<=b, Aeq=beq in your example are not stable under small perturbations of the equation data. For example, if I were to perturb b(3) to -epsilon, where epsilon>0 is any arbitrarily small positive number, the equations no longer have any solutions. A symbolic solver would possibly be able to identify the solution, as you have done analytically, but I'm afraid I don't know where you can find code for that.

Raphaël Boudreault

Hi Matt J,
I may be doing something wrong, but in one iteration of my algorithm, I need to find the vertices of the polyhedron given by :
A =
-1 0 0
0 -1 0
0 0 -1
b =
0
0
0
Aeq =
-1 1 -1
0 -1.5 0
0 0 0.5
beq =
-1
0
0
By hand, I find easily that there is a unique solution to these equations, i.e. (1, 0, 0). But with lcon2vert, I get an empty vertices set!
Thanks,
Raphael

Kavit Tolia

Matt J

Hi Shubham,

That is a numerically unstable approach because small floating point errors in your plane equation can cause large changes in the plane's intersection with the tetrahedron. The better approach is to first find all the vertices of the tetrahedron. Then, plug them into the plane equation to see which ones lie in that plane.

SHUBHAM SHUKLA

Thanks for your solution.
I wish to use this code for finding the vertices of the tetrahedron (as given in the above example) by giving the equality constraints Aeq*x=beq as one of the planes forming the tetrahedron. I am assuming it should give the vertex lying on that plane of the tetrahedron.
However, its implementation is not giving the entire set of vertices.

Can anyone suggest some solution for this application?

Taniya Seth

@Matt

Thank you so much for your response! I never noticed that message! Thanks a lot.

Matt J

@Taniya,

Your polyhedron is unbounded, and therefore outside the scope of the tool. Did you see the following in the error message?

Non-bounding constraints detected. (Consider box constraints on variables.)

If so, that is one option to consider.

Taniya Seth

Hi, I may be going wrong with something really trivial. But I have the following values
A = [10000 20000; 100 75];
b = [100000; 500];
When I run the command V = lcon2vert(A,b,[],[])
I get the following:
Error in lcon2vert (line 254)
Zt=con2vert(AAA,bbb,TOL,checkbounds);
Can anybody suggest a get-around? I'm sorry I may sound really dumb.

Emanuele

Matt J

The bug has been fixed.

Matt J

Just as an FYI, I found a small bug in lcon2vert which is on my TODO list to fix. It only affects the paticular case when the equality constraints have only 1 solution. The code neglects to check if the inequality constraints are satisified by that unique point.

Matt J

Hi Matthew (Thirkettle),

Thanks a lot for your feedback and your kind rating. The issue is not the robustness of the code, but rather of your description of the polytope. Due to finite precision issues, no computer can be relied upon to correctly recognize a polytope of zero volume from a set of inequalities, . A similar example to yours is the following. Even though the polytope should contain the single point x=[1/100/pi,1]', my computer determines that it does not. Whether you see the same on your machine may depend on CPU differences, but the point is that it has nothing to do with lcon2vert.

>> A=[100*pi,-1;-pi,+.01]; b=[0;0]; x=[1/100/pi,1]';
>> all(A*x<=b)

ans =

logical

0

As the file description above says, "it is important that the inequalities A*x<=b still be chosen to correspond with a region of non-zero N-dimensional volume." This is the reason why.

If your polytope has non-zero volume, but is very thin, you could try passing a smaller TOL argument to lcon2vert

[V,nr,nre]=lcon2vert(A,b,Aeq,beq,TOL)

but my recommendation is that you express some of your variables in different units, so that it is not so thin. If it has zero volume, you must find a way to express it as the intersection of a non-zero volume and a set of equalities.

Matthew Thirkettle

Hi Matt:

Great code. In my problem, I am trying to find the vertices of a polytope defined by inequality that can potentially become thin. Unfortunately, your code is not robust to this. Do you have any suggestions on how to solve this? In particular, consider the inequality constraints:

A = [1 0 ; 0 1; -1 0 ; 0 -1], b = [1;1;-1;-1].

There is only one point in the polytope defined by {x : Ax <= b}, namely x = [1;1]. Your code returns V = [] if we set V = lcon2vert(A,b) and returns V = [1;1] if we set V = lcon2vert([],[],A,b).

Now suppose we set b(3) = -1 + e and b(4) = -1 + e for e small. This defines a tiny cube with vertices [1,1], [1 - e, 1], [1, 1-e], [1-e,1-e]. Provided e is sufficiently large (say > 1e-6), V = lcon2vert(A,b) returns the correct answer. However, for small e we get the same problem as above.

Matt J

Mathieu,

My suggestion would have been | Facet * vertex - bound | <=epsilon where epsilon is some tolerance parameter. Note that vert2lcon returns the facets normalized to unit vectors so | Facet * vertex - bound | will represent a Euclidean distance. If the tolerance is too hard to choose to get the correct vertices, that would probably indicate something degenerate about your polyhedron.

Mathieu

Hi Matt !

I have a question: once I find the facts of a polyhedron, how can I know which vertices of the polyhedron lie on a particular facet ?

I tried to check if | Facet * vertex - bound | == 0, but I was wondering if there was a better way, as doing this not always gives exactly 0 for a vertex that lies on the facet I want. Maybe by combining it with convhulln ?

Thanks,
Mathieu

Mathieu

Hi Matt !

I have a question: once I find the facts of a polyhedron, how can I know which vertices of the polyhedron lie on a particular facet ?

I tried to check if | Facet * vertex - bound | == 0, but I was wondering if there was a better way, as doing this not always gives exactly 0 for a vertex that lies on the facet I want.

Thanks,
Mathieu

zhang baoqiang

ok, thanks again!

Matt J

Thanks, zhang. The code isn't based on any journal reference. It was largely self-derived with inspiration from Michael Kleder's code.

http://www.mathworks.com/matlabcentral/fileexchange/7894-con2vert-constraints-to-vertices

My best advice is just to cite the link to this FEX page,

http://www.mathworks.com/matlabcentral/fileexchange/30892-analyze-n-dimensional-polyhedra-in-terms-of-vertices-or--in-equalities

zhang baoqiang

Thank you very much for your quick replies! And I was sorry for my careless mistake...You supply a great code!
I want to quote your code in my studies, could you give me some references about it?

Matt J

Hi again Zhang,

I looked more closely at your data and ran it through lcon2vert myself. I do not get the error message you posted. Rather, lcon2vert reports that the problem is infeasible,

>> lcon2vert(A,b,Aeq,beq)

ans =

[]

It is easy to see from your equality constraints, that this is the correct result. Your first equality constraint is x1+x2=9 whereas the 3rd equality constraint contradicts the first with x1+x2=10.

Matt J

Hi zhang,

lcon2vert has determined that your A,b,Aeq,beq data violates the assumption that your polyhedron is bounded. You might try adding some upper bounds to the requirement x>=0 that you already have.

zhang baoqiang

HI Matt J,
thanks very much for great work you've done!

however, something is wrong when I give A,b,Aeq,beq like this:

A =

-1 0 0 0 0 0 0 0
0 -1 0 0 0 0 0 0
0 0 -1 0 0 0 0 0
0 0 0 -1 0 0 0 0
0 0 0 0 -1 0 0 0
0 0 0 0 0 -1 0 0
0 0 0 0 0 0 -1 0
0 0 0 0 0 0 0 -1

b =

0
0
0
0
0
0
0
0

Aeq =

1 1 0 0 0 0 0 0
1 1 1 0 0 0 0 0
1 1 0 0 0 0 0 0
0 0 1 0 0 0 1 1
0 0 1 0 0 0 1 1
0 0 0 1 1 0 0 0
0 0 0 0 0 1 1 1
0 0 0 0 0 1 1 0
0 1 1 0 0 0 0 1
0 0 1 0 0 0 0 1

beq =

9
10
10
9
10
8
9
9
10
9

Error using lcon2vert>con2vert (line 308)
Non-bounding constraints detected. (Consider box constraints on variables.)

Error in lcon2vert (line 234)
[Zt,nr]=con2vert(AAA,bbb,TOL,checkbounds);

Error in test (line 8)
[V,nr,nre] = lcon2vert(A,b);

I can't understand the error, is there something wrong in my data? Could you give me some guidance?

Thank you in advance

Matt J

Hi mirshad,

I'm afraid I don't know of a systematic way to check intersection that doesn't involve actually computing the intersection.

However, one way to save computation might be to first test the vertices of each polyhedron to see if a vertex lies inside the other polyhedron. This is a sufficient, but not necessary condition for intersection. However, it might hit often enough to make it a worthwhile test. Similary, you could check if the vertices of one polyhedron all lie on the wrong side of one of the bounding planes of the other polyhedron. This is a sufficient, but not necessary condition for non-intersection of the two polyhedra.

mirshad

Hi Matt,
I only want to check if two convex polyhedrons that have created by delaunayTriangulation command intersect each other or not and do not want to find polyhedron resulting from their intersection. Would you give me some guidance?

nilay kant

Thanks a lot. Amazing Work!!

Matt J

Hi nilay,

The basic idea is that if polytope 1 is given by inequalities A1*x<=b1 and polytope 2 is given by A2*x<=b2, then their intersection is given by the combined system [A1;A2]*x<=[b1;b2]. In other words, once you have the polytopes in algebraic form, the intersection is obtained just be concatenating inequalities.

nilay kant

Hii Matt,
Thanks a lot for developing this code. I have a small query. I still cant understand how you compute the intersection of two convex polytope?
Thanx in advance

Matt J

@nilay,
Yes, it sounds like intersectionHull would do what you want. Your "2D slice" is presumably given by 2 linear equalities in 4 variables. You can use those equalities to specify this region for intersectionHull.

nilay kant

Hi,
I need to take a 2-D slice of conved hull formed by scattered points in 4-D. Can this be done through this?
Thanx

Jon Dattorro

finds vertices of a permutation polyhedron.
nice job.

Joerg

Great, thank you very much!

Nidal Kochrad

Hassan Naseri

Really useful and straightforward tool. Thanks indeed.

Matt J

@squirrel,

The vert2lcon parent function takes the N-dimensional input vertices and determines the minimal dimension, D<=N, of the polyhedron and the D-dimensional affine set in which it lies. It does this using QR analysis. It then transforms the vertex data into D-dimensional data and reformulates the problem in R^D where the polyhedron is solid. The R^D data is then passed to the vert2con subfunction. Because the transformed polyhedron is solid, you can use convhulln on it and find the facets of the polyhedron's hull. It is then a simple matter to fit hyperplanes to all the facets and then transform everything back to the original space R^N.

squirrel

Dear Matt J,

thank you very much for great work you've done!

Could you PLEASE describe mathematical rationale, on which functions (vert2lcon, vert2con) from vert2lcon.m file are based on?
Could you provide a short description to the main steps of these functions?

Matt J

Dear Christophe,

No, changing tolerances won't help. Your problem is ill-posed. Because the ideal triangles only touch at a single point, small perturbations of the input data can make the triangles intersect at multiple points or not at all.

This is why, in the Description section above, I mentioned "It is to be emphasized that A,b must define a solid region". Your intersection point is not solid in R^2, so the inequalities [A1;A2],[b1;b2] are not a legal way of expressing it.

Christophe Lauwerys

Dear Matt,

tried calculating the intersection between a two triangles touching each other in a point but got following error.

Assume can be solved by setting tolerances correct?

Thanks

Christophe

>> [A1,b1,Aeq1,beq1]=vert2lcon([2,0;0,2;0,0]);
>> [A2,b2,Aeq2,beq2]=vert2lcon([1,1;2,1;1,2]);
>> V = lcon2vert([A1;A2],[b1;b2],[],[])
Something's wrong. We should have found a recession vector (bb<0).

Matt J

@Bhimavarapu,
The code doesn't generate any plots, but there are other utilities on the FEX that you could use for plotting polyhedra, e.g.,

As for "not getting the vertices", if you mean you are getting empty output V=[], it means your polyhedron appears empty to the code. Hard to say more without seeing what you're doing.

Bhimavarapu Reddy

Dear Sir,
Even I am giving A,b I am not getting the vertices and also not getting the plot

Matt J

Cong,

Points on the boundary of the polytope are expected to violate the inequalities by small amounts due to finite precision arithmetic. When I compute the violations for your A,b data with the code below, I find that they are indeed very small O(1e-13),

>> V=lcon2vert(A,b);
>> slacks=bsxfun(@minus, b,A*V');
>> violations = slacks(slacks<0)

violations =

1.0e-13 *

-0.1066
-0.0355
-0.5684
-0.0711
-0.0355
-0.0355
-0.5684
-0.0355
-0.2842
-0.5684
-0.1066
-0.2487
-0.3375

Cong

I just did a test about your program, and I found that an extreme vertex that obtained does not satisfy all inequalities. Note that I just test inequalities case (3D):
A=
-32.7680000000000 17.6400000000000 -3.20000000000000
-3.37500000000000 6.25000000000000 -1.50000000000000
17.5760000000000 2.56000000000000 2.60000000000000
79.5070000000000 10.8900000000000 4.30000000000000
32.7680000000000 -17.6400000000000 3.20000000000000
3.37500000000000 -6.25000000000000 1.50000000000000
-17.5760000000000 -2.56000000000000 -2.60000000000000
-79.5070000000000 -10.8900000000000 -4.30000000000000
b=
5
18
68
244
2
-11
-61
-219
One of the extreme vertices: [1.5530 4.9322 9.7233]
In order to test if this three values satisfy all inequalities, I multiply A by [1.5530 4.9322 9.7233]' to compare the new set of b values with the previous b values, and I found not all of new b values less than the original b values, which means they are not satisfied.

Ivan Bogun

Matt J

Hi Richard,
You might already have a modification that you like, but I encourage you to look at QLCON2VERT, in my release yesterday. It allows you to supply an x0 as you asked, but it also works when there are equality constraints, i.e. x0 only has to be in the relative interior. Also, QLCON2VERT will skip the boundedness check by default, which can save a lot of time if you know a priori that the set is bounded. Finally, LCON2VERT itself should be faster now for certain problems. It now checks first whether the origin is an interior point, something it can do very quickly, before performing more complicated searches.

Richard Katzwer

Disregard the previous question, I now see that it is easily modifiable. Simply
1. pass an extra argument x0 to con2vert
2. check if x0 in the polyhedron
2a. If yes, continue as normal with s = x0;
2b. Otherwise, set s = slackfun(pinv(A)*b) and proceed as if there were no x0.

Richard Katzwer

Great software, Matt J. I've found it to be far more robust than con2vert.

Quick question, though. Suppose I have a set of linear constraints [A,b], and I *know* a point x0 that satisfies A*x0 <= b. Is there a way to modify lcon2vert so that it doesn't have to work to find an interior point and can use a user-supplied one?

Matt J

Hi Andreas,

Thanks for your feedback.

Your example doesn't demonstrate a bug, as far as I can see. The rows of the matrix returned by lcon2vert are the vertices of the polyhedron. The vertices are not guaranteed to be listed in any particular order.

So, in your example, there is no reason that it would return eye(3), if that's what you were expecting -- only some permutation of the rows of eye(3).

Andreas

Hi Matt,

first of all thanks for the great work!

It seems that the example with V=eye(3) does not work. Dimensions are swapped when using the output of vert2lcon with lcon2vert. This is what I get:

>> V=eye(3)

V =

1 0 0
0 1 0
0 0 1

>> [A,b,Aeq,beq]=vert2lcon(V)

A =

0.4082 -0.8165 0.4082
0.4082 0.4082 -0.8165
-0.8165 0.4082 0.4082

b =

0.4082
0.4082
0.4082

Aeq =

0.5774 0.5774 0.5774

beq =

0.5774

>> lcon2vert(A,b,Aeq,beq)

ans =

1.0000 0.0000 0.0000
0.0000 0.0000 1.0000
-0.0000 1.0000 0.0000

I am using R2012a. Is there something that can be done about this "bug"?

Matt J

@Sandro: that would be a problem with the way you generate your constraints. It might be worth elaborating on how you do that in a post to the Newsgroup, or to Answers.

Sandro

Hi Matt,

The problem is that by changing input parameters I got different sets of inequality constraints, but I don't know beforehand whether some will degenerate to equality.

How do you think?

Matt J

Hi Sandro,

The inequality constraints x2<=0, -x2<=0 are equivalent to x2=0. You should use the Aeq input argument to express this.

Sandro

Hi Matt,

I just wonder whether lcon2vert can handle this case: 0<=x1<=1, x2<=0, -x2<=0, in R^2.

It seems to fail in this case. Any suggestion?

Best,
Sandro

Matt J

If you cannot upgrade, just replace the line with

[trash,E]=max(P);

Xu

My bad. Perhaps my Matlab version (2007) is too old.

Xu

Syntax error in line 89..... Could you please check that?

Matt J

Hi Cong,
No. There is no connection to Nelder-Mead that is apparent to me.

Cong

Hi Matt,
Thanks for your explanation. I was wondering that if the algorithm is related to the simplex method which proposed by Nelder and Mead?

Matt J

Hi Cong,
Thanks for your feedback. There is no formal reference that I can give you, I'm afraid. I largely reverse-engineered Michael Kleder's code. Approximately speaking, though, it works like this: every bounded polyhedron has a dual polyhedron. The face normals of the original polyhedron are the vertices of the dual polyedron and vice versa. Because we are given the face normals of the original polyhedron, we know the vertices of the dual. The code can use CONVHULLN to compute the dual face normals from the dual vertices. The dual face normals are the vertices of the original polyhedron.

Cong

Hi Matt,

I am sorry about my mistake, and you are right. I should only expand upper values or decrease lower values to test this program, and it give me the results quickly. In order to have a better understanding of this program, could you provide some references to the algorithm used?
Best regards,
Cong

Matt J

Hi Cong,

I disagree that your second data set defines a non-empty polyhedron. If it were non-empty, then I should be able to sum together any 2 of your inequalities and obtain another true inequality. But if I add your 7th and 17th inequality together, I obtain

0*x(1)+0*x(2) <= -9.5583e+005

which clearly cannot be satisfied by any point [x(1), x(2)].

Cong

Matt, I agree with you that my previous data does not exist a polytope. However, I think it failed to give a polytope when I use another data which definitely exists a polyhedron:
A=[0.1882 0.0936
0.2080 0.0853
0.2435 0.0393
0.2488 -0.0170
0.2482 -0.0211
0.2444 -0.0366
0.2287 -0.0668
0.2005 -0.0892
0.1939 -0.0918
0.1342 -0.0913
-0.1882 -0.0936
-0.2080 -0.0853
-0.2435 -0.0393
-0.2488 0.0170
-0.2482 0.0211
-0.2444 0.0366
-0.2287 0.0668
-0.2005 0.0892
-0.1939 0.0918
-0.1342 0.0913];

b=[-7495.00
-1048038.00
-11224.49
-11567.01
-10924.00
-11159.33
-968495.00
-9129.12
-7411.05
-6410.33
9718.79
1358994.00
14616.50
1512609.00
14285.24
14592.97
12664.93
11938.08
9650.65
8312.30];

If I do not let it run sufficient time which would be quite long, then the error is "Unable to locate a point near the interior of the feasible region". Thus, is it possible if I just let it run and it will provide a polyhedron sooner or later?

Best regards,
Cong

Matt J

Cong - it appears that the polyhedron is empty and therefore the "Initializers" cannot find a point inside it. I intend to modify the code soon so that it will bail out with V=[] when this occurs.

Cong

Hi Matt,
Thanks for providing the code.
However, the program will run a long time(unknown) if I use the data:
A=[0.188 0.094
0.208 0.085
0.243 0.039
0.249 -0.017
0.248 -0.021
0.244 -0.037
0.229 -0.067
0.201 -0.089
0.194 -0.092
0.134 -0.091
-0.188 -0.094
-0.208 -0.085
-0.243 -0.039
-0.249 0.017
-0.248 0.021
-0.244 0.037
-0.229 0.067
-0.201 0.089
-0.194 0.092
-0.134 0.091];
b=[-8495.00
-11480.40
-13224.50
-11567.00
-10924.00
-11159.30
-12684.95
-9129.12
-7411.05
-7410.33
9718.79
11589.90
17616.50
15126.10
14285.20
14593.00
16664.90
11938.10
9650.65
8312.30];
In addition, I find the "Initializer2" always spend much time on this calculation. Do you have any idea to make it more efficient?

Best regards,
Cong

Deepanshu Vasal

Matt J

@Deepanshu - actually the 2nd case doesn't appear to work for me, either. I think I've found the problem, but it will take me a few days until I can fix it.

Deepanshu Vasal

Matt, Thanks for your reply. Having equality constraint was not the problem. If I take
A = - [
0 0 0 0 0 1 -2 1 ;
0 0 1 -2 1 0 0 0 ;
eye(8);
-eye(8)];
and
b = - [zeros(2,1); -10*ones(16,1)];

it does not work but if I reduce redundant dimensions and put

A = - [
0 0 0 1 -2 1 ;
1 -2 1 0 0 0 ;
eye(6);
-eye(6)];
and
b = - [zeros(2,1); -10*ones(12,1)];

it works.

but I used lrs package at http://cgm.cs.mcgill.ca/~avis/C/lrs.html

and I got vertices for my problem.

Anyway I appreciate your effort and your such timely response. good luck

Deepanshu

Matt J

@Deepanshu - Are you sure the first 8 constraints are not meant to be equality constraints?

In any case, I answered you in an earlier email, but here's the recap: the region described by your inequality constraints appears to be non-solid. Such inequality constraints are numerically-unstable, which is why all MATLAB Optimization Toolbox functions ask you to represent your linear constraints in the form:

Aineq*x<=bineq
Aeq*x=beq

where Aineq*x=bineq describes a SOLID region.

LCON2VERT works similarly. You must rewrite your non-solid region as above and invoke lcon2vert in 4-argument form:

lcon2vert(Aineq,bineq,Aeq,beq)

Deepanshu Vasal

Hi Matt,

I am trying to use your code with my dataset, but I am getting singular matrix warning with the vertices generated are out of bounds (Nan, inf).
My data set is as follows :-

A = - [

0 1 0 -2 0 1 0 0 ;
0 0 0 0 0 1 -2 1 ;
0 0 1 -2 1 0 0 0 ;

-1 1 1 -1 0 0 0 0 ;
0 0 0 -1 1 1 -1 0 ;

0 0 0 1 -1 0 -1 1 ;
1 -1 0 -1 1 0 0 0 ;
0 0 1 -1 0 -1 1 0 ;
eye(8);
-eye(8)];

b = - [zeros(8,1); -10*ones(16,1)];

is there something that cannot be handled by this code? I see that the initial point taken is fine (satisfies the constraints.) but cannot figure out where is the problem in the code.
If you come across any problem, please let me know. Also do you know the reference for the algorithm used in this code?

Thanks so much.
Deepanshu
Thanks so much
Deepanshu

Matt J

@david: You need to tell me how I can reproduce the error.

david

Hi Matt,
I'll be very thankful if you helps me with this problem: When I applied your function i get the following error :
"??? Maximum recursion limit of 500 reached. Use set(0,'RecursionLimit',N)
to change the limit. Be aware that exceeding your available stack space can
crash MATLAB and/or your computer."

What i have to do in this case ?

Matt J

Hi Sandro,

I'm not having any such difficult in R2011a, so perhaps you need to upgrade. Here is the test that I ran, and no errors resulted

>> [a,b,ae,be]=vert2lcon(eye(20));
>> V=lcon2vert(a,b,ae,be);
>> norm(V'*V-eye(20))

ans =

5.7516e-015

Sandro

hi Matt,

Do you have any exposure of how high the dim could be? I am testing something like x1+x2+...+xn=1, xi>=0, i=1,2,...n. I randomly generate some points within this region and also give the vertices like (0,0,0..,0), (1,0,0,..0), (0,1,0,...0)...to help checking.

It looks like n=10 works, but n=20 is stuck due to qullmx function called by convexhulln.

Is there any possible we can have function like qullmx but more powerful in the sense of dim and speed?

Sandro

Harshavardhan Sundar

Hi Matt, Thanks again (Officially)....the latest of your code lcon2vert and vert2lcon are working absolutely fine....

Matt J

@Harshavardhan - I've just updated the package with what I believe is a more robust version of LCON2VERT. At the very least, it successfully handles your example data set. Also, it will more reliably report cases where it cannot find a solution.

Matt J

@Harshavardhan - This is a resend of a previous comment which seems to have gotten lost by the FEX.

The data set you've provided is quite intriguing! I've traced the problem you're seeing to Michael Kleder's CON2VERT code (my lcon2vert is mainly a wrapper for this).

The problem is that, in order to work, CON2VERT needs to find an initial point c inside the polytope. The simple method that it uses seems to have had a good rate of success historical success, but lo and behold, your particular A,b data set has defeated it!

I know this is the problem however, because

(i) I've checked that the initial point c that it generates does not satisfy A*c<b

(ii) When I manually change CON2VERT, hardcoding the initial c to a point that I know is in your polytope, I get the correct vertices.

I've alerted Michael Kleder to this via the CON2VERT webpage, but he doesn't appear to have been monitoring the page for several years. I doubt he'll be responding with any input.

I have some ideas of my own to improve the initialization scheme, but they will take some time to develop.

In the meantime, if you're in a situation where you already know an initial point, c, inside your polytope, a workaround for you might be to modify CON2VERT so that it accepts this c is an additional input argument.

Beyond that, I don't know what to say. You have difficult data...

Matt J

@Harshavardhan - one additional note. In the meantime, you should also modify this line of code

if ef ~= 1

to this

if ~all(A*c < b)

This will at least warn you properly that the function is about to fail.

Hi Matt. First of all thanks a lot for the wonderful piece of code. The code works fine for most of the cases. However there are a few occasions where I get erroneous vertices. One such example is below. My A matrix is :
A=[1 0 0
0 1 0
0 0 1
-1 0 0
0 -1 0
0 0 -1
-0.387696185212551 -0.295276274520099 0.703237014010810
5.73326513862449 -0.422448628998633 -0.113292931747631
-3.79255233000727 -0.360645958414035 -0.673892379691217
-1.27188203250127 0.544750863219514 -0.791324376598499
-2.30540199426269 -2.50522755054822 -0.291771858420232
0.616578245027929 -0.654271120193835 0.574406987898971
-1.49508806389885 -2.68877782397954 -0.989225453798092
5.34556895341194 -0.717724903518731 0.589944082263178
-3.40485614479472 -0.0653696838939366 -1.37712939370203
-0.884185847288714 0.840027137739613 -1.49456139060931
-1.91770580905014 -2.20995127602812 -0.995008872431042
0.228882059815378 -0.949547394713933 1.27764400190978
-1.10739187868630 -2.39350154945944 -1.69246246780890
1.94071280861722 -0.783094587412668 -0.787185311438849
4.46138310612323 0.122302234220881 -0.904617308346131
3.42786314436180 -2.92767617954685 -0.405064790167863
5.11668689359657 0.231822491195202 -0.687699919646602
4.23817707472564 -3.11122645297817 -1.10251838554572
-2.52067029750601 -0.905396821633549 0.117431996907282
1.48715033574458 -2.14458159213419 0.382120521270985
-3.17597408497934 -1.01491707860787 -0.0994853917922465
2.29746426610842 -2.32813186556550 -0.315333074106875
-1.03351996176143 -3.04997841376774 0.499552518178267
-0.655303787473336 -0.109520256974320 -0.216917388699529
-0.223206031397585 -3.23352868719905 -0.197901077199593
-1.68882374923476 -3.15949867074206 0.282635129478739
-0.810313930363841 0.183550273431316 0.697453595377860
-0.878509818870922 -3.34304894417337 -0.414818465899121]
b=[6
4
2
0
0
0
-0.154147860211922
16.6671436561829
-8.83179250583593
-1.48693603688105
-8.92116339978305
0.423928365940654
-8.95704617639833
16.5129957959710
-8.67764464562401
-1.33278817666913
-8.76701553957113
0.269780505728732
-8.80289831618640
7.83535115034699
15.1802076193019
7.74598025639987
16.2432152902423
7.71009747978459
-7.34485646895489
-0.0893708939471196
-8.40786413989528
-0.125253670562397
-7.43422736290201
-1.06300767094039
-7.47011013951728
-8.49723503384240
0.0358827766152754
-8.53311781045768]

I used the plotregion command available in MATLAB central to plot the region. I have reasons to believe that the plotregion command is plotting the region accurately. However for this particular example the vertices I get (graphically) from the plotregion command are not coinciding with the vertices obtained from lcon2vert. I also tried con2vert by Michael Kleder with little success.

Timothy Stratton

Nice very powerful!

Sandro

Hi Matt,

my concern of using consolidator is because I think function unique has pretty high requirement of two vectors to be the same. When run vert2lcon on a large data set, i have seen a large A with many rows are the same up to 6-8 digital number. This enlarges the size of A and it may not be necessary for the user.

"A and b are produced by Michael Kleder's functions which seem to normalize it this way already."

I doubt about it.

Best,
Sandro

Matt J

Hi Sandro. Thanks for the feedback. I'll think about your recommendations. Could you explain the advantage of consolidator over unique? I chose a normalization for Aeq and beq that led to nice whole integers where possible. A and b are produced by Michael Kleder's functions which seem to normalize it this way already.

Regarding your rating. It's fine. The FEX only scores according to your most recent rating. This is to allow people to update their ratings if the author makes improvements to the submission.

Sandro

Sorry Matt, I am not sure how to remove the first rating. Two comments: "unique" function may be better to be replaced by "consolidator"; It may be good to normalize A,b,etc before the final output(I am not quite sure about its necessity.)

Matt J

Ah, ok.

Sandro

It is definitely a good job!!! I may mess up the rate submission above.

Sandro

MATLAB Release Compatibility
Created with R2010b
Compatible with any release
Platform Compatibility
Windows macOS Linux

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!