#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:

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 15 months ago by gh-kliem

  • Branch set to public/29661
  • Commit set to ad162d70516726db61edec848c82cdef99863064
  • Status changed from new to needs_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 15 months ago by jipilab

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

Last edited 15 months ago by jipilab (previous) (diff)

comment:3 follow-up: Changed 15 months 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 15 months ago by 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.

comment:5 follow-up: Changed 15 months ago by gh-kliem

  • Branch changed from public/29661 to public/29661-reb
  • Commit changed from ad162d70516726db61edec848c82cdef99863064 to b05df006ca52e7943b73df4b63124aad5b9654f8
  • Status changed from needs_work to needs_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 15 months ago by git

  • Commit changed from b05df006ca52e7943b73df4b63124aad5b9654f8 to 9a690e4b96331cf9e66128e8326d088d81fdf347

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

9a690e4sorting entries by appearance year

comment:7 in reply to: ↑ 5 Changed 15 months ago by 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?

comment:8 follow-up: Changed 15 months 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 15 months ago by 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...

comment:10 Changed 15 months ago by git

  • Commit changed from 9a690e4b96331cf9e66128e8326d088d81fdf347 to 00103d0a6a3aeb0c1512a241316b8904d825de8e

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

00103d0fix empty hyperplane arrangement

comment:11 follow-up: Changed 15 months ago by 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
Last edited 15 months ago by tscrim (previous) (diff)

comment:12 Changed 15 months ago by git

  • Commit changed from 00103d0a6a3aeb0c1512a241316b8904d825de8e to f18fd2ab164f6abf95fd4e75219e17c85a5b75d6

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

f18fd2asmall improvments

comment:13 in reply to: ↑ 11 ; follow-up: Changed 15 months 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 15 months ago by 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.

comment:15 Changed 15 months ago by git

  • Commit changed from f18fd2ab164f6abf95fd4e75219e17c85a5b75d6 to 59dc0740b561563f81e8bfb0731ea5e90245da92

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

59dc074backslash and alignment

comment:16 Changed 15 months ago by git

  • Commit changed from 59dc0740b561563f81e8bfb0731ea5e90245da92 to a2da77530b455059013f3ab6201df939d662f921

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

a2da775alignment

comment:17 Changed 15 months 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 15 months ago by git

  • Commit changed from a2da77530b455059013f3ab6201df939d662f921 to 930217fb8b5979ee798da8fd8542048bc1e1a90f

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

930217falignment again

comment:19 Changed 15 months ago by jipilab

  • Keywords hyperplane added; hyerplane removed

comment:20 Changed 15 months ago by jipilab

  • 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 15 months ago by tscrim

  • Status changed from needs_review to positive_review

Yep, I am happy with it. Thanks.

comment:22 Changed 15 months ago by gh-kliem

Thank you.

comment:23 Changed 14 months ago by vbraun

  • Branch changed from public/29661-reb to 930217fb8b5979ee798da8fd8542048bc1e1a90f
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.