c++ - 加权外积的向量化

Question

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

具体来说，我有一个 Eigen::VectorXd(N) w 和一个 Eigen::MatrixXd(M,N) 点。我想计算 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()，但对于前者，如果 指向 和 MatrixXcd 代替的结果(w 必须仍然是实数，如果结果必须是自伴随的)。

另外，请注意以下(来自您的原始代码)不会将所有条目设置为零:

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

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

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