Question about singular matrix inversion

General Usage
1

Hello :grinning:,
I was wondering why inversion works (when it shouldn’t) with floats, but returns an error with integers? Is it a rounding issue?

julia> inv([0.23 0.23; 0.3 0.3])
2×2 Matrix{Float64}:
 -3.3777e17   2.58957e17
  3.3777e17  -2.58957e17

julia> inv([1 1; 2 2])
ERROR: SingularException(2)
Stacktrace:
  [1] checknonsingular
    @ ~/.julia/juliaup/julia-1.11.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/factorization.jl:69 [inlined]
  [2] _check_lu_success
    @ ~/.julia/juliaup/julia-1.11.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/lu.jl:84 [inlined]
  [3] #lu!#182
    @ ~/.julia/juliaup/julia-1.11.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/lu.jl:92 [inlined]
  [4] lu!
    @ ~/.julia/juliaup/julia-1.11.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/lu.jl:90 [inlined]
  [5] lu!
    @ ~/.julia/juliaup/julia-1.11.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/lu.jl:89 [inlined]
  [6] _lu
    @ ~/.julia/juliaup/julia-1.11.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/lu.jl:347 [inlined]
  [7] lu(::Matrix{Int64}; kwargs::@Kwargs{})
    @ LinearAlgebra ~/.julia/juliaup/julia-1.11.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/lu.jl:341
  [8] lu
    @ ~/.julia/juliaup/julia-1.11.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/lu.jl:341 [inlined]
  [9] inv(A::Matrix{Int64})
    @ LinearAlgebra ~/.julia/juliaup/julia-1.11.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/dense.jl:993
 [10] top-level scope
    @ ~/Documents/julia_code/DLI_project/blend_functions.jl:313

Merci!,
fdekerm

2

Correct. You can watch what happens in more detail by breaking down the inversion into it’s constituent steps. Start by computing the LU factorization of the matrix:

ulia> A = [0.23 0.23; 0.3 0.3]
2×2 Matrix{Float64}:
 0.23  0.23
 0.3   0.3

julia> (; P, L, U) = lu(A)
LU{Float64, Matrix{Float64}, Vector{Int64}}
L factor:
2×2 Matrix{Float64}:
 1.0       0.0
 0.766667  1.0
U factor:
2×2 Matrix{Float64}:
 0.3  0.3
 0.0  2.96059e-18

Look at the last diagonal element in U: it should be zero, but we got 2.96e-18. Compared to the other numbers here, 2.96e-18 is really close to zero; in fact, it’s closer to zero than the actual number hiding behind 0.3 is to the real number 0.3. (Try big"0.3" - 0.3 to see this.) But since it’s not exactly zero, the remaining steps of the inversion will succeed and give a bogus answer:

julia> Ainv = U \ (L \ P)
2×2 Matrix{Float64}:
 -3.3777e17   2.58957e17
  3.3777e17  -2.58957e17

The reason your second example fails as expected is not due to the number type; in fact, the type must be converted to float before the factorization can proceed anyway. The difference is simply that for these particular numbers, the algorithm used for LU factorization happens to return an exact instead of approximate zero pivot:

julia> (; P, L, U) = lu([1.0 1.0; 2.0 2.0]; check=false)
Failed factorization of type LU{Float64, Matrix{Float64}, Vector{Int64}}

julia> U
2×2 Matrix{Float64}:
 2.0  2.0
 0.0  0.0
4 Likes
3

The reason for this is that 1.0 and 2.0 and 1.0/2.0 are represented exactly as floating-point numbers, so Gaussian elimination on [1 1; 2 2] incurs no roundoff errors.

PS. On my machine, inv([0.23 0.23; 0.3 0.3]) happens to gives a SingularException too.

2 Likes