Why does Julia fails to solve linear system systematically?

Question

The problem Ax=b for square A is solved by the \ function. With that in mind, I've tried to do the following:

A = rand(1:4,3,3)
x = fill(1.0, 3)
b = A * x 
A\b

For some reason, the code seems to works at times. But sometimes it returns me the following error:

LinearAlgebra.SingularException(3)

Stacktrace:
 [1] checknonsingular
   @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/factorization.jl:19 [inlined]
 [2] checknonsingular
   @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/factorization.jl:21 [inlined]
 [3] #lu!#136
   @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:85 [inlined]
 [4] #lu#140
   @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:273 [inlined]
 [5] lu (repeats 2 times)
   @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:272 [inlined]
 [6] \(A::Matrix{Int64}, B::Vector{Float64})
   @ LinearAlgebra /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/generic.jl:1136
 [7] top-level scope
   @ In[208]:4
 [8] eval
   @ ./boot.jl:360 [inlined]
 [9] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1116

So, I tried to understand what is happening, and executed the code 10000000 times and found out that it failed 10% of the times it was executed.

using Printf

i = 0
test = 10000000
for x in 1:test    
    try
        A = rand(1:4,3,3)
        x = fill(1.0, 3)
        b = A * x 
        A\b
    catch 
        i = i+1        
    end
end

fail_percentage = (i/test)*100

@printf "this code has failed in %.2f%%" fail_percentage

Can someone explain me what is happening here?

Answer 1

The error is explicit: LinearAlgebra.SingularException. This is not a failure of Julia, but a property of a system of equations.

There is no single solution if the matrix A is singular - either an infinite amount of solutions if the system is homogeneous, or none in the general case. Seems you have empirically calculated the probability of generating a singular system using the properties you tested (dimensions of A and x, x filled with 1s, A filled between 1 and 4).

In case, like OP, you are looking to skip singular matrices, you need to ensure determinant of A is not 0. You can either use the built in function to check and skip such matrices, or generate them by noting that the determinant is itself an equation, and so, say for your 3x3 example with no 0 entries, you choose 8 numbers, you can calculate what the 9th one cannot be to ensure the determinant is non-zero. If you allow for 0s you need to check all possibilities.