Pairwise differences between two matrices in Eigen

Question

In matlab/octave pairwise distances between matrices as required for e.g. k-means are calculated by one function call (see cvKmeans.m), to distFunc(Codebook, X) with as arguments two matrices of dimensions KxD.

In Eigen this can be done for a matrix and one vector by using broadcasting, as explained on eigen.tuxfamily.org:

 (m.colwise() - v).colwise().squaredNorm().minCoeff(&index);

However, in this case v is not just a vector, but a matrix. What's the equivalent oneliner in Eigen to calculate such pairwise (Euclidean) distances across all entries between two matrices?

Answer 1

I think the appropriate solution is to abstract this functionality into a function. That function may well be templated; and it may well use a loop - the loop will be really short, after all. Many matrix operations are implemented using loops - that's not a problem.

For example, given your example of...

MatrixXd p0(2, 4);
p0 <<
    1, 23, 6, 9,
    3, 11, 7, 2;

MatrixXd p1(2, 2);
p1 <<
    2, 20,
    3, 10;

then we can construct a matrix D such that D(i,j) = |p₀(i) - p₁(j)|²

MatrixXd D(p0.cols(), p0.rows());
for (int i = 0; i < p1.cols(); i++)
    D.col(i) = (p0.colwise() - p1.col(i)).colwise().squaredNorm().transpose();

I think this is fine - we can use some broadcasting to avoid 2 levels of nesting: we iterate over p₁'s points, but not over p₀'s points, nor over their dimensions.

However, you can make a oneliner if you observe that |p₀(i) - p₁(j)|² = |p₀(i)|² + |p₁(j)|² - 2 p₀(i)^T p₁(j). In particular, the last component is just matrix multiplication, so D = -2 p₀^T p₁ + ...

The blank left to be filled is composed of a component that only depends on the row; and a component that only depends on the column: these can be expressed using rowwise and columnwise operations.

The final "oneliner" is then:

D = ( (p0.transpose() * p1 * -2
      ).colwise() + p0.colwise().squaredNorm().transpose()
    ).rowwise() + p1.colwise().squaredNorm();

You could also replace the rowwise/colwise trickery with an (outer) product with a 1 vector.

Both methods result in the following (squared) distances:

  1 410
505  10
 32 205
 50 185

You'd have to benchmark which is fastest, but I wouldn't be surprised to see the loop win, and I expect that's more readable too.