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++

Krp

Element-wise multiply, then sum? In MATLAB it would be sum(A .* B). Eigen provides those operations, but I don't know the exact names of the calls. – Ghee

nice one! That should work. Thanks! – Krp

Can use an "Eigen::Map" to reshape the matrices into vectors, then take their inner product. – Drover

@Zedd: If my comment provides the hint you needed to make working code, please answer your own question showing future visitors how it is done using Eigen functions. – Ghee

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

Yemen

One curious question: By appending diagonal to the dot product operation, would Eigen only calculate the (ith row of A) dot product (ith col of B) internally so that it would avoid unnecessary calculation? For these three implementation which one would be faster? (Since the 1st and 2nd would avoid unnecessary calculation). – Lizethlizette

As written in the answer they all perform lazy evaluation, so the last one does compute the diagonal entries only. They are expected to generate roughly the same code. You can check on compiler explorer. – Yemen

Will (A.transpose() * B).diagonal(); have the right asymptotic performance? Meanwhile, how come A.colwise().dot(B.colwise()) doesn't work? – Prurigo

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.

Lizethlizette

@ggael, please feel free to comment. – Lizethlizette

I performed a very similar test but used 16x3 single precision matrices and performed each operation a million times. Using clang++ -O3 I had times 175, 129, 131 milliseconds which is similar to your results. – Cuttle

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

Drover

Now I see that you want a vector as the result? That is, are you after a \vect{result} where each component is A(i,)*B(:,i)' ? – Drover

Yeah, thats right. I'm not sure how this would translate into the code you have above. – Krp