Opened 2 years ago
Last modified 5 weeks ago
#19112 needs_review enhancement
Add a function "isometry" to the quadratic forms package.
Reported by:  tgaona  Owned by:  

Priority:  minor  Milestone:  sage7.1 
Component:  quadratic forms  Keywords:  isometry 
Cc:  annahaensch, tgaona  Merged in:  
Authors:  Tyler Gaona  Reviewers:  Simon Brandhorst 
Report Upstream:  N/A  Work issues:  
Branch:  u/tgaona/ticket/19112 (Commits)  Commit:  2b394f41bd18eb22c6558a1441444692b79caec7 
Dependencies:  Stopgaps: 
Description
Adds a function "isometry" that returns an isometry from one quadratic form to another, provided that it exists.
Change History (53)
comment:1 Changed 2 years ago by
 Branch set to u/tgaona/ticket/19112
comment:2 Changed 2 years ago by
 Commit set to 9e211fecd4b525ad30276d37d442a68ce9bec8f0
 Status changed from new to needs_review
comment:3 Changed 2 years ago by
 Milestone changed from sage6.9 to sage6.10
 Status changed from needs_review to needs_info
 Work issues needs implementation deleted
This looks like it reinvents the alreadyexisting method is_globally_equivalent_to
. What does your code add?
comment:4 Changed 2 years ago by
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.
comment:5 Changed 2 years ago by
Right, I think it's better to mention this more clearly in the documentation (I always get confused when quadraticform people talk about isomorphic/isometric/equivalent, I never know what they mean).
comment:6 Changed 2 years ago by
That being said, I think you are really overcomplicating 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:
 implement an affine version of
qfsolve()
to solve equations of the formQ(x) = C
whereQ
is a quadratic form andC
is a constant. Do this in a separate Sage ticket.  use this new function to implement
isometry
on this ticket.
This will have several major advantages:
 your code will be a lot faster.
 it will no longer run forever if there is isometry (which is unacceptable).
 you can remove the assumption that your quadratic form is positivedefinite.
comment:7 Changed 2 years ago by
Also, your branch is based on an old version of Sage. You should really develop from the latest beta version (currently 6.10.beta2
).
comment:8 Changed 2 years ago by
Please follow http://doc.sagemath.org/html/en/developer/coding_basics.html#documentationstrings for the correct formatting of docstrings.
comment:9 Changed 2 years ago by
 Status changed from needs_info to needs_work
comment:10 Changed 2 years ago by
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.
 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.
comment:11 Changed 2 years ago by
 Dependencies set to #18669
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.
comment:12 Changed 2 years ago by
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.
comment:13 Changed 2 years ago by
 Dependencies changed from #18669 to #18669 #19533
Thanks for the clear explanation. I've opened up a separate ticket for this, and will address the various changes you've suggested.
comment:14 Changed 2 years ago by
 Commit changed from 9e211fecd4b525ad30276d37d442a68ce9bec8f0 to afcd21b014a9c543b02c0b957c6d2c028d60b17f
comment:15 Changed 2 years ago by
 Commit changed from afcd21b014a9c543b02c0b957c6d2c028d60b17f to 8a040f3e2bdc5b312bd14b1fe36e85dda87ead91
Branch pushed to git repo; I updated commit sha1. New commits:
8a040f3  Merge branch 'develop' into ticket/19112

comment:16 Changed 2 years ago by
 Commit changed from 8a040f3e2bdc5b312bd14b1fe36e85dda87ead91 to 8f8ce1575dcfb53fee524a9e546d297d22aa82a0
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

comment:17 Changed 2 years ago by
 Status changed from needs_work to needs_review
comment:18 Changed 2 years ago by
Can you review #18669?
comment:19 Changed 2 years ago by
The patchbot complains about a TAB character.
comment:20 Changed 2 years ago by
 Commit changed from 8f8ce1575dcfb53fee524a9e546d297d22aa82a0 to d2b3dde11edb22ba199d4666b99ad753865c5322
Branch pushed to git repo; I updated commit sha1. New commits:
d2b3dde  Removed errant TAB character in quadratic_form.py

comment:21 Changed 2 years ago by
 Dependencies #18669 #19533 deleted
comment:22 Changed 2 years ago by
I haven't looked closely at the code, but some cleanup should be done:
 the function
diagonal_isometry
should have doctests  remove commented code like
#print ...
comment:23 followup: ↓ 31 Changed 2 years ago by
Why does the final term need to be a special case? I don't like special cases unless they are justified.
comment:24 Changed 2 years ago by
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.
comment:25 Changed 2 years ago by
Do you really need the copy
in the function isometry
?
comment:26 Changed 2 years ago by
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
comment:27 Changed 2 years ago by
 Status changed from needs_review to needs_work
comment:28 Changed 2 years ago by
Can you move the helper functions out of the diagonal_isometry
function and also doctest them?
comment:29 Changed 2 years ago by
 Commit changed from d2b3dde11edb22ba199d4666b99ad753865c5322 to 4b22c5155459c95fe2fac874c49d8bc7843f08c5
