Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Don't subtype AbstractQ <: AbstractMatrix #46196

Merged
merged 27 commits into from Feb 16, 2023
Merged

Don't subtype AbstractQ <: AbstractMatrix #46196

merged 27 commits into from Feb 16, 2023

Conversation

dkarrasch
Copy link
Member

This is a technically breaking experiment. I believe there is very good reason to not subtype AbstractQ as AbstractMatrix. The main reason is... it isn't! Q's are very different from container-like matrices, they are more like container-backed, function-based linear operators. I believe we've had avoidable issues in the past due to automatically allowing broadcast and indexing, see, for instance, #38972 (comment), and method ambiguities, see the deletions in special.jl in this PR. While this may seem weird to some people, let me point out that we have nice precedence for this: AbstractRotation and Factorization, which both behave like matrices in many ways, but also don't subtype AbstractMatrix! Here, I went even one step further and introduced a minimal AdjointQ <: AbstractQ type to avoid making Q' an AbstractMatrix again. In fact, I believe the same would be reasonable for (R::AbstractRotation)', i.e., to introduce a minimal AdjointRotation type to get these things out of the matrix dispatch (see the disambiguation methods at the end of givens.jl!).

I'm not sure such a change in the type system is acceptable at all, but to the very least, we should be able to detect very bad usage of AbstractQ objects in a AbstractMatrix context (especially like in broadcasting) via a nanosoldier run, but I don't believe this will be widespread, more like by accident.

@dkarrasch dkarrasch added kind:breaking This change will break code domain:linear algebra Linear algebra kind:speculative Whether the change will be implemented is speculative labels Jul 27, 2022
@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.

kellertuer pushed a commit to JuliaManifolds/Manopt.jl that referenced this pull request Jul 29, 2022
This came up in a nanosoldier run in JuliaLang/julia#46196. Transposes of Q's don't have specialized functions and fall back to `AbstractMatrix` methods, in particular in multiplication. Generic multiplication, however, is defined in terms of elementwise `getindex`, which is expensive and very slow in contrast to dense or structured matrices.
@dkarrasch
Copy link
Member Author

