fcmaes - a Python 3 gradient-free optimization library

Find the Optimal Subset(s)

This tutorial

shows how to find a subset of a given set that optimizes a chosen property.
shows how to map continuous decision variables to subset selection based on numpy arrays.
provides examples for both single- and multi-objective optimization.

Motivation

This tutorial was inspired by this reddit post. The author shared a good open source solution based on CVXPY and GLPK_MI, a mixed-integer convex optimizer. Here we look at two questions:

Can this problem be solved using continuous optimization?
Are there any advantages over using CVXPY?

The first advantage is obvious. Continuous optimization is not limited to convex solution landscapes, so the approach is more general. But even when the problem is convex, as in the example below, there are still good reasons to consider it:

GLPK_MI delivers only one solution for a given problem instance, even if valid alternatives exist.
Executing GLPK_MI again always results in the same solution.
GLPK_MI struggles with large problem instances.
GLPK_MI cannot compute the pareto-front if multiple competing objectives are defined.
GLPK_MI does not use modern many-core CPU architectures.

In contrast, continuous optimization:

Computes multiple alternative solutions utilizing parallel threads.
Can handle huge problem instances.
Supports multiple objectives and computes a set of non-dominated solutions.
Size and complexity of the code and performance are similar to GLPK_MI.

Implementation

The example code is at subset.py. First we define the value of a specific subset selection. The goal is to identify the subset of transactions that matches a set of payments, so we compare the two sums.

class transaction_value():

    def __init__(self, transactions, payments):
        self.transactions = transactions
        self.sum_payments = sum(payments)

    def __call__(self, selection):
        return abs(sum(self.transactions[selection]) - self.sum_payments)

The function returns the difference between the sum of the selected transactions and the sum of the payments. To use continuous optimization, we represent the problem as a Python fitness class. This class can be reused for other selection-to-value mappings.

class fitness():

    def __init__(self, selection_value, dim):
        self.selection_value = selection_value
        self.bounds = Bounds([0]*dim, [1.99999999]*dim)

    # all decision variables are in the [0,2[ interval and mapped to a boolean array.
    def __call__(self, x):
        return self.selection_value(x.astype(int).astype(bool))

