1

我有一个矩阵X(10000, 800)。我想计算 gramm 矩阵K(10000,10000)，其中K(i,j)= exp(-(X(i,:)-X(j,:))^2).

首先我使用双循环，但后来它永远挂起。然后我尝试了这个：

[N d] = size(X);
aa = repmat(X',[1 N]);
bb = repmat(reshape(X',1,[]),[N 1]);
K = reshape((aa-bb).^2, [N*N d]);
K = reshape(sum(D,2),[N N]);

但是后来它使用了很多额外的空间，我很快就用完了内存。是否有任何有效的矢量化方法。我确信一定有一些东西，因为这对于许多内核 svm 以及图像处理来说都是一个相当标准的中间步骤。

4

2 回答 2

3

使用pdist2或 pdist。请注意，来自 Matlab 的 pdist2 很慢...... Matlab余弦距离很慢

代码：

X = rand(100, 3);
K = squareform(pdist(X, 'euclidean'));
K = exp(-K.^2);

我会为更一般的情况写这个，你有两个矩阵并且你想找到所有的距离。 (x-y)^2 = x'x - 2x'y + y'y如果要计算 Gram 矩阵，则需要所有差异组合。

X = rand(100, 3);
Y = rand(50, 3);
A = sum(X .* X, 2);
B = -2 *X * Y';
C = sum(Y .* Y, 2);
K = bsxfun(@plus, A, B);
K = bsxfun(@plus, K, C);
K = exp(-K);

编辑：速度比较

代码

% http://stackoverflow.com/questions/13109826/compute-a-gramm-matrix-in-matlab-without-loops/24407122#24407122
function time_gramm()
% I have a matrix X(10000, 800). I want to compute gramm matrix K(10000,10000), where K(i,j)= exp(-(X(i,:)-X(j,:))^2).
X = rand(100, 800);

%% The straight-forward pdist solution.
tic;
K = squareform(pdist(X, 'euclidean'));
K1 = exp(-K .^2);
t1 = toc;
fprintf('pdist took \t%d seconds\n', t1);

%% The vectorized solution
tic;
A = sum(X .* X, 2);
B = -2 * X * X';
K = bsxfun(@plus, A, B);
K = bsxfun(@plus, K, A');
K2 = exp(-K);
t2 = toc;
fprintf('Vectorized solution took \t%d seconds.\n', t2);

%% The not-so-smart triple-loop solution
tic;
N = size(X, 1);
K3 = zeros(N, N);
for i=1:N
    %     fprintf('Running outer loop for i= %d\n', i);
    for j=1:N
        xij = X(i,:) - X(j,:);
        xij = norm(xij, 2);
        xij = xij ^ 2;
        K3(i,j) = -xij;
        %         d = X(i,:) - X(j,:); % Alternative way, twice as fast but still
        %         orders of magnitude slower than the other solutions.
        %         K3(i,j) = exp(-d * d');
    end
end
K3 = exp(K3);
t3 = toc;
fprintf('Triple-nested loop took \t%d seconds\n', t3);
%% Assert results are the same...
assert(all(abs(K1(:) - K2(:)) < 1e-6 ));
assert(all(abs(K1(:) - K3(:)) < 1e-6 ));
end

结果

我用 N=100 运行了上面的代码

pdist took  8.600000e-03 seconds
Vectorized solution took    3.916000e-03 seconds.
Triple-nested loop took     2.699330e-01 seconds

请注意，在问题的第 100 个请求大小处，另一个答案 ( O(m^2 n)) 中建议的代码的性能要慢两个数量级。到那时，我插入了 100k 作为X矩阵的大小，这比我愿意等待的时间要长得多。

在全尺寸问题 ( X = rand(10000, 800)) 上的表现是这样的：

pdist took  5.470632e+01 seconds
Vectorized solution took    1.141894e+01 seconds.

评论

矢量化解决方案耗时 11 秒，Matlab 的 pdist 耗时 55 秒，而另一个示例中建议的手动解决方案从未完成。

于 2014-06-25T11:13:26.887 回答
0

为什么不使用简单的公式呢？对于元素K(i, j) = exp( sum_{k=0}^{n} (X(i, k) - X(j, k))^2。所以，这是两个外循环 fori和j一个内循环 for k。时间复杂度为O(m^2 n)，其中有m行和n列X。空间复杂度是因为除了和矩阵O(1)之外，您不使用更多空间来计算答案。XK

你试过这个吗，真的这么慢吗？

于 2012-10-28T15:10:43.053 回答