# Rank is wrong for Rational matrices

**URL:** <https://discourse.julialang.org/t/rank-is-wrong-for-rational-matrices/19380>\
**Category:** General Usage\
**Created:** [January 8, 2019, 10:46am UTC](https://discourse.julialang.org/t/rank-is-wrong-for-rational-matrices/19380 "2019-01-08T10:46:30Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Jean\_Michel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jean_michel/32/8282_2.png) [@Jean\_Michel](https://discourse.julialang.org/u/Jean_Michel)\
**Post date:** [January 8, 2019, 10:46am UTC](https://discourse.julialang.org/t/rank-is-wrong-for-rational-matrices/19380/1 "2019-01-08T10:46:30Z")

</div>

consider the matrix:

`julia> m=[1//(n+m) for n in 1:11, m in 1:11];`

it can be inverted exactly with no problem:

```julia
julia> one(m)==inv(m)*m
true

```

however:

```julia
julia> rank(m)
10

```

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).

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 8, 2019, 12:19pm UTC](https://discourse.julialang.org/t/rank-is-wrong-for-rational-matrices/19380/2 "2019-01-08T12:19:05Z")

</div>

Are you using this in a pedagogical setting?

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.

---

<div class="post-metadata">

**Author:** ![Jean\_Michel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jean_michel/32/8282_2.png) [@Jean\_Michel](https://discourse.julialang.org/u/Jean_Michel)\
**Post date:** [January 8, 2019, 1:33pm UTC](https://discourse.julialang.org/t/rank-is-wrong-for-rational-matrices/19380/3 "2019-01-08T13:33:08Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 8, 2019, 1:37pm UTC](https://discourse.julialang.org/t/rank-is-wrong-for-rational-matrices/19380/4 "2019-01-08T13:37:00Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Jean\_Michel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jean_michel/32/8282_2.png) [@Jean\_Michel](https://discourse.julialang.org/u/Jean_Michel)\
**Post date:** [January 8, 2019, 1:41pm UTC](https://discourse.julialang.org/t/rank-is-wrong-for-rational-matrices/19380/5 "2019-01-08T13:41:28Z")

</div>

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…

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 8, 2019, 1:51pm UTC](https://discourse.julialang.org/t/rank-is-wrong-for-rational-matrices/19380/6 "2019-01-08T13:51:37Z")

</div>

> [@Jean\_Michel](#):
>
> A proper algorithm would work for any field, like Rational Fractions, Algebraic Numbers, etc…

You could make a PR then for this “proper algorithm”.

---

<div class="post-metadata">

**Author:** ![Jean\_Michel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jean_michel/32/8282_2.png) [@Jean\_Michel](https://discourse.julialang.org/u/Jean_Michel)\
**Post date:** [January 8, 2019, 1:54pm UTC](https://discourse.julialang.org/t/rank-is-wrong-for-rational-matrices/19380/7 "2019-01-08T13:54:21Z")

</div>

As I said, a proper algorithm would be simply to compute the echelon form of the matrix.

---

<div class="post-metadata">

**Author:** ![Nosferican](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nosferican/32/9275_2.png) [@Nosferican](https://discourse.julialang.org/u/Nosferican)\
**Post date:** [January 8, 2019, 8:25pm UTC](https://discourse.julialang.org/t/rank-is-wrong-for-rational-matrices/19380/8 "2019-01-08T20:25:40Z")

</div>

> **[GitHub - blegat/RowEchelon.jl: Small package containing the rref fonction for...](https://github.com/blegat/RowEchelon.jl)**
>
> Small package containing the rref fonction for computing the reduced row echelon form of the matrix A - GitHub - blegat/RowEchelon.jl: Small package containing the rref fonction for computing the r...

---

<div class="post-metadata">

**Author:** ![Jean\_Michel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jean_michel/32/8282_2.png) [@Jean\_Michel](https://discourse.julialang.org/u/Jean_Michel)\
**Post date:** [January 8, 2019, 11:29pm UTC](https://discourse.julialang.org/t/rank-is-wrong-for-rational-matrices/19380/9 "2019-01-08T23:29:36Z")

</div>

The code there has a suspicious ε which may not make sense in every field. I was thinking of something like this:

```julia
" 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.

---

<div class="post-metadata">

**Author:** ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)\
**Post date:** [January 9, 2019, 3:19am UTC](https://discourse.julialang.org/t/rank-is-wrong-for-rational-matrices/19380/10 "2019-01-09T03:19:20Z")

</div>

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](http://nemocas.org/) where you may be happier to see that  
`rank(hilbert(MatrixSpace(QQ,n,n))) == n`.

---

<div class="post-metadata">

**Author:** ![Jean\_Michel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jean_michel/32/8282_2.png) [@Jean\_Michel](https://discourse.julialang.org/u/Jean_Michel)\
**Post date:** [January 9, 2019, 10:54pm UTC](https://discourse.julialang.org/t/rank-is-wrong-for-rational-matrices/19380/11 "2019-01-09T22:54:19Z")

</div>

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}}`.

---

<div class="post-metadata">

**Author:** ![abulak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abulak/32/28314_2.png) [@abulak](https://discourse.julialang.org/u/abulak)\
**Post date:** [January 10, 2019, 1:10am UTC](https://discourse.julialang.org/t/rank-is-wrong-for-rational-matrices/19380/12 "2019-01-10T01:10:16Z")

</div>

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
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.

Speaking of hostility:

> [@Jean\_Michel](#):
>
> As a mathematician, I do not think I have to take refuge in a ghetto as Nemo is currently.

well done! 😉

---

<div class="post-metadata">

**Author:** ![simonbyrne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonbyrne/32/19_2.png) [@simonbyrne](https://discourse.julialang.org/u/simonbyrne)\
**Post date:** [January 10, 2019, 1:29am UTC](https://discourse.julialang.org/t/rank-is-wrong-for-rational-matrices/19380/13 "2019-01-10T01:29:26Z")

</div>

> [@Jean\_Michel](#):
>
> As I said, a proper algorithm would be simply to compute the echelon form of the matrix.

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.

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [January 10, 2019, 3:02am UTC](https://discourse.julialang.org/t/rank-is-wrong-for-rational-matrices/19380/14 "2019-01-10T03:02:09Z")

</div>

It appears like the same problem applies to integer matrices:

```julia
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](http://pretty.structures.free.fr/talks/(7.1.1)Solving.pdf)?

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?

---

<div class="post-metadata">

**Author:** ![andreasnoack](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasnoack/32/27_2.png) [@andreasnoack](https://discourse.julialang.org/u/andreasnoack)\
**Post date:** [January 10, 2019, 7:45am UTC](https://discourse.julialang.org/t/rank-is-wrong-for-rational-matrices/19380/15 "2019-01-10T07:45:36Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Jean\_Michel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jean_michel/32/8282_2.png) [@Jean\_Michel](https://discourse.julialang.org/u/Jean_Michel)\
**Post date:** [January 10, 2019, 9:04am UTC](https://discourse.julialang.org/t/rank-is-wrong-for-rational-matrices/19380/16 "2019-01-10T09:04:22Z")

</div>

> [@andreasnoack](#):
>
> As people have already pointed out, the current `rank` function is for computing the floating point rank and when 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.

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 `BigInt`s should never be treated as floats). Has there been any thoughts along these lines among developers/designers of the language?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 10, 2019, 9:44am UTC](https://discourse.julialang.org/t/rank-is-wrong-for-rational-matrices/19380/17 "2019-01-10T09:44:54Z")

</div>

> [@Jean\_Michel](#):
>
> This is an important problem for me.

In that case, you should perhaps consider making a PR to get the discussion going with something concrete.

---

<div class="post-metadata">

**Author:** ![Jean\_Michel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jean_michel/32/8282_2.png) [@Jean\_Michel](https://discourse.julialang.org/u/Jean_Michel)\
**Post date:** [January 10, 2019, 10:31am UTC](https://discourse.julialang.org/t/rank-is-wrong-for-rational-matrices/19380/18 "2019-01-10T10:31:56Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![chakravala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chakravala/32/6832_2.png) [@chakravala](https://discourse.julialang.org/u/chakravala)\
**Post date:** [January 10, 2019, 10:45am UTC](https://discourse.julialang.org/t/rank-is-wrong-for-rational-matrices/19380/19 "2019-01-10T10:45:10Z")

</div>

> [@foobar\_lv2](#):
>
> It appears like the same problem applies to integer matrices:

You could try Hermite Normal Form

> **[Hermite normal form](https://en.m.wikipedia.org/wiki/Hermite_normal_form)**
>
> In linear algebra, the Hermite normal form is an analogue of reduced echelon form for matrices over the integers Z. Just as reduced echelon form can be used to solve problems about the solution to the linear system Ax=b where x is in Rn, the Hermite normal form can solve problems about the solution to the linear system Ax=b where this time x is restricted to have integer coordinates only. Other applications of the Hermite normal form include integer programming, cryptography, and abstract a Vario...

The Hecke.jl library supports it

[https://github.com/thofma/Hecke.jl/search?&q=hnf](https://github.com/thofma/Hecke.jl/search?&q=hnf)

As stated in the references [http://nemocas.org/links.html](http://nemocas.org/links.html)

---

<div class="post-metadata">

**Author:** ![thofma](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thofma/32/1691_2.png) [@thofma](https://discourse.julialang.org/u/thofma)\
**Post date:** [January 10, 2019, 10:58am UTC](https://discourse.julialang.org/t/rank-is-wrong-for-rational-matrices/19380/20 "2019-01-10T10:58:05Z")

</div>

No need for Hecke, [AbstractAlgebra.jl](https://github.com/Nemocas/AbstractAlgebra.jl) (pure julia without C dependencies) is enough:

```julia
julia> using AbstractAlgebra

julia> m = matrix(QQ, [1//(n+m) for n in 1:11, m in 1:11]);

julia> rank(m)
11

```

It also has a generic HNF implementation.