是否有优化的方法来计算 Eigen++ 中对称“A”的“x^T A x”？

Question

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

Answer 1

你多久计算一次？如果经常出现不同的 x，那么它可能会给出一点 speed-up 来计算 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 的总和。第二项，如果小心使用三角结构，使用原始表达式的一半操作。

如果三角形 matrix-vector 产品必须在 Eigen 过程之外实现，这可能会破坏通用 [=28] 中 BLAS-like 块操作的 efficiency/optimizations =] 产品。最后，算术运算次数的减少可能不会有任何改善。

Answer 2

对于小矩阵，自己写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;
    }
}

是否有优化的方法来计算 Eigen++ 中对称“A”的“x^T A x”？

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

matrix

vector

linear-algebra

numerical-methods

eigen3