The reason is that rank begins by converting its argument to floating-point.
It seems to me that given a field (Rationals, but it could be rational fractions or any other field),
the rank should be computed by computing the echelon form of the matrix rather that trying to
convert to floats (which would not make sense for rational fractions, for instance).
Usually, rank is a thin wrapper over some decomposition, typically SVD, and just checks the singular values. Regardless of the algorithm, rationals can very quickly blow up even for small matrices.
The point of Rationals, for instance Rational(BigInt) is to make exact computations (for instance in case of blowup). But mathematicians in Group theory (my field) or number theory are interested in exact computations more generally.
I am not an expert, but I guess there would be no harm in implementing some LinearAlgebra functions for Rational{T}. However, you could still easily get overflows except for T === BigInt, so it is unclear whether one should just convert to that immediately for all T, otherwise it’s not so useful.
To get back to my example of Hilbert matrices, they can be inverted without trouble using Rational{Int128} up to size 25. You get an overflow (an explicit error) trying to invert size 26. However, the rank is wrong starting from size 11. The rank is also wrong using BigInt. This is an unfortunate lack of consistency. My claim is that the algorithm to compute the rank is wrong for all fields excepted “approximate Real numbers”. A proper algorithm would work for any field, like Rational Fractions, Algebraic Numbers, etc…
The code there has a suspicious ε which may not make sense in every field. I was thinking of something like this:
" returns: echelon form of m, indices of linearly independent rows of m"
function echelon!(m::Matrix)
rk=0
inds=collect(axes(m,1))
for k in axes(m,2)
j=findfirst(x->!iszero(x),m[rk+1:end,k])
if j===nothing continue end
j+=rk
rk+=1
row=m[j,:]
m[j,:].=m[rk,:]
m[rk,:].=inv(row[k]).*row
inds[[j,rk]]=inds[[rk,j]]
for j in axes(m,1)
if rk!=j && !iszero(m[j,k]) m[j,:].-=m[j,k].*m[rk,:] end
end
end
m,inds[1:rk]
end
echelon(m::Matrix)=echelon!(copy(m))
I am sure there are many optimizations of the above code possible since I am novice to Julia. For instance, it could be more efficient to compute a column echelon form instead of a row echelon form.
The standard LinearAlgebra package is indeed focused on floating-point, so that one might better say that the type promotions there are not adequately documented, rather than that “rank is wrong.”
For serious treatment of rationals (and other abstract-algebraic or number-theoretic matters), it is appropriate to use other packages such as Nemo.jl where you may be happier to see that rank(hilbert(MatrixSpace(QQ,n,n))) == n.
I do not buy these arguments. I think Julia, even though it was developed for the purpose of numerical analysis, is able to be a good programming language in general. I am sure the Julia developers are not deliberately hostile, as you seem to be, to people using Julia for broader mathematics than numerical analysis.
As a mathematician, I do not think I have to take refuge in a ghetto as Nemo is currently. And on the point
currently debated, I am sure that even hard-core numerical analysts would like to see an exact answer if they ask for the rank of a Matrix{Rational{BigInt}}.
Almost nobody in numerical world does computations with Rationals, and this is why “rational path-code” is not as complete and as well maintained. Mainly because
julia> A = rand(1:100, 100); B = rand(1:100, 100); C = A.//B; sum(C)
ERROR: OverflowError: 140428563049948639 * 83 overflowed for type Int64
...
you are not able to compute with machine Integers a sum of 100 rationals with small num/denum. As for your example: it’s not that developers don’t care, maybe nobody actually needed this. And if you care about rank, please submit a PR.
Is this the best option though? For floats there are a variety of rank-revealing factorisations (typically involving some sort of pivoting): typically you can modify these factorizations to work with rationals, which I would have thought would be more efficient. Also I imagine you would want to minimise the risk of overflow.
It appears like the same problem applies to integer matrices:
julia> m=[1//(n+m) for n in 1:11, m in 1:11];
julia> p = lcm([n+m for n=1:11 for m=1:11])
232792560
julia> mi = Int.(p*m); Mi = BigInt.(p*m);
julia> rank(mi)
10
julia> rank(Mi)
ERROR: MethodError: no method matching svdvals!(::Array{BigFloat,2})
Re echeleon form: There are some pitfalls involving large intermediate values, and very naive calculation requires (inexact / rational) division or is in a suboptimal complexity class in the number of input bits. I vaguely remember that sufficiently nasty examples with sufficiently naive variants of Gauss can even get exponential runtime over Rational{BigInt}. Try this?
Presumably nemo is capable of handling rank over PIDs like integers and fields like rationals. Maybe nemo has a good generic implementation with compatible license or friendly authors that can just be included in base?
As people have already pointed out, the current rank function is for computing the floating point rank and it can be hard to combine convenient promotion while correctly limiting the input, e.g. we can’t just require the elements to be <:AbstractFloat since that would exclude complex numbers and quaternions based on floating point numbers.
While our current support for rationals is far from complete, it would be fine to make rank compute the correct rank for rational matrices. However, I don’t think reducing to echelon is an efficient way of doing it. I suspect that we could simply use the LU and check for zeros in the diagonal of U. However, we’d still have to figure out how to dispatch to this version. I’m not sure how to do that yet. Worst case would be to rename rank to numericalrank and use rank only for exact computations but I wouldn’t like that.
This is an important problem for me. To make the julia library friendly to other mathematics than floating-point numbers, it is needed to be able to dispatch on “approximate numbers” (Floats, Complex of Floats, Quaternions of Floats) versus “something else” (Rationals, user-defined types like algebraic numbers, rational fractions, whatever). The grey zone will be Integers, which are often treated as Floats, but not always, by Julia (I think BigInts should never be treated as floats). Has there been any thoughts along these lines among developers/designers of the language?
Do you mean an issue or a PR? Anyway, I have no precise proposal yet, and I think the current discussion
is OK to get the opinion of various people on the topic, including developers/committers who drop by.
