Implementation of CMA-ES (Covariance Matrix Adaptation - Evolution Strategy) by Ruby. It requires some gems, narray, gsl, gnuplot.

cma_es.rb

require "gnuplot"
require "./octavize"
require "./function"

class CmaEs
  def initialize
    @objective_function = Proc.new do |x|
      rastrigin10 x
    end

    # 探索変数の次元数
    @N = 10

    # テスト用初期平均ベクトル
    @xmean = NVector.random(@N)
    @sigma = 0.5
    @stopfitness = 1E-10
    @stopeval = 1E3 * @N**2

    # 戦略パラメータ：選択
    @lambda = 4 + (100 * log(@N)).truncate
    @mu = @lambda / 2.0
    @weights = NVector.new_fill(@mu.to_i, log(@mu + 1.0 / 2.0)) - NVector.by_log(@mu)
    @mu = @mu.truncate
    @weights = @weights / @weights.sum.to_f
    @mueff = @weights.sum**2 / @weights.to_a.inject(0.0){|sum, w| sum += w**2}

    # 戦略パラメータ：適応
    @cc = (4.0 + @mueff / @N) / (@N + 4.0 + 2.0 * @mueff / @N)
    @cs = (@mueff + 2.0) / (@N + @mueff + 5.0)
    @c1 = 2.0 / ((@N + 1.3)**2.0 + @mueff)
    @cmu = 2.0 * (@mueff - 2.0 + 1.0 / @mueff) / ((@N + 2.0)**2.0 + 2.0 * @mueff / 2.0)
    @damps = 1.0 + 2.0 * [0.0, sqrt((@mueff - 1.0) / (@N + 1.0)) - 1.0].max + @cs

    # 初期動的（内的）戦略パラメータと定数
    @pc = @ps = NVector.float(@N)
    @b = NMatrix.float(@N, @N).I
    @d = NVector.float(@N).fill(1.0)
    @c = @b * @d.each_power(2.0).to_diag * @b.transpose
    @invsqrtC = @b * @d.each_power(-1.0).to_diag * @b.transpose
    @eigeneval = 0.0;
    @chiN = @N**0.5 * (1.0 - 1.0 / (4.0 * @N) + 1.0 / (21.0 * @N**2.0))

    # @xmean = NVector.ref NArray.load("./xmean.dat")
    # @d_rand = NMatrix.ref NArray.load("./d_rand.dat")

    @path = Array.new(@N).map{|a| a = []}
  end

  attr_accessor :xmean, :path

  def gen_offspring
    # λ個の子孫の生成と評価
    @arx = NMatrix.float(@lambda, @xmean.total)
    @arfitness = NVector.float(@lambda)
    @lambda.times do |k|
      # @arx[k, 0..-1] = @xmean + @sigma * @b * @d_rand.col(k)
      @arx[k, 0..-1] = @xmean + @sigma * @b * @d.collect{|v| v *= randn}
      @arfitness[k] = @objective_function.call(@arx.col(k))
      @counteval += 1
    end
  end

  def update_mean
    # 適応度によるソートと重み付き平均による平均の更新
    @arindex = @arfitness.sort_index
    @arfitness = @arfitness.sort
    @xold = @xmean
    @xmean = @arx[@arindex[0..@mu-1], 0..-1] * @weights
  end

  def update_path
    # 集積：進化経路のアップデート
    @ps = (1 - @cs) * @ps +
      sqrt(@cs * (2 - @cs) * @mueff) * @invsqrtC * (@xmean - @xold) / @sigma
    @hsig = @ps.each_power(2).sum / (1.0 - (1.0 - @cs)**(2.0 * @counteval.to_f / @lambda.to_f)) / @N < (2.0 + 4.0 / (@N + 1.0)) ? 1 : 0
    @pc = (1 - @cc) * @pc +
      @hsig * sqrt(@cc * (2 - @cc) * @mueff) * (@xmean - @xold) / @sigma
  end

  def adapt_covmat
    # 共分散行列Cの適応
    artmp = (1 / @sigma) * (@arx[@arindex[0..@mu-1], 0..-1] - @xold.repmat(@mu))
    @c = (1 - @c1 - @cmu) * @c +
      @c1 * (@pc ** @pc + (1 - @hsig) * @cc * (2 - @cc) * @c) +
      @cmu * artmp * @weights.to_diag * artmp.transpose
  end

  def adapt_sigma
    # ステップサイズσの適応
    @sigma = @sigma * exp((@cs / @damps) * (@ps.norm / @chiN - 1.0))
  end

  def update_eigen
    # 行列aとbをcからアップデート
    if @counteval - @eigeneval > @lambda / (@c1 + @cmu) / @N / 10.0
      @eigeneval = @counteval
      # @c = @c.triu + @c.triu(1).transpose
      val, vec = @c.eigen
      @b = vec
      @d = val.each_power(0.5)
      @invsqrtC = @b * @d.to_diag.each_power(-1) * @b.transpose
    end
  end

  def generation_loop
    # 世代ループ
    @counteval = 0
    while @counteval < @stopeval
      gen_offspring
      update_mean
      update_path
      adapt_covmat
      adapt_sigma
      update_eigen

      if @arfitness[0] <= @stopfitness || @d.max > 1e7 * @d.min
        break
      end

      # 経過表示
      print "count: #{@counteval}\t-> "
      print sprintf("fit: %.5f\t%.5f\t%.5f\t%.5f\t", @arfitness[0], @sigma * sqrt(@c.diag.max), @d.max / @d.min, @lambda)
      print "xmean: "
      @xmean.puts true

      # パス記録
      @path.each_with_index{|p, i| p << @xmean[i]}
    end
  end
