快速且近乎“完美”的划分问题解决方案

• 下载 Java 源代码 - 54.7 KB
• 下载 C# 源代码 - 150.4 KB
• 下载 C++ 源代码 - 7.9 MB

目录

• 引言
• 算法
  • 贪心算法术语
  • 贪心算法
  • 动态搜索树
• 最终排列
• 实现
• 测试
• 结论
• 致谢

简介

一般而言，划分问题是指判断给定一个包含 N 个正整数的集合是否可以被划分为 k 个子集，使得每个子集中数字的和相等。在现实生活中，子集和的绝对相等通常是无法实现的，因此我们将采用划分问题的优化版本，即划分多重集为 k 个子集，使得子集和之间的差值最小化。一个典型的应用场景是将一组具有不同执行时间的独立工作分配给多个并发线程。如果每个线程运行在独立的 CPU 核心上，我们希望确保每个 CPU 节点的负载均匀。

划分问题是 NP 完全问题，而优化版本是 NP 难问题。在实际应用中，这意味着当 N 大于约 20 且 k 大于约 3 时，完全枚举方法将变得不堪重负。

尽管有实际应用中可以高效解决优化版本的说法，但在我工作中遇到这个问题时，却找不到任何可靠的现成解决方案。这迫使我不得不自己寻找解决方案。

算法

术语

我们将通过子集中数字的最大和与最小和之间的差值来评估划分的质量。我将这个差值称为“绝对不准确度”，或简称为“不准确度”。相应地，最大和差值占平均和的百分比将被称为“相对不准确度”。

在“完美”划分的使用上存在一些不一致。在一些出版物中，“完美”表示对基本问题的陈述进行划分，而在另一些出版物中，则同时表示基本问题和优化问题。为了更确定，我们将使用“完美”（不带引号）一词来表示不准确度为零的划分。使用“完美”（带引号）一词，我们将表示优化问题的解决方案，即具有可能最小不准确度的划分。

贪心算法

已接受的贪心算法

众所周知的启发式方法，孩子们通常这样解决问题。它按降序遍历数字，将每个数字分配给和较小的子集。这种方法的运行时间为 O(n log n)。

替代方案

另一种启发式方法基于一个简单的想法：子集中已分配的数字越多，就越容易将剩余的小数字均匀分配，而不是大数字。因此，我们应该尝试从最大的数字开始填充每个子集，然后到最小的数字。

算法如下：首先，确定所有子集中数字的平均和。然后，按降序遍历数字，将每个数字分配给第一个子集，直到其和不超出某个阈值，该阈值定义为平均和加上 delta。如果添加下一个数字超出此阈值，则按降序遍历剩余数字，寻找第一个满足此限制条件的数字。当第一个子集不再需要这样的数字时，过程将为下一个子集重复。

在这种形式下，该算法存在一个严重的缺点：由于子集是根据接近和的唯一标准依次填充的，最后一个子集可能会留空，尤其是在 N 小于约 2k 的情况下，这对应于最大的可能不准确度。为了消除这个缺陷，在整个周期完成后，将未填充的数字按照贪心方法分配到子集中。但这只在数字数量不超过 k/2 的情况下进行，这是一个折衷的方法。否则，整个过程将再次与 delta 增加其初始值一起重复。然而，应该指出的是，这种调整仍然不能保证在 N≈k 的极端情况下填充所有子集。

分配初始 delta 值是一个重要环节。如果它太小，迭代次数（以及总时间）会不成比例地增加，而太大的值可能会导致糟糕（高不准确度）的划分。由于没有计算此值的数学标准，它根据测试选择为输入数字中最小值的 1/20，但不小于 1。这个启发式值，连同上述完成迭代的条件，在所有测试测量中都保证了单次传递中的令人满意的不准确度。

我将此算法称为“顺序填充贪心”（Sequential stuffing Greedy, SGA）。在一次迭代中，其复杂度也为 O(n log n)。总的来说，它在结果质量和速度方面与贪心算法竞争，但启发式常数和条件的存在，以及填充所有子集的保证不足，是该方法的薄弱环节。

简单

为了进行比较，提出一个更简单且重要的是更快的算法很有趣。它按降序遍历数字，将每个数字分配给下一个子集，并在到达最后一个子集或第一个子集时更改子集遍历的顺序。这种方法的运行时间为 O(n)。我称之为“无条件贪心”（Unconditional Greedy）算法，UGA（尽管“愚蠢贪心”也合适）。

这种方法确实是最快的，但它的划分质量通常是可以预见的最低（如图所示）。

动态搜索树

