What exactly does the rank function compute? I thought it would give the number of llinearly independent rows in a matrix, but evidently not.
Here cov is a 138x138 matrix.
1|julia> rank(cov[1:20,1:20])
20
1|julia> rank(cov)
13
Shouldn’t that be impossible?
Indeed, this looks weird.
At first I thought that you may have introduced some problems by using the name of the function cov from the standard Statistics library, but unless you import it, you are free to call your matrix cov.
Do you think it is possible to share you matrix somehow? Do you generate it through a code? In general, it is really useful to share a minimum working example (including the imports/using of the packages) because who knows where you take your function rank from. I certainly cannot reproduce the problem with some randomly generated singular matrix.
julia> using LinearAlgebra: rank
julia> A = rand(13,138);
julia> X = A'*A
138×138 Matrix{Float64}:
6.37017 3.87032 3.7639 2.93673 3.90856 2.89665 … 3.58424 4.81491 2.77375 3.95746 4.24206
3.87032 4.4881 3.17356 3.16758 3.26719 2.35837 3.2946 4.02962 2.37569 3.14432 2.57209
3.7639 3.17356 4.64649 3.36158 2.97343 2.63004 4.15317 4.36145 3.35719 2.9598 4.12072
2.93673 3.16758 3.36158 3.80284 2.74161 2.42185 3.45357 3.46843 2.79067 2.18312 2.65219
3.90856 3.26719 2.97343 2.74161 3.84898 2.41936 3.0891 3.70304 2.32024 2.98657 3.03344
2.89665 2.35837 2.63004 2.42185 2.41936 2.40715 … 2.71721 3.39881 2.12739 2.21862 2.54849
5.21516 4.57535 4.43188 4.21679 3.69402 3.22894 4.32022 5.20635 3.58125 4.07824 4.01366
⋮ ⋮ ⋱ ⋮
4.42637 4.1324 3.98606 3.16609 3.45972 2.20833 3.65279 4.18168 2.93238 3.29186 3.46011
3.58424 3.2946 4.15317 3.45357 3.0891 2.71721 4.12401 4.17301 3.20623 2.70909 3.51342
4.81491 4.02962 4.36145 3.46843 3.70304 3.39881 4.17301 5.53489 3.24436 3.89468 4.29279
2.77375 2.37569 3.35719 2.79067 2.32024 2.12739 … 3.20623 3.24436 2.90164 2.27987 2.93651
3.95746 3.14432 2.9598 2.18312 2.98657 2.21862 2.70909 3.89468 2.27987 3.79522 3.11997
4.24206 2.57209 4.12072 2.65219 3.03344 2.54849 3.51342 4.29279 2.93651 3.11997 4.53387
julia> rank(X)
13
julia> rank(X[1:20,1:20])
13
What’s a good way to share a 138x138 matrix?
Btw, could this matrix just be too large for the LinearAlgebra package?
@unhandyandy: rank(A) computes the number of singular values of A that are above a given tolerance.
help?> rank
search: rank RankDeficientException lowrankupdate lowrankupdate! lowrankdowndate
rank(A::AbstractMatrix; atol::Real=0, rtol::Real=atol>0 ? 0 : n*ϵ)
rank(A::AbstractMatrix, rtol::Real)
Compute the rank of a matrix by counting how many singular values of A have
magnitude greater than max(atol, rtol*σ₁) where σ₁ is A's largest singular
value. atol and rtol are the absolute and relative tolerances, respectively.
The default relative tolerance is n*ϵ, where n is the size of the smallest
dimension of A, and ϵ is the eps of the element type of A.
| John_Gibson
September 6 |
- | - |
`
`
Try computing the singular values of your 20 x 20 submatrix. I’ll bet seven of them are very close to zero., meaning your matrix is effectively rank 13 in 64-bit floating point artithmetic.
@zdenek_hurak: Did you mean to type
A = rand(138,138)instead ofA = rand(13,138)?
No, he used A to form A’*A, which he knew would have rank 13.
Here’s a simple example illustrating the behavior you see with your cov matrix
julia> A = randn(5,5);
julia> A[1:5,1:3] = 1e-32*randn(5,3);
julia> rank(A)
2
julia> rank(A[1:3,1:3])
3
The issue is that the submatrix is well-conditioned but the full matrix is not:
julia> svdvals(A)
5-element Vector{Float64}:
2.8403134593841486
1.883573969215261
1.5494894381007906e-16
9.045455141730465e-17
1.0078154590482648e-32
julia> svdvals(A[1:3,1:3])
3-element Vector{Float64}:
3.9530622919506925e-32
1.2513920009022788e-32
1.9253723883320212e-33
julia> cond(A[1:3,1:3])
20.531416758164323
julia> cond(A)
2.8182872507893575e32
Sorry about earlier edits. I misread your original post and didn’t see the issue correcty.
I see, thank you.
Perfect. Thanks. I have composed even a tiny bit more explicit example that demonstrates the phenomenon:
julia> A = [1e-16 0 0; 0 1e-16 0; 0 0 1.0]
3×3 Matrix{Float64}:
1.0e-16 0.0 0.0
0.0 1.0e-16 0.0
0.0 0.0 1.0
julia> rank(A)
1
julia> rank(A[1:2,1:2])
2
I confess I have learnt something. Another confirmation that one shouldn’t blindly rely on the (discontinuous) rank function – checking the singular values is a good practice.
That’s even better!