Potential issue with LRS

[worc4021]

Dear Benoit,

first off: This is not really an issue this repository directly! I have also written an interface to the lrs library but for Matlab. (https://github.com/worc4021/GeoCalcLib.git) It calls the lrs_set_row_mp function to construct the data internally. Everything works fine until you try calculate a facet enumeration of a set of points that do not have the origin in their convex hull. If you also communicate with the lrs internal methods by calling their set_row functions then you might experience similar issues. I couldn't really work out how you construct the internal dic and dat types. If you did work out how to resolve this issue in any way please let me know, as I spent ages reverse engineering the text file based lrs executable which led me nowhere.

All the best, Manuel

[blegat]

Hey, I am glad you share your experience. I first wrote a wrapper for CDD trying not to modify CDD's source code but I gave up ^^ CDD, LRS, ... are, as you say, text file based and trying to use them without text files is quite challenging. So forked cddlib and I added functions to CDD to set the matrices and other fields and it works great.

LRS is even more text file based than CDD. It does not even store the whole output. It outputs it incrementally (which wouldn't be possible for CDD since it uses a different algorithm). This design of LRS has the advantage that the set of vertices does not have to fit in memory. I plan to support something like that in Polyhedra in the future.

With that in mind, I first wanted to use LRS in a non-text file based interface. I didn't try to use lrs_set_row_mp and forked it immediately. In lrs_main, it does a few things.

1. First it calls lrs_alloc_dat.
2. Then it calls lrs_read_dat to read the dimension of the matrix.
3. Then it calls lrs_alloc_dic using this information to allocate the matrix.
4. Then it calls lrs_read_dic to populate the matrix.
5. Then it computes the vertices and rays and output them as soon as they are found.

I created two new functions: lrs_init_dat and lrs_init_dic_vec to be able to initialize the structures from Julia without using files. I tried to initialize them and then extract the matrix from Julia and it works, the code is in this package.

What I still need to do is to add a function for the 5th step in my fork that store the result in a matrix rather than storing it in a file.

What do you think ?

[worc4021]

Hi Benoit,

Thank you for your response.

You might find it helpful to look at what I've done in my implementation. It also works without any intermediate text files. The main procedure is

1. Create the necessary data (in particular A,b for a vertex enumeration or V and R for a facet enumeration) from the data we supply from outside in the GMP format mpq_t which is the rational format in order to not lose precision and also because what the LRS does internally is really difficult to keep track of. (This is essentially one huge row major mpq_t array)
2. Using GMP inbuilt functions i extract the numerator and denominator of each row in the construct.
3. I use the the lrs_set_row_mp to feed it into the lrs_dic and lrs_dat structure.
4. I mimic what is outlined in vedemo.c and chdemo.c to compute the the vertex enumeration or facet enumeration and extract vertex by vertex (facet by facet) into my GMPmat format.

This all works perfectly fine as long as I don't try to go from a set of points that does not contain the origin in its convex hull to a facet description. I have contacted David Avi in order to figure out whether this is a problem with the LRS or I am just missing a flag. LLDB is only so helpful when trying to work out what flags are set in the lrs call and which ones in my interface.

Using something like the GMPmat format I use has the advantage that the data is comprehensible (to some extent) since comparing the dictionaries of two different runs can be extremely intransparent. I agree that the dictionary description might be a lot more efficient in terms of memory.

When I created the testcase that reproduces my issue without anywhere relying on Matlab I found that the functions are pretty easily plug&play. It should be almost trivial to get the data from julia into the GMP format and then back into the Float64 when it's done.

[blegat]

The Julia BigInt type uses GMP so I do not need to translate it to GMP :-P I am not sure to have understood your last paragraph: Do you have an example that does not work even if you do not use your MATLAB wrapper and directly use LRS ? The issue when the origin might be that a convex set need to contain the origin for its duality/polarity to work. For example "the polar of the polar of a closed convex set containing the origin is itself". (see here for example). This is why if you have a convex set that is not a cone, it is lifted (we add one dimension that is should be equal to one).

• Ax >= b becomes [-b A][z; x] >= 0 with z = 1
• A ray r becomes [0; r] and a vertex v becomes [1;v].

Once you have a cone, there is no issue with duality/polarity and the fact that the origin is inside the initial convex set shouldn't matter. If you have not lifted your convex set, LRS thinks that you have given a cone and hence behave unexpectingly. That would explain the behaviour that you explain.

[blegat]

I have just looked at chdemo.c and vedemo.c, that is very helpful. I thought that was meant to be used only from lrs_main but seeing this demo, its seems that as you said, the cleaner way is to use lrs_set_row. I think I will delete my function lrs_init_dic_vec and I will not add a function for the 5th step but instead get each vertices/rays on Julia. That will be much cleaner and less intrusive.

[worc4021]

I agree. The point where I get lost is that the lrs_main algorithm performs its homogenisation that does the lifting itself, but when you ask it to do it by setting Q->homogeneous it does not.

The way I worked around this is by shifting the entire set of points into its barycentre and then shifting the facet representation back again.

It is just important to notice that the LRS library seems to ignore the setting of the Q->homogeneous flag if it is set manually.

[blegat]

It seemed to me that Q->homogeneous meant the opposite.

• If Q->homogeneous is set, that means that there is no vertex, it is a cone so no lift needed.
• If Q->polytope is set, that means that there is no ray, it is a polytope, lift needed.
• If none is set, there are vertices and rays, lift needed.

I don't know what is the optimization that it does for Q->polytope though.

[worc4021]

This is what you would think. As far as I can tell the choice of either of them does not affect anything.

[blegat]

I have just read very carefully lrs_set_row_mp and here is what I understand: Q->hull is very important to be set but Q->homogeneous and Q->polytope are computed automatically so they should not be set (then why is it set in chdemo.c ??????).

Let D be the unlifted dimension.

When you create the dat, you give Q->n = D + 1. If hull = 0, d = Q->n-1 = D and if hull = 1, d = Q->n = D+1 (weird that d != D). Then the matrix A has dimension Q->m x d+1.

If hull = 0, it is the H-representation b + Ax >= 0. A[i][0] has b_i and A[i][1:d] has A_i. num and den are supposed to have length d+1=D+1 and contain b_i followed by A_i. If all b_i are 0 then it detects automatically that Q->homogeneous should be true. Q->polytope is set to false since it cannot be detected from the H-representation.

If hull = 1, it is the V-representation with vertices and rays. A[i][0] contains 0 and A[i][1:d] is num. num and den are supposed to have length d=D+1 and contain 1 followed by a vertex or 0 followed by a ray. If A[i][1] is always 0 then it detects automatically that Q->homogeneous should be true. If A[i][0] is always 1 then it detects automatically that Q->polytope is should be true.

[blegat]

Can you give me your example that fails ?

[worc4021]

The way I did this is that I followed the description of 'Building the data matrix' given in http://cgm.cs.mcgill.ca/%7Eavis/C/lrslib/lrslib.html, this to me explains that the matrix I feed in is basically [type,data] where type is 1 for a vertex and 0 for a ray and data is the vertices and rays stacked up.

I have sent this to David Avi and he sent me this:

I suspect that when you are doing a facet enumeration you do not have Q->hull=1 set.

lrs by default does vertex/ray enumeration. For facet enumeration lrs lifts the vertices to a cone in one higher dimension by adding a column of zeroes in column 1 (homogenization) then essentially both types of enumeration are the same and it just does vertex/ray enumeration. If V-representation is in the input file it sets Q->hull=1 to enable this.

In your file, if you delete the 'V-representation' line and run lrs you should get:

Totals: vertices=16 rays=2 bases=16 integer_vertices=0 vertices+rays=18 Note! Duplicate rays may be present Dictionary Cache: max size= 10 misses= 4/15 Tree Depth= 13 lrs:lrslib v.6.0 2015.10.10(32bit,lrsgmp.h) *0.015u 0.000s 4020Kb 1025 flts 0 swaps 0 blks-in 0 baks-out

which is correct when the input is interpreted as inequalities. Now if you insert a column of zeroes and set n=4 you get:

Totals: vertices=1 rays=30 bases=28 integer_vertices=1 vertices+rays=31 Dictionary Cache: max size= 10 misses= 18/27 Tree Depth= 27 lrs:lrslib v.6.0 2015.10.10(32bit,lrsgmp.h) 0.015u 0.000s 4028Kb 1027 flts 0 swaps 0 blks-in 0 baks-out

Which now corresponds to what you get from lrs from the exp.ext file, except that 'rays' are 'facets':

Volume= 3953390694893980748961675022025/1267650600228229401496703205376 Totals: facets=30 bases=28 Dictionary Cache: max size= 10 misses= 18/27 Tree Depth= 27 lrs:lrslib v.6.0 2015.10.10(32bit,lrsgmp.h) *0.000u 0.015s 4040Kb 1030 flts 0 swaps 0 blks-in 0 baks-out

I actually do set the Q->hull flag, so in theory the LRS algorithm should do the lift itself. From what I see it might be easier and quicker to resolve this by explicitly adding the additional column and using the vertex enumeration instead to obtain a facet description.

The stimulus that does produce this malfunction is in my testcase.c.. Basically its a bunch of vertices that are [sin(t)+5,cos(t)+5], i.e. they lie on the unit circle so facet description has exactly one less face than there are vertices and they also look like [sin,cos] and the right hand side should be shifted.

[worc4021]

So for me this works. I hope it does for you too!

I've actually desperately been waiting for the Polyhedra package to work in julia. So best of luck!

[blegat]

I have implemented the change of representation and it seems to work but I haven't tested it yet on your example.

You added a column of zero and didn't use Q->hull and it worked ?

[worc4021]

Yes exactly! What then comes out is one 'vertex' at the origin, [1 0 0 0 ...] and the rest are all rays that have the form [0,b,-A]. But I avoided using the internal homogenisation which didn't work for me.

[worc4021]

I have news:

I internally I used the GMP rational format mpq_t whereas the LRS library uses the integer format mpz_t, the output of the facet enumeration is given as a + b_x + c_y >= 0. Since it comes out as purely integer a, b and c are horribly large, which makes direct conversion to double (remember this was for an interface to Matlab) unwise. So I converted them into 1 + b/a_x + c/a_y >= 0, which is exactly the same if a is positive, but since a is negative if the convex hull of the set of points does not contain the origin it flipped the orientation of the inequality. All inequality reduction algorithms applied after would produce only nonsense, but that is not really a malfunction. Rather me being an idiot and forgetting that case very early when I adopted the rational format over the integer one.

So the good news is that you hopefully wouldn't make such a silly mistake and will produce sensible results with your wrapper.

[blegat]

That nasty bug, I'm glad you figured it out. Fortunately, I do not need any conversion so I keep the inequalities as is. I only need to divide when it is a vertex. I need to transform

den x_1 x_2 ... x_d

to

x_1/den ... x_d/den

but normally den is positive.

[worc4021]

This is exactly how I introduced this error: The output of the LRS for a vertex enumeration is [lcd v1 v2...] as you say and the lowest common denominator is always positive. The vertex enumeration was the first thing I implemented back then, and when I added the facet enumeration I was already extracting the data that way and ignored the case.

Are you using a lifted version then? That is how I solved it, although now it seems redundant, but it's already working.

[blegat]

No I don't add the zero column manually but I haven't tested your example yet