A bug in rank() function in package LinearAlgebra

[sunnsi]

julia> using LinearAlgebra
julia> df=ones(200,2);
julia> rank(df)
1
julia> df=ones(500,2);
julia> rank(df)
2
julia> versioninfo()
Julia Version 1.1.0
Commit 80516ca202 (2019-01-21 21:24 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-4790 CPU @ 3.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, haswell)
Environment:
  JULIAPATH = /home/sunsh/software/julia-1.1.0/bin

[andreasnoack]

See JuliaLang/julia#23523 and the linked issues. It's a bit subtle to determine the correct tolerance to use in this function. I did some simulations that suggested that the rounding errors were mainly influenced by the largersmaller dimension but that doesn't seem to hold for the matrix of ones:

julia> [prod(svdvals(ones(10^i,2))) for i in 1:6]
6-element Array{Float64,1}:
 4.213000162292042e-15
 2.9535326035845185e-14
 3.6564529777304436e-12
 1.0315077471517604e-11
 2.451559578683964e-10
 1.0685103408661952e-8

and I'm getting somewhat similar results for random matrices now. As mentioned in the linked issues, the problem with using the larger dimension is that ones(n,1) would have rank zero for sufficiently large n. It would be great if someone could derive a better tolerance as a function of m and n.

[eulerkochy]

Hey, notice this interesting observation. @andreasnoack

Why 287 ? It was the first number where rank of df was changing from 1 to 2.

[eulerkochy]

Check JuliaLang/julia#31103 It contains the heuristic I employed in resolving this issue
@StefanKarpinski

[simonbyrne]

There seems to be something up with LAPACKs algorithm, since:

julia> svdvals(ones(2^10,2))
2-element Array{Float64,1}:
 45.25483399593904
  4.24193471892493e-15

julia> svdvals(ones(2,2^10))
2-element Array{Float64,1}:
 45.254833995939535
  4.770749899131662e-13

Unfortunately the LAPACK docs aren't much help: all they say is that the error is a "modestly growing function of m and n"

[higham]

(4.770749899131662e-13/45.254833995939535)/2^10
=
1.0295e-17

so I don't think there is anything amiss.
To check correctness of the computed SVD it's best to check that U and V
are close to orthogonal and that the residual norm(A-U*S*V)/norm(A)
is small - then you know you have a small backward error.

[simonbyrne]

@higham I think our question is: what is an appropriate scaling factor for the error of the singular values, relative to the largest singular value. In other words, what is the approximate form of the function p(m,n) mentioned in the LAPACK docs.

Is max(m,n) a reasonable choice?

[higham]

Worst-case error constant is cubic in the dimensions.
Square root of that is more realistic: http://eprints.maths.manchester.ac.uk/2673/.
Beyond that, any choice is heuristic.