向量化 octave/matlab 代码

Question

以下是八度编码（kmeans 的一部分）

centroidSum = zeros(K);
valueSum = zeros(K, n);
for i = 1 : m  
  for j = 1 : K    
    if(idx(i) == j)
      centroidSum(j) = centroidSum(j) + 1;
      valueSum(j, :) = valueSum(j, :) + X(i, :);
    end
  end
end

代码有效，是否可以对代码进行矢量化？ 不用 if 语句很容易对代码进行矢量化， 但是我们如何用 if 语句对代码进行向量化呢？

Answer 1

我假设代码的目的是计算n 维空间中一组m 数据点的子集的质心，其中点存储在矩阵X 中（点x 坐标) 并且向量 idx 为每个数据点指定该点所属的子集 (1 ... K)。那么一个部分向量化就是：

centroid = zeros(K, n)
for j = 1 : K
   centroid(j, :) = mean(X(idx == j, :));
end

if 被索引消除，特别是逻辑索引：idx == j 给出一个布尔数组，指示哪些数据点属于子集j。

我认为也有可能摆脱第二个 for 循环，但这会导致代码非常复杂、难以理解。

Answer 2

简介及解决方案代码

这可能是一种基于 -

的完全矢量化方法

• accumarray：用于累积和计算 valueSum。这还介绍了一种如何使用 accumarray on a 2D matrix along a certain direction 的技术，而这在直接的方式中是不可能的。
• bsxfun：用于计算所有列的线性索引，以匹配来自idx 的行索引。

这是实现 -

%// Store no. of columns in X for frequent usage later on
ncols = size(X,2); 

%// Find indices in idx that are within [1:k] range, call them as labels
%// Also, find their locations in that range array, call those as pos
[pos,id] = ismember(idx,1:K);
labels = id(pos);
%// OR with bsxfun: [pos,labels] = find(bsxfun(@eq,idx(:),1:K));

%// Find all labels, i.e. across all columns of X
all_labels = bsxfun(@plus,labels(:),[0:ncols-1]*K);

%// Get truncated X corresponding to all indices matches across all columns
X_cut = X(pos,:);

%// Accumulate summations within each column based on the labels.
%// Note that accumarray doesn't accept matrices, so we were required
%// to create all_labels that had same labels within each column and
%// offsetted at constant intervals from consecutive columns
acc1 = accumarray(all_labels(:),X_cut(:));

%// Regularise accumulated array and reshape back to a 2D array version
acc1_reg2D = [acc1 ; zeros(K*ncols - numel(acc1),1)]; 
valueSum = reshape(acc1_reg2D,[],ncols);
centroidSum = histc(labels,1:K); %// Get labels counts as centroid sums

基准代码

%// Datasize parameters
K = 5000;
n = 5000;
m = 5000;

idx = randi(9,1,m);
X = rand(m,n);

disp('----------------------------- With Original Approach')
tic
centroidSum1 = zeros(K,1);
valueSum1 = zeros(K, n);
for i = 1 : m  
  for j = 1 : K    
    if(idx(i) == j)
      centroidSum1(j) = centroidSum1(j) + 1;
      valueSum1(j, :) = valueSum1(j, :) + X(i, :);
    end
  end
end
toc, clear valueSum1 centroidSum1

disp('----------------------------- With Proposed Approach')
tic
%// ... Code from earlied mentioned section
toc

运行时结果

----------------------------- With Original Approach
Elapsed time is 1.235412 seconds.
----------------------------- With Proposed Approach
Elapsed time is 0.379133 seconds.

Answer 3

不确定它的运行时性能，但这是一个非卷积向量化实现：

b = idx == 1:K;
centroids = (b' * X) ./ sum(b)';

Answer 4

向量化计算会在性能上产生巨大差异。基准测试

1. 原代码，
2. 来自 A. Donda 的部分矢量化和
3. 来自 Tom 的完整矢量化，

给了我以下结果：

Original Code: Elapsed time is 1.327877 seconds.

Partial Vectorization: Elapsed time is 0.630767 seconds.

Full Vectorization: Elapsed time is 0.021129 seconds.

这里的基准代码：

%// Datasize parameters
K = 5000;
n = 5000;
m = 5000;

idx = randi(9,1,m);
X = rand(m,n);

fprintf('\nOriginal Code: ')
tic
centroidSum1 = zeros(K,1);
valueSum1 = zeros(K, n);
for i = 1 : m  
  for j = 1 : K    
    if(idx(i) == j)
      centroidSum1(j) = centroidSum1(j) + 1;
      valueSum1(j, :) = valueSum1(j, :) + X(i, :);
    end
  end
end
centroids = valueSum1 ./ centroidSum1;
toc, clear valueSum1 centroidSum1 centroids

fprintf('\nPartial Vectorization: ')
tic
centroids = zeros(K,n);
for k = 1:K
    centroids(k,:) = mean( X(idx == k, :) );
end
toc, clear centroids

fprintf('\nFull Vectorization: ')
tic
centroids = zeros(K,n);
b = idx == 1:K;
centroids = (b * X) ./ sum(b)';
toc

注意，我在原始代码中添加了额外的一行，以元素方式将 valueSum1 除以 centroidSum1，以使每种代码的输出相同。

最后，我知道这不是严格意义上的“答案”，但是我没有足够的声誉来添加评论，而且我认为基准测试数据对任何正在学习 MATLAB（如我自己）并需要掌握矢量化的一些额外动力。