共分散行列適応進化戦略(CMA-ES)

MetaTrader 5 — トレーディング | YYYY年06月DD日, 11:41

Andrey Dik

内容

はじめに

最適化の世界には多くのアルゴリズムが存在しますが、本記事では自動売買ロボットの最適化問題を解くために、最も強力な手法を探します。CMA-ES (Covariance Matrix Adaptation Evolution Strategy)は、そのような希少な例の一つであり、数学的厳密性と生物学的直観が融合し、単なる探索アルゴリズムではなく、問題構造そのものを捉えるアルゴリズムへと発展しています。

CMA-ESの歴史は1990年代後半のドイツの研究室に始まります。Nikolaus HansenとAndreas Ostermeierは、「解を単に探索するのではなく、問題の幾何構造に適応する最適化アルゴリズムを構築できるか」という根本的な問いを立てました。従来の進化的アルゴリズムは、親個体の周囲の球状の領域から子個体を生成していました。この方法は単純な関数には有効でしたが、条件の悪い複雑な問題では非効率であることが明らかになりました。ここで、この興味深いアルゴリズムを見ていきます。

アルゴリズムの実装

不規則な形状の島で宝探しをする場面を想像してください。通常の方法では、島が円形であるかのようにあらゆる方向へ均等に探索します。一方CMA-ESは、徐々に島の形状を学習し、宝が見つかる確率の高い方向へ探索を移行していきます。さらに、成功した経路を記憶し、その記憶をもとに将来の探索を計画します。

CMA-ESは、一見単純な方程式x_k ~ N(m, σ²C)に基づいています。しかし、この単純さの背後には深い数学的構造が存在します。この式の各記号は探索状態に関する重要な情報を持っています。mは最適解に対する現在の最良推定、σは既知領域からどれだけ大胆に離れるかの尺度、Cは問題の幾何構造を表現する共分散行列です。本記事における唯一の正当化された変更は、正規分布を冪分布に置き換えることです。すなわち、修正式 xₖ ~ PowerDist(m, σ²C) を用います。この変更は探索空間における挙動を変化させ（より広い「ジャンプ」を許容します）、しかしアルゴリズムの本質的な適応性は維持されます。

共分散行列Cはアルゴリズムの核心です。初期状態では単位行列として、球状分布を表しています。しかし反復ごとに進化し、改善が速い方向には伸び、進展が遅い方向には縮みます。結果として、分布は球から楕円へ、さらに細長い楕円体へと変形し、最終的には目的関数の等高線に沿う形へと整列していきます。

CMA-ESの主要な革新は進化経路の概念です。これはアルゴリズムの遺伝的記憶のようなものであり、成功した点の位置だけでなく、そこに到達するまでの過程そのものを記憶します。この経路は成功した連続的ステップの情報を蓄積し、最も有望な探索方向を示すベクトルを形成します。もう一つの進化経路はステップサイズ制御を担います。連続するステップが同じ方向に相関していれば、アルゴリズムは正しい方向に進んでいると判断しステップサイズを増加させます。逆にステップがランダムで無相関であれば、最適解付近に到達している可能性があるためステップサイズを減少させ、より慎重に探索します。

CMA-ESの背後にある本質的原理は、生物学的比喩ではなく最大尤度推定です。アルゴリズムは常に「成功した点を最も高い確率で生成する分布パラメータは何か？」という問いを解いています。これにより進化的最適化はヒューリスティックから統計的基盤を持つ手法へと変換されます。

共分散行列Cの更新は2つの成分から構成されます。1つは進化経路に基づくランク1更新、もう1つは選択された集団全体の情報を用いるランクμ更新です。ランク1更新は長期的傾向の安定的な捕捉を可能にし、ランクμ更新は新しい世代の情報に迅速に適応する役割を果たします。

CMA-ESでは改良されたヘヴィサイド関数が用いられ、アルゴリズムの停滞検出機構として重要な役割を果たします。この関数は進化経路の長さを期待値と比較し、短すぎる場合には「よろめき」の兆候と判断します。 hsig = 0（無効化条件）では、アルゴリズムがランダムに漂い、ステップが相殺されている状態を意味し、ステップサイズが過大である可能性があります。この場合ランク1更新の影響は減少し、集団全体の情報により依存します。 hsig = 1（有効化条件）では、探索が一方向に進展しており、ステップ間に相関が存在するため、現在のステップサイズは適切であると判断されます。

CMA-ESの最も深い性質の一つは、探索空間の変換に対する不変性です。このアルゴリズムは関数f(x)とf(Ax + b)を同等に扱うことができ、ここでAは非特異行列、bはシフトベクトルです。これは座標系の選択に依存しないことを意味します。この不変性は偶然ではなく、最大尤度原理と共分散行列の適応から自然に導かれる性質です。アルゴリズムは問題に対して自然な座標系を自動的に発見し、その軸を関数の主方向に一致させます。

理論的な美しさは実用性と両立する必要があります。CMA-ESは共分散行列の記憶にO(n²)、固有値分解にO(n³)の計算量を必要とし、数百次元程度までの問題に適用可能です。高次元問題に対しては、sep-CMA-ES（対角近似）、VkD-CMA（可変次元）、LM-CMA（リミテッドメモリ）といった拡張が提案されています。これらの詳細な実装については本稿の範囲外ですが、今後の議論で再び取り上げる可能性があります。

cmaesアルゴリズム

図1：CMA-ESアルゴリズムの図解

図は世代を通じた進化を示しており、アルゴリズムの3つの状態スナップショット（世代1、5、15）を通して、集団が最適解へと収束していく様子を表しています。

共分散行列の適応では、青い楕円が分布形状の変化を示します。世代1では球状（単位行列）、世代5では伸長・回転し、世代15では局所的な幾何構造に精密に適合しています。
アルゴリズムの構成要素として、赤点は全個体（λ子孫）、緑点は上位解（μ親）、黒点mは集団平均、破線円は目的関数の等高線を示します。

図の下部には主要な数式が示されており、冪分布への修正に関する注記も付されています。この図は、アルゴリズムが探索戦略を関数の地形に適応させ、探索空間を徐々に絞り込みながら最適解へ接近していく様子を明確に示しています。

それでは、擬似コードを準備します。

初期化

アルゴリズムパラメータの設定
- 集団サイズ(λ) = 50
- 親個体数(μ) = 25
- ランク1更新の学習率(c1) = 0.01
- ランクμ更新の学習率(cμ) = 0.8
- ステップサイズ減衰係数 = 0.6
- 初期ステップサイズ(σ) = 0.3
再結合重みの計算
- μ個の親それぞれについて、重みを log(μ + 0.5) − log(i + 1)として計算する
- 重みの合計が1になるよう正規化する
- μの有効選択質量(μ_eff)を計算する
戦略パラメータの初期化
- 進化経路の学習率（cs、cc）を計算する
- 標準正規分布ベクトルの期待ノルム(chiN)を計算する
- 固有値分解の実行間隔を決定する
データ構造の初期化
- C共分散行列 = 単位行列
- 固有ベクトル行列B = 単位行列
- 固有値ベクトルD = 単位ベクトル
- 進化経路pc、ps = ゼロベクトル
- 集団平均m = 探索空間内のランダム点
初期集団の生成
- 探索空間内でλ個のランダム点を生成する

進化の基本サイクル

停止条件が満たされるまで繰り返します。

ステップ1：子個体の生成

各λ個の子個体について：
- 標準正規分布からランダムベクトルzを生成する
- 変換を適用する：y = B × D × z
- 子個体を生成する：x_k = m + σ × y
- 探索空間の制約をx_kに適用する

ステップ2：評価と更新

