Opened 3 years ago
Closed 2 years ago
#29661 closed enhancement (fixed)
Some optimizations for method regions of hyperplane arrangements
Reported by:  ghkliem  Owned by:  

Priority:  major  Milestone:  sage9.2 
Component:  geometry  Keywords:  hyperplane arrangements, regions 
Cc:  JeanPhilippe Labbé, Laith Rastanawi, Sophia Elia  Merged in:  
Authors:  Jonathan Kliem  Reviewers:  JeanPhilippe Labbé, Travis Scrimshaw 
Report Upstream:  N/A  Work issues:  
Branch:  930217f (Commits, GitHub, GitLab)  Commit:  930217fb8b5979ee798da8fd8542048bc1e1a90f 
Dependencies:  Stopgaps: 
Description
There are some easy optimzations that speed up the computation of regions of hyperplane arrangements:
 To see if an hyperplane intersects a polyhedron nontrivially 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 3 years ago by
Branch:  → public/29661 

Commit:  → ad162d70516726db61edec848c82cdef99863064 
Status:  new → needs_review 
comment:2 Changed 3 years ago by
Status:  needs_review → needs_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:
+ 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.
comment:3 followup: 4 Changed 3 years ago by
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 fulldimensional 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 Changed 3 years ago by
Replying to ghkliem:
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 fulldimensional 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 followup: 7 Changed 3 years ago by
Branch:  public/29661 → public/29661reb 

Commit:  ad162d70516726db61edec848c82cdef99863064 → b05df006ca52e7943b73df4b63124aad5b9654f8 
Status:  needs_work → needs_review 
comment:6 Changed 3 years ago by
Commit:  b05df006ca52e7943b73df4b63124aad5b9654f8 → 9a690e4b96331cf9e66128e8326d088d81fdf347 

Branch pushed to git repo; I updated commit sha1. New commits:
9a690e4  sorting entries by appearance year

comment:7 Changed 3 years ago by
Replying to ghkliem:
Btw, is there a way to obtain a ndimensional 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 followup: 9 Changed 3 years ago by
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 Changed 3 years ago by
Replying to ghkliem:
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 3 years ago by
Commit:  9a690e4b96331cf9e66128e8326d088d81fdf347 → 00103d0a6a3aeb0c1512a241316b8904d825de8e 

Branch pushed to git repo; I updated commit sha1. New commits:
00103d0  fix empty hyperplane arrangement

comment:11 followup: 13 Changed 3 years ago by
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:12 Changed 3 years ago by
Commit:  00103d0a6a3aeb0c1512a241316b8904d825de8e → f18fd2ab164f6abf95fd4e75219e17c85a5b75d6 

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

comment:13 followup: 14 Changed 3 years ago by
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 anyclause, 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 Changed 3 years ago by
Replying to ghkliem:
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 anyclause, 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 3 years ago by
Commit:  f18fd2ab164f6abf95fd4e75219e17c85a5b75d6 → 59dc0740b561563f81e8bfb0731ea5e90245da92 

Branch pushed to git repo; I updated commit sha1. New commits:
59dc074  backslash and alignment

comment:16 Changed 3 years ago by
Commit:  59dc0740b561563f81e8bfb0731ea5e90245da92 → a2da77530b455059013f3ab6201df939d662f921 

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

comment:17 Changed 3 years ago by
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 3 years ago by
Commit:  a2da77530b455059013f3ab6201df939d662f921 → 930217fb8b5979ee798da8fd8542048bc1e1a90f 

Branch pushed to git repo; I updated commit sha1. New commits:
930217f  alignment again

comment:19 Changed 3 years ago by
Keywords:  hyperplane added; hyerplane removed 

comment:20 Changed 3 years ago by
Reviewers:  → JeanPhilippe Labbé, Travis Scrimshaw 

It looks good to me. If Travis agrees, I suggest to give positive review.
comment:21 Changed 3 years ago by
Status:  needs_review → positive_review 

Yep, I am happy with it. Thanks.
comment:23 Changed 2 years ago by
Branch:  public/29661reb → 930217fb8b5979ee798da8fd8542048bc1e1a90f 

Resolution:  → fixed 
Status:  positive_review → closed 
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:
optimize region computation of hyperplane arrangements