Column-wise dot product in Eigen C++

Question

Is there an easy way to evaluate the column wise dot product of 2 matrices (lets call them A and B, of type Eigen::MatrixXd) that have dimensions mxn, without evaluating A*B or without having to resort to for loops? The resulting vector would need to have dimensions of 1xn or nx1. Also, I'm trying to do this with Eigen in C++

Answer 1

There are many ways to achieve this, all performing lazy evaluation:

res = (A.array() * B.array()).colwise().sum();
res = (A.cwiseProduct(B)).colwise().sum();

And my favorite:

res = (A.transpose() * B).diagonal();

Answer 2

I did experiment based on @ggael's answer.

MatrixXd A = MatrixXd::Random(600000,30);
MatrixXd B = MatrixXd::Random(600000,30);

MatrixXd res;
clock_t start, end;
start = clock();
res.noalias() = (A * B.transpose()).diagonal();
end = clock();
cout << "Dur 1 : " << (end - start) / (double)CLOCKS_PER_SEC << endl;

MatrixXd res2;
start = clock();
res2 = (A.array() * B.array()).rowwise().sum();
end = clock();
cout << "Dur 2 : " << (end - start) / (double)CLOCKS_PER_SEC << endl;

MatrixXd res3;
start = clock();
res3 = (A.cwiseProduct(B)).rowwise().sum();
end = clock();
cout << "Dur 3 : " << (end - start) / (double)CLOCKS_PER_SEC << endl;

And the output is:

Dur 1 : 10.442
Dur 2 : 8.415
Dur 3 : 7.576

Seems that the diagonal() solution is the slowest one. The cwiseProduct one is the fastest. And the memory usage is the same.

Answer 3

Here is how I'd do it with an Eigen::Map (assuming real matrices, can extend to complex via taking the adjoint), where rows and cols denote the number of rows/columns:

#include <Eigen/Dense>
#include <iostream>

int main()
{
    Eigen::MatrixXd A(2, 2);
    Eigen::MatrixXd B(2, 2);
    A << 1, 2, 3, 4;
    B << 5, 6, 7, 8;

    int rows = 2, cols = 2;

    Eigen::VectorXd vA = Eigen::Map<Eigen::VectorXd>(
                             const_cast<double *>(A.data()), rows * cols, 1);
    Eigen::VectorXd vB = Eigen::Map<Eigen::VectorXd>(
                             const_cast<double *>(B.data()), rows * cols, 1);

    double inner_prod = (vA.transpose() * vB).sum();

    std::cout << inner_prod << std::endl;
}