Efficient Eigen Matrix From Function

Question

I'm trying to build a matrix from a kernel, such that A(i,j) = f(i,j) where i,j are both vectors (hence I build A from two matrices x,y which each row corresponds to a point/vector). My current function looks similar to this:

Eigen::MatrixXd get_kernel_matrix(const Eigen::MatrixXd& x, const Eigen::MatrixXd& y, double(&kernel)(const Eigen::VectorXd&)) {
    Eigen::MatrixXd res (x.rows(), y.rows());
    for(int i = 0; i < res.rows() ; i++) {
        for(int j = 0; j < res.cols(); j++) {
            res(i, j) = kernel(x.row(i), y.row(j));
            }
        }
    }
    return res;
}

Along with some logic for the diagonal (which would in my case likely cause division by zero).

Is there a more efficient/idiometric way to do this? In some of my tests it appears that Matlab code beats the speed of my C++/Eigen implementation (I'm guessing due to vectorization).

I've looked through a considerable amount of documentation (such as the unaryExpr function), but can't seem to find what I'm looking for.

Thanks for any help.

Answer 1

You can use NullaryExpr with an appropriate lambda to remove your for loops:

MatrixXd res = MatrixXd::NullaryExpr(x.rows(), y.rows(),
               [&x,&y,&kernel](int i,int j) { return kernel(x.row(i), y.row(j)); });

Here is a working self-contained example reproducing a matrix product:

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;

double my_kernel(const MatrixXd::ConstRowXpr &x, const MatrixXd::ConstRowXpr &y) {
  return x.dot(y);
}

template<typename Kernel>
MatrixXd apply_kernel(const MatrixXd& x, const MatrixXd& y, Kernel kernel) {
  return MatrixXd::NullaryExpr(x.rows(), y.rows(),
                [&x,&y,&kernel](int i,int j) { return kernel(x.row(i), y.row(j)); });
}

int main()
{
  int n = 10;
  MatrixXd X = MatrixXd::Random(n,n);
  MatrixXd Y = MatrixXd::Random(n,n);
  MatrixXd R = apply_kernel(X,Y,std::ptr_fun(my_kernel));
  std::cout << R << "\n\n";
  std::cout << X*Y.transpose() << "\n\n";
}

If you don't want to make apply_kernel a template function, you can use std::function to pass the kernel.