如何避免循环以减少此代码的计算时间？

Question

如何避免循环以减少此代码（我上一个问题的一个解决方案）的计算时间：

我希望找到A(1:3,:)的列向量，其在M(4,:)中对应的值不是单元格X的其中一个向量的一部分（显然也不等于这些向量之一）。如果X非常大，我希望找到一种快速的解决方案。

M = [1007  1007  4044  1007  4044  1007  5002 5002 5002 622 622;
      552   552   300   552   300   552   431  431  431 124 124; 
     2010  2010  1113  2010  1113  2010  1100 1100 1100  88  88;
        7    12    25    15    12    30     2   10   55  32  12];

这里直接取 A：

A = [1007  4044  5002  622;
      552   300   431  124;
     2010  1113  1100   88];

A 包含 M（1:3，：） 的唯一列向量

X = {[2 5 68 44],[2 10 55 9 17],[1 55 6 7 8 9],[32 12]};

[~, ~, subs] = unique(M(1:3,:)','rows');

A4 = accumarray(subs(:),M(4,:).',[],@(x) {x});

%// getting a mask of which columns we want
idxC(length(A4)) = false;
for ii = 1:length(A4)
    idxC(ii) = ~any(cellfun(@(x) all(ismember(A4{ii},x)), X));
end

显示我们想要的列

out = A(:,idxC)

结果：

>> out

out =

    1007        4044
     552         300
    2010        1113

由于[2;10;55]包含在X{2} = [2 10 55 9 17]中，因此列向量[5002;431;1100]被排除。

由于[32 12] = X{4}，因此列向量[622;124;88]被排除。

另一个例子：使用相同的X

    M = [1007  4044  1007  4044  1007  5002 5002 5002 622 622  1007  1007  1007;
          552   300   552   300   552   431  431  431 124 124   552    11    11; 
         2010  1113  2010  1113  2010  1100 1100 1100  88  88  2010    20    20;
           12    25    15    12    30     2   10   55  32  12     7    12     7];

X = {[2 5 68 44],[2 10 55 9 17],[1 55 6 7 8 9],[32 12]};

A = [1007  4044  5002  622  1077;
      552   300   431  124    11;
     2010  1113  1100   88    20];

结果: (带有scmg答案)

如果按照第一行对 A 进行排序，我得到了正确的结果：

out =

         1007        1007        4044
           11         552         300
           20        2010        1113

如果我不对矩阵A进行排序，结果会是错误的：（错误的结果）

out =

        4044        5002         622
         300         431         124
        1113        1100          88

应该消除列向量A(:,4) = [622;124;88]，因为[32 12] = X{4}。

应该消除列向量[5002;431;1100]，因为[2;10;55]包含在X{2} = [2 10 55 9 17]中。

- bzak

---

你可以解释一下你得出输出结果的逻辑吗？这将节省我们从代码中推导出逻辑所需的时间。 - Luis Mendo

@LuisMendo：我收到了两个关于我的问题的答案。scmg的回复给出了正确的输出，就像例子中一样，但如果X非常大，则需要大量计算时间。Ben Voigt开发的逻辑很有趣，但输出结果是错误的，我找不出原因！我的问题中输入是M、A和X，输出是out = A（：，idxC）。 - bzak

@LuisMendo：我希望找到A(1:3,:)的列向量，其在M(4,:)中对应的值不是单元格X中任意一个向量的一部分（并且显然不等于这些向量中的任何一个）。如果X非常大，我希望快速解决此问题。 - bzak

只是澄清一下：您的意思是“在M（4，:)中对应的值不属于单元格X的同一向量”，对吗？ - Luis Mendo

@LuisMendo：是的，就是单元格X的同一向量。 - bzak

也许在你的问题中加入这个解释，这样你就可以得到更多的帮助。 - Luis Mendo

Answer 1

在这种情况下，你不应该试图消除循环。向量化实际上会严重影响性能。

特别地，给你的匿名lambda函数取一个名称。

issubset = @(x) all(ismember(A4{ii},x))

这种写法非常低效，因为它没有短路。请使用循环代替。

同样适用于：

any(cellfun(issubset, X))

请使用类似于以下的方法：

idxC = true(size(A4));
NX = numel(X);
for ii = 1:length(A4)
    for jj = 1:NX
        xj = X{jj};
        issubset = true;
        for A4i=A4{ii}
            if ~ismember(A4i, xj)
                issubset = false;
                break;
            end;
        end;
        if issubset
            idxC(ii) = false;
            break;
        end;
    end;
end;

两个break语句，尤其是第二个，触发了一个提前退出，可能为您节省了大量的计算。

- Ben Voigt

---

是的，那些应该是括号而不是大括号。 - Ben Voigt

你的回答很快，但是它给出了错误的结果，我认为你的代码有误。 - bzak

“for A4i=A4{ii}” 是什么意思？循环在哪里进行？我找不到你代码中的错误！ - bzak

@bzak：有两个快捷方式。第一个是，如果 A4{ii} 中的任何元素不在 X{jj} 中找到，则不要测试 A4{ii} 的其余部分，重新开始下一个 jj。另外，如果 A4{ii} 的所有元素都在任何 X{jj} 中找到，则不要测试 jj 的其余值，已经删除该 A4{ii}。 - Ben Voigt

但是为什么你的代码没有给出预期的结果，而另一个答案却给出了正确的结果呢？ - bzak

Answer 2

本·沃特的回答非常好，但是for A4i = A4{ii}这一行可能会导致问题：对于列向量来说，这种方式的for循环是不起作用的。

%row vector
for i = 1:3
    disp('foo');
end

    foo
    foo
    foo

%column vector
for i = (1:3).'
    disp('foo');
end

    foo

尝试使用A4i = A4{ii}.'，它可以完成您的工作！

现在，如果我们查看输出：

A(:,idxC) =

    4044        5002
     300         431
    1113        1100

正如您所见，最终结果并非我们期望的那样。

只要unique执行了某种排序，子字符串就不是按照 A 中遇到的顺序编号的，而是按照在已排序的 C 中遇到的顺序编号的:

subs =

 2
 2
 3
 2
 3
 2
 4
 4
 4
 1
 1

因此，您应该通过由unique给出的矩阵而不是A来获取最终输出。

输入

[C, ~, subs] = unique(M(1:3,:)','rows'); 
%% rather than [~, ~, subs] = unique(M(1:3,:)','rows');

然后，要获得最终输出，请输入

>> out = C(idxC,:).'
out =

        1007        4044
         552         300
        2010        1113

- Ikaros

---

我批准这个澄清的答案。如何使循环在行向量和列向量上工作？也许像A4{ii}(:)'这样的东西？ - Ben Voigt

转置操作是 .' 而不是 '！仅使用 ' 会进行复共轭转置，这是一种不同的操作并会导致错误的结果。 - hbaderts

你尝试过在大数据上使用A(:,idxC)这种方法吗？因为即使是scmg代码在这个例子上也给出了错误的答案，让我觉得你的例子有问题。 - Ikaros

@HamtaroWarrior：我在我的大数据上尝试了A(:,idxC)和C(idxC,:)两种方法，但发现结果错误！而且奇怪的是，只有scmg代码能够给我正确的结果，无论是在示例还是真实数据中！ - bzak

好的，我会去搜索它。 - Ikaros

Answer 3

第一步

这一节中列出了一种被认为是快速且直接的解决我们问题的方法。请注意，由于A是从M中考虑前三行唯一列的矩阵，因此在此处跳过作为输入，因为我们将使用解决方案代码在内部生成它。这也保留在下一个方法/步骤中。以下是实现方式：

function out = shot1_func(M,X)

%// Get unique columns and corresponding subscripts
[unqrows, ~, subs_idx] = unique(M(1:3,:)','rows');
unqcols = unqrows.'; %//'

counts = accumarray(subs_idx(:),1); %// Counts of each unique subs_idx

%// Modify each cell of X based on their relevance with the fourth row of M
X1 = cellfun(@(x) subs_idx(ismember(M(4,:),x)),X,'Uni',0);

lensX = cellfun('length',X1); %// Cell element count of X1

Xn = vertcat(X1{:}); %// Numeric array version of X
N = max(subs_idx);   %// Number of unique subs_idx

%// Finally, get decision mask to select the correst columns from unqcols
sums = cumsum(bsxfun(@eq,Xn,1:N),1);
cumsums_at_shifts = sums(cumsum(lensX),:);

mask1 = any(bsxfun(@eq,diff(cumsums_at_shifts,[],1),counts(:).'),1); %//'
decision_mask = mask1 | cumsums_at_shifts(1,:) == counts(:).';    %//'
out = unqcols(:,~decision_mask);

return

---

第二步

前面提到的方法可能会在以下位置出现瓶颈：

cellfun(@(x) subs_idx(ismember(M4,x)),X,'Uni',0)

因此，为了保持良好的性能，可以将整个过程分成两个阶段。第一阶段可以处理M的第四行中未重复的X单元格，这可以使用向量化方法实现；另一个阶段解决其余X单元格，使用较慢的cellfun方法。

因此，代码会变得有点臃肿，但希望能够具有更好的性能。最终的实现将类似于：

%// Get unique columns and corresponding subscripts
[unqrows, ~, subs_idx] = unique(M(1:3,:)','rows')
unqcols = unqrows.' %//'
counts = accumarray(subs_idx,1);

%// Form ID array for X
lX = cellfun('length',X)
X_id = zeros(1,sum(lX))
X_id([1 cumsum(lX(1:end-1)) + 1]) = 1
X_id = cumsum(X_id)

Xr = cellfun(@(x) x(:).',X,'Uni',0); %//'# Convert to cells of row vectors
X1 = [Xr{:}]                         %// Get numeric array version

%// Detect cells that are to be processed by part1 (vectorized code)
[valid,idx1] = ismember(M(4,:),X1)
p1v = ~ismember(1:max(X_id),unique(X_id(accumarray(idx1(valid).',1)>1))) %//'

X_part1 = Xr(p1v)
X_part2 = Xr(~p1v)

%// Get decision masks from first and second passes and thus the final output
N = size(unqcols,2);
dm1 = first_pass(X_part1,M(4,:),subs_idx,counts,N)
dm2 = second_pass(X_part2,M(4,:),subs_idx,counts)
out = unqcols(:,~dm1 & ~dm2)

相关函数 -

function decision_mask = first_pass(X,M4,subs_idx,counts,N)

lensX = cellfun('length',X)'; %//'# Get X cells lengths
X1 = [X{:}];                  %// Extract cell data from X

%// Finally, get the decision mask
vals = changem(X1,subs_idx,M4) .* ismember(X1,M4);

sums = cumsum(bsxfun(@eq,vals(:),1:N),1);
cumsums_at_shifts = sums(cumsum(lensX),:);
mask1 = any(bsxfun(@eq,diff(cumsums_at_shifts,[],1),counts(:).'),1); %//'
decision_mask = mask1 | cumsums_at_shifts(1,:) == counts(:).';    %//'
return


function decision_mask = second_pass(X,M4,subs_idx,counts)

%// Modify each cell of X based on their relevance with the fourth row of M
X1 = cellfun(@(x) subs_idx(ismember(M4,x)),X,'Uni',0);

lensX = cellfun('length',X1); %// Cell element count of X1

Xn = vertcat(X1{:}); %// Numeric array version of X
N = max(subs_idx);   %// Number of unique subs_idx

%// Finally, get decision mask to select the correst columns from unqcols
sums = cumsum(bsxfun(@eq,Xn,1:N),1);
cumsums_at_shifts = sums(cumsum(lensX),:);

mask1 = any(bsxfun(@eq,diff(cumsums_at_shifts,[],1),counts(:).'),1); %//'
decision_mask = mask1 | cumsums_at_shifts(1,:) == counts(:).';       %//'

return

---

验证

本部分列出了用于验证输出的代码。以下是验证shot＃1代码的代码-

%// Setup inputs and output
load('matrice_data.mat');   %// Load input data
X = cellfun(@(x) unique(x).',X,'Uni',0); %// Consider X's unique elements
out = shot1_func(M,X); %// output with Shot#1 function

%// Accumulate fourth row data from M based on the uniqueness from first 3 rows
[unqrows, ~, subs] = unique(M(1:3,:)','rows');    %//'
unqcols = unqrows.';                              %//'
M4 = accumarray(subs(:),M(4,:).',[],@(x) {x});    %//'
M4 = cellfun(@(x) unique(x),M4,'Uni',0);

%// Find out cells in M4 that correspond to unique columns unqcols
[unqcols_idx,~] = find(pdist2(unqcols.',out.')==0);

%// Finally, verify output
for ii = 1:numel(unqcols_idx)
    for jj = 1:numel(X)
        if all(ismember(M4{unqcols_idx(ii)},X{jj}))
            error('Error: Wrong output!')
        end
    end
end
disp('Success!')

- Divakar

---

让我们在聊天室中继续此讨论。 - Divakar

我将真实的M和X数据保存在一个文件中，然后直接将您的答案应用于这些数据（输入：M和X），得到的输出正好是矩阵A。但是，对于示例数据，我得到了正确的结果，我觉得这非常奇怪！！！ - bzak

@bzak 请查看修改内容，看看是否适用于您的实际情况？ - Divakar

非常感谢Divakar为帮助我所付出的努力。应用于示例的代码给出了正确的结果，但是对于我的真实数据，我仍然得到错误的结果！这里有没有一种方法可以将我的真实数据发送给您？ - bzak

你找到了与matrice_Correct_result相同的结果吗？ - bzak

Answer 4

也许你可以使用2次cellfun:

idxC = cellfun(@(a) ~any(cellfun(@(x) all(ismember(a,x)), X)), A4, 'un', 0);
idxC = cell2mat(idxC);
out = A(:,idxC)

- scmg

---

是的，起初您应该避免使用cellfun，而像其他答案一样使用for循环，直到一切都正确为止，然后才能逐部分尝试矢量化。我只想指出，您可以连续使用2次cellfun，但正确性取决于您的实际问题，在其中必须自行进行调整。 - scmg