<p>Note: the original answer was completely redefined, because I misinterpreted the definition of the problem.</p>

<p>The evaluation of the kernel distance between <code>Xi</code> and <code>Xj</code> is presented below. Two codes that implement the algorithm are presented. The first code is inefficient but can be easily related to the definition of the kernel distance. The second code is much more efficient, but may not be as clear due to several vectorisation tricks.</p>

<p>The code assumes the following interpretation of the problem: </p>

<ol>
<li><code>Xi</code> and <code>Xj</code> are 2 data sets that contain 425 and 4 points, respectively. Each point belongs to <code>R^3</code> (real vector space with dimension 3).</li>
<li>The kernel distance between two data sets is calculated according to the definition given in the article by J.M. Phillips and S. Venkatasubramanian "A Gentle Introduction to the Kernel Distance" that can be found at the following <a href="https://arxiv.org/pdf/1103.1625v2.pdf" rel="nofollow noreferrer">link</a>. The definition is also provided below:</li>
</ol>

<hr>

<p><a href="https://i.stack.imgur.com/4TBRo.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/4TBRo.png" alt="enter image description here"></a></p>

<hr>

<p>The most straightforward implementation of the algorithm:</p>

<pre><code>% Initialisation.
clear;
clc;

% Construct Xi.
Xi = [randn(425, 1) randn(425, 1) randn(425, 1)];

% Definition of Xj.
Xj = [0.1 0.2 0.3; 0 0 0; -0.1 -0.1 -0.2; 1 -8 4];

% Convert to cell arrays.
Xi = mat2cell(Xi, ones(1, length(Xi(:, 1))), 3);
Xj = mat2cell(Xj, ones(1, length(Xj(:, 1))), 3);

% First, construct the kernel function for the evaluation of individual
% points in Xi and Xj
omega = 150;
a = 2;
kerFunction = @(xi, xj) exp(sum(abs(xi - xj).^a)/(omega^2));

kerDist = 0;
for i = 1 : length(Xj)
    for j = 1 : length(Xj)
        kerDist = kerDist + kerFunction(Xj{i}, Xj{j});
    end
end
for i = 1 : length(Xi)
    for j = 1 : length(Xi)
        kerDist = kerDist + kerFunction(Xi{i}, Xi{j});
    end
end
for i = 1 : length(Xi)
    for j = 1 : length(Xj)
        kerDist = kerDist - 2*kerFunction(Xi{i}, Xj{j});
    end
end
</code></pre>

<hr>

<p>A more efficient implementation of the algorithm is presented below:</p>

<pre><code>clear;

% Define constants.
omega = 150;
a = 2;

% Definition of Xi.
Xi = [randn(425, 1) randn(425, 1) randn(425, 1)];

% Definition of Xj.
Xj = [0.1 0.2 0.3; 0 0 0; -0.1 -0.1 -0.2; 1 -8 4];

% Definition of the characteristics of the data sets.
numPointsXj = length(Xj(:, 1));
numPointsXi = length(Xi(:, 1));

% Define a handle function for the definition of indices for the
% vectorisation of the kernel function.
hdlRepIdxPermutation = @(numPoints, numMatrixRep) ...
    repmat( ...
    (1 : numPoints : numPoints*(numMatrixRep - 1) + 1)', ...
    1, numPoints ...
    ) + ...
    repmat(0 : (numPoints - 1), numMatrixRep, 1);

tic

% Calculate the term that corresponds to K(p, p') in the definition of the
% kernal distance.
repXiRight = repmat(Xi, numPointsXi, 1);
leftIdxPermutationXi = hdlRepIdxPermutation(numPointsXi, numPointsXi);
repXiLeft = repXiRight(leftIdxPermutationXi(:), :);

kerDistComp1 = sum(exp(sum(abs(repXiLeft - repXiRight).^a, 2)/(omega^2)));

% Calculate the term that corresponds to K(q, q') in the definition of the
% kernal distance.
repXjRight = repmat(Xj, numPointsXj, 1);
leftIdxPermutationXj = hdlRepIdxPermutation(numPointsXj, numPointsXj);
repXjLeft = repXjRight(leftIdxPermutationXj(:), :);

kerDistComp2 = sum(exp(sum(abs(repXjLeft - repXjRight).^a, 2)/(omega^2)));

% Calculate the term that corresponds to K(p, q) in the definition of the
% kernal distance.
repXjRight = repmat(Xj, numPointsXi, 1);
repXiLeft = repmat(Xi, numPointsXj, 1);
leftIdxPermutationXi = hdlRepIdxPermutation(numPointsXi, numPointsXj);
repXiLeft = repXiLeft(leftIdxPermutationXi(:), :);
kerDistComp3 = -2*sum(exp(sum(abs(repXiLeft - repXjRight).^a, 2)/(omega^2)));

kerDist = kerDistComp1 + kerDistComp2 + kerDistComp3;

toc

disp(kerDist);
</code></pre>
Gaussian Kernel confirmation
<p>I'm trying to implement this kernel function<a href="https://i.stack.imgur.com/Uqzqu.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/Uqzqu.png" alt="enter image description here"></a></p>

<p>which is also known as radial basis function. Suppose that <code>a = 2</code>, <code>b = 1</code> and <code>σ = 150</code>. </p>

<ul>
<li>Xi is a 425x3 matrix</li>
<li>Xj is a 4x3 matrix</li>
</ul>

<p>I've came up with this code but I'm not sure that is correct. Can you help me?</p>

<pre><code>kS = exp( - (pdist2(Xj,Xi).^2) / (sigma^2) )
</code></pre>