How to vectorize the evaluation of a kernel function.

I have a kernal function which is defined for . And now I have to compute a matrix for m points and n points , where K is given by . It is direct when using two for loop. But how can I vectorize the evaluation? For example, I tried

k_fun = @(x, y) 1 / norm(x - y);
d = 2;  % Make d = 1 if you want it runs correctly.
m = 100;
n = 100;
x_points = rand(m, d);
y_points = rand(n, d);
% The following code is the two for loop version.
K = zeros(m, n);
for i = 1 : m
    for j = 1 : n
        K(i, j) = k_fun(x_points(i,:), y_points(j, :));
    end
end
% The folloing code works when d = 1, but when d > 1 it failes.
K = k_fun(x_points, y_points');
% When d > 1, the error is "Arrays have incompatible sizes for this
% operation."

When , it gives the result I want, But for , it failes. How can I improve it?

[Jan]

"It is direct when using two for loop" - then it would be useful, if you post this solution, which would uniquely define, what you want to achieve.

"it fails" - for me it fails due to the unfound function "ran()". Should this be "rand()"? Post the error message, if you mention an error.

[Jingyu]

@Jan Thanks for pointing out my mistake. I correct the typos in the code. And the code I posted is only correct when you modify . When m and n are large, two for loop for this evaluation is expensive. You can compare the behavior when . What I want to know is that is there any treatment like the case .

[Jan]

@Jingyu: What is norm(x - y) if x is [M x 2] and y is [N x 2]? If N~=M it does not work for d=1 also.

Note, that you need the elementwise division ./ to divide by an array. What does the error message tell you?

[Jingyu]

@Jan You must make sure x and y are the vectors of the same length. If you given me an x and an y, then actually I want to compute a matrix K whose elements is given by where k is the kernel function. I have told you the code will occur error. I just want to know, for example, how to vectorize when the dimension is larger than 1.

[Jan]

@Jingyu: "I have told you the code will occur error" - yes, you did. Please insert the error message also in future questions.

While your code is vectorized already, you let the readers guess, what you want to achieve. All we know, is that your kernal function is "special" and the not working code.

[Jingyu]

@Jan Thanks for your advice. I have changed the code.

Answer 1

[Matt J]

K=1./pdist2(x_points,y_points);

[Jingyu]

Thanks for your answer. But the kernel k in my code is so special. Is there any method for a general k?

[Torsten]

Then supply the special kernel for which you need the implementation.

[Jingyu]

@Torsten I mean for any kernel. Do we have a method? The treatment is independent of the kernel function,

[Torsten]

I thought you wanted optimized code.

Can you vectorize your kernel function such that it accepts vectors or matrices as inputs for x and y ?

[Matt J]

Is there any method for a general k?

No, there is no universal vectorization method. Each specific k must be considered individually.

[Jingyu]

@Torsten Yes. The main cost of my function is the two for loop so I want to optimize it. The kernel function accepts inputs only for vectors.

[Jingyu]

@Matt J Thanks! Then, if the kernel function has the form of , do we have some method to handle it?

[Torsten]

It's not possible to optimize the call to your kernel function if it only accepts two vectors (one for x and one for y) . Or did you mean two matrices of vectors of the same length (one matrix for the x vectors and the other for the y-vectors) ?

[Jingyu]

@Torsten Sorry. The latter (two matrices of vectors of the same length).

[Torsten]

Sorry. The latter (two matrices of vectors of the same length).

If your function does so, you can call it using

k_fun(x_points, y_points);

because x_points and y_points are matrices of vectors of the same length.

But the kernel function you wrote does not accept matrices of vectors of the same length.

[Matt J]

Then, if the kernel function has the form of , do we have some method to handle it?

It depends on and how "special" it is.

[Jingyu]

@Torsten Oh. My fault. Then it only accepts two vectors...Does it mean I must write two for loops?

[Torsten]

Does it mean I must write two for loops?

Yes.

[Jingyu]

@Torsten Thanks for your answer！