python实现CMA-ES

CMA-ES（Covariance Matrix Adaptation Evolution Strategy）是一种基于进化算法的优化方法，用于寻找非线性、高维、非凸函数的最优解。下面是Python实现CMA-ES的示例代码：

import numpy as np
from scipy.stats import multivariate_normal

class CMAES:
    def __init__(self, dim, mu, sigma):
        self.dim = dim
        self.mu = mu
        self.sigma = sigma
        self.lambda_ = int(4 + np.floor(3 * np.log(dim)))  # population size
        self.weights = np.log(self.lambda_ + 1 / 2) - np.log(np.arange(1, self.lambda_ + 1))
        self.weights /= np.sum(self.weights)
        self.mu_eff = np.sum(self.weights) ** 2 / np.sum(self.weights ** 2)  # variance-effective size of mu
        self.c_sigma = (self.mu_eff + 2) / (dim + self.mu_eff + 5)  # learning rate for sigma control
        self.d_sigma = 1 + 2 * max(0, np.sqrt((self.mu_eff - 1) / (dim + 1)) - 1) + self.c_sigma  # damping parameter for sigma control
        self.pc = np.zeros(dim)  # evolution path for C
        self.ps = np.zeros(dim)  # evolution path for sigma
        self.B = np.eye(dim)  # transformation matrix
        self.D = np.ones(dim)  # diagonal matrix
        self.C = np.eye(dim)  # covariance matrix
        self.sigma_hist = []

    def ask(self):
        self.z = np.random.randn(self.lambda_, self.dim)
        self.y = np.dot(self.z, self.B * np.diag(self.D)) * self.sigma + self.mu
        return self.y

    def tell(self, x, fit):
        idx = np.argsort(fit)
        z = np.dot(self.B.T, self.z[idx, :].T).T / np.sqrt(self.D)
        self.pc = (1 - self.c_sigma) * self.pc + np.sqrt(self.c_sigma * (2 - self.c_sigma) * self.mu_eff) * np.sum(self.weights[:, None] * z[:self.mu, :], axis=0)
        self.ps = (1 - self.c_sigma) * self.ps + np.sqrt(self.c_sigma * (2 - self.c_sigma)) * np.dot(self.B, np.dot(np.diag(1 / np.sqrt(self.D)), np.mean(z[:self.mu, :], axis=0)))
        hsig = np.linalg.norm(self.ps) / np.sqrt(1 - (1 - self.c_sigma) ** (2 * self.iteration)) / self.d_sigma < 1.4 + 2 / (self.dim + 1)
        self.sigma *= np.exp((self.c_sigma / self.d_sigma) * (np.linalg.norm(self.ps) / self.sigma - 1))
        if hsig:
            self.C = (1 - self.c_sigma) * self.C + self.c_sigma * np.dot(self.pc[:, None], self.pc[None, :]) + self.c_sigma * (1 - np.sum(self.weights ** 2)) * self.C
        else:
            self.C = (1 - self.c_sigma) * self.C + self.c_sigma * np.dot(self.pc[:, None], self.pc[None, :])
        self.D, self.B = np.linalg.eigh(self.C)
        self.D = np.sqrt(self.D)
        self.sigma_hist.append(self.sigma)
        self.iteration += 1

    def result(self):
        return self.mu, self.sigma, self.sigma_hist

使用示例：

def sphere(x):
    return np.sum(x ** 2)

es = CMAES(dim=10, mu=5, sigma=1)
for i in range(100):
    x = es.ask()
    fit = np.array([sphere(x[i]) for i in range(len(x))])
    es.tell(x, fit)
print(es.result())

此示例演示了使用CMA-ES来优化10维球函数的最小值。