贪心算法的划分质量就像掷骰子。为了获得可靠的“完美”结果，我们必须使用蛮力。它通过动态搜索树来实现，树中的每个节点都分配了一个数字和一个子集。同样，所有数字都预先按降序排序。节点依次搜索剩余未分配的数字，并决定是将每个数字添加到分配的子集，还是添加到下一个子集。这两种选择都会创建一个新节点；这是“动态”一词的一种含义。拥有最后一个分配的子集的节点是最终的，即树的“底部”。它通过填充所有剩余未填充的数字来组装最后一个子集。

因此，树的深度受初始整数数量的限制，但每个节点的分支数量不是预先确定的。我们只能说它受当前未分配整数数量的限制。

这种解决方案考虑了加法的交换律，但不控制不同子集的组合频率。换句话说，组合 {n1, n2} 和 {n2, n1} 被保证排除，但两个子集 {n1, n2} 和 {n3, n4} 可能会在不同路径上以不同顺序重复出现。这种冗余是算法相对简单的结果。它在很大程度上被“预先完成”调整（见下文）所抵消。

此算法的复杂度为 O(n!)，非常耗时。为了使执行时间可接受，必须确保四个关键要求得到满足：速度最大化、自收紧精度选择、预先完成和时间限制。

加速

实现最高可能速度的明显条件是避免在迭代期间在堆上分配内存。解决方案是使用递归过程在堆栈上构建动态树。每个特定的树路径保留在堆栈上，直到到达最低点或直到它不再通过精度标准，然后让位给下一个。这是“动态”一词的另一种含义。在迭代期间，仅检查数字是否属于特定子集，所有子集的实际动态内存分配是在所有迭代结束时执行的。作为递归调用深度的输入整数的数量，对于实际任务是完全可以接受的。

自收紧精度选择

这一点与 SGA 方法相同：选择子集的输入数字，如果它们的和不超过某个上限。不同之处在于这次，我们要寻找所有未分配的数字。正如常规图所示，这种过滤会大大减少可能的树路径数量，并相应地减少执行时间。

再次出现关于此阈值大小或 delta 值的问题。这次，它最初被指定为 1（如果平均和超过最大数字，否则，该值将是它们的差值加上 2）。然而，尤其是在 N 数量较少的情况下，可能不存在满足此条件的组合。在这种情况下，将应用与之前相同的重复完整搜索方法，这次 delta 加倍。具有极高不准确度的路径稀疏且短，因此可以快速执行此类迭代。然而，delta 的连续加倍很快会导致总路径长度显著增加。为了解决这个矛盾，delta 是动态的：对于每个最终变体，检查其实际值，如果它小于当前指定值，则它成为当前值，并存储已实现的结果作为备用。否则，结果将被忽略。因此，在迭代过程中，delta 可以减小，呈现“自沉”状态。

此外，这次不仅使用了“delta 上限”（平均和的最大允许超限），还使用了“delta 下限”（平均和的最大允许低限）。这两个值是同步分配的，因此初始可接受的不准确度为 2。概括地说，是根据动态允许的不准确度进行选择。

预先完成

枚举的限制基于这样一个事实：随着 N 和 k 的增长，“完美”组合的数量也会增加。我们满足于找到的第一个。如果平均和没有小数部分，我们可以期望完美的划分，但不一定。找到的第一个此类结果会中断进一步的搜索。否则，最好的可能是一个不准确度为 1 的划分，并且同样，第一个此类结果会取消执行。

局限性

总的来说，我们无法确定当前找到的解决方案是否是“完美的”。不准确度为 1 的划分可能根本不存在，或者需要很长时间来搜索，所以我们别无选择，只能限制总执行时间（或迭代次数，本质上是一样的），并对在此期间找到的最佳结果感到满意。

具体来说，迭代次数限制为启发式值 1,000,000。该值在 168 个测试测量中只有 6 次给出的结果比某些贪心算法差，而算法执行时间对于分配 1000 个整数的最困难情况不超过约 0.5 秒（有关更多详细信息，请参阅测试）。稍微增加此限制会改善划分质量。例如，将限制加倍会将最差测试尝试次数减少到 3 次，但会使运行时间加倍。当限制低于 100 万时，DST 分区的质量会明显下降。因此，这个值可以被认为是最佳折衷。毕竟，迭代次数可以通过作为参数传递给类构造函数的乘数来管理。默认值为乘数值 1。

最终排列

在运气好的情况下，贪心算法可以提供一个不错的结果（见测试）。而且，由于完成它们所需的时间以微秒到数百微秒计，因此没有理由不去尝试运气。

