Define new dot method for abstract vectors

[dalum]

#22374 undefines dot between matrices. The concerns raised in #22220 was that dot(::Matrix, ::Matrix) is presently the Euclidean product. This behaviour will be fixed/broken by #22374 and the Euclidean product between matrices will have to be computed by calling vecdot.

There are two main arguments for adding methods between matrices and numbers to dot as in this PR:

• numpy defines dot in a similar way (except for the conjugation) between arrays and numbers and arrays and arrays
• this allows exploitation of the nested calls to dot between the entries in vectors, where the entries are arrays or numbers to write sums of products as dot products of vectors

(I don't know how to make #22220 point to this new branch, so I'm making a new PR. Sorry for the noise.)

[andreasnoack]

@eveydee Sorry for obstructing your PRs. It is great that you contribute. However, I think we should go for the complete separation of dot and vecdot and simply define

dot(x::AbstractVector, y::AbstractVector) = sum(xx'yy for (xx,yy) in zip(x,y)

instead of having dot call vecdot. I believe that would allow everybody to do what they'd like to do with either dot or vecdot.

[dalum]

@andreasnoack I'm really glad to receive feedback and hope that I contribute more than cause trouble.

I think yours is a much cleaner solution, but it looks like it comes with a performance cost which is quite significant for normal short vectors:

julia> dott(x::AbstractVector, y::AbstractVector) = sum(xx'yy for (xx,yy) in zip(x,y))
julia> a = 1:4
julia> @benchmark dot(a, a)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     22.109 ns (0.00% GC)
  median time:      22.961 ns (0.00% GC)
  mean time:        23.265 ns (0.00% GC)
  maximum time:     91.116 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     997

julia> @benchmark dott(a, a)
BenchmarkTools.Trial: 
  memory estimate:  128 bytes
  allocs estimate:  3
  --------------
  minimum time:     64.720 ns (0.00% GC)
  median time:      68.734 ns (0.00% GC)
  mean time:        84.562 ns (15.93% GC)
  maximum time:     3.641 μs (97.08% GC)
  --------------
  samples:          10000
  evals/sample:     980

I would imagine this is the reason dot was defined to be equal to vecdot in the first place. 😕

[stevengj]

Presumably @andreasnoack's suggestion was mainly for the semantics. We'd probably still want specialized methods for specific types (e.g. arrays of numbers, cases that would call BLAS, arrays of vectors, arrays of matrices, etcetera), and the fallback method should probably be implemented by explicit loops rather than summing a generator.

[andreasnoack]

@andreasnoack I'm really glad to receive feedback and hope that I contribute more than cause trouble.

Good to hear that.

As @stevengj suggests, I'm mainly after the semantics here.

It would be great if we had two argument mapreduce. I think we could do something like

dot(x::AbstractVector{<:Number}, y::AbstractVector{<:Number}) = vecdot(x, y)
function dot(x::AbstractVector, y::AbstractVector)
  if length(x) != length(y)
    throw(ArgumentError("something"))
  end
  return isa(first(x)'first(y)) ? vecdot(x,y) : sum(xx'yy for (xx,yy) in zip(x,y))
end

and avoid the overhead in the most important cases. The second version would also be fairly fast in the Number case but it has a branch and can't handle zero length arguments. We could try to get the zero argument version working in more cases but it wouldn't work for Vector{Vector} anyway so I'm not sure it is worth it.

[dalum]

I have tried to optimise the dot method to the best of my knowledge, by studying the implementations of foldl, reduce, etc. and benchmarking different approaches. I'm not entirely sure if calling @inbounds is safe, but it does make it slightly faster on my machine. Curiously, the BLAS methods are a bit slower for short vectors of numbers, but have a cross-over around vector lengths of 16.

[StefanKarpinski]

Curiously, the BLAS methods are a bit slower for short vectors of numbers, but have a cross-over around vector lengths of 16.

This is typical: BLAS calls have a fair amount of overhead, but are highly optimized, so the larger the input, the better they fare; the flip side is that they do less well for smaller inputs.