Skip to content

Conversation

@dkarrasch
Copy link
Member

This is a careful manual rebase of #35420. I'd say half of it was already done in other PRs, so here is what remained to be done. There was one concern by @andreasnoack that some implementations may now be open to array types that (i) are invariant under 3-arg similar and (ii) have slow elementwise indexing. If that is still a concern, then this PR should be reviewed in light of that. This needs a pkgeval run to see if it creates many ambiguities, but, as you can see, there was not that much to be resolved. I wonder if this [c/sh]ould be made a v1.9 milestone, so packages could iron out ambiguities in the beta phase?

Closes #35420. Fixes JuliaLang/LinearAlgebra.jl#713. Fixes JuliaLang/LinearAlgebra.jl#620. Fixes JuliaLang/LinearAlgebra.jl#587.

Comment on lines -351 to -355
# There is no fallback solver for Bunch-Kaufman so we'll have to promote to same element type
function ldiv!(B::BunchKaufman{T}, R::StridedVecOrMat{S}) where {T,S}
TS = promote_type(T,S)
return ldiv!(convert(BunchKaufman{TS}, B), convert(AbstractArray{TS}, R))
end
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is an untested method that doesn't even follow the ldiv! contract: when the types TS and S are different, then R gets copied. Usually, we don't promote eltypes in in-place functions.

@dkarrasch
Copy link
Member Author

@nanosoldier runtests(ALL, vs = ":master")

@nanosoldier
Copy link
Collaborator

Your package evaluation job has completed - possible new issues were detected. A full report can be found here.

@dkarrasch
Copy link
Member Author

Wow, what a clean nanosoldier run! ArrayLayouts.jl "fails" due to unexpected passes. 🤣 BlockArrays.jl fails due to

MethodError: ldiv!(::LU{Float64, Matrix{Float64}, Vector{Int64}}, ::PseudoBlockMatrix{Float64, Matrix{Float64}, Tuple{BlockedUnitRange{RangeCumsum{Int64, UnitRange{Int64}}}, BlockedUnitRange{Vector{Int64}}}}) is ambiguous.
  
  Candidates:
    ldiv!(A::Factorization, x::LayoutMatrix)
      @ ArrayLayouts ~/.julia/packages/ArrayLayouts/DYIpy/src/ldiv.jl:135
    ldiv!(A::LU, B::AbstractVecOrMat)
      @ LinearAlgebra /opt/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/lu.jl:426
  
  Possible fix, define
    ldiv!(::LU, ::LayoutMatrix)

