所有可能组合的最快解决方案，在 k>2 和 n 大的情况下从 n 个可能的元素中取出 k 个元素
Fastest solution for all possible combinations, taking k elements out of n possible with k>2 and n large
我正在使用 MATLAB 从 n 个可能的元素中找出 k 个元素的所有可能组合。我偶然发现了 , but unfortunately it does not solve my problem. Of course, neither does nchoosek，因为我的 n 大约是 100。
事实是，我不需要同时使用所有可能的组合。我会解释我需要什么，因为可能有更简单的方法来达到预期的结果。我有一个 100 行和 25 列的矩阵 M。
将 M 的子矩阵视为由 M 的所有列和仅行的子集组成的矩阵。我有一个函数 f，它可以应用于任何结果为 -1 或 1 的矩阵。例如，您可以将函数视为 sign(det(A))，其中 A是任意矩阵（具体函数与这部分问题无关）。
我想知道 M 的最大行数是多少，由这些行组成的子矩阵 A 是这样的f(A) = 1。请注意，如果 f(M) = 1，我就完成了。但是，如果不是这种情况，那么我需要开始组合行，从所有具有 99 行的组合开始，然后采用具有 98 行的组合，依此类推。
到目前为止，我的实现与 nchoosek 有关，它在 M 只有几行时有效。然而，现在我正在处理一个相对更大的数据集，事情就卡住了。你们有没有想到一种无需使用上述功能即可实现此功能的方法？任何帮助将不胜感激。
这是我的最小工作示例，它适用于小 obs_tot 但当我尝试使用更大的数字时失败：
value = -1; obs_tot = 100; n_rows = 25;
mat = randi(obs_tot,n_rows);
while value == -1
posibles = nchoosek(1:obs_tot,i);
[num_tries,num_obs] = size(possibles);
num_try = 1;
while value == 0 && num_try <= num_tries
check = mat(possibles(num_try,:),:);
value = sign(det(check));
num_try = num_try + 1;
end
i = i - 1;
end
obs_used = possibles(num_try-1,:)';
前言
正如您在问题中注意到的那样，最好不要同时拥有 nchoosek 到 return 所有可能的组合，而是将它们一一枚举以免爆炸n变大时的内存。所以像：
enumerator = CombinationEnumerator(k, n);
while(enumerator.MoveNext())
currentCombination = enumerator.Current;
...
end
这是像 Matlab class 这样的枚举器的实现。它基于 C# / .NET 中的 classic IEnumerator<T> 接口，并模仿 nchoosek 中的子函数 combs (展开的方式):
%
% PURPOSE:
%
% Enumerates all combinations of length 'k' in a set of length 'n'.
%
% USAGE:
%
% enumerator = CombinaisonEnumerator(k, n);
% while(enumerator.MoveNext())
% currentCombination = enumerator.Current;
% ...
% end
%
%% ---
classdef CombinaisonEnumerator < handle
properties (Dependent) % NB: Matlab R2013b bug => Dependent must be declared before their get/set !
Current; % Gets the current element.
end
methods
function [enumerator] = CombinaisonEnumerator(k, n)
% Creates a new combinations enumerator.
if (~isscalar(n) || (n < 1) || (~isreal(n)) || (n ~= round(n))), error('`n` must be a scalar positive integer.'); end
if (~isscalar(k) || (k < 0) || (~isreal(k)) || (k ~= round(k))), error('`k` must be a scalar positive or null integer.'); end
if (k > n), error('`k` must be less or equal than `n`'); end
enumerator.k = k;
enumerator.n = n;
enumerator.v = 1:n;
enumerator.Reset();
end
function [b] = MoveNext(enumerator)
% Advances the enumerator to the next element of the collection.
if (~enumerator.isOkNext),
b = false; return;
end
if (enumerator.isInVoid)
if (enumerator.k == enumerator.n),
enumerator.isInVoid = false;
enumerator.current = enumerator.v;
elseif (enumerator.k == 1)
enumerator.isInVoid = false;
enumerator.index = 1;
enumerator.current = enumerator.v(enumerator.index);
else
enumerator.isInVoid = false;
enumerator.index = 1;
enumerator.recursion = CombinaisonEnumerator(enumerator.k - 1, enumerator.n - enumerator.index);
enumerator.recursion.v = enumerator.v((enumerator.index + 1):end); % adapt v (todo: should use private constructor)
enumerator.recursion.MoveNext();
enumerator.current = [enumerator.v(enumerator.index) enumerator.recursion.Current];
end
else
if (enumerator.k == enumerator.n),
enumerator.isInVoid = true;
enumerator.isOkNext = false;
elseif (enumerator.k == 1)
enumerator.index = enumerator.index + 1;
if (enumerator.index <= enumerator.n)
enumerator.current = enumerator.v(enumerator.index);
else
enumerator.isInVoid = true;
enumerator.isOkNext = false;
end
else
if (enumerator.recursion.MoveNext())
enumerator.current = [enumerator.v(enumerator.index) enumerator.recursion.Current];
else
enumerator.index = enumerator.index + 1;
if (enumerator.index <= (enumerator.n - enumerator.k + 1))
enumerator.recursion = CombinaisonEnumerator(enumerator.k - 1, enumerator.n - enumerator.index);
enumerator.recursion.v = enumerator.v((enumerator.index + 1):end); % adapt v (todo: should use private constructor)
enumerator.recursion.MoveNext();
enumerator.current = [enumerator.v(enumerator.index) enumerator.recursion.Current];
else
enumerator.isInVoid = true;
enumerator.isOkNext = false;
end
end
end
end
b = enumerator.isOkNext;
end
function [] = Reset(enumerator)
% Sets the enumerator to its initial position, which is before the first element.
enumerator.isInVoid = true;
enumerator.isOkNext = (enumerator.k > 0);
end
function [c] = get.Current(enumerator)
if (enumerator.isInVoid), error('Enumerator is positioned (before/after) the (first/last) element.'); end
c = enumerator.current;
end
end
properties (GetAccess=private, SetAccess=private)
k = [];
n = [];
v = [];
index = [];
recursion = [];
current = [];
isOkNext = false;
isInVoid = true;
end
end
我们可以通过命令 window 测试实现是否正常，如下所示：
>> e = CombinaisonEnumerator(3, 6);
>> while(e.MoveNext()), fprintf(1, '%s\n', num2str(e.Current)); end
return符合以下 n!/(k!*(n-k)!) 组合的预期：
1 2 3
1 2 4
1 2 5
1 2 6
1 3 4
1 3 5
1 3 6
1 4 5
1 4 6
1 5 6
2 3 4
2 3 5
2 3 6
2 4 5
2 4 6
2 5 6
3 4 5
3 4 6
3 5 6
4 5 6
可以进一步优化此枚举器的实现以提高速度，或者通过以更适合您的情况的顺序枚举组合（例如，先测试某些组合而不是其他组合）……好吧，至少它有效！ :)
问题解决
现在解决您的问题真的很简单：
n = 100;
m = 25;
matrix = rand(n, m);
k = n;
cont = true;
while(cont && (k >= 1))
e = CombinationEnumerator(k, n);
while(cont && e.MoveNext());
cont = f(matrix(e.Current(:), :)) ~= 1;
end
if (cont), k = k - 1; end
end
我正在使用 MATLAB 从 n 个可能的元素中找出 k 个元素的所有可能组合。我偶然发现了 nchoosek，因为我的 n 大约是 100。
事实是，我不需要同时使用所有可能的组合。我会解释我需要什么，因为可能有更简单的方法来达到预期的结果。我有一个 100 行和 25 列的矩阵 M。
将 M 的子矩阵视为由 M 的所有列和仅行的子集组成的矩阵。我有一个函数 f，它可以应用于任何结果为 -1 或 1 的矩阵。例如，您可以将函数视为 sign(det(A))，其中 A是任意矩阵（具体函数与这部分问题无关）。
我想知道 M 的最大行数是多少，由这些行组成的子矩阵 A 是这样的f(A) = 1。请注意，如果 f(M) = 1，我就完成了。但是，如果不是这种情况，那么我需要开始组合行，从所有具有 99 行的组合开始，然后采用具有 98 行的组合，依此类推。
到目前为止，我的实现与 nchoosek 有关，它在 M 只有几行时有效。然而，现在我正在处理一个相对更大的数据集，事情就卡住了。你们有没有想到一种无需使用上述功能即可实现此功能的方法？任何帮助将不胜感激。
这是我的最小工作示例，它适用于小 obs_tot 但当我尝试使用更大的数字时失败：
value = -1; obs_tot = 100; n_rows = 25;
mat = randi(obs_tot,n_rows);
while value == -1
posibles = nchoosek(1:obs_tot,i);
[num_tries,num_obs] = size(possibles);
num_try = 1;
while value == 0 && num_try <= num_tries
check = mat(possibles(num_try,:),:);
value = sign(det(check));
num_try = num_try + 1;
end
i = i - 1;
end
obs_used = possibles(num_try-1,:)';
前言
正如您在问题中注意到的那样，最好不要同时拥有 nchoosek 到 return 所有可能的组合，而是将它们一一枚举以免爆炸n变大时的内存。所以像：
enumerator = CombinationEnumerator(k, n);
while(enumerator.MoveNext())
currentCombination = enumerator.Current;
...
end
这是像 Matlab class 这样的枚举器的实现。它基于 C# / .NET 中的 classic IEnumerator<T> 接口，并模仿 nchoosek 中的子函数 combs (展开的方式):
%
% PURPOSE:
%
% Enumerates all combinations of length 'k' in a set of length 'n'.
%
% USAGE:
%
% enumerator = CombinaisonEnumerator(k, n);
% while(enumerator.MoveNext())
% currentCombination = enumerator.Current;
% ...
% end
%
%% ---
classdef CombinaisonEnumerator < handle
properties (Dependent) % NB: Matlab R2013b bug => Dependent must be declared before their get/set !
Current; % Gets the current element.
end
methods
function [enumerator] = CombinaisonEnumerator(k, n)
% Creates a new combinations enumerator.
if (~isscalar(n) || (n < 1) || (~isreal(n)) || (n ~= round(n))), error('`n` must be a scalar positive integer.'); end
if (~isscalar(k) || (k < 0) || (~isreal(k)) || (k ~= round(k))), error('`k` must be a scalar positive or null integer.'); end
if (k > n), error('`k` must be less or equal than `n`'); end
enumerator.k = k;
enumerator.n = n;
enumerator.v = 1:n;
enumerator.Reset();
end
function [b] = MoveNext(enumerator)
% Advances the enumerator to the next element of the collection.
if (~enumerator.isOkNext),
b = false; return;
end
if (enumerator.isInVoid)
if (enumerator.k == enumerator.n),
enumerator.isInVoid = false;
enumerator.current = enumerator.v;
elseif (enumerator.k == 1)
enumerator.isInVoid = false;
enumerator.index = 1;
enumerator.current = enumerator.v(enumerator.index);
else
enumerator.isInVoid = false;
enumerator.index = 1;
enumerator.recursion = CombinaisonEnumerator(enumerator.k - 1, enumerator.n - enumerator.index);
enumerator.recursion.v = enumerator.v((enumerator.index + 1):end); % adapt v (todo: should use private constructor)
enumerator.recursion.MoveNext();
enumerator.current = [enumerator.v(enumerator.index) enumerator.recursion.Current];
end
else
if (enumerator.k == enumerator.n),
enumerator.isInVoid = true;
enumerator.isOkNext = false;
elseif (enumerator.k == 1)
enumerator.index = enumerator.index + 1;
if (enumerator.index <= enumerator.n)
enumerator.current = enumerator.v(enumerator.index);
else
enumerator.isInVoid = true;
enumerator.isOkNext = false;
end
else
if (enumerator.recursion.MoveNext())
enumerator.current = [enumerator.v(enumerator.index) enumerator.recursion.Current];
else
enumerator.index = enumerator.index + 1;
if (enumerator.index <= (enumerator.n - enumerator.k + 1))
enumerator.recursion = CombinaisonEnumerator(enumerator.k - 1, enumerator.n - enumerator.index);
enumerator.recursion.v = enumerator.v((enumerator.index + 1):end); % adapt v (todo: should use private constructor)
enumerator.recursion.MoveNext();
enumerator.current = [enumerator.v(enumerator.index) enumerator.recursion.Current];
else
enumerator.isInVoid = true;
enumerator.isOkNext = false;
end
end
end
end
b = enumerator.isOkNext;
end
function [] = Reset(enumerator)
% Sets the enumerator to its initial position, which is before the first element.
enumerator.isInVoid = true;
enumerator.isOkNext = (enumerator.k > 0);
end
function [c] = get.Current(enumerator)
if (enumerator.isInVoid), error('Enumerator is positioned (before/after) the (first/last) element.'); end
c = enumerator.current;
end
end
properties (GetAccess=private, SetAccess=private)
k = [];
n = [];
v = [];
index = [];
recursion = [];
current = [];
isOkNext = false;
isInVoid = true;
end
end
我们可以通过命令 window 测试实现是否正常，如下所示：
>> e = CombinaisonEnumerator(3, 6);
>> while(e.MoveNext()), fprintf(1, '%s\n', num2str(e.Current)); end
return符合以下 n!/(k!*(n-k)!) 组合的预期：
1 2 3
1 2 4
1 2 5
1 2 6
1 3 4
1 3 5
1 3 6
1 4 5
1 4 6
1 5 6
2 3 4
2 3 5
2 3 6
2 4 5
2 4 6
2 5 6
3 4 5
3 4 6
3 5 6
4 5 6
可以进一步优化此枚举器的实现以提高速度，或者通过以更适合您的情况的顺序枚举组合（例如，先测试某些组合而不是其他组合）……好吧，至少它有效！ :)
问题解决
现在解决您的问题真的很简单：
n = 100;
m = 25;
matrix = rand(n, m);
k = n;
cont = true;
while(cont && (k >= 1))
e = CombinationEnumerator(k, n);
while(cont && e.MoveNext());
cont = f(matrix(e.Current(:), :)) ~= 1;
end
if (cont), k = k - 1; end
end