Matrix inverse over rational function field is slow

[awerick11]

I came across a matrix (over rational_function_field(QQ, :q)) where taking its inverse makes Julia hang. While it hangs, memory usage increases seemingly unbounded. After an interruption, futher matrix computations make Julia abort and throw a malloc error, (malloc_consolidate(): unaligned fastbin chunk detected).

I know that the matrix is invertible, as it squares to the identity.

The matrix is below.

[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 (-q^8+q^6+2*q^4+q^2+1)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (2*q^7)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (-2*q^6)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (-2*q^6)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (2*q^5)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (-2*q^4)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 0; 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 (2*q^7)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (q^8-q^6+2*q^4+q^2+1)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (2*q^5)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (2*q^5)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (-2*q^4)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (2*q^3)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 (-2*q^6)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (2*q^5)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (-2*q^4)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (q^8+q^6+q^2+1)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (2*q^3)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (-2*q^2)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 (-2*q^6)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (2*q^5)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (q^8+q^6+q^2+1)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (-2*q^4)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (2*q^3)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (-2*q^2)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 (2*q^5)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (-2*q^4)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (2*q^3)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (2*q^3)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (q^8+q^6+2*q^4-q^2+1)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (2*q)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0; 0 0 0 0 0 (-2*q^4)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (2*q^3)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (-2*q^2)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (-2*q^2)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (2*q)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 (q^8+q^6+2*q^4+q^2-1)//(q^8+q^6+2*q^4+q^2+1) 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]

[joschmitt]

First of all, I should point out the difference between Julia matrices and AbstractAlgebra matrices. You can read the details here: https://docs.oscar-system.org/stable/General/faq/#Frequently-Asked-Questions (scroll to the question "Why do you have your own matrix types, and why do they not support the exact same commands as Julia matrices?"). The point is that the matrix should be constructed via:

julia> F, q = rational_function_field(QQ, :q);

julia> M = matrix(F, [1 0 0 0 0 ...])

(Maybe you did this, I am not sure.)

I don't know what you mean by "Julia hangs". The computation did not terminate quickly for me either and I aborted it after around 3 minutes (and 12 GB RAM). But I don't think this is a bug; computing over this kind of field is not easy.

However, I would recommend to use Nemo.jl. This builds on top of AbstractAlgebra, but adds much more efficient implementations, for example, for polynomials over the rationals. With Nemo, computing the inverse of this matrix takes less than a second.

[awerick11]

I did construct it as an AbstractAlgebra matrix. The matrix I put in the post is simply the output of print().

I believe I let it run for at least an hour without termination. I suspected that this is a bug as I've been able to do much larger computations in shorter times (e.g. multiplying 1216 x 1216 matrices over the same field), though of course the particular matrices are very critical to the results. I don't actually have any hard evidence that there's a bug, though; it's believable that the performance is sensitive by orders of magnitude to the input matrix.

Thank you for the tip to use Nemo! I think I misread its docs saying that there's specialized functionality for series over Q as being only about Q.

[thofma]

I agree with @joschmitt, that this is not a bug, but it is just slow. But I think the strategy of the algorithm is just wrong. We should probably always go via the polynomial ring:

julia> d = mapreduce(denominator, lcm, A) # A is your matrix
q^8 + q^6 + 2*q^4 + q^2 + 1

julia> @time map_entries(numerator, d*A) |> pseudo_inv;
38.115909 seconds (846.10 M allocations: 28.331 GiB, 26.24% gc time, 1.98% compilation time)