lazier, less-jazzy linalg internals

[Sacha0]

This pull request completes the first two major blocks of work towards JuliaLang/LinearAlgebra.jl#57. Specifically,

1. The first commit introduces Adjoint and Transpose types, basic functionality for those types, and an extensive test suite for that functionality.

2. The remaining twenty-one commits then strip all A[ct]_(mul|ldiv|rdiv)_B[ct][!] definitions from Base. More precisely, these commits transform each such definition to a corresponding *///\/mul!/ldiv!/rdiv! definition over Adjoint/Transpose, and then place the original signature in base/deprecated.jl as a shim that calls the transformed definition.

   (In a vain attempt to minimize erosion of my sanity in this process, I made the preceding transformations in the most safe and mechanically simple manner I could: Each method like A_ldiv_Bc!(A::TypeA, B::TypeB) = $METHBODY tranforms in place to ldiv!(A::TypeA, adjB::Adjoint{<:Any,<:TypeB}) = (B = adjB.parent; $METHBODY) with shim A_ldiv_Bc!(A::TypeA, B::TypeB) = ldiv!(A, Adjoint(B)) placed in base/deprecated.jl. I also de-meta-fied definitions where that made life simpler.)

All intermediate commits pass the linear algebra and sparse test suites locally.

Collapsing all e.g. A[ct]_mul_B[ct] methods into a single e.g. * function yields the mother of all ambiguity resolution nightmares. So while I fixed all ambiguities necessary to get each commit to pass tests, I expect CI's general method ambiguity filter to detonate on contact. Ambiguities resolved.

(Updated.) Next steps, whether in this or subsequent pull requests:

3. Strip all A[ct]_(mul|ldiv|rdiv)_B[ct][!] definitions from stdlib/. done here
4. Strip all explicit A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls from base/. done here
5. Strip all explicit A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls from test/. done here
6. Strip all explicit A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls from stdlib/. done here
7. Transition linalg to Adjoint/Transpose externally. done in transition linalg to Adjoint/Transpose externally #25083
8. Strip all .' instances from base/. done in remove special lowering for and deprecate .' #25125
9. Strip all .' instances from test/. done in remove special lowering for and deprecate .' #25125
10. Strip all .' instances from stdlib/. done in remove special lowering for and deprecate .' #25125
11. Remove special A[t]_(mul|ldiv|rdiv)_B[t] lowering. done in remove special lowering for and deprecate .' #25125
12. Deprecate .'. done in remove special lowering for and deprecate .' #25125
13. Rewrite all semantically eager ' calls in base/. done in update context-independent and remove context-dependent lowering of ' #25148
14. Rewrite all semantically eager ' calls in test/. done in update context-independent and remove context-dependent lowering of ' #25148
15. Rewrite all semantically eager ' calls in stdlib/. done in update context-independent and remove context-dependent lowering of ' #25148
16. Change ' lowering from adjoint to Adjoint. done in update context-independent and remove context-dependent lowering of ' #25148
17. Remove special A[c]_(mul|ldiv|rdiv)_B[c] lowering. done in update context-independent and remove context-dependent lowering of ' #25148
18. Deprecate all the things. done in sunset linalg jazz #25217
19. Clean up, docs, news. partial sunset linalg jazz #25217
20. Sleep.

Best!

[Sacha0]

(CI appears happy apart from the global method ambiguity check and stdlib/suitesparse, which I had forgotten.)

[Sacha0]

Fixed the one non-ambiguity CI failure. Satisfying the ambiguity check...

julia> length(Test.detect_ambiguities(Core, Base; imported=true, recursive=true, ambiguous_bottom=false))
387

... might take a few days.

[Sacha0]

Unsurprisingly, the vast majority of the ambiguities come from transferring the broad fallbacks for A[ct]_(mul|ldiv|rdiv)_B[ct](a, b) from base/operators.jl (for example Ac_mul_B(a, b) = adjoint(a) * b) to the corresponding consolidation functions; the number of ambiguities falls from 300+ to 99 on dropping the relevant commit (cbea711). The resolution signatures for the 99 ambiguities without that commit seem mostly reasonable, e.g. ambiguous pair