end

cma_es = CmaEs.new
cma_es.generation_loop
path = cma_es.path

x = (0..path[0].size - 1).to_a
y = {}
path.each_with_index do |xs, i|
  y["dim: #{i}"] = xs
end

def draw_chart(x, y)
  Gnuplot.open do |gp|
    Gnuplot::Plot.new(gp) do |plot|
      y.each do |name, value|
        if x.size == value.size
          plot.data << Gnuplot::DataSet.new([x, value]) do |ds|
            ds.with = "lines"
            ds.title = name
          end
        end
      end
    end
  end
end

draw_chart(x, y)

function.rb

require "./octavize"

def rosenbrock(x)
  if x.total < 2
    raise 'dimension must be greater one'
  end
  100 * ((x[0..-2].each_power(2) - x[1..-1]).each_power(2)).sum + (x[0..-2].each_minus(1).each_power(2)).sum
end

def rastrigin(x)
  n = x.total
  if n < 2
    raise 'dimension must be greater one'
  end
  b = x.to_a.inject(0.0){|sum, val| sum += val**2 - 10.0 * cos(2 * PI * val)}
  10 * n + b
end

def rastrigin10(x)
  n = x.total
  if n < 2
    raise 'dimension must be greater one'
  end
  scale = NVector.range((0..n - 1)).each_times(1.0 / (n - 1)).collect{|v| v = 10.0**v}
  10.0 * n + ((scale.times(x)).each_power(2.0) - 10.0 * (2.0 * PI * (scale.times(x))).each_cos).sum
end

def griewank(x)
  1.0 + x.each_power(2).sum / 4000.0 - x.total.times.inject(1.0){|prod, i| prod *= cos(x[i] / sqrt(i + 1))}
end

octavize.rb

require "narray"
require "gsl"
include Math

def randn
  sqrt(-2 * log(rand)) * cos(2 * PI * rand)
end