全子個体の適応度を評価する
集団をソートする
- 子個体を適応度の降順にソートする
- 必要に応じて最良解を更新する
集団平均の更新
- 以前の平均をm_oldとして保存する
- 上位μ個の重み付き平均として新しい平均を計算する：m_new = Σ(w_i × x_i), i = 1..μ
進化経路の更新
- 平均の移動量を計算する：y_w = (m_new - m_old) / σ
- C^(-1/2) × y_wを固有値分解を用いて計算する
- ステップサイズ用進化経路：ps = (1 - cs) × ps + √(cs × (2 - cs) × μ_eff) × C^(-1/2) × y_wでpsを更新する
- ヘヴィサイド関数による停滞条件をチェックする
- 共分散行列用進化経路：pc = (1 - cc) × pc + stagnation_indicator × √(cc × (2 - cc) × μ_eff) × y_wによってpcを更新する
共分散行列の更新
- pcの外積からランク1更新を構築する
- 上位個体の偏差からランクμ更新を構築する
- Cを更新する：C = (1 - c1 - cμ) × C + c1 × (pc × pc^T) + cμ × Σ(w_i × y_i × y_i^T)
- 行列の対称性を保証する
ステップサイズの更新
- ps経路長に基づき適応係数を計算する
- σを更新する：σ = σ × exp((cs/damps) × (||ps||/chiN - 1))
- 数値安定性のためσを適切な範囲に制限する
固有値分解（必要時）
- 一定回数の反復ごとに：
  - C = B × D² × Bᵀのヤコビ分解を実行する
  - 固有値と固有ベクトルを降順にソートする
  - 行列の正定性を保証する
共分散行列のチェックと補正
- 定期的にCの正定性をチェックする
- 必要に応じて対角成分に微小値を加える
- 行列の対称性を保証する

ここから、C_AOクラスを継承したC_AO_CMAESを設計し、CMA-ES最適化アルゴリズムを実装します。

公開フィールドとメソッド

SetParams()：params配列からCMA-ESのパラメータを設定する
Init ()：アルゴリズムを初期化し、各変数の最小値・最大値・ステップサイズおよびエポック数を受け取る
Moving ()：新しい解を生成するメインループを実装する
Revision ()：新しい解の評価とアルゴリズム状態の更新をおこなう

非公開フィールド

CMA-ESパラメータ：CMA-ESアルゴリズムの特定のパラメータを保持する変数群：親個体数(mu)、ステップサイズ(sigma)、学習率（learningRateC1、learningRateCMu）、減衰係数(stepSizeDamping)、およびアルゴリズムの動作に必要なその他のパラメータ。
CMA-ESデータ構造：再結合重み、共分散行列covMatrix、固有ベクトル行列B、固有値ベクトルD、進化経路（pc、ps）を格納する配列群、およびその他の補助データ。mu_eff、counteval、eigeneval、chiN、eigenIntervalといった変数は、アルゴリズムの進行を制御・管理するために使用される。
性能最適化用変数：アルゴリズム実行中の計算を高速化するための値であり、キャッシュされた学習率、減衰係数、および繰り返し使用される事前計算済み定数などを含む。
補助配列：各種線形代数演算（サンプリング、変換、選択、共分散更新など）の過程で一時的にベクトルや行列を保存するために使用される構造。
補助メソッド：CMA-ESアルゴリズムの各ステップを実装する関数群であり、分布の初期化、分布パラメータの更新、固有値分解の計算、集団のソート、再結合重みの計算、平均ベクトルの更新、共分散行列の正定性チェックおよび強制的な正定性の保証などをおこなう。

全体として、C_AO_CMAESクラスはCMA-ESアルゴリズムのロジックをカプセル化しており、初期化、パラメータ調整、および最適化実行のためのインターフェースを提供します。このクラスには、アルゴリズムの実装に必要なパラメータ、データ配列、および各種メソッドが含まれています。また、公開フィールドと非公開フィールドを分離することで、クラス内部コンポーネントへのアクセス制御を実現しています。

//————————————————————————————————————————————————————————————————————
class C_AO_CMAES : public C_AO
{
  public: //----------------------------------------------------------
  ~C_AO_CMAES () { }
  C_AO_CMAES ()
  {
    ao_name = "CMAES";
    ao_desc = "Covariance Matrix Adaptation Evolution Strategy";
    ao_link = "https://www.mql5.com/ja/articles/18227";

    // Default parameters
    popSize         = 50;          // Default population size (lambda)
    mu              = 25;          // Number of parents (half of the population)
    learningRateC1  = 0.01;        // Learning rate for rank-1 update
    learningRateCMu = 0.8;         // Learning rate for rank-μ update
    stepSizeDamping = 0.6;         // Damping for step size

    // Create and initialize the parameters array
    ArrayResize (params, 5);
    params [0].name = "popSize";         params [0].val = popSize;
    params [1].name = "mu";              params [1].val = mu;
    params [2].name = "learningRateC1";  params [2].val = learningRateC1;
    params [3].name = "learningRateCMu"; params [3].val = learningRateCMu;
    params [4].name = "stepSizeDamping"; params [4].val = stepSizeDamping;
  }

  void SetParams ()
  {
    popSize         = (int)params [0].val;
    mu              = (int)params [1].val;
    learningRateC1  = params      [2].val;
    learningRateCMu = params      [3].val;
    stepSizeDamping = params      [4].val;
  }

  bool Init (const double &rangeMinP  [],  // minimum values
             const double &rangeMaxP  [],  // maximum values
             const double &rangeStepP [],  // step size
             const int     epochsP = 0);

  void Moving   ();
  void Revision ();

  //------------------------------------------------------------------
  private: //---------------------------------------------------------
  // CMA-ES specific parameters
  int mu;                  // Number of parents (selected points)
  double sigma;            // Step size
  double learningRateC1;   // Learning rate for rank-1 update
  double learningRateCMu;  // Learning rate for rank-μ update
  double stepSizeDamping;  // Damping factor for step size update

  // CMA-ES specific data structures
  double weights   [];      // Recombination weights
  double covMatrix [];      // Covariance matrix (stored as a one-dimensional array)
  double B  [];             // C eigenvectors
  double D  [];             // C eigenvalues (square roots)
  double pc [];             // Evolution path for C
  double ps [];             // Evolution path for step-size control
  double mu_eff;            // Effective selection mass
  int    counteval;         // Counter of function evaluations since the last decomposition
  int    eigeneval;         // Generation counter when decomposition was performed
  double chiN;              // Expected norm N(0,I)
  int    eigenInterval;     // Interval for eigenvalue decomposition

  // Variables for performance optimization
  double cs;                // Learning rate for the sigma path
  double cc;                // Learning rate for rank-1 path 
  double damps;             // Damping for sigma
  double hsig_threshold;    // Threshold for the Heaviside function

  // Auxiliary arrays
  double y_vec     [];         // Mutation vector
  double arindex   [];         // Array of indices for sorting
  double arfitness [];         // Fitness array for sorting
  double temp_vec  [];         // Temporary vector for matrix operations
  double invsqrtC_times_yw []; // Temporary storage for C^(-1/2) * y_w

  // Caching variables for Box-Muller
  double cached_normal;
  bool   has_cached;

  // Auxiliary methods
  void   InitDistribution          ();
  void   UpdateDistribution        ();
  void   ComputeEigendecomposition ();
  double GetChiN                   ();
  void   SortPopulation            ();
  void   ComputeWeights            ();
  void   UpdateMean                ();
  bool   CheckPositiveDefinite     ();
  void   EnforcePositiveDefinite   ();
};
//————————————————————————————————————————————————————————————————————

C_AO_CMAESクラスのInitメソッドは、CMA-ESアルゴリズムを初期化します。このメソッドは、各パラメータの最小値・最大値、ステップサイズ、およびエポック数を入力として受け取ります。まず、標準的な初期化を実行するためにStandardInitメソッドが呼び出されます。続いて、ボックスミュラー法による乱数生成にキャッシュされた値が存在しないことを示すため、has_cachedフラグがfalseに初期化されます。その後、ステップサイズsigmaが初期値0.3（探索範囲の30%）に設定されます。さらに、再結合重みを計算するためにComputeWeightsメソッドが呼び出され、標準正規分布N(0, I)の期待ノルムを計算するためにGetChiNメソッドが実行されます。

次に、mu_eff（有効選択質量）および座標数(coords)に基づいて、cs、cc、damps、hsig_thresholdといったパラメータが計算されます。これらのパラメータは、CMA-ES戦略を実行中に適応させるために使用されます。また、固有値分解を実行する頻度を決定するeigenIntervalが計算・設定されます。このパラメータは、アルゴリズムの性能を最適化するために調整されます。

