I need to calculate a the error function V

V = Σi Σj X[i] X[j] σ[i][j]

Where σ[i][j] is a given matrix and I need a relatively fast solution. To do this I want to create another matrix Y where

Y[i][j] = X[i]*X[j]

So I that I can simply sum over Y * σ. Is there a good way to implement this with numpy functions?

So far I have tried to meshgrid(X, X) and then apply np.prod to each row, however this did yield the expected output and would have required a for-loop in python.