How to compute the sum of squares of outer products of two matrices minus a common matrix in Matlab?

You can use:

X2 = X.^2;
Y2 = Y.^2;
S2 = S.^2;
b = sum(sum(X2.' * Y2 - 2 * (X.' * Y ) .* S + n * S2));

Given your example

b=0;
for i = 1:n
   b = b  + sum(sum((X(i,:).' * Y(i,:) - S).^2));
end

We can first bring the summation out of the loop:

b=0;
for i = 1:n
  b = b  + (X(i,:).' * Y(i,:) - S).^2;
end
b=sum(b(:))

Knowing that we can write (a - b)^2 as a^2 - 2*a*b + b^2

b=0;
for i = 1:n
  b = b  + (X(i,:).' * Y(i,:)).^2 - 2.* (X(i,:).' * Y(i,:)) .*S + S.^2;
end
b=sum(b(:))

And we know that (a * b) ^ 2 is the same as a^2 * b^2:

X2 = X.^2;
Y2 = Y.^2;
S2 = S.^2;
b=0;
for i = 1:n
  b = b  + (X2(i,:).' * Y2(i,:)) - 2.* (X(i,:).' * Y(i,:)) .*S + S2;
end
b=sum(b(:))

Now we can compute each term separately:

 b = sum(sum(X2.' * Y2 - 2 * (X.' * Y ) .* S + n * S2));

Here is the result of a test in Octave that compares my method and two other methods provided by @AndrasDeak and the original loop based solution for inputs of size 500*500:

===rahnema1 (B)===
Elapsed time is 0.0984299 seconds.

===Andras Deak (B2)===
Elapsed time is 7.86407 seconds.

===Andras Deak (B3)===
Elapsed time is 2.99158 seconds.

===Loop solution===
Elapsed time is 2.20357 seconds


n=500;
X= rand(n);
Y= rand(n);
S= rand(n);

disp('===rahnema1 (B)===')
tic
    X2 = X.^2;
    Y2 = Y.^2;
    S2 = S.^2;
    b=sum(sum(X2.' * Y2 - 2 * (X.' * Y ) .* S + n * S2));
toc
disp('===Andras Deak (B2)===')
tic
    b2 = sum(reshape((permute(reshape(X, [n, 1, n]).*Y, [3,2,1]) - S).^2, 1, []));
toc
disp('===Andras Deak (B3)===')
tic
    b3 = sum(reshape((reshape(X, [n, 1, n]).*Y - reshape(S.', [1, n, n])).^2, 1, []));
toc
tic
    b=0;
    for i = 1:n
      b = b  + sum(sum((X(i,:)' * Y(i,:) - S).^2));
    end
toc

You probably can't spare time complexity, but you can make use of vectorization to get rid of the loop and make use of low-level code and caching as much as possible. Whether it's actually faster depends on your dimensions, so you need to do some timing tests to see if this is worth it:

% dummy data
n = 3;
X = rand(n);
Y = rand(n);
S = rand(n);

% vectorize
b2 = sum(reshape((permute(reshape(X, [n, 1, n]).*Y, [3,2,1]) - S).^2, 1, []));

% check
b - b2 % close to machine epsilon i.e. zero

What happens is that we insert a new singleton dimension in one of the arrays, ending up with an array of size [n, 1, n] against one with [n, n], the latter being implicitly the same as [n, n, 1]. The overlapping first index corresponds to the i in your loop, the remaining two indices correspond to the matrix indices of the dyadic product you have for each i. Then we permute the indices in order to put the "i" index last, so that we can again broadcast the result with S of (implicit) size [n, n, 1]. What we then have is a matrix of size [n, n, n] where the first two indices are matrix indices in your original and the last one corresponds to i. We then just have to take the square and sum each term (instead of summing twice I reshaped the array to a row and summed once).

A slight variation of the above transposes S instead of the 3d array which might be faster (again, you should time it):

b3 = sum(reshape((reshape(X, [n, 1, n]).*Y - reshape(S.', [1, n, n])).^2, 1, []));

In terms of performance, reshape is free (it only reinterprets data, it doesn't copy) but permute/transpose will often lead to a perforance hit when data gets copied.