是否可以优化此 Matlab 代码以使用 k-means 的质心进行矢量量化？答案

【问题标题】：Is it possible to optimize this Matlab code for doing vector quantization with centroids from k-means?是否可以优化此 Matlab 代码以使用 k-means 的质心进行矢量量化？
【发布时间】：2011-08-09 20:21:32
【问题描述】：

我使用大小为 4000x300（4000 个质心，每个质心有 300 个特征）的 k-means 创建了一个码本。使用密码本，然后我想标记一个输入向量（用于稍后进行分箱）。输入向量的大小为 Nx300，其中 N 是我收到的输入实例的总数。

为了计算标签，我为每个输入向量计算最近的质心。为此，我将每个输入向量与所有质心进行比较，并选择距离最小的质心。那么标签就是那个质心的索引。

我当前的 Matlab 代码如下所示：

function labels = assign_labels(centroids, X)
labels = zeros(size(X, 1), 1);

% for each X, calculate the distance from each centroid
for i = 1:size(X, 1)
    % distance of X_i from all j centroids is: sum((X_i - centroid_j)^2)
    % note: we leave off the sqrt as an optimization
    distances = sum(bsxfun(@minus, centroids, X(i, :)) .^ 2, 2);
    [value, label] = min(distances);
    labels(i) = label;
end

但是，这段代码仍然相当慢（出于我的目的），我希望有办法进一步优化代码。

一个明显的问题是有一个 for 循环，它是 Matlab 良好性能的祸根。我一直在尝试想出一种方法来摆脱它，但没有运气（我研究过将 arrayfun 与 bsxfun 结合使用，但还没有让它起作用）。或者，如果有人知道任何其他方法可以加快速度，我将不胜感激。

更新

经过一番搜索，我找不到使用 Matlab 的好解决方案，因此我决定查看 Python 的 scikits.learn 包中用于 'euclidean_distance'（缩短）的内容：

 XX = sum(X * X, axis=1)[:, newaxis]
 YY = Y.copy()
 YY **= 2
 YY = sum(YY, axis=1)[newaxis, :]
 distances = XX + YY
 distances -= 2 * dot(X, Y.T)
 distances = maximum(distances, 0)

它使用欧几里得距离 ((x-y)^2 -> x^2 + y^2 - 2xy) 的二项式形式，据我所知，它通常运行得更快。我完全未经测试的 Matlab 翻译是：

 XX = sum(data .* data, 2);
 YY = sum(center .^ 2, 2);
 [val, ~] = max(XX + YY - 2*data*center');

【问题讨论】：

相关：pdist2 equivalent in MATLAB version 7

标签： optimization matlab vector k-means quantization

【解决方案1】：

使用以下函数计算您的距离。你应该会看到一个数量级的加速

两个矩阵 A 和 B 将列作为维度，将行作为每个点。 A 是您的质心矩阵。 B 是您的数据点矩阵。

function D=getSim(A,B)
    Qa=repmat(dot(A,A,2),1,size(B,1));
    Qb=repmat(dot(B,B,2),1,size(A,1));
    D=Qa+Qb'-2*A*B';

【讨论】：

+1 给你。这确实更快。想起来，我不知道我为什么不这样做，而不是使用cellfun。我最近使用cellfun 的方式太多了，应该将它从我的工具包中退役，至少暂时:)

【解决方案2】：

您可以通过转换为单元格并使用cellfun来对其进行矢量化：

[nRows,nCols]=size(X);
XCell=num2cell(X,2);
dist=reshape(cell2mat(cellfun(@(x)(sum(bsxfun(@minus,centroids,x).^2,2)),XCell,'UniformOutput',false)),nRows,nRows);
[~,labels]=min(dist);

说明：

我们将X 的每一行分配给第二行中自己的单元格
这块@(x)(sum(bsxfun(@minus,centroids,x).^2,2))是一个匿名函数，和你的distances=...行一样，使用cell2mat，我们将它应用到X的每一行。
标签就是每列最小行的索引。

【讨论】：

看起来 num2cell(X) 将 X 的每个元素转换为自己的向量。但是，对于 (x-centroids)，我们希望将 X 的每一行与每个质心相减。那么它应该改为：XCell=num2cell(X, 2)？
@Abe，你是对的。虽然它肯定会返回相同的答案，但在将其应用于每个元素时会产生不必要的函数调用开销。我已经解决了。但请记住bsxfun 和cellfun 通常只是编写循环的简洁方式，并且不一定要更快（有时是，但并非总是如此）。为与您相同维度的矩阵计时您的循环和 cellfun 代码，它们甚至在 93 秒左右时相当不错，仅相差十分之一。
太好了，谢谢！是的，我希望 Matlab 的内部循环可能比他们的一般 for 循环更好，但我也没有看到很大的改进。此外，for 循环允许“parfor”。最初想使用 GPU，但 bsxfun 不允许匿名函数，因此无法传递“质心”。我确实设法找到的一项改进是其他人在 SO 上的帖子循环了质心而不是数据。我的数据中有 4000 个质心和 1e5 到 1e6 个特征向量。所以通过在质心上循环，我可以用矩阵数学加快速度。
我之前的评论中没有空格，但我使用了这里的标签分配代码：stackoverflow.com/questions/1373516/matlabk-means-clustering/…

【解决方案3】：

对于真正的矩阵实现，您可以考虑尝试以下方式：

  P2 = kron(centroids, ones(size(X,1),1));
  Q2 = kron(ones(size(centroids,1),1), X);

  distances = reshape(sum((Q2-P2).^2,2), size(X,1), size(centroids,1));

注意这假设数据组织为 [x1 y1 ...; x2 y2 ...;...]

【讨论】：

这会占用大量内存——任何真正的矢量化版本都会如此。 @yoda 已经提供了（较慢的）替代方案。你也可以用 repmat 替换 kron - 这只是我已经在我自己的一个项目中使用的一个版本
嗯。我要试试这个，因为我目前的方法需要> 4天。不过，我对 kron 并不熟悉，所以我将尝试使用上面的 repmats 重写它。如果有机会，您介意验证一下吗？
没关系。在纸上，看起来如果我将两个矩阵都设为 3 维（重复数据矩阵 size(centroids, 1) 次，并为每个质心重复每个质心 size(data, 1) 次，我可以减去，平方，总和，然后对它们进行最小化。但是，我想不出确保第二个矩阵的好方法。
当然 - 你能告诉我N 有多大吗？正如我所提到的，内存使用量会很快变大——考虑到 300 维的测量空间，我认为它最多只能处理几百个左右的 N。
好的，很高兴知道。我现在正在阅读克朗，但我以前从未见过它以这种方式使用过。至于N，它相当大，所以这可能是一个问题。我的数据一般是 100,000，我的质心在 4,000 左右。我可以更改数据的大小，因为我正在分页，但如果它太小，我担心会产生高成本。

【解决方案4】：

您可以使用比蛮力更有效的最近邻搜索算法。最流行的方法是 Kd-Tree。 O(log(n)) 平均查询时间而不是 O(n) 蛮力复杂度。关于 Kd-Trees 的 Maltab 实现，可以看看here

【讨论】：