How to optimize a function involving max

Question

I am having problems minimizing a simple if slightly idiosyncratic function. I have scipy.optimize.minimize but I can't get consistent results. Here is the full code:

from math import log, exp, sqrt
from bisect import bisect_left
from scipy.optimize import minimize
from scipy.optimize import Bounds
import numpy as np


def new_inflection(x0, x1):
    return log((exp(x0)+exp(x1) + sqrt(exp(2*x0)+6*exp(x0+x1)+exp(2*x1)))/2)


def make_pairs(points):
    new_points = []
    for i in range(len(points)):
        for j in range(i+1, len(points)):
            new_point = new_inflection(points[i], points[j])
            new_points.append(new_point)
    return new_points


def find_closest_number(numbers, query):
    index = bisect_left(numbers, query)
    if index == 0:
        return numbers[0]
    if index == len(numbers):
        return numbers[-1]
    before = numbers[index - 1]
    after = numbers[index]
    if after - query < query - before:
        return after
    else:
        return before

    
def max_distance(target_points):
    pair_points = make_pairs(target_points)
    target_points = sorted(target_points)
    dists = []
    return max(abs(point - find_closest_number(target_points, point)) for point in pair_points)



num_points = 20
points = np.random.rand(num_points)*10
print("Starting score:", max_distance(points))
bounds = Bounds([0]*num_points, [num_points] * num_points)
res = minimize(max_distance, points, bounds = bounds, options={'maxiter': 100}, method="SLSQP")
print([round(x,2) for x in res.x])
print(res)

Every time I run it I get quite different results. This is despite the output saying Optimization terminated successfully. An example output:

 message: Optimization terminated successfully
 success: True
  status: 0
     fun: 0.4277378933292031
       x: [ 5.710e+00  1.963e+00 ...  1.479e+00  6.775e+00]
     nit: 15
     jac: [ 0.000e+00  0.000e+00 ...  0.000e+00  0.000e+00]
    nfev: 364
    njev: 15

Sometimes I get a result as low as 0.40 and other times as high as 0.51.

Is there any way to optimize this function properly in Python?

Answer 1

The problem here is that you are exploring immense non convex search space by trying to optimize 20 variables at the same time, so usual optimization methods such as gradient descent related etc. will potentially be trapped in a local minima. In your case, as you are starting from a random coordinates each time, you end in a different local minima each time.

If the problem has no analytical solution, there is no optimizer (at least that I know) that can solve this kind of problem and guarantee you that you are in the global minimum, but that doesn't mean that we can get a very good approximations.

Approach One

Try different types of algorithms more suited for non convex high dimension space optimization, within scipy you have a global optimizer function:

from scipy import optimize
optimize.shgo(eggholder, bounds)

In my case that was very slow, but maybe it can help you.

Update: The global optimizer basinhopping seems to give good results faster:

from scipy.optimize import basinhopping
res = basinhopping(max_distance, points, minimizer_kwargs={"method": "SLSQP", "bounds": bounds}, niter=100)
print([round(x,2) for x in res.x])
print(res)

With this optimizer you can also reach a fitness of 3.9 consistently.

Approach Two

I would give a try to genetic algorithms using pygad, here is my trial and it reaches 3.9 fitness and consistently bellow 4.1 (though I change the sign for the optimizer).

import pygad

def fitness_func(solution, solution_idx):
    return -max_distance(solution)

fitness_function = fitness_func

num_generations = 500
num_parents_mating = 5

sol_per_pop = 100
num_genes = len(points)

init_range_low = 0
init_range_high = 20

parent_selection_type = "sss"
keep_parents = 5

crossover_type = "single_point"

mutation_type = "random"
mutation_percent_genes = 10

ga_instance = pygad.GA(num_generations=num_generations,
                       num_parents_mating=num_parents_mating,
                       fitness_func=fitness_function,
                       sol_per_pop=sol_per_pop,
                       num_genes=num_genes,
                       gene_space = {'low': 0, 'high': 20},
                       init_range_low=init_range_low,
                       init_range_high=init_range_high,
                       parent_selection_type=parent_selection_type,
                       keep_parents=keep_parents,
                       crossover_type=crossover_type,
                       mutation_type=mutation_type,
                       mutation_percent_genes=mutation_percent_genes)

ga_instance.run()

To show the solution:

solution, solution_fitness, solution_idx = ga_instance.best_solution()
print("Parameters of the best solution : {solution}".format(solution=solution))
print("Fitness value of the best solution = {solution_fitness}".format(solution_fitness=solution_fitness))
# plots the optimization process
ga_instance.plot_fitness(title="PyGAD with Adaptive Mutation")

Approach Three

Probably the worst idea but also maybe is everything you need, Just iterate many times the algorithm you already have, with different initial conditions saving the best solution.

Update As the OP commented, somewhat similar to this is what the algorithm basinhopping is doing (with a bit more complexity) but it archives very good results, see approach one.

Conclusion

There is very little you can do for having consistent results in the optimization algorithm in a problem like this one (if not fixing the seed, which I strongly discourage), but at least you can choose an algorithm that is more suited for the task maximizing the search space, so you get closer and closer to the global minimum.

On the other hand you should notice that there may be multiple/infinite solutions with the same fitness, but that can look very different, specially if there are hidden symmetries in the problem, this problem for instance seems invariant under permutations, so it would make sense to sort the solution list. Also in this case little changes in one element of the 20 numbers doesn't necessarily change the fitness.