続いて、CMA-ES専用の配列に対してメモリが確保されます。これには、共分散行列covMatrix、固有ベクトル行列B、固有値ベクトルD、進化経路pcおよびps、さらにy_vec、arindex、arfitness、temp_vec、invsqrtC_times_ywなどの補助配列が含まれます。これらの配列は、問題の次元数を変更する必要がある場合にのみ再生成されます。

その後、進化経路配列pcおよびpsはゼロで初期化され、共分散行列covMatrixと固有ベクトル行列Bは単位行列として初期化されます。また、固有値配列Dはすべて1で初期化されます。続いて、初期分布を設定するためにInitDistributionメソッドが呼び出されます。その後、countevalおよびeigenevalの評価カウンタがリセットされます。さらに、適応度値の再計算が必要であることを示すため、revisionフラグがtrueに設定されます。初期化が正常に完了した場合、このメソッドはtrueを返します。

要約すると、InitメソッドはCMA-ESアルゴリズムの動作に必要なすべての初期化処理を実行します。これには、各種パラメータの設定、配列用メモリの確保、および初期値の設定が含まれます。

//————————————————————————————————————————————————————————————————————
bool C_AO_CMAES::Init (const double &rangeMinP [],  // minimum values
                       const double &rangeMaxP [],  // maximum values
                       const double &rangeStepP [], // step size
                       const int     epochsP = 0)   // number of epochs
{
  if (!StandardInit (rangeMinP, rangeMaxP, rangeStepP)) return false;

  //------------------------------------------------------------------
  // Initialize Box-Muller caching
  has_cached = false;

  // Initialize CMA-ES specific variables
  sigma          = 0.3;  // Initial step size (30% of search range)

  // Calculate the effective mass of the variance selection
  ComputeWeights ();

  // Expected norm N(0,I)
  chiN = GetChiN ();

  // Calculate and save strategy parameters
  cs = (mu_eff + 2.0) / (coords + mu_eff + 5.0);
  cc = (4.0 + mu_eff / coords) / (coords + 4.0 + 2.0 * mu_eff / coords);
  damps = 1.0 + 2.0 * MathMax (0.0, MathSqrt ((mu_eff - 1.0) / (coords + 1.0)) - 1.0) + cs;
  hsig_threshold = 1.4 + 2.0 / (coords + 1.0);

  // Set the eigenvalue decomposition interval - tuning for performance
  eigenInterval = (int)(coords / (10.0 * MathSqrt (learningRateC1 + learningRateCMu)));
  eigenInterval = MathMax (1, eigenInterval);

  // Allocate arrays only once
  ArrayResize (covMatrix, coords * coords);
  ArrayResize (B, coords * coords);
  ArrayResize (D, coords);
  ArrayResize (pc, coords);
  ArrayResize (ps, coords);
  ArrayResize (y_vec, coords);
  ArrayResize (arindex, popSize);
  ArrayResize (arfitness, popSize);
  ArrayResize (temp_vec, coords);
  ArrayResize (invsqrtC_times_yw, coords);

  // Initialize evolution paths with zeros
  ArrayInitialize (pc, 0);
  ArrayInitialize (ps, 0);

  // Fast initialization of the identity covariance matrix and decomposition
  ArrayInitialize (covMatrix, 0.0);
  ArrayInitialize (B, 0.0);

  for (int i = 0; i < coords; i++)
  {
    covMatrix [i * coords + i] = 1.0;
    B [i * coords + i] = 1.0;
    D [i] = 1.0;
  }

  // Initialize initial distribution
  InitDistribution ();

  // Reset calculation counters
  counteval = 0;
  eigeneval = 0;

  // Forced fitness recalculation
  revision = true;

  return true;
}
//————————————————————————————————————————————————————————————————————

C_AO_CMAESクラスのMovingメソッドは、集団内に新しい子個体（解）を生成する役割を担います。 y = B * D * zの変換：集団全体を処理するループの内部には、各個体のすべての座標を走査する入れ子ループがあります。各座標iにおいて、y_vec[i]の値が計算されます。これは行列積によって実行されます。実際には、乱数ベクトルzと行列B × Dの積を計算しています。ここで、Bは共分散行列の固有ベクトル、Dは共分散行列の固有値、そしてzはu.PowerDistribution関数によって生成された乱数です。

子個体の生成：各子個体kの各座標iについて、a[k].c[i]の値が計算されます。これは、現在の平均値cB[i]（現在の分布の中心）に、スケーリングされたベクトルy_vec[i]を加算することでおこなわれます。スケーリングは、y_vec[i] にステップサイズsigmaを乗算することで実現されます。したがって、a [k].c [i] = cB [i] + sigma * y_vec [i]が得られます。

u.SeInDiSp関数は、生成された子個体の座標値が rangeMin[i]、rangeMax[i]、rangeStep[i] によって定義された有効範囲内に収まることを保証します。この関数は、値を境界で切り詰めたり、境界で反射させたり、あるいは最も近い有効なステップ値へ量子化したりすることができます。すべての子個体が生成されると、revisionフラグがtrueに設定されます。これは、生成された子個体の適応度（品質）を再計算する必要があることを示します。

結果として、このメソッドは現在の分布（平均値cBおよびBとDによって表現される共分散構造）とステップサイズsigmaに基づいて、新しい解集団を生成します。探索空間を調査するために乱数を用いて現在の平均周辺に小さな変動を作り出し、新たな候補解を生成します。また、SeInDiSp関数によって新しい解が許容範囲内に維持され、revisionフラグによってその後に適応度評価が実行されることが保証されます。

//————————————————————————————————————————————————————————————————————
void C_AO_CMAES::Moving ()
{
  // Generate new lambda descendants
  for (int k = 0; k < popSize; k++)
  {

    // Apply the transformation y = B*D*z
    for (int i = 0; i < coords; i++)
    {
      y_vec [i] = 0.0;
      for (int j = 0; j < coords; j++)
      {
        y_vec [i] += B [i * coords + j] * D [j] * u.PowerDistribution (0.0, -8, 8, 20);
      }

      // Generate a descendant: x_k = m + σ * y
      a [k].c [i] = cB [i] + sigma * y_vec [i];
      a [k].c [i] = u.SeInDiSp (a [k].c [i], rangeMin [i], rangeMax [i], rangeStep [i]);
    }
  }

  // Fitness recalculation mark
  revision = true;
}
//————————————————————————————————————————————————————————————————————

C_AO_CMAESクラスのRevisionメソッドは、新しい集団の適応度評価が完了した後に、アルゴリズムの各種パラメータを更新する役割を担います。まず、このメソッドはrevisionフラグの値を確認します。続いて、SortPopulation()メソッドが呼び出され、現在の集団が適応度に基づいてソートされます。その結果、最も優れた個体が集団の先頭に配置されます。次に、UpdateDistribution()メソッドが実行されます。このメソッドは、集団内の最良個体の情報を利用して、CMA-ESの分布パラメータを更新します。更新対象には、戦略ベクトルの平均値cB、共分散行列、sigmaステップサイズ、および進化経路が含まれます。分布の更新はCMA-ESにおける最も重要な処理の一つです。なぜなら、この更新によってアルゴリズムは得られた探索結果を学習し、それに基づいて探索戦略を適応的に変化させることができるからです。

結果として、RevisionメソッドはCMA-ESの主要な適応サイクルを構成します。このメソッドは、集団のソート、最良個体に基づく分布パラメータの更新、そして計算カウンタの増加をおこない、アルゴリズムが探索の進行に応じて継続的に学習・適応できるようにします。

//————————————————————————————————————————————————————————————————————
void C_AO_CMAES::Revision ()
{
  if (!revision) return;

  revision = false;

  // Sort the population by fitness
  SortPopulation ();

  // Update distribution parameters based on selected individuals
  UpdateDistribution ();

  // Update the calculation counter
  counteval++;
}
//————————————————————————————————————————————————————————————————————

C_AO_CMAESクラスのInitDistributionメソッドは、解探索の初期分布を初期化することを目的としています。このメソッドは、探索空間における各座標ごとに初期平均値cBを決定します。各i座標について、cB[i]の値はu.RNDfromCI()関数を用いてrangeMin[i]からrangeMax[i]の範囲内でランダムに設定され、各変数の有効範囲の中心付近に初期平均値が配置されます。外側のループは集団全体（popSize個体）を反復し、その内側のループは各個体の全座標を反復します。各子個体について座標が生成され、まずu.RNDfromCIによって指定範囲から一様にランダム生成され、その後SeInDiSpによってサンプリングステップに従って丸められます。

