Inverse of a p-adic matrix fails. Is this a bug?

[Simon Brandhorst]

The following may be a bug or me not understanding p-adic floating point computations:

sage: R = Qp(2,type='floating-point',print_mode='terse')
sage: M = Matrix(R,4,[0, 0, 1, 1, 2^20, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1])
sage: M.det()
1048575
sage: M.inverse()
---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
....

ZeroDivisionError: input matrix must be nonsingular

sage: M = Matrix(R,4,[0, 0, 1, 1, 2^19, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1])
sage: M.inverse()
[     1048575            2            0            1]
[274878955520       524288            1      1048575]
[      524289       524287            1       524287]
[274877382656       524289      1048575       524289]

Works.

[Simon King]

Hi!

I am afraid I don't know a definitive answer to your question.
But at least I find the following a bit strange:
sage: d = M.det()
sage: d == -1
True
sage: d == 0
False

So, the inversion of M should work, I think. Let us do compute
the inverse "manually" (i.e., by applying Gaussian elimination
to the augmented matrix) to see where it fails:

sage: N = M.augment(M.parent()(1))
sage: N.swap_rows(0,3)
sage: N.add_multiple_of_row(1,0,-N[1,0])
sage: N.add_multiple_of_row(2,0,-N[2,0])
sage: N.add_multiple_of_row(2,1,-N[2,1])
sage: N.rescale_row(2,~N[2,2])
sage: N.add_multiple_of_row(3,2,-N[3,2])
sage: N.add_multiple_of_row(1,2,-N[1,2])
sage: N.add_multiple_of_row(0,2,-N[0,2])
sage: N.rescale_row(3,~N[3,3])
sage: N.add_multiple_of_row(2,3,-N[2,3])
sage: N.add_multiple_of_row(0,3,-N[0,3])
sage: N[:,:4] == 1
True

Therefore, N[:,4:] should be inverse to M. But it isn't:

sage: N[:,4:]*M == 1
False
sage: M*N[:,4:] == 1
False

I suppose that there is some p-adic rounding going on. Which would also
explain why there is a division-by-zero error (a number which is
p-adically close to zero may have been rounded to zero).

Best regards,
Simon

[Nils Bruin]

I found this over CC and RR as well: the generic matrix inversion in sage (which is probably what is getting invoked) looks like it's based on straight-up row reduction: No pivoting (which it wouldn't know how to do in the generic case). For CC and RR there's PLU decomposition in mpmath which does use partial pivoting. It would be nice to have "natively" in sage and would be easy to program and would be nice to have for non-square matrices. The code is easily adapted (see https://git.sagemath.org/sage.git/commit?id=f2ccf44acc673458aba3cf92d7692c7a92a505b6 for a POC).

Partial pivoting for p-adic PLU decomposition would be straightforward to implement as well, and very worthwhile to have. It's well-known that matrix inversion is tricky business numerically, even for well-conditioned matrices, but with pivoting one should get reasonable results in a lot more cases. p-adics might be a little better behaved, but fundamentally the same numerical considerations play a role.

[Nils Bruin]

On Tuesday, January 30, 2018 at 3:50:45 PM UTC, Simon Brandhorst wrote:

The following may be a bug or me not understanding p-adic floating point computations:

sage: R = Qp(2,type='floating-point',print_mode='terse')
sage: M = Matrix(R,4,[0, 0, 1, 1, 2^20, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1])
sage: M.det()
1048575
sage: M.inverse()
---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)

So, for this particular example:

sage: P=Matrix(R,4,[0,0,0,1,0,1,0,0,0,0,1,0,1,0,0,0])
sage: M* (P*(M*P)^(-1))
[            1             0             0             0]
[1099510579200             1             0       1048576]
[            0             0             1             0]
[            0             0             0             1]

Which is indeed congruent to the identity matrix modulo 2^20 (which is about what you could hope for)

The key is that the inversion routine should be figuring out a suitable permutation on the fly.


 

[Nils Bruin]

On Tuesday, January 30, 2018 at 3:50:45 PM UTC, Simon Brandhorst wrote:

sage: M = Matrix(R,4,[0, 0, 1, 1, 2^19, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1])

sage: M.inverse()
[     1048575            2            0            1]
[274878955520       524288            1      1048575]
[      524289       524287            1       524287]
[274877382656       524289      1048575       524289]

Works.

 

Only in the sense that you don't get an error:

sage: min([a.valuation() for a in (M*M^(-1)-M^0).list()])
1

that's quite far from what it should be. With appropriate pivoting of course:

sage: min([a.valuation() for a in (M*(P*(M*P)^(-1))-M^0).list()])
20

 

[John Cremona]

Here is one way (I am not saying it is always a good strategy):

sage: R = Qp(2,type='floating-point',print_mode='terse')

sage: M = Matrix(R,4,[0, 0, 1, 1, 2^20, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1])

sage: Mi = (M.change_ring(QQ)^(-1)).change_ring(R)

sage: Mi

[1048575       0       0       1]

[      0       0       1 1048575]

[      1 1048575       1 1048575]

[1048576       1 1048575       1]

sage: M*Mi

[      1       0       0       0]

[      0       1       0       0]

[1048576       0       1       0]

[1048576       0       0       1]

sage: Mi*M

[1 0 0 0]

[0 1 0 0]

[0 0 1 0]

[0 0 0 1]

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

[Simon Brandhorst]

Thank you for your replies.

For the application I had in mind it seems that 'floating-point' numbers are not suitable.

I was doing computations in

sage: R = Zp(2,type='fixed-mod')

sage: M = Matrix(R,4,[0, 0, 1, 1, 2^19, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1])

sage: M.inverse().base_ring()
2-adic Field with floating precision 20
sage: (M.adjoint()*M.det().inverse_of_unit()).base_ring()
2-adic Ring of fixed modulus 2^20

Unfortunately the syntax is a little cumbersome. Is there a shortcut to this?
Invert a matrix without exending the base_ring (or throw an error)

[Julian Rüth]

Hi Simon,

Linear Algebra is known to have issues over the p-adics in Sage. The generic algorithms used are mostly unaware of precision issues. Xavier Caruso has been working on improving Sage's Linear Algebra capabilities for example at https://trac.sagemath.org/ticket/23505 and https://trac.sagemath.org/ticket/23450. You might want to raise this issue again on the sage-padics list (https://groups.google.com/forum/#!forum/sage-padics)

julian

[Kiran Kedlaya]

Would it make sense to implement a "generic pivot" method, which takes as an argument a function that does the pivot selection (that would then be handled separately over RR/CC vs p-adics)?

Kiran