RFC: Don't subtype `AbstractQ <: AbstractMatrix`

[dkarrasch]

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, https://github.com/JuliaLang/julia/issues/38972#issuecomment-1066120782, 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]

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

[nanosoldier]

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

[dkarrasch]

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

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

[dkarrasch]

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:

• [x] ArrayLayouts.jl
• [x] BandedMatrices.jl: looks like a bit more work
• [x] GenericLinearAlgebra.jl
• [x] packages like JOLI.jl and LinearMaps.jl, perhaps also LinearOperators.jl need to accept AbstractQs in their matrix-based operator types, should be straightforward
• [x] 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]

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

[andreasnoack]

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]

I think that most of the matrix abstraction still works via the explicitly provided methods, like size, axes, multiplication and division, similar to the Linear[Maps/Operators].jl packages. Unless I'm missing something, this is "only" killing automatic iteration and broadcasting. If there are more matrix aspects that we want to preserve, that's totally fine, and should be included, as long as it can be realized efficiently.

[dkarrasch]

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

[maleadt]

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

[KlausC]

IMO that is a good idea, which solves a minor issue with AbstractArray. I would like to go a step further:

You are right that AbstractQ "is not" an AbstractMatrix, once we had a clear understanding what an AbstractArray is. So I think it would be more important to define that, for example as: "An AbstractArray is the representation of a linear form, which supports indexed access to its canonical elements with getindex and setindex! functions with an effort of O(1) (independent of the dimension and size)." With this definition, neither an AbstractQ nor a SparseMatrixCSC would be subtypes of AbstractArray; nevertheless all of them could be subtypes of a common, supertype, say, AbstractLinearSomething, to be defined.

That could also be used to handle the AbstractArray fallback"*, which has not been treated seriously up to now, if we would also define the "wrapper" containers (like Adjoint, Hermitian, ...) as subtypes of that new AbstractLinearSomething.

*) There are some algorithmes in LinearAlgebra and other packages, which take AbstractMatrix input arguments and silently assume, that there is a fast getindex. That has fatal performance impacts, if the assumption is not true, like for AbstractQ or UpperTriangular{...,Adjoint{...,SparseMatrixCSC}} or similar.

[dkarrasch]

@KlausC Your definition makes lots of sense to me. The role of fast getindex (and perhaps also setindex!) is sometimes downplayed a bit in the community, in the understandable hope to squeeze in FFTs or Linear[Maps/Operators] in the AbstractMatrix hierarchy. The challenge with the new supertype is that I see a few examples of types in the ecosystem that would seem natural candidates, but which likely won't subtype to AbstractLinearSomething, e.g., because they have their own larger type hierarchy including own supertypes. Anyway, do you think your proposal is a showstopper for this PR and should be replaced by some more ambitious action, or should we move forward with the direction taken here, and consider introducing abstract supertypes etc later?

[KlausC]

is a showstopper for this PR ? No, I only wanted to point out, that there is still a way to go for a satisfactory treatment of AbstractArray / AbstractMatrix.

[dlfivefifty]

1. I agree AbstractQ is not an AbstractMatrix since it is sometimes rectangular and sometimes square.
2. SparseMatrixCSC is definitely an AbstractMatrix. I do not agree that one should assume a fast getindex from an abstract type. If you need that feature you should add a trait rather than trying to use an abstract type to determine these behaviour. I.e. we already have IndexLinear and IndexCartesian rather than making new supertypes.
3. If one is going through this effort it is important that there's a crystal clear description of the AbstractQ interface in the docs.

[StefanKarpinski]

What would be nice might be adding a new abstraction called AbstractLinear or something like that and having AbstractMatrix subtype that and then have AbstractQ <: AbstractLinear. However, I'm not sure that we can arrange for an inheritance structure like that. The difference would be that AbstractLinear can be applied as a linear transform but need not have fast scalar indexing.

[dkarrasch]

We could easily introduce another abstract type. I'm just not sure how widely it will be adopted. Just like with AbstractMatrix, which is already subtyping AbstractArray, there are a few other "candidates" that already have their own supertypes (I'm thinking of SciMLLinearOperators [not sure about the naming]). Some others might actually start to subtype, like AbstractFFT, Linear[Maps/Operators] and alike. OTOH, this can be added in another PR, and perhaps requires more discussion: naming, interface, ...

I wonder how we should proceed with this. I have PRs open in the fundamental linalg packages. Their tests pass now and on this PR branch. Shall we merge the downstream PRs, release versions, rerun nanosoldier and see if something is missing? At least, I don't see any objections against this change.

[dkarrasch]

