Last active
February 10, 2017 16:21
-
-
Save yukoba/8fe869f6a019bd097f20b4c6aef954cd to your computer and use it in GitHub Desktop.
最適化アルゴリズムCMA-ESの実装。MATLAB版をScalaに移植。
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
package jp.yukoba | |
import breeze.linalg.eigSym.EigSym | |
import breeze.linalg.{*, Axis, DenseMatrix, DenseVector, argsort, diag, eigSym, max, min, norm, sum} | |
import breeze.numerics.{log, pow, sqrt} | |
import breeze.stats.distributions.Gaussian | |
/** | |
* 最適化アルゴリズム CMA-ES の実装。 | |
* https://www.lri.fr/~hansen/purecmaes.m を Scala に移植。 | |
* https://arxiv.org/abs/1604.00772 も参考にしています。 | |
* | |
* 依存ライブラリ:Breeze 0.13 | |
*/ | |
object Cmaes { | |
def minimize(feval: DenseMatrix[Double] => DenseVector[Double], | |
N: Int, | |
xmean1: Option[DenseVector[Double]] = None, | |
sigma1: Option[Double] = None, | |
stopfitness1: Option[Double] = None, | |
stopeval1: Option[Int] = None, | |
lambdaFactor: Int = 1): Unit = { | |
var sigma = sigma1.getOrElse(1.0) | |
val stopfitness = stopfitness1.getOrElse(1e-10) | |
var xmean = xmean1.getOrElse(DenseVector.rand(N, Gaussian(0, 1))) | |
val stopeval = stopeval1.getOrElse(1000 * N * N) | |
// Strategy parameter setting: Selection | |
val lam = 4 + (3 * math.log(N)).toInt * lambdaFactor | |
val mu = lam / 2 | |
val weights1 = log((lam + 1) / 2d) - log(DenseVector.range(0, mu) + 1) | |
val weights = weights1 / sum(weights1) | |
val mueff = 1 / sum(weights *:* weights) | |
// Strategy parameter setting: Adaptation | |
val cc = (4 + mueff / N) / (N + 4 + 2 * mueff / N) | |
val cs = (mueff + 2) / (N + mueff + 5) | |
val alpha_cov = 2 | |
val c1 = alpha_cov / ((N + 1.3) * (N + 1.3) + mueff) | |
val cmu = math.min(1 - c1, alpha_cov * (mueff - 2 + 1d / mueff) / ((N + 2) * (N + 2) + alpha_cov * mueff / 2d)) | |
val damps = 1 + 2 * math.max(0, math.sqrt((mueff - 1) / (N + 1d)) - 1) + cs | |
// Initialize dynamic (internal) strategy parameters and constants | |
var pc = DenseVector.zeros[Double](N) | |
var ps = DenseVector.zeros[Double](N) | |
var B = DenseMatrix.eye[Double](N) | |
var D = DenseVector.ones[Double](N) | |
var C = B * diag(pow(D, 2)) * B.t | |
var invsqrtC = B * diag(pow(D, -1)) * B.t | |
var eigenval = 0 | |
val chiN = math.sqrt(N) * (1 - 1d / (4 * N) + 1d / (21 + N * N)) | |
var counteval = 0 | |
while (counteval < stopeval) { | |
// Generate and evaluate lambda offspring | |
val tmp = DenseMatrix.rand(N, lam, Gaussian(0, 1)).apply(::, *) *:* D | |
val arx = (sigma * (B * tmp)).apply(::, *) + xmean | |
val arfitnessOrig = feval(arx) | |
counteval += lam | |
// Sort by fitness and compute weighted mean into xmean | |
val arindex = argsort(arfitnessOrig) | |
val arfitness = arfitnessOrig(arindex) | |
val xold = xmean | |
xmean = arx(::, arindex.take(mu)).toDenseMatrix * weights | |
// Cumulation: Update evolution paths | |
ps = (1 - cs) * ps + math.sqrt(cs * (2 - cs) * mueff) / sigma * invsqrtC * (xmean - xold) | |
val hsig = if (sum(pow(ps, 2)) / (1 - math.pow(1 - cs, 2d * counteval / lam)) / N < 2 + 4d / (N + 1)) 1 else 0 | |
pc = (1 - cc) * pc + hsig * math.sqrt(cc * (2 - cc) * mueff) / sigma * (xmean - xold) | |
// Adapt covariance matrix C | |
val artmp = (1 / sigma) * (arx(::, arindex.take(mu)).toDenseMatrix.apply(::, *) - xold) | |
C = (1 - c1 - cmu) * C + | |
c1 * (pc * pc.t + (1 - hsig) * cc * (2 - cc) * C) + | |
cmu * (artmp * diag(weights) * artmp.t) | |
// Adapt step size sigma | |
sigma *= math.exp((cs / damps) * (norm(ps) / chiN - 1)) | |
// Update B and D from C | |
if (counteval - eigenval > lam / (c1 + cmu) / N / 10) { | |
eigenval = counteval | |
val EigSym(lambda, evs) = eigSym(C) | |
B = evs | |
D = sqrt(lambda) | |
invsqrtC = B * diag(pow(D, -1)) * B.t | |
} | |
println(s"$counteval: ${arfitness(0)}") | |
if (arfitness(0) <= stopfitness || max(D) > 1e7 * min(D)) { | |
counteval = stopeval | |
} | |
} | |
} | |
} | |
object CmaesTest extends App { | |
// fsphere | |
// Cmaes.minimize(x => sum(pow(x, 2), Axis._0).t, 12) | |
// Rosebrock | |
Cmaes.minimize(x => | |
((sum(pow(pow(x(0 to -2, ::), 2) - x(1 to -1, ::), 2), Axis._0) *:* 100d) + | |
sum(pow(x(0 to -2, ::) - 1d, 2), Axis._0)).t, 12, | |
lambdaFactor = 5) | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment