Fix mapreduce on AdjOrTrans

[SobhanMP]

see JuliaSparse/SparseArrays.jl#244. Lots of operations are defined in term of _mapreduce i.e. sum. We think it makes sense that transpose "forwards" the map reduce to the parent as to not use the abstractarray implementation if the parent has a better implementation. While this is not as important for dense arrays, for special arrays like sparse matrices, this has grave consequences. We think that it makes sense that this should be added to LinearAlgebra.

(didn't run tests, waiting for CI)
@dkarrasch

Fixes JuliaLang/LinearAlgebra.jl#951.

[dkarrasch]

@tkf I'm not familiar with the transducer stuff. I was expecting sum(A') to be much slower than sum(A) for big matrices due to an unfortunate memory reading pattern (as I'm seeing in a simple hand-written for-loop case, iterating over eachindex(A)). But I'm seeing essentially identical runtimes for both cases, even though they take different code paths: sum(A) takes some mapreduce(..., ::IndexLinear, A) route, whereas sum(A') enters the mapfoldl mechanism. To me surprisingly, the approach taken in this PR then doesn't give any benefit on dense matrices.

As stated in the OP, it certainly helps in the sparse case, or just as much, in cases of matrices that have "non-symmetric" reduction methods implemented, whose transposes/adjoints then would go inefficient mapreduce(..., ::IndexCartesian, A) routes.

[martinholters]

I was about to complain that this would change the order in which the reduction happens (while op has to be associative, we don't necessarily expect it to be commutative in general, I think). But apparently, that's already the case:

julia> mapreduce(identity, (x, y) -> 10x+y, copy([1 2; 3 4]'))
1234

julia> mapreduce(identity, (x, y) -> 10x+y, [1 2; 3 4]')
1324

This looks like a bug to me, but code out there might be relying on this behaviour.

Another thing I've noticed:

julia> mapreduce(string, *, copy([1 2; 3 4]'))
"1234"

julia> mapreduce(string, *, [1 2; 3 4]')
ERROR: MethodError: no method matching adjoint(::String)

Why are we applying adjoint to the result of applying string here? OTOH, this works:

julia> reduce(*, map(string, [1 2; 3 4]'))
"1234"

I won't hold up this PR for not fixing things it wasn't meant to fix, but while you're at it, maybe you could?

[dkarrasch]

Actually, I just found that we do have code that does this: 🤦

julia/stdlib/LinearAlgebra/src/adjtrans.jl

Lines 382 to 385 in fa3981b

	Base._mapreduce_dim(f, op, init::Base._InitialValue, A::Transpose, dims::Colon) =
	transpose(Base._mapreduce_dim(_sandwich(transpose, f), _sandwich(transpose, op), init, parent(A), dims))
	Base._mapreduce_dim(f, op, init::Base._InitialValue, A::Adjoint, dims::Colon) =
	adjoint(Base._mapreduce_dim(_sandwich(adjoint, f), _sandwich(adjoint, op), init, parent(A), dims))

The _sandwich pre- and postcomposes op/f with transpose or adjoint, so when you read the parent array, you first account for the implied, say, adjoint, then you apply f to it. Then the adjoints of the sandwiched f and op annihilate, so next we reduce via op the adjoint of the reduced value and the f value. After the reduction, we apply the postcomposed adjoint, which may get annihilated by the outermost adjoint. Perhaps the design is easiest inspected by this example (note the return types):

julia> A = [rand(ComplexF64, 2, 2) for _ in 1:3, _ in 1:3]
3×3 Matrix{Matrix{ComplexF64}}:
 [0.058704+0.652161im 0.733817+0.57694im; 0.805281+0.489643im 0.380885+0.726077im]   …  [0.26835+0.213657im 0.69784+0.0284174im; 0.799889+0.0295452im 0.672902+0.40377im]
 [0.188201+0.998513im 0.954231+0.504686im; 0.0905644+0.143256im 0.867253+0.42945im]     [0.0246882+0.259391im 0.433305+0.917742im; 0.385457+0.934409im 0.0770505+0.746884im]
 [0.771951+0.966885im 0.91714+0.325295im; 0.490785+0.325342im 0.999969+0.335089im]      [0.402092+0.956559im 0.925539+0.988305im; 0.31698+0.379112im 0.804623+0.107361im]

julia> B = A'
3×3 adjoint(::Matrix{Matrix{ComplexF64}}) with eltype LinearAlgebra.Adjoint{ComplexF64, Matrix{ComplexF64}}:
 [0.058704-0.652161im 0.805281-0.489643im; 0.733817-0.57694im 0.380885-0.726077im]  …  [0.771951-0.966885im 0.490785-0.325342im; 0.91714-0.325295im 0.999969-0.335089im]
 [0.571035-0.567911im 0.628074-0.53731im; 0.293206-0.539971im 0.693246-0.26058im]      [0.00273056-0.293979im 0.0554113-0.127508im; 0.184346-0.705327im 0.398427-0.488137im]
 [0.26835-0.213657im 0.799889-0.0295452im; 0.69784-0.0284174im 0.672902-0.40377im]     [0.402092-0.956559im 0.31698-0.379112im; 0.925539-0.988305im 0.804623-0.107361im]

julia> C = copy(B)
3×3 Matrix{Matrix{ComplexF64}}:
 [0.058704-0.652161im 0.805281-0.489643im; 0.733817-0.57694im 0.380885-0.726077im]  …  [0.771951-0.966885im 0.490785-0.325342im; 0.91714-0.325295im 0.999969-0.335089im]
 [0.571035-0.567911im 0.628074-0.53731im; 0.293206-0.539971im 0.693246-0.26058im]      [0.00273056-0.293979im 0.0554113-0.127508im; 0.184346-0.705327im 0.398427-0.488137im]
 [0.26835-0.213657im 0.799889-0.0295452im; 0.69784-0.0284174im 0.672902-0.40377im]     [0.402092-0.956559im 0.31698-0.379112im; 0.925539-0.988305im 0.804623-0.107361im]

julia> sum(B)
2×2 adjoint(::Matrix{ComplexF64}) with eltype ComplexF64:
 2.53179-5.62525im  4.53982-2.97494im
 5.67989-4.6595im   5.53009-4.18912im

julia> sum(C)
2×2 Matrix{ComplexF64}:
 2.53179-5.62525im  4.53982-2.97494im
 5.67989-4.6595im   5.53009-4.18912im

Your first examples indeed indicate that we implicitly assume commutativity of op when we pass to the parent array.

julia> prod(B)
2×2 adjoint(::Matrix{ComplexF64}) with eltype ComplexF64:
 2.64354-6.88695im  7.20743-8.63406im
 7.83816-6.16369im  14.3462-5.03128im

julia> prod(C)
2×2 Matrix{ComplexF64}:
 11.3249-11.4946im  7.58097-1.75261im
 15.9034-12.9864im  9.83046-1.17683im

Ouch! Perhaps we should restrict that behavior to known commutative cases.

[dkarrasch]

@martinholters Can you please take another look when you find the time. I'm not sure we can fix your mapreduce-string example, because that fails in a broader context, see:

julia> ["a" "b"; "c" "d"]'
2×2 adjoint(::Matrix{String}) with eltype Union{}:
Error showing value of type LinearAlgebra.Adjoint{Union{}, Matrix{String}}:
ERROR: MethodError: no method matching adjoint(::String)

[martinholters]

Regarding the original intention of this PR and the current implementation, I cannot offer any helpful feedback. But I'm glad to see the bug I've discovered is fixed.

Regarding the string case, I'd like to point out that this works:

julia> reduce(*, map(string, [1 2; 3 4]'))
"1234"

But as I've said before, I don't want to force this PR into fixing things it wasn't meant to fix. So feel free to ignore this case.

[codecov]

Codecov Report

Merging #46605 (81eb6ef) into master (293031b) will decrease coverage by 0.00%.
The diff coverage is n/a.

❗ Current head 81eb6ef differs from pull request most recent head d6deee2. Consider uploading reports for the commit d6deee2 to get more accurate results

@@            Coverage Diff             @@
##           master   JuliaLang/julia#46605      +/-   ##
==========================================
- Coverage   93.66%   93.66%   -0.01%     
==========================================
  Files         382      386       +4     
  Lines       85268    85567     +299     
==========================================
+ Hits        79868    80147     +279     
- Misses       5400     5420      +20     

Impacted Files	Coverage Δ
compiler/compiler.jl	62.50% <0.00%> (-25.74%)	⬇️
...usr/share/julia/stdlib/v1.9/Profile/src/Profile.jl	68.18% <0.00%> (-19.67%)	⬇️
compiler/typelattice.jl	84.72% <0.00%> (-4.45%)	⬇️
methodshow.jl	87.17% <0.00%> (-3.87%)	⬇️
...e/julia/stdlib/v1.9/NetworkOptions/src/ca_roots.jl	70.00% <0.00%> (-3.34%)	⬇️
...9/SparseArrays/src/solvers/lib/x86_64-linux-gnu.jl	96.82% <0.00%> (-3.18%)	⬇️
float.jl	96.87% <0.00%> (-2.25%)	⬇️
...usr/share/julia/stdlib/v1.9/Unicode/src/Unicode.jl	98.27% <0.00%> (-1.73%)	⬇️
strings/substring.jl	97.16% <0.00%> (-1.42%)	⬇️
...ed/usr/share/julia/stdlib/v1.9/Pkg/src/GitTools.jl	68.83% <0.00%> (-1.30%)	⬇️
... and 102 more

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

[dkarrasch]

Regarding the string case, I'd like to point out that this works:

I'm not sure we can do much about that. Taking adjoints of strings is apparently not intended to work, and the above combination is one that works "by chance", whereas moving the adjoint by one bracket fails for the same reason: reduce(*, map(string, [1 2; 3 4])'). I assume this could be made to work by defining a generic fallback:

adjoint(x) = x
tranpose(x) = x

which I'm not going to propose. A similar generalization of basic operators got reverted a few months ago after causing lots of hassle 😅 .

[dkarrasch]

With the next commit, the following now "works":

julia> mapreduce(string, *, [1 2; 3 4]')
"1324"

julia> mapreduce(string, *, copy([1 2; 3 4]'))
"1234"

Apparently, we don't necessarily need that the eltype of the matrix is such that * commutes, but that * commutes for the return types of f. I guess we should be conservative and restrict to f::Identity for *? Then the fast path would apply to prod(A'), which is perhaps our main concern.

[dkarrasch]

Alright, tests are passing. Thanks everyone for suggestions, more test cases and comments!