fitrsvm doesn't vectorize my custom kernel

I'm working on an iterative function approximation algorithm. As part of that algorithm I use fitrsvm to update a kernel regression each iteration. This all works great when I use MATLAB's built in kernels. However, once I venture outside of MATLAB's kernels and start using my own the runtime of fitrsvm increases by nearly 4 times. That is, the runtime increases when I call:

m = fitrsvm(X, Y, 'KernelFunction','my_kernel');

To try to understand why I placed a break point in my_kernel. By doing this it appears that fitrsvm is not taking advantage of vectorization. The arguments being passed into my_kernel by fitrsvm behave as if it is in a loop similar to what I have below:

%how my_kernel appears to be called in fitrsvm based on the arguments passed to my_kernel
%though I admit I can't actually see the code since it seems to dispatch to a C file
for i = 1:size(X,1):
    for j = 1:size(X,1):
        G(i,j) = feval('my_kernel',X(i,:), X(j,:))
    end
end

This same inefficient behavior also occurs when I predict the responses:

response = predict(m,Z)

I've managed to get a considerable speed up by replacing predict with the following:

response = m.Bias + m.Alpha' * feval('my_kernel',m.SupportVectors,Z);

Has anybody else noticed this behavior? Any suggestions? Am I missing a parameter maybe? I'm working on 2020a.

Answer 1

[Mario Andrés Muñoz]

I had this annoying problem too. I reached the conclusion that the best way to solve it is to use a 'precompiled' kernel and hack the kernel function as follows:

function g = precomp(u,v)

global K

g = K(u,v);

end

Where u and v correspond to the indexes of the observations. To achieve such indexing, I replaced the training data with the indexes as follows:

global K

K = myvectorizedkernelfcn(X,X);
Xidx = 1:size(X,1); % This corressponds to the number of rows
mdl = fitrsvm(Xidx',Y,'KernelFunction','precomp');

I hope this helps. Cheers.

[Mark Rucker]

That's actually a really good idea. I'd been "pre-computing" answers outside of the svm by only calling predict one time on all x, but doing it inside of the kernel call I think would be more flexible. Plus it is nice to know I'm not crazy and that somebody else has noticed this. Thanks!

[Mario Andrés Muñoz]

Not a problem. On further investigation, there are a couple of pointers more:

1. K must include both training and testing data, otherwise the indexing is not going to work.
2. It seems that for some reason, the code first checks whether the Kernel of (u,u) is always equal to 1 a vector at the time, after which it does run in parallel. You can further debug this function by changing it to this:

function g = precomp(u,v)

global K

g = K(u,v);
try
    disp([u v K(u,v)])
catch
    disp('ME')
end

end

Cheers

[Mark Rucker]

Regarding your first point. Thanks, yeah I understand that about the index approach.

Regarding your second point. You're right, I hadn't looked that closely at the values being passed in. After doing more testing tonight I also see what you say. It looks like it first it passes in every vector with itself. Once that is done it it then passes in every prediction vector one at a time to compute against the fitted vectors. So still not quite fully vectorized but better than vector by vector like I had originally thought.

[Mark Rucker]

Alright just finished all my testing and I take it back. In my case pre-computing my custom kernel is slower and takes more memory than pre-computing all solutions for my model. That is I'm doing this:

%assume that X is all possible values that could be found in my population
%assume that x,y are my training samples from my population X

m = fitrsvm(x,y,'KernelFunction','my_kernel');
Y = predict(m,X);

%this means that I can n;ow use x_indexes to directly lookup the solution
answers = Y(x_indexes)