slow matrix multiplication for transpose(real) * complex

[jishnub]

Cross-posting from this discourse thread that didn't generate a lot of discussion:

On v1.7 (and also nightly)

julia> ZC = rand(1800, 1680); M = rand(ComplexF64, 1800, 1800);

julia> @time ZC' * M;
  9.169792 seconds (3 allocations: 46.143 MiB)

julia> @time permutedims(ZC) * M;
  0.877853 seconds (6 allocations: 115.357 MiB, 1.55% gc time)

julia> @which ZC' * M
*(A::AbstractMatrix, B::AbstractMatrix) in LinearAlgebra at /home/jishnu/Downloads/julia/julia-1.7.0/share/julia/stdlib/v1.7/LinearAlgebra/src/matmul.jl:151

julia> @which permutedims(ZC) * M
*(A::StridedMatrix{var"#s859"} where var"#s859"<:Union{Float32, Float64}, B::StridedMatrix{var"#s858"} where var"#s858"<:Union{Float32, Float64, ComplexF32, ComplexF64}) in LinearAlgebra at /home/jishnu/Downloads/julia/julia-1.7.0/share/julia/stdlib/v1.7/LinearAlgebra/src/matmul.jl:158

This performance difference doesn't seem to exist for real matrices:

julia> ZC = rand(1800, 1680); M = rand(1800, 1800);

julia> @time ZC' * M;
  0.261731 seconds (3 allocations: 23.071 MiB, 2.57% gc time)

julia> @time permutedims(ZC) * M;
  0.278280 seconds (4 allocations: 46.143 MiB)

This difference seems to arise because ZC' * M for a complex M calls

@inline function mul!(C::AbstractMatrix, adjA::Adjoint{<:Any,<:AbstractVecOrMat}, B::AbstractVecOrMat,
                 alpha::Number, beta::Number)
    A = adjA.parent
    return generic_matmatmul!(C, 'C', 'N', A, B, MulAddMul(alpha, beta))
end

and generic_matmatmul! is presumably slow. For a real M, it calls

@inline function mul!(C::StridedMatrix{T}, adjA::Adjoint{<:Any,<:StridedVecOrMat{T}}, B::StridedVecOrMat{T},
                 alpha::Real, beta::Real) where {T<:BlasReal}
    A = adjA.parent
    return mul!(C, transpose(A), B, alpha, beta)
end

Would it be possible to add a similar method for complex matrices to avoid falling back to the generic method?

[stevengj]

We call BLAS for complex A' * B in mul!: https://github.com/JuliaLang/julia/blob/6893f21a896e51989c40b0b6de5cb2512b8afc97/stdlib/LinearAlgebra/src/matmul.jl#L410-L418
which is called by the 3-argument mul!: https://github.com/JuliaLang/julia/blob/6893f21a896e51989c40b0b6de5cb2512b8afc97/stdlib/LinearAlgebra/src/matmul.jl#L254-L256

[stevengj]

Oh, I see, the issue is that you have a complex ZC' multiplied by a real M. The reason why this is slower is that that BLAS does not provide an accelerated matmul for this case.

In contrast, without the conjugate-transpose (note that ZC' is the conjugate transpose), there is a trick to exploit BLAS for complex * real matrix: https://github.com/JuliaLang/julia/blob/6893f21a896e51989c40b0b6de5cb2512b8afc97/stdlib/LinearAlgebra/src/matmul.jl#L151-L163
The only way to employ this method for ZC' * M would be to materialize the adjoint via Matrix(ZC') * M or similar.

(Note that permutedims(ZC) * M is not equivalent to ZC' * M because it lacks the conjugation. )

[stevengj]

I don't really see any way of fixing this issue short of writing our own BLAS (e.g. using https://github.com/JuliaLinearAlgebra/Octavian.jl or similar).

[mcabbott]

Another possibility might be to copy, as is done for mismatched element types:

julia> A = rand(2000, 2000); A32 = Float32.(A); C = similar(A);

julia> @time A * A;
  0.079667 seconds (2 allocations: 30.518 MiB)

julia> @time A * A32;  # twice the memory
  0.072415 seconds (4 allocations: 61.035 MiB, 7.09% gc time)

julia> @time mul!(C, A, A32);
  3.763957 seconds (3 allocations: 25.688 KiB)

[stevengj]

Yes, we could create a copy for *, but we probably wouldn't want to do that for mul!.

[jishnub]

Actually in the OP it's transpose(real) * complex, where the conjugation in the adjoint is somewhat superfluous (I merely used ' as the notation is convenient). In this case the transpose is equivalent to permutedims. Creating a copy seems reasonable if there's no alternative.

Edit: Although for real ZC and complex M, we obtain ZC' * M = (M' * ZC)', so the operation adjoint(real) * complex is equivalent to adjoint(adjoint(complex) * real), although I'm not sure if this is how it is evaluated internally.

[stevengj]

Right, I mis-read it, it's adjoint(real) * complex == transpose(real) * complex. That variant also seems incompatible with BLAS, unfortunately.