Will we ever do this? Will we include this in v1.9? If some authority gives a thumb up, I can go ahead and release package versions that have an adoption mechanism included (except for GenericLinearAlgebra.jl, for which I don't have permissions; and we would need to bump SparseArrays.jl in the sysimage after the corresponding PR is merged there), and then run nanosoldier again. Or should we first merge this, then adjust the packages and then tidy up related issues? A nanosoldier run doesn't make sense without including "fixes" in the fundamental linalg packages, because many packages will fail for the same reasons in those packages.

[dkarrasch]

@ViralBShah @StefanKarpinski @andreasnoack @dlfivefifty @stevengj Any opinions on this piece? I have functional PRs in BandedMatrices.jl, ArrayLayouts.jl, GenericLinearAlgebra.jl and SparseArrays.jl. If we want this, we have two options:

1. Merge this PR and update the downstream PRs with the specific DEV-version number.
2. Release downstream patches that should work across different versions (including this PR), run pkgeval and see if we missed something, fix things downstream, re-run and reiterate until this doesn't cause any failures.

I'd really love to get this into v1.9 (same with the AdjointFactorization PR), so we can make glue fixes in the beta phase.

[stevengj]

I'm not that convinced that the benefits here are worth the potential breakage…

[dkarrasch]

Thanks for the feedback @stevengj, much appreciated to move towards a conclusion. There are at least two aspects to this PR which are breaking, though not equally severely.

1. Give up AbstractQ <: AbstractMatrix. Since subtypes of AbstractQ are densely covered with their own methods, any intended usage is covered anyway, and what is "inherited" from AbstractMatrix is what we wouldn't want as usage: broadcasting, implicit iteration over elements. All issues arising from this change are presumably fixed by the many small PRs listed in the beginning of this PR.
2. Introduce AdjointQ <: AbstractQ and make adjoint(Q::AbstractQ) = AdjointQ(Q) instead of adjoint(Q::AbstractQ) = Adjoint(Q), because the latter subtypes AbstractMatrix, which we didn't quite like for Q itself. This is breaking because specialized methods with ::Adjoint{<:Any,<:SomeQType} are no longer hit and need "retyping" to ::AdjointQ{<:Any,<:SomeQType}. The relevant changes for the registered packages affected by this change are already in the pipeline, but I understand that there may be heavy unregistered linear algebra-related packages that are not covered by PkgEval.

Would you support keeping 1. but abandoning 2. for the time being (until v2)?

[stevengj]

I'm guess I'm not clear on what we're gaining, besides avoiding accidentally using Q in methods where it is slow, that is worth breakage here.

[dkarrasch]

Here are some issues with our current design:

• Accidentally using Q in methods where it is slow, as you say. If you look at julia v1.6, then products of structured matrices (that have "their own" fast multiplication methods) with Q's are all slow. The reason is that Q's also have their own specialized multiplication code, and it requires tons of methods to disambiguate the situation, and in this situation it means to always densify the structured matrix and then to use the Q code. On my machine, the benchmark code from https://github.com/JuliaLang/julia/pull/44615#issuecomment-1068080855 yields the following results:

133.588 μs (3 allocations: 705.70 KiB)
  140.725 μs (3 allocations: 705.70 KiB)
  309.615 ms (270002 allocations: 448.38 MiB)
  316.566 ms (270004 allocations: 448.38 MiB)
  311.354 ms (270003 allocations: 448.38 MiB)
  984.662 ms (808236 allocations: 1.31 GiB)
  312.473 ms (270002 allocations: 448.38 MiB)
  958.825 ms (808235 allocations: 1.31 GiB)
  309.983 ms (270002 allocations: 448.38 MiB)
  948.879 ms (808235 allocations: 1.31 GiB)
  132.578 μs (3 allocations: 705.70 KiB)
  140.487 μs (3 allocations: 705.70 KiB)
  144.079 μs (3 allocations: 705.70 KiB)
  304.931 ms (270004 allocations: 448.38 MiB)
  310.116 ms (270003 allocations: 448.38 MiB)
  964.747 ms (808236 allocations: 1.31 GiB)
  308.169 ms (270002 allocations: 448.38 MiB)
  967.538 ms (808235 allocations: 1.31 GiB)
  302.516 ms (270002 allocations: 448.38 MiB)
  954.196 ms (808235 allocations: 1.31 GiB)

Missing methods yield a performance drop of more than three orders of magnitude (and I believe that the overhead scales squared with the dimension, because it needs to compute n^2 elements in matrixification). Any package that defines custom matrices that have preferential multiplication code when multiplied with "generic"/"abstract" matrices need to be taught that they need to step back when they meet a Q. The situation is particularly bad when, as in v1.6+ up until a few days ago's nightly, we have on the one hand

*(A::StridedMatrix{T} where T, Q::LinearAlgebra.AbstractQ)
     @ LinearAlgebra /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/qr.jl:775

and on the other hand

*(D::Diagonal, A::AbstractMatrix) in LinearAlgebra at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/diagonal.jl:249

Without explicit disambiguation (for every single special matrix type out there), this is a performance nightmare.

• Even after breaking the subtyping AbstractQ <: AbstractMatrix, the exact same story applies to Adjoint(AbstractQ).
• If a package defines a new AbstractQ type, it is at risk that either Q or Q' pose performance challenges. For instance, this PR here deletes the following method:

*(D::Diagonal, adjQ::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) =
    rmul!(Array{promote_type(eltype(D), eltype(adjQ))}(D), adjQ)

But what is happening to adjoints of other Q-types? The disambiguation of special matrix types and Q's takes up more than fifty lines of code alone in special.jl, not counting the quoted D*Q' method above and perhaps a few more.

• Currently, AbstractQ "pretends to be" an abstract supertype, but it's not used as such. It doesn't provide an interface, and it barely provides methods that new concrete subtypes would fall back to. Thus, you can find that every package with its custom AbstractQ subtype defines the adjoint method, even within the LinearAlgebra stdlib. In fact, the Q from LQ didn't subtype to AbstractQ, and initially I was confused why. The reason is that AbstractQ meant to summarize types that are flexible in size when multiplied from the left or their adjoint from the right. It was more like a type union for those Q-types, excluding others.

So, the benefits of this PR are, in particular:

1. a clear interface/abstract type with a variety of fallbacks that allow for convenient subtyping and minimal effort to provide a fully functional AbstractQ subtype;
2. a clear distinction between "abstract matrices" and Q's (and their adjoints!), where Qs always have precedence in multiplication, no matter how sparse/structured the matrix is, thus avoiding the need for disambiguation between Q's and special/custom matrices.

Addendum: currently, you can write Q\A, and now guess which method is called!

julia> @which Q\rand(3,3)
\(A::AbstractMatrix, B::AbstractVecOrMat)
     @ LinearAlgebra /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/generic.jl:1096

Ouch! With the current state of this PR, this call would at least yield a method error, but I think we should add methods that translate solving into multiplication by the adjoint. We have that encoded in inv already, and ArrayLayouts.jl has it as well, so I think that makes sense for AbstractQ in the stdlib, too.

[stevengj]

LGTM.

[dkarrasch]

Since it's very early in the v1.10 release cycle, we have plenty of time, so I'll release patches for the fundamental linalg packages that would pick up the changes and should therefore work on both current master and this PR's branch. Then I'll run nanosoldier and fix all issues before merging this.

[dlfivefifty]

What would be nice might be adding a new abstraction called AbstractLinear

It would then also make sense to have UniformScaling <: AbstractLinear. There is a nice similarity that both UniformScaling and AbstractQ can act on vectors of different sizes.

[dkarrasch]

@nanosoldier runtests()

[dlfivefifty]

What about views? I.e. will this no longer work:

julia> view(qr(randn(5,5)).Q, :, 1:3)
5×3 view(::LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}, :, 1:3) with eltype Float64:
 -0.297051  -0.673535   0.544172
  0.552457   0.283495   0.75936
  0.74138   -0.590937  -0.316325
  0.113867   0.177987  -0.163323
 -0.209624  -0.291712   0.0226731

[dkarrasch]

Currently probably not. But what should a view on a function-based AbstractLinear object be? Probably any operation downstream will treat it like an AbstractArray and start to get the entries one by one, which will allocate a huge amount of memory even though you thought everything should be allocation-free, since it's a view. I think this is the educational part of this endeavour here: don't think of Qs as of arrays.

[nanosoldier]

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

[dkarrasch]

From the nanosoldier run, here is another common pattern of failure with this PR:

MethodError: no method matching lmul!(::LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}, ::SparseArrays.SparseMatrixCSC{Float64, Int64})

ERROR: LoadError: MethodError: no method matching lmul!(::LinearAlgebra.AdjointQ{Float64, LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}}, ::SparseArrays.SparseMatrixCSC{Float64, Int64})

On v1.8 and v1.9, the calls Q * S and Q'S, S sparse, go through generic_matmatmul!, which means that each entry of the matrix representation gets computed during the matrix-multiplication loop, and the (typically) dense result gets written into a sparse matrix. 😱 Perhaps we should override the fallback behaviour and make S dense and then multiply into it.