If (in Julia) we compute the QR factorization of a rank deficient matrix like A=[1 2 3;4 5 6;7 8 9], some of the diagonal entries of the R matrix will be very small. However, when doing this numerically in Julia, how small must these diagonal entries be to be considered zero (thus rank deficient)? I am trying to find a formula for this rather than an absolute number.
1 Answers
If your original matrix A can be thought of as exactly rank deficient, but with some additive noise, the small-bug-nonzero elements of the diagonal of R that should be ignored will be about the same size as the noise in A.
For example,
julia> d = Diagonal(-5*log.(rand(5)))
5×5 Diagonal{Float64, Vector{Float64}}:
5.56654 ⋅ ⋅ ⋅ ⋅
⋅ 17.6294 ⋅ ⋅ ⋅
⋅ ⋅ 2.20874 ⋅ ⋅
⋅ ⋅ ⋅ 2.91583 ⋅
⋅ ⋅ ⋅ ⋅ 26.5664
julia> A = [d zeros(5,3); zeros(3,8)]
8×8 SparseArrays.SparseMatrixCSC{Float64, Int64} with 5 stored entries:
5.56654 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ 17.6294 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ 2.20874 ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 2.91583 ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ 26.5664 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
julia> diag(qr(A).R)
8-element SparseArrays.SparseVector{Float64, Int64} with 5 stored entries:
[1] = 5.56654
[2] = 17.6294
[3] = 2.20874
[4] = 2.91583
[5] = 26.5664
julia> diag(qr(A+rand(8,8)*1e-6).R)
8-element Vector{Float64}:
-5.566543507821689
-17.629411781406233
-2.20873750628613
-2.915828986361413
-26.56641073082019
-1.1241491098747558e-6
-6.122032725632311e-7
-2.3444092119580006e-7
Related