Precision loss in p-adic linear algebra

[simonbrandhorst]

Maybe the following example is not a bug ... but the solution computed is suboptimal
The following computations can be done modulo 29^2.
They should be possible with padics?

julia> k = PadicField(29,2)
Field of 29-adic numbers

julia> A = k[0 29 1; 29 1 0; 1 0 0]
[       O(29^2)   29^1 + O(29^3)   29^0 + O(29^2)]
[29^1 + O(29^3)   29^0 + O(29^2)          O(29^2)]
[29^0 + O(29^2)          O(29^2)          O(29^2)]

julia> A = k[0 29 1; 29 1 0; 1 0 0];

julia> A
[       O(29^2)   29^1 + O(29^3)   29^0 + O(29^2)]
[29^1 + O(29^3)   29^0 + O(29^2)          O(29^2)]
[29^0 + O(29^2)          O(29^2)          O(29^2)]

julia> det(A) # precision loss
28*29^0 + O(29^1)

julia> inv(A)   # precision loss
[       O(29^0)          O(29^0)             O(29^0)]
[       O(29^2)   29^0 + O(29^2)   28*29^1 + O(29^2)]
[29^0 + O(29^1)          O(29^1)             O(29^1)]

julia> pseudo_inv(A)[1]  # not as bad but still not good
[          O(29^2)             O(29^2)          O(29^2)]
[          O(29^2)   28*29^0 + O(29^1)   29^1 + O(29^2)]
[28*29^0 + O(29^1)             O(29^2)          O(29^2)]

julia> b = k[1; 1; 1;];

julia> x1 = solve(A,b;side=:right)  #precision loss
[                 O(29^0)]
[29^0 + 28*29^1 + O(29^2)]
[          29^0 + O(29^1)]

julia> A*x1
[       O(29^0)]
[29^0 + O(29^1)]
[       O(29^0)]

julia> A*x1 == b
true

julia> b
[29^0 + O(29^2)]
[29^0 + O(29^2)]
[29^0 + O(29^2)]

julia> x2 = inv(A^2)*A * b  # this is the best possible and desired solution
[          29^0 + O(29^2)]
[29^0 + 28*29^1 + O(29^2)]
[29^0 + 28*29^1 + O(29^2)]

julia> A*x2
[29^0 + O(29^2)]
[29^0 + O(29^2)]
[29^0 + O(29^2)]

julia> A*x2 == b
true

Similar, but with wierd printing

julia> k = PadicField(29,2)
Field of 29-adic numbers

julia> K,toK = completion(F,prime_ideals_over(maximal_order(F),29)[1], 2)
(Eisenstein extension with defining polynomial x + 26*29^1 + O(29^2) over Unramified extension of 29-adic numbers of degree 1, Map: F -> eisenstein extension with defining polynomial x + 26*29^1 + O(29^2) over QQ_29^1)

julia> A = K[0 29 1; 29 1 0; 1 0 0];

julia> b = K[1; 1; 1;];

julia> A*solve(A,b;side=:right)  # not really a solution
[1]
[1]
[0]

In some larger examples solve complains that no solution exists, although it clearly does.

[simonbrandhorst]

We can go a bit further and create an error

julia> k = PadicField(29,2)
Field of 29-adic numbers

julia> A = k[29^2 29 1; 29 1 0; 1 0 0];

julia> A[1,1] = setprecision!(A[1,1],3)
29^2 + O(29^3)

julia> b = k[1; 1; 1;];

julia> A
[29^2 + O(29^3)   29^1 + O(29^3)   29^0 + O(29^2)]
[29^1 + O(29^3)   29^0 + O(29^2)          O(29^2)]
[29^0 + O(29^2)          O(29^2)          O(29^2)]

julia> inv(A)
ERROR: Singular matrix in inv
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] inv(M::AbstractAlgebra.Generic.MatSpaceElem{PadicFieldElem})
   @ AbstractAlgebra ~/.julia/packages/AbstractAlgebra/5q3jO/src/Matrix.jl:3439
 [3] top-level scope
   @ REPL[303]:1

julia> det(A)
O(29^2)

[simonbrandhorst]

I guess in my case I want to call a division free algorithm ... but AA sees the padics as a field.

[thofma]

Yes, unfortunately we don't have any useful linear algebra over inexact fields. Also, the whole model for p-adics, unramified extenions and arbitrary local fields is currently remodelled.

In the meantime, you might make use of AbstractAlgebra._solve_ff and AbstractAlgebra.det_df. They are supposed to be division free. (I have used det_df successfully in the past for local fields, don't know how useful _solve_ff will be).

[simonbrandhorst]

julia> AbstractAlgebra._solve_ff(A,b)
ERROR: System not solvable in _solve_ff

[fieker]

In this case, where the elementary divisors are all 1, the ring-approach together with the "Book" fixes will do the job neatly:

A = map_entries(maximal_order(base_ring(A)), A)
s, u, v = snf_with_transform(A)

julia> v*u*b
[          29^0 + O(29^2)]
[29^0 + 28*29^1 + O(29^2)]
[29^0 + 28*29^1 + O(29^2)]

Da haette noch duth "s geteilt" werden muessen, da kommt das der Praezisionsverlusst rein

[simonbrandhorst]

@fieker interesting, I wasn't aware of maximal_order to give the ring of p-adic integers.
However after conversion, I still get precision loss in my examples.
I think I will just give up and use the residue rings of the number field.

[fieker]

I have a patched version of hnf/ snf that will produce transformation matrices without precisio loss, precisely for this purpose. It's still under discussion...