RFC: Performance improvement for Diagonal' * Vector|Matrix

[georgemarrows]

Added tests
Only addresses diagonal part of #21286, not sparse matrices

julia> include("test.jl")
  0.000014 seconds (6 allocations: 78.359 KiB)
  0.000012 seconds (6 allocations: 78.359 KiB)
  0.000012 seconds (6 allocations: 78.359 KiB)

from

begin
	D = Diagonal(randn(10000))
	v = randn(10000)
	@time D * v
	@time D' * v
	@time D.' * v
end

[georgemarrows]

First Julia PR - please be gentle :-)
Potential problem, so far as I can see:

• tests A'*B, not Ac/t_mul_B directly
• no performance tests
• lack of generality? - Diagonal' * Vector only

[andreasnoack]

Our (c)transpose is recursive so it should be ctranspose.(A.diag) .* B. It would be great if you could make the same fix for Ac/t_mul_B!.

[georgemarrows]

Ah, for block diagonal matrices? New PR tomorrow covering both points.

[andreasnoack]

Yes. You can find a couple of issues where the details have been discussed. Please just update this PR instead of opening a new one.

[georgemarrows]

Fixed for block matrices and a new test added.

For Ac/t_mul_B!, I notice that A_mul_B! is actually even slower:

begin
	D = Diagonal(randn(10000))
	v = randn(10000)
	vv = similar(v)
	@time A_mul_B!(vv, D, v)
	@time Ac_mul_B!(vv, D, v)
	@time At_mul_B!(vv, D, v)

	@time vv .= D * v
	@time vv .= Ac_mul_B(D, v)
	@time vv .= At_mul_B(D, v)
	nothing
end

 1.570805 seconds (4 allocations: 160 bytes)
  1.095022 seconds (4 allocations: 160 bytes)
  1.093878 seconds (4 allocations: 160 bytes)
  0.000030 seconds (6 allocations: 78.359 KiB)
  0.000019 seconds (6 allocations: 78.359 KiB)
  0.000017 seconds (6 allocations: 78.359 KiB)

Should l do a change for that too?
Are the vv .= Ac_mul_B(D, v) style implementations ok? They're way faster but allocate more.
Should block diagonals work here too?

BTW should the following work? I wanted to use it in writing my test for transposes.

julia> D = Diagonal([[1 2; 3 4], [1 2; 3 4]])
2×2 Diagonal{Array{Int64,2}}:
 [1 2; 3 4]      ⋅     
     ⋅       [1 2; 3 4]

julia> full(D)
ERROR: MethodError: no method matching zero(::Type{Array{Int64,2}})
Closest candidates are:
  zero(::Type{Base.LibGit2.GitHash}) at libgit2/oid.jl:106
  zero(::Type{Base.Pkg.Resolve.VersionWeights.VWPreBuildItem}) at pkg/resolve/versionweight.jl:80
  zero(::Type{Base.Pkg.Resolve.VersionWeights.VWPreBuild}) at pkg/resolve/versionweight.jl:120
  ...
Stacktrace:
 [1] diagm(::Array{Array{Int64,2},1}, ::Int64) at ./linalg/dense.jl:251
 [2] full(::Diagonal{Array{Int64,2}}) at ./linalg/diagonal.jl:56

[georgemarrows]

Now with perf fixes for Ac/t_mul_B for Vector and Diagonal. All tests pass. Will push again with squashed commits if you're happy with this.

[andreasnoack]

Couldn't this just be vout .= A.diag .* vin and thereby avoiding the allocation?

[georgemarrows]

Indeed, and then you don't need to have separate non-! versions. I tried that briefly before but missed that I was measuring compilation time in there too :-(

Notes on performance here: https://gist.github.com/georgemarrows/678876779e7292a954a344aa3addeaf1

[andreasnoack]

Did you confirm that this version still fixes #21286? I'm not sure that D::Diagonal'*B::Matrix dispatches to the two argument Ac_mul_B!s in the current version.

[georgemarrows]

Latest push covers D' * Matrix as well as D' * Vector, and adds tests for this combination too. It doesn't address sparse matrices, which were also raised in #21286 - I think that's better as a separate PR (and perhaps a separate bug report).

I ran tests/linalg/diagonal.jl, and this gist shows that performance is hugely increased in both Vector and Matrix cases so that transpose and ctranspose cases are now similar to vanilla multiplication.

[georgemarrows]

Nudging for review when time is available...

[Sacha0]

LGTM from a brief review! Thanks for the great first PR @georgemarrows! :)

[georgemarrows]

Thanks @Sacha0. I've rebased and squashed commits.
make test all good locally, and we still have the performance improvements in this gist.
PR title changed to "RFC:"

[Sacha0]

Additional thoughts @andreasnoack? Mergeworthy? Best!

[tkelman]

we really shouldn't be syntax-highlighting in...

[georgemarrows]

Thanks @andreasnoack, I appreciate your guidance and reviews while I thrashed towards a solution.

[andreasnoack]

You are welcome. Thanks for the contribution.