Branch pushed to git repo; I updated commit sha1. New commits:
4b22c51  Added doctesting to helper methods of isometry

comment:30 Changed 2 years ago by
 Commit changed from 4b22c5155459c95fe2fac874c49d8bc7843f08c5 to 582ae4500eb3652446c6f829310298953f8dbd02
Branch pushed to git repo; I updated commit sha1. New commits:
582ae45  Minor bugfixes/refactoring on isometry

comment:31 in reply to: ↑ 23 ; followup: ↓ 34 Changed 2 years ago by
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 1dimensional quadratic forms, and the final iteration of the algorithm operates on two 1dimensional 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 nonsingular 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 functionisometry
?
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 3dimensional 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

comment:32 Changed 2 years ago by
 Status changed from needs_work to needs_review
comment:33 Changed 2 years ago by
comment:34 in reply to: ↑ 31 Changed 2 years ago by
Replying to tgaona:
This block finds a pair of vectors
v
andw
such thatQ(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 withv = [1, 0, 0]
(for a 3dimensional form) and findsw
by callingF.solve(Q(v))
. If this pair doesn't work, it finds a newv
and starts over. I get newv
'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.
comment:35 Changed 2 years ago by
 Status changed from needs_review to needs_info
comment:36 Changed 23 months ago by
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.
comment:37 Changed 23 months ago by
 Commit changed from 582ae4500eb3652446c6f829310298953f8dbd02 to 17d12c6faa5f18a9f822e30239bf201e2fe76c5f
comment:38 Changed 22 months ago by
 Status changed from needs_info to needs_review
comment:39 Changed 22 months ago by
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.
comment:40 Changed 22 months ago by
 Status changed from needs_review to needs_work
Merge conflicts with Sage 7.0
comment:41 Changed 22 months ago by
 Commit changed from 17d12c6faa5f18a9f822e30239bf201e2fe76c5f to 56b5225a8caec33e4b4366a90d1686fc0a5f6768
Branch pushed to git repo; I updated commit sha1. New commits:
56b5225  Merge branch 'develop' of git://trac.sagemath.org/sage into ticket/19112

comment:42 Changed 21 months ago by
 Commit changed from 56b5225a8caec33e4b4366a90d1686fc0a5f6768 to a56cb960a859fd2cb318d00bc06eacebb8cccb91
comment:43 Changed 21 months ago by
 Commit changed from a56cb960a859fd2cb318d00bc06eacebb8cccb91 to 642df1071e10e3a771abd7d02ceff017c0a3ef30
Branch pushed to git repo; I updated commit sha1. New commits:
642df10  Fixed typo in documentation example for _diagonal_isometry

comment:44 Changed 21 months ago by
 Status changed from needs_work to needs_review
comment:45 Changed 21 months ago by
 Milestone changed from sage6.10 to sage7.1
 Status changed from needs_review to needs_work
Why this added line?
+ v = parilist[i]
comment:46 Changed 21 months ago by
I am still missing some highlevel explanation of the algorithm. For example, what is FM
?
comment:47 followup: ↓ 48 Changed 21 months ago by
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:
 In your documentation for
_diagonal_isometry()
you call the functionisometry()
(of course it works out to be the same thing here, since the forms are diagonal to begin with.
 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.
comment:48 in reply to: ↑ 47 Changed 21 months ago by
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.
comment:49 Changed 21 months ago by
 Commit changed from 642df1071e10e3a771abd7d02ceff017c0a3ef30 to 70ed4293b680ec8859683a8e1c3b7791c1c67e73
comment:50 Changed 21 months ago by
 Status changed from needs_work to needs_review
comment:51 Changed 15 months ago by
 Commit changed from 70ed4293b680ec8859683a8e1c3b7791c1c67e73 to 0b8841492ca2e5714879cadacb6076f595927bd7
comment:52 Changed 15 months ago by
 Commit changed from 0b8841492ca2e5714879cadacb6076f595927bd7 to 2b394f41bd18eb22c6558a1441444692b79caec7
Branch pushed to git repo; I updated commit sha1. New commits:
2b394f4  Removed trailing whitespace.

comment:53 Changed 5 weeks ago by
 Reviewers set to Simon Brandhorst
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 giveis_rationally_isometric
a new argumentreturn_matrix
. You can model analogous behavior as inis_globally_equivalent_to
.
 In the current implementation you do not check the
base_ring
isZZ
orQQ
. This should be done. Or does it work more generally? For example the code does not work over the padics.  You do not need a new method
_compute_gram_matrix_from_basis
. Instead you can just useQ(T).Gram_matrix_rational()
. HereQ
is a quadratic form andT
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
New commits:
Implements a method to compute isometries between positive definite quadratic forms.
Adds a flag parameter to short_vector_list_up_to_length to pass to PARI's