6
  • z - 双精度矩阵,大小 Nx2;
  • x - 双精度矩阵,大小 Nx2;

sup = x(i, :);

phi(1, i) = {@(z) exp(-g * sum((z - sup(ones([size(z, 1) 1]),:)) .^ 2, 2))};

这是逻辑回归的径向基函数 (RBF)。这是公式:

在此处输入图像描述

我需要你的建议,我可以优化这个公式吗?因为它调用了数百万次,而且需要很多时间......

4

2 回答 2

6

在您最近的编辑中,您似乎引入了一些语法错误,但我想我理解您想要做什么(从第一个版本开始)。

考虑使用高效的BSXFUN函数,而不是使用REPMAT或索引来重复向量x(i,:)以匹配 的行:z

rbf(:,i) = exp( -g .* sum(bsxfun(@minus,z,x(i,:)).^2,2) );

上面显然循环了 x 的每一行


您可以更进一步,并使用PDIST2z计算和中每对行之间的欧几里德距离x

%# some random data
X = rand(10,2);
Z = rand(10,2);
g = 0.5;

%# one-line solution
rbf = exp(-g .* pdist2(Z,X,'euclidean').^2);

现在矩阵中的每个值:rbf(i,j)对应于和之间的函数z(i,:)x(j,:)


编辑:

我对不同的方法进行了计时,这是我使用的代码:

%# some random data
N = 5000;
X = rand(N,2);
Z = rand(N,2);
g = 0.5;

%# PDIST2
tic
rbf1 = exp(-g .* pdist2(Z,X,'euclidean').^2);
toc

%# BSXFUN+loop
tic
rbf2 = zeros(N,N);
for j=1:N
    rbf2(:,j) = exp( -g .* sum(bsxfun(@minus,Z,X(j,:)).^2,2) );
end
toc

%# REPMAT+loop
tic
rbf3 = zeros(N,N);
for j=1:N
    rbf3(:,j) = exp( -g .* sum((Z-repmat(X(j,:),[N 1])).^2,2) );
end
toc

%# check if results are equal
all( abs(rbf1(:)-rbf2(:)) < 1e-15 )
all( abs(rbf2(:)-rbf3(:)) < 1e-15 )

结果:

Elapsed time is 2.108313 seconds.     # PDIST2
Elapsed time is 1.975865 seconds.     # BSXFUN
Elapsed time is 2.706201 seconds.     # REPMAT
于 2011-08-08T23:46:34.713 回答
3

Amro 提到了一些非常好的方法。但是可以通过重塑其中一个矩阵来进一步利用 bsxfun。

>> type r.m

N = 5000;
X = rand(N,2);
Z = rand(N,2);
g = 0.5;

%BSXFUN+loop
tic
rbf2 = zeros(N,N);
for j=1:N
    rbf2(:,j) = exp( -g .* sum(bsxfun(@minus,Z,X(j,:)).^2,2) );
end
toc

tic
diffs = bsxfun(@minus, reshape(X', [1, 2, N]), Z);
dist = reshape(sum(diffs.^2, 2), [N, N]);
rbf3 = exp(-g .* dist);
toc

>> r
Elapsed time is 2.235527 seconds.
Elapsed time is 0.877833 seconds.
>> r
Elapsed time is 2.253943 seconds.
Elapsed time is 1.047295 seconds.
>> r
Elapsed time is 2.234132 seconds.
Elapsed time is 0.856302 seconds.
>> max(abs(rbf2(:) - rbf3(:)))

ans =

     0

您想从 Z 的每一行中减去 X 的每一行。当其中一个是向量而另一个是矩阵时,这通常是直截了当的。但是如果它们都是矩阵,我们可以通过确保卷中的每个矩阵只包含一个向量来做到这一点。这里我选择了 X,但是 Z 可以和 X 互换使用。

于 2011-08-09T03:43:44.640 回答