fcmaes: Gradient-Free Optimization Library

The library provides several gradient-free optimization algorithms tailored for different problem scales and types:

• BiteOpt (C++ only): A robust single-objective optimizer. If you are unsure which algorithm to use, try this first. Use weighted penalties for constraints.
• MO-DE (Python/C++): Multi-objective optimization merging Differential Evolution and NSGA concepts. Supports constraints and mixed-integer problems.
• CMA-ES (Python/C++): Efficient for a low number of decision variables (< 500).
• CR-FM-NES (Python/C++): Fast Moving Natural Evolution Strategy. Recommended for high-dimensional problems (> 100 variables).
• PGPE (C++): Parameter Exploring Policy Gradients. Best for very high-dimensional problems (> 1000 variables) and machine learning tasks.
• Differential Evolution (Python/C++): Supports parallel function evaluation and mixed-integer problems.
• Quality Diversity (QD) Methods: Includes CVT-map-elites and the Diversifier meta-algorithm for exploring diverse niches in the search space.

When implementing objective functions in C++ to be called from Python (e.g., via ctypes) for optimization, follow these patterns to avoid performance bottlenecks:

• Avoid CLI overhead: Do not use a command-line interface to call the objective function from Python. This creates massive overhead, especially in parallel settings. Instead, use a shared library with a ctypes-based API.
• Ensure Reentrancy: Implement the objective function as reentrant so multiple invocations can run safely in parallel.
  • Avoid global variables: All shared state should be member variables of a singleton class instance.
  • Instance per call: Each reentrant API call should create its own singleton instance to avoid sharing state between parallel calls.
• Handle Discrete Variables correctly: Do not hide discrete decision variables from the optimizer. Use integer decision variables to generate discrete values so optimizers like fcmaes-MODE or DE can adapt their algorithms accordingly.
• Separate Initialization: Move expensive initialization code out of the objective function loop so it runs only once.

When choosing an algorithm for high-dimensional optimization in EvoJax, consider the following properties:

Performance Heuristics

• Low Optimization Budget: PGPE (and its C++ wrapper FPGPEC) tends to lead in the early stages of optimization.
• High Optimization Budget: CR-FM-NES (and its C++ wrapper FCRFMC) shows a higher improvement rate and often takes the lead as the budget increases.
• Strategy Tip: For complex tasks, consider a phased approach: use PGPE in the first phase, then initialize CR-FM-NES with the resulting population to leverage its superior late-stage convergence.

When using modecpp.retry, constraints are treated as a special kind of objective. To implement this:

1. Define the number of constraints in your fitness class using self.ncon.
2. In your objective function, return the constraint violation (e.g., a count of failed operations) as the last element of the returned array.
3. Filter the results after optimization to ensure only valid solutions (where the constraint value is 0) are kept.

For single-objective optimization, constraints are typically handled by assigning them a very high weight in a weighted sum of objectives.

def retry_modecpp(fit, retry_num = 32, popsize = 48, max_eval = 500000, 
                  nsga_update = True, logger = logger(), workers=mp.cpu_count()):
    xf, yf = modecpp.retry(fit.fun, fit.nobj, fit.ncon, fit.bounds, retry_num, popsize, 
                  max_evaluations = 960000, nsga_update = nsga_update, logger = logger, workers=workers)
    xs = []; ys = []
    for i in range(len(yf)):
        if yf[i][-1] == 0: # filter valid solutions (constraint == 0)
            xs.append(xf[i])
            ys.append(yf[i][:fit.nobj])
    return np.array(xs), np.array(ys)

In applications like laser fusion reactor design, points (lasers) must be arranged symmetrically, and weights (intensities) must be non-negative and ideally consist of only a few distinct values to simplify hardware requirements.

To achieve this, you can reduce the search space by:

1. Generating Orbits: Instead of optimizing all $N$ points, optimize a smaller set of 'reference points' and use symmetry operations (permutation and mirroring) to generate the full set of 12 symmetric points per reference point.
2. Reducing Decision Variables: This significantly shrinks the dimensionality of the optimization problem, speeding up convergence.
3. Weight Constraints: The optimization naturally tends toward solutions with a small number of distinct weights (e.g., only 2 or 3 different intensity levels), which is highly practical for physical implementation.

def generate_orbit(A, B, C):
    # apply 12 symmetries (permutations and sign changes)
    cart = []
    for perm in [(A, B, C), (B, C, A), (C, A, B)]:
        for signs in [(1, 1, 1), (1, -1, -1), (-1, 1, -1), (-1, -1, 1)]:
            x, y, z = [sign * coord for sign, coord in zip(signs, perm)]
            cart.append((x, y, z))
    return cartesian_to_spherical(cart)

# The fitness function then uses these generated orbits
def fit(x):
    points = generate_orbits(N, x[:2*n])
    weights = x[2*n:]
    return t_design_error(points, l_max, weights)

When implementing or verifying Spherical T-Designs (specifically for computing Legendre functions and Spherical Harmonics), the following mathematical steps are required:

1. Legendre Functions: Compute $\cos(\theta)$ for the Legendre functions.
2. Normalization: Compute the normalization factor.
3. Associated Legendre Functions: Compute the associated Legendre function for order $m$ and degree $l$.
   • For $m = 0$, no extra trigonometric factor is needed.
   • For $m > 0$, compute both the cosine and sine components.
4. Spherical Harmonics (SH) Calculation:
   • Step 1: Compute the SH values at all points (unweighted).
   • Step 2: Multiply each $Y_{l,m}$ by the corresponding weight and sum them up (element-wise multiply by $w_i$ and sum).
5. Coordinate Assumptions: Assume points are provided as pts[:,0]=theta and pts[:,1]=phi.
6. Weighting: For a given $l$, create a multiplier $1 / (2l + 1)$ for $l=0, ext{까지 } l_{max}$.

Follow these conventions for consistent optimization:

General Rules:

• Minimization: All objectives are minimized. Negate quantities that must be maximized.
• Constraint Format: For MODE/moretry, use [objective_0, ..., objective_m, constraint_0, ..., constraint_k]. A constraint is feasible if constraint <= 0.
• Constraint Conversion: Convert h(x) = 0 to abs(h(x)) - tolerance <= 0.
• Penalty Method: For scalar optimizers, use decoding/repair for hard structure and a scaled penalty for violations: objective + rho * sum(max(0, g_i(x))**2). Normalize terms before tuning rho.
• Rejection: Return a large finite rejection value for invalid simulations.

Variable Encoding:

• Integers: Encode as real coordinates; round, clamp, and repair in the objective. Optionally set DE or MODE's integer mask.
• Categories: Encode as an index into a fixed list. Do not interpret distances between unrelated categories physically.
• Complex Structures: Decode random keys or repair permutations/schedules/subsets deterministically.