Potentially erring paths in matrix logarithms

[jishnub]

julia> using JET

julia> using LinearAlgebra

julia> @report_opt log(UpperTriangular(rand(4,4)))
# [some lines redacted for brevity]
│┌ log_quasitriu(A0::UpperTriangular{Float64, Matrix{Float64}}) @ LinearAlgebra /cache/build/builder-amdci4-5/julialang/julia-release-1-dot-11/usr/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:1711
││┌ _log_quasitriu!(A0::UpperTriangular{Float64, Matrix{Float64}}, A::UpperTriangular{ComplexF64, Matrix{ComplexF64}}) @ LinearAlgebra /cache/build/builder-amdci4-5/julialang/julia-release-1-dot-11/usr/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:1732
│││┌ _sqrt_pow_diag_quasitriu!(A::UpperTriangular{ComplexF64, Matrix{…}}, A0::UpperTriangular{Float64, Matrix{…}}, s::Int64) @ LinearAlgebra /cache/build/builder-amdci4-5/julialang/julia-release-1-dot-11/usr/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:1895
││││┌ _sqrt_pow_diag_block_2x2!(A::SubArray{…}, A0::SubArray{…}, s::Int64) @ LinearAlgebra /cache/build/builder-amdci4-5/julialang/julia-release-1-dot-11/usr/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:1939
│││││┌ _sqrt_real_2x2!(R::SubArray{…}, A::SubArray{…}) @ LinearAlgebra /cache/build/builder-amdci4-5/julialang/julia-release-1-dot-11/usr/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:2376
││││││┌ _real_sqrt(θ::ComplexF64, μ::Float64) @ LinearAlgebra /cache/build/builder-amdci4-5/julialang/julia-release-1-dot-11/usr/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:2388
│││││││┌ >=(x::ComplexF64, y::Int64) @ Base ./operators.jl:426
││││││││┌ <=(x::Int64, y::ComplexF64) @ Base ./operators.jl:402
│││││││││┌ <(x::Int64, y::ComplexF64) @ Base ./operators.jl:353
││││││││││ runtime dispatch detected: isless(x::Int64, y::ComplexF64)

Should the comparison be using the real part of the complex number? Perhaps someone familiar with this part of the code could comment on this.

Also,

julia> @report_opt log(rand(4,4))
# [some lines redacted for brevity]
┌ log(A::Matrix{Float64}) @ LinearAlgebra /cache/build/builder-amdci4-5/julialang/julia-release-1-dot-11/usr/share/julia/stdlib/v1.11/LinearAlgebra/src/dense.jl:844
│┌ log(A::UpperTriangular{ComplexF64, Matrix{ComplexF64}}) @ LinearAlgebra /cache/build/builder-amdci4-5/julialang/julia-release-1-dot-11/usr/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:1692
││┌ log_quasitriu(A0::UpperTriangular{ComplexF64, Matrix{ComplexF64}}) @ LinearAlgebra /cache/build/builder-amdci4-5/julialang/julia-release-1-dot-11/usr/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:1711
│││┌ _log_quasitriu!(A0::UpperTriangular{ComplexF64, Matrix{ComplexF64}}, A::UpperTriangular{ComplexF64, Matrix{ComplexF64}}) @ LinearAlgebra /cache/build/builder-amdci4-5/julialang/julia-release-1-dot-11/usr/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:1764
││││┌ _log_diag_quasitriu!(A::Matrix{ComplexF64}, A0::UpperTriangular{ComplexF64, Matrix{ComplexF64}}) @ LinearAlgebra /cache/build/builder-amdci4-5/julialang/julia-release-1-dot-11/usr/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:2041
│││││┌ _log_diag_block_2x2!(A::SubArray{…}, A0::SubArray{…}) @ LinearAlgebra /cache/build/builder-amdci4-5/julialang/julia-release-1-dot-11/usr/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:2060
││││││ runtime dispatch detected: LinearAlgebra.atan(%429::Float64, %127::ComplexF64)
│││││└────────────────────

atan isn't defined for this combination, and I wonder if this should be called for the real part again?

Unfortunately, I haven't been able to come up with examples of matrices that hit these code paths, but it would be good to fix these.

[Seelengrab]

There's definitely some erroring paths in log. I wasn't able to produce the errors from your investigation with JET, but I was able to produce this one:

julia> using LinearAlgebra

julia> [NaN 0.0; 0.0 0.0] |> UpperTriangular
2×2 UpperTriangular{Float64, Matrix{Float64}}:
 NaN    0.0
    ⋅   0.0

julia> [NaN 0.0; 0.0 0.0] |> UpperTriangular |> log
ERROR: InexactError: Float64(NaN + NaN*im)
Stacktrace:
  [1] Real
    @ ./complex.jl:44 [inlined]
  [2] convert
    @ ./number.jl:7 [inlined]
  [3] setindex!
    @ ./array.jl:979 [inlined]
  [4] setindex!
    @ ~/julia/usr/share/julia/stdlib/v1.12/LinearAlgebra/src/triangular.jl:286 [inlined]
  [5] _pow_superdiag_quasitriu!(A::UpperTriangular{Float64, Matrix{…}}, A0::UpperTriangular{Float64, Matrix{…}}, p::Float64)
    @ LinearAlgebra ~/julia/usr/share/julia/stdlib/v1.12/LinearAlgebra/src/triangular.jl:2088
  [6] _log_quasitriu!(A0::UpperTriangular{Float64, Matrix{Float64}}, A::UpperTriangular{Float64, Matrix{Float64}})
    @ LinearAlgebra ~/julia/usr/share/julia/stdlib/v1.12/LinearAlgebra/src/triangular.jl:1816
  [7] log_quasitriu(A0::UpperTriangular{Float64, Matrix{Float64}})
    @ LinearAlgebra ~/julia/usr/share/julia/stdlib/v1.12/LinearAlgebra/src/triangular.jl:1798
  [8] log
    @ ~/julia/usr/share/julia/stdlib/v1.12/LinearAlgebra/src/triangular.jl:1779 [inlined]
  [9] |>(x::UpperTriangular{Float64, Matrix{Float64}}, f::typeof(log))
    @ Base ./operators.jl:967
 [10] top-level scope
    @ REPL[38]:1
Some type information was truncated. Use `show(err)` to see complete types.

Being unfamiliar with the math involved here and the docs not mentioning anything about NaN in the input matrix, I have no idea whether that is intended though :/ At the very least, it's a hard-to-diagnose/reconcile error.

Here's the fuzzing setup that found this:

julia> using Supposition, LinearAlgebra

julia> square_mats = Data.SquareMatrices(Data.Floats{Float64}(); max_size=20) # upcoming Possibility in the next release
Supposition.Data.SquareMatrices(Supposition.Data.Vectors(Supposition.Data.Floats{Float64}(; nans=true, infs=true); min_size=0, max_size=400); min_size=0, max_size=20)

julia> mats = map(UpperTriangular, square_mats);

julia> @check db=false (m=mats) -> log(m) isa UpperTriangular
Encountered an error
  Context: ##SuppositionAnon#646

  Arguments:
      m::UpperTriangular{Float64, Matrix{Float64}} = [NaN 0.0; 0.0 0.0]

  Exception:
    Message:
      InexactError: Float64(NaN + NaN*im)

If I catch that InexactError, the test passes even when I crank it up to 1_000_000 examples. The same goes for disabling NaN in the generation of data.