Given an integer $n$, and 2 real sequences $\{a_1, \dots, a_n\}$ and $\{b_1, \dots, b_n\}$, with $a_i$, $b_i$ > 0, for all $i$. For a fixed $m < n$ let $\{P_1, \dots, P_m\}$ be a partition of the set $\{1, \dots, n\}$ as in $P_1 \cup \dots \cup P_m$ = $\{1, \dots, n\}$, with the $P_i$'s pairwise disjoint. I wish to find a partition of size $m$ that solves
$$\max_{P=\{P_1, ..., P_m\}}\sum_{j=1}^{m}\frac{(\sum_{i \in P_j}a_i)^2}{\sum_{i \in P_j}b_i}$$
I am really looking for an algorithm which solves the problem in polynomial time, a brute-force solution is not feasible, as it would involve the a Bell number of order $(n, k)$, with $n$ over 1e6 for realistic cases.
I would be happy to prove that the partition is monotonic in increasing values of $a/b$, in the sense that a partition expressed in the indices of the two sequences $a, b$, sorted by increasing values of $a/b$ will contain monotonic, increasing sets of integers in $1, ..., n$. I believe this is the case - can someone provide a proof?
If so, the brute-force search could be improved to an order $n \choose m-1$ algorithm, still long, but a significant savings.
The script below solves the problem by brute-force. For example, a sample run with $n = 12$, $m = 3$, gives an optimal partition of (expressed in indices of the sorted sequence $a/b$):
[[0, 1, 2, 3, 4, 5, 6, 7], [8, 9], [10, 11]]
which is monotonic, as claimed.
import numpy as np
import multiprocessing
import concurrent.futures
from functools import partial
from itertools import chain, islice
# n
NUM_POINTS = 12
# m
PARTITION_SIZE = 4
rng = np.random.RandomState(55)
def knuth_partition(ns, m):
def visit(n, a):
ps = [[] for i in range(m)]
for j in range(n):
ps[a[j + 1]].append(ns[j])
return ps
def f(mu, nu, sigma, n, a):
if mu == 2:
yield visit(n, a)
else:
for v in f(mu - 1, nu - 1, (mu + sigma) % 2, n, a):
yield v
if nu == mu + 1:
a[mu] = mu - 1
yield visit(n, a)
while a[nu] > 0:
a[nu] = a[nu] - 1
yield visit(n, a)
elif nu > mu + 1:
if (mu + sigma) % 2 == 1:
a[nu - 1] = mu - 1
else:
a[mu] = mu - 1
if (a[nu] + sigma) % 2 == 1:
for v in b(mu, nu - 1, 0, n, a):
yield v
else:
for v in f(mu, nu - 1, 0, n, a):
yield v
while a[nu] > 0:
a[nu] = a[nu] - 1
if (a[nu] + sigma) % 2 == 1:
for v in b(mu, nu - 1, 0, n, a):
yield v
else:
for v in f(mu, nu - 1, 0, n, a):
yield v
def b(mu, nu, sigma, n, a):
if nu == mu + 1:
while a[nu] < mu - 1:
yield visit(n, a)
a[nu] = a[nu] + 1
yield visit(n, a)
a[mu] = 0
elif nu > mu + 1:
if (a[nu] + sigma) % 2 == 1:
for v in f(mu, nu - 1, 0, n, a):
yield v
else:
for v in b(mu, nu - 1, 0, n, a):
yield v
while a[nu] < mu - 1:
a[nu] = a[nu] + 1
if (a[nu] + sigma) % 2 == 1:
for v in f(mu, nu - 1, 0, n, a):
yield v
else:
for v in b(mu, nu - 1, 0, n, a):
yield v
if (mu + sigma) % 2 == 1:
a[nu - 1] = 0
else:
a[mu] = 0
if mu == 2:
yield visit(n, a)
else:
for v in b(mu - 1, nu - 1, (mu + sigma) % 2, n, a):
yield v
n = len(ns)
a = [0] * (n + 1)
for j in range(1, m + 1):
a[n - m + j] = j - 1
return f(m, n, 0, n, a)
def Bell_n_k(n, k):
''' Number of partitions of {1,...,n} into
k subsets, a restricted Bell number
'''
if (n == 0 or k == 0 or k > n):
return 0
if (k == 1 or k == n):
return 1
return (k * Bell_n_k(n - 1, k) +
Bell_n_k(n - 1, k - 1))
def slice_partitions(partitions):
# Have to consume it; can't split work on generator
partitions = list(partitions)
num_partitions = len(partitions)
bin_ends = list(range(0,num_partitions,int(num_partitions/NUM_WORKERS)))
bin_ends = bin_ends + [num_partitions] if num_partitions/NUM_WORKERS else bin_ends
islice_on = list(zip(bin_ends[:-1], bin_ends[1:]))
rng.shuffle(partitions)
slices = [list(islice(partitions, *ind)) for ind in islice_on]
return slices
def reduce(return_values, fn):
return fn(return_values, key=lambda x: x[0])
class SimpleTask(object):
def __init__(self, a, b):
self.a = a
self.b = b
def __call__(self):
time.sleep(1)
return '{self.a} * {self.b} = {product}'.format(self=self, product=self.a * self.b)
def __str__(self):
return '{self.a} * {self.b}'.format(self=self)
class Task(object):
def __init__(self, a, b, partition):
self.partition = partition
self.task = partial(Task._task, a, b)
def __call__(self):
return self.task(self.partition)
@staticmethod
def _task(a, b, partitions, report_each=1000):
max_sum = float('-inf')
arg_max = -1
for ind,part in enumerate(partitions):
val = 0
part_val = [0] * len(part)
part_vertex = [0] * len(part)
for part_ind, p in enumerate(part):
part_sum = sum(a[p])**2/sum(b[p])
part_vertex[part_ind] = part_sum
part_val[part_ind] = part_sum
val += part_sum
if val > max_sum:
max_sum = val
arg_max = part
max_part_vertex = part_vertex
# if not ind%report_each:
# print('Percent complete: {:.{prec}f}'.
# format(100*len(slices)*ind/num_partitions, prec=2))
return (max_sum, arg_max, max_part_vertex)
class Worker(multiprocessing.Process):
def __init__(self, task_queue, result_queue):
multiprocessing.Process.__init__(self)
self.task_queue = task_queue
self.result_queue = result_queue
def run(self):
proc_name = self.name
while True:
task = self.task_queue.get()
if task is None:
# print('Exiting: {}'.format(proc_name))
self.task_queue.task_done()
break
result = task()
self.task_queue.task_done()
self.result_queue.put(result)
NUM_WORKERS = multiprocessing.cpu_count() - 1
INT_LIST= range(0, NUM_POINTS)
num_partitions = Bell_n_k(NUM_POINTS, PARTITION_SIZE)
partitions = knuth_partition(INT_LIST, PARTITION_SIZE)
slices = slice_partitions(partitions)
while True:
a0 = rng.uniform(low=-0.0, high=100.0, size=NUM_POINTS)
b0 = rng.uniform(low=-0.0, high=100.0, size=NUM_POINTS)
# sort by increasing a/b, to check claim
ind = np.argsort(a0/b0)
(a,b) = (seq[ind] for seq in (a0,b0))
tasks = multiprocessing.JoinableQueue()
results = multiprocessing.Queue()
workers = [Worker(tasks, results) for i in range(NUM_WORKERS)]
num_slices = len(slices) # should be the same as NUM_WORKERS
for worker in workers:
worker.start()
for i,slice in enumerate(slices):
tasks.put(Task(a, b, slice))
for i in range(NUM_WORKERS):
tasks.put(None)
tasks.join()
allResults = list()
slices_left = num_slices
while not results.empty():
result = results.get()
allResults.append(result)
# print('result: {!r}'.format(result))
slices_left -= 1
r_max = reduce(allResults, max)
c = a/b
part = r_max[1]
endpoints = [(a[-1], b[0]) for a,b in zip(part[:-1], part[1:])]
d = [(c[r]-c[l]) for l,r in endpoints]
r = [(c[r]-c[l])/c[l] for l,r in endpoints]
all_diffs = np.concatenate([[np.nan], np.diff(c)])
all_rets = np.concatenate([np.diff(c), [np.nan]]) / c
max_diffs = sorted(all_diffs)[-(PARTITION_SIZE-1):]
max_rets = sorted(all_rets)[-(PARTITION_SIZE-1):]
print('TRIAL: {} : max: {:4.6f} pttion: {!r}'.format(i, *r_max[:-1]))
# print('TRIAL: {} : max: {:4.6f} {!r} {!r}'.format(i, *r_max[:-1],
# [float(x)
# for x in ['{0:0.2f}'.format(i)
# for i in r_max[2]]], prec=2))
try:
assert all(np.diff(list(chain.from_iterable(r_max[1]))) == 1)
except AssertionError as e:
import pdb
pdb.set_trace()