Add a function "isometry" to the quadratic forms package.

[264bb60e-7393-4731-bda3-074974467a5a]

Adds a function "isometry" that returns an isometry from one rational quadratic form to another, provided that it exists.

CC: @annahaensch @sagetrac-tgaona

Component: quadratic forms

Keywords: isometry

Author: Tyler Gaona, Simon Brandhorst

Branch: ffa5295

Reviewer: Simon Brandhorst

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

[264bb60e-7393-4731-bda3-074974467a5a]

Branch: u/tgaona/ticket/19112

[264bb60e-7393-4731-bda3-074974467a5a]

New commits:

0d32f4b	Implements a method to compute isometries between positive definite quadratic forms.
9e211fe	Adds a flag parameter to short_vector_list_up_to_length to pass to PARI's

[264bb60e-7393-4731-bda3-074974467a5a]

Commit: 9e211fe

[jdemeyer]

Changed work issues from needs implementation to none

[jdemeyer]

comment:3

This looks like it reinvents the already-existing method is_globally_equivalent_to. What does your code add?

[264bb60e-7393-4731-bda3-074974467a5a]

comment:4

Hi Jeroen,

Sorry for the slow response. This method should be able to return isometries for forms over fields, where I believe is_globally_equivalent_to just works over the ring of integers. For example:

sage: Q = DiagonalQuadraticForm(QQ, [1, 5])
sage: F = QuadraticForm(QQ, 2, [1, 12, 81])
sage: Q.is_globally_equivalent_to(F)
False
sage: Q.is_rationally_isometric(F)
True
sage: T = Q.isometry(F)
sage: T
[  1  -2]
[  0 1/3]
sage: Q.Gram_matrix() == T.transpose() * F.Gram_matrix() * T
True

So Q is equivalent to F over the rationals, but is_globally_equivalent_to doesn't recognize this. This method was intended to complement the is_rationally_isometric method, so perhaps it should be refactored to work in a similar manner as is_globally_equivalent_to, i.e, adding an optional flag to is_rationally_isometric to return the transition matrix for the isometry.

[jdemeyer]

comment:5

Right, I think it's better to mention this more clearly in the documentation (I always get confused when quadratic-form people talk about isomorphic/isometric/equivalent, I never know what they mean).

[jdemeyer]

comment:6

That being said, I think you are really over-complicating the algorithm.

Suppose one of the forms is given by the diagonal 5*x0^2 + .... Then really all you need to do in the first step is to find a vector v of "length" 5. This can be done with some algorithm based on PARI's qfsolve(): the only issue is that qfsolve() finds a projective point, but you really want an affine point.

So I suggest to do the following:

1. implement an affine version of qfsolve() to solve equations of the form Q(x) = C where Q is a quadratic form and C is a constant. Do this in a separate Sage ticket.
2. use this new function to implement isometry on this ticket.

This will have several major advantages:

1. your code will be a lot faster.
2. it will no longer run forever if there is isometry (which is unacceptable).
3. you can remove the assumption that your quadratic form is positive-definite.

[jdemeyer]

comment:7

Also, your branch is based on an old version of Sage. You should really develop from the latest beta version (currently 6.10.beta2).

[jdemeyer]

comment:8

Please follow http://doc.sagemath.org/html/en/developer/coding_basics.html#documentation-strings for the correct formatting of docstrings.

[264bb60e-7393-4731-bda3-074974467a5a]

comment:10

I will make it more clear in the documentation that the method returns a transformation matrix, as opposed to saying it returns an isometry.

As to your second point, I agree that a cleaner, more efficient replacement for the first step of the algorithm is very desirable. However, I have looked at the source for PARI's qfsolve() and it isn't clear to me how to adapt it to return an affine vector x such that Q(x) = C. I would appreciate it if you could offer more guidance here, otherwise, I'm not sure how to proceed.

2. it will no longer run forever if there is isometry (which is unacceptable).

I can fix this by throwing an exception when is_rationally_isometric() returns false for the two forms. The reason I didn't include this initially is that I was working off of Sage 6.8 but is_rationally_isometric() was added in 6.9.

I will fix the formatting in the docstring.

[jdemeyer]

Dependencies: #18669

[jdemeyer]

comment:11

Given that your new functionality relies on rational_diagonal_form(), it would be a lot better in fact to base your patch on top of #18669.

[jdemeyer]

comment:12

About solving quadratic forms:

Suppose you want to solve Q(x) = c. Consider the quadratic form Q(x) - c*z^2 = 0, where z is an extra variable. Find a solution (x,z) to this quadratic form using qfsolve().

Case 1: If z != 0, then Q(x/z) = c.

Case 2: We found a solution Q(x) = 0. Let e be any vector such that B(x,e) != 0, where B is the bilinear form corresponding to Q. To find e, just try all unit vectors (0,..0,1,0...0). Let a = (c - Q(e))/B(x,e) and let y = e + a*x. Then Q(y) = c.

