Best way to compute matrix quadratic form in Julia

Question

Given a Vector denoted v and Matrix denoted M, what is the fastest method for computing the matrix quadratic form 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 v'*M*v returns a vector containing one element, not a scalar. I was not expecting this behaviour, although have read enough github issue pages to suspect there is a good reason for this behaviour that I'm just not smart enough to see. So, obviously (v'*M*v)[1] will do the job, just wondering if there is a better method...

Answer 1

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

Answer 2

How about a double for? Would save allocations of 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 @inbounds and perhaps @simds would optimize.

Answer 3

Just a note to point out that n(n-1)/2 operations are redundant in the above double loop; if v is a vector on n elements and M is a nxn matrix, the computation of quadratic form v'Mv requires the following operations only:

sum_i=1:n(v[i]^2 * M[i, i]) + sum_j=1:n-1, i=j+1:n(2 * (v[i] * v[j] * M[i, j]))