u.RNDfromCI 関数は、指定された区間内で一様分布する乱数を生成します。u.RNDfromCI関数は、与えられた区間内で一様分布に従う乱数を生成します。一方、u.SeInDiSp関数は、生成された各個体の座標がrangeMin[i]、rangeMax[i]、rangeStep[i]で定義される許容範囲内に収まることを保証します。

このメソッドは、探索空間全体に一様に分布したランダム解を用いて初期集団を生成します。また分布の中心cBもランダムに選択されます。これにより、CMA-ESの後続反復のための出発点が設定され、アルゴリズムはその後の探索過程で最適解へ向けて戦略を適応させていきます。

//+------------------------------------------------------------------+
//| Initialize search distribution                                   |
//+------------------------------------------------------------------+
void C_AO_CMAES::InitDistribution ()
{
  // Set the initial mean to the center of the search space
  for (int i = 0; i < coords; i++)
  {
    cB [i] = u.RNDfromCI (rangeMin [i], rangeMax [i]);
  }

  for (int k = 0; k < popSize; k++)
  {
    for (int i = 0; i < coords; i++)
    {
      // Generate a uniformly distributed point
      a [k].c [i] = u.RNDfromCI (rangeMin [i], rangeMax [i]);
      a [k].c [i] = u.SeInDiSp (a [k].c [i], rangeMin [i], rangeMax [i], rangeStep [i]);
    }
  }
}
//————————————————————————————————————————————————————————————————————

C_AO_CMAESクラスのGetChiNメソッドは、与えられた座標数(coords)に対して、N(0,I)から生成されるランダムベクトルの期待ノルムを計算します。この結果は、探索空間の次元数に基づいて近似的に求められます。

//+------------------------------------------------------------------+
//| Calculate the expected norm N(0,I)                               |
//+------------------------------------------------------------------+
double C_AO_CMAES::GetChiN ()
{
  double n = (double)coords;
  return MathSqrt (n) * (1.0 - 1.0 / (4.0 * n) + 1.0 / (21.0 * n * n));
}
//————————————————————————————————————————————————————————————————————

C_AO_CMAESクラスのSortPopulationメソッドは、目的関数値に基づいて解の集団をソートし、最良解を特定するために設計されています。まず、このメソッドは各個体のfitness値を集団からarfitness補助配列へコピーします。同時に、arindex配列にも各個体のインデックスが格納されます。これは、fitnessでソートした後でも、ソート結果と元の個体との対応関係を復元できるようにするためです。

ソートには挿入ソートアルゴリズムが使用されます。arfitness配列を順に走査し、各要素をソート済み部分の適切な位置へ挿入しながら、大きい要素を右側へシフトしていきます。重要なのは、このソートが降順（arfitness[j] < tempFitness条件）でおこなわれる点であり、目的はfitnessの最大化です。arfitnessと同時にarindexも同期して並べ替えられ、個体の元インデックス情報が追跡されます。

ソート完了後、最良個体のfitnessが現在記録されている最良解よりも優れているかがチェックされ、改善されている場合にはfBが更新されます。その際、最良解は集団内の最良個体の座標で置き換えられます。このように、SortPopulationメソッドは単に集団をfitness順に並べ替えるだけでなく、現在時点での最良解を継続的に追跡・更新する役割も担っています。

//+------------------------------------------------------------------+
//| Sort the population by fitness                                   |
//+------------------------------------------------------------------+
void C_AO_CMAES::SortPopulation ()
{
  // Copy fitness values and indices
  for (int i = 0; i < popSize; i++)
  {
    arindex [i] = i;
    arfitness [i] = a [i].f;
  }

  for (int i = 1; i < popSize; i++)
  {
    double tempFitness = arfitness [i];
    double tempIndex = arindex [i];
    int j = i - 1;

    // Sort in descending order (for maximization)
    while (j >= 0 && arfitness [j] < tempFitness)
    {
      arfitness [j + 1] = arfitness [j];
      arindex [j + 1] = arindex [j];
      j--;
    }

    arfitness [j + 1] = tempFitness;
    arindex [j + 1] = tempIndex;
  }

  // Update the best solution if necessary
  if (arfitness [0] > fB)
  {
    fB = arfitness [0];
    int bestIdx = (int)arindex [0];
    ArrayCopy (cB, a [bestIdx].c, 0, 0, coords);
  }
}
//————————————————————————————————————————————————————————————————————

C_AO_CMAESクラスのUpdateMeanメソッドは、最良個体群の重み付き再結合によって集団平均cBを更新します。各座標について、新しい平均値は最良個体群における各個体の座標値を重み付きで加算した総和として計算されます。ここで使用される重みはweights配列で定義されています。

//+------------------------------------------------------------------+
//| Update the mean using weighted recombination                     |
//+------------------------------------------------------------------+
void C_AO_CMAES::UpdateMean ()
{
  // Weighted recombination: m^(g+1) = Σ w_i * x_{i:λ}^(g+1)
  for (int j = 0; j < coords; j++)
  {
    double meanSum = 0.0;
    for (int i = 0; i < mu; i++)
    {
      int idx = (int)arindex [i];
      meanSum += weights [i] * a [idx].c [j];
    }
    cB [j] = meanSum;
  }
}
//————————————————————————————————————————————————————————————————————

C_AO_CMAESクラスのComputeWeightsメソッドは、重み付き再結合に用いるweights配列を計算します。まず、配列weightsをサイズmuで確保し、重み計算に使用するlog(mu+0.5)を事前計算します。次に、i=0からmu-1までについて、各重みをlog(mu+0.5)−log(i+1)として設定します。その後、すべての重みの総和を計算し、合計が1になるように各重みを正規化します。さらに、重みの二乗和も計算されます。これにより、mu_eff（有効選択質量）が算出され、アルゴリズムにおける更新速度や分散調整のチューニングに利用されます。

このメソッドは、最良解の寄与度を考慮した最適な重みを生成し、CMA-ESにおける再結合プロセスの基盤を提供します。

//+------------------------------------------------------------------+
//| Calculate weighted recombination weights                         |
//+------------------------------------------------------------------+
void C_AO_CMAES::ComputeWeights ()
{
  // Allocate the array of weights
  ArrayResize (weights, mu);

  // Pre-calculate log(mu + 0.5)
  double log_mu_plus_half = MathLog (mu + 0.5);

  // Calculate positive weights
  double sum = 0.0;
  for (int i = 0; i < mu; i++)
  {
    weights [i] = log_mu_plus_half - MathLog (i + 1);
    sum += weights [i];
  }

  // Normalize weights
  double sum_weights = 0.0;
  double sum_squares = 0.0;
  for (int i = 0; i < mu; i++)
  {
    weights [i] /= sum;
    sum_weights += weights [i];
    sum_squares += weights [i] * weights [i];
  }

  // Calculate the effective mass of the variance selection
  mu_eff = sum_weights * sum_weights / sum_squares;
}
//————————————————————————————————————————————————————————————————————

C_AO_CMAESクラスのUpdateDistributionメソッドは、CMA-ESアルゴリズムにおける分布パラメータ、すなわち共分散行列covMatrixおよびステップサイズsigmaを更新します。まず、十分な世代が経過している場合(counteval−eigeneval>eigenInterval)、固有値分解が実行され、現在の平均cBがoldMean一時配列に保存されます。その後、UpdateMeanメソッドが呼び出され新しい平均が計算されます。次に、新旧平均の差をsigmaで割った値を計算し、さらにC^(-1/2)との積を求めます。この計算はBおよびD分解を用いて段階的におこなわれ、その結果はinvsqrtC_times_ywに格納されます。

続いて、このinvsqrtC_times_ywを用いてステップサイズ用の進化経路psを更新します。その後、psノルムと期待長さに基づいて進捗が「停滞しているか」を判定します。さらに、pcおよびy_wと進捗指標hsigを用いて、共分散行列更新用の進化経路pcを更新します。その後、まずランク1更新行列c1aを計算し、続いて旧平均との差分（sigmaで正規化された値）とweightsに基づいてランクμ更新行列cmuを構築します。これらを用いて共分散行列covMatrixを更新し、進捗が停滞している場合にはc1が調整されます。また、共分散行列の正定性を維持するためにEnforcePositiveDefiniteが定期的に呼び出されます。さらに、ステップサイズsigmaはpsノルムに基づいて更新され、数値安定性のためmin_sigmaとmax_sigmaの範囲内に制限されます。

