Open jishnub opened 3 months ago
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.
In this comment, I am just analysing the first case log(UpperTriangular(rand(4,4))
.
UpperTriangular(rand(4,4))
has only nonnegative elements.A
, i.e., _log_quasitriu!(A0::UpperTriangular{Float64, Matrix{Float64}}, A::UpperTriangular{ComplexF64, Matrix{ComplexF64}})
is chosen in the function log_quasitriu
https://github.com/JuliaLang/julia/blob/fdc109088e7a25c7b412546634f2db70a2d584ad/stdlib/LinearAlgebra/src/triangular.jl#L1920isreal(A0)
is true, istriu(A0)
is true, and any(x -> real(x) < zero(real(T)), diag(A0))
is false, since A0
only has nonnegative elements as described above. Hence true && (!true || !false)
is true. Therefore, the chosen code branch is https://github.com/JuliaLang/julia/blob/fdc109088e7a25c7b412546634f2db70a2d584ad/stdlib/LinearAlgebra/src/triangular.jl#L1924 This cannot result in A
being complex since A0
is real valued.@report_opt
chooses an incorrect branch leading to the above confusion. Note: I can reproduce the reported incorrect behaviour of @report_opt
on Version 1.10.5 (2024-08-27)
.In this comment, I analyse log(rand(4,4))
. Here too, the path
││││┌ _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
is invalid. This is because in _log_diag_quasitriu!
with eltype(A)
and eltype(A0)
as ComplexF64
, the check in
https://github.com/JuliaLang/julia/blob/e217f938d291854caaca4e5b81bbeb27c3e985ba/stdlib/LinearAlgebra/src/triangular.jl#L2254 will always be true
since the schur-decomposition output matrix A
is upper triangular for complex matrices (not quasi triangular) and the else
path
https://github.com/JuliaLang/julia/blob/e217f938d291854caaca4e5b81bbeb27c3e985ba/stdlib/LinearAlgebra/src/triangular.jl#L2279-L2282 will never be taken for complex inputs.
Both provided paths cannot be reached in normal operation. Hence, I suggest to review my analyses and, if acceptable, close this issue.
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,
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.