加权外积的向量化

Question

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

具体来说，我有一个 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());
}

期待看到可以做什么！

Answer 1

以下应该有效：

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

加权外积的向量化

Vectorization of weighted outer product

c++

vectorization

eigen