Opened 2 years ago

Closed 2 years ago

#29661 closed enhancement (fixed)

Some optimizations for method regions of hyperplane arrangements

Reported by: gh-kliem Owned by:
Priority: major Milestone: sage-9.2
Component: geometry Keywords: hyperplane arrangements, regions
Cc: Jean-Philippe Labbé, Laith Rastanawi, Sophia Elia Merged in:
Authors: Jonathan Kliem Reviewers: Jean-Philippe Labbé, Travis Scrimshaw
Report Upstream: N/A Work issues:
Branch: 930217f (Commits, GitHub, GitLab) Commit: 930217fb8b5979ee798da8fd8542048bc1e1a90f
Dependencies: Stopgaps:

Status badges

Description

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

Change History (23)

comment:1 Changed 2 years ago by gh-kliem

Branch: public/29661
Commit: ad162d70516726db61edec848c82cdef99863064
Status: newneeds_review

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:

ad162d7optimize region computation of hyperplane arrangements

comment:2 Changed 2 years ago by Jean-Philippe Labbé

Status: needs_reviewneeds_work

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:

!#diff
+                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.

Version 0, edited 2 years ago by Jean-Philippe Labbé (next)

comment:3 Changed 2 years ago by 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.

comment:4 in reply to:  3 Changed 2 years ago by Jean-Philippe Labbé

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.

comment:5 Changed 2 years ago by gh-kliem

Branch: public/29661public/29661-reb
Commit: ad162d70516726db61edec848c82cdef99863064b05df006ca52e7943b73df4b63124aad5b9654f8
Status: needs_workneeds_review

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


New commits:

79e5f80optimize region computation of hyperplane arrangements
61da5b6readability of code
b05df00added example

comment:6 Changed 2 years ago by git

Commit: b05df006ca52e7943b73df4b63124aad5b9654f89a690e4b96331cf9e66128e8326d088d81fdf347

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

9a690e4sorting entries by appearance year

comment:7 in reply to:  5 Changed 2 years ago by Jean-Philippe Labbé

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?

comment:8 Changed 2 years ago by 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)

comment:9 in reply to:  8 Changed 2 years ago by Jean-Philippe Labbé

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

comment:10 Changed 2 years ago by git

Commit: 9a690e4b96331cf9e66128e8326d088d81fdf34700103d0a6a3aeb0c1512a241316b8904d825de8e

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

00103d0fix empty hyperplane arrangement

comment:11 Changed 2 years ago by Travis Scrimshaw

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
Last edited 2 years ago by Travis Scrimshaw (previous) (diff)

comment:12 Changed 2 years ago by git

Commit: 00103d0a6a3aeb0c1512a241316b8904d825de8ef18fd2ab164f6abf95fd4e75219e17c85a5b75d6

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

f18fd2asmall improvments

comment:13 in reply to:  11 ; Changed 2 years ago by 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

comment:14 in reply to:  13 Changed 2 years ago by Travis Scrimshaw

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.

comment:15 Changed 2 years ago by git

Commit: f18fd2ab164f6abf95fd4e75219e17c85a5b75d659dc0740b561563f81e8bfb0731ea5e90245da92

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

59dc074backslash and alignment

comment:16 Changed 2 years ago by git

Commit: 59dc0740b561563f81e8bfb0731ea5e90245da92a2da77530b455059013f3ab6201df939d662f921

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

a2da775alignment

comment:17 Changed 2 years ago by 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)):

comment:18 Changed 2 years ago by git

Commit: a2da77530b455059013f3ab6201df939d662f921930217fb8b5979ee798da8fd8542048bc1e1a90f

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

930217falignment again

comment:19 Changed 2 years ago by Jean-Philippe Labbé

Keywords: hyperplane added; hyerplane removed

comment:20 Changed 2 years ago by Jean-Philippe Labbé

Reviewers: Jean-Philippe Labbé, Travis Scrimshaw

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

comment:21 Changed 2 years ago by Travis Scrimshaw

Status: needs_reviewpositive_review

Yep, I am happy with it. Thanks.

comment:22 Changed 2 years ago by gh-kliem

Thank you.

comment:23 Changed 2 years ago by Volker Braun

Branch: public/29661-reb930217fb8b5979ee798da8fd8542048bc1e1a90f
Resolution: fixed
Status: positive_reviewclosed
Note: See TracTickets for help on using tickets.