Simplistic failure case for AL-CMA-ES

[TGlas]

I want to solve the following problem: place $n$ points inside a 2D ball of radius $r$ so that the sum of (squared) distances to the nearest neighbor of each point is maximized. The following code is a direct adaptation of the
constraint handling notebook:

import cma
import numpy as np

# problem instance: place n points in 2D into a ball of radius r
n = 3
r = 3

# sum of smallest squared distance to nearest neighbor
def objective(x):
	total = 0
	for i in range(n):
		best = 1e100
		for j in range(n):
			if j == i: continue
			dist = np.linalg.norm(x[2*i:2*i+2] - x[2*j:2*j+2])
			best = min(best, dist**2)
		total += best
	return -total

# all points must be inside of the (closed) ball
def constraints(x):
	c = np.zeros(n)
	for i in range(n):
		c[i] = np.linalg.norm(x[2*i:2*i+2]) - r
	return c

x0 = np.random.randn(2*n)
print("initial:", constraints(x0), objective(x0))

x, es = cma.fmin_con2(objective, x0, 1, constraints)
print("optimized:", constraints(x), objective(x))

xx = r * np.array([1,0,-0.5,-np.sqrt(0.75),-0.5,np.sqrt(0.75)])
print("analytic:", constraints(xx), objective(xx))

For $n=3$ the problem has an analytic solution with an objective function value of $81$ (points placed on an equilateral triangle in a circle of radius 3). However, in many runs, AL-CMA-ES returns a much larger value, corresponding to a solution that is hopelessly infeasible. It is worth pointing out that going from squared distances to normal distances "solves" the problem; this the result is very close to the analytical optimum of $9 \cdot \sqrt{3}$. Also, changing the radius to $r=5$ seems to solve the problem.

Note that the problem is hard in the sense that the case distinction that is implicit in the minimum computation is active at the optimum. However, it is unclear to me why the problem is solved with the norm, while squaring the norm results in a breakdown of AL-CMA-ES. Also, it is off that the stability of the algorithm depends on the radius value.

This should probably not be tagged as a bug, but rather as a problem related to stability and parameter tuning. It would be nice to make the algorithm succeed on this simple problem.

[nikohansen]

FTR, when solving such a problem, I tend to try setting $r=\max_i ||p_i||$ and adding a mild penalty on the objective like $r^2 / 10$ to prevent divergence. EDIT: we may need $r^{2.5}$ here?

On the main point: I agree that this should not lead to a failure in practice. It should be relatively simple to get more insights by looking at the plots of the constraints logger.

First idea: the problem may well be that the objective is not convex: because the objective "decays" quadratically (grows negatively) but the constraint grows linearly, the Lagrangian coefficient changes may never be able catch the quickly growing difference between the objective and constraints gradients with increasing $\sigma$? If this is correct, the effect should go away with a constraint like $|p_i| - r + 1_{|p_i| > r} \times ||p_i| - r|^{2.2}$.

[nikohansen]

Indeed, I did a quick check, when I add

max(((np.linalg.norm(x[2*i:2*i+2]) - r)**2, 0))

in the constraints computation to the term np.linalg.norm(x[2*i:2*i+2]) - r, the divergence is gone (too).

[nikohansen]

FTR, a larger radius does not help when x0 and sigma0 are also changed proportionally. I suspect that being just long enough not too far from the feasible domain is enough to have the Lagrangian coefficients catch up.

We have a race between increasing $\sigma$ and the increasing coefficients. In principle, if we wait long enough, the $\mu_i$ coefficients should win the race, however the 'tolfacupx' termination gets in the way.

How can we prevent this in practice? Possible candidates are

• better initialization in particular of $\mu_i$ (and/or $\lambda_i$)
• convexify the constraints by default like above
• faster adaptation of $\mu_i$, e.g. based on the infeasibility ratio
• do not increase $\lambda_i$ beyond $\mu_i g_i$ or at all when all solutions are infeasible
• set 'stall_sigma_change_on_divergence_iterations': 2 to mitigate divergence (seems to have little effect and is problematic on a dynamically changing fitness function)
• ...?

[TGlas]

First idea: the problem may well be that the objective is not convex: because the objective "decays" quadratically (grows negatively) but the constraint grows linearly, the Lagrangian coefficient changes may never be able catch the quickly growing difference between the objective and constraints gradients with increasing σ ? If this is correct, the effect should go away with a constraint like | p i | − r + 1 | p i | > r × | | p i | − r | 2.2 .

Ideally, the quadratic gets linearized close to the optimum, but the problem is of course that we may be very far away.

We have a race between increasing σ and the increasing coefficients. In principle, if we wait long enough, the μ i coefficients should win the race, however the 'tolfacupx' termination gets in the way.

The race was exactly my hypothesis when I wrote that it may be a matter of tuning. However, tuning would likely only shift the failure point away, in contrast to offering a principled solution. On the other hand, convexification also has assumptions. Maybe the adaptation process itself needs a different rule, i.e., change the parameters only in case of sufficient evidence?

[nikohansen]

Ideally, the quadratic gets linearized close to the optimum

If the constraint is locally linear adding the quadratic term doesn't change that, I think.

Maybe the adaptation process itself needs a different rule, i.e., change the parameters only in case of sufficient evidence?

AFAICS, the adaptation is too slow (not too fast), so waiting for more evidence should make the problem worse.

If I "turn off" 'tolfacupx': 1e33, an unsuccessful run ends up in all zeros and looks like

A successful run looks like

It looks like the $\mu_i$ adaptation is

• either (way) too slow ($\times10^2$ in 200 iterations) compared to the $\lambda_i$ adaptation ($\times10^7$ in 50 iterations), at least in the beginning to get stably close to the feasible domain
• or somewhat badly initialized. EDIT: the problem seem to be rather the divergence after initialization which makes the (good) initial values moot.

This may be fixable.

Code:

import matplotlib.pyplot as plt
plt.figure(np.random.randint(1e8))

es.objective_function.al.loggers.plot()

[nikohansen]

1. The main "first problem" is that the step-size increases fast. This is an algorithm feature we probably don't want to give up. It seems desirable but difficult to moderate this behavior depending on diverging fitness.
2. AFAICS, only convexification of the constraints (or the objective) addresses the "root problem".
3. The adaptation rate of $\mu_i$ should most likely depend on the feasibility ratio. I can't see any good reason why it should not. This alleviates the issue, but in my first implementation attempts I still see initial divergence.