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

[tgaona]

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

[tgaona]

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

[jdemeyer]

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

[tgaona]

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]

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]

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]

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]

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

[tgaona]

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]

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]

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.

[tgaona]

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

[git]

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

[git]

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

​8a040f3

Merge branch 'develop' into ticket/19112

[git]

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

[jdemeyer]

Can you review #18669?

[jdemeyer]

The patchbot complains about a TAB character.

[git]

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

​d2b3dde

Removed errant TAB character in quadratic_form.py

[jdemeyer]

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]

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

[jdemeyer]

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]

Do you really need the copy in the function isometry?

[jdemeyer]

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]

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

[git]

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

​4b22c51

Added doctesting to helper methods of isometry

[git]

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

​582ae45

Minor bugfixes/refactoring on isometry

[tgaona]

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]

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]

Replying to 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.

[annahaensch]

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.

[git]

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]

Sorry, but 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]

Merge conflicts with Sage 7.0

[git]

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

​56b5225

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

[git]

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

[git]

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

​642df10

Fixed typo in documentation example for _diagonal_isometry

[jdemeyer]

Why this added line?

+            v = parilist[i]

[jdemeyer]

I am still missing some high-level explanation of the algorithm. For example, what is FM?

[annahaensch]

I understand the algorithm, I don't think it's so complicated. It just looks for a vector w in W of length corresponding to, V[0][0], the upper left most entry in V. Then it simply extends the basis of W to contain w, preforms Gram Schmidt (fixing w), and then picks off <Q(w)> as an orthogonal component of W and <V[0][0]> as an orthogonal component of V. I feel satisfied the the algorithm does what it claims to do, and that it is mathematically correct.

A few notes on the documentation:

1. In your documentation for _diagonal_isometry() you call the function isometry() (of course it works out to be the same thing here, since the forms are diagonal to begin with.

2. Your documentation for _gram_schmidt() is a bit misleading, since it looks like you are treating the columns of a Gram matrix as a set of vectors to be orthogonalized, rather, do the process on some sort of (preferably basis) matrix.

In response to comment 45, I'm also not sure why v=parilist[i] got added. Is there a reason you're touching short_vector_list_up_to_length()? I think something may have happened on the merge.

[jdemeyer]

Replying to annahaensch:

I understand the algorithm

The fact that you understand it, is not sufficient. It should be documented in the code such that everybody can understand it.

[git]

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

​2b8f694	removed old line from quadratic_form__automorphisms
​70ed429	Fixed documentation errors and added more comments to _diagonal_isometry

[git]

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

​44c29ef	Merge branch 'master' into isometry2
​0b88414	Changed 'from quadratic_form' to 'from sage.quadratic_forms.quadratic_form'

[git]

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

​2b394f4

Removed trailing whitespace.

[sbrandhorst]

Thank you for writing this code. It is well documented and easy to read. The basic algorithm works. Good job! After computing some examples I noticed that coefficients tend to explode and examples beyond n=4 take ages. Is this expected? Is there an algorithm avoiding coefficient blowups?

I can review it. Here are some comments:

• use imperative in the docstrings: "Computes something". Should rather be: "Compute something"
• There is no need to introduce a new method isometry(). Instead give is_rationally_isometric a new argument return_matrix. You can model analogous behavior as in is_globally_equivalent_to.

• In the current implementation you do not check the base_ring is ZZ or QQ. This should be done. Or does it work more generally? For example the code does not work over the p-adics.
• You do not need a new method_compute_gram_matrix_from_basis. Instead you can just use Q(T).Gram_matrix_rational(). Here Q is a quadratic form and T a transformation matrix.
• make sure degenerate quadratic forms work as well and doctest it.
• add a random doctest. For example something in the direction of:

  sage: n = 4
  sage: G = matrix.random(QQ,n) #random
  sage: G = G + G.transpose()
  sage: Q = QuadraticForm(QQ,G)
  sage: T = matrix.random(QQ,n)
  sage: S = Q.isometry(Q(T))
  sage: Q == Q(S)
  True
  

[dimpase]

the patchbot results look OK (apart from pyflake plugin, which might tell you how to improve the code).

[embray]

Tickets still needing working or clarification should be moved to the next release milestone at the soonest (please feel free to revert if you think the ticket is close to being resolved).

[embray]

Ticket retargeted after milestone closed

[mkoeppe]

Batch modifying tickets that will likely not be ready for 9.1, based on a review of the ticket title, branch/review status, and last modification date.