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: | jipilab, gh-LaisRast, selia | 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: |
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
- 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 follow-up: ↓ 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 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
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 follow-up: ↓ 7 Changed 2 years ago by
- Branch changed from public/29661 to public/29661-reb
- 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 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 follow-up: ↓ 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 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
- Commit changed from 9a690e4b96331cf9e66128e8326d088d81fdf347 to 00103d0a6a3aeb0c1512a241316b8904d825de8e
Branch pushed to git repo; I updated commit sha1. New commits:
00103d0 | fix empty hyperplane arrangement
|
comment:11 follow-up: ↓ 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 ; follow-up: ↓ 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 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
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
- 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 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
- 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/29661-reb 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