decompositions for Rational matrices

[mlubin]

Currently, lu, chol, (\), etc. convert an Array{Rational{Int64},2} to Array{Float64,2} before doing any computation. These should actually compute the decompositions using rational arithmetic. Textbook implementations would likely be fine to start with, since the only possible numerical issues are integer overflow.

Of course it's not always possible to use rational arithmetic. For QR it's pretty reasonable to just convert to floating-point.

[StefanKarpinski]

That would be very cool.

[johnmyleswhite]

Indeed, it would be pretty awesome.

[dmbates]

Isn't there a problem with the square roots? A Cholesky decomposition begins by taking the square root of the [1,1] element. Unless I am misremembering (always a possibility as my math classes were many, many years ago) that takes you out of the rational numbers right away. You may be able to avoid the square root in the LDLt form of the decomposition.

[mlubin]

That's right. It seems that one must use the LDLt form.

[toivoh]

On thing to look out for is that the numerators and denominators would probably blow up pretty quickly.
Unless using bigints, you would have to be very careful to deal with integer overflow.

[StefanKarpinski]

That was my initial reaction too.

[mlubin]

For the use case I was thinking of, all of the elements in the matrix have small numerators and denominators, so I wouldn't expect overflow to be an issue with Int64. To some extent the user is responsible for selecting the right precision, especially since the base integer type must be explicitly stated to use rational matrices.

I can also think of some quick sanity checks that would detect many cases of overflow. For example, when eliminating a variable in Gaussian elimination, one could check something of the form a - (a/b)*b == 0//1.

[andreasnoack]

This one is almost completed. You can compute the lu in rational arithmetic and therefore solve linear systems of equations, which I think is the main goal for this issue.

The only remaining part is the ldlt factorization of symmetric or hermitian matrices. I guess these need pivoting to be stable for floating point arithmetic, i.e. what we and LAPACK call Bunch-Kaufman. The pivoting in this case is a little tedious and I wonder if we really need our own version of this factorization.

[mlubin]

Having LU factorizations is a big step forward. For positive definite rational matrices, one can actually implement LDLt without pivoting, though I wouldn't say that this is essential.

[jiahao]

Cross-ref: JuliaLang/julia#5381, JuliaLang/julia#5430, JuliaLang/julia#5526, #61

[mlubin]

We've made a lot of progress on this. Only thing missing I'd say is LDLt for symmetric positive definite matrices.

[ViralBShah]

Much of this works now. Should we still leave this open?