add orthogonal/orthonormal feature to affine_hull of polyhedra

[mo271]

As mentioned in #16045 comment:22 it would be desirable to get the affine hull of a polyhderon, while preserving the polyhedron isometrically (instead of applying an affine transformation, that might distort lenght, volumes, angels, etc.).

I propose to add an option "algorithm" for the .affine_hull method for polyhedra. This can then be used to correctly calculate volume of non full-diemsional polyhedra, see #16045.

CC: @videlec

Component: geometry

Keywords: polyhedron

Author: Moritz Firsching

Branch: 2ca1a84

Reviewer: Matthias Koeppe, Vincent Delecroix

Issue created by migration from https://trac.sagemath.org/ticket/22804

[mo271]

Branch: u/moritz/affine_hull

[videlec]

comment:2

Two quick comments

1) The keyword algorithm is very misleading. "Algorithm" refers to a method to compute something not the result itself. Moreover, both constructions are projections. You should find something more coherent! And also describe in the documentation what is the result.

2) I like better the more explicit

P.vertices()[0].vector() == P.ambient_space().zero()

---

New commits:

91f239b

adding isometry feature to affine_hull

[videlec]

Commit: 91f239b

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

43dd759

change name of parameter to 'force_isometry'

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 91f239b to 43dd759

[mo271]

comment:4

Thanks for your comments, Vincent! I hope the new commit addresses all the proposed changes.

I was not quite happy with the choice of "algorithm" for a name of the parameter. Now I named the parameter "force_isometry".

[videlec]

comment:5

I don't understand why this step

    #if the self is full-dimensional, return it immediately
    if self.ambient_dim() == self.dim():
        return self

is inside the if force_isometry.

[mo271]

comment:6

Replying to @videlec:

I don't understand why this step

    #if the self is full-dimensional, return it immediately
    if self.ambient_dim() == self.dim():
        return self

is inside the if force_isometry.

The intent was not to disturb the original method too much.

Howeer , I will gladly move it above the if force_isometry. (I don't know why the old method didn't include this check..)

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 43dd759 to 7afd9aa

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

7afd9aa

move dimenion check to top of method

[mkoeppe]

comment:8

This looks like a nice interface, but let me just say that for volume computations it's overkill to compute this (which will lead to a double description computation in the slow "field" backend), when all one needs is to compute a scaling factor for the volume...

Another comment: I don't think this is a suitable test:

if self.base_ring() in [ZZ,QQ]

because you could as well have some number field that is too small.

sage: D = polytopes.dodecahedron()
sage: F = D.faces(2)[0].as_polyhedron()
sage: F.affine_hull(force_isometry=True)
TypeError: QR decomposition unable to compute square roots in Number Field in sqrt5 with defining polynomial x^2 - 5

[mo271]

comment:9

I agree that for volume computations, this is an overkill. (see more on that below)