一般に、UpdateDistributionメソッドは進化経路と集団データという探索履歴に基づいて分布パラメータ（共分散行列およびステップサイズ）を調整し、目的関数の地形に適応できるようにする役割を担います。

//+------------------------------------------------------------------+ 
//| Update distribution parameters                                   |
//+------------------------------------------------------------------+
void C_AO_CMAES::UpdateDistribution ()
{
  // Check the necessity of eigenvalue decomposition
  if (counteval - eigeneval > eigenInterval)
  {
    ComputeEigendecomposition ();
    eigeneval = counteval;
  }

  // Preserve the old average
  double oldMean [];
  ArrayResize (oldMean, coords);
  ArrayCopy (oldMean, cB, 0, 0, coords);

  // Update average
  UpdateMean ();

  // Calculate the mean shift
  double y_w [];
  ArrayResize (y_w, coords);
  for (int j = 0; j < coords; j++)
  {
    y_w [j] = (cB [j] - oldMean [j]) / sigma;
  }

  // Calculate C^(-1/2) * y_w
  // Step 1: B^T * y_w
  ArrayInitialize (temp_vec, 0.0);
  for (int i = 0; i < coords; i++)
  {
    for (int j = 0; j < coords; j++)
    {
      temp_vec [i] += B [j * coords + i] * y_w [j];
    }
  }

  // Step 2: D^(-1) * (B^T * y_w)
  for (int i = 0; i < coords; i++)
  {
    temp_vec [i] /= D [i];
  }

  // Step 3: B * D^(-1) * B^T * y_w
  ArrayInitialize (invsqrtC_times_yw, 0.0);
  for (int i = 0; i < coords; i++)
  {
    for (int j = 0; j < coords; j++)
    {
      invsqrtC_times_yw [i] += B [i * coords + j] * temp_vec [j];
    }
  }

  // Update the evolution path for sigma
  double norm_ps_squared = 0.0;
  for (int i = 0; i < coords; i++)
  {
    ps [i] = (1.0 - cs) * ps [i] + MathSqrt (cs * (2.0 - cs) * mu_eff) * invsqrtC_times_yw [i];
    norm_ps_squared += ps [i] * ps [i];
  }

  // Heaviside function
  double norm_ps = MathSqrt (norm_ps_squared);
  double expected_length = MathSqrt (1.0 - MathPow (1.0 - cs, 2.0 * counteval)) * chiN;
  bool hsig = norm_ps / expected_length < hsig_threshold;

  // Update the evolution path for C
  double delta_hsig = hsig ? 1.0 : 0.0;
  for (int i = 0; i < coords; i++)
  {
    pc [i] = (1.0 - cc) * pc [i] + delta_hsig * MathSqrt (cc * (2.0 - cc) * mu_eff) * y_w [i];
  }

  // Prepare rank-1 update
  double c1a [];
  ArrayResize (c1a, coords * coords);
  for (int i = 0; i < coords; i++)
  {
    for (int j = 0; j <= i; j++)
    {
      c1a [i * coords + j] = c1a [j * coords + i] = pc [i] * pc [j];
    }
  }

  // Prepare rank-μ update
  double cmu [];
  ArrayResize (cmu, coords * coords);
  ArrayInitialize (cmu, 0.0);

  for (int k = 0; k < mu; k++)
  {
    int idx = (int)arindex [k];

    // Calculate y_i = (x_i - m_old) / sigma
    for (int i = 0; i < coords; i++)
    {
      temp_vec [i] = (a [idx].c [i] - oldMean [i]) / sigma;
    }

    // Add the weighted outer product
    for (int i = 0; i < coords; i++)
    {
      for (int j = 0; j <= i; j++)
      {
        double update = weights [k] * temp_vec [i] * temp_vec [j];
        cmu [i * coords + j] += update;
        if (i != j) cmu [j * coords + i] += update;
      }
    }
  }

  // Update C the covariance matrix
  double c1 = learningRateC1;
  double cmu_rate = learningRateCMu;

  // Adjust c1 if hsig is 'false' (progress stalled)
  if (!hsig)
  {
    c1 *= (1.0 - (1.0 - delta_hsig) * cc * (2.0 - cc));
  }

  double one_minus_c1_cmu = 1.0 - c1 - cmu_rate;

  // Update C with rank-1 and rank-μ updates
  for (int i = 0; i < coords; i++)
  {
    for (int j = 0; j <= i; j++)
    {
      covMatrix [i * coords + j] = one_minus_c1_cmu * covMatrix [i * coords + j] +
                                   c1 * c1a [i * coords + j] +
                                   cmu_rate * cmu [i * coords + j];

      // Maintain symmetry
      if (i != j)
      {
        covMatrix [j * coords + i] = covMatrix [i * coords + j];
      }
    }
  }

  // Ensure positive definiteness
  if (counteval % (10 * eigenInterval) == 0)
  {
    EnforcePositiveDefinite ();
  }

  //  Update the sigma step size
  double exponent = (cs / damps) * (norm_ps / chiN - 1.0);
  sigma *= MathExp (exponent);

  // Limit sigma for numerical stability
  double min_sigma = 1e-16;
  double max_eigenvalue = D [0]; // D sorted in descending order
  double max_sigma = 1e4 * MathMax (1.0, MathSqrt (max_eigenvalue));

  if (sigma < min_sigma) sigma = min_sigma;
  else
    if (sigma > max_sigma) sigma = max_sigma;
}
//————————————————————————————————————————————————————————————————————

C_AO_CMAESクラスのComputeEigendecompositionメソッドは、共分散行列covMatrixの対称分解を固有値および固有ベクトルへと変換するために、改良されたJacobi法を用いて分解を実行します。まず、このメソッドは元のcovMatrixを変更しないよう、そのコピーを作成します。初期状態では、Bは単位行列として設定され、初期基底を表します。その後、最大反復回数に達するか、または非対角要素がtolerance未満になるまで反復計算がおこなわれます。

非対角要素の最大値を求める：次に、非対角要素のうち絶対値が最大の要素を探索し、回転操作に使用する要素を選択します。もし最大値がtolerance未満であれば、反復処理は終了します。続いて、ヤコビ回転を実行するための角度phiを計算し、その回転を適用してC_copyの要素を更新します。同時に、回転に従ってB行列（固有ベクトル）の各列も更新されます。反復終了後、対角成分の平方根を計算し（数値安定性のため最小値1e−14を保証）、固有値を得ます。その後、固有値と対応する固有ベクトルを降順に並び替え、後続の計算で使用できるように整列します。このようにして本メソッドは、共分散行列の固有値および固有ベクトルを高精度に算出する役割を果たします。

