Some optimizations for method regions of hyperplane arrangements

[gh-kliem]

There are some easy optimzations that speed up the computation of regions of hyperplane arrangements:

• To see if an hyperplane intersects a polyhedron non-trivially can be done by quickly checking the evaluation of the hyperplane normal on the Vrepresentation. This is much faster than computing two new polyhedra every time.
• In case of a linear hyperplane arrangement, we can only consider the upper half space with respect to the first hyperplane and then append to each region the negative.

Before this ticket:

sage: from itertools import product
sage: def zero_one(d):
....:     for x in product([0,1], repeat=d):
....:         if any(y for y in x):
....:             yield [0] + list(x)
....:          
sage: K.<x,y,z> = HyperplaneArrangements(QQ)
sage: A = K(*zero_one(3))
sage: %time len(A.regions())
CPU times: user 141 ms, sys: 7.6 ms, total: 148 ms
Wall time: 146 ms
32
sage: K.<x,y,z,w> = HyperplaneArrangements(QQ)
sage: A = K(*zero_one(4))
sage: %time len(A.regions())
CPU times: user 1.64 s, sys: 18.2 ms, total: 1.66 s
Wall time: 1.66 s
370
sage: K.<x,y,z,w,r> = HyperplaneArrangements(QQ)
sage: A = K(*zero_one(5))
sage: %time len(A.regions())
CPU times: user 2min 2s, sys: 298 ms, total: 2min 2s
Wall time: 2min 2s
11292

With this ticket:

sage: K.<x,y,z,w,r> = HyperplaneArrangements(QQ)
sage: A = K(*zero_one(5))
sage: %time len(A.regions())
CPU times: user 22.7 s, sys: 73.1 ms, total: 22.8 s
Wall time: 22.8 s
11292
sage: K.<x,y,z> = HyperplaneArrangements(QQ)
sage: A = K(*zero_one(3))
sage: %time len(A.regions())
CPU times: user 46.4 ms, sys: 8.03 ms, total: 54.4 ms
Wall time: 53 ms
32
sage: K.<x,y,z,w> = HyperplaneArrangements(QQ)
sage: A = K(*zero_one(4))
sage: %time len(A.regions())
CPU times: user 1.05 s, sys: 4.01 ms, total: 1.05 s
Wall time: 1.05 s
370
sage: K.<x,y,z,w,r> = HyperplaneArrangements(QQ)
sage: A = K(*zero_one(5))
sage: %time len(A.regions())
CPU times: user 22.2 s, sys: 31.8 ms, total: 22.2 s
Wall time: 22.2 s
11292
sage: %time y = tuple(-P for P in A.regions())
CPU times: user 10.7 s, sys: 124 ms, total: 10.8 s
Wall time: 10.8 s

The last timing indicates that half of the time in this example is taken just to create the object and cannot be improved (in this method).

[gh-kliem]

Note that this ticket is motivated by a new implementation in polymake by Lars Kastner and Marta Panizzut, see ​https://arxiv.org/pdf/2003.13548.pdf.

The example to benchmark is taken from their paper.

---

New commits:

​ad162d7

optimize region computation of hyperplane arrangements

[jipilab]

One comment:

Probably a comment like

+                # Determine if all vertices lie on one side of the hyperplane.
+                splits = False
+                direction = 0

Would complete the following one:

+                if not splits:
+                    # All vertices lie in one closed halfspace of the hyperplane.

I am a bit perturbed by the for loop followed by an if statement that ultimately modify the value of splits. Is there a way to fusion them and keep the efficiency provided by the break?

But now that I look at it more, probably it is the best way to do so. Hm.

Another thing: if the algorithm you wrote is described in the article you mention, it should be added in the documentation of the function.

[gh-kliem]

No, it's not their algorithm. It's exactly the algorithm that Volker Braun implemented back in 2013 (according to git blame). It's just that I looked at the paper and their timings and compared it to Sage. Sage was faster, but I realized that with the above two simple observations we can make it even faster. I would expect their algorithm to perform better asymptotically and/or on some other classes of hyperplane arrangements. Their algorithm computes first a full-dimensional cell and then works it's way through the region by "reflecting" on the hyperplanes. So I would guess that they pay a big penalty for that, which only pays of in large examples.

I came up with both ideas optimizing Sage's approach myself. I don't see anything like this mentioned in their paper, but I will ask them about that.

[jipilab]

Replying to gh-kliem:

No, it's not their algorithm. It's exactly the algorithm that Volker Braun implemented back in 2013 (according to git blame). It's just that I looked at the paper and their timings and compared it to Sage. Sage was faster, but I realized that with the above two simple observations we can make it even faster. I would expect their algorithm to perform better asymptotically and/or on some other classes of hyperplane arrangements. Their algorithm computes first a full-dimensional cell and then works it's way through the region by "reflecting" on the hyperplanes. So I would guess that they pay a big penalty for that, which only pays of in large examples.

I came up with both ideas optimizing Sage's approach myself. I don't see anything like this mentioned in their paper, but I will ask them about that.

Ok, fine. Just making sure.

[gh-kliem]

Btw, is there a way to obtain a n-dimensional hyperplane arrangement without concretely specifying the generators?

---

New commits:

​79e5f80	optimize region computation of hyperplane arrangements
​61da5b6	readability of code
​b05df00	added example

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​9a690e4

sorting entries by appearance year

[jipilab]

Replying to gh-kliem:

Btw, is there a way to obtain a n-dimensional hyperplane arrangement without concretely specifying the generators?

