2

我有以下代码：

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

原则上，匿名函数f可以有不同的形式。我需要创建两个数组。

首先，我需要n x s x s使用通用元素创建一个大小数组

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

在哪里i=1,...n同时j,k=1,...,s。我可以轻松做的是使用 构造矩阵bsxfun，例如

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

但是然后我需要3D通过将这两个矩阵的每一列逐元素相乘来将它们组合成数组。

在第二步中，我需要创建一个 size 的数组n x 3 x s x s，它从一侧看起来像下面的矩阵

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

其中i=1,...s, j=1,...s. 我可以用类似的东西循环两个额外的维度

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

有没有办法避免循环？

在第三步中，假设在获得数组out1（从第一步的输出）之后，我想创建一个out3维度为 的新数组n x n x s x s，其中包含主对角线上的原始数组out1，即out3(i,i,s,s) = out1(i, s, s)and out3(i,j,s,s)=0for all i~=j。是否有某种替代方法diag来创建“对角线阵列”？或者，如果我创建n x n x s x s零数组，有没有办法放在out1主对角线上？

4

1 回答 1

2

代码

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
于 2014-10-27T17:33:37.067 回答