CMA-ES 演化策略
CMA-ES（covariance matrix adaptation evolution strategy，Hansen 與 Ostermeier 1996～2001）是連續參數黑箱最佳化最常用的方法之一。它和遺傳演算法一樣，每一代產生一群候選解、留下好的，但不用交配與突變，而是維護一個多變量常態分布，並從好的樣本「學到」這個分布應該往哪裡移動、拉成什麼形狀。
演算法 Algorithm
第 \(g\) 代的分布是 \(\mathcal N(\mathbf m, \sigma^2 \mathbf C)\)：平均 \(\mathbf m\) 是目前的猜測、\(\sigma\) 是整體的步長、共變異矩陣 \(\mathbf C\) 決定形狀與方向。每一代：
- 取樣：\(\mathbf x_k = \mathbf m + \sigma \mathbf y_k\)，\(\mathbf y_k \sim \mathcal N(\mathbf 0, \mathbf C)\)，\(k = 1, \ldots, \lambda\)（預設 \(\lambda = 4 + \lfloor 3 \ln n \rfloor\)）。
- 排序：計算 \(f(\mathbf x_k)\)，只用名次，不用數值大小。
- 更新平均：前 \(\mu = \lambda/2\) 名的加權平均，名次越前權重越大： $$\mathbf m \leftarrow \mathbf m + \sigma \mathbf y_w, \qquad \mathbf y_w = \sum_{i=1}^{\mu} w_i \mathbf y_{i:\lambda}$$
- 更新共變異矩陣：讓「成功的步伐」以後更容易被抽到： $$\mathbf C \leftarrow (1 - c_1 - c_\mu)\mathbf C + c_1\, \mathbf p_c \mathbf p_c^\mathsf T + c_\mu \sum_{i=1}^{\mu} w_i\, \mathbf y_{i:\lambda} \mathbf y_{i:\lambda}^\mathsf T$$ 第二項（rank-1）用演化路徑 \(\mathbf p_c\)，也就是最近幾代 \(\mathbf y_w\) 的指數加權累積，記住「一直往同一個方向走」的資訊；第三項（rank-μ）用這一代的好樣本。
- 更新步長（累積步長調整）：另一條路徑 \(\mathbf p_\sigma\) 是 \(\mathbf C^{-1/2}\mathbf y_w\) 的累積，比較它和隨機走動時的期望長度： $$\sigma \leftarrow \sigma \exp\!\left(\frac{c_\sigma}{d_\sigma}\left(\frac{\lVert \mathbf p_\sigma \rVert}{\mathbb E\lVert \mathcal N(\mathbf 0, \mathbf I)\rVert} - 1\right)\right)$$ 連續幾代都往同一方向（路徑長）表示步子太小，就放大；方向互相抵消（路徑短）表示在最小值附近來回，就縮小。
所有學習率 \(c_1, c_\mu, c_\sigma, d_\sigma\) 都有依維度 \(n\) 推導出的預設值（見範例程式），使用者通常只需要給初始點與初始步長，幾乎不用調參。這是它和 GA、PSO 很大的不同。
學到的 \(\mathbf C\) 會接近函數在最小值附近 Hessian 矩陣的反矩陣（差一個常數倍），所以 CMA-ES 等於在不算任何導數的情況下，學會了類似牛頓法的「把問題轉成圓形」。另外它有兩個重要的不變性：只用名次，所以對 \(f\) 做任何遞增轉換（例如取對數）結果都一樣；座標軸旋轉也不影響它的表現。
動手看 Try It
下方動畫在網頁裡執行二維 CMA-ES（每代 6 個樣本）。在 Rosenbrock 函數上，可以看到橢圓沿著彎曲的山谷拉長、跟著轉彎；在 Rastrigin 函數上，常常停在離起點最近的某個區域最小值：
和 GA、PSO 比較 Benchmark
三個 10 維的測試函數，每個演算法都只能評估 10,000 次，各跑 5 個亂數種子。GA 與 PSO 使用常見的設定（見程式），CMA-ES 用預設值：
# CMA-ES 和遺傳演算法（GA）、粒子群（PSO）的比較：10 維的三個測試函數，各演算法都只能評估 10,000 次，每個情況跑 5 個亂數種子
import math
import numpy as np
class CMAES:
"""(μ/μ_w, λ)-CMA-ES，參數用 Hansen（2016）教學文的預設值"""
def __init__(self, x0, sigma0, rng):
n = self.n = len(x0)
self.lam = 4 + int(3 * math.log(n)); self.mu = self.lam // 2 # 每代樣本數 λ、選前 μ 個
w = np.array([math.log(self.mu + 0.5) - math.log(i + 1) for i in range(self.mu)])
self.w = w / w.sum(); me = self.mueff = 1 / np.sum(self.w ** 2)
self.cs = (me + 2) / (n + me + 5) # 步長的學習率與阻尼
self.ds = 1 + 2 * max(0.0, math.sqrt((me - 1) / (n + 1)) - 1) + self.cs
self.cc = (4 + me / n) / (n + 4 + 2 * me / n) # 共變異矩陣的學習率
self.c1 = 2 / ((n + 1.3) ** 2 + me)
self.cmu = min(1 - self.c1, 2 * (me - 2 + 1 / me) / ((n + 2) ** 2 + me))
self.chin = math.sqrt(n) * (1 - 1 / (4 * n) + 1 / (21 * n * n)) # E‖N(0, I)‖
self.m, self.sigma, self.rng = np.array(x0, float), sigma0, rng
self.C, self.ps, self.pc, self.g = np.eye(n), np.zeros(n), np.zeros(n), 0
self.D, self.B = np.ones(n), np.eye(n)
def ask(self):
Z = np.array([[float(self.rng.normal()) for _ in range(self.n)] for _ in range(self.lam)])
self.Y = Z @ np.diag(self.D) @ self.B.T # y ~ N(0, C)
return self.m + self.sigma * self.Y
def tell(self, fits):
n, ys = self.n, self.Y[np.argsort(fits, kind="stable")[:self.mu]]
yw = self.w @ ys
self.m = self.m + self.sigma * yw # 1. 平均移到好樣本的加權平均
self.ps = (1 - self.cs) * self.ps + math.sqrt(self.cs * (2 - self.cs) * self.mueff) * (self.B @ ((self.B.T @ yw) / self.D))
self.g += 1
hs = float(np.linalg.norm(self.ps) / math.sqrt(1 - (1 - self.cs) ** (2 * self.g)) < (1.4 + 2 / (n + 1)) * self.chin)
self.pc = (1 - self.cc) * self.pc + hs * math.sqrt(self.cc * (2 - self.cc) * self.mueff) * yw
self.C = ((1 - self.c1 - self.cmu) * self.C # 2. 共變異矩陣：rank-1（演化路徑）＋ rank-μ
+ self.c1 * (np.outer(self.pc, self.pc) + (1 - hs) * self.cc * (2 - self.cc) * self.C)
+ self.cmu * sum(wi * np.outer(y, y) for wi, y in zip(self.w, ys)))
self.C = (self.C + self.C.T) / 2
self.sigma *= math.exp((self.cs / self.ds) * (np.linalg.norm(self.ps) / self.chin - 1)) # 3. 步長：路徑太長就放大
d2, self.B = np.linalg.eigh(self.C); self.D = np.sqrt(np.maximum(d2, 1e-300))
def cma(f, n, budget, seed):
rng = np.random.default_rng(seed)
es = CMAES(rng.uniform(-5, 5, n), 2.0, rng)
best, ev = math.inf, 0
while ev + es.lam <= budget:
X = es.ask(); fits = np.array([f(x) for x in X]); ev += es.lam
best = min(best, fits.min()); es.tell(fits)
return best
def ga(f, n, budget, seed, pop=50):
"""實數編碼 GA：錦標賽選擇、BLX-0.5 交配（機率 0.9）、高斯突變（每個基因 1/n）、保留最好的 1 個"""
rng = np.random.default_rng(100 + seed); lo, hi = -5.0, 5.0
P = rng.uniform(lo, hi, (pop, n)); F = np.array([f(x) for x in P]); ev = pop
while ev + pop - 1 <= budget:
def pick():
i, j = rng.integers(pop, size=2)
return P[i] if F[i] < F[j] else P[j]
kids = [P[int(np.argmin(F))].copy()]
while len(kids) < pop:
a, b = pick(), pick()
if rng.random() < 0.9:
l, h = np.minimum(a, b), np.maximum(a, b); d = h - l
c = rng.uniform(l - 0.5 * d, h + 0.5 * d)
else:
c = a.copy()
c = np.clip(c + (rng.random(n) < 1 / n) * rng.normal(0, 0.1 * (hi - lo), n), lo, hi)
kids.append(c)
P = np.array(kids); F = np.r_[F.min(), [f(x) for x in P[1:]]]; ev += pop - 1
return float(F.min())
def pso(f, n, budget, seed, swarm=30, w=0.7298, c=1.49618):
"""標準 PSO：慣性 0.7298、c₁ = c₂ = 1.49618，速度限制在 0.2 倍範圍"""
rng = np.random.default_rng(200 + seed); lo, hi, vm = -5.0, 5.0, 0.2 * 10
X = rng.uniform(lo, hi, (swarm, n)); V = rng.uniform(-0.2, 0.2, (swarm, n)) * (hi - lo)
F = np.array([f(x) for x in X]); ev = swarm; Pb, Fb = X.copy(), F.copy()
while ev + swarm <= budget:
g = Pb[int(np.argmin(Fb))]
V = np.clip(w * V + c * rng.random((swarm, n)) * (Pb - X) + c * rng.random((swarm, n)) * (g - X), -vm, vm)
X = np.clip(X + V, lo, hi); F = np.array([f(x) for x in X]); ev += swarm
Pb[F < Fb], Fb[F < Fb] = X[F < Fb], F[F < Fb]
return float(Fb.min())
n = 10
Q, _ = np.linalg.qr(np.random.default_rng(0).normal(size=(n, n))) # 隨機旋轉
wts = 1e6 ** (np.arange(n) / (n - 1))
funcs = {
"Rosenbrock（彎曲的狹長山谷）": lambda x: float(np.sum(100 * (x[1:] - x[:-1] ** 2) ** 2 + (1 - x[:-1]) ** 2)),
"旋轉橢球（條件數 10⁶）": lambda x: float(np.sum(wts * (Q @ x) ** 2)),
"Rastrigin（很多區域最小值）": lambda x: float(10 * n + np.sum(x ** 2 - 10 * np.cos(2 * np.pi * x))),
}
print("10 維、每次最多評估 10,000 次、5 個種子；表中是找到的最小值（真正的最小值都是 0）")
for name, f in funcs.items():
print(name)
for alg, run in (("CMA-ES", cma), ("GA", ga), ("PSO", pso)):
r = np.array([run(f, n, 10000, s) for s in range(5)])
print(f" {alg:6s} 中位數 {np.median(r):10.3g} 最好 {r.min():10.3g} 最差 {r.max():10.3g}")⬇ 下載 cmaes_benchmark.py（Python 3＋numpy）
實際執行的輸出：
10 維、每次最多評估 10,000 次、5 個種子；表中是找到的最小值（真正的最小值都是 0）
Rosenbrock（彎曲的狹長山谷）
CMA-ES 中位數 6.18e-30 最好 0 最差 1.28e-28
GA 中位數 13.3 最好 10.3 最差 14
PSO 中位數 4.44 最好 3.51 最差 9.14
旋轉橢球（條件數 10⁶）
CMA-ES 中位數 1.77e-35 最好 5.8e-37 最差 3.71e-34
GA 中位數 5.56e+03 最好 1.95e+03 最差 8.39e+03
PSO 中位數 422 最好 220 最差 3.45e+03
Rastrigin（很多區域最小值）
CMA-ES 中位數 15.9 最好 5.97 最差 27.9
GA 中位數 19.9 最好 13.4 最差 23.4
PSO 中位數 5.97 最好 2.98 最差 16.9- Rosenbrock 與旋轉橢球：CMA-ES 都收斂到數值精度（\(10^{-28}\) 以下），GA、PSO 還差很遠。這兩個函數的困難在於變數彼此相關、各方向陡峭程度差很多，正是學習共變異矩陣能解決的。GA 的交配與突變都是沿著座標軸進行，旋轉後的問題對它特別難。
- Rastrigin：PSO 最好。這個函數有非常多規則排列的區域最小值，CMA-ES 的步長縮小得太快，就困在附近的坑裡。實務上的解法是重新啟動並加大族群（IPOP-CMA-ES），大族群會「看到」整體的碗形趨勢；本頁沒有示範。
沒有一種演算法在所有問題上都最好（no free lunch）。比較時也要注意，GA 與 PSO 的表現很依賴參數設定，換一組參數，結果可能不同。
應用：調 PID PID Tuning
把三個增益取對數、縮放到 \([0, 1]^3\)（\(K_p\) 在 0.1～31.6、\(K_i\) 與 \(K_d\) 在 0.01～10），從範圍中央、步長 0.3 出發，最小化 ITAE。模擬使用取樣 0.01 s 的數位 PID：微分只作用在量測值（避免設定值改變時的微分衝擊）並加一階濾波，控制量限制在 ±10，飽和時停止積分。每代 7 個樣本，評估 294 次後（5 個種子的中位數那一個）：
| \(K_p\) | \(K_i\) | \(K_d\) | ITAE | 超越量 | 2% 安定時間 | |
|---|---|---|---|---|---|---|
| Ziegler–Nichols | 4.80 | 2.65 | 2.18 | 5.20 | 52.7% | 9.66 s |
| CMA-ES 294 次 | 6.48 | 0.89 | 4.68 | 1.05 | 2.6% | 4.08 s |
最佳化的結果用了較大的 \(K_d\) 與小得多的 \(K_i\)：Ziegler–Nichols 法則的目標是「四分之一衰減」，本來就容許很大的超越量，只適合當起點。完整的程式與其他方法的比較在貝氏最佳化一頁。
注意：最佳化只對「模擬的受控體＋這個成本函數」是最好的。模型有誤差、有雜訊、有負載變化時，最佳增益不一定最好，實務上常在成本裡加入控制量、雜訊放大，或同時對好幾個模型（強健性）評估。