class NArray
  def each_power(n)
    self.collect do |e|
      e = e**n
      e.finite? ? e : 0.0
    end
  end

  def each_times(n)
    self.collect{|e| e *= n}
  end

  def each_minus(n)
    self.collect{|e| e -= n}
  end

  def each_cos
    self.collect{|e| cos(e)}
  end

  def repmat(m, n)
    col, row = self.shape
    out = NMatrix.float(col * m, row * n)
    n.times do |i|
      m.times do |j|
        out[col * j..col * j + col - 1, row * i..row * i + row - 1] = self
      end
    end
    out
  end

  def puts(t = false)
    print "#{self.shape} "
    print "\n" if not t
    self.to_a.each do |a|
      if a.instance_of?(Array)
        a.each do |b|
          print sprintf("% .5f", b) + "\t"
        end
      else
        print sprintf("% .5f", a) + "\t"
      end
      print "\n" if not t
    end
    print "\n"
  end

  def self.load(file, split = " ")
    mat = []
    File.open(file).each do |row|
      if row[0] != "#"
        r = row.split(split).map{|v| v.to_f}
        if not r.empty?
          val = r.size == 1 ? r[0].to_f : r
          mat << val
        end
      end
    end
    NArray.to_na(mat)
  end

  def save(file_name, delim = " ")
    File.open(file_name, "w") do |file|
      file.puts "# shape: #{self.shape}"
      self.to_a.each do |a|
        if a.instance_of?(Array)
          a.each do |b|
            file.print "#{b}#{delim}"
          end
        else
          file.print "#{a}#{delim}"
        end
        file.print "\n"
      end
    end
  end
end

class NMatrix
  def self.rand(cols, rows)
    NMatrix.float(cols, rows).randomn
  end

  def col(n)
    NVector.to_na(self[n, 0..-1].flatten.to_a)
  end

  def row(n)
    NVector.to_na(self[0..-1, n].flatten.to_a)
  end

  def repmat(m, n)
    col, row = self.shape
    out = NMatrix.float(col * m, row * n)
    n.times do |i|
      m.times do |j|
        out[col * j..col * j + col - 1, row * i..row * i + row - 1] = self
      end
    end
    out
  end

  def triu(k = 0)
    mat = self.clone
    mat.shape[0].times do |i|
      (i + 1 - k..mat.shape[1] - 1).each do |j|
        mat[i, j] = 0
      end
    end
    mat
  end

  def eigen
    # tmp = self.to_gm.eigen_symmv
    # [NVector.ref(tmp[0].to_na), NMatrix.ref(tmp[1].to_na)]
    u, v, s = self.to_gm.SV_decomp
    [NVector.ref(s.to_na), NMatrix.ref(v.to_na)]
  end

  def diag
    cols, rows = self.shape
    if cols == rows
      vec = NVector.float(cols)
      cols.times do |i|
        vec[i] = self[i, i]
      end
    end
    vec
  end

end

class NVector
  def self.random(dim)
    vec = NVector.float(dim)
    vec.each do |val|
      val = rand * 2.0 - 1.0
    end
    vec
  end

  def self.new_fill(dim, val)
    NVector.float(dim).fill(val)
  end

  def to_diag
    mat = NMatrix.float(self.total, self.total)
    self.total.times do |i|
      mat[i, i] = self[i]
    end
    mat
  end

  def part_power(range, n)
    x = self[range].each_power(n)
    self[range] = x
    self
  end

  def times(vector)
    i = 0
    self.collect{|v| v *= vector[i]; i += 1; v}
  end

  def repmat(n)
    dim = self.total
    out = NMatrix.float(n, dim)
    n.times do |i|
      out[i, 0..dim - 1] = self
    end
    out
  end

  def **(other)
    dim = self.total
    mat = NMatrix.float(dim, dim)
    dim.times do |i|
      dim.times do |j|
        mat[j, i] = self[i] * other[j]
      end
    end
    mat
  end

  def norm
    sqrt(self.to_a.inject(0.0){|sum, v| sum += v**2})
  end

  def self.by_log(dim)
    dim.to_i.times.inject(NVector.float(dim)){|v, i| v[i] = log(i + 1); v}
  end

  def self.range(range)
    vec = NVector.float(range.size)
    range.each_with_index do |r, i|
      vec[i] = r
    end
    vec
  end
end