在Eigen ++中，是否存在优化的方法来为对称`A`计算`x ^ T A x`？

问题描述

鉴于对称矩阵A和向量x，我经常需要计算x^T * A * x。我可以使用x.transpose() * (A * x)来实现Eigen ++ 3.x，但这并不能利用x双方相同且A是对称的信息。有没有更有效的方法来计算这个？

解决方法

您多久计算一次？如果经常使用不同的x，则可能会加快计算A的Cholesky或LDLT分解的速度，并使用三角形矩阵与向量的乘积仅需要一半的乘法。

或者甚至更简单，如果分解A=L+D+L.T，其中L严格是下三角形，而D是对角线，那么

x.T*A*x = x.T*D*x + 2*x.T*(L*x)

其中第一项是d[k]*x[k]**2上的总和。如果仔细使用三角形结构，第二项将使用原始表达式的一半运算。

如果必须在Eigen过程之外实施三角矩阵向量乘积，则这可能会破坏通用矩阵向量乘积中类似BLAS的块运算的效率/优化。最后，减少算术运算的数量可能没有任何改善。

,

对于小型矩阵，我自己编写for循环似乎比依赖Eigen的代码要快。对于大型矩阵，使用.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;
    }
}