组合解决方案包括所有四种算法，它们交替运行：UGA-GA-SGA-DST。每次轮换，最佳结果被保留为最终结果。具有极低不准确度（0 或 1）的结果会中断进一步计算。UGA 不会从链中排除，主要是为了在 N/k ≤ 2/3 的退化情况下获得极快的响应。

开发人员可以选择“快速”选项（不调用 DST 方法），此时划分质量不能保证，但执行时间不超过某个非常低的水平（粗略地说，在 2.5 GHz 的机器上约 0.1 毫秒）。

实现

下面的代码是用 C++ 编写的。如果没有文档注释和计时/测试插入，大约需要 500 行。唯一使用的数据容器是 std::vector。我没有使用 C++11（及更高版本）的特性，如 auto 说明符、lambda、仿函数和类型别名代替 typedef，主要是为了不限制编译器版本。在代码中，“bin”等同于“subset”。为了防止有关类与结构的问题：在对象仅存在于容器中的情况下，我使用结构。这只是一个约定问题。

该解决方案在一个名为 effPartition 的类中实现。输入数字用标识符标记，生成的划分是一组子集，而子集本身是数字 ID 的集合。

最初，给出了类型定义

typedef unsigned short  numb_id;    // type of number's ID
typedef unsigned int    numb_val;   // type of number's value
typedef unsigned char   ss_id;      // type of subset's ID
typedef unsigned long   ss_sum;     // type of subset's sum

输入数字由公共 IdNumber 结构表示

struct IdNumber
{
    friend class  Partition;
    friend struct Subset;
    friend class  IdNumbers;
    friend struct Result;
private:
    ss_id    BinIInd;   // bin's inverse index: count_of_bins-1 for the first and 0 for the last one
    numb_id  ID;        // number's ID
    numb_val Val;       // number's value

    ...
public:
    // Constructor by number's ID and value
    inline IdNumber(numb_id id, numb_val val) : BinIInd(0), Id(id), Val(val) {}

    // Returns number's value
    inline numb_val Value() const { return Val; }
 
    // Returns number's ID
    inline numb_val ID() const { return Id; }

    ...
};

一个公共的输入数字集合只是一个 IdNumber 的向量，带有 7 个额外的 private 方法

class IdNumbers : public std::vector<IdNumber>
{
    ...
}

最后一个公共结构 Subset 定义了一个划分数据

struct Subset
{
    friend class Partition;
    friend class effPartition;
private:
    ss_id    id;       // subset's ID
    ss_sum   sumVal;   // sum of number values
    numb_ids numbIDs;  // subset number IDs container

    // Takes into account (charges) number's value
    //  @it: number's iterator
    //  @binInd: bin's inverse index
    void TakeVal(const numb_it& it, ss_id binInd)
    {
        sumVal += it->Val;
        it->BinIInd = binInd;
    }

    // Eliminates (discharges) number's value
    //  @it: number's iterator
    void effPartition::Subset::ElimVal(const numb_it& it)
    {
        sumVal -= it->Val;
        it->BinIInd = 0;
    }

    ...

public:
    // Gets subset's ID
    inline ss_id ID() const { return id; }

    // Gets sum of subset numbers
    inline ss_sum Sum() const { return sumVal; }

    // Gets number IDs container
    inline const numb_ids& NumbIDs() const { return numbIDs; }
};

effPartition 类有三个构造函数，除了 IdNumbers、numb_val 的向量或 numb_val 的 C 数组。在后两种情况下，标识符会自动分配给整数，按照它们的顺序号。类似于 RAII idiom，划分在构造函数中执行。实际上，所有工作都由一个特殊的私有类 Partition 的对象完成，该对象在 effPartition 构造函数中创建，并在完成时自动释放所有临时变量。因此，在 effPartition 对象初始化后，它只保留最终结果作为子集向量，以及它的不准确度，并允许公共访问。

为简洁起见，我不会给出贪心算法的明显代码。这是仅包含 DST 方法的 Partition 类

class Partition
{
    static const UINT minCallLimit = 1000000;   // the minimum DSTree() invokes limit
    static const BYTE limFlag = 0x1;            // flag of DSTree() completion by limit
    static const BYTE perfFlag = 0x2;           // flag of DSTree() completion by 'perfect' result

