Suppose we are given a finite list of real numbers \(a_1,\ldots,a_n\). We want to divide the index set \(\{1,\ldots,n\}\) into \(r\) non-empty pieces \[ \pi=\{p_1,\ldots,p_r\} \] so that the sums over the pieces are as close to each other as possible. For a block \(p\), put \[ S(p)=\sum_{i\in p}a_i . \] The discrepancy of the partition is \[ \nu(\pi)=\max_{p\in\pi}S(p)-\min_{p\in\pi}S(p). \] Thus, for fixed \(r\), the problem is to compute \[ \nu_{n,r}(a_1,\ldots,a_n)= \min_{\pi\in\Pi_{n,r}} \left( \max_{p\in\pi}\sum_{i\in p}a_i - \min_{p\in\pi}\sum_{i\in p}a_i \right), \] where \(\Pi_{n,r}\) is the set of partitions of \(\{1,\ldots,n\}\) into exactly \(r\) non-empty blocks.
I will write an exact algorithm. It is not meant to be a fast heuristic for large inputs. The problem is already hard for \(r=2\), so one should not expect a polynomial-time miracle.
function intervals(m, t):
if |m| < t:
return empty set
if t = 1:
return { [w(m), w(m)] with block m }
first = least element of m
result = empty set
for each subset B of m containing first:
R = m \ B
if |R| >= t - 1:
for each interval [L,U] in intervals(R, t - 1):
s = w(B)
add [min(s,L), max(s,U)] with block B to result
return Pareto-minimal intervals in result
answer:
choose [L,U] in intervals({1,...,n}, r) minimizing U - LThe algorithm is a dynamic programming procedure over subsets. A
subset of the index set is represented by a bit mask \(m\), and \(w(m)\) denotes the sum of the entries whose
indices occur in \(m\). For every pair
\((m,t)\), the function
intervals(m,t) stores the Pareto-minimal intervals \([L,U]\) that can occur when the subset
\(m\) is partitioned into exactly \(t\) non-empty blocks, with all block sums
lying between \(L\) and \(U\). The case \(t=1\) is immediate: the only possible block
is the whole subset, hence the only interval is \([w(m),w(m)]\). For \(t>1\), the algorithm chooses the block
containing the smallest element of \(m\), recursively partitions the complement
into \(t-1\) blocks, and updates the
resulting interval by adjoining the sum of the chosen block. Requiring
the chosen block to contain the smallest available element removes the
permutation symmetry among the blocks, so each unordered partition is
generated once. After each step, dominated intervals are discarded: if
\([L,U]\subseteq
[L',U']\), then \([L',U']\) can never lead to a
better range and is removed. The final answer is the interval of least
width among intervals(fullMask,r).
Let us start with base classes:
import scala.collection.mutable
import java.lang.Long.{bitCount, numberOfTrailingZeros}
case class Node(lo: Double, hi: Double, blocks: List[Long]) {
def width: Double = hi - lo
}
case class Answer(value: Double, parts: Vector[Vector[Int]], sums: Vector[Double])import scala.collection.mutable
import java.lang.Long.{bitCount, numberOfTrailingZeros}
class Node
class AnswerNow, the main solver:
def solve(x: IndexedSeq[Double], r: Int, eps: Double = 1e-10): Answer = {
require(1 <= r && r <= x.length)
require(x.length <= 62)
val n = x.length
val sums = mutable.HashMap[Long, Double](0L -> 0.0)
val memo = mutable.HashMap.empty[(Long, Int), Vector[Node]]
def sum(m: Long): Double =
sums.getOrElseUpdate(
m,
sum(m & (m - 1L)) + x(numberOfTrailingZeros(m))
)
def subsets(m: Long): Iterator[Long] =
Iterator.iterate(m)(s => (s - 1L) & m).takeWhile(_ != 0L) ++ Iterator.single(0L)
def dominates(a: Node, b: Node): Boolean =
a.lo >= b.lo - eps &&
a.hi <= b.hi + eps &&
(a.lo > b.lo + eps || a.hi < b.hi - eps)
def pareto(v: Vector[Node]): Vector[Node] =
v.filter(b => !v.exists(a => dominates(a, b)))
def intervals(m: Long, t: Int): Vector[Node] =
memo.getOrElseUpdate(
(m, t),
{
if (bitCount(m) < t) Vector.empty
else if (t == 1) {
val s = sum(m)
Vector(Node(s, s, m :: Nil))
}
else {
val first = m & -m
pareto(
subsets(m ^ first).flatMap { sub =>
val block = sub | first
val rem = m ^ block
if (bitCount(rem) < t - 1) Iterator.empty
else {
val s = sum(block)
intervals(rem, t - 1).iterator.map { z =>
Node(math.min(s, z.lo), math.max(s, z.hi), block :: z.blocks)
}
}
}.toVector
)
}
}
)
val best = intervals((1L << n) - 1L, r).minBy(_.width)
Answer(
value = best.width,
parts = best.blocks
.map(m => (0 until n).filter(i => ((m >>> i) & 1L) == 1L).toVector)
.toVector,
sums = best.blocks.map(sum).toVector
)
}def solve(x: IndexedSeq[Double], r: Int, eps: Double): AnswerThe only slightly non-obvious line is the Pareto reduction. Suppose one construction gives the interval \([L,U]\), and another gives \([L',U']\). If \[ L\geq L' \quad\text{and}\quad U\leq U', \] then \([L,U]\subseteq [L',U']\). The larger interval cannot be better for minimizing \(U-L\). So the larger interval may be discarded.
Let us test the function on a small example.
val a = solve(Vector(1.0, 0.45, 0.40, 0.40, 0.10, 0.10), 4)
a.valueval a: Answer = Answer(0.55,Vector(Vector(0), Vector(1), Vector(2, 5), Vector(3, 4)),Vector(1.0, 0.45, 0.5, 0.5))
val res0: Double = 0.55The value is the minimum possible discrepancy among all partitions into \(4\) blocks. The returned partition is encoded by zero-based indices.
a.parts.zip(a.sums).zipWithIndex.foreach {
case ((p, s), j) => println(s"p$j = $p ; sum = $s")
}p0 = Vector(0) ; sum = 1.0
p1 = Vector(1) ; sum = 0.45
p2 = Vector(2, 5) ; sum = 0.5
p3 = Vector(3, 4) ; sum = 0.5Here is a reproducible random example.
val rng = new scala.util.Random(2026L)
val x = Vector.fill(12)(rng.nextDouble())
val b = solve(x, 3)val rng: scala.util.Random = scala.util.Random@6f430ea8
val x: scala.collection.immutable.Vector[Double] = Vector(0.6186327814734252, 0.849209341093155, 0.7903839677394705, 0.5296053232288097, 0.8416668624371814,
0.14829601053373576, 0.8437822755735771, 0.5847572137985457, 0.1667271317594733, 0.7157842860912638, 0.5392902493576879, 0.27949326840295363)
val b: Answer = Answer(0.007453088971694388,Vector(Vector(0, 4, 6), Vector(1, 5, 7, 9), Vector(2, 3, 8, 10, 11)),Vector(2.304081919484184, 2.2980468515167005, 2.305499940488395))b.valueval res2: Double = 0.007453088971694388The three block sums are:
b.sumsval res3: Vector[Double] = Vector(2.304081919484184, 2.2980468515167005, 2.305499940488395)Their range is:
b.sums.max - b.sums.minval res4: Double = 0.007453088971694388which agrees with the value returned by solve.
The complexity is exponential. The memoized state is a pair \((m,t)\), where \(m\) is a subset mask and \(t\) is the number of blocks still required. There are at most \(r2^n\) such states, and each state may enumerate many submasks. Thus the algorithm is exact, but it is suitable only for moderately sized inputs.
The important point is that this function is not the greedy equal-pile procedure. The greedy procedure gives a fast balancing heuristic. The present function computes the actual optimum
\[ \min_{\pi\in\Pi_{n,r}} \left( \max_{p\in\pi}\sum_{i\in p}a_i - \min_{p\in\pi}\sum_{i\in p}a_i \right). \]