The trick is simple. We use a vector of continuous decision variables in the [0,2[ interval and convert it with x.astype(int).astype(bool) into a boolean vector. That boolean vector is then used for the selection via transactions[selection].

We run the optimization with parallel retry. We do not use retry.minimize here because we want more than the single best solution. We want a list of good alternatives. retry.retry uses all available CPU cores, mp.cpu_count(), but the number of parallel workers is configurable.

    def optimize(fitness, opt, num_retries = 32):
        store = retry.Store(wrapper(fitness), fitness.bounds)
        retry.retry(store, opt.minimize, num_retries)
        xs = store.get_xs()
        ys = store.get_ys()
        ...

Finally, we choose a concrete problem instance and an optimizer. For this task, the BiteOpt algorithm combined with parallel retry performs best.

    transactions= rng.integers(100, 2500, 1000) / 100
    payments = rng.integers(10, 50, 100)
    selection_value = transaction_value(transactions, payments)
    fit = fitness(selection_value, len(transactions))
    opt = Bite_cpp(50000, popsize=500)
    optimize(fit, opt, num_retries=32)

Unlike De_cpp, Bite_cpp does not need a configuration parameter that marks variables with discrete values.

Results

Executing subset.py produces a result in less than one second on an AMD 5950x 16 core CPU for a superset of 1000 transactions:

...
0.86 353695 411273.0 0.0 [0.404784861806257, 0.3056760325932503, ...]
1 ) Optimal Objective value:  0.0
[0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 1 0 0 0
 0 0 0 1 1 0 0 0 0 0 0 0 0 1 0 0 0 1 1 0 0 1 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0
 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 1 1 0 0 0 0 0 1 0 0 0 0 0 0 1 1 0 0 0 0 1 0
 0 0 0 0 0 0 0 0 0 0 1 1 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 1 0 0
 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0
 1 0 0 0 1 1 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 1 0 0 0 0
 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1
 1 0 0 0 0 0 1 1 0 0 0 1 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 1 1 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 1 0 0 0 0
 0 1 0 0 1 1 0 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 1 0 0 0 1 1 0 0 0 0 0 0 0 0
 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
 0 0 0 0 1 0 1 1 0 0 0 0 0 1 0 0 0 1 0 0 1 0 1 1 0 0 0 0 0 0 1 1 0 0 0 0 0
 1 0 0 0 0 0 0 0 0 1 0 1 1 0 1 0 0 0 1 1 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1 0
 0 0 1 0 0 0 0 0 1 0 0 1 0 1 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0
 0 1 0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 1 0 1 0 0 0 0 1 0 1 0 1 0 0 0 0 0 1 1 0
 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 1
 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 1 0 1 0 0 1 1 0 0 0 0
 1 0 1 1 0 0 1 0 0 0 0 0 0 1 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 1 0 0 1 0
 0 1 1 0 1 0 1 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 1 0 0 0 1 0 1 1 0 0 1 0 0 1 0 0 0 0 1 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0
 1 0 1 0 1 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 1 1 0 0 0 0 0 1
 0 1 0 0 0 0 1 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 1 1 0 1 0 1 1 1 0 0 0 0 0 1
 0 0 1 0 1 1 0 0 0 0 0 0 1 1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0
 0 0 0 1 0 1 0 0 0 0 0 1 0 1 1 1 1 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 0 0
 0 0 0 1 1 1 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0 1 1 0 0 1 0 0 0 0 0 0 1 0 0 0 0
 1 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0
 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 1 0 1 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0 0
 0]
2 ) Optimal Objective value:  0.0
[0 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0

Excercise

Compare the performance of different optimization algorithms from fcmaes.optimize, such as de_cma, Cma_cpp, De_cpp, Da_cpp, Csma_cpp, Bite_cpp and Crfmnes_cpp.

Multi-Objective Optimization

The example code is at subset_mo.py. Now add a second objective. We want to maximize the minimal selected transaction, and we want to compute a set of non-dominated solutions, a pareto front. transaction_value only needs a small modification to return the second objective:

class transaction_value():

    def __init__(self, transactions, payments):
        self.transactions = transactions
        self.sum_payments = sum(payments)

    def __call__(self, selection):
        trs = self.transactions[selection]
        return abs(sum(trs) - self.sum_payments), -min(trs)

As optimizer we choose modecpp with NSGA-II population update by setting nsga_update = True.

    def optimize(fitness, num_retries = 32):
        nobj = 2
        ncon = 0
        xs, ys = modecpp.retry(mode.wrapper(fitness, nobj), nobj, ncon,
                               fit.bounds, num_retries=num_retries, popsize = 500,
                               max_evaluations = 100000, nsga_update = True,
                               logger = logger(), workers=32)
    ...

Results

Executing subset_mo.py takes more time than the single-objective case. It needs almost 30 seconds on an AMD 5950x 16 core CPU:

...
retries = 32: time = 27.0 i = 32000
[(0.0, -2.46), (4.547473508864641e-13, -2.65), (1.3642420526593924e-12, -2.97), (2.720000000001164, -2.99), (5.199999999998909, -3.15), (136.39999999999736, -3.17)]
1 ) Optimal Objective values:  [ 0.   -2.46]
[0 0 0 0 1 1 0 0 1 0 0 1 1 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 1 1 1 0 0 0 0 1 1
 0 0 0 0 0 1 1 0 0 0 0 1 0 0 0 0 0 0 1 0 0 1 0 1 1 0 0 0 0 0 1 1 0 1 0 1 0
...

Conclusion

Continuous optimization, especially the BiteOpt algorithm used with parallel retry, is well suited for optimization problems that involve selecting a subset.
Many alternative solutions can be computed in a single run.
Performance on a modern many-core CPU is comparable to CVXPY/GLPK_MI.
Adding another objective and computing the pareto-front is easy.
This approach also handles noisy fitness or a rugged solution landscape well.