Hecke.jl
Hecke.jl copied to clipboard
thofma
Precision loss in p-adic linear algebra
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.
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)
I guess in my case I want to call a division free algorithm ... but AA sees the padics as a field.
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).
julia> AbstractAlgebra._solve_ff(A,b)
ERROR: System not solvable in _solve_ff
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
@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.
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...