Given rational coefficients, lu returns rational arrays for dense inputs, floating point arrays for rational inputs

[henselman-petrusek]

Given a rational input, the lu function seems to return rational arrays for dense inputs and floating point arrays for sparse inputs. Is there any way around this inconsistency?

julia> using LinearAlgebra
julia> using SparseArrays

julia> A = Rational.([3 1; 1 3])
2×2 Array{Rational{Int64},2}:
 3//1  1//1
 1//1  3//1

julia> F = lu(A)
LU{Rational{Int64},Array{Rational{Int64},2}}
L factor:
2×2 Array{Rational{Int64},2}:
 1//1  0//1
 1//3  1//1
U factor:
2×2 Array{Rational{Int64},2}:
 3//1  1//1
 0//1  8//3

julia> G = lu(sparse(A))
SuiteSparse.UMFPACK.UmfpackLU{Float64,Int64}
L factor:
2×2 SparseMatrixCSC{Float64,Int64} with 3 stored entries:
  [1, 1]  =  1.0
  [2, 1]  =  0.333333
  [2, 2]  =  1.0
U factor:
2×2 SparseMatrixCSC{Float64,Int64} with 3 stored entries:
  [1, 1]  =  0.75
  [1, 2]  =  0.25
  [2, 2]  =  0.666667

[stevengj]

We don't have a Julia-native sparse LU function.

Honestly, though, I'm not seeing why you would use the sparse-matrix routines here. They are only beneficial for large arrays, where Rational calculations are probably not useful — you will quickly overflow typemax(Int64).

[henselman-petrusek]

Thank you!

There's a branch of applied mathematics (topological data analysis) that deals with very large, very sparse matrices, for which this is relevant problem is relevant. One big hurdle for the field as a whole is a lack of base-level linear algebra functions available to work over the rationals. To give some sense of scale, the arrays typically have hundreds of millions of rows and columns, with ~3 nonempty entries per row and column. In numerical experiments overflow has never been an issue, even for the largest arrays. As noted in a related issue, there are practical reasons for using standard Julia sparse arrays for this type of work, even though there exist packages like Hecke that offer sparse matrix algorithms over general fields. We have tried running standard LU on our rational arrays, but conversion to float produces numerical errors above tolerance.

Overflowing typemax(Int64) does sound like a reasonable concern in general, but shouldn't that throw an error when it occurs? Do you think there's any circumstance under which a Julia-native LU algorithm could be contributed to LinearAlgebra? Thanks again!

[stevengj]

Overflowing typemax(Int64) does sound like a reasonable concern in general, but shouldn't that throw an error when it occurs?

Only if you use a checked integer type like https://github.com/JeffreySarnoff/SaferIntegers.jl, which slows things down a bit but is probably worth it here.

Do you think there's any circumstance under which a Julia-native LU algorithm could be contributed to LinearAlgebra?

They would be very useful for many purposes besides rational arithmetic. e.g. arithmetic with dual numbers, intervals, quaternions, etcetera. But it's a fairly big undertaking to re-implement the sparse-matrix libraries like UMFPACK from scratch. (It's the sort of thing that is most likely to happen in a package first, not in Base.)

[henselman-petrusek]

This is very helpful to know, thank you. We've gone back and forth about whether to implement our own algorithms, but if Base doesn't have this functionality on the horizon that clarifies the to do list considerably. We certainly won't be re-implementing UMFPACK in its entirety, more like LU and a select few other factorizations, so it'll be up to some other good samaritans to help build such a package. I agree, it would be very useful.

Shall I close the issue or leave open to anyone else who wants to ask/contribute? Thanks again!

[stevengj]

I would close it. If you want to open a more general issue to track pure-Julia sparse direct solvers, that would be better — anything we do here would not be specific to rationals, probably.

[jmichel7]

Note that rational arithmetic is checked for overflow in Julia.