//+------------------------------------------------------------------+
//| Calculate the eigenvalue decomposition using the Jacobi method   |
//+------------------------------------------------------------------+
void C_AO_CMAES::ComputeEigendecomposition ()
{
  // Create a copy of the covariance matrix for decomposition
  double C_copy [];
  ArrayResize (C_copy, coords * coords);
  ArrayCopy (C_copy, covMatrix);

  // Initialize B as the identity matrix
  for (int i = 0; i < coords; i++)
  {
    for (int j = 0; j < coords; j++)
    {
      B [i * coords + j] = (i == j) ? 1.0 : 0.0;
    }
  }

  // Improved Jacobi eigenvalue decomposition
  int max_iterations = 10; //50 * coords;
  double tolerance   = 0.01; //1e-14 * coords * coords;

  for (int iter = 0; iter < max_iterations; iter++)
  {
    // Find the largest off-diagonal element
    double max_val = 0.0;
    int p = 0, q = 1;

    for (int i = 0; i < coords - 1; i++)
    {
      for (int j = i + 1; j < coords; j++)
      {
        double val = MathAbs (C_copy [i * coords + j]);
        if (val > max_val)
        {
          max_val = val;
          p = i;
          q = j;
        }
      }
    }

    // Check convergence
    if (max_val < tolerance) break;

    // Calculate the rotation angle
    double app = C_copy [p * coords + p];
    double aqq = C_copy [q * coords + q];
    double apq = C_copy [p * coords + q];

    double phi = 0.5 * MathArctan (2.0 * apq / (aqq - app + 1e-14));
    double c = MathCos (phi);
    double s = MathSin (phi);

    // Update the matrix elements
    double app_new = c * c * app - 2 * c * s * apq + s * s * aqq;
    double aqq_new = s * s * app + 2 * c * s * apq + c * c * aqq;

    C_copy [p * coords + p] = app_new;
    C_copy [q * coords + q] = aqq_new;
    C_copy [p * coords + q] = C_copy [q * coords + p] = 0.0;

    // Update other elements in p and q rows/columns
    for (int i = 0; i < coords; i++)
    {
      if (i != p && i != q)
      {
        double aip = C_copy [i * coords + p];
        double aiq = C_copy [i * coords + q];
        C_copy [i * coords + p] = C_copy [p * coords + i] = c * aip - s * aiq;
        C_copy [i * coords + q] = C_copy [q * coords + i] = s * aip + c * aiq;
      }
    }

    // Update eigenvectors
    for (int i = 0; i < coords; i++)
    {
      double bip = B [i * coords + p];
      double biq = B [i * coords + q];
      B [i * coords + p] = c * bip - s * biq;
      B [i * coords + q] = s * bip + c * biq;
    }
  }

  // Extract eigenvalues and ensure positivity
  double min_eigenvalue = 1e-14;
  for (int i = 0; i < coords; i++)
  {
    D [i] = MathSqrt (MathMax (min_eigenvalue, C_copy [i * coords + i]));
  }

  // Sort eigenvalues and vectors in descending order
  for (int i = 0; i < coords - 1; i++)
  {
    int max_idx = i;
    for (int j = i + 1; j < coords; j++)
    {
      if (D [j] > D [max_idx]) max_idx = j;
    }

    if (max_idx != i)
    {
      // Exchange eigenvalues
      double temp = D [i];
      D [i] = D [max_idx];
      D [max_idx] = temp;

      // Exchange eigenvectors
      for (int k = 0; k < coords; k++)
      {
        temp = B [k * coords + i];
        B [k * coords + i] = B [k * coords + max_idx];
        B [k * coords + max_idx] = temp;
      }
    }
  }
}
//————————————————————————————————————————————————————————————————————

C_AO_CMAESクラスのCheckPositiveDefiniteメソッドは、共分散行列が正定値であるかどうかを確認します。まず、対角成分の正値性を素早くチェックし、さらにD配列（降順にソートされた固有値）の最小固有値が小さな正の値1e−14より大きいかを比較します。これら2つの条件を満たす場合、メソッドはtrueを返します。

//+------------------------------------------------------------------+
//| Check the positive definiteness of the covariance matrix         |
//+------------------------------------------------------------------+
bool C_AO_CMAES::CheckPositiveDefinite ()
{
  // Quick check: all diagonal elements should be positive
  for (int i = 0; i < coords; i++)
  {
    if (covMatrix [i * coords + i] <= 0) return false;
  }

  // Check if the minimum eigenvalue is positive
  double min_eigenvalue = D [coords - 1]; // D is sorted in descending order
  return min_eigenvalue > 1e-14;
}
//————————————————————————————————————————————————————————————————————

C_AO_CMAESクラスのEnforcePositiveDefiniteメソッドは、共分散行列covMatrixの正定値性を保証するために設計されています。まず、対角要素の最小値を求め、不足している「補正値」を対角成分に加算することで、最小対角要素が1e−10になるよう調整します。次に、各オフ対角要素を対応する転置要素と平均化することで、行列の対称性を最大化します。

それでも行列が正定値でない場合は、ComputeEigendecompositionメソッドを用いて固有値分解を実行します。その後、1e−10の平方根より小さいすべての固有値を√1e−10に置き換え、正の値を保証します。続いて、修正された固有値Dと固有ベクトルBを用いてcovMatrixを再構成し、補正済みの正定値行列を生成します。このようにして本メソッドは、共分散行列が常に正定値であり続けることを保証し、CMA-ESアルゴリズムにおける後続の計算が安定して実行できるようにします。

//+------------------------------------------------------------------+
//| Ensuring positive definiteness of the covariance matrix          |
//+------------------------------------------------------------------+
void C_AO_CMAES::EnforcePositiveDefinite ()
{
  // Method 1: Adding a small value to the diagonal
  double min_diag = 1e308; // A very large number
  for (int i = 0; i < coords; i++)
  {
    if (covMatrix [i * coords + i] < min_diag)
    {
      min_diag = covMatrix [i * coords + i];
    }
  }

  if (min_diag < 1e-10)
  {
    double correction = 1e-10 - min_diag;
    for (int i = 0; i < coords; i++)
    {
      covMatrix [i * coords + i] += correction;
    }
  }

  // Method 2: Ensuring symmetry
  for (int i = 0; i < coords; i++)
  {
    for (int j = i + 1; j < coords; j++)
    {
      double avg = (covMatrix [i * coords + j] + covMatrix [j * coords + i]) * 0.5;
      covMatrix [i * coords + j] = covMatrix [j * coords + i] = avg;
    }
  }

  // If still not positive definite, perform the factorization and fix
  if (!CheckPositiveDefinite ())
  {
    ComputeEigendecomposition ();

    double min_eigenvalue = 1e-10;
    for (int i = 0; i < coords; i++)
    {
      if (D [i] < MathSqrt (min_eigenvalue))
      {
        D [i] = MathSqrt (min_eigenvalue);
      }
    }

    // Reconstruct C = B * D^2 * B^T
    ArrayInitialize (covMatrix, 0.0);
    for (int i = 0; i < coords; i++)
    {
      for (int j = 0; j <= i; j++)
      {
        double sum = 0.0;
        for (int k = 0; k < coords; k++)
        {
          sum += B [i * coords + k] * D [k] * D [k] * B [j * coords + k];
        }
        covMatrix [i * coords + j] = covMatrix [j * coords + i] = sum;
      }
    }
  }
}
//————————————————————————————————————————————————————————————————————

テスト結果

さて、いよいよ結果を見ていきます。アルゴリズムは、中規模の問題に対しては非常に良好かつ高速に収束することがすぐに分かります。また、設定によっては低次元問題でも良好な性能を示します。しかし、多次元の高難度問題になると、計算時間が不釣り合いに増大するため、実質的に結果評価からは除外されています。なぜこのような現象が起きるのでしょうか。従来、MQL5における行列計算は深刻な計算負荷の問題を抱えています。行列演算を手動で実装した場合、高次元問題に対してアルゴリズムの性能は著しく低下し、実用性が大きく制限されます。しかし、組み込みの行列処理クラスが導入されたことで状況は大きく改善されました。現在では、計算負荷の高い処理についてはプラットフォーム標準の行列演算実装を使用することが極めて重要となっています。

CMAES|Covariance Matrix Adaptation Evolution Strategy|50.0|25.0|0.01|0.8|1.0|
=============================
5 Hilly's; Func runs:10000; result:0.7625797883550075
25 Hilly's; Func runs:10000; result:0.7208874560706138
500 Hilly's; Func runs:10000; result:0.0
=============================
5 Forest's; Func runs:10000; result:0.8205636421348295
25 Forest's; Func runs:10000; result:0.7961602627346933
500 Forest's; Func runs:10000; result:0.0
=============================
5 Megacity's; Func runs:10000; result:0.7584615384615383
25 Megacity's; Func runs:10000; result:0.49076923076923074
500 Megacity's; Func runs:10000; result:0.0
=============================
All score:4.34942 (48.33%)

可視化では、アルゴリズムが長距離ジャンプと関数の極値への集中の両方を行い、それらを徹底的に探索していることが確認できます。

Hilly

Hillyテスト関数でのCMA-ES

Forest

Forestテスト関数でのCMA-ES

Megacity

Megacity関数でのCMA-ES

テスト結果に基づくと、CMA-ESアルゴリズムは集団ベース最適化アルゴリズムのランキング表で第38位に位置しています。