It would be great if you could implement this on a separate ticket.

[264bb60e-7393-4731-bda3-074974467a5a]

Changed dependencies from #18669 to #18669 #19533

[264bb60e-7393-4731-bda3-074974467a5a]

comment:13

Thanks for the clear explanation. I've opened up a separate ticket for this, and will address the various changes you've suggested.

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

Changed commit from 9e211fe to afcd21b

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

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

51ba9f8	Adds isometry method for quadratic forms.
afcd21b	Merge branch 'research' into ticket/19112

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

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

8a040f3

Merge branch 'develop' into ticket/19112

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

Changed commit from afcd21b to 8a040f3

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

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

3ee588d	Use PARI to diagonalize quadratic forms
2689325	Merge tag '6.10.beta1' into t/18669/use_pari_to_compute_rational_diagonal_form__
45cbc2d	Further fixes to rational diagonal form
ee778c8	Implements a method to compute isometries between positive definite quadratic forms.
27b2190	Adds a flag parameter to short_vector_list_up_to_length to pass to PARI's
92de9ef	Adds isometry method for quadratic forms.
1b1a574	Merge branch 'ticket/19533' into ticket/19112
2f8ebc7	Bug fixes in 'isometry' method and refactoring to utilize the 'solve' method for quadratic forms
8f8ce15	Merge branch 'u/tgaona/ticket/19112' of git://trac.sagemath.org/sage into ticket/19112

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

Changed commit from 8a040f3 to 8f8ce15

[jdemeyer]

comment:18

Can you review #18669?

[jdemeyer]

comment:19

The patchbot complains about a TAB character.

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

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

d2b3dde

Removed errant TAB character in quadratic_form.py

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

Changed commit from 8f8ce15 to d2b3dde

[jdemeyer]

Changed dependencies from #18669 #19533 to none

[jdemeyer]

comment:22

I haven't looked closely at the code, but some clean-up should be done:

1. the function diagonal_isometry should have doctests
2. remove commented code like #print ...

[jdemeyer]

comment:23

Why does the final term need to be a special case? I don't like special cases unless they are justified.

[jdemeyer]

comment:24

Are you sure it's worth to check is_diagonal? Is there significant gain?

Are you sure it's worth to check is_rationally_isometric? I would just try to compute the isometry and let it fail if it doesn't exist.

[jdemeyer]

comment:25

Do you really need the copy in the function isometry?

[jdemeyer]

comment:26

Can you explain this block of code:

+            # Find a vector w such that Q(v) = F(w) where v = [1, ..., 0]
+            # v, w = vectors_of_common_length_dev(Q, F, q_basis, f_basis, i)
+            v = vector([0] * (n - i))
+            index = 0;
+            while True:
+                v[index] = v[index] + 1
+                index = (index + 1) % (n - i)
+                c = Q(v)
+                try:
+                    w = F.solve(c)
+                    #print("Find vectors {0} and {1} such that Q(v) = F(w)".format(v, w))        
+                    if not zero_row(f_basis, w, i) and not zero_row(q_basis, v, i):
+                        break
+                except ArithmeticError:
+                    # No solution found, try another vector.
+                    pass

[jdemeyer]

comment:28

Can you move the helper functions out of the diagonal_isometry function and also doctest them?

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

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

4b22c51

Added doctesting to helper methods of isometry

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

Changed commit from d2b3dde to 4b22c51

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

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

582ae45

Minor bugfixes/refactoring on isometry

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

Changed commit from 4b22c51 to 582ae45

[264bb60e-7393-4731-bda3-074974467a5a]

comment:31

Replying to @jdemeyer:

Why does the final term need to be a special case? I don't like special cases unless they are justified.

It was originally necessary because the short_vector_list_up_to_length() function didn't work for 1-dimensional quadratic forms, and the final iteration of the algorithm operates on two 1-dimensional forms. However, since that's been replaced by the PARI method, it's no longer necessary to have it as a special case, so thanks for pointing that out.

Are you sure it's worth to check is_diagonal? Is there significant gain?

I suppose not, and simplifies the code to remove it, so I will.

Are you sure it's worth to check is_rationally_isometric?

I believe this is necessary. The loop that looks for two vectors such that the modified bases including them will be non-singular will try vectors indefinitely until it finds a satisfactory pair. If the forms aren't equivalent there's no guarantee such a pair will be found.

Do you really need the copy in the function isometry?

Nope. That's a relic from an early version when isometry and diagonal_isometry were one function. I'll remove it, thanks for catching that.

Can you explain this block of code:

