对于特征++中的对称"A"，是否有一种优化的方法来计算"x^T A x"

Question

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

Answer 1

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

或者可能更简单，如果你分解A=L+D+L.T，其中L是严格的下三角，D是对角线，那么

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

其中第一项是CCD_ 12上的和。第二个术语，如果仔细使用三角形结构，则使用原始表达式的一半运算。

如果三角矩阵向量积必须在Eigen过程之外实现，这可能会破坏通用矩阵向量积中类似BLAS的块运算的效率/优化。最终，算术运算次数的减少可能不会带来任何改进。

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;
}
}