SingularException when dividing lu on zero-dimensional matrices of some generic types

[lukas-weber]

Many parts of julia linear algebra handle empty matrices gracefully, which is helpful when writing codes where they appear as a corner case. However, dividing LU-decompositions of generic numbers will yield a SingularException instead.

using LinearAlgebra
import ForwardDiff: Dual

a = fill(Dual(1.0, 1.0), 0, 0)
v = fill(Dual(1.0, 1.0), 0, 10) 

a \ v # return 0x0 matrix as expected

lu(a) \ v # throws SingularException(0)

ERROR: SingularException(0)
Stacktrace:
 [1] generic_trimatdiv!(C::Matrix{Dual{…}}, uploc::Char, isunitc::Char, tfun::typeof(identity), A::Matrix{Dual{…}}, B::Matrix{Dual{…}})
   @ LinearAlgebra ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/triangular.jl:1226
 [2] _ldiv!
   @ ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/triangular.jl:833 [inlined]
 [3] ldiv!
   @ ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/triangular.jl:838 [inlined]
 [4] ldiv!(A::LU{Dual{Nothing, Float64, 1}, Matrix{Dual{Nothing, Float64, 1}}, Vector{Int64}}, B::Matrix{Dual{Nothing, Float64, 1}})
   @ LinearAlgebra ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:498
 [5] ldiv(F::LU{Dual{Nothing, Float64, 1}, Matrix{Dual{Nothing, Float64, 1}}, Vector{Int64}}, B::Matrix{Dual{Nothing, Float64, 1}})
   @ LinearAlgebra ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/LinearAlgebra.jl:646
 [6] \(F::LU{Dual{Nothing, Float64, 1}, Matrix{Dual{Nothing, Float64, 1}}, Vector{Int64}}, B::Matrix{Dual{Nothing, Float64, 1}})
   @ LinearAlgebra ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/LinearAlgebra.jl:621
 [7] top-level scope
   @ REPL[6]:1

julia> versioninfo()
Julia Version 1.12.0-DEV.90
Commit b3b27360672 (2024-02-27 13:57 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 32 × Intel(R) Xeon(R) Gold 6244 CPU @ 3.60GHz
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, cascadelake)
Threads: 1 default, 0 interactive, 1 GC (on 32 virtual cores)
Environment:
  LD_LIBRARY_PATH = /mnt/sw/nix/store/pmwk60bp5k4qr8vsg411p7vzhr502d83-openblas-0.3.23/lib:/mnt/home/lweber/app/cpython/prefix/lib:/mnt/home/lweber/app/cpython/prefix/lib:/mnt/home/lweber/app/cpython/prefix/lib:/mnt/home/lweber/app/cpython/prefix/lib:/mnt/home/lweber/app/cpython/prefix/lib:/mnt/home/lweber/app/cpython/prefix/lib:/mnt/home/lweber/app/cpython/prefix/lib:/cm/shared/apps/slurm/current/lib64:/mnt/home/lweber/app/cpython/prefix/lib:/mnt/home/lweber/app/cpython/prefix/lib:/mnt/home/lweber/app/cpython/prefix/lib:/mnt/home/lweber/app/cpython/prefix/lib:/mnt/home/lweber/app/cpython/prefix/lib:/mnt/home/lweber/app/cpython/prefix/lib:/mnt/home/lweber/app/cpython/prefix/lib

It is worth noting that doing the same with complex Dual numbers on 1.10.1 will lead to segfaults as well. In the nightly, I was not able to reproduce that so I assume it has been fixed.

using LinearAlgebra
import ForwardDiff: Dual

a = fill(complex(Dual(1.0, 1.0)), 0, 0)
v = fill(complex(Dual(1.0, 1.0)), 0, 10) 

a \ v # return 0x0 matrix as expected

lu(a) \ v # throws SingularException(0)

ERROR: SingularException(0)
Stacktrace:
 [1] generic_trimatdiv!(C::Matrix{Complex{Dual{…}}}, uploc::Char, isunitc::Char, tfun::typeof(identity), A::Matrix{Complex{Dual{…}}}, B::Matrix{Complex{Dual{…}}})
   @ LinearAlgebra ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/triangular.jl:1226
 [2] _ldiv!
   @ ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/triangular.jl:833 [inlined]
 [3] ldiv!
   @ ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/triangular.jl:838 [inlined]
 [4] ldiv!(A::LU{Complex{Dual{Nothing, Float64, 1}}, Matrix{Complex{Dual{Nothing, Float64, 1}}}, Vector{Int64}}, B::Matrix{Complex{Dual{Nothing, Float64, 1}}})
   @ LinearAlgebra ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:498
 [5] ldiv(F::LU{Complex{Dual{Nothing, Float64, 1}}, Matrix{Complex{Dual{Nothing, Float64, 1}}}, Vector{Int64}}, B::Matrix{Complex{Dual{Nothing, Float64, 1}}})
   @ LinearAlgebra ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/LinearAlgebra.jl:646
 [6] \(F::LU{Complex{Dual{Nothing, Float64, 1}}, Matrix{Complex{Dual{Nothing, Float64, 1}}}, Vector{Int64}}, B::Matrix{Complex{Dual{Nothing, Float64, 1}}})
   @ LinearAlgebra ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/LinearAlgebra.jl:621
 [7] top-level scope
   @ REPL[6]:1
Some type information was truncated. Use `show(err)` to see complete types.

[aravindh-krishnamoorthy]

The line causing the exception on the latest master 71f68b4c is as follows:

julia/stdlib/LinearAlgebra/src/triangular.jl

Line 1226 in 71f68b4

iszero(amm) && throw(SingularException(m))

[dkarrasch]

@lukas-weber Could you please post the full stacktraces and a copy-pastable code example for the complex Dual case. One possibility is to catch this case earlier at the lu/factorization level, another one is to catch it at the generic triangular solve level. That would be easier to judge given the stacktraces.

[aravindh-krishnamoorthy]

Perhaps it's better to solve it in triangular solve as it will also be solved for the vanilla L \ b when used independently of lu? What I saw at that time was that in the other "graceful" cases, it works because the loop variables are not accessed outside of the loop over their size. But, in this case, amm is accessed before the loop over m, perhaps for efficiency reasons.

[lukas-weber]

@lukas-weber Could you please post the full stacktraces and a copy-pastable code example for the complex Dual case. One possibility is to catch this case earlier at the lu/factorization level, another one is to catch it at the generic triangular solve level. That would be easier to judge given the stacktraces.

They are added now. Sorry for the inconvenience!