a random-sampling heuristic for optimization that usually works better than i'd like it to

rsh.py

#!/usr/bin/env python3

import numpy as np

# get this here: https://github.com/imh/hipsterplot/blob/master/hipsterplot.py
from hipsterplot import plot

def heuristic(costs, deltas, center_cost):
    d = np.linalg.norm(deltas, axis=1)
    c0 = costs[:, 0] - center_cost
    c1 = costs[:, 1] - center_cost
    l0 = np.abs(3 * c1 + c0)
    l1 = np.abs(c1 + 3 * c0)
    peak = np.maximum(l0, l1) / (2 * d)
    return costs / peak[:, None]

def populate(f, center, popsize, sigma):
    costs, deltas = [], []
    for p in range(popsize):
        delta = np.random.normal(scale=sigma, size=center.shape)
        cost_pos = f(center + delta)
        cost_neg = f(center - delta)
        costs.append((cost_pos, cost_neg))
        deltas.append(delta)
    return np.array(costs, float), np.array(deltas, float)

def minimize(objective, init, iterations,
             sigma=0.1, popsize=None, step_size=1.0,
             true_objective=None):
    if popsize is None:
        # it's still better to provide one yourself.
        popsize = int(np.sqrt(len(init)))

    center = np.array(init, float, copy=True)
    center_cost = objective(center)
    history = []

    def track():
        if true_objective is None:
            cost = center_cost
        else:
            cost = true_objective(center)
        history.append(cost)

    track()

    for i in range(iterations):
        costs, deltas = populate(objective, center, popsize, sigma)
        costs = heuristic(costs, deltas, center_cost)

        flat_costs = costs[:, 0] - costs[:, 1]
        step = np.average(deltas / sigma * flat_costs[:, None], axis=0)
        center -= step_size * step
        center_cost = objective(center)

        track()

    return center, history

problem_size = 100
iterations = 1000
rotation, _ = np.linalg.qr(np.random.randn(problem_size, problem_size))
init = np.random.uniform(-5, 5, problem_size)

def ellipsoid_problem(x):
    return np.sum(10**(6 * np.linspace(0, 1, len(x))) * np.square(x))

def rotated_problem(x):
    return ellipsoid_problem(rotation @ x)

def noisy_problem(x):
    multiplicative_noise = np.random.uniform(0.707, 1.414)
    additive_noise = np.abs(np.random.normal(scale=50000))
    return rotated_problem(x) * multiplicative_noise + additive_noise

objective, true_objective = noisy_problem, rotated_problem
optimized, history = minimize(objective, init, iterations,
                              step_size=0.25, sigma=0.3,
                              true_objective=true_objective)

print(" " * 11 + "plot of log10-losses over time")
plot(np.log10(history), num_y_chars=23)

print("loss, before optimization: {:9.6f}".format(true_objective(init)))
print("loss,  after optimization: {:9.6f}".format(true_objective(optimized)))