You mean for example something like a library of hyperplane arrangements?

There is

sage: p = polytopes.simplex()
sage: p.hyperplane_arrangement()
Arrangement of 5 hyperplanes of dimension 4 and rank 4

Do you mean whether hyperplane arrangements accept other types of input?

[gh-kliem]

Thanks. That's not exactly what I meant. But it helped.

I was looking for something like

HyperplaneArrangements(ZZ, dim=5)

that just figures out some meaningful generators, but I guess one extra line does the job (as is illustrated in the method hyperplane_arrangement):

names = tuple('t' + str(i) for i in range(5))
HyperplaneArrangements(ZZ, names)

[jipilab]

Replying to gh-kliem:

Thanks. That's not exactly what I meant. But it helped.

I was looking for something like

HyperplaneArrangements(ZZ, dim=5)

that just figures out some meaningful generators, but I guess one extra line does the job (as is illustrated in the method hyperplane_arrangement):

names = tuple('t' + str(i) for i in range(5))
HyperplaneArrangements(ZZ, names)

Ahhh. Yes, I got you. Yes... I like this method too, it is better for internal usage...

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​00103d0

fix empty hyperplane arrangement

[tscrim]

Two minor comments:

-Example 6 of [KP2020]::
+Example 6 of [KP2020]_::

in order to get the link.

To avoid some repeated calls:

+                    region_lines = region.lines()
                     if direction == 0:
                         # In this case all vertices lie on the hyperplane and we must
                         # check if rays are contained in one closed halfspace given by the hyperplane.
                         valuations = tuple(ieq[1:]*ray[:] for ray in region.rays())
-                        if region.lines():
-                            valuations += tuple(ieq[1:]*line[:] for line in region.lines())
-                            valuations += tuple(-ieq[1:]*line[:] for line in region.lines())
+                        if region_lines:
+                            valuations += tuple(ieq[1:]*line[:] for line in region_lines)
+                            valuations += tuple(-ieq[1:]*line[:] for line in region_lines)
                         if any(x > 0 for x in valuations) and any(x < 0 for x in valuations):
                             splits = True
                     else:
                         # In this case, at least one of the vertices is not on the hyperplane.
                         # So we check if any ray or line pokes the hyperplane.
-                        if any(ieq[1:]*r[:]*direction < 0 for r in region.rays()) or \
-                                any(ieq[1:]*l[:] != 0 for l in region.lines()):
+                        if any(ieq[1:]*r[:]*direction < 0 for r in region.rays()) or
+                               any(ieq[1:]*l[:] != 0 for l in region_lines):
                             splits = True

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​f18fd2a

small improvments

[gh-kliem]

Thank you. I applied that.

The backslash in the if clause is needed though. I also don't understand why I should adjust the alignment the way you suggested. It looks nicer, but it suggests that the second clause is part of the any-clause, doesn't it?

Replying to tscrim:

Two minor comments:

-Example 6 of [KP2020]::
+Example 6 of [KP2020]_::

in order to get the link.

To avoid some repeated calls:

+                    region_lines = region.lines()
                     if direction == 0:
                         # In this case all vertices lie on the hyperplane and we must
                         # check if rays are contained in one closed halfspace given by the hyperplane.
                         valuations = tuple(ieq[1:]*ray[:] for ray in region.rays())
-                        if region.lines():
-                            valuations += tuple(ieq[1:]*line[:] for line in region.lines())
-                            valuations += tuple(-ieq[1:]*line[:] for line in region.lines())
+                        if region_lines:
+                            valuations += tuple(ieq[1:]*line[:] for line in region_lines)
+                            valuations += tuple(-ieq[1:]*line[:] for line in region_lines)
                         if any(x > 0 for x in valuations) and any(x < 0 for x in valuations):
                             splits = True
                     else:
                         # In this case, at least one of the vertices is not on the hyperplane.
                         # So we check if any ray or line pokes the hyperplane.
-                        if any(ieq[1:]*r[:]*direction < 0 for r in region.rays()) or \
-                                any(ieq[1:]*l[:] != 0 for l in region.lines()):
+                        if any(ieq[1:]*r[:]*direction < 0 for r in region.rays()) or
+                               any(ieq[1:]*l[:] != 0 for l in region_lines):
                             splits = True

[tscrim]

Replying to gh-kliem:

The backslash in the if clause is needed though.

Ah right. I miscounted the number of parentheses.

I also don't understand why I should adjust the alignment the way you suggested. It looks nicer, but it suggests that the second clause is part of the any-clause, doesn't it?

See the previous point. I agree with not changing the alignment. If I was coding this, I would have done it like this:

                        if (any(ieq[1:]*r[:]*direction < 0 for r in region.rays())
                            or any(ieq[1:]*l[:] != 0 for l in region_lines)):

However, that is my style and my general desire to avoid the \ so it generates an error if things are not properly matched. So don't feel you have to change it.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​59dc074

backslash and alignment

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​a2da775

alignment

[gh-kliem]

I never thought about using parenthesis to avoid backslash. I like that.

The reason that I often do something like or or + at the end of the line, is that it makes it easier to align both lines alike:

if (any(ieq[1:]*r[:]*direction < 0 for r in region.rays()) or
    any(ieq[1:]*l[:]          != 0 for l in region_lines)):

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​930217f

alignment again

[jipilab]

It looks good to me. If Travis agrees, I suggest to give positive review.

[tscrim]

Yep, I am happy with it. Thanks.

[gh-kliem]

Thank you.