#	AO	詳細	Hilly			Hilly ファイナル	Forest			Forest ファイナル	Megacity（離散）			Megacity ファイナル	ファイナル結果	％ MAX
#	AO	詳細	10 p (5 F)	50p(25F)	1000p(500F)	Hilly ファイナル	10 p (5 F)	50p(25F)	1000p(500F)	Forest ファイナル	10 p (5 F)	50p(25F)	1000p(500F)	Megacity ファイナル	ファイナル結果	％ MAX
1	ANS	近隣検索	0.94948	0.84776	0.43857	2.23581	1.00000	0.92334	0.39988	2.32323	0.70923	0.63477	0.23091	1.57491	6.134	68.15
2	CLA	コードロックアルゴリズム(joo)	0.95345	0.87107	0.37590	2.20042	0.98942	0.91709	0.31642	2.22294	0.79692	0.69385	0.19303	1.68380	6.107	67.86
3	AMOm	動物移動最適化M	0.90358	0.84317	0.46284	2.20959	0.99001	0.92436	0.46598	2.38034	0.56769	0.59132	0.23773	1.39675	5.987	66.52
4	(P+O)ES	(P+O)進化戦略	0.92256	0.88101	0.40021	2.20379	0.97750	0.87490	0.31945	2.17185	0.67385	0.62985	0.18634	1.49003	5.866	65.17
5	CTA	彗星の尾アルゴリズム(joo)	0.95346	0.86319	0.27770	2.09435	0.99794	0.85740	0.33949	2.19484	0.88769	0.56431	0.10512	1.55712	5.846	64.96
6	TETA	時間進化移動アルゴリズム(joo)	0.91362	0.82349	0.31990	2.05701	0.97096	0.89532	0.29324	2.15952	0.73462	0.68569	0.16021	1.58052	5.797	64.41
7	SDSm	確率的拡散探索M	0.93066	0.85445	0.39476	2.17988	0.99983	0.89244	0.19619	2.08846	0.72333	0.61100	0.10670	1.44103	5.709	63.44
8	BOAm	ビリヤード最適化アルゴリズムM	0.95757	0.82599	0.25235	2.03590	1.00000	0.90036	0.30502	2.20538	0.73538	0.52523	0.09563	1.35625	5.598	62.19
9	AAm	アーチェリーアルゴリズムM	0.91744	0.70876	0.42160	2.04780	0.92527	0.75802	0.35328	2.03657	0.67385	0.55200	0.23738	1.46323	5.548	61.64
10	ESG	社会集団の進化(joo)	0.99906	0.79654	0.35056	2.14616	1.00000	0.82863	0.13102	1.95965	0.82333	0.55300	0.04725	1.42358	5.529	61.44
11	SIA	等方的焼きなまし(joo)	0.95784	0.84264	0.41465	2.21513	0.98239	0.79586	0.20507	1.98332	0.68667	0.49300	0.09053	1.27020	5.469	60.76
12	BBO	生物地理学に基づく最適化	0.94912	0.69456	0.35031	1.99399	0.93820	0.67365	0.25682	1.86867	0.74615	0.48277	0.17369	1.40261	5.265	58.50
13	ACS	人工協調探索	0.75547	0.74744	0.30407	1.80698	1.00000	0.88861	0.22413	2.11274	0.69077	0.48185	0.13322	1.30583	5.226	58.06
14	DA	弁証法的アルゴリズム	0.86183	0.70033	0.33724	1.89940	0.98163	0.72772	0.28718	1.99653	0.70308	0.45292	0.16367	1.31967	5.216	57.95
15	BHAm	ブラックホールアルゴリズムM	0.75236	0.76675	0.34583	1.86493	0.93593	0.80152	0.27177	2.00923	0.65077	0.51646	0.15472	1.32195	5.196	57.73
16	ASO	無政府社会の最適化	0.84872	0.74646	0.31465	1.90983	0.96148	0.79150	0.23803	1.99101	0.57077	0.54062	0.16614	1.27752	5.178	57.54
17	RFO	ロイヤルフラッシュ最適化(joo)	0.83361	0.73742	0.34629	1.91733	0.89424	0.73824	0.24098	1.87346	0.63154	0.50292	0.16421	1.29867	5.089	56.55
18	AOSm	原子軌道探索M	0.80232	0.70449	0.31021	1.81702	0.85660	0.69451	0.21996	1.77107	0.74615	0.52862	0.14358	1.41835	5.006	55.63
19	TSEA	亀の甲羅進化アルゴリズム(joo)	0.96798	0.64480	0.29672	1.90949	0.99449	0.61981	0.22708	1.84139	0.69077	0.42646	0.13598	1.25322	5.004	55.60
20	DE	差分進化	0.95044	0.61674	0.30308	1.87026	0.95317	0.78896	0.16652	1.90865	0.78667	0.36033	0.02953	1.17653	4.955	55.06
21	SRA	レストラン経営達人アルゴリズム(joo)	0.96883	0.63455	0.29217	1.89555	0.94637	0.55506	0.19124	1.69267	0.74923	0.44031	0.12526	1.31480	4.903	54.48
22	CRO	化学反応の最適化	0.94629	0.66112	0.29853	1.90593	0.87906	0.58422	0.21146	1.67473	0.75846	0.42646	0.12686	1.31178	4.892	54.36
23	BIO	血液型遺伝最適化(joo)	0.81568	0.65336	0.30877	1.77781	0.89937	0.65319	0.21760	1.77016	0.67846	0.47631	0.13902	1.29378	4.842	53.80
24	BSA	鳥群アルゴリズム	0.89306	0.64900	0.26250	1.80455	0.92420	0.71121	0.24939	1.88479	0.69385	0.32615	0.10012	1.12012	4.809	53.44
25	HS	ハーモニーサーチ	0.86509	0.68782	0.32527	1.87818	0.99999	0.68002	0.09590	1.77592	0.62000	0.42267	0.05458	1.09725	4.751	52.79
26	SSG	苗木の播種と成長	0.77839	0.64925	0.39543	1.82308	0.85973	0.62467	0.17429	1.65869	0.64667	0.44133	0.10598	1.19398	4.676	51.95
27	BCOm	細菌走化性最適化M	0.75953	0.62268	0.31483	1.69704	0.89378	0.61339	0.22542	1.73259	0.65385	0.42092	0.14435	1.21912	4.649	51.65
28	ABO	アフリカ水牛の最適化	0.83337	0.62247	0.29964	1.75548	0.92170	0.58618	0.19723	1.70511	0.61000	0.43154	0.13225	1.17378	4.634	51.49
29	(PO)ES	(PO)進化戦略	0.79025	0.62647	0.42935	1.84606	0.87616	0.60943	0.19591	1.68151	0.59000	0.37933	0.11322	1.08255	4.610	51.22
30	FBA	フラクタルベースのアルゴリズム	0.79000	0.65134	0.28965	1.73099	0.87158	0.56823	0.18877	1.62858	0.61077	0.46062	0.12398	1.19537	4.555	50.61
31	TSm	タブーサーチM	0.87795	0.61431	0.29104	1.78330	0.92885	0.51844	0.19054	1.63783	0.61077	0.38215	0.12157	1.11449	4.536	50.40
32	BSO	ブレインストーミング最適化	0.93736	0.57616	0.29688	1.81041	0.93131	0.55866	0.23537	1.72534	0.55231	0.29077	0.11914	0.96222	4.498	49.98
33	WOAm	ウェール最適化アルゴリズムM	0.84521	0.56298	0.26263	1.67081	0.93100	0.52278	0.16365	1.61743	0.66308	0.41138	0.11357	1.18803	4.476	49.74
34	AEFA	人工電界アルゴリズム	0.87700	0.61753	0.25235	1.74688	0.92729	0.72698	0.18064	1.83490	0.66615	0.11631	0.09508	0.87754	4.459	49.55
35	AEO	人工生態系ベースの最適化アルゴリズム	0.91380	0.46713	0.26470	1.64563	0.90223	0.43705	0.21400	1.55327	0.66154	0.30800	0.28563	1.25517	4.454	49.49
36	CAm	ラクダアルゴリズムM	0.78684	0.56042	0.35133	1.69859	0.82772	0.56041	0.24336	1.63149	0.64846	0.33092	0.13418	1.11356	4.444	49.37
37	ACOm	アリコロニー最適化M	0.88190	0.66127	0.30377	1.84693	0.85873	0.58680	0.15051	1.59604	0.59667	0.37333	0.02472	0.99472	4.438	49.31
38	CMAES	共分散行列適応進化戦略	0.76258	0.72089	0.00000	1.48347	0.82056	0.79616	0.00000	1.61672	0.75846	0.49077	0.00000	1.24923	4.349	48.33
39	BFO-GA	細菌採餌最適化 - ga	0.89150	0.55111	0.31529	1.75790	0.96982	0.39612	0.06305	1.42899	0.72667	0.27500	0.03525	1.03692	4.224	46.93
40	SOA	シンプル最適化アルゴリズム	0.91520	0.46976	0.27089	1.65585	0.89675	0.37401	0.16984	1.44060	0.69538	0.28031	0.10852	1.08422	4.181	46.45
41	ABHA	人工蜂の巣アルゴリズム	0.84131	0.54227	0.26304	1.64663	0.87858	0.47779	0.17181	1.52818	0.50923	0.33877	0.10397	0.95197	4.127	45.85
42	ACMO	大気雲モデルの最適化	0.90321	0.48546	0.30403	1.69270	0.80268	0.37857	0.19178	1.37303	0.62308	0.24400	0.10795	0.97503	4.041	44.90
43	ADAMm	適応モーメント推定M	0.88635	0.44766	0.26613	1.60014	0.84497	0.38493	0.16889	1.39880	0.66154	0.27046	0.10594	1.03794	4.037	44.85
44	CGO	カオスゲーム最適化	0.57256	0.37158	0.32018	1.26432	0.61176	0.61931	0.62161	1.85267	0.37538	0.21923	0.19028	0.78490	3.902	43.35
45	CROm	サンゴ礁最適化M	0.78512	0.46032	0.25958	1.50502	0.86688	0.35297	0.16267	1.38252	0.63231	0.26738	0.10734	1.00703	3.895	43.27
	RW	ランダムウォーク	0.48754	0.32159	0.25781	1.06694	0.37554	0.21944	0.15877	0.75375	0.27969	0.14917	0.09847	0.52734	2.348	26.09

