Rational Linear Algebra in Julia

[David Zhang]

Does Julia have any support for linear algebra on rational-valued matrices? Currently, most linear algebra functionality seems to be provided by BLAS/LAPACK, which only support operations on floating-point matrices.

I am interested in this because I would like to use Julia to derive high-order numerical integration methods, whose coefficients are the solutions of rational linear system.

(Continuation from http://stackoverflow.com/questions/20985783/rational-matrix-division-in-julia)

[Tim Holy]

If you look in base/linalg/matmul.jl, you'll find functions like
generic_matmatmul that do multiplication on arbitrary matrices.
Division/factorization is a completely different story, I think. You could
presumably implement a gaussian-elimination algorithm (or something fancier)
in pure Julia, and that would give you the ability to work with Rational
matrices.

--Tim

[Viral Shah]

There is an open issue on this, and it shouldn't be too difficult to put something in place.

https://github.com/JuliaLang/julia/issues/3014

-viral

[Alasdair McAndrew]

Since Julia plays so nicely with Python, why not just - at least in the short term - create a package that interfaces with SymPy (http://sympy.org/en/index.html)?  This will certainly give you what you want, and a great deal more besides.

[Tim Holy]

On Wednesday, January 08, 2014 06:07:43 PM Viral Shah wrote:
> There is an open issue on this, and it shouldn't be too difficult to put
> something in place.

But to be realistic, it requires a motivated volunteer to do it, and it seems
unlikely that this would just happen to be sitting at the very top of anyone
else's priority list right now. David?

--Tim

[j verzani]

There is SymPy.jl. The matrix bit could use some ironing out. By default, it uses an Array of symbolic objects, rather than a symbolic array. This means taking an inverse is bit cumbersome:

using SymPy

r = sympy.Rational 

m = [r(1,2) r(2,3); r(3,4) r(4,5)] ## Array of Sym objects

inv(m) ## converts to float, so no good

m = convert(SymMatrix, m)     ## SymMatrix object

m[:inv]()

Might be better to just interface directly with sympy through PyCall for this.