(*(A::AbstractArray{T,2} where T, adjB::Base.LinAlg.Adjoint{#s336,#s335} where #s335<:(Union{Hermitian{T,S}, Hermitian{Complex{T},S}, Symmetric{T,S}} where S where T<:Real) where #s336) in Base.LinAlg at linalg/symmetric.jl:332, *(adjA::Base.LinAlg.Adjoint{#s336,#s335} where #s335<:Base.LinAlg.AbstractTriangular where #s336, B::AbstractArray{T,2} where T) in Base.LinAlg at linalg/triangular.jl:1700)

requiring resolution signature

Tuple{typeof(*),Base.LinAlg.Adjoint{#s336,#s335} where #s335<:Base.LinAlg.AbstractTriangular where #s336,Base.LinAlg.Adjoint{#s336,#s335} where #s335<:(Union{Hermitian{T,S}, Hermitian{Complex{T},S}, Symmetric{T,S}} where S where T<:Real) where #s336}

But the resolution signatures for the 300+ ambiguities with that commit get a bit abstract, e.g. ambiguous pair

(*(a::Adjoint, b::Adjoint) in Base.LinAlg at linalg/adjtrans.jl:105, *(adjA::Adjoint{#s336,#s335} where #s335<:Base.LinAlg.AbstractTriangular where #s336, B::AbstractArray{T,2} where T) in Base.LinAlg at linalg/triangular.jl:1700)

requiring resolution signature

Tuple{typeof(*),Adjoint{#s336,#s335} where #s335<:Base.LinAlg.AbstractTriangular where #s336,Adjoint}

I begin to lean towards simply not transferring those (overly broad?) fallbacks. Thoughts? Thanks!

[Sacha0]

Update, and a fundamental design point (opt-in versus opt-out laziness):

Transferring the broad fallbacks for A[ct]_(mul|ldiv|rdiv)_B[ct][!] from base/operators.jl to the corresponding *///\/mul!/ldiv!/rdiv! functions seems practically untenable due to the resulting ambiguities: Resolving one such ambiguity often introduces O(1)-O(10) more, so resolving both the initial 300-some and all downstream ambiguities likely would require substantially more than 300-some methods. Moreover, the resulting collection of methods would be brittle / difficult to work with at best.

So instead I added the shims for those broad fallbacks, e.g. Ac_mul_B(a, b) = *(Adjoint(a), b), without transferring the associated methods to the corresponding functions. This changed dispatch, causing O(100) errors (mostly MethodErrors, many differing only in operand eltypes). Fortunately those errors were relatively straightforward to fix, and the present version of this branch passes the linalg, sparse, and stdlib test suites locally. This leaves about 100 (relatively reasonable) ambiguities to fix.

What are the implications of not transferring those fallbacks?

Formerly, due to those broad fallbacks, most A[ct]_(mul|ldiv|rdiv)_B[ct] calls that did not directly hit appropriate A[ct]_(mul|ldiv|rdiv)_B[ct] specializations would eagerly adjoint/transpose the operands and then hit a decent specialization (instead of directly hitting a relatively slow generic AbstractArray method). The chief advantage of transferring those broad fallbacks was retaining that behavior. But there's another way to achieve similar behavior in practice:

This pull request presently makes lazy Adjoint/Transpose opt-out by type, in that if one does not want to define specializations over Adjoint{...,MyType}/Transpose{...,MyType}, one must define Adjoint(A::MyType) = adjoint(A) and likewise for Transpose. But laziness could instead be opt-in by type, and after some thought that might be much better. That is, rather than having the generic Adjoint(A) method construct an Adjoint-wrapped A::MyType, instead have the eager fallback Adjoint(A) = adjoint(A), and if one wants to opt-in to laziness, one must define Adjoint(A::MyType) = _Adjoint(A) (where the generic _Adjoint(...) method creates an Adjoint-wrapped A::MyType) and specializations over Adjoint{...,MyType}. Then (1) behavior like the nice former A[ct]_(mul|ldiv|rdiv)_B[ct] fallback behavior described above happens mostly by default, (2) we achieve that behavior while avoiding the combinatorial ambiguity explosion that transferring the broad fallbacks involves, (3) the eager->lazy adjoint/transpose transition becomes less breaking overall and easier for users, and (4) after the transition, users need pay attention to lazy vs. eager adjoint/transpose primarily only for optimization (rather than of necessity).

Given those upsides and the problems associated with transferring the broad fallbacks, I presently lean heavily towards not transferring the broad fallbacks and making laziness opt-in. By default I will follow that approach in this and subsequent PRs and flip the opt-in switch at some point in that process. Best!

[Sacha0]

This pull request now replaces/rewrites all A[ct]_(mul|ldiv|rdiv)_B[ct][!] definitions and (explicit) calls in base/, test/, and stdlib/, completing items 1-6 in the OP's updated task list.

@andreasnoack (thanks again!) and I discussed how to best move forward earlier today. The OP's task list has been updated accordingly, and tomorrow morning PT I plan to give this pull request another pass, make any necessary touchups, merge, and then approach tasks 7+ in a followup pull request. Best!

[vchuravy]

Why does this need to be @pure

[Sacha0]

This definition was preexisting for RowVectors :).

[vtjnash]

This is very strongly not a pure function and terrible things will happen.

[Sacha0]

As above, this definition was preexisting for RowVectors :).

[Jutho]

Thanks this is great. However, I thought this PR would come together with the parser changes that no longer parse A'*B as Ac_mul_B. Am I mistaken that it cannot be the plan that the following no longer works:

julia> struct MyType{T<:Number}
           a::T
       end

julia> Base.:*(x::MyType,y::MyType) = MyType(x.a*y.a)

julia> Base.adjoint(x::MyType) = x

julia> x=MyType(3)
MyType{Int64}(3)

julia> x'*x
ERROR: MethodError: no method matching *(::Base.LinAlg.Adjoint{Any,MyType{Int64}}, ::MyType{Int64})
Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...) at operators.jl:469
  *(::AbstractArray, ::Number) at arraymath.jl:55
  *(::AbstractArray{T,2} where T, ::Base.LinAlg.Transpose{#s344,#s343} where #s343<:RowVector where #s344) at linalg/rowvector.jl:221
  ...
Stacktrace:
 [1] Ac_mul_B(::MyType{Int64}, ::MyType{Int64}) at ./deprecated.jl:2968
 [2] top-level scope

Could we please have adjoint and tranpose in the fallback methods? Or are custom types supposed to overload Adjoint and Transpose? It feels strange to overload type constructors.

Or is that what #25083 is about?

[Sacha0]

Or is that what #25083 is about?

Yes, #25083 and downstream PRs continue this work :).

[Jutho]

Ok thanks. I don't see the parser changes there either yet, but I guess those are to come then? Thanks again for all the work!

[Sacha0]

Please see the OP for a task list. Best!