    const float  avrSum;        // average sum value among subsets
    const bool   isAvrFract;    // true if avrSum is not an integer (has a nonzero fractional part)
    const ULLONG callLimit;     // current DSTree() invokes limit
    ULLONG  callCnt;            // counter of SGreedy() iterations or DSTree() invokes
    ss_sum  minSum;             // current minimum sum among subsets
    ss_sum  maxSum;             // current maximum sum among subsets
    ss_sum  sumDiff;            // current maximum difference between subset sums (inaccuracy)
    ss_sum  lastSumDiff;        // last best inaccuracy: used in DSTree()
    ss_sum  standbySumDiff;     // temporarily stored inaccuracy: used in DSTree()
    BYTE    complete;           // holder of DSTree() completion flags (signs)
    effPartition::Result&    finalResult;    // final result
    effPartition::Result     currResult;     // current result
    effPartition::IdNumbers& numbers;        // input numbers
    effPartition::IdNumbers  standbyNumbers; // numbers with the best unsaved maximum sum difference:
                                             // used in DSTree()
    ...
    
    // Performs 'Dynamic Search Tree' ('perfect') partition
    //  @currnit: number's iterator, from which the cycle continues
    //  @invInd: current inverse index of subset
    void DSTree(effPartition::numb_it currnit, ss_id invInd)
    {
        if(complete)    return;
        effPartition::numb_it nit;
        const ss_it sit = currResult.Bins.end() - 1 - invInd;    // curent bins iterator

        if(++callCnt == callLimit)    RaiseComplFlag(limFlag);

        if(invInd) {        // not last bin
            for(nit=currnit; nit!=numbers.end(); nit++)    // loop through the numbers starting from the current one
                if(!nit->BinIInd && nit->Val + sit->sumVal < maxSum) {   // accuracy selection
                    sit->TakeVal(nit, invInd);             // charge number's value
                    if(nit+1 != numbers.end())             // checkup just to avoid blank recursive invoke
                        DSTree(nit+1, invInd);             // try to fit next number to the same bin
                    if(sit->sumVal > minSum)               // bin is full?
                        DSTree(numbers.begin(), invInd-1); // try to fit unallocated numbers to the next bin
                    sit->ElimVal(nit);                     // discharge number's value
                }
        }
        else {                // last bin
            // accumulate sum for the last bin
            for(nit=numbers.begin(); nit!=numbers.end(); nit++)
                if(!nit->BinIInd)              // zero invIndex means that number belongs to the last bin
                    sit->sumVal += nit->Val;
            if(SetRange(currResult)) {         // is inaccuracy better than the previous one?
                standbyNumbers.Copy(numbers);  // keep current numbers as the standby one
                lastSumDiff = standbySumDiff = sumDiff;    // keep standbySumDiff for the 
                                                           // next standby selection
                if(IsResultPerfect())        RaiseComplFlag(perfFlag);
            }
            else if(currResult.SumDiff < standbySumDiff) {  // should we keep current result as standby?
                standbyNumbers.Copy(numbers);   // keep current numbers as the standby one
                standbySumDiff = currResult.SumDiff;
            }
            sit->sumVal = 0;   // clear last bin sum
        }
    }
    ...
};

完整的代码（包含测试和计时过程）已附上。

模板

编写一个通用版本似乎很诱人，因为只需定义数字的加法、减法和比较运算即可。然而，在这种情况下，“预先完成”要求将变得不可行。实际上，在大多数 N 小于 50 的情况下，正是这个要求提供了非常短的算法工作时间。此外，对 float 数字（以及其他非整数类型）的算术运算速度稍慢。因此，在划分 float 时，更有效的方法是将它们向下取整为整数并使用此解决方案：结果的质量不会改变，而执行时间保持在可接受的水平。

测试

测试在配备 Intel Core i5 2.5 GHz 的笔记本电脑上进行。代码以最大化速度选项编译。

将 10 到 1000 个范围在 50-500 之间的随机整数分配给 2 到 10（前 4 个系列）和 10 到 100（后 2 个系列）之间的不同数量的子集。每个划分都使用不同的 RNG 种子值执行 3 次，然后对结果进行平均，并绘制成不准确度和运行时间图。

每个算法都独立测试。因此，总共进行了 168 次测量。

结果摘要

根据下表中的统计数据，排除了 N=10 的系列，因为它们是相当低效的替代方案。

	UGA		GA		IGA		DSTree
平均相对不准确度，%	19.23		3.04		3.41		1.13
在贪心算法中最佳案例的数量	10	(7%)	61	(43%)	78	(55%)
在所有算法中最佳案例的数量	4	(3%)	10	(7%)	3	(2%)	133	(94%)
不准确度为 0 的案例数量	0		0		0		12
不准确度为 1 的案例数量	0		2		2		29

最佳案例总数占总案例数的比例超过 100%，因为有时两个算法显示相同的最佳结果，这种情况被计算在两者之内。

不考虑 N=10 的系列，DST 算法在 16 个可能的案例中有 12 个完美地分配了整数。然而，必须反复强调，所有整数的平均和（无小数部分）并不一定意味着存在不准确度为零的划分。

