Hi, this is probably dumb, but I want to know if there is some standard way of doing the dot product. I know that
a'b
will give a correct answer, though I do not know if there is some standard when extending it to matrix multiplication (I mean if there is a used operation besides using LinearAlgebra). Some sources cite \cdot but apparently that is no longer defined on the general library.
If you import LinearAlgebra you should have \cdot. I wasn’t aware that a’b works without that package? (But that may be because I import LinearAlgebra whenever I need linear algebra stuff such as dot products)
The standard way is to use the dot function from the LinearAlgebra standard library.
Julia 1.11 is somewhat helpful about finding it.
julia> dot
ERROR: UndefVarError: `dot` not defined in `Main`
Suggestion: check for spelling errors or missing imports.
Hint: a global variable of this name may be made accessible by importing LinearAlgebra in the current active module Main
julia> using LinearAlgebra
help?> dot
search: dot do cot det stdout ldlt acot sort @doc edit cotd coth
dot(x, y)
x â‹… y
Compute the dot product between two vectors. For complex vectors, the first vector is conjugated.
dot also works on arbitrary iterable objects, including arrays of any dimension, as long as dot is defined
on the elements.
[...]
I find this definition in wiki
Algebraically, the dot product is the sum of the products of the corresponding entries of the two sequences of numbers.
Translating into julia, it is
function dot_product(x, y) return sum(x .* y) end
But due to the fact that we haven’t restricted the type of args, this function may have some other (pleasantly) unexpected usage. Therefore, I would prefer a more appropriate name for it, see this.
My larger point is: there is no foolproof method. And we should not hope for this.
As an example, if sparsity exists, we may sometimes need to write our own code, see this.
In general, the adjoint (') or transpose are versatile, not restricted to only offering a scalar output (which in math sense, inner product). I think these are already enough.
If there is anyone who knows more about the underlying performance issue, just correct me 
I think you’re just spreading confusion here needlessly.
The dot(x,y) function from the LinearAlgebra package (and its synonym x â‹… y) is the standard interface for inner products in Julia. It has methods for inner products of matrices as well. And it is also extended to methods for other types of objects (e.g. sparse vectors) in other packages like SparseArrays.
a'b or a' * b is equivalent to calling dot(a, b) if these are two “column vectors” (1d arrays) of scalars [*see below].
It’s not clear what “extension” to matrix multiplication you are asking for. A'B for two matrices computes the matrix product \overline{A^T} B. dot(A, B) computes the Frobenius inner product. These are two different things, both of which are useful depending on what you are doing.
Numpy muddied the waters even further because numpy.dot is ordinary matrix multiplication (A * B in Julia) for two 2d arrays (and is a type of tensor contraction more generally). Naming things is hard.
PS. Note that sum(x .* y), besides being inefficient, is also not an inner product if the elements aren’t real: a true inner product needs a complex conjugation, which is what dot(x, y) does. (There was reasonable disagreement about whether dot in Julia should be unconjugated, and we should have some other name like inner(x,y) for the true inner product, but nowadays Julia is past that bikeshedding.)
[*] There’s some subtlety here because linear-algebra notation sometimes glosses over the isomorphism between a 1-row matrix and a linear form or “dual” vector, or between a 1 \times 1 matrix and a scalar. As a programming language, Julia had to define this more explicitly; see also the talk Taking Vector Transposes Seriously by @jiahao.
I watched that video. I prefer the current (julia’s) behavior, comparing to other languages.
Yes, dot(x, y) calls BLAS, which is specialized.
Looks like sum(x .* y) will call mapreduce and therefore resolves into standard + and * operations (perhaps calling standard C functions?).
Moreover, it appears that dot(x, A, y) also outperforms x' * A * y, at least in terms of allocations, sometimes also in time (can someone do a test?).
It first computes x .* y, allocating a new array, and then calls sum on this new array. This is much less efficient than doing a single pass over the arrays and accumulating the sum without allocating anything, as dot does.
Is this dispatch as expected? I don’t see the special code owing to SparseArrays.jl.
julia> using LinearAlgebra, SparseArrays
julia> x = sprand(100, .1);
julia> A = rand(100, 100);
julia> @which dot(x, A, x)
dot(x::AbstractVector, A::AbstractMatrix, y::AbstractVector)
@ LinearAlgebra K:\julia-1.11.4\share\julia\stdlib\v1.11\LinearAlgebra\src\generic.jl:931
julia> x = Vector(x); # as a comparison
julia> @which dot(x, A, x) # the same method
dot(x::AbstractVector, A::AbstractMatrix, y::AbstractVector)
@ LinearAlgebra K:\julia-1.11.4\share\julia\stdlib\v1.11\LinearAlgebra\src\generic.jl:931
Looks like the 3-arg dot(x, A, y) doesn’t have an optimized method for sparse vectors yet.
Not related to OP’s question, but I am wondering what is the best way to do this. Is transpose(x)y the best?
what is “this”? just dot(x, y) ?
“this” is the dot product without conjugating the first vector, I guess.
What is “best” depends on to what extent the corresponding method is optimized.
For regular numerical operations, dot will call BLAS, and is therefore “the best”.
You can write transpose(x) * y or x' * y resembling the textbook, but you don’t have to.
dot(x, y) is considered “standard”.
On the other hand, transpose and ' are more versatile so that they can reach more abundant results. e.g., you can do [rand(3, 3) for _ in 1:4]' * rand(4), which is not an inner product but is also meaningful itself. If we look into this operation, we find that it essentially calls
sum(uu*vv for (uu, vv) in zip(u, v))
inside *(u::AdjOrTransAbsVec, v::AbstractVector) = _dot_nonrecursive(u, v).
Do you have a particular case with complex vectors where you specifically need to avoid conjugation?
If the vectors are real-valued, you should use x' * y, not transpose(x) * y.
For vectors with a real element type, both x' * y and transpose(x) * y are equivalent: both call dot(x, y).
Sorry for the confusion. Yes, I want to multiply two vector element-wise and sum the result, like the following
res = 0.0im
for i = eachindex(x)
res += x[i]*y[i]
end
When x is real, I can use dot, but sometimes I need this for complex vectors.
The above code is fine, but a little lengthy, so I’m wondering if there is a better solution.
Surprisingly, I found the functions cdotu and zdotu exist in LAPACK, which perform this operation.
I really hope these functions can be included, maybe through extending the current dot function. (I had trouble with the dot function before because the name for me suggests only element-wise product. I feel inner implies the conjugation involved.)
If you use transpose(x) * y (as suggested above) it will do what you want. It calls dot(x, y) for real arrays, calls BLAS.dotu(x, y) for BLAS-compatible complex types, and otherwise uses a sum(xx*yy for (xx, yy) in zip(x, y)) fallback.
Great! Thanks for the information. Now I know what to do when I need dot product without conjugation.