vectorization - Generating arrays using bsxfun with anonymous function and for elementwise subtractions - MATLAB

Question

I have the following code:

n = 10000; 
s = 100;
Z = rand(n, 2);
x = rand(s, 1);
y = rand(s, 1);
fun = @(a) exp(a);

In principle, the anonymous function f can have a different form. I need to create two arrays.

First, I need to create an array of size n x s x s with generic elements

fun(Z(i, 1) - x(j)) * fun(Z(i, 2) - y(k))

where i=1,...n while j,k=1,...,s. What I can easily do, is to construct matrices using bsxfun, e.g.

bsxfun(@(x, y) fun(x - y), Z(:, 1), x');
bsxfun(@(x, y) fun(x - y), Z(:, 2), y');

But then I would need to combine them into 3D array by multiplying element-wise each column of those two matrices.

---

In the second step, I need to create an array of size n x 3 x s x s, which would look from one side as the following matrix

[ones(n, 1), Z(:, 1) - x(i), Z(:, 2) - y(j);]

where i=1,...s, j=1,...s. I could loop over the two extra dimensions with something like

A = [ones(n, 1), Z(:, 1) - x(1), Z(:, 2) - y(1)];
for i = 1:s
    for j = 1:s
        A(:, :, i, j) = [ones(n, 1), Z(:, 1) - x(i), Z(:, 2) - y(j);];
    end
end

Is there a way to avoid loops?

---

In the third step, suppose that after obtaining array out1 (output from first step), I want to create a new array out3 of dimension n x n x s x s, which contains the original array out1 on the main diagonal, i.e. out3(i,i,s,s) = out1(i, s, s) and out3(i,j,s,s)=0 for all i~=j. Is there some kind of alternative of diag for creating "diagonal arrays"? Alternatively, if I create n x n x s x s array of zeros, is there a way to put out1 on the main diagonal?

Answer 1

Code

exp_Z_x = exp(bsxfun(@minus,Z(:,1),x.')); %//'
exp_Z_y = exp(bsxfun(@minus,Z(:,2),y.')); %//'
out1 = bsxfun(@times,exp_Z_x,permute(exp_Z_y,[1 3 2]));

Z1 = [ones(n,1) Z(:,1) Z(:,2)];
X1 = permute([ zeros(s,1) x zeros(s,1)],[3 2 1]);
Y1 = permute([ zeros(s,1) zeros(s,1) y],[4 2 3 1]);
out2 = bsxfun(@minus,bsxfun(@minus,Z1,X1),Y1);

out3 = zeros(n,n,s,s); %// out3(n,n,s,s) = 0; could be used for performance
out3(bsxfun(@plus,[1:n+1:n*n]',[0:s*s-1]*n*n)) = out1; %//'

%// out1, out2 and out3 are the desired outputs