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

Cmaes.scala

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)
}