这个问题是由非常具体的组合优化问题推动的，其中搜索空间被定义为具有多重性的向量未排序离散值集的置换子集的空间。
我正在寻找有效（足够快，矢量化或任何其他更聪明的解决方案）功能，它能够以下列方式找到子集的索引：
t = [1 1 3 2 2 2 3 ]
是所有可能值的未排序向量，包括其多重性。
item = [2 3 1; 2 1 2; 3 1 1; 1 3 3]
是矢量t的置换子集的列表。
我需要找到对应于向量t的子集项的相应索引列表。所以，对于上面提到的例子，我们有：
item =
2 3 1
2 1 2
3 1 1
1 3 3
t =
1 1 3 2 2 2 3
ind = item2ind(item,t)
ind =
4 3 1
4 1 5
3 1 2
1 3 7
因此，对于item = [2 3 1]，我们得到ind = [4 3 1]，这意味着：
项目的第一个值“2”对应于位置“4”的t处的第一个值“2”，项目的第二个值“3”对应于位置“3”的t处的第一个值“3”和第三个值“1” “at item对应于位于”1“的t处的第一个值”1“。
在案例项= [2 1 2]中，我们得到ind = [4 1 5]，这意味着：
项目的第一个值“2”对应于位置“4”的t处的第一个值“2”，项目的第二个值“1”对应于位置“1”的t处的第一个值“1”，第三个值“ 2“在项目对应于位于”5“的t处的第二个（!!!）值”1“。
对于
item = [1 1 1]
不存在任何解决方案，因为向量t只包含两个“1”。
我当前版本的函数“item2ind”是非常简单的串行代码，通过将“for”改为“parfor”循环可以简单地并行化：
function ind = item2ind(item,t)
[nlp,N] = size(item);
ind = zeros(nlp,N);
for i = 1:nlp
auxitem = item(i,:);
auxt = t;
for j = 1:N
I = find(auxitem(j) == auxt,1,'first');
if ~isempty(I)
auxt(I) = 0;
ind(i,j) = I;
else
error('Incompatible content of item and t.');
end
end
end
end
但我需要一些更聪明的东西......而且更快:)
更大输入数据的测试用例：
t = 1:10; % 10 unique values at vector t
t = repmat(t,1,5); % unsorted vector t with multiplicity of all unique values 5
nlp = 100000; % number of item rows
[~,p] = sort(rand(nlp,length(t)),2); % 100000 random permutations
item = t(p); % transform permutations to items
item = item(:,1:30); % transform item to shorter subset
tic;ind = item2ind(item,t);toc % runing and timing of the original function
tic;ind_ = item2ind_new(item,t);toc % runing and timing of the new function
isequal(ind,ind_) % comparison of solutions
为了实现代码的矢量化，我假设不会出现错误情况。它应该首先丢弃，我将在下面介绍一个简单的过程。
方法首先，让我们计算t中所有元素的索引：
t = t(:);
mct = max(accumarray(t,1));
G = accumarray(t,1:length(t),[],@(x) {sort(x)});
G = cellfun(@(x) padarray(x.',[0 mct-length(x)],0,'post'), G, 'UniformOutput', false);
G = vertcat(G{:});
说明：将输入放入列向量形状后，我们使用t计算accumarray中每个可能值的最大出现次数。现在，我们形成所有数字的所有索引的数组。它形成一个单元阵列，因为每个值的出现次数可能不同。为了形成矩阵，我们将每个数组独立地填充到最大长度（命名为mct）。然后我们可以将单元阵列转换为矩阵。在这一步，我们有：
G =
1 11 21 31 41
2 12 22 32 42
3 13 23 33 43
4 14 24 34 44
5 15 25 35 45
6 16 26 36 46
7 17 27 37 47
8 18 28 38 48
9 19 29 39 49
10 20 30 40 50
现在，我们处理item。为此，让我们弄清楚如何在向量内创建值的累积和。例如，如果我有：
A = [1 1 3 2 2 2 3];
然后我想得到：
B = [1 2 1 1 2 3 2];
由于隐式扩展，我们可以将它放在一行：
B = diag(cumsum(A==A'));
就这么简单。语法A==A'扩展为矩阵，其中每个元素都是A(i)==A(j)。仅在一个维度上生成累积和并且取对角线给出了良好的结果：每个列出现在一个值上的累积总和。
要使用2-D的item这个技巧，我们应该使用3D数组。我们打电话给m=size(item,1)和n=size(item,2)。所以：
C = cumsum(reshape(item,m,1,n)==item,3);
所有累积出现的（大）3D矩阵。最后一件事是沿着尺寸2和3选择对角线上的列：
ia = C(sub2ind(size(C),repelem((1:m).',1,n),repelem(1:n,m,1),repelem(1:n,m,1)));
现在，使用所有这些矩阵，索引很容易：
ind = G(sub2ind(size(G),item,ia));
最后，让我们回顾一下函数的代码：
function ind = item2ind_new(item,t)
t = t(:);
[m,n] = size(item);
mct = max(accumarray(t,1));
G = accumarray(t,1:length(t),[],@(x) {sort(x)});
G = cellfun(@(x) padarray(x.',[0 mct-length(x)],0,'post'), G, 'UniformOutput', false);
G = vertcat(G{:});
C = cumsum(reshape(item,m,1,n)==item,3);
ia = C(sub2ind(size(C),repelem((1:m).',1,n),repelem(1:n,m,1),repelem(1:n,m,1)));
ind = G(sub2ind(size(G),item,ia));
结果在旧的4核上运行提供的脚本，我得到：
Elapsed time is 4.317914 seconds.
Elapsed time is 0.556803 seconds.
ans =
logical
1
加速是次要的（超过8倍），以及内存消耗（使用矩阵C）。我想这个部分可以做一些改进以节省更多内存。
编辑为了生成ia，此过程可能会导致内存丢失。节省内存的一种方法是使用for循环直接生成此数组：
ia = zeros(size(item));
for i=unique(t(:)).'
ia = ia+cumsum(item==i, 2).*(item==i);
end
在所有情况下，当你有ia时，很容易测试item与t相比是否有错误：
any(ind(:)==0)
然后，获得错误项目（作为掩码）的简单解决方案
min(ind,[],2)==0