A_mul_B and all that jazz

[jiahao]

The purpose of this issue is to discuss the future of the dozens of functions in Julia that are slight variations of matrix multiplication, matrix division and backslash. I presume that the point was to have fine-grained control over the use of matrix variables and temporaries when needed. However, I don't understand why the availability of these functions is so inconsistent: is this intentional design or a consequence of organic growth with methods being implemented as needed?

List of all the functions (existing and plausibly existing in the future)

Here are the 3x3x3x2=54 possible functions that can be constructed out of A, its transpose At and its Hermitian conjugate Ac; and similarly for B, Bt and Bc; mul, rdiv and ldiv; and mutating and non-mutating variants. (### means that this function does not exist.)

A_mul_B #deprecated synonym of A*B
A_mul_B!
A_mul_Bt
A_mul_Bt!
A_mul_Bc
A_mul_Bc!

At_mul_B
At_mul_B!
At_mul_Bt
At_mul_Bt!
At_mul_Bc   ###
At_mul_Bc!  ###

Ac_mul_B
Ac_mul_B!
Ac_mul_Bt #deprecated
Ac_mul_Bt!  ###
Ac_mul_Bc
Ac_mul_Bc!


A_ldiv_B    ### synonym of A\B
A_ldiv_B!
A_ldiv_Bt
A_ldiv_Bt!  ###
A_ldiv_Bc
A_ldiv_Bc!  ###

At_ldiv_B
At_ldiv_B!  ###
At_ldiv_Bt
At_ldiv_Bt! ###
At_ldiv_Bc  ###
At_ldiv_Bc! ###

Ac_ldiv_B
Ac_ldiv_B!  ###
Ac_ldiv_Bt  ###
Ac_ldiv_Bt! ###
Ac_ldiv_Bc
Ac_ldiv_Bc! ###


A_rdiv_B    ### synonym of A/B
A_rdiv_B!   ###
A_rdiv_Bt
A_rdiv_Bt!  ###
A_rdiv_Bc
A_rdiv_Bc!  ###

At_rdiv_B
At_rdiv_B!  ###
At_rdiv_Bt
At_rdiv_Bt! ###
At_rdiv_Bc  ###
At_rdiv_Bc! ###

Ac_rdiv_B
Ac_rdiv_B!  ###
Ac_rdiv_Bt  ###
Ac_rdiv_Bt! ###
Ac_rdiv_Bc
Ac_rdiv_Bc! ###

I'm rocking this boat because it's come up as a natural consequence of systematically making methods available for special matrix types.

Specific issues

1. Documentation. The existence, motivation and use of these functions is sparsely documented. Not all the functions even have defined help strings.

2 Growth and proliferation. Do we foresee an eventual future where all 54-3=51 functions (omitting the synonyms for *, / and \) become defined and have methods implemented? There would a great many methods that would have to be implemented to cover all the possible types of matrices that currently exist.

• Can we eliminate any of these functions on the basis of being useless? There are no functions defined that mix Ac
  and Bt, for example; however, this is a somewhat rare but entirely plausible combination to operate on. Do we care?
• How many of these functions can be implemented with just a fallback method?
• Having more things in Base without reason is not ideal, but conversely these functions form a coherent grouping for which it doesn't make sense to have some but not others.
• How many of these functions can be defined trivially in terms of lazy array views (#5003) and rewriting rules (Label QRCompactWY.T as Triangular #122)?

cc: @JeffBezanson @StefanKarpinski @andreasnoackjensen @dmbates @lindahua

[johnmyleswhite]

If we get rid of these functions, am I right to assume we'll be pushing information into the type system using transpose types and conjugate types?

[jiahao]

I think that would depend quite strongly on how #42 is resolved.

[dmbates]

Because I develop iterative methods that require a lot of numerical linear algebra within each iteration, I want to be careful about reusing existing storage when possible. Hence I find these mutating functions to be worthwhile. I could always write them as direct calls to functions in the BLAS or LAPACK modules but this makes the code harder to read for those who didn't grow up needing to know Lapack calling sequences as I did.

However, I would be amenable to hiding all these functions in the Base namespace and only exporting the non-mutating versions that are called when expressions like C = A'*B are parsed. The use of mutating functions to save memory allocations is a "consenting adults" activity in my opinion. If you are going to engage in it you should be aware of the potential problems.

[ViralBShah]

I think only export non-mutating versions is the right choice. What about the transposed multiply operators? Should those also not be exported by default?

I would be ok with removing the c and t combinations too. It wasn't fun dealing with all of those in a recent PR.

For the stuff with lazy array views, we should wait until we get there, or else we will see a performance degradation. At that point, I would expect that we can get rid of the transposed multiply and divide operators.

[johnmyleswhite]

To me, it seems like it would be better to wait until we have a general tool for expressing mutating operations before removing these functions.

[ViralBShah]

I am only giving a +1 to the proposal to remove the Ac_mul_Bt where you are mixing complex and real transposed matrices, and not exporting all these in Base by default.

[johnmyleswhite]

Ok. That makes sense to me, although I'd think there's something to be said for keeping an eyesore around as a reminder that we need a cleaner solution.

[ViralBShah]

There is some value to eyesore, but keeping it around longer means more people use it thinking it is a supported API and makes it difficult to remove in the future.

[johnmyleswhite]

Very true.

[jiahao]

Time for an @Embarrassing macro?

[andreasnoack]

Yesterday I found out that

julia> n=1000;
julia> A=randn(n,n);A=A+A';Aeig=eigfact(Symmetric(A));
julia> V=Aeig[:vectors];D=diagm(Aeig[:values]);
julia> @time V*D*V';
elapsed time: 0.263635665 seconds (24000512 bytes allocated)
julia> @time V*(D*V');
elapsed time: 0.235679586 seconds (16000272 bytes allocated)
julia> @time (V*D)*V';
elapsed time: 0.237066579 seconds (16000272 bytes allocated)

because

julia> @which V*D*V'
*(a,b,c) at operators.jl:71

causing the transpose to be calculated before the multiplication.

I think that @johnmyleswhite's proposal sounds like the right one. It could eliminate all the A_x_B functions and only require mul!,*,ldiv!,rdiv!,\ and /. The rest is dispatch. I think all the function bodies would be unaffected. Should we try this change? I can do the coding, but I would like to know it could fly before spending the time.

[johnmyleswhite]

I believe the problem with the transpose type (which somebody else proposed before me and I just brought up again) is that we're not sure how to implement it for higher-order tensors.

[Jutho]

Another suggestion would be to follow more closely the BLAS approach and have a single function with additional arguments 'N', 'T' and 'C' specifying the multiplication pattern. This of course only shifts the problem to branching within the method, and probably calling all of the above functions anyway, but the advantage is that only a single function needs to be exported and documented.

Occasionally I regret that the mutating multiplication functions in Julia do not more closely resemble the BLAS design, i.e. not only allow to store the result of the multiplication in a given array, but actually add the result to it. I guess there has been some thought in the BLAS design, and it would be great if Julia could generalise this functionality (i.e. matrices with no stride 1, different types of dense and sparse matrices, ...) without sacrificing the existing functionality.

I know that gemm and gemv are available within Base.LinAlg.BLAS , but if you define a new type hierarchy, say AbstractLinearOperator, for sparse linear operators for use in an iterative linear solver or eigensolver, then you want to implement a method A_mul_B! for its concrete subtypes. If you would allow to actually add the result of the matrix vector multiplication to a given vector in this method, then you could also define a SumOfLinearOperators <: AbstractLinearOperator type, define + for any AbstractLinearOperator to return a SumOfLinearOperators which just keeps track of the different AbstractLinearOperator terms, and implements A_mul_B! for an input vector by calling A_mul_B! sequentially on all the terms, thereby adding the result and not having to create any copies.