Kernel Matrix

library(kernlab)
N = 10000
p = 5
b = 1000
X = matrix(rnorm(N*p), ncol = p)
center = X[sample(1:N, b),]
sigma = 1
kernel_X <- kernelMatrix(rbfdot(sigma=1/(2*sigma^2)), X, center)

library(Rcpp)
library(RcppArmadillo)
library(inline)
settings <- getPlugin("Rcpp")
settings$env$PKG_CXXFLAGS <- paste('-fopenmp', settings$env$PKG_CXXFLAGS)
settings$env$PKG_LIBS <- paste('-fopenmp -lgomp', settings$env$PKG_LIBS)
do.call(Sys.setenv, settings$env)
sourceCpp(code = '
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <omp.h>
using namespace Rcpp;
using namespace arma;
// [[Rcpp::export]]

NumericMatrix kernelMatrix_cpp(NumericMatrix Xr, NumericMatrix Centerr, double sigma) {
    omp_set_num_threads(omp_get_max_threads());
    uword n = Xr.nrow(), b = Centerr.nrow(), row_index, col_index;
    mat X(Xr.begin(), n, Xr.ncol(), false);
    mat Center(Centerr.begin(), b, Centerr.ncol(), false);
    mat KerX(n, b);
    #pragma omp parallel private(row_index, col_index)
    for (row_index = 0; row_index < n; row_index++)
    {
        #pragma omp for nowait
        for (col_index = 0; col_index < b; col_index++)
        {
            KerX(row_index, col_index) = exp(sum(square(X.row(row_index) - Center.row(col_index))) / (-2.0 * sigma * sigma));
        }
    }
    return wrap(KerX);
}')
kernel_X_cpp = kernelMatrix_cpp(X, center, sigma)
all.equal(kernel_X@.Data, kernel_X_cpp)
sum(abs(kernel_X@.Data - kernel_X_cpp) < 1e-12)

library(rbenchmark)
benchmark(kernelMatrix_cpp(X, center, sigma), kernelMatrix(rbfdot(sigma=1/(2*sigma^2)), X, center), columns=c("test", "replications","elapsed", "relative"), replications=50, order="relative")