+            # Find a vector w such that Q(v) = F(w) where v = [1, ..., 0]
+            # v, w = vectors_of_common_length_dev(Q, F, q_basis, f_basis, i)
+            v = vector([0] * (n - i))
+            index = 0;
+            while True:
+                v[index] = v[index] + 1
+                index = (index + 1) % (n - i)
+                c = Q(v)
+                try:
+                    w = F.solve(c)
+                    #print("Find vectors {0} and {1} such that Q(v) = F(w)".format(v, w))        
+                    if not zero_row(f_basis, w, i) and not zero_row(q_basis, v, i):
+                        break
+                except ArithmeticError:
+                    # No solution found, try another vector.
+                    pass

This block finds a pair of vectors v and w such that Q(v) == F(v). These vectors will represent a linear combination of the vectors in the basis for each quadratic form. It's necessary that modifying the bases to include these vectors not produce a basis whose matrix is singular. So essentially this loop just looks for a pair of vectors satsifying these properties. It starts with v = [1, 0, 0] (for a 3-dimensional form) and finds w by calling F.solve(Q(v)). If this pair doesn't work, it finds a new v and starts over. I get new v's by incrementing each term in the vector so the first few vectors that are generated are: [1, 0, 0], [1, 1, 0], [1, 1, 1], [2, 1, 1]...

Also, I realized that matrices have an is_singular function, so I'm getting rid of the zero_row function.

---

New commits:

582ae45

Minor bugfixes/refactoring on isometry

[jdemeyer]

comment:33

More comments later, but already this: you make changes to short_vector_list_up_to_length() which I think do not belong to this ticket. In fact, it looks more likely that they belong to #14868. So I suggest you undo those changes here and move them to a new branch on #14868.

[jdemeyer]

comment:34

Replying to @sagetrac-tgaona:

This block finds a pair of vectors v and w such that Q(v) == F(v). These vectors will represent a linear combination of the vectors in the basis for each quadratic form. It's necessary that modifying the bases to include these vectors not produce a basis whose matrix is singular. So essentially this loop just looks for a pair of vectors satsifying these properties. It starts with v = [1, 0, 0] (for a 3-dimensional form) and finds w by calling F.solve(Q(v)). If this pair doesn't work, it finds a new v and starts over. I get new v's by incrementing each term in the vector so the first few vectors that are generated are: [1, 0, 0], [1, 1, 0], [1, 1, 1], [2, 1, 1]...

Sorry, I still don't understand the algorithm. You really need more comments in the code. In particular, you need to explain the purpose of the variables i, q_basis, f_basis, qb and fb.

My feeling is that the current algorithm is too complicated (I don't think you need to loop for v), but I cannot really say how to improve it since I don't understand it.

[4d3919ea-19e6-493c-bf47-650e092876a7]

comment:36

I have looked over tgaona's algorithm, and it makes sense to me. The block of code mentioned in comment 26 is a bit confusing, but I believe that it is necessary since its eventually necessary to invert, and therefore avoid the possibility of having a singular change of basis matrix.

Another source of confusion, I believe, is in the description of the helper function modify basis. I do believe that it does precisely what it needs to do in the context of the code, but the description is not quite correct, and a little bit misleading, since it is always a function acting on a submatrix of the original basis matrix. A better description might read "Given a lattice L with basis matrix M and a vector, v=(v_1,...,v_n) of length n, this function extends the basis {b_1,...,b_n} of an underlying nxn orthogonal component of L to contain the vector v_1b_1+...+v_nb_n."

The functions diagonal_isometry, compute_gram_matrix_from_basis, modify_basis, and graham_schmidt are all very particular to this parent function isometry, so (I believe) the standard Sage convention is to begin those with an underscore. Also, I believe a better name for the parent function would be explicit_isometry just to differentiate this from the binary function returning True or False.

Other than that, I'll concede that there may very well be a faster way to run this algorithm, but as tgaona has it, it's certainly correct.

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

Changed commit from 582ae45 to 17d12c6

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

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

6d9be4f	removed changes to 'short_vector_list_up_to_length()'
17d12c6	Prefixed helper functions with an underscore and reworded documentation.

[jdemeyer]

comment:39

Sorry, but [comment:34] hasn't been addressed yet. The code is simply too hard to understand and therefore impossible to review for me.

And why did you rename isometry -> explicit_isometry? That's much harder to find.

[jdemeyer]

comment:40

Merge conflicts with Sage 7.0

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

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

56b5225

Merge branch 'develop' of git://trac.sagemath.org/sage into ticket/19112

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

Changed commit from 17d12c6 to 56b5225

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

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

132da46	attempt to simplify algorithm
24e01fe	Simplified algorithm and code for 'isometry' and its helpers
a56cb96	Merge branch 'u/tgaona/ticket/19112' of git://trac.sagemath.org/sage into isometry2

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

Changed commit from 56b5225 to a56cb96

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

Changed commit from a56cb96 to 642df10