worc4021 / GeoCalcLib

Interface for MATLAB to use LRS library and derivative functionality.
GNU General Public License v2.0
7 stars 6 forks source link

Large facet enumeration problem #3

Open cja222 opened 3 years ago

cja222 commented 3 years ago

Hello,

I am trying to solve large facet enumeration problems with GeoCalcLib. It works great for fairly large problems -- I have a V-representation cone in 18 dimensions with 30 rays. The resulting H-representation has about 50,000 facets. As a check, I can solve a linear program over the H and V representations separately and they give the same answer.

However, I would like to add more rays to the V-representation and perform facet enumeration on, say, 40 rays in 18 dimensions. When I do this, GeoCalcLib/lrs produces an H-representation with about 4,000,000 facets. However, in my test of solving a linear program over the H and V representations, each gives a different result. In general, do you know if LRS struggles with very large outputs? Or could there be numerical precision issues at this scale - might you recommend changing any settings for large problems?

Thanks!

worc4021 commented 3 years ago

Hi,

To be honest I never used the LRS algorithm for such huge problems. I used it for academically small 2-5 dimensional problems. So I cannot give you a satisfactory response from first hand experience. However, the version we use here does its computation using the GMP library, that means that it is meant to be lossless for any rational representable polyhedron, it will just get excessively slow reallocating more and more memory to hold the representation. We only go back to floating point precision, where we could be incurring precision loss when we go from the mex into Matlab. So if I had to predict what would happen if we solve a linear program on the halfspace representation of a particular polyhedron, I would expect that it pretty much matches the solution of a linear program on the vertex representation. Do you see the minimum itself vary much or only the solution?

If you want more details about the scaling ability of the algorithm I'm afraid you will have to ask David Avis directly.

In response to previous issues I created a more modern interface to the lrs library but never got around to documenting it, potentially more recent versions of the LRS library perform better, which might be useful for your 18 dimensional analysis. Whether it changes the numerical accuracy I couldn't comment on.

Good luck! Manuel

cja222 commented 3 years ago

Hi Manuel,

Thanks for getting back to me. For the "smaller" problem with 30 rays, the linear program results exactly match for V and H-representations, but for the larger problem with 40 rays they are completely different solutions for the V and H representations. Actually I tried sequentially increasing the number of rays from 30 to 40 and at a specific point the linear program results stop matching. Do you know what happens if memory is maxed out during the LRS run - does GMP abort, or could it continue with a loss in precision or some other issue? I am running on a laptop with 16 GB RAM, which is entirely used throughout the run.

I'll try getting the most recent LRS library set up, maybe that would make a difference.

worc4021 commented 3 years ago

Hello,

sorry for the long delay. For some reason my emails with the notifications about this aren't coming through..

I would expect that if you eventually run out of memory that your mex file just crashes and takes matlab out with it because it fails some assertion. You could easily notice that if you run top while the enumeration is running. You should notice it consuming more and more memory, but before that happens it would go and spill into your hard drive, i.e. I would not expect that to be your issue.

It's been a few years now since I worked with polyhedra in anger during my phd, but I remember that there exist polyhedra that have a rational H representation but no rational V representation and vice versa. So an H representation that you can potentially hold on a finite precision machine but its V representation does not have a rational and therefore computer cannot capture it. I had been under the impression that these were reasonably rare and don't exist in low dimensions. But since you're working on much higher dimensions I wonder whether you are running into those. In that case I have no idea what the LRS would find when it enumerates, as it will certainly find something. But the graph of what it finds is not the dual to the graph you had put in.

Where does your problem arise? I've always started from an H representation that came from some for of inequality constraints and had to compute some guarantees with those. You're starting from the rays and vertices themselves..?

cja222 commented 3 years ago

No worries! I'm actually just finishing up my PhD now so I'm a bit preoccupied anyways. I'm working with a V representation cone - the origin and a set of rays. My goal is essentially to compute the "central direction" inside the cone that is furthest from the boundaries of the cone, which can be computed by solving a linear program given the H-representation (with an additional constraint added so that the LP is bounded). Performing facet enumeration then solving the LP works just fine even for large problems (up to say 35 rays in 18 dimension). But there is a point where if I include more rays, say 40, LRS still gives and output and I can compute a "central direction", but something is goes wrong because this direction does not lie within the original V representation of the cone.

Of course there could be a bug somewhere in my code, but it's strange that it works for most examples I've tried. So what you're describing sounds plausible, if the H representation that is output does not correspond the V representation. Do you have any references that discuss cases like this or any reasoning for why this might occur?

worc4021 commented 3 years ago

That's interesting! I would have expected something more complex to result to find the central ray, but I guess that is influenced by your chosen norm.. Very cool!

I don't know who would discuss computational issues relating to these polyhedra, but the basics you'll find in Ziegler - Lectures on Polytopes. Bear in mind that the book doesn't care about the computational aspects themselves. But they result as a natural consequence. The first hint at the issue I'm thinking about is around Prop 2.17 about rational polytopes, the actual examples for this are later e.g. Example 6.21. But again, no mention of what you'd find numerically.

You'd probably end up having to define something like a metric on the combinatorial makeup of polyhedra and try to find out what a perturbation in that space would yield in the realisation of the polytopes. I.e. if you have a nonrational polytope and you try to find a rational approximate, does it resemble the nonrational one? If so, in what sense?! All the enumeration methods rely on the graphs of polyhedra; actually simplex based linear solvers do too!! Thinking about it this way, it's probably way less common to run into those by chance. Otherwise there would be way more literature about simplex searches in nonrational polytopes. If you cannot represent a polytope in a computer then searching for a vertex in it makes no sense.. :/