最后，在这些测试中，组合解决方案的平均相对不准确度为 1.004%。我必须强调，这个统计数据是纯粹相对的。

结果详情

图中显示了平均划分绝对不准确度的分布

最后两个案例是“硬”案例，更属于“集装箱和货车”的任务。

这是平均执行时间

所有算法都能正确处理单个子集或初始数字相等的退化情况。

运行时间：托管代码 vs 非托管代码

将 30 个范围在 50-500 之间的随机数分配给 4、6 和 8 个子集，分别通过 C++、C# 和 Java 实现。通过启动每个程序，effPartition 对象被创建 3 次，每次使用不同数量的子集。3 次启动之间的结果平均值显示在表中

lang	子 子集	UGA			GA			SGA			DSTree
		时间，毫秒	次	不准确度	时间，毫秒	次	不准确度	时间，毫秒	次	不准确度	时间，毫秒	次	不准确度
C++	4	0.0003	1	72	0.0004	1	38	0.0006	1	39	0.006	1	0
C#		1.4	4491		1.9	4762		1.2	1873		2.8	491
Java		0.68	2184		0.22	551		0.44	686		0.81	143
C++	6	0.0003	1	205	0.0005	1	31	0.0008	1	45	0.01	1	0
C#		0.005	18		0.007	16		0.013	17		0.047	5
Java		0.036	125		0.074	158		0.10	136		0.30	30
C++	8	0.0003	1	142	0.0006	1	55	0.0009	1	38	44.9	1	2
C#		0.002	6		0.005	10		0.011	12		104.0	2
Java		0.03	101		0.08	139		0.12	126		100.8	2

粗体字体显示了每个子集执行时间的相对关系，其中 C++ 被设为 1。

第一个显而易见的问题是，在 .NET 和 Java 的情况下，第一次创建 effPartition 对象和随后的对象之间，执行时间差异如此之大。很明显，CLR 在调用秒表测量的函数之前，会在堆上分配一个数据段。显然，垃圾回收器从未被调用。特殊之处在于，所有分区方法都积极引用初始数字的集合。并且所有方法（除了 UGreedy）也积极引用一个临时的 Result 对象（但实际上只在最后一次动态分配内存）。快速实验表明，“迟疑”方法首次调用的最可能原因是对象在 CPU 缓存（C#）或缓存映射（Java）中驻留。在这种情况下，C# 和 Java 中首次和后续分区方法调用运行时间差异可以用这些缓存策略效率的差异来解释。

第二个问题是，在 8 个子集的情况下，托管 DST 方法实现比非托管实现慢了近一倍，而在 4 个和 6 个子集的情况下则不然。为了获得这个结果，该方法被调用了二百万次，delta 分别加倍为 2 和 4，而对于 4 个和 6 个子集，分别只需要 247 次和 487 次调用。在如此大量的递归调用下，初始化托管数据内存所需的时间变得无关紧要。最有可能的是，托管和非托管代码之间运行时间上的微小差异可以用两个因素来解释：与索引（C# 和 Java）相比，使用迭代器（C++）可以更快地遍历集合，以及直接访问集合元素（C++）与通过引用访问（C# 和 Java）。

对托管内存特性影响的详细分析超出了主题范围。主要的是，即使是 C# 和 Java 解决方案的执行时间，在最困难的 N=1000 和 k≤100 的测试中也不超过 0.7 秒。

结论

回到我最初的应用任务，它涉及到基因组文库的计算机模拟处理。每个染色体都有一个独特的大小，直接与其处理时间成正比，并且独立处理。下表显示了鼠标基因组并行处理总时间为 10 分钟时，每个 4、6 和 8 个核心的最大空闲时间。使用“完美”解决方案选项（DST 方法），在前两种情况下，所有处理器几乎同时关闭。在最后一种情况下，空闲时间约为半分钟；但必须考虑到，将 22 个差异很大的大整数划分为 8 个子集的更好划分可能根本不存在。

子集	4		6		8
	'快速'	'完美'	'快速'	'完美'	'快速'	'完美'
相对不准确度，%	3.03	0.016	11.52	0.23	20.52	5.66
运行时间，毫秒	0.007	22	0.009	35	0.01	110
总计 10 分钟；差异时间，秒	18.2	0.1	69.1	1.4	123.1	33.9

不能断言这个组合解决方案总是给出“完美”的结果。然而，在实际意义上，它是。至少，它是质量和速度之间最佳的折衷。

致谢

我感谢Oleg Tolmachev（英国格林威治大学）为改进本文的语言所提供的宝贵帮助。