从可能的n个元素中取出k(其中k>2且n较大)个元素的最快组合解决方案

3
我正在使用MATLAB查找n个可能元素中的所有k个元素的所有可能组合。我偶然发现了this question,但不幸的是它没有解决我的问题。当然,nchoosek也不能解决我的问题,因为我的n大约是100。 事实上,我并不需要同时得到所有可能的组合。我将解释我需要什么,因为可能有更简单的方法来实现所需的结果。我有一个由100行和25列组成的矩阵M。 将M的子矩阵视为由M的所有列和只有一部分行组成的矩阵。我有一个可以应用于任何矩阵的函数f,它给出-1或1的结果。例如,您可以将该函数视为sign(det(A)),其中A是任何矩阵(精确的函数对此问题的这一部分无关紧要)。 我想知道在子矩阵 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,:)';

1
所以,你说你在处理更多行的nchoosek时遇到了一些问题。我猜你和我们一样明白,这是因为组合数变得更大了。然而,我想在这里问一下,你预计需要测试多少个“i”?我尝试估算了一些最坏情况下的组合数,比如从100行中选择50行。可能的组合数可以计算为100!/(50!* 50!),这给出了约1.0 * 10 ^ 29个组合的最大数量。从1到1.0 * 10 ^ 29的循环会非常耗费时间。 - patrik
是的,我明白那是个问题。我之前以为大概10-20次迭代就足够了,但我无法确定。现在我在思考一个很好的估算方法。例如,抽样一次观察并检查条件是否成立。然后再抽样另一次观察,将其添加并再次检查,如果条件不符合则再抽样另一次观察并重试,如此进行多次而不会导致过程崩溃。 - MathUser
йӮЈд№Ҳиҝҷж„Ҹе‘ізқҖжӮЁдёҚжғід»Һ100дёӘеҸҜиғҪзҡ„иЎҢдёӯиҺ·еҸ–жүҖжңүеҸҜиғҪзҡ„kиЎҢз»„еҗҲпјҢжҳҜиҝҷж ·еҗ—пјҹиҝҳжҳҜжҲ‘зҺ°еңЁзҗҶи§Јй”ҷдәҶпјҹ - patrik
1
不要对获得nchoosek的真正替代品抱有太大希望。虽然有一些比nchoosek更快的实现,但是由此产生的矩阵会非常大,以至于这并不重要。如果无法在数学上简化,那么您应该使用启发式方法。我假设您已经验证了您的问题不能通过简单的贪心方法解决?为什么不从大小为k=1的组合开始随机生成,并且每找到一个f(A)=1就增加k呢?(例如:comb=randperm(n,k)。)这样,您至少有一点机会找到一个大的k - knedlsepp
@scmg 在第4段中所写的内容由于计算复杂度很快就会变得无法解决。即使在k=10时,问题已经有了10^13种可能的组合,这意味着遍历所有组合不是一个选项。然而,OP表明他有一些方法可以减少组合数。但是,除非我知道如何做到这一点,否则很难给出适当的答案。 - patrik
显示剩余2条评论
1个回答

1

前言

正如您在问题中所注意到的那样,最好不要让 nchoosek 一次返回所有可能的组合,而是按顺序逐个枚举它们,以避免当 n 变大时内存溢出。因此,类似于:

enumerator = CombinationEnumerator(k, n);
while(enumerator.MoveNext())
    currentCombination = enumerator.Current;
    ...
end
这是一个Matlab类的枚举器实现,与C# / .NET中经典的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
我们可以从命令窗口测试实现是否正常,方法如下:
>> e = CombinaisonEnumerator(3, 6);
>> while(e.MoveNext()), fprintf(1, '%s\n', num2str(e.Current)); end

这将按预期返回以下 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 

这是一个绝对惊人的解决方案,适用于更大的N,我很惊讶它没有更多的赞。不过有一个问题。我正在尝试使用您自己的示例来实现此功能。在您的问题解决部分中,行cont = f(matrix(e.Current(:), :) ~= 1;缺少一个括号,导致我无法运行代码。我尝试在基本上所有地方添加它,但没有任何地方有效。我不确定在这种情况下f()是什么意思,也不想破坏它。您能澄清一下这行代码应该做什么或者括号放在哪里吗(我很惭愧地请求这个澄清)。 - chainhomelow
@chainhomelow 感谢您的点赞和提供有关缺少括号的信息(我已编辑我的帖子以修复此问题)。该示例仅为伪代码,以展示如何在特定情况下使用我的枚举器。函数“f”未定义(您必须将其替换为您想要的任何函数)。@MathUser 在他的问题中只指定它返回“-1”或“1”,并且当它返回“1”时,代码应停止。 - CitizenInsane
啊,谢谢你的澄清。我以前没有使用过用户构建的类,担心它与其结构有关。感谢你的时间。 - chainhomelow

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,