b = np.arange(9).reshape(3, 3)
np.diag(np.dot(np.dot(a, b), a.T))
np.einsum('ij,ji->i', np.dot(a, b), a.T)
np.einsum('ij,ij->i', np.dot(a, b), a)
Is there a numpy/scipy dot product, calculating only the diagonal entries of the result?