While computing Smith normal form, I got bitten

[Saul Schleimer]

I sent the following to John Cremona - he suggested I post it here.

Dear John -

  I was computing Smith normal forms (over ZZ) and I ran into an issue.  I boiled it down to this:

sage: version()
'Sage Version 5.12, Release Date: 2013-10-07'
sage: M = matrix([[1,0],[0,1]])
sage: M
[1 0]
[0 1]
sage: M.parent()
Full MatrixSpace of 2 by 2 dense matrices over Integer Ring
sage: M.inverse()
[1 0]
[0 1]
sage: M.inverse().parent()
Full MatrixSpace of 2 by 2 dense matrices over Rational Field

Note that the base ring has changed.  This caught me out: I was computing Smith normal forms in two ways and they were different.  Eventually I realized that (obviously) the Smith normal form is sensitive to the base ring!

Does the base ring have to change here?  I guess that there is a subtle algorithm to quickly compute inverses, that requires extending the ring.  But once we've got the answer, couldn't we do a fast check to see if the answer is actually in the original base ring?

That is - if M is the original matrix and N is the computed inverse then add this to the very end:

if N.parent() == M.parent():
    return N

try:

    R = M.parent().base_ring()

    P = N.change_ring(R)

    assert M*P == P*M == 1

    return P

except:

    return N

all the best,

saul

[Edit]
John pointed me to the code for M.__invert__() so I suppose that I am suggesting placing some version of the above code after the try and before the return, around line 11.

"try:
                        return (~self.lift()).change_ring(self.base_ring())"

Of course, in this particular place in the code the first two lines of my suggestion are pointless. 

best,

saul

[David Loeffler]

Hi Saul,

This is a consistent design decision across Sage: when you invert a
ring element, it automatically returns the inverse as an element of
the fraction field. E.g.

sage: a = ZZ(1)
sage: a.parent()
Integer Ring
sage: b = ~a
sage: b.parent()
Rational Field

Same with polynomials, power series, etc. Generally the idea is that
the parent of the result of an operation should depend only on the
parents of the operands (not on their actual values). I don't know
exactly what the arguments for and against this are, but if it's to be
changed we would certainly want to change it consistently across all
Sage parents (not just for matrices alone).

David

> --
> You received this message because you are subscribed to the Google Groups
> "sage-support" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sage-support...@googlegroups.com.
> To post to this group, send email to sage-s...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sage-support.
> For more options, visit https://groups.google.com/groups/opt_out.

[Saul Schleimer]

Hi David -

This is a consistent design decision across Sage: when you invert a
ring element, it automatically returns the inverse as an element of
the fraction field.

It seems to me that this violates the principle of "least surprise":  if I have
a unit in a ring, and I invert it, I can reasonably expect answer to be a unit,
in that ring...

Here is (a reconstruction of) what happened to me:  I was
preparing a question on Smith normal forms for my class
on homology.  I computed

sage: A = matrix(ZZ, [[-6, -26, -82], [0, 4, 3], [1, 0, 7], [0, 2, 5], [2, 10, 30]])
sage: D, U, V = A.smith_form()
sage: U*A*V == D
True
sage: D == (U.inverse()*D).smith_form()[0]
False

I was thrown by this.  The theorem for Smith normal form tells me that U lies
in GL(5,ZZ) so multiplying by U.inverse() cannot change the Smith normal
form...  Indeed, I print U and U.inverse() and both are integer matrices...  

Generally the idea is that the parent of the result of an operation should 

depend only on the parents of the operands (not on their actual values). 

I don't know exactly what the arguments for and against this are,

Ah, ok.  I'm interested to know the for and the against here.

best,

saul

[Simon King]

Hi Saul,

On 2014-01-23, Saul Schleimer <sau...@gmail.com> wrote:
> It seems to me that this violates the principle of "least surprise": if I
> have a unit in a ring, and I invert it, I can reasonably expect answer to
> be a unit, in that ring...

Yes and no.

If you have an integral domain, you can reasonably expect the answer to be an
element of the fraction field: 1/1 naturally belongs to QQ, and it would be
very confusing if the type of 1/n would depend on n: Think of writing a program
that needs special cases all over the place.

But in your example, it is a matrix algebra, hence, no integral domain.
In this case, it might indeed be reasonable to expect the inverse of a
unit to stay in the same ring.