@nanosoldier runtests(["AKNS", "ASE", "AbstractAlgebra", "ArrayLayouts", "BasicInterpolators", "BifurcationInference", "BosonSampling", "BugReporting", "COPT", "CairoMakie", "ControlSystems", "CountdownNumbers", "CrystalInfoFramework", "DataDeps", "DifferentiableTrajectoryOptimization", "DiffinDiffsBase", "DrelTools", "FastMarching", "FileIO", "FilesystemDatastructures", "FiniteDifferences", "FlameGraphs", "Folds", "FundamentalsNumericalComputation", "GeneNetworkAPI", "GenericLinearAlgebra", "GeoParquet", "GeometricIntegrators", "GraphMLDatasets", "GraphNeuralNetworks", "IncrementalPruning", "InfiniteArrays", "InformationGeometry", "Inpaintings", "InteractBase", "IntervalLinearAlgebra", "IscaTools", "JOLI", "JSONLines", "KernelEstimator", "KissMCMC", "KrylovMethods", "LIBSVM", "Lehmann", "LineSearches", "LinearRationalExpectations", "LoggingExtras", "LowRankArithmetic", "Lux", "MPSGE", "ManifoldLearning", "Manopt", "MatrixMarket", "Metalhead", "MultiscaleGraphSignalTransforms", "MultivariatePolynomials", "NLSolversBase", "NOMAD", "Nemo", "NiceNumbers", "NonconvexMultistart", "Optim", "OptimKit", "OptimTestProblems", "OptimizationAlgorithms", "Oscar", "ParameterSpacePartitions", "Pathfinder", "Pfam", "Pidfile", "PolaronPathIntegrals", "Polynomials", "PoreMatMod", "ProbNumDiffEq", "Probably", "ProfileView", "ProxSDP", "Pseudospectra", "QXTns", "QuadraticToBinary", "QuantumTomography", "RandomQuantum", "ReinforcementLearningExperiments", "Relief", "Retriever", "RetroCap", "Ripserer", "RungeKutta", "SDPSymmetryReduction", "SemiseparableMatrices", "SparseIR", "StateSpaceEcon", "StochasticPrograms", "StochasticRounding", "StrBase", "SuiteSparseMatrixCollection", "SunAsAStar", "SymbolicRegression", "SymbolicWedderburn", "Symbolics", "TopoPlots", "TuringGLM", "YaoBase", "YaoBlocks"], 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

dkarrasch commented Jul 31, 2022

This is currently failing only in SparseArrays.jl and should be easily fixable (rename Adjoint{..., <:AbstractQ} to AdjointQ{...,<:AbstractQ}). Otherwise there are only very few packages that require adjustments:

  • ArrayLayouts.jl
  • BandedMatrices.jl: looks like a bit more work
  • GenericLinearAlgebra.jl
  • packages like JOLI.jl and LinearMaps.jl, perhaps also LinearOperators.jl need to accept AbstractQs in their matrix-based operator types, should be straightforward
  • SparseArrays.jl

Otherwise, this has already lead to the detection of "suboptimal" performance usage, see the list of PRs above.

Nanosoldier doesn't run GPU packages, so we may have missed breakage in that part of the ecosystem. @maleadt Would you help me detect packages that would require adjustments, or ping somebody?

Given that breakage appears to be very limited, I think this PR has some chance.

@ViralBShah @KristofferC How does it work? Is SparseArrays.jl still tied to some Julia version, or would earlier Julia versions potentially pick up (breaking) changes like this? I think they wouldn't, because older SparseArrays.jl "versions" are baked into the sysimage, right? If SparseArrays is tied to Julia, controlling the breakage would be quite easy via a VERSION-dependent branching like this

@static if VERSION > v"1.9....DEV..."
    AdjointQtype = AdjointQ
else
    AdjointQtype = Adjoint
end

and then replace the current Adjoint in method signatures by AdjointQtype?

@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.

@andreasnoack
Copy link
Member

This is probably the right thing to do so great that you go through all the work required to make this possible. Part of me thinks it a bit sad that we can't use the matrix abstraction for something like this but, as you point out, AbstractMatrix effectively assumes fast scalar indexing. A natural continuation would be to drop the subtyping for SparseMatrixCSC since it's pretty similar but the costs are not as extreme as for the Qs so the trade-off balance is different.

@dkarrasch
Copy link
Member Author

@nanosoldier runtests(["YaoBlocks", "JOLI", "HarmonicBalance", "StochasticRounding", "ReachabilityAnalysis", "IntervalLinearAlgebra", "Pathfinder", "MultivariateStats", "Pajarito", "NumericalRepresentationTheory", "YaoBase", "DataAssimilationBenchmarks", "SBMLToolkit"])

@nanosoldier
Copy link
Collaborator

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

@dkarrasch dkarrasch changed the title RFC: Don't subtype AbstractQ <: AbstractMatrix Don't subtype AbstractQ <: AbstractMatrix Feb 16, 2023
@dkarrasch
Copy link
Member Author

All green, let's go!

@dkarrasch dkarrasch merged commit d7cbf39 into master Feb 16, 2023
@dkarrasch dkarrasch deleted the dk/q branch February 16, 2023 10:01
@andreasnoack
Copy link
Member

Small comment related to the git history and also the ability to review: IMO it's generally best to split refactoring and functional changes into separate PRs. It's hard to review a diff when functions are also moved to new files.

@dkarrasch
Copy link
Member Author

Thanks for the feedback @andreasnoack. I should have started by first collecting all the existing methods in a new file, and then change the code.

@maleadt
Copy link
Member

maleadt commented Apr 14, 2023

I just noticed that this broke CUDA.jl (and probably other GPU back-ends too). Specifically, our CUSOLVER wrappers testing *(Adjoint{Float32, CuMatrix{Float32}},QRPackedQ{Float32, CuMatrix{Float32}, CuVector{Float32}}), which used to dispatch to

function *(adjA::Adjoint{<:Any,<:StridedVecOrMat}, Q::AbstractQ)
A = adjA.parent
TAQ = promote_type(eltype(A), eltype(Q))
Ac = similar(A, TAQ, (size(A, 2), size(A, 1)))
adjoint!(Ac, A)
return rmul!(Ac, convert(AbstractMatrix{TAQ}, Q))
now end up in
function (*)(A::AbstractMatrix, Q::AbstractQ)
T = promote_type(eltype(A), eltype(Q))
return rmul!(copy_similar(A, T), convert(AbstractQ{T}, Q))
end
where the call to convert loses our GPU array data because of
QRPackedQ{T}(Q::QRPackedQ) where {T} = QRPackedQ(convert(AbstractMatrix{T}, Q.factors), convert(Vector{T}, Q.τ))

Notice the hard-coded Vector, resulting in:

julia> Q
15×15 LinearAlgebra.QRPackedQ{Float32, CuArray{Float32, 2, CUDA.Mem.DeviceBuffer}, CuArray{Float32, 1, CUDA.Mem.DeviceBuffer}}

julia> LinearAlgebra.QRPackedQ{T}(Q)
15×15 LinearAlgebra.QRPackedQ{Float32, CuArray{Float32, 2, CUDA.Mem.DeviceBuffer}, Vector{Float32}}

@dkarrasch
Copy link
Member Author

What was the definition of convert(AbstractMatrix{TAQ}, Q) before? As far as I can tell, the two * methods are the same: for A::Adjoint the copy_similar(A, T) expands to the similar(...); adjoint!(...) code, so we need to catch the convert(AbstractQ, Q) call and make it do the same as the convert(AbstractMatrix, Q) call for Q::QRPackedQ{<:Any,<:CuArray}.

@dkarrasch
Copy link
Member Author

Oh, or would it already help to convert Q.τ to AbstractVector{T}?

@maleadt
Copy link
Member

maleadt commented Apr 14, 2023

What was the definition of convert(AbstractMatrix{TAQ}, Q) before?

It used to fall back to convert(::Type{T}, a::T) where {T<:AbstractArray} = a, so it was a no-op whereas now it does an actual conversion.

I'd also be happy to add a no-op definition for convert(AbstractQ{T}, Q::QRPackedQ{Float32, <:CuArray}), if you think that is the right thing to do.

@maleadt
Copy link
Member

maleadt commented Apr 14, 2023

Yeah, that works:

julia> typeof(convert(LinearAlgebra.AbstractQ{T}, Q))
LinearAlgebra.QRPackedQ{Float32, CuArray{Float32, 2, CUDA.Mem.DeviceBuffer}, Vector{Float32}}

julia> LinearAlgebra.QRPackedQ{T}(Q::LinearAlgebra.QRPackedQ) where {T} = LinearAlgebra.QRPackedQ(convert(AbstractMatrix{T}, Q.factors), convert(AbstractVector{T}, Q.τ))

julia> typeof(convert(LinearAlgebra.AbstractQ{T}, Q))
LinearAlgebra.QRPackedQ{Float32, CuArray{Float32, 2, CUDA.Mem.DeviceBuffer}, CuArray{Float32, 1, CUDA.Mem.DeviceBuffer}}

@dkarrasch
Copy link
Member Author

And for equal eltypes this will be a no-op as well, right? While you're at it, could you please tell me which convert method gets called? Maybe I need to review the dispatch, because I'd want the above method to be used (for equal eltypes).

@dkarrasch
Copy link
Member Author

Nevermind, I found it.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
domain:linear algebra Linear algebra kind:breaking This change will break code kind:speculative Whether the change will be implemented is speculative linalg triage
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

10 participants