For other purposes it might still be useful to have. At the moment, it is not possible it is cumbersome to plot a regular tetrahedron, which looks like a regular tetrahedron. With this feature this is really easy. (Same for other polytopes that are most easily defined as non full-dimensional polytopes. Therefore I would like to keep this new feature.

I propose to add two new ways to obtain the affine hull:

• orthonormal, which should be simply what has been called isometry above
• orthogonal, a scaled version, that uses gram_schmitt. From that the scaling factor can be obtained.

Now I wonder where to store or return the scaling factor.

Would the following be a reasonable way to do it: have a method like

def affine_hull(self, mode= "default", return_scaling_factor=False)

Here mode can take values "default", "orthogonal", "orthonormal". If return_scaling_factor is True, then instead of returning simply the polyhedron, return a pair (polyhedron, scaling factor).

I am not quite happy with the names "mode" and "return_scaling_factor", I am sure we can find better alternatives.

What do you think?

[mkoeppe]

comment:10

How about returning the map that sends the result back to the space, instead of returning the factor? Then document how the volume transforms, with an example.

Instead of mode, I guess one could just use two boolean keyword parameters.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

c9425eb	adding isometry feature to affine_hull
8d23a22	change name of parameter to 'force_isometry'
83e5277	move dimension check to top of method
a7e3a0b	introduce parameter 'orthogonal', 'orthonormal' and 'as_affine_map'
f0cac6b	typos and one more doctest

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 7afd9aa to f0cac6b

[mo271]

comment:12

Thank you very much for your quick reply, mkoeppe!

I have now indeed introduced two booleans, orthogonal and orthonormal. "Returning the map" means in this case returning an affine map, I guess. I went ahead and return now a "linear_transformation" together with a vector. I guess there is some support of affine groups, but I am not sure if this should be used instead: http://doc.sagemath.org/html/en/reference/groups/sage/groups/affine_gps/group_element.html.

In any case, the method is all set to use it with the volume computations now.

As a side effect, a user will know be able to plot nice 3d-plots of the regular Tetrahedron as well as the Permutahedron. (Once this ticket gets a positive review, I suggest to alter the plot routine to enforce an orthonormal transformation by default)

[videlec]

comment:15

• a orthogonal -> an orthogonal
• a orthonormal -> an orthonormal
• orhtonormal -> orthonormal
• I think that the fact that orthonormal=True/orthogonal=True will silently modify the base ring should be mentioned in the INPUT. One also might actually want to prevent this kind of behavior. What about an additional boolean keyword extend that would be inserted as follows in the code

try:
    A = M.gram_schmidt(orthonormal=orthonormal)[0]
except TypeError:
+    if not extend:
+        raise
    M = matrix(AA, M)
    A = M.gram_schmidt(orthonormal=orthonormal)[0]

• you should test this with polyhedron over embedded number fields like sage: K.<sqrt3> = QuadraticField(3)

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

6304fef	adding isometry feature to affine_hull
c8638fe	change name of parameter to 'force_isometry'
6e6c68b	move dimension check to top of method
bef4cf9	introduce parameter 'orthogonal', 'orthonormal' and 'as_affine_map'
caa089b	typos and one more doctest
ded4ac6	added base_extend parameter, plus a few typos

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from f0cac6b to ded4ac6

[mo271]

comment:17

Thank you Vincent for the comments! I have implemented all of them. I called the parameter base_extend instead of extend.

Also its rebased on the current beta.

[videlec]

comment:18

Please change to extend. It is what is mostly used in Sage

sage: V=QQ^2
sage: H=V.endomorphism_ring()([[0,-1],[1,0]])
sage: H.eigenvalues()
[-1*I, 1*I]
sage: H.eigenvalues(extend=False)
[]

Or

sage: 3.sqrt()
sqrt(3)
sage: 3.sqrt(extend=False)
Traceback (most recent call last)
...
ValueError: square root of 3 not an integer

If you want to change conventions, that is the purpose of another ticket.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from ded4ac6 to 87c86b5

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

87c86b5

changed 'base_extend' to 'extend'

[mo271]

comment:20

Replying to @videlec:

Please change to extend. It is what is mostly used in Sage

done; I didn't mean to alter the conventions.

[videlec]

comment:21

Could you update the NOTE as the base ring can not be changed by default. Or it can still be the case the base ring changes silently?

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 87c86b5 to 2ca1a84

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

0aaab2b	adding isometry feature to affine_hull
17ac7c3	change name of parameter to 'force_isometry'
777ce7e	move dimension check to top of method
836be61	introduce parameter 'orthogonal', 'orthonormal' and 'as_affine_map'
1679534	typos and one more doctest
1f64046	added base_extend parameter, plus a few typos
ee003eb	changed 'base_extend' to 'extend'
2ca1a84	removed NOTE and made sure base ring stays fixed

[mo271]

comment:23

Thanks again, for looking at the ticket, Vincent!

I have now updated the code to make sure that the base ring is never changed (and removed the note).

Before that it was never extended, but sometimes "simplified", i.e. restricted to a smaller base ring, if you were lucky. See the following example:

sage: K.<sqrt2> = QuadraticField(2)
sage: sage: P = Polyhedron([2*[K.zero()],2*[sqrt2]]); P
A 1-dimensional polyhedron in (Number Field in sqrt2 with defining polynomial x^2 - 2)^2 defined as the convex hull of 2 vertices
sage: P.vertices()
(A vertex at (0, 0), A vertex at (sqrt2, sqrt2))

previously we had:

A = P.affine_hull(orthonormal=True); A
A 1-dimensional polyhedron in ZZ^1 defined as the convex hull of 2 vertices
sage: A.vertices()
(A vertex at (0), A vertex at (2))

So the base_ring was changed from K to ZZ.

Now I changed this to:

sage: A = P.affine_hull(orthonormal=True); A
A 1-dimensional polyhedron in (Number Field in sqrt2 with defining polynomial x^2 - 2)^1 defined as the convex hull of 2 vertices
sage: A.vertices()
(A vertex at (0), A vertex at (2))

also, I rebased on the current beta...

I think there are advantages to both, simplifying the base ring and also to keeping it. It could also be controlled by the extend parameter, but the name doesn't quite capture what is going on. But in general I don't think this should be done here, but there should be a general method, which tries to find the smallest base ring possible. (This might be usefull for demoting rational polytopes to integer ones, but also from AA to some smaller field.)

[videlec]

comment:24

I feel it is good enough. I just launched the patchbot to see how it goes.

Concerning ring/field simplification there should be something at the level of fields. Such method would also be useful on matrices for example.

[videlec]

Reviewer: Matthias Koeppe, Vincent Delecroix

[videlec]

comment:25

I think that the ticket is good enough (and patchbot is happy)

[vbraun]

Changed branch from u/moritz/affine_hull to 2ca1a84

[mo271]

Changed commit from 2ca1a84 to none

[mo271]

comment:27

Thanks for the review, Vincent! I will continue to work on #16045.