Opened 2 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:  jipilab, ghLaisRast, selia  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 2 years ago by
 Branch set to public/29661
 Commit set to ad162d70516726db61edec848c82cdef99863064
 Status changed from new to needs_review
comment:2 Changed 2 years ago by
 Status changed from needs_review to 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 2 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 in reply to: ↑ 3 Changed 2 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 2 years ago by
 Branch changed from public/29661 to public/29661reb
 Commit changed from ad162d70516726db61edec848c82cdef99863064 to b05df006ca52e7943b73df4b63124aad5b9654f8
 Status changed from needs_work to needs_review
comment:6 Changed 2 years ago by
 Commit changed from b05df006ca52e7943b73df4b63124aad5b9654f8 to 9a690e4b96331cf9e66128e8326d088d81fdf347
Branch pushed to git repo; I updated commit sha1. New commits:
9a690e4  sorting entries by appearance year

comment:7 in reply to: ↑ 5 Changed 2 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 2 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 in reply to: ↑ 8 Changed 2 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 2 years ago by
 Commit changed from 9a690e4b96331cf9e66128e8326d088d81fdf347 to 00103d0a6a3aeb0c1512a241316b8904d825de8e
Branch pushed to git repo; I updated commit sha1. New commits:
00103d0  fix empty hyperplane arrangement

comment:11 followup: ↓ 13 Changed 2 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 2 years ago by
 Commit changed from 00103d0a6a3aeb0c1512a241316b8904d825de8e to f18fd2ab164f6abf95fd4e75219e17c85a5b75d6
Branch pushed to git repo; I updated commit sha1. New commits:
f18fd2a  small improvments

comment:13 in reply to: ↑ 11 ; followup: ↓ 14 Changed 2 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 in reply to: ↑ 13 Changed 2 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 2 years ago by
 Commit changed from f18fd2ab164f6abf95fd4e75219e17c85a5b75d6 to 59dc0740b561563f81e8bfb0731ea5e90245da92
Branch pushed to git repo; I updated commit sha1. New commits:
59dc074  backslash and alignment

comment:16 Changed 2 years ago by
 Commit changed from 59dc0740b561563f81e8bfb0731ea5e90245da92 to a2da77530b455059013f3ab6201df939d662f921
Branch pushed to git repo; I updated commit sha1. New commits:
a2da775  alignment

comment:17 Changed 2 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 2 years ago by
 Commit changed from a2da77530b455059013f3ab6201df939d662f921 to 930217fb8b5979ee798da8fd8542048bc1e1a90f
Branch pushed to git repo; I updated commit sha1. New commits:
930217f  alignment again

comment:19 Changed 2 years ago by
 Keywords hyperplane added; hyerplane removed
comment:20 Changed 2 years ago by
 Reviewers set to JeanPhilippe Labbé, Travis Scrimshaw
It looks good to me. If Travis agrees, I suggest to give positive review.
comment:21 Changed 2 years ago by
 Status changed from needs_review to positive_review
Yep, I am happy with it. Thanks.
comment:22 Changed 2 years ago by
Thank you.
comment:23 Changed 2 years ago by
 Branch changed from public/29661reb to 930217fb8b5979ee798da8fd8542048bc1e1a90f
 Resolution set to fixed
 Status changed from positive_review to 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