Векторизация взвешенного внешнего продукта

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.transpose();

Обратите внимание, что следующее не устанавливает все записи в ноль:

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

Если вы хотите это, вам нужно написать:

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