如何解决在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;
}
}
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。