eigen - 使用 Eigen 的指数平均

Question

考虑以下代码。

const int N = 100;
const float alpha = 0.9;

Eigen::MatrixXf myVec = Eigen::MatrixXf::Random(N,1);
Eigen::MatrixXf symmetricMatrix(N, N);
for(int i=0; i<N; i++)
    for(int j=0; j<=i; j++)
        symmetricMatrix(i,j) = symmetricMatrix(j,i) =   i+j;

symmetricMatrix *= alpha;
symmetricMatrix += ((1-alpha)*myVec*myVec.adjoint());

它本质上实现了指数平均。
我知道最后一行可以通过以下方式进行优化。

symmetricMatrix_copy.selfadjointView<Eigen::Upper>().rankUpdate(myVec, 1-alpha);

我想知道我是否可以以有效的方式组合最后两行。
简而言之，我想计算 A = alpha*A+(1-alpha)*(x*x') .

Answer 1

最重要的是，您应该声明 myVec如 Eigen::VectorXf ，如果它保证是一个向量。并确保您使用 -O3 -march=native -DNDEBUG 进行编译.

您可以尝试这些替代方案(我使用 A 和 v 来节省输入)，哪个最快可能取决于您的问题大小和您的 CPU:

A.noalias() = alpha * A + (1.0f-alpha)*v*v.adjoint();
A.noalias() = alpha * A + (1.0f-alpha)*v.lazyProduct(v.adjoint());
A.noalias() = alpha * A + ((1.0f-alpha)*v).lazyProduct(v.adjoint());
A.noalias() = alpha * A + v.lazyProduct((1.0f-alpha)*v.adjoint());

A.triangularView<Eigen::Upper>() = alpha * A + (1.0f-alpha)*v*v.adjoint();
// or any `lazyProduct` as above.

不幸的是， .noalias()和 .triangularView()暂时无法合并。

你也可以考虑计算这个:

A.selfadjointView<Eigen::Upper>().rankUpdate(v, (1.0f-alpha)/alpha);

和每个 N迭代扩展您的 A矩阵由 pow(alpha, N)