`lu(A::AbstractMatrix{Taylor1{T}})` method not properly inferred in Julia v1.8

[lbenet]

lu(A::AbstractMatrix{Taylor1{T}}) method here is not properly inferred in Julia v1.8 (currently nightly), which is the reason for this error in the tests. (The error is silent due to this line.)

The following illustrates the problem:

julia> versioninfo()
Julia Version 1.8.0-DEV.1221
Commit 98b485e064 (2022-01-05 22:00 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin19.6.0)
  CPU: Intel(R) Core(TM) i9-9880H CPU @ 2.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, skylake)

julia> using LinearAlgebra, SparseArrays

julia> t = Taylor1(5)
 1.0 t + 𝒪(t⁶)

julia> a = Diagonal(rand(0:10,3)) + rand(3, 3);

julia> b = Taylor1.(a, 5);

julia> @which lu(b)
lu(A::StridedMatrix{T} where T) in LinearAlgebra at /usr/local/julia/usr/share/julia/stdlib/v1.8/LinearAlgebra/src/lu.jl:79

For comparison, in Julia v1.7.1, the last line yields

julia> @which lu(b)
lu(A::AbstractArray{Taylor1{T}, 2}; check) where T<:Number in TaylorSeries at /Users/benet/.julia/dev/TaylorSeries/src/arithmetic.jl:629

Refs: #281, #284, and older #221, #223

@dkarrasch @daanhb Maybe you have a suggestion for this...

[lbenet]

A way to solve this is by extending the method to the signature lu(A::Union{AbstractMatrix{Taylor1{T}},StridedMatrix{Taylor1{T}}}; check...).

Still I don't understand why the method is not inferred at once....

[lbenet]

Or simply replacing AbstracMatrix{Taylor1{T}} by StridedMatrix{Taylor1{T}}...

[daanhb]

Hmm, I'm not sure what changed in v1.8, but indeed there is a potential for ambiguities when defining lu for your own matrix type. I think you may want to define both lu(A::AbstractMatrix{Taylor1{T}} and lu(A::StridedMatrix{Taylor1{T}} here, since LinearAlgebra defines fallbacks for both. Like your first suggestion. It is probably best to include the full signature, with the optional second argument including its default value. (Maybe try a HermOrSym too...)

I don't think this is an inference issue, it seems like a genuine ambiguity problem (which may have invisibly existed before v1.8).

@dkarrasch I did not think this through, but would it be an idea to add an indirection in LinearAlgebra? Something like lu(A) -> _lu(A) -> lu!(A). We could dispatch internally on _lu rather than on lu itself. This "frees up" lu a little more for specialization in other packages. Alternatively, the logic to choose between copy_oftype here and copy_to_array here could be in a generic subroutine, rather than in the main definition of lu itself. It would be good to have only one (generic) definition of lu in LinearAlgebra, because I fear that anyone wanting to specialize for a particular element type will run into the problem of this issue.

[dkarrasch]

Hmm, I'm not sure what changed in v1.8

On master we have JuliaLang/julia#40831 included, but not in v1.7.x.

Couldn't we remove the lu(::StridedMatrix, ...) method from LinearAlgebra.jl? If the eltype is not a BlasFloat (and it's not HermOrSym), then it would fall back to the lu(A::AbstractMatrix{T}, ...) method, which should give the exact same result.

[dkarrasch]

I'm working on mildly cleaning up JuliaLang/julia#40831, so we should continue to discuss there. I hope we can find solutions that would prevent the need for downstream fixes.

[daanhb]

Hmm, I'm not sure what changed in v1.8

On master we have JuliaLang/julia#40831 included, but not in v1.7.x.

Of course, that's it. I mistakenly thought that went into 1.7.

[lbenet]

Thanks a lot for looking up this @dkarrasch @daanhb! As proposed, the discussion will follow in JuliaLang/julia#43700.

(Maybe try a HermOrSym too...)

Just to answer @daanhb's comment, you are right, using #296 for a Symmetric{Taylor1{T}} matrix still uses the lu method of LinearAlgebra:

julia> a = Symmetric(rand(3,3));

julia> b = Symmetric(Taylor1.(a, 5))
3×3 Symmetric{Taylor1{Float64}, Matrix{Taylor1{Float64}}}:
  0.19349955324833334 + 𝒪(t⁶)  …    0.3400375943006272 + 𝒪(t⁶)
  0.47521570481828634 + 𝒪(t⁶)       0.5766819299947761 + 𝒪(t⁶)
   0.3400375943006272 + 𝒪(t⁶)      0.07220722813704294 + 𝒪(t⁶)

julia> @which lu(b)
lu(A::AbstractMatrix{T}) where T in LinearAlgebra at /usr/local/julia/julia-1.7.1/share/julia/stdlib/v1.7/LinearAlgebra/src/lu.jl:277

[dkarrasch]

The symmetric case should be taken care of by @daanhb' proposal which is now also part of JuliaLang/julia#43700.

[lbenet]

Fixed by JuliaLang/julia#43700.