make (*)(::Adjoint{Vector}, ::Vector) non-recursive

[dkarrasch]

Closes #35174.

These changes don't break any tests, so I believe this was more oversight than deep intention. With this PR, I get

julia> using LinearAlgebra

julia> σ = [[0  1
                    1  0],
                   [0 -im
                    im 0],
                   [1  0
                    0  -1]]
3-element Array{Array{Complex{Int64},2},1}:
 [0 + 0im 1 + 0im; 1 + 0im 0 + 0im]
 [0 + 0im 0 - 1im; 0 + 1im 0 + 0im]
 [1 + 0im 0 + 0im; 0 + 0im -1 + 0im]

julia> [1,2,3]'* σ
2×2 Array{Complex{Int64},2}:
 3+0im   1-2im
 1+2im  -3+0im

which used to throw before.

[dkarrasch]

Currently, we don't have any tests for arrays of arrays related to the issue here.

[ericphanson]

I don’t think the old behavior is really coherent; as I see it, the fact that *(::Adjoint, ::Matrix) dispatched to dot was a leaking implementation detail. And I don’t see how one could formulate a coherent system for array-of-array operations from which you would get the rule that *(::Adjoint{Vector{T}}, Vector{Matrix{T}}) where T <: Number gives a number.

To me, the general rule should be that e.g. *(::Matrix{Matrix{T}}, ::Matrix{Matrix{T}}) where T <: Number should dispatch (semantically at least) to a generic matmul that uses * and + for Matrix{T}, and similarly for other linear algebra operations. In other words, it should just lower one level at a time by using the operations as defined for the lower level.

The recursive dot fits this rule if you define dot(u,v) = sum( dot(u[i], v[i]) for i in eachindex(u,v)) and dot for numbers as conjugated multiplication. And the fix here, as I understand it, fits this rule too.

[simeonschaub]

I would propose a small change to @ericphanson's proposal, to change that to (*)(u::Adjoint{Vector}, v::Vector) = sum(u[k]'v[k] for k in eachindex(u,v)), which would work for cases like this without any other changes to dot for matrices and numbers and would be consistent with the recursive behavior of Adjoint everywhere else.

[dkarrasch]

Oops, I overlooked one error. Fixing it is instructive, as the failing test really sticks out of the surrounding tests. I found that it was changed in #27401, but the thread is very long. I'm pinging @ranocha to check what was the rationale there.

[dkarrasch]

@simeonschaub I think that the recursive transposition behavior is coded into the getindex for Adjoints.

[simeonschaub]

Oops, yes, of course! Then I definitely agree with the implementation here

[dkarrasch]

Actually, currently we have

(*)(u::TransposeAbsVec, v::AdjointAbsVec{<:Any,<:TransposeAbsVec}) =
    sum(uu*vv for (uu, vv) in zip(u, v))
(*)(u::AdjointAbsVec,   v::AdjointAbsVec{<:Any,<:TransposeAbsVec}) =
    sum(uu*vv for (uu, vv) in zip(u, v))
(*)(u::TransposeAbsVec, v::TransposeAbsVec{<:Any,<:AdjointAbsVec}) =
    sum(uu*vv for (uu, vv) in zip(u, v))
(*)(u::AdjointAbsVec,   v::TransposeAbsVec{<:Any,<:AdjointAbsVec}) =
    sum(uu*vv for (uu, vv) in zip(u, v))

Should we simplify the *(u::AdjOrTransAbsVec, v::AbstractVector) method to the same

sum(uu*vv for (uu, vv) in zip(u, v))

without bounds checking and require_one_based_indexing? Feels strange to be so strict in one case and so relaxed in so many others...

[ranocha]

Oops, I overlooked one error. Fixing it is instructive, as the failing test really sticks out of the surrounding tests. I found that it was changed in #27401, but the thread is very long. I'm pinging @ranocha to check what was the rationale there.

The intention of #27401 was to make dot really an inner/scalar product for various things and use the behavior of the old vecdot (and similarly for norm). Hence, dot(Xv2, Xv2) ≈ norm(Xv2)^2 is the desired behavior. The test

julia/stdlib/LinearAlgebra/test/matmul.jl

Line 553 in cc60547

@test Xv2'*Xv2 ≈ norm(Xv2)^2

reflects that since *(::Adjoint, ::Matrix) dispatched to dot, which is not necessarily an optimal choice, as mentioned by @ericphanson above. I didn't change *(::Adjoint, ::Matrix), only dot in #27401. Hence, I think the rational behind this PR is fine.

[codecov-io]

Codecov Report

Merging #35257 into master will decrease coverage by 0.00%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master   #35257      +/-   ##
==========================================
- Coverage   86.79%   86.79%   -0.01%     
==========================================
  Files         344      344              
  Lines       63715    63716       +1     
==========================================
- Hits        55303    55301       -2     
- Misses       8412     8415       +3     

Impacted Files	Coverage Δ
base/generator.jl	100.00% <ø> (ø)
base/broadcast.jl	85.40% <100.00%> (+0.03%)	⬆️
base/compiler/optimize.jl	83.90% <100.00%> (ø)
stdlib/Base64/src/encode.jl	84.76% <0.00%> (-3.81%)	⬇️
base/grisu/fastprecision.jl	84.90% <0.00%> (-1.89%)	⬇️
base/missing.jl	87.77% <0.00%> (-0.56%)	⬇️
base/task.jl	82.89% <0.00%> (-0.38%)	⬇️
stdlib/SparseArrays/src/linalg.jl	96.33% <0.00%> (-0.32%)	⬇️
base/bitset.jl	91.74% <0.00%> (ø)
stdlib/SparseArrays/src/sparsematrix.jl	90.23% <0.00%> (+0.05%)	⬆️
... and 3 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update f449765...3a4f267. Read the comment docs.

[marius311]

Hi guys, I ran into a few issues with this PR as I was upgrading my code to 1.5:

1. I don't think the dimentionality checking should be removed. The 1.5 behavior now is just bizarre:

   # julia 1.4
   julia> [1, 2, 3]' * [1., 2.]
   ERROR: DimensionMismatch("first array has length 3 which does not match the length of the second, 2.")
   
   # julia 1.5rc
   julia> [1, 2, 3]' * [1., 2.]
   5.0

   and lots of dot implementations in LinearAlgebra/src/generic.jl have dimension checking.

2. E.g. (*)(::Adjoint{Vector{Int}}, Vector{Float}) uses slightly different summation code rather than the existing https://github.com/JuliaLang/julia/blob/release-1.4/stdlib/LinearAlgebra/src/generic.jl#L895-L897. I've benchmarked the two for simple builtin types (as you may have too) and I don't see any performance regressions, but are we sure there's no types this might cause an issue?

3. The fact that mixed-but-still-scalar eltypes go down different codepaths introduces what feels like an unnecessary breaking API change because whereas previously it sufficed to override dot(x,y) for your custom types, now you have to do both dot(::MyVector, ::MyVector) and (*)(::Adjoint{MyVector}, ::MyVector}. A solution might be

   *(u::AdjointAbsVec{<:Number}, v::AbstractVector{<:Number})

   instead of

   *(u::AdjointAbsVec{T}, v::AbstractVector{T}) where {T <: Number}

   in line 270, but its not totally clear.

My impression is that while the original issue certainly deserves to be solved, maybe this PR needed to be iterated on a little more before merging? Would it be possible to fix it quickly before 1.5, or maybe revert and iterate on it a bit? For me certainly the breaking API change, while I think technically allowed by Julia's semver rules as they pertain to stdlibs, seems worrisome.

[KristofferC]

cc @dkarrasch