Missing Ac_mul_B(StridedMatrix,SparseMatrixCSC) causes dispatch to generic matmul

[mes9]

I'm running the following code which takes about 200 seconds on Windows and Linux, while only .3 seconds on a Mac

srand(1)
X=randn(20000,30)
W=spdiagm(rand(20000))
@time res1=X'W*X

My versioninfo():

Julia Version 0.6.0
Commit 903644385b* (2017-06-19 13:05 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i5-6200U CPU @ 2.30GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, skylake)

[andreasnoack]

I don't have a Windows machine handy so could you try timing @time res1=(X'W)*X and @time res1=X'*(W*X) as well?

[mes9]

@time res1=(X'W)*X is slow, taking about the same amount of time, but @time res1=X'*(W*X) is fast! taking less than .01 seconds.

Is there any rule of thumb for how I should be grouping the matrices or can it be optimized in the code?

[mauro3]

On Linux with 0.6 I also get a slow @time res1=X'W*X.
The expressions get expanded to different things:

julia> expand(:(X'W*X))                                                                                                                
:(Ac_mul_B(X, W) * X)                                                                                                                  

julia> expand(:((X'*W)*X))                                                                                                             
:(Ac_mul_B(X, W) * X)                                                                                                                  

julia> expand(:(X'*(W*X)))                                                                                                             
:(Ac_mul_B(X, W * X))                                                                                                                  

julia> expand(:(X'*W*X))                                                                                                               
:(ctranspose(X) * W * X)                                                                                                               

What do those show on the Mac?

[andreasnoack]

It is also slow on my Mac. @mes9 Are you sure the expressions are exactly the same on Windows and Linux. I get that X'*W*X is fast and X'W*X is slow so the problem seems to be that Ac_mul_B(::Matrix,::SparseMatrixCSC) hits our generic fallback via https://github.com/JuliaLang/julia/blob/dcf39a1ddaa788feb231f06611f849c18dcece76/base/linalg/matmul.jl#L218 whereas X'*W*X=(*)(ctranspose(X),W,X)=(ctranspose(X)*W)*W and (*)(::Matrix,::SparseMatrixCSC) has an optimized version for sparse matrices. The solution is to define Ac_mul_B(::StridedMatrix,::SparseMatrixCSC) in the sparse matrix code.

[mes9]

yes, if I include the * between X' and W it is fast, otherwise it is slow. I agree, the issue just appears to be sparse matrices. When I make W just a regular diagonal matrix it takes about a second whether or not the * is included

[andreasnoack]

This will be fixed by JuliaLang/julia#24045

[dkarrasch]

This has been resolved, I think. All timings are equal and far below 1 second.