Opened 3 years ago

Closed 2 years ago

# Some optimizations for method regions of hyperplane arrangements

Reported by: Owned by: gh-kliem major sage-9.2 geometry hyperplane arrangements, regions Jean-Philippe Labbé, Laith Rastanawi, Sophia Elia Jonathan Kliem Jean-Philippe Labbé, Travis Scrimshaw N/A 930217f 930217fb8b5979ee798da8fd8542048bc1e1a90f

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

### comment:1 Changed 3 years ago by gh-kliem

Branch: → public/29661 → ad162d70516726db61edec848c82cdef99863064 new → 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:

 ​ad162d7 `optimize region computation of hyperplane arrangements`

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

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.

Last edited 3 years ago by Jean-Philippe Labbé (previous) (diff)

### comment:3 follow-up:  4 Changed 3 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 3 years ago by Jean-Philippe Labbé

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 3 years ago by gh-kliem

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

New commits:

 ​79e5f80 `optimize region computation of hyperplane arrangements` ​61da5b6 `readability of code` ​b05df00 `added example`

### comment:6 Changed 3 years ago by git

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

 ​9a690e4 `sorting entries by appearance year`

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

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 3 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 3 years ago by Jean-Philippe Labbé

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 git

Commit: 9a690e4b96331cf9e66128e8326d088d81fdf347 → 00103d0a6a3aeb0c1512a241316b8904d825de8e

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

 ​00103d0 `fix empty hyperplane arrangement`

### comment:11 follow-up:  13 Changed 3 years ago by Travis Scrimshaw

```-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 3 years ago by Travis Scrimshaw (previous) (diff)

### comment:12 Changed 3 years ago by git

Commit: 00103d0a6a3aeb0c1512a241316b8904d825de8e → 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 3 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?

```-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 3 years ago by Travis Scrimshaw

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 3 years ago by git

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 git

Commit: 59dc0740b561563f81e8bfb0731ea5e90245da92 → a2da77530b455059013f3ab6201df939d662f921

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

 ​a2da775 `alignment`

### comment:17 Changed 3 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 3 years ago by git

Commit: a2da77530b455059013f3ab6201df939d662f921 → 930217fb8b5979ee798da8fd8542048bc1e1a90f

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

 ​930217f `alignment again`

### comment:20 Changed 3 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 3 years ago by Travis Scrimshaw

Status: needs_review → positive_review

Yep, I am happy with it. Thanks.

Thank you.

### comment:23 Changed 2 years ago by Volker Braun

Branch: public/29661-reb → 930217fb8b5979ee798da8fd8542048bc1e1a90f → fixed positive_review → closed
Note: See TracTickets for help on using tickets.