加快矩阵输入

Question

我的 matlab 计算速度有点问题。我能够在 matlab 中编写代码来运行小矩阵的计算，但它使用嵌套的 for 循环，并且对于我正在使用的大型数据集，matlab 无法完成计算。

注意：我对 Matlab 不是很熟悉，所以虽然程序可以运行，但效率极低。

简而言之，我正在尝试创建一个矩阵，其条目描述一组独特位置之间的关系。作为一个具体的例子，我们从这个矩阵开始：

B = 

5873 4 1

5873 7 1

5873 1 1

2819 8 2

2819 1 2

9771 4 3

9771 2 3

9771 5 3

9771 6 3

5548 7 4

其中第三列是唯一的位置标识符，第二列是恰好在该位置中的“段”的编号。请注意，多个段可以落入不同的位置。

我想要创建一个描述不同位置之间关系的矩阵。具体来说，对于位置 i 和 j，我希望新矩阵的 (i,j) 条目是 i 和 j 共有的段数除以 i 和 j 组合的段总数。

目前我的代码如下所示：

C = zeros(max(B.data(:,3)), max(B.data(:,3)));

for i = 1:max(B.data(:,3))

  for j = 1:max(B.data(:,3))

    vi = B.data(:,3) == i;
    vj = B.data(:,3) == j;
    C(i,j) = numel(intersect(B.data(vi,2), B.data(vj,2))) / numel(union(B.data(vi,2), B.data(vj,2)));
  end

end

但它非常非常慢。有人对加快计算有什么建议吗？

非常感谢！！

Answer 1

分组和循环的方法

1. 位置根据段进行分组（单元格数组groups 在下面的代码中）。带有自定义函数的accumarray 用于该任务。
2. 大小等于最大位置的方阵C 被初始化为零。对于每个组，属于该组的所有位置对在C 中的条目都增加了1。
3. 矩阵C 已归一化。

代码：

groups = accumarray(B.data(:,2), B.data(:,3), [], @(x) {x});  %// step 1
C = zeros(max(B.data(:,3)));                                  %// step 2
for n = 1:numel(groups);
    ind = groups{n};
    C(ind,ind) = C(ind,ind)+1;
end
d = diag(C);                                                  %// step 3
C = C./(bsxfun(@plus,d,d.')-C);

使用 3D 数组的矢量化方法

这种方法需要大量内存；并且根据下面@Divakar 的 cmets，对于非常大的数据集，它并不比循环方法快：

1. 它构建一个 3D logical 数组 T 使得 T(m,n,s) 是 1 当且仅当位置 m 和 n 共享段 s。这可以通过bsxfun 高效完成。
2. 将T 沿第三个维度求和，得到的C 与之前的方法相同。
3. 矩阵C 的归一化方式与之前相同。

代码：

T = full(sparse(B.data(:,3), B.data(:,2), 1));               %// step 1
T = bsxfun(@and, permute(T, [1 3 2]), permute(T, [3 1 2]));
C = sum(T, 3);                                               %// step 2
d = diag(C);                                                 %// step 3
C = C./(bsxfun(@plus,d,d.')-C);

---

基准测试

感谢@Divakar 提供以下基准测试和解释：

这里有一些运行时比较 loop-based 方法和后者 vectorized 之一 -

***** Parameters: No. of rows in B ~= 10000 and No. of groups = 10 *****
---------------------------------- With loopy approach
Elapsed time is 0.242635 seconds.
---------------------------------- With vectorized approach
Elapsed time is 0.022174 seconds.
 
 
***** Parameters: No. of rows in B ~= 10000 and No. of groups = 100 *****
---------------------------------- With loopy approach
Elapsed time is 0.318268 seconds.
---------------------------------- With vectorized approach
Elapsed time is 0.451242 seconds.
 
 
***** Parameters: No. of rows in B ~= 100000 and No. of groups = 100 *****
---------------------------------- With loopy approach
Elapsed time is 1.173182 seconds.
---------------------------------- With vectorized approach
Elapsed time is 0.464339 seconds.
 
 
***** Parameters: No. of rows in B ~= 100000 and No. of groups = 1000 *****
---------------------------------- With loopy approach
Elapsed time is 10.310780 seconds.
---------------------------------- With vectorized approach
Elapsed time is 54.216923 seconds.

可以注意到，对于 B 中的相同行数，组数的增加意味着矢量化方法的性能会显着下降。

因此，在为问题案例选择适当的方法时，请牢记这一点。

Answer 2

这可能是一种使用单循环的方法 -

%// Store column-2 and 3 values in separate variables for easy access
col2 = B(:,2);
col3 = B(:,3);

%// Start and end indices of the groups w.r.t column-3 of B
ends = [find(diff(col3)) ; numel(col3)];
starts = [1 ; ends(1:end-1)+1];

ngrp = max(col3); %// number of groups
intersectM = zeros(ngrp); %// initialize 2D matrix to store intersect counts
for k = 1:ngrp-1

    %// Matches for each group with respect to each other group 
    %// but one less after each iteration
    matches = ismember(col2(starts(k+1):end),col2(starts(k):ends(k)));
    %// OR matches =
    %// any(bsxfun(@eq,col2(starts(k):ends(k)),col2(starts(k+1):end).'),1)'

    %// Store intersect counts
    intersectM(k,k+1:end) = accumarray(col3(starts(k+1):end)-k,matches,[]);
end

grp_counts = histc(col3,1:ngrp); %// counts of group elements

%// Peform numel(intersect)/numel(union)
C = intersectM./(bsxfun(@plus,grp_counts,grp_counts') - intersectM); %//'

%// Since the output is diagonal symmetric, so just copy over the upper
%// diagonal elements into the lower diagonal places
C = triu(C)' + C; %//'
C(1:ngrp+1:end) = 1; %// Set diagonal elements as ones