Generando arreglos usando bsxfun con función anónima y para restas de elementos-MATLAB

vectorization

Question

Tengo el siguiente código:

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

En principio, la función anónima f puede tener una forma diferente. Necesito crear dos matrices.

Primero, necesito crear una matriz de tamaño nxsxs con elementos genéricos

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

donde i=1,...n mientras j,k=1,...,s . Lo que puedo hacer fácilmente, es construir matrices usando bsxfun , por ejemplo

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

Pero luego necesitaría combinarlos en una matriz 3D multiplicando elementos por cada columna de esas dos matrices.

En el segundo paso, necesito crear una matriz de tamaño nx 3 xsxs , que se vería desde un lado como la siguiente matriz

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

donde i=1,...s , j=1,...s . Podría recorrer las dos dimensiones extra con algo como

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

¿Hay alguna manera de evitar los bucles?

En el tercer paso, supongamos que después de obtener array out1 (salida desde el primer paso), quiero crear una nueva matriz out3 de dimensión nxnxsxs , que contiene la matriz original out1 en la diagonal principal, es decir, out3(i,i,s,s) = out1(i, s, s) y out3(i,j,s,s)=0 para todo i~=j . ¿Hay algún tipo de alternativa de diag para crear "arreglos diagonales"? De forma alternativa, si creo nxnxsxs array de ceros, ¿hay alguna manera de poner out1 en la diagonal principal?

Answer 1

Código

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