まとめ

CMA-ESの持つ、目的関数の局所的な幾何構造に適応して探索戦略を調整する能力は、構造が未知で条件の悪い複雑な問題に対してほぼ不可欠なものとなっています。冪分布の裾は重いため、アルゴリズムは長距離ジャンプできるようになります。これは、局所最適から脱出する上で重要な役割を果たします。一方で、その計算量はメモリO(n²)、計算時間O(n³)と非常に大きく、高次元問題への適用性を著しく制限します。特にn>100のような高次元領域では、その計算コストが得られる利点を上回るため、実用性が大きく低下します。多次元関数における実行時間は急激に増加し、大規模問題では事実上使用不可能になります。

それにもかかわらず、CMA-ESはノイズを含む関数、非連続な目的関数、多峰性の地形に対しても有効に機能します。この柔軟性の本質は、アルゴリズムが関数構造について事前の仮定を置かず、蓄積された経験に基づいて探索戦略を適応させながら慎重に探索するという基本思想にあります。CMA-ESは進化計算の見取り図そのものを変えたアルゴリズムであり、生物学的に着想された手法が厳密な数学的基盤を持ち得ること、そしてアルゴリズムの適応が単なるパラメータ調整ではなく「問題構造そのものの学習」であることを示しています。

このアルゴリズムは中規模問題に対して非常に優れた性能を発揮し、その分野では有力な選択肢となる専門的ツールです。数学的な美しさと実用的な効率性を兼ね備えていますが、計算コストの制約を考慮した慎重な運用が必要です。なお、近年MQL5開発者によって高速な行列演算機能が導入されており、本実装では従来型の逐次計算が使用されています。

図2：対応するテストに応じたアルゴリズムのカラーグラデーション

チャート

図3：アルゴリズムテスト結果のヒストグラム（0から100のスケール、高いほど良い）100は理論上の最大値であり、アーカイブには評価表を計算するためのスクリプトがあります。

CMAESの長所と短所

長所

中次元関数において良好な収束性を示す

短所

低次元関数では結果のばらつきが大きい
大規模な問題ではリソース消費が極めて大きい

この記事には、最新版のアルゴリズムコードを含むアーカイブが添付されています。記事の著者は、古典的なアルゴリズムの説明の絶対的な正確さについて責任を負いません。探索性能を向上させるために、それらの多くに変更が加えられています。記事に示された結論と判断は、実験結果に基づいています。

github：https://github.com/JQSakaJoo/Population-optimization-algorithms-MQL5

記事で使用されているプログラム

#	名前	種類	詳細
1	#C_AO.mqh	インクルード	集団最適化アルゴリズムの親クラス
2	#C_AO_enum.mqh	インクルード	集団最適化アルゴリズムの列挙
3	TestFunctions.mqh	インクルード	テスト関数のライブラリ
4	TestStandFunctions.mqh	インクルード	テストスタンド関数ライブラリ
5	Utilities.mqh	インクルード	補助関数のライブラリ
6	CalculationTestResults.mqh	インクルード	比較表の結果を計算するスクリプト
7	Testing AOs.mq5	スクリプト	すべての集団最適化アルゴリズムの統一テストスタンド
8	Simple use of population optimization algorithms.mq5	スクリプト	可視化せずに集団最適化アルゴリズムを使用する簡単な例
9	Test_AO_CMAES.mq5	スクリプト	CMAESテストベンチ

MetaQuotes Ltdによってロシア語から翻訳されました。
元の記事: https://www.mql5.com/ru/articles/18227

添付されたファイル |

ZIPをダウンロード

CMAES.zip (230.9 KB)

警告: これらの資料についてのすべての権利はMetaQuotes Ltd.が保有しています。これらの資料の全部または一部の複製や再プリントは禁じられています。

この記事はサイトのユーザーによって執筆されたものであり、著者の個人的な見解を反映しています。MetaQuotes Ltdは、提示された情報の正確性や、記載されているソリューション、戦略、または推奨事項の使用によって生じたいかなる結果についても責任を負いません。

Andrey Dik

I WILL CONSIDER PROPOSALS FOR THE PUBLICATION OF A BOOK (TEXTBOOK) ON OPTIMIZATION ALGORITHMS.
A group for communication on optimization and free product testing://t.me/+vazsAAcney4zYmZi
Attention! My Telegram doppelgangers have appeared, my real nickname is @JQS_aka_Joo
My github with optimization algorithms: https://github.com/JQSakaJoo/Population-optimization-algorithms-MQL5
All my publications: https://www.mql5.com/en/users/joo/publications
I have been developing systems based on machine learning technologies since 2007 and in the field of artificial
intelligence, optimization and forecasting.

この著者による他の記事

ディスカッションに移動

機械学習を用いたフラクタル市場構造入門

本記事では、金融時系列を自己相似的なフラクタル構造という観点から考察します。市場の価格変動が自己相似フラクタルとして捉えられる可能性を支持する類似性が多数存在することから、このような構造の予測可能性の地平線について考えることができます。

市場シミュレーション（第24回）：SQL入門(VII)

前回の記事では、SQLに必要な導入を完了しました。SQLについて何を説明したいのかは、十分に明確にできたと思います。これは、市場のリプレイ/シミュレーションシステムの構築を見に来る人であれば誰でも、そこで何が起きているのかを少なくともある程度イメージできるようにするためのものでした。重要な点は、SQLが完全に処理できることをわざわざプログラミングする意味はないということです。

エラー 146 (「トレードコンテキストビジー」) と、その対処方法

この記事では、MT4において複数のEAの衝突をさける方法を扱います。ターミナルの操作、MQL4の基本的な使い方がわかる人にとって、役に立つでしょう。

初級から中級まで：関数ポインタ

プログラミングにおいて「ポインタ」という概念について聞いたことがある方は多いと思います。しかし、MQL5でもこの種のデータを利用できることをご存じでしょうか。もちろん、これは実行時の挙動を適切に制御し、不可解な動作を回避できるよう慎重に扱う必要があります。ただし、この機能は非常に限定的な用途に向けたものであり、特定のタスクに特化しているため、MQL5におけるポインタの概念や使用方法について語られる機会はあまり多くありません。