I think the discussion whether or not the algebraic structure (aka
"Parent",which means "object in a subcategory of the category of sets)
containing the result of arithmetic operations should depend
on the operands or *only* on the parents containing the operands took
place before I joined Sage. However, I can give some reasons for letting
the resulting parent only depend on the operand's parents:

- Sage's coercion system ("automatic conversion into a suitable parent
when doing arithmetic") relies on morphisms, which I think is the right
thing to do for providing a mathematical fundament to coercion (which goes
way beyond the "automatic type conversion" known from the C
programming language). Hence, an arithmetic operation such as *
involving parents R and S is modelled as a morphism (in an appropriate
category) RxS --> T, for another parent T (multiplication of elements
of R and S happens after mapping both elements to T). Of course, this
will not always be possible, there may only be a partial map or there
may be no canonical choice of a morphism. In that case, Sage would
normally refuse to do the multiplication (but I guess the user can
force it by doing explicit rather than implicit conversions).

Long paragraph in a nutshell: Sage's arithmetic system largely relies
on morphisms, and hence it would be against some of Sage's fundaments
to make the result of an arithmetic operation live in different parents,
depending on the operands.

- Imagine you have a huge number of arithemtic operations to do, with
operands that all belong to a fixed pair of parents. Now think what
happens if for *each* pair of operands you have to do (expensive)
tests before being able to choose the appropriate algorithm: It's a
waste of time! Instead, you do the expensive test *once* (namely only
for the fixed pair of parents), resulting in a single action, and
then you feed this single action with a zillion of operand pairs.
That's a lot more efficient than creation zillions of actions,
namely one for each operand pair!

Example for the second point:

sage: R = ZZ[x]
sage: S = QQ
sage: A = R.get_action(S, operator.mul)
sage: A
Right scalar multiplication by Rational Field on Univariate Polynomial
Ring in x over Integer Ring
sage: A.codomain()
Univariate Polynomial Ring in x over Rational Field

Hence, you have one action, taking a rational number and a polynomial
over the integers, and returning a polynomial over the rationals. And
this single action is in fact what Sage uses internally, for all
elements of R and S. So, when you do
sage: r = R.random_element(); r
x^2 - 8*x + 2
sage: s = S.random_element(); s
-1/2
sage: r*s
-1/2*x^2 + 4*x - 1

then internally Sage just looks up A (it is cached) and does
sage: A(r,s)
-1/2*x^2 + 4*x - 1

Best regards,
Simon

[Saul Schleimer]

Hi Simon -

  Nice to hear from you.  This discussion of morphisms is pretty convincing.  I tried it out, but I think I am doing something wrong.

sage: Id = matrix([[1,0],[0,1]])
sage: type(Id)
<type 'sage.matrix.matrix_integer_dense.Matrix_integer_dense'>
sage: A = Mat.get_action(Mat, operator.mul); A
Left action by Full MatrixSpace of 2 by 2 dense matrices over Integer Ring on Full MatrixSpace of 2 by 2 dense matrices over Integer Ring
sage: A.codomain()

Full MatrixSpace of 2 by 2 dense matrices over Integer Ring

sage: A = Mat.get_action(Mat, operator.inv); A
sage: A.codomain()
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-414-ef5880878bbf> in <module>()
----> 1 A.codomain()

AttributeError: 'NoneType' object has no attribute 'codomain'

What am I doing wrong?

best,

saul

[Saul Schleimer]

Dear all -

Here is a much lighter way make my "problem" go away. 
Suppose that A is an m by n matrix (number of rows by number of columns). 

A = matrix(ZZ, [[-6, -26, -82], [0, 4, 3], [1, 0, 7], [0, 2, 5], [2, 10, 30]])

sage: type(A)
<type 'sage.matrix.matrix_integer_dense.Matrix_integer_dense'>

Currently we get this:

sage: D, U, V = A.smith_form()

sage: [type(x) for x in A.smith_form()]
[sage.matrix.matrix_integer_dense.Matrix_integer_dense,
 sage.matrix.matrix_integer_dense.Matrix_integer_dense,
 sage.matrix.matrix_integer_dense.Matrix_integer_dense]

So U (say) is a matrix.  Now, the whole point of U's existence is to be a product of
elementary matrices, that do row operations to A, trying to turn it into D.  So if we
change smith_form to return U as an element of GL(m,ZZ) (and return V as an
element of GL(n.ZZ)) then I am good to go.  Just to be clear what I am talking about:

sage: G = GL(5,ZZ)
sage: P = G.generators()[0]
sage: type(P)
<class 'sage.groups.matrix_gps.group_element.LinearMatrixGroup_gap_with_category.element_class'>
sage: type(P.inverse())
<class 'sage.groups.matrix_gps.group_element.LinearMatrixGroup_gap_with_category.element_class'>

This is good...

sage: type(A)
<type 'sage.matrix.matrix_integer_dense.Matrix_integer_dense'>
sage: type(P*A)
<type 'sage.matrix.matrix_integer_dense.Matrix_integer_dense'>
sage: A.smith_form()[0] == (P*A).smith_form()[0]
True

and I am happy.  Of course, in general U and V lie in GL(m, R) and GL(n,R)
where R is the parent of the entries of A.

best,

saul

[Simon King]

Hi Saul,

On 2014-01-24, Saul Schleimer <sau...@gmail.com> wrote:
> Nice to hear from you. This discussion of morphisms is pretty
> convincing. I tried it out, but I think I am doing something wrong.
>

> ...

>
> sage: A = Mat.get_action(Mat, operator.inv); A
> sage: A.codomain()
> ---------------------------------------------------------------------------
> AttributeError Traceback (most recent call last)
><ipython-input-414-ef5880878bbf> in <module>()
> ----> 1 A.codomain()
>
> AttributeError: 'NoneType' object has no attribute 'codomain'
>
> What am I doing wrong?

Nothing. What I said was an attempt to explain why Sage generally uses a
"per parent" and not "per element" approach towards arithmetic/algebra:
It makes the coercion system nicer, and for consistency one might prefer
to have a "per parent" approach even when coercion is actually not
directly involved.

Here, we don't have an action (with two parents/two operands), since
operator.inv is a unary operator. Hence, A in your example is None, and
no coercion of the *matrices* is involved. However, in the underlying
inversion algorithm, I suppose one sometimes has to divide by the
*coefficients* of the matrix, and thus inversion of an integer matrix
naturally leads to a matrix whose coefficients live in the fraction
field of ZZ, even if in the end the coefficients have denominator one.

Inversion is implemented (as usual in Python) by the `__invert__`
method. In the case of matrices, one reads in the documentation:
"""
Return this inverse of this matrix, as a matrix over the fraction
field.
...
"""

In principle it would be legal to have something that is not a morphism
(hence, choose the result's parent based an a "per element" decision),
since it is not part of the coercion system.

I guess, after all it was a conscious design decision that may not be
*forced* by the coercion system, but certainly *suggested* by the coercion
system.

Best regards,
Simon