Why this formula is written with this Matlab code?

Question

In section 4 of this paper, the following formula is presented to zero-center the kernel matrix K of size p:

And this is the code I have corresponding to the above formula:

K = K - (1/p)*(K*ones(p,1))*ones(1,p) - (1/p)*ones(p,1)*(ones(1,p)*K) + (1/p^2)*sum(sum(K));

I'm confused about the relation of the code to the actual formula in the paper. In particular about the last two members - (1/p)*ones(p,1)*(ones(1,p)*K) and (1/p^2)*sum(sum(K)).

Can someone please explain it?

Answer 1

Well, the code is not completely correct. The third term includes a K which is not there in the presented formula from the paper. Additionally, the last term does not multiply with e e^T. However, the latter can be omitted in this case since MATLAB will automatically add the scalar to all the elements in the matrix. The same applies for the third term, so it can be omitted there as well.
Here is a correct version of the line with the above mentioned simplification:

K = K - (1/p)*(K*ones(p,1))*ones(1,p) - 1/p + (1/p^2)*sum(sum(K))

---

We can simplify further by making only one call to ones, since ones(p,1)*ones(1,p) gives you the same result as ones(p). Furthermore sum(sum(K)) can be replaced with sum(K(:)).
It would look like this:

K = K - (1/p)*K*ones(p) - 1/p + (1/p^2)*sum(K(:))

---

Now we can compare this to what would be the one-to-one implementation of the formula. Therefore, we will use e = ones(p,1) to represent e. To get e^T, you can just transpose e with .'. So the formula can be written as follows:

K = K - (1/p)*K*e*e.' - (1/p)*e*e.' + ((e.'*K*e)/p^2)*e*e.'

Note that e.'*K*e just calculates the sum of all elements in K, which is equal to sum(K(:)). This is valid because e = ones(p,1).

Let's generate some sample data and compare the results:

rng(8);             % make it reproducible
p = 3;              % size of matrix
K = randi(10,p);    % generate random matrix
e = ones(p,1);      % generate e-vector

K1 = K - (1/p)*(K*ones(p,1))*ones(1,p) - 1/p + (1/p^2)*sum(sum(K))
K2 = K - (1/p)*K*ones(p) - 1/p + (1/p^2)*sum(K(:))
K3 = K - (1/p)*K*e*e.' - (1/p)*e*e.' + ((e.'*K*e)/p^2)*e*e.'
K4 = K - (1/p)*K*e*e.' - (1/p)*e*e.' + sum(K(:))/p^2

Here is the result:

K1 =
    8.0000    5.0000    4.0000
    9.6667    2.6667    4.6667
    9.3333    1.3333    6.3333
K2 =
    8.0000    5.0000    4.0000
    9.6667    2.6667    4.6667
    9.3333    1.3333    6.3333
K3 =
    8.0000    5.0000    4.0000
    9.6667    2.6667    4.6667
    9.3333    1.3333    6.3333
K4 =
    8.0000    5.0000    4.0000
    9.6667    2.6667    4.6667
    9.3333    1.3333    6.3333