Is there an optimised way to calculate `x^T A x` for symmetric `A` in Eigen++?

Question

Given a symmetric matrix A and a vector x, I often need to calculate x^T * A * x. I can do it Eigen++ 3.x with x.transpose() * (A * x), but that doesn't exploit the information that x is the same on both sides and A is symmetric. Is there a more efficient way to calculate this?

Answer 1

For small matrices, writing the for loop myself seems to be faster than relying on Eigen's code. For large matrices I got good results using .selfAdjointView:

double xAx_symmetric(const Eigen::MatrixXd& A, const Eigen::VectorXd& x)
{
    const auto dim = A.rows();
    if (dim < 15) {
        double sum = 0;
        for (Eigen::Index i = 0; i < dim; ++i) {
            const auto x_i = x[i];
            sum += A(i, i) * x_i * x_i;
            for (Eigen::Index j = 0; j < i; ++j) {
                sum += 2 * A(j, i) * x_i * x[j];
            }
        }
        return sum;
    }
    else {
        return x.transpose() * A.selfadjointView<Eigen::Upper>() * x;
    }
}