最適化アルゴリズムCMA-ESの実装。MATLAB版をPython 3に移植。

cmaes.py

#############################################################
# 最適化アルゴリズム CMA-ES の実装。
# https://www.lri.fr/~hansen/purecmaes.m を Python 3 に移植。
# https://arxiv.org/abs/1604.00772 も参考にしています。
#############################################################

import numpy as np


def minimize(feval, N=12, xmean=None, sigma=1.0, stopfitness=1e-10, stopeval=None):
    xmean = np.random.randn(N, 1) if xmean is None else xmean
    stopeval = 1e3 * N ** 2 if stopeval is None else stopeval

    # Strategy parameter setting: Selection
    lam = 4 + int(3 * np.log(N))
    lam *= 5  # 最小値に到達する確率を上げるために、独自に増やしました
    mu = int(lam / 2)
    weights = np.log((lam + 1) / 2) - np.log(np.arange(mu) + 1)
    weights /= weights.sum()
    mueff = (weights.sum() ** 2) / (weights ** 2).sum()

    # Strategy parameter setting: Adaptation
    cc = (4 + mueff / N) / (N + 4 + 2 * mueff / N)  # (56)
    cs = (mueff + 2) / (N + mueff + 5)
    alpha_cov = 2
    c1 = alpha_cov / ((N + 1.3) ** 2 + mueff)  # (57)
    cmu = min(1 - c1, alpha_cov * (mueff - 2 + 1 / mueff) / ((N + 2) ** 2 + alpha_cov * mueff / 2))  # (58)
    damps = 1 + 2 * max(0, np.sqrt((mueff - 1) / (N + 1)) - 1) + cs  # (55)

    # Initialize dynamic (internal) strategy parameters and constants
    pc = np.zeros([N, 1])
    ps = np.zeros([N, 1])
    B = np.eye(N)
    D = np.ones([N, 1])
    C = B @ np.diag((D ** 2).flatten()) @ B.T
    invsqrtC = B @ np.diag((D ** -1).flatten()) @ B.T
    eigenval = 0
    chiN = N ** 0.5 * (1 - 1 / (4 * N) + 1 / (21 + N ** 2))

    counteval = 0
    while counteval < stopeval:
        # Generate and evaluate lambda offspring
        arx = xmean + sigma * (B @ (D * np.random.randn(N, lam)))
        arfitness = feval(arx)
        counteval += lam

        # Sort by fitness and compute weighted mean into xmean
        arindex = np.argsort(arfitness)
        arfitness = arfitness[arindex]
        xold = xmean
        xmean = (arx[:, arindex[:mu]] @ weights).reshape([-1, 1])

        # Cumulation: Update evolution paths
        ps = (1 - cs) * ps + np.sqrt(cs * (2 - cs) * mueff) * invsqrtC @ (xmean - xold) / sigma
        hsig = (ps ** 2).sum() / (1 - (1 - cs) ** (2 * counteval / lam)) / N < 2 + 4 / (N + 1)
        pc = (1 - cc) * pc + hsig * np.sqrt(cc * (2 - cc) * mueff) * (xmean - xold) / sigma

        # Adapt covariance matrix C
        artmp = (1 / sigma) * (arx[:, arindex[:mu]] - xold)
        C = (1 - c1 - cmu) * C + \
            c1 * (pc @ pc.T + (1 - hsig) * cc * (2 - cc) * C) + \
            cmu * artmp @ np.diag(weights) @ artmp.T

        # Adapt step size sigma
        sigma *= np.exp((cs / damps) * (np.linalg.norm(ps) / chiN - 1))

        # Update B and D from C
        if counteval - eigenval > lam / (c1 + cmu) / N / 10:
            eigenval = counteval
            # C = np.triu(C) + np.triu(C, 1).T
            D, B = np.linalg.eigh(C, 'U')
            D = np.sqrt(D.real).reshape([N, 1])
            invsqrtC = B @ np.diag((D ** -1).flatten()) @ B.T

        if arfitness[0] <= stopfitness or D.max() > 1e7 * D.min():
            break

        if np.allclose(arfitness[0], arfitness[int(1 + 0.7 * lam)], 1e-10, 1e-10):
            sigma *= np.exp(0.2 + cs / damps)

        print("%d: %g" % (counteval, arfitness[0]))


# 最小化すべき評価関数。x.shape == (N, lam)
# minimize(lambda x: (x ** 2).sum(axis=0))  # fsphere
minimize(lambda x: 100 * ((x[:-1] ** 2 - x[1:]) ** 2).sum(axis=0) + ((x[:-1] - 1) ** 2).sum(axis=0))  # Rosebrock