How to compute gaussian kernel matrix efficiently?

Hi,

I have a matrix X(10000, 800). I want to compute gramm matrix K(10000,10000), where K(i,j)= exp(-(X(i,:)-X(j,:))^2).

First i used double for loop, but then it just hangs forever. Then I tried this:

[N d] = size(X); aa = repmat(X',[1 N]); bb = repmat(reshape(X',1,[]),[N 1]); K = reshape((aa-bb).^2, [N*N d]); K = reshape(sum(D,2),[N N]); But then it uses a lot of extra space and I run out of memory very soon. Is there any efficient vectorized method for this. I am sure there must be something as this is quite a standard intermediate step for many kernel svms and also in image processing.

[Matt J]

Your expression for K(i,j) does not evaluate to a scalar. Are you sure you don't want something like

exp(-norm( X(i,:) - X(j,:) ))^2);

Also, please format your code so it's more readable.

[Shahid Mahmood]

can you explain the whole procedure in detail to compute a kernel matrix in matlab

Answer 1

[Matt J]

Assuming you really want exp(-norm( X(i,:) - X(j,:) ))^2), then one way is

nsq=sum(X.^2,2);

K=bsxfun(@minus,nsq,(2*X)*X.');
K=bsxfun(@plus,nsq.',K);
K=exp(-K);

[Jeon]

How I can modify the code when I want to involve 'sigma', that is, I want to calculate 'exp(-norm(X1(:,i)-X2(:,j))^2/(2*sigma^2));' instead? Thanks

[Matt J]

Just pre-normalize X1 and X2 by sigma.

[Farzan Zaheer]

I am working on Kernel LMS, and I am having issues with the implementation of Kernel. I want to know what exactly is "X2" here. I am implementing the Kernel using recursion.

I am using the following statement,

for n=2:K-1
     Kernel(n)=exp(-0.5*(dist(x(:,2:n),x(:,n)')/ker_bw^2));
end

where ker_bw is the kernel bandwidth/sigma and x is input of (1000,1) and I have reshaped the input x as

x = [x(1:end-1),x(2:end)];

as mentioned in the research paper I am following. Any help will be highly appreciated.

Answer 2

[Image Analyst]

If you have the Image Processing Toolbox, why not use fspecial()?

h = fspecial('gaussian', hsize, sigma) returns a rotationally symmetric Gaussian lowpass filter of size hsize with standard deviation sigma (positive). hsize can be a vector specifying the number of rows and columns in h, or it can be a scalar, in which case h is a square matrix. The default value for hsize is [3 3]; the default value for sigma is 0.5.

[Mohammad Alwardat]

Hi Image Analyst,

How can I study the similarity between 2 vectors x and y using Gaussian kernel similarity algorithm? And how can I determine the parameter sigma?

Thanks♥

[amel kaouane]

i have the same problem, don't know to get the parameter sigma

[Image Analyst]

@amel kaouane it comes from your mind. You think up some sigma that might work, assign it like

sigma = 5;
hsize = 21;
h = fspecial('gaussian', hsize, sigma)

h = 21×21

0.0001 0.0002 0.0003 0.0003 0.0005 0.0006 0.0007 0.0008 0.0009 0.0009 0.0009 0.0009 0.0009 0.0008 0.0007 0.0006 0.0005 0.0003 0.0003 0.0002 0.0001 0.0002 0.0003 0.0004 0.0005 0.0007 0.0008 0.0010 0.0011 0.0012 0.0013 0.0014 0.0013 0.0012 0.0011 0.0010 0.0008 0.0007 0.0005 0.0004 0.0003 0.0002 0.0003 0.0004 0.0005 0.0007 0.0009 0.0012 0.0014 0.0016 0.0018 0.0019 0.0019 0.0019 0.0018 0.0016 0.0014 0.0012 0.0009 0.0007 0.0005 0.0004 0.0003 0.0003 0.0005 0.0007 0.0010 0.0012 0.0016 0.0019 0.0021 0.0024 0.0025 0.0026 0.0025 0.0024 0.0021 0.0019 0.0016 0.0012 0.0010 0.0007 0.0005 0.0003 0.0005 0.0007 0.0009 0.0012 0.0016 0.0020 0.0024 0.0028 0.0031 0.0033 0.0033 0.0033 0.0031 0.0028 0.0024 0.0020 0.0016 0.0012 0.0009 0.0007 0.0005 0.0006 0.0008 0.0012 0.0016 0.0020 0.0025 0.0030 0.0035 0.0038 0.0041 0.0042 0.0041 0.0038 0.0035 0.0030 0.0025 0.0020 0.0016 0.0012 0.0008 0.0006 0.0007 0.0010 0.0014 0.0019 0.0024 0.0030 0.0036 0.0042 0.0046 0.0049 0.0050 0.0049 0.0046 0.0042 0.0036 0.0030 0.0024 0.0019 0.0014 0.0010 0.0007 0.0008 0.0011 0.0016 0.0021 0.0028 0.0035 0.0042 0.0048 0.0053 0.0056 0.0057 0.0056 0.0053 0.0048 0.0042 0.0035 0.0028 0.0021 0.0016 0.0011 0.0008 0.0009 0.0012 0.0018 0.0024 0.0031 0.0038 0.0046 0.0053 0.0058 0.0062 0.0063 0.0062 0.0058 0.0053 0.0046 0.0038 0.0031 0.0024 0.0018 0.0012 0.0009 0.0009 0.0013 0.0019 0.0025 0.0033 0.0041 0.0049 0.0056 0.0062 0.0066 0.0067 0.0066 0.0062 0.0056 0.0049 0.0041 0.0033 0.0025 0.0019 0.0013 0.0009

imshow(h, []);
axis('on', 'image')

If you don't like 5 for sigma then just try others until you get one that you like. It's not like I can tell you the perfect value of sigma because it really depends on your situation and image.

[amel kaouane]

am looking to get similarity between two time series by using this gaussian kernel, i think it's not the same situation, right?!