R functional programming: using apply family to calculate kernel matrix for gaussian processes

Question

My Minimal Working Example

Here is my minimal working example (notice I've simplified the mathematics). Suppose I have a function of two variables where x and y are two vectors of the same dimension.

kernel_func <- function(x, y){
    return(sum((x - y)^2))
}

I also have two matrices with different number of rows, but same number of columns.

• A matrix X of dimension n times d
• A matrix Y of dimension m times d

Now I would like to obtain a matrix, call it K, whose i,j element is calculated by passing the ith row of X as the first argument to kernel_func and the jth row of Y as the second argument. That is

kernel_func(X[i, ], Y[j, ])

How can I write a concise piece of code that does this, hopefully using apply, lapply, mapply or similar?

Stupid MWE

Create the two matrices with the same number of columns

X = matrix(1:9, nrow=3, ncol=3)
Y = matrix(1:12, nrow=4, ncol=3)

Initiate the K matrix with zeros

K <- matrix(0, nrow(X), nrow(Y))

Use a double loop to create the matrix

for (i in rep(1:nrow(X), 4)){
    for (j in 1:nrow(Y)) {
        K[i, j] = kernel_func(X[i, ], Y[j, ])
    }
}

Answer 1

Here are two approaches to make it:

• Solution with apply():

K <- t(apply(X, 1, function(p) apply(Y, 1, function(q) kernel_func(p,q))))

• Solution with expand.grid():

K <- matrix(apply(expand.grid(1:nrow(X),1:nrow(Y)),1, 
                  function(k) kernel_func(X[k[1],],Y[k[2],])),nrow = nrow(X))

Output

> K
     [,1] [,2] [,3] [,4]
[1,]    5   14   29   50
[2,]    2    5   14   29
[3,]    5    2    5   14