which points to a general pattern for all Factorizations that admit a generic ldiv! solver (in which we relaxed the StridedVecOrMat restriction). ExtendableSparse.jl can be easily fixed by splitting \(symm_ext::Symmetric{Tm, ExtendableSparseMatrix{Tm, Ti}}, B::AbstractVecOrMat) where {Tm, Ti} into two methods, one for matrices and one for vectors; and same for Hermitian. InfiniteLinearAlgebra fails due to the same reason as ArrayLayouts.jl (actually it's the same ArrayLayouts.jl issue with Cholesky. So, it seems to be all about the ldiv!(::Factorization, ::LayoutMatrix) in ArrayLayouts.jl. Finally, I'm not sure what to do with OptimKit.jl. The computation doesn't seem to converge, so perhaps @Jutho could take a look at this PR if there is something relevant?

@Jutho
Copy link
Contributor

Jutho commented Oct 10, 2022

I don't think the OptimKit.jl failure is relevant. I run some optimisation tests using quadratic problems with random matrices, and do not even set the seed appropriately. So I guess it can happen that a matrix happens to have bad convergence properties for GD, and convergence is thus terribly slow. I should probably make the tests of OptimKit a bit more robust.

@dkarrasch
Copy link
Member Author

@dlfivefifty I think this can be merged once we find a solution for the ambiguity in ArrayLayouts.jl, see #47107 (comment). Do you have any idea how to fix it? A suboptimal fix would be to introduce specific methods for all factorizations from LinearAlgebra, but that wouldn't extend to other factorizations from the general ecosystem.

OTOH, I noticed that ArrayLayouts.jl has about 1k method ambiguities (counted by Aqua.jl), so I wonder how it can be working at all. 😅

@dkarrasch
Copy link
Member Author

Hm:

               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.8.2 (2022-09-29)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> using LinearAlgebra

julia> methods(ldiv!, (Factorization, AbstractMatrix), LinearAlgebra)
# 19 methods for generic function "ldiv!":
[1] ldiv!(A::LinearAlgebra.QRCompactWY{T, M, C} where {M<:AbstractMatrix{T}, C<:AbstractMatrix{T}}, B::StridedMatrix{T}) where T<:Union{Float32, Float64, ComplexF32, ComplexF64} in LinearAlgebra at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/qr.jl:876
[2] ldiv!(A::QRPivoted{T, S, C} where {S<:AbstractMatrix{T}, C<:AbstractVector{T}}, B::StridedVecOrMat{T}) where T<:Union{Float32, Float64, ComplexF32, ComplexF64} in LinearAlgebra at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/qr.jl:920
[3] ldiv!(A::QRPivoted, B::StridedMatrix{T} where T) in LinearAlgebra at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/qr.jl:989
[4] ldiv!(A::LU{T, <:StridedMatrix{T} where T}, B::StridedVecOrMat{T}) where T<:Union{Float32, Float64, ComplexF32, ComplexF64} in LinearAlgebra at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/lu.jl:404
[5] ldiv!(A::LU{<:Any, <:StridedMatrix{T} where T}, B::StridedVecOrMat) in LinearAlgebra at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/lu.jl:407
[6] ldiv!(A::LU{T, Tridiagonal{T, V}}, B::AbstractVecOrMat) where {T, V} in LinearAlgebra at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/lu.jl:598
[7] ldiv!(C::CholeskyPivoted{T}, B::StridedMatrix{T}) where T<:Union{Float32, Float64, ComplexF32, ComplexF64} in LinearAlgebra at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/cholesky.jl:592
[8] ldiv!(C::CholeskyPivoted, B::AbstractMatrix) in LinearAlgebra at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/cholesky.jl:615
[9] ldiv!(A::QR{T, S, C} where {S<:AbstractMatrix{T}, C<:AbstractVector{T}}, B::StridedMatrix{T}) where T in LinearAlgebra at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/qr.jl:970
[10] ldiv!(A::SVD{T, Tr, M} where {Tr, M<:(AbstractArray{T})}, B::StridedVecOrMat) where T in LinearAlgebra at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/svd.jl:251
[11] ldiv!(C::Cholesky{T, <:StridedMatrix{T} where T}, B::StridedVecOrMat{T}) where T<:Union{Float32, Float64, ComplexF32, ComplexF64} in LinearAlgebra at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/cholesky.jl:578
[12] ldiv!(C::Cholesky, B::StridedVecOrMat) in LinearAlgebra at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/cholesky.jl:581
[13] ldiv!(A::LQ, B::StridedVecOrMat) in LinearAlgebra at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/lq.jl:337
[14] ldiv!(S::LDLt{<:Any, <:SymTridiagonal}, B::AbstractVecOrMat) in LinearAlgebra at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/ldlt.jl:174
[15] ldiv!(F::Hessenberg{<:Complex, <:Any, <:AbstractMatrix{<:Real}}, B::AbstractVecOrMat{<:Complex}) in LinearAlgebra at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/hessenberg.jl:662
[16] ldiv!(F::Hessenberg, B::AbstractVecOrMat) in LinearAlgebra at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/hessenberg.jl:646
[17] ldiv!(B::BunchKaufman{T}, R::StridedVecOrMat{T}) where T<:Union{Float32, Float64} in LinearAlgebra at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/bunchkaufman.jl:329
[18] ldiv!(B::BunchKaufman{T}, R::StridedVecOrMat{T}) where T<:Union{ComplexF32, ComplexF64} in LinearAlgebra at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/bunchkaufman.jl:336
[19] ldiv!(B::BunchKaufman{T}, R::StridedVecOrMat{S}) where {T, S} in LinearAlgebra at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/bunchkaufman.jl:352

So in v1.8 we have

[6] ldiv!(A::LU{T, Tridiagonal{T, V}}, B::AbstractVecOrMat)
[8] ldiv!(C::CholeskyPivoted, B::AbstractMatrix)
[14] ldiv!(S::LDLt{<:Any, <:SymTridiagonal}, B::AbstractVecOrMat)
[16] ldiv!(F::Hessenberg, B::AbstractVecOrMat)

that are already causing ambiguities with ArrayLayouts.jl. Apparently, they are not called anywhere with a LayoutMatrix in the registered ecosystem, otherwise that would throw or tests would fail. I guess we only have two options: restrict to StridedArrays here or introduce a few disambiguating ldiv! methods in ArrayLayouts.jl.

@dkarrasch
Copy link
Member Author

@nanosoldier runtests(["ArrayLayouts", "BlockArrays", "ExtendableSparse", "InfiniteLinearAlgebra", "NeXLCore", "OptimKit", "Singular"], vs = ":master")

@nanosoldier
Copy link
Collaborator

Your package evaluation job has completed - no new issues were detected. A full report can be found here.

@dkarrasch
Copy link
Member Author

@nanosoldier runtests(["ArrayLayouts", "BlockArrays", "ExtendableSparse", "InfiniteLinearAlgebra"], vs = ":master")

@nanosoldier
Copy link
Collaborator

Your package evaluation job has completed - no new issues were detected. A full report can be found here.

@dkarrasch
Copy link
Member Author

Alright, this seems to fix broken tests instead of creating new issues, so this looks promising. I'll put PRs together to adjust the tests in ArrayLayouts.jl and BandedMatrices.jl.

@dkarrasch
Copy link
Member Author

CI failures seem unrelated.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

linear algebra Linear algebra

Projects

None yet

5 participants