加权外积的矢量化

Vectorization of weighted outer product

本文关键字：矢量化 加权      更新时间：2023-10-16

我希望加快近似加权协方差的计算。

具体来说，我有一个Eigen::VectorXd(N) w和一个Eigen::MatrixXd(M,N) points.我想计算w(i)*points.col(i)*(points.col(i).transpose())的总和。

我正在使用 for 循环，但想看看我是否可以走得更快：

Eigen::VectorXd w = Eigen::VectorXd(N) ;
Eigen::MatrixXd points = Eigen::MatrixXd(M,N) ;
Eigen::MatrixXd tempMatrix = Eigen::MatrixXd(M,M) ;
for (int i=0; i < N ; i++){
tempMatrix += w(i)*points.col(i)*(points.col(i).transpose());
}

期待看到能做些什么！

以下方法应该有效：

Eigen::MatrixXd tempMatrix; // not necessary to pre-allocate
// assigning the product allocates tempMatrix if needed
// noalias() tells Eigen that no factor on the right aliases with tempMatrix
tempMatrix.noalias() = points * w.asDiagonal() * points.adjoint();

或直接：

Eigen::MatrixXd tempMatrix = points * w.asDiagonal() * points.adjoint();

如果M真的很大，那么只计算一端并复制它(如果需要(可能会快得多：

Eigen::MatrixXd tempMatrix(M,M);
tempMatrix.triangularView<Eigen::Upper>() = points * w.asDiagonal() * points.adjoint();
tempMatrix.triangularView<Eigen::StrictlyLower>() = tempMatrix.adjoint();

请注意，对于非复杂标量，.adjoint()等效于.transpose()，但对于前者，如果points代码也可以工作，而结果则MatrixXcd相反(如果结果必须是自伴随的，则w必须是真实的(。

另请注意，以下内容(来自原始代码(并未将所有条目设置为零：

Eigen::MatrixXd tempMatrix = Eigen::MatrixXd(M,M);

如果你想要这个，你需要写：

Eigen::MatrixXd tempMatrix = Eigen::MatrixXd::Zero(M,M);
相关文章：