Best way to calculate matrix quadratic form in Julia

Question

Given Vector, designated vand Matrixdesignated as Mwhat is the fastest way to calculate the quadratic form of a matrix in Julia, i.e. v'Mv? What about the most elegant?

Note: I would like the return value to be a scalar. Interestingly, if v = rand(3)and M = rand(3, 3), then it v'*M*vreturns a vector containing one element, not a scalar. I did not expect this behavior, although I read enough pages of the github problem to suspect that there is a good reason for this behavior that I am not smart enough to see. So, it’s obvious that it (v'*M*v)[1]will do the job just by wondering if there is a better method ...

Answer 1

One parameter that returns a scalar dot(v, M*v).

Answer 2

How about a double for? Saves the allocation of the intermediate vector. For completeness, in code:

vMv{T}(v::Vector{T},M::Matrix{T}) = begin
    size(M) == (size(v,1),size(v,1)) || throw("Vector/Matrix size mismatch")
    res = zero(T)
    for i = 1:size(v,1)
        for j = 1:size(v,1)
            res += v[i]*v[j]*M[i,j]
        end
    end
    res
end

Adding some